Macroscopic properties of heterogeneous materials, including thermal and electrical conductivity, elastic moduli, and fluid permeability, depend crucially on their microstructure. Optimal design of such materials with desired properties, as well as characterization of their microstructure, are fundamental problems that have been studied for a long time. Accurate characterization of the microstructure is possible if several $n$-point correlation functions that describe their statistical properties can be computed. While the limits $n=1$, corresponding to phase fractions, and $n=2$, that represents such two-point correlation functions as the radial distribution function, have yielded useful information and insights and have been utilized for reconstruction of models of heterogeneous materials, in many cases higher-order correlation functions are required in order to develop deeper understanding of materials' properties, as well as obtaining accurate estimates for them. We describe an algorithm for computing third-order correlation functions. We test the accuracy of the algorithm for a model of fully penetrable disks, i.e., the so-called cherry-pit model in the limit of the penetrability parameter $\ensuremath{\lambda}=0$, for which an exact expression for the three-point probability function ${S}_{3}({\mathbf{x}}_{1},{\mathbf{x}}_{2},{\mathbf{x}}_{3})$ is known. Excellent agreement between the computed results and the theoretical predictions is demonstrated. We then report on the results of extensive computations for several types of heterogeneous materials and the analysis of their three-point correlation functions. They include the probability and cluster functions ${S}_{3}$ and ${C}_{3}$, as well as the three-point surface and surface-surface-void correlation functions ${F}_{sss}$ and ${F}_{ssv}$, for a variety of two- and three-dimensional disordered materials and media.