Journal of Microscopy, Vol. 258, Issue 3 2015, pp. 173–178
doi: 10.1111/jmi.12234
Received 30 September 2014; accepted 23 January 2015
Reconstruction of random heterogeneous media F. BALLANI & D. STOYAN Institut f¨ur Stochastik, Technische Universit¨at Bergakademie Freiberg, Freiberg, Germany
Key words. Random heterogeneous media, simulated annealing, stochastic reconstruction.
Summary Stochastic reconstruction is a technique to generate samples of random structures with prescribed distributional properties in the sense that certain of their statistical summary characteristics match target values or forms. This technique can be used to produce structures of any wanted size for further statistical analysis starting from small samples, which may be even only lower dimensional, for instance, when three-dimensional imaging techniques are not available. In this review we explain the main ideas of stochastic reconstruction, concentrating on the most important case of digitized binary media and with a particular emphasis on stereological reconstruction.
Introduction Stochastic reconstruction is a technique to generate samples of random structures, typically in the sense that the obtained samples, in some respects, are very similar to a given real sample. This includes the generation of similar and possibly larger planar patterns of a random heterogeneous medium starting from a planar pattern, likewise that of a spatial pattern of a random heterogeneous medium from a spatial or planar pattern, as well as the reconstruction of point patterns. In this review, we explain the main ideas of this technique, where we restrict to the special case that samples of a three-dimensional medium consisting of two phases (implying its heterogeneity) have to be reconstructed. Throughout the paper we assume that the random medium is statistically homogeneous (spatially stationary) and isotropic (see, e.g. Torquato, 2002). We expect some knowledge of fundamental ideas of spatial statistics as, e.g. presented in Chiu et al. (2013). In the analysis and characterization of random heterogeneous media it is a fundamental task to make inference on an appropriate three-dimensional representation of a real-world specimen in a computer (typically a pixel image), for instance to get information on topological properties and on characteristics of connectivity, flow and transport properties, etc. Such ¨ Stochastik, Technische Universit¨at Correspondence to: F. Ballani, Institut fur Bergakademie Freiberg, D-09596 Freiberg, Germany. Tel: +49 (0)3731 39-2796; fax: +49 (0)3731 39-3598; e-mail: [email protected] C 2015 The Authors C 2015 Royal Microscopical Society Journal of Microscopy
virtual structures must be large enough to obtain reliable estimates of complicated characteristics. However, very often they are not directly available: samples obtained by modern imaging methods such as X-ray computed tomography may be too small and data from planar sections like scanning electron microscope images, which are often the only alternative in case three-dimensional imaging techniques are not available, are by their nature not sufficient to estimate spatial topological characteristics. Stochastic reconstruction is able to produce the three-dimensional samples needed, starting from (small) samples, which may be even only lower dimensional, resulting from planar (or even linear) sections. It can yield many samples of arbitrary size in full dimensionality. The classical approach to generate samples of threedimensional sets as described in Chiu et al. (2013) consists in fitting a model (like a Boolean model or a thresholded Gaussian random field) to the data, where the model parameters come somehow from the data, typically from estimated summary characteristics. Then the model is simulated according to its properties. For the nonspecialist this approach has the disadvantage to use complicated statistical techniques and the need to believe in a model she/he does perhaps not understand. Modern stochastic reconstruction is more general and (presumably) user-friendly since it does not use explicit models but only models implicitly determined by the algorithm. The user has only to understand the general idea and some functional summary characteristics we explain in the next section. The general idea is simple: construct a three-dimensional structure which has summary characteristics that match exactly (or in very good approximation) the prescribed ones, which come from statistics of the originally given sample. It is even possible to build-in explicitly the given data in a conditional approach. Though in every simulation of stochastic models given statistical data there is a reconstructive aspect (in the sense that ‘reconstruction means to build up a complete structure of which one has only a few parts or only partial evidence’, see Hornby, 1974) the term ‘stochastic reconstruction’ is used for this special method – maybe a clever marketing brand. According to Torquato (2002, section 8.2) an origin of stochastic reconstruction may be seen in ideas of statistics of Gaussian random fields. However, the modern reconstruction
174
F. BALLANI AND D. STOYAN
approach does not use such models and can be understood without deeper knowledge of probability theory. This modern approach goes back to a paper by Farmer in 1992, the Ph.D. thesis of Deutsch in 1992 and several subsequent publications, e.g. Deutsch & Cockerham (1994) (see the references in Hazlett, 1997), and two papers of Yeong & Torquato (1998a,b). Since the year 1992 there have been published many papers on various methods of reconstruction, about variants of the classical procedure, about the use of various summary characteristics and directed to various structural properties of the random media to reconstruct. A lot of references can be found in Torquato (2002), Chiu et al. (2013) and Li (2014). In the following, we first explain some spatial summary characteristics which are commonly used in stochastic reconstruction of two-phase structures. Then we sketch the reconstruction algorithm for three-dimensional statistically homogeneous (stationary) and isotropic two-phase structure in case it is represented as a binary pixel image, where 0-pixels stand for pore, and 1-pixels for solid (or for two different phases). Since stochastic reconstruction is a computerintensive method we also list some practical aspects which might be helpful to implement it in a fairly fast version. Finally, we discuss in some detail the stereological problem where a planar section is given and based on that a spatial sample has to be reconstructed, and consider even an example, which yields interesting insights. Summary characteristics The structure to be reconstructed is characterized by various summary characteristics. These are estimated from the original sample, based on scattering methods, X-ray computed tomography or microscopy. The choice of the summary characteristics depends on the available data, on the nature of the information aimed to get by the reconstruction, and on the associated computational costs. The main characteristic of, say, the solid phase Z of a statistically homogeneous structure is volume fraction VV (or porosity 1 − VV ), which is usually a priori matched exactly since only structures with the prescribed volume fraction VV are admitted. The further characteristics are usually of a functional nature. Most popular and important is the covariance C (r ) (or twopoint probability function) of the solid phase Z , which is the probability that two points a distance r apart are both covered by Z . (Note that in a two-phase structure the covariance of one phase unambiguously determines that of the other phase: the covariance of the complementary phase of Z equals 1 − 2C (0) + C (r ).) Using C (r ) in stochastic reconstruction is usually considered a must since it characterizes the second-order behaviour of a structure, in particular, besides VV = C (0), it also determines the specific surface area SV since its slope at the origin is related to SV via SV = −4C (0) (Chiu et al., 2013, p. 219). The function C (r ) has also the advantage
C(r)
Z
r Hd(r)
r L(r) r
Hl(r) r Fig. 1. Planar section together with a schematic illustration of events contributing to summary characteristics available from planar sections, for covariance C (r ), chord length distribution function L (r ) and linear and discoidal contact distribution function Hl (r ) and Hd (r ), see the text for further details.
that it may be determined even from small-angle scattering experiments. Further functional characteristics (that are often incorrectly called ‘correlation functions’) are, e.g. the chord-length distribution function, various contact distribution (or, equivalently, empty space) functions and the two-point cluster function, see Chiu et al. (2013) for a detailed introduction. The chord length distribution function L (r ) (not to be confused with the popular L -function derived from Ripley’s K function for random point patterns, see Chiu et al., 2013) related to the solid phase Z is the portion of Z -chords with length less or equal than r , where a Z -chord is a segment between intersection points of a line with the two-phase interface fully contained in Z , see also Figure 1. Hence, L (r ) contains information on connectedness and correlation along a lineal path. The same information is carried by the linear contact distribution function Hl (r ) of the complement of Z as well as by the lineal path function (see Torquato, 2002, p. 44). Particularly apparent for L (r ), Figure 1 suggests that the estimation of L (r ) (but also that of any other summary characteristic depending on some distance r ) might be heavily influenced by edge effects. Often a corrected estimation accounting for such edge ¨ effects is meaningful, see Ohser & Mucklich (2000, p. 102) for L (r ) and Chiu et al. (2013) for general ideas and further references. It should be noted that all summary characteristics considered so far are of a one-dimensional nature. This has C 2015 The Authors C 2015 Royal Microscopical Society, 258, 173–178 Journal of Microscopy
RECONSTRUCTION OF RANDOM HETEROGENEOUS MEDIA
the advantage that they can be determined even from onedimensional samples, and the disadvantage that they represent only simpler spatial properties of the structure. Nevertheless, Li (2014) states that the ‘lineal’ summary characteristics are percolation-sensitive descriptors. As a summary characteristic of a two-dimensional nature, but still available from planar sections, is the discoidal contact distribution function Hd (r ). The function Hd (r ) for the solid phase Z is the probability that a random two-dimensional disc of radius r with centre outside Z hits Z . Thus, Hd (r ) characterizes the planar expansion of regions in the void phase. An analogous function is the quadratic contact distribution function Hq (r ), where the disc is replaced by a square. Finally, if three-dimensional data are available, also summary characteristics with a truly three-dimensional nature might be used. One numerical characteristic of this type is the specific Euler number NV , carrying information on connectivity. Sometimes three-point correlation functions are used for stochastic reconstruction. A further functional summary characteristic leading to good reconstruction results (see Jiao et al., 2009) is the two-point cluster function (denoted by C 2 (r ) in Torquato, 2002, and by C 1 (r ) in Chiu et al., 2013). For the solid phase Z this is the probability that two points a distance r apart are both covered by Z and belong to the same cluster/connectivity component, i.e. there is a connecting path fully covered by Z . With the exception of the case of particle structures with well-separated particles, the determination of C 2 (r ) is rather difficult and time consuming. Simpler three-dimensional summary characteristics are contact distributions with three-dimensional structuring elements. The best known example is the spherical contact distribution function Hs (r ). The function Hs (r ) for the solid phase Z is the probability that a random ball of radius r with centre outside Z hits Z . Hence, in 2D Hs (r ) is equal to Hd (r ). If the ball is replaced by a cube, the cubical contact distribution Hc (r ) is obtained. Reconstruction by simulated annealing In what follows we consider the most important case that the reconstruction is to be done on a cubic lattice. Many other reconstruction schemes are at least briefly discussed in Li (2014). The aim is to obtain a three-dimensional binary pixel structure with summary characteristics VV and Fˆ 1 (r ), . . . , Fˆ s (r ) which are as close as possible to the prescribed F1 (r ), . . . , Fs (r ), the latter coming from the available data. In order to approximate statistical homogeneity usually periodic boundary conditions are assumed. In least-squares thinking the following ‘energy’ E = α1
(F1 (r ) − Fˆ 1 (r ))2 + . . . (r )
+ αs
(Fs (r ) − Fˆ s (r ))2 → min (r )
C 2015 The Authors C 2015 Royal Microscopical Society, 258, 173–178 Journal of Microscopy
(1)
175
is minimized. Usually, for convenience, an optimal pixel structure is searched for only among those structures which have the prescribed fraction VV of 1-pixels. The αi in Eq. (1) are weights, often αi ≡ 1 is used. Of course, in the latter case the resulting energy depends mainly on those Fi which have the highest magnitude. A more equitable influence of the Fi can thus be obtained if the αi are taken inversely proportional to the magnitude of the prescribed Fi . On the other hand, certain Fi might get higher weights αi if they are subjectively most important in the sense that primarily these are intended to be recovered by the reconstructed structure. The summation over r goes over suitable interpixel distances. In general, the employed set of summary characteristics F1 (r ), . . . , Fs (r ) does not uniquely characterize a pixel structure, hence, an optimal pixel structure is not unique. The today standard minimization method, usually preferable even to methods with a higher potential of parallelization, such as genetic algorithms, is simulated annealing. However, in the case that only certain summary characteristics are used also other optimization techniques, such as, e.g. a gradientbased method (Fullwood et al., 2008) might be applied, see Li (2014) and the references therein. In its most basic form simulated annealing works as follows. 1. In the beginning the pixels are randomly marked as 0 (void) and 1 (solid), such that volume fraction is VV . 2. The starting configuration is then stepwise changed, at iteration step k a 0- and a 1-pixel are randomly chosen and their marks are interchanged. The change slightly modifies the functional summary characteristics Fˆ i (r ), which describe the actual pixel structure. Therefore, also the energy E defined by (1) changes from E k to E k+1 , while volume fraction VV is preserved. 3. A pixel change is accepted with probability p k given by the Metropolis rule E k+1 − E k , (2) p k = min 1, exp − Tk i.e. the modification of the pixel image is always accepted if energy is decreased, otherwise it is accepted only with probability p k . 4. This kind of update is continued until a very small energy E min is achieved. In the Metropolis rule (2) Tk is a control parameter, called ‘temperature,. It decreases with increasing k with the goal that a global minimum is achieved as quickly as possible. Several choices of the temperature/annealing schedule are sketched in Deutsch & Cockerham (1994) and Talukdar et al. (2002). Cost savings might be achieved with the simpler form 1 if E k > E k+1 pk = 0 otherwise, of accepting only ‘downhill’ changes (Torquato, 2002).
176
F. BALLANI AND D. STOYAN
The energy threshold E min in the stopping rule of the reconstruction algorithm indicates how good (resp. accurate) the reconstructed pixel structure is. A too small value of E min might, however, lead to very long run times. Jiao et al. (2008) propose to determine E min from the just acceptable number of misplaced 1-pixels and give an explicit rule for the case that F1 (r ) = C (r ) is the only employed summary characteristic. The general idea is to study how the energy changes if in the available pixel structure (the data) one, two or more 1-pixels are, respectively, swapped with 0-pixels. Even if no clear relationship between energy and the number of swapped pixel pairs can be detected, such a preliminary study can give valuable hints for an acceptable range of E min . Practical aspects and modifications As for many optimization methods also for simulated annealing the general statement applies that the ‘parameters used in optimization [...] are more critical than the method itself. [...] If the parameters are chosen inappropriately, then the optimization method will not reach the minimization efficiently [...] Still there is plenty of room for efficiency improvements [...]’ (Li, 2014). A very important point is to calculate the Fˆ i (r ) in a way that uses the fact that in each step the pixel structure is changed only slightly and thus one can reuse the Fˆ i (r ) measured in time step k instead of recomputing them fully at k + 1. This is considered a basic requirement and is discussed in detail in Yeong & Torquato (1998a) for C (r ) and the linealpath function from where the ideas directly carry over to the updates of L (r ) and Hl (r ), and, accordingly, also to Hd (r ) and Hq (r ). An efficient update of Hs (r ) and Hc (r ) is, due to their three-dimensional nature, much more involved but still possible. A further important aspect of efficiency is that of the random pixel-selecting process. In case the media to be reconstructed contain large convex domains of one type of pixels, the basic pixel-selecting process in later stages of the simulated annealing algorithm very often selects pixels inside such domains. However, these trials are mostly rejected since they usually do not decrease energy, and, hence, convergence of the algorithm becomes very slow. Jiao et al. (2008) propose to replace the basic pixel-selecting process in later stages by biasing the selecting process towards pixels on the surface of such clusters, i.e. towards pixels having a small number of equally marked nearest neighbour pixels. Although the minimization scheme of simulated annealing is an inherently sequential algorithm since any decision about accepting or rejecting the proposed modification of the current state influences all subsequent states of the structure, there are several ways to accelerate it by parallelizing. Often, depending on the used characteristics, parallelization during the computation of their updates is possible. If, for instance, a characteristic of a lower-dimensional nature like C (r ), L (r )
or Hd (r ) is used then the computations for its update can be separated with respect to the involved pixel lines or layers and, hence, be parallelized. Numerical experiments show that this makes the computations faster with a factor 3 to 4. Furthermore, already Ouenes & Saad (1993) proposed to compute trials in later stages of the algorithm parallely since then the acceptance rate is very low and nearly all proposed pixel swaps would be rejected. Another idea towards a faster convergence of the simulated annealing algorithm is to start with an initial pixel structure somehow better than the purely random one in step 1. Talukdar et al. (2002) propose to save time by first generating a sample of a thresholded Gaussian random field with the target volume fraction and covariance and then continuing with simulated annealing in order to get good matches with further target summary characteristics, see also Jiang et al. (2013) for a recent application in this journal. However, the effect of this kind of precomputing is surely case-dependent and not well explored in the literature, leaving enough potential for further investigations. For instance, instead of using a thresholded Gaussian random field for the initial guess any other model, e.g. a Boolean (germ-grain) model (Chiu et al., 2013) or, as in Politis et al. (2008), a random sphere packing might take this part over. There are many further modifications of the basic algorithm. Of special interest may be the method in the spirit of germ-grain or particle models as described in Singh et al. (2006). There the grains are empirically given and first the germs are modelled and then the grains are placed. Subsequently, the position and orientation of the grains are optimized. We also mention the method described in Li (2014, p. 449), where ellipsoids are packed. Stereological reconstruction In relation to a classical topic of this journal, we discuss below the special case of stereological reconstruction, i.e. reconstruction based on planar sections. The data come here from a section plane, which is in the following assumed to be a two-dimensional binary pixel image. This pattern is analysed with the statistical methods for planar random sets as in Chiu et al. (2013), which yield area fraction A A (of 1-pixels), C (r ), L (r ), Hl (r ) and Hd (r ) (or Hq (r )). These characteristics can be used without further operations for three-dimensional reconstruction, since VV = A A and the named functional summary characteristics are valid also in three dimensions. A three-dimensional reconstruction should be carried out conditionally in order to make maximal use of the given information, i.e. in the top layer of the three-dimensional pixel structure the given pixel values are fixed during the whole simulation, which is achieved by applying the random pixelselecting process to all layers except the top layer. This conditioning implies that the values of pixels in layers near to the top layer are, as it has to be, spatially correlated to the values C 2015 The Authors C 2015 Royal Microscopical Society, 258, 173–178 Journal of Microscopy
RECONSTRUCTION OF RANDOM HETEROGENEOUS MEDIA
C 2015 The Authors C 2015 Royal Microscopical Society, 258, 173–178 Journal of Microscopy
1.0 0.8 0.6 0.4 0.2 0.0
at their neighbours in the top layer, see Figure 2. Stereological reconstruction first appeared in Hazlett (1997) and is a particular case of what Deutsch & Cockerham (1994) introduce as ‘conditioning to local data’. We demonstrate stereological reconstruction by a small example. As reference structure we consider a three-dimensional sample of size 100 × 100 × 100 pixels which was obtained by digitizing a part of a large random configuration of nonoverlapping balls within a cubic box of side length 1 cm. In stereological thinking the top layer plays the role of the available planar data. Figure 3 shows both the reference structure and structures of the same size obtained by stereological reconstruction using different sets of summary characteristics and keeping always the top layer of the reference structure fixed. The volume fraction of the reference structure is VV = 0.3608 whereas the area fraction of the top layer is A A = 0.3727, hence, the volume fraction of the reconstructed structures was taken to be 0.3727, which deviates only little from that of the reference structure. Obviously, the reconstruction does not lead to digital balls and with increasing distance from the section plane the irregularity increases. Nevertheless, the reconstruction yields realistic estimates of three-dimensional summary characteristics that cannot be estimated stereologically. An example is the cubical contact distribution function Hc (r ). Figure 4 shows that the reconstruction using the quadratic contact
Fig. 3. Reference structure (top left) and stereologically reconstructed structures using C (r ) related to the solid phase (top right), C (r ) and Hq (r ) related to the solid phase (bottom left), C (r ) related to the solid phase and L (r ) related to both solid and void phase (bottom right).
H c (r)
Fig. 2. First three top layers of the stereologically reconstructed structures using C (r ) related to the solid phase (left column), C (r ) and Hq (r ) related to the solid phase (middle column), C (r ) related to the solid phase and L (r ) related to both solid and the void phase (right column). The solid phase is represented by white pixels.
177
0
2
4
6
8
10
12
r [pixel length] Fig. 4. Cubical contact distribution Hc (r ) for the reference structure (black), for the reconstruction using C (r ) related to the solid phase (red), for the reconstruction using C (r ) and Hq (r ) related to the solid phase (green) and for the reconstruction using C (r ) related to the solid phase and L (r ) related to both solid and void phase (blue).
178
F. BALLANI AND D. STOYAN
distribution Hq (r ) reproduces very well Hc (r ) of the reference sample. It is not surprising that just the approach using Hq (r ) is best for the purpose of determining Hc (r ). We tried also to determine the specific Euler number NV . (Note that the determination of the spatial particle number is one of the big problems of stereology.) In our small example the specific Euler number of the reference structure was estimated as 122 cm−3 , which is close to the particle number per unit volume. For all three reconstructions we obtained specific Euler numbers in the order of −200 cm−3 . This result cannot surprise as the reconstruction process generates particles with holes and mutual contacts. When we dilated the reconstructed structures with digital balls of radius one pixel length and then eroded by digital balls of radius six pixel lengths, we obtained positive specific Euler numbers and for the reconstruction on the basis of C (r ) and Hq (r ) a value of 134 cm−3 . Probably much better results can be obtained when applying stochastic reconstruction to serial sections. Then the data from all section planes, which have only little distances, are used in conditioning. Conclusions Stochastic reconstruction is today an established method aimed to generate samples of three-dimensional random heterogeneous media with prescribed properties for further analysis, based on a stochastic optimization scheme such as simulated annealing. For the user it has the advantage that no explicit model assumptions are necessary, they are hidden in the algorithm itself. Although classical models of spatial structures such as Boolean models or thresholded Gaussian random fields start from clear structural assumptions and need mathematical work to obtain their summary characteristics, in a stochastic reconstruction framework some summary characteristics are prescribed and spatial structures are generated in such a way that it has just these given summary characteristics. The only requirement for the employed summary characteristics is that they are available from the given data. In the proper sense of ‘reconstruction’ it is also possible to build-in explicitly the given data in a conditional approach. The method of stochastic reconstruction is computerintensive and needs clever programming. Even with most of the proposed acceleration methods, running times are usually quite long and are growing exponentially with increasing size. It is thus helpful to use the potential of modern parallel computer architectures. Experiences from frequent use of stochastic reconstruction will lead to improvements in the application of the method. It is still a topic of further research to cope with the question which summary characteristics are suitable for
which kind of random heterogeneous media and aimed threedimensional information. Also the extent in which methods of mathematical morphology may be used in improving reconstructed structures will be a topic of future research.
References Chiu, S.N., Stoyan, D., Kendall, W.S. & Mecke, J. (2013) Stochastic Geometry and Its Applications. 3rd edn. Wiley, Chichester. Deutsch, C.V. & Cockerham, P.W. (1994) Practical considerations in the application of simulated annealing to stochastic simulation. Math. Geol. 26(1), 67–82. Fullwood, D.T., Kalidindi, S.R., Niezgoda, S.R., Fast, A. & Hampson, N. (2008) Gradient-based microstructure reconstructions from distributions using fast Fourier transforms. Mater. Sci. Eng. A 494, 68–72. Hazlett, R.D. (1997) Statistical characterization and stochastic modeling of pore networks in relation to fluid flow. Math. Geol. 29, 801–822. Hornby, A.S. (1974) Oxford Advanced Learner’s Dictionary of Current English. Oxford University Press. Oxford. Jiang, Z., Chen, W. & Burkhart, C. (2013) Efficient 3D porous microstructure reconstruction via Gaussian random field and hybrid optimization. J. Microsc. 252(2), 135–148. Jiao, Y., Stillinger, F.H. & Torquato, S. (2008) Modeling heterogeneous materials via two-point correlation functions. II. Algorithmic details and applications. Phys. Rev. E 77, 031135. Jiao, Y., Stillinger, F.H. & Torquato, S. (2009) A superior descriptor of random textures and its predictive capacity. Proc. Nat. Acad. Sci. 106, 17634–17639. Li, D. (2014) Review of structure representation and reconstruction on mesoscale and microscale. JOM 66(3), 444–454. ¨ Ohser, J. & Mucklich, F. (2000) Statistical Analysis of Microstructures in Materials Science. Wiley, Chichester. Ouenes, A. & Saad, N. (1993) A new, fast parallel simulated annealing algorithm for reservoir characterization. SPE paper 26419 presented at the 68th Annual Technical Conference and Exhibition of the Soc. Petroleum Engineers (Houston), pp. 19–29. Politis, M.G., Kikkinides, E.S., Kainourgiakis, M.E. & Stubos, A.K. (2008) A hybrid process-based and stochastic reconstruction method of porous media. Micropor. Mesopor. Mat. 110, 92–99. Singh, H., Gokhale, A.M., Mao, Y. and Spowart, J.E. (2006) Computer simulations of realistic microstructures of discontinuously reinforced aluminum alloy (DRA) composites. Acta Mater. 54, 2131–2143. Talukdar, M.S., Torsaeter, O., Ioannidis, M.A. & Howard, J.J. (2002) Stochastic reconstruction, 3D characterization and network modeling of chalk. J. Petrol. Sci. Eng. 35, 1–21 Torquato, S. (2002) Random Heterogeneous Materials: Microstructure and Macroscopic Properties. Springer, New York. Yeong, C.L.Y. & Torquato, S. (1998a) Reconstructing random media. Phys. Rev. E 57(1), 495–506. Yeong, C.L.Y. & Torquato, S. (1998b) Reconstructing random media. II. Three-dimensional media from two-dimensional cuts. Phys. Rev. E 58(1), 224–233.
C 2015 The Authors C 2015 Royal Microscopical Society, 258, 173–178 Journal of Microscopy
doi: 10.1111/jmi.12234
Received 30 September 2014; accepted 23 January 2015
Reconstruction of random heterogeneous media F. BALLANI & D. STOYAN Institut f¨ur Stochastik, Technische Universit¨at Bergakademie Freiberg, Freiberg, Germany
Key words. Random heterogeneous media, simulated annealing, stochastic reconstruction.
Summary Stochastic reconstruction is a technique to generate samples of random structures with prescribed distributional properties in the sense that certain of their statistical summary characteristics match target values or forms. This technique can be used to produce structures of any wanted size for further statistical analysis starting from small samples, which may be even only lower dimensional, for instance, when three-dimensional imaging techniques are not available. In this review we explain the main ideas of stochastic reconstruction, concentrating on the most important case of digitized binary media and with a particular emphasis on stereological reconstruction.
Introduction Stochastic reconstruction is a technique to generate samples of random structures, typically in the sense that the obtained samples, in some respects, are very similar to a given real sample. This includes the generation of similar and possibly larger planar patterns of a random heterogeneous medium starting from a planar pattern, likewise that of a spatial pattern of a random heterogeneous medium from a spatial or planar pattern, as well as the reconstruction of point patterns. In this review, we explain the main ideas of this technique, where we restrict to the special case that samples of a three-dimensional medium consisting of two phases (implying its heterogeneity) have to be reconstructed. Throughout the paper we assume that the random medium is statistically homogeneous (spatially stationary) and isotropic (see, e.g. Torquato, 2002). We expect some knowledge of fundamental ideas of spatial statistics as, e.g. presented in Chiu et al. (2013). In the analysis and characterization of random heterogeneous media it is a fundamental task to make inference on an appropriate three-dimensional representation of a real-world specimen in a computer (typically a pixel image), for instance to get information on topological properties and on characteristics of connectivity, flow and transport properties, etc. Such ¨ Stochastik, Technische Universit¨at Correspondence to: F. Ballani, Institut fur Bergakademie Freiberg, D-09596 Freiberg, Germany. Tel: +49 (0)3731 39-2796; fax: +49 (0)3731 39-3598; e-mail: [email protected] C 2015 The Authors C 2015 Royal Microscopical Society Journal of Microscopy
virtual structures must be large enough to obtain reliable estimates of complicated characteristics. However, very often they are not directly available: samples obtained by modern imaging methods such as X-ray computed tomography may be too small and data from planar sections like scanning electron microscope images, which are often the only alternative in case three-dimensional imaging techniques are not available, are by their nature not sufficient to estimate spatial topological characteristics. Stochastic reconstruction is able to produce the three-dimensional samples needed, starting from (small) samples, which may be even only lower dimensional, resulting from planar (or even linear) sections. It can yield many samples of arbitrary size in full dimensionality. The classical approach to generate samples of threedimensional sets as described in Chiu et al. (2013) consists in fitting a model (like a Boolean model or a thresholded Gaussian random field) to the data, where the model parameters come somehow from the data, typically from estimated summary characteristics. Then the model is simulated according to its properties. For the nonspecialist this approach has the disadvantage to use complicated statistical techniques and the need to believe in a model she/he does perhaps not understand. Modern stochastic reconstruction is more general and (presumably) user-friendly since it does not use explicit models but only models implicitly determined by the algorithm. The user has only to understand the general idea and some functional summary characteristics we explain in the next section. The general idea is simple: construct a three-dimensional structure which has summary characteristics that match exactly (or in very good approximation) the prescribed ones, which come from statistics of the originally given sample. It is even possible to build-in explicitly the given data in a conditional approach. Though in every simulation of stochastic models given statistical data there is a reconstructive aspect (in the sense that ‘reconstruction means to build up a complete structure of which one has only a few parts or only partial evidence’, see Hornby, 1974) the term ‘stochastic reconstruction’ is used for this special method – maybe a clever marketing brand. According to Torquato (2002, section 8.2) an origin of stochastic reconstruction may be seen in ideas of statistics of Gaussian random fields. However, the modern reconstruction
174
F. BALLANI AND D. STOYAN
approach does not use such models and can be understood without deeper knowledge of probability theory. This modern approach goes back to a paper by Farmer in 1992, the Ph.D. thesis of Deutsch in 1992 and several subsequent publications, e.g. Deutsch & Cockerham (1994) (see the references in Hazlett, 1997), and two papers of Yeong & Torquato (1998a,b). Since the year 1992 there have been published many papers on various methods of reconstruction, about variants of the classical procedure, about the use of various summary characteristics and directed to various structural properties of the random media to reconstruct. A lot of references can be found in Torquato (2002), Chiu et al. (2013) and Li (2014). In the following, we first explain some spatial summary characteristics which are commonly used in stochastic reconstruction of two-phase structures. Then we sketch the reconstruction algorithm for three-dimensional statistically homogeneous (stationary) and isotropic two-phase structure in case it is represented as a binary pixel image, where 0-pixels stand for pore, and 1-pixels for solid (or for two different phases). Since stochastic reconstruction is a computerintensive method we also list some practical aspects which might be helpful to implement it in a fairly fast version. Finally, we discuss in some detail the stereological problem where a planar section is given and based on that a spatial sample has to be reconstructed, and consider even an example, which yields interesting insights. Summary characteristics The structure to be reconstructed is characterized by various summary characteristics. These are estimated from the original sample, based on scattering methods, X-ray computed tomography or microscopy. The choice of the summary characteristics depends on the available data, on the nature of the information aimed to get by the reconstruction, and on the associated computational costs. The main characteristic of, say, the solid phase Z of a statistically homogeneous structure is volume fraction VV (or porosity 1 − VV ), which is usually a priori matched exactly since only structures with the prescribed volume fraction VV are admitted. The further characteristics are usually of a functional nature. Most popular and important is the covariance C (r ) (or twopoint probability function) of the solid phase Z , which is the probability that two points a distance r apart are both covered by Z . (Note that in a two-phase structure the covariance of one phase unambiguously determines that of the other phase: the covariance of the complementary phase of Z equals 1 − 2C (0) + C (r ).) Using C (r ) in stochastic reconstruction is usually considered a must since it characterizes the second-order behaviour of a structure, in particular, besides VV = C (0), it also determines the specific surface area SV since its slope at the origin is related to SV via SV = −4C (0) (Chiu et al., 2013, p. 219). The function C (r ) has also the advantage
C(r)
Z
r Hd(r)
r L(r) r
Hl(r) r Fig. 1. Planar section together with a schematic illustration of events contributing to summary characteristics available from planar sections, for covariance C (r ), chord length distribution function L (r ) and linear and discoidal contact distribution function Hl (r ) and Hd (r ), see the text for further details.
that it may be determined even from small-angle scattering experiments. Further functional characteristics (that are often incorrectly called ‘correlation functions’) are, e.g. the chord-length distribution function, various contact distribution (or, equivalently, empty space) functions and the two-point cluster function, see Chiu et al. (2013) for a detailed introduction. The chord length distribution function L (r ) (not to be confused with the popular L -function derived from Ripley’s K function for random point patterns, see Chiu et al., 2013) related to the solid phase Z is the portion of Z -chords with length less or equal than r , where a Z -chord is a segment between intersection points of a line with the two-phase interface fully contained in Z , see also Figure 1. Hence, L (r ) contains information on connectedness and correlation along a lineal path. The same information is carried by the linear contact distribution function Hl (r ) of the complement of Z as well as by the lineal path function (see Torquato, 2002, p. 44). Particularly apparent for L (r ), Figure 1 suggests that the estimation of L (r ) (but also that of any other summary characteristic depending on some distance r ) might be heavily influenced by edge effects. Often a corrected estimation accounting for such edge ¨ effects is meaningful, see Ohser & Mucklich (2000, p. 102) for L (r ) and Chiu et al. (2013) for general ideas and further references. It should be noted that all summary characteristics considered so far are of a one-dimensional nature. This has C 2015 The Authors C 2015 Royal Microscopical Society, 258, 173–178 Journal of Microscopy
RECONSTRUCTION OF RANDOM HETEROGENEOUS MEDIA
the advantage that they can be determined even from onedimensional samples, and the disadvantage that they represent only simpler spatial properties of the structure. Nevertheless, Li (2014) states that the ‘lineal’ summary characteristics are percolation-sensitive descriptors. As a summary characteristic of a two-dimensional nature, but still available from planar sections, is the discoidal contact distribution function Hd (r ). The function Hd (r ) for the solid phase Z is the probability that a random two-dimensional disc of radius r with centre outside Z hits Z . Thus, Hd (r ) characterizes the planar expansion of regions in the void phase. An analogous function is the quadratic contact distribution function Hq (r ), where the disc is replaced by a square. Finally, if three-dimensional data are available, also summary characteristics with a truly three-dimensional nature might be used. One numerical characteristic of this type is the specific Euler number NV , carrying information on connectivity. Sometimes three-point correlation functions are used for stochastic reconstruction. A further functional summary characteristic leading to good reconstruction results (see Jiao et al., 2009) is the two-point cluster function (denoted by C 2 (r ) in Torquato, 2002, and by C 1 (r ) in Chiu et al., 2013). For the solid phase Z this is the probability that two points a distance r apart are both covered by Z and belong to the same cluster/connectivity component, i.e. there is a connecting path fully covered by Z . With the exception of the case of particle structures with well-separated particles, the determination of C 2 (r ) is rather difficult and time consuming. Simpler three-dimensional summary characteristics are contact distributions with three-dimensional structuring elements. The best known example is the spherical contact distribution function Hs (r ). The function Hs (r ) for the solid phase Z is the probability that a random ball of radius r with centre outside Z hits Z . Hence, in 2D Hs (r ) is equal to Hd (r ). If the ball is replaced by a cube, the cubical contact distribution Hc (r ) is obtained. Reconstruction by simulated annealing In what follows we consider the most important case that the reconstruction is to be done on a cubic lattice. Many other reconstruction schemes are at least briefly discussed in Li (2014). The aim is to obtain a three-dimensional binary pixel structure with summary characteristics VV and Fˆ 1 (r ), . . . , Fˆ s (r ) which are as close as possible to the prescribed F1 (r ), . . . , Fs (r ), the latter coming from the available data. In order to approximate statistical homogeneity usually periodic boundary conditions are assumed. In least-squares thinking the following ‘energy’ E = α1
(F1 (r ) − Fˆ 1 (r ))2 + . . . (r )
+ αs
(Fs (r ) − Fˆ s (r ))2 → min (r )
C 2015 The Authors C 2015 Royal Microscopical Society, 258, 173–178 Journal of Microscopy
(1)
175
is minimized. Usually, for convenience, an optimal pixel structure is searched for only among those structures which have the prescribed fraction VV of 1-pixels. The αi in Eq. (1) are weights, often αi ≡ 1 is used. Of course, in the latter case the resulting energy depends mainly on those Fi which have the highest magnitude. A more equitable influence of the Fi can thus be obtained if the αi are taken inversely proportional to the magnitude of the prescribed Fi . On the other hand, certain Fi might get higher weights αi if they are subjectively most important in the sense that primarily these are intended to be recovered by the reconstructed structure. The summation over r goes over suitable interpixel distances. In general, the employed set of summary characteristics F1 (r ), . . . , Fs (r ) does not uniquely characterize a pixel structure, hence, an optimal pixel structure is not unique. The today standard minimization method, usually preferable even to methods with a higher potential of parallelization, such as genetic algorithms, is simulated annealing. However, in the case that only certain summary characteristics are used also other optimization techniques, such as, e.g. a gradientbased method (Fullwood et al., 2008) might be applied, see Li (2014) and the references therein. In its most basic form simulated annealing works as follows. 1. In the beginning the pixels are randomly marked as 0 (void) and 1 (solid), such that volume fraction is VV . 2. The starting configuration is then stepwise changed, at iteration step k a 0- and a 1-pixel are randomly chosen and their marks are interchanged. The change slightly modifies the functional summary characteristics Fˆ i (r ), which describe the actual pixel structure. Therefore, also the energy E defined by (1) changes from E k to E k+1 , while volume fraction VV is preserved. 3. A pixel change is accepted with probability p k given by the Metropolis rule E k+1 − E k , (2) p k = min 1, exp − Tk i.e. the modification of the pixel image is always accepted if energy is decreased, otherwise it is accepted only with probability p k . 4. This kind of update is continued until a very small energy E min is achieved. In the Metropolis rule (2) Tk is a control parameter, called ‘temperature,. It decreases with increasing k with the goal that a global minimum is achieved as quickly as possible. Several choices of the temperature/annealing schedule are sketched in Deutsch & Cockerham (1994) and Talukdar et al. (2002). Cost savings might be achieved with the simpler form 1 if E k > E k+1 pk = 0 otherwise, of accepting only ‘downhill’ changes (Torquato, 2002).
176
F. BALLANI AND D. STOYAN
The energy threshold E min in the stopping rule of the reconstruction algorithm indicates how good (resp. accurate) the reconstructed pixel structure is. A too small value of E min might, however, lead to very long run times. Jiao et al. (2008) propose to determine E min from the just acceptable number of misplaced 1-pixels and give an explicit rule for the case that F1 (r ) = C (r ) is the only employed summary characteristic. The general idea is to study how the energy changes if in the available pixel structure (the data) one, two or more 1-pixels are, respectively, swapped with 0-pixels. Even if no clear relationship between energy and the number of swapped pixel pairs can be detected, such a preliminary study can give valuable hints for an acceptable range of E min . Practical aspects and modifications As for many optimization methods also for simulated annealing the general statement applies that the ‘parameters used in optimization [...] are more critical than the method itself. [...] If the parameters are chosen inappropriately, then the optimization method will not reach the minimization efficiently [...] Still there is plenty of room for efficiency improvements [...]’ (Li, 2014). A very important point is to calculate the Fˆ i (r ) in a way that uses the fact that in each step the pixel structure is changed only slightly and thus one can reuse the Fˆ i (r ) measured in time step k instead of recomputing them fully at k + 1. This is considered a basic requirement and is discussed in detail in Yeong & Torquato (1998a) for C (r ) and the linealpath function from where the ideas directly carry over to the updates of L (r ) and Hl (r ), and, accordingly, also to Hd (r ) and Hq (r ). An efficient update of Hs (r ) and Hc (r ) is, due to their three-dimensional nature, much more involved but still possible. A further important aspect of efficiency is that of the random pixel-selecting process. In case the media to be reconstructed contain large convex domains of one type of pixels, the basic pixel-selecting process in later stages of the simulated annealing algorithm very often selects pixels inside such domains. However, these trials are mostly rejected since they usually do not decrease energy, and, hence, convergence of the algorithm becomes very slow. Jiao et al. (2008) propose to replace the basic pixel-selecting process in later stages by biasing the selecting process towards pixels on the surface of such clusters, i.e. towards pixels having a small number of equally marked nearest neighbour pixels. Although the minimization scheme of simulated annealing is an inherently sequential algorithm since any decision about accepting or rejecting the proposed modification of the current state influences all subsequent states of the structure, there are several ways to accelerate it by parallelizing. Often, depending on the used characteristics, parallelization during the computation of their updates is possible. If, for instance, a characteristic of a lower-dimensional nature like C (r ), L (r )
or Hd (r ) is used then the computations for its update can be separated with respect to the involved pixel lines or layers and, hence, be parallelized. Numerical experiments show that this makes the computations faster with a factor 3 to 4. Furthermore, already Ouenes & Saad (1993) proposed to compute trials in later stages of the algorithm parallely since then the acceptance rate is very low and nearly all proposed pixel swaps would be rejected. Another idea towards a faster convergence of the simulated annealing algorithm is to start with an initial pixel structure somehow better than the purely random one in step 1. Talukdar et al. (2002) propose to save time by first generating a sample of a thresholded Gaussian random field with the target volume fraction and covariance and then continuing with simulated annealing in order to get good matches with further target summary characteristics, see also Jiang et al. (2013) for a recent application in this journal. However, the effect of this kind of precomputing is surely case-dependent and not well explored in the literature, leaving enough potential for further investigations. For instance, instead of using a thresholded Gaussian random field for the initial guess any other model, e.g. a Boolean (germ-grain) model (Chiu et al., 2013) or, as in Politis et al. (2008), a random sphere packing might take this part over. There are many further modifications of the basic algorithm. Of special interest may be the method in the spirit of germ-grain or particle models as described in Singh et al. (2006). There the grains are empirically given and first the germs are modelled and then the grains are placed. Subsequently, the position and orientation of the grains are optimized. We also mention the method described in Li (2014, p. 449), where ellipsoids are packed. Stereological reconstruction In relation to a classical topic of this journal, we discuss below the special case of stereological reconstruction, i.e. reconstruction based on planar sections. The data come here from a section plane, which is in the following assumed to be a two-dimensional binary pixel image. This pattern is analysed with the statistical methods for planar random sets as in Chiu et al. (2013), which yield area fraction A A (of 1-pixels), C (r ), L (r ), Hl (r ) and Hd (r ) (or Hq (r )). These characteristics can be used without further operations for three-dimensional reconstruction, since VV = A A and the named functional summary characteristics are valid also in three dimensions. A three-dimensional reconstruction should be carried out conditionally in order to make maximal use of the given information, i.e. in the top layer of the three-dimensional pixel structure the given pixel values are fixed during the whole simulation, which is achieved by applying the random pixelselecting process to all layers except the top layer. This conditioning implies that the values of pixels in layers near to the top layer are, as it has to be, spatially correlated to the values C 2015 The Authors C 2015 Royal Microscopical Society, 258, 173–178 Journal of Microscopy
RECONSTRUCTION OF RANDOM HETEROGENEOUS MEDIA
C 2015 The Authors C 2015 Royal Microscopical Society, 258, 173–178 Journal of Microscopy
1.0 0.8 0.6 0.4 0.2 0.0
at their neighbours in the top layer, see Figure 2. Stereological reconstruction first appeared in Hazlett (1997) and is a particular case of what Deutsch & Cockerham (1994) introduce as ‘conditioning to local data’. We demonstrate stereological reconstruction by a small example. As reference structure we consider a three-dimensional sample of size 100 × 100 × 100 pixels which was obtained by digitizing a part of a large random configuration of nonoverlapping balls within a cubic box of side length 1 cm. In stereological thinking the top layer plays the role of the available planar data. Figure 3 shows both the reference structure and structures of the same size obtained by stereological reconstruction using different sets of summary characteristics and keeping always the top layer of the reference structure fixed. The volume fraction of the reference structure is VV = 0.3608 whereas the area fraction of the top layer is A A = 0.3727, hence, the volume fraction of the reconstructed structures was taken to be 0.3727, which deviates only little from that of the reference structure. Obviously, the reconstruction does not lead to digital balls and with increasing distance from the section plane the irregularity increases. Nevertheless, the reconstruction yields realistic estimates of three-dimensional summary characteristics that cannot be estimated stereologically. An example is the cubical contact distribution function Hc (r ). Figure 4 shows that the reconstruction using the quadratic contact
Fig. 3. Reference structure (top left) and stereologically reconstructed structures using C (r ) related to the solid phase (top right), C (r ) and Hq (r ) related to the solid phase (bottom left), C (r ) related to the solid phase and L (r ) related to both solid and void phase (bottom right).
H c (r)
Fig. 2. First three top layers of the stereologically reconstructed structures using C (r ) related to the solid phase (left column), C (r ) and Hq (r ) related to the solid phase (middle column), C (r ) related to the solid phase and L (r ) related to both solid and the void phase (right column). The solid phase is represented by white pixels.
177
0
2
4
6
8
10
12
r [pixel length] Fig. 4. Cubical contact distribution Hc (r ) for the reference structure (black), for the reconstruction using C (r ) related to the solid phase (red), for the reconstruction using C (r ) and Hq (r ) related to the solid phase (green) and for the reconstruction using C (r ) related to the solid phase and L (r ) related to both solid and void phase (blue).
178
F. BALLANI AND D. STOYAN
distribution Hq (r ) reproduces very well Hc (r ) of the reference sample. It is not surprising that just the approach using Hq (r ) is best for the purpose of determining Hc (r ). We tried also to determine the specific Euler number NV . (Note that the determination of the spatial particle number is one of the big problems of stereology.) In our small example the specific Euler number of the reference structure was estimated as 122 cm−3 , which is close to the particle number per unit volume. For all three reconstructions we obtained specific Euler numbers in the order of −200 cm−3 . This result cannot surprise as the reconstruction process generates particles with holes and mutual contacts. When we dilated the reconstructed structures with digital balls of radius one pixel length and then eroded by digital balls of radius six pixel lengths, we obtained positive specific Euler numbers and for the reconstruction on the basis of C (r ) and Hq (r ) a value of 134 cm−3 . Probably much better results can be obtained when applying stochastic reconstruction to serial sections. Then the data from all section planes, which have only little distances, are used in conditioning. Conclusions Stochastic reconstruction is today an established method aimed to generate samples of three-dimensional random heterogeneous media with prescribed properties for further analysis, based on a stochastic optimization scheme such as simulated annealing. For the user it has the advantage that no explicit model assumptions are necessary, they are hidden in the algorithm itself. Although classical models of spatial structures such as Boolean models or thresholded Gaussian random fields start from clear structural assumptions and need mathematical work to obtain their summary characteristics, in a stochastic reconstruction framework some summary characteristics are prescribed and spatial structures are generated in such a way that it has just these given summary characteristics. The only requirement for the employed summary characteristics is that they are available from the given data. In the proper sense of ‘reconstruction’ it is also possible to build-in explicitly the given data in a conditional approach. The method of stochastic reconstruction is computerintensive and needs clever programming. Even with most of the proposed acceleration methods, running times are usually quite long and are growing exponentially with increasing size. It is thus helpful to use the potential of modern parallel computer architectures. Experiences from frequent use of stochastic reconstruction will lead to improvements in the application of the method. It is still a topic of further research to cope with the question which summary characteristics are suitable for
which kind of random heterogeneous media and aimed threedimensional information. Also the extent in which methods of mathematical morphology may be used in improving reconstructed structures will be a topic of future research.
References Chiu, S.N., Stoyan, D., Kendall, W.S. & Mecke, J. (2013) Stochastic Geometry and Its Applications. 3rd edn. Wiley, Chichester. Deutsch, C.V. & Cockerham, P.W. (1994) Practical considerations in the application of simulated annealing to stochastic simulation. Math. Geol. 26(1), 67–82. Fullwood, D.T., Kalidindi, S.R., Niezgoda, S.R., Fast, A. & Hampson, N. (2008) Gradient-based microstructure reconstructions from distributions using fast Fourier transforms. Mater. Sci. Eng. A 494, 68–72. Hazlett, R.D. (1997) Statistical characterization and stochastic modeling of pore networks in relation to fluid flow. Math. Geol. 29, 801–822. Hornby, A.S. (1974) Oxford Advanced Learner’s Dictionary of Current English. Oxford University Press. Oxford. Jiang, Z., Chen, W. & Burkhart, C. (2013) Efficient 3D porous microstructure reconstruction via Gaussian random field and hybrid optimization. J. Microsc. 252(2), 135–148. Jiao, Y., Stillinger, F.H. & Torquato, S. (2008) Modeling heterogeneous materials via two-point correlation functions. II. Algorithmic details and applications. Phys. Rev. E 77, 031135. Jiao, Y., Stillinger, F.H. & Torquato, S. (2009) A superior descriptor of random textures and its predictive capacity. Proc. Nat. Acad. Sci. 106, 17634–17639. Li, D. (2014) Review of structure representation and reconstruction on mesoscale and microscale. JOM 66(3), 444–454. ¨ Ohser, J. & Mucklich, F. (2000) Statistical Analysis of Microstructures in Materials Science. Wiley, Chichester. Ouenes, A. & Saad, N. (1993) A new, fast parallel simulated annealing algorithm for reservoir characterization. SPE paper 26419 presented at the 68th Annual Technical Conference and Exhibition of the Soc. Petroleum Engineers (Houston), pp. 19–29. Politis, M.G., Kikkinides, E.S., Kainourgiakis, M.E. & Stubos, A.K. (2008) A hybrid process-based and stochastic reconstruction method of porous media. Micropor. Mesopor. Mat. 110, 92–99. Singh, H., Gokhale, A.M., Mao, Y. and Spowart, J.E. (2006) Computer simulations of realistic microstructures of discontinuously reinforced aluminum alloy (DRA) composites. Acta Mater. 54, 2131–2143. Talukdar, M.S., Torsaeter, O., Ioannidis, M.A. & Howard, J.J. (2002) Stochastic reconstruction, 3D characterization and network modeling of chalk. J. Petrol. Sci. Eng. 35, 1–21 Torquato, S. (2002) Random Heterogeneous Materials: Microstructure and Macroscopic Properties. Springer, New York. Yeong, C.L.Y. & Torquato, S. (1998a) Reconstructing random media. Phys. Rev. E 57(1), 495–506. Yeong, C.L.Y. & Torquato, S. (1998b) Reconstructing random media. II. Three-dimensional media from two-dimensional cuts. Phys. Rev. E 58(1), 224–233.
C 2015 The Authors C 2015 Royal Microscopical Society, 258, 173–178 Journal of Microscopy
