Adaptive phase-retrieval stochastic reconstruction with correlation functions: 3D images from 2D cuts

Precise characterization of three-dimensional heterogeneous media is indispensable in finding the relationships between structure and macroscopic physical properties (permeability, conductivity, and others). The most widely used experimental methods (electronic and optical microscopy) provide high-resolution bi-dimensional images of the samples of interest. However, 3D material inner microstructure registration is needed to apply numerous modeling tools. Numerous research areas search for cheap and robust methods to obtain"full"3D information about the structure of the studied sample from its 2D cuts. In this work, we develop a dynamic phase-retrieval stochastic reconstruction algorithm that can create 3D replicas from 2D original images - DDTF. The DDTF is free of artifacts characteristic of previously proposed phase-retrieval techniques. While based on a two-point $S_2$ correlation function, any correlation function or other morphological metrics can be accounted for during the reconstruction, thus, paving the way to the hybridization of different reconstruction techniques. In this work, we use two-point probability and surface-surface functions for optimization. To test DDTF, we performed reconstructions for three binary porous media samples of different genesis: sandstone, carbonate, and ceramic. Based on computed permeability and connectivity ($C_2$ and $L_2$ correlation functions), we have shown that the proposed technique in terms of accuracy is comparable to the classic simulated annealing-based reconstruction method but is computationally very effective. Our findings open the possibility of utilizing DDTF to produce fast or crude replicas further polished by other reconstruction techniques such as simulated annealing or process-based methods.


I. INTRODUCTION
Precise characterization of three-dimensional heterogeneous media is indispensable in finding the relationships between structure and macroscopic physical properties (permeability, conductivity, and other properties). [1][2][3].The most widely used experimental methods (electronic and optical microscopy) provide high-resolution bidimensional images of the samples of interest. However, 3D material inner microstructure registration is needed to apply modeling tools. Obtaining quality 3D information is technically more complicated and requires expensive equipment and qualified operators. Thus, porous media reconstruction from thin sections [4,5] is actual and has many petroleum and food engineering [6,7], material and soil sciences, [8,9], electrochemistry [10], and * cherkasov.am@phystech.edu medical applications [11]. Effective reconstruction procedures allow the creation of high-quality digital models for further analysis and calculation of physical properties. This is true for various multiphase materials, including porous media with a binary structure consisting of pores and solids. A versatile set of methods is available in the literature for stochastic reconstructions, including creating 3D replicas from 2D slices. Truncated Gaussian random fields [12,13] are a relatively computationally inexpensive approach but do not always allow the preservation of morphological properties of porous media. Simulated annealing (SA) using correlation functions [14,15] allows obtaining high-quality reconstructions, but it is a long-duration process even when accelerated (e.g., hierarchical simulated annealing) [16,17]). Note that hierarchical annealing schemes with correlation function rescaling were proposed recently to overcome this computational disadvantage [18]. Multiple-point statistics (MPS) [6,19,20] reconstruction was shown to reproduce mor-phology quite well compared to other methods. However, it is computationally expensive and may produce repetitive textures based on the number of points used for statistics. Quite surprisingly, some MPS methods are even slower and less accurate than SA based on twopoint statistics [21]. Deep learning approaches [22][23][24] are getting popular and show great promise, but even more computationally expensive than MPS and SA, with lower accuracies not balanced by massive training times. Process-based methods are very efficient and accurate, but their applicability is limited to granular porous media [25,26]. Phase-retrieval [27] is a fast algorithm that allows reconstructing volume from its two-point statistics [3] (or two-point probability correlation function). A method for approximating three-dimensional two-point statistics by a set of two-dimensional ones has been developed by Hasanabadi et al [28], and a three-dimensional image has been reconstructed from a bi-dimensional image. This approach has high performance and retains the basic microstructural properties of the media. However, the resulting image is rather noisy: artifacts like stripes and individual misplaced voxels appeared in the reconstructed sample. While this had only a minor impact on simulated electrical conductivity, other physical properties, especially permeability, are very sensitive to such noise [29]. The originally described phase-retrieval method [27] used two-point statistics between each possible combination of pairs of pixels/voxels. This way, it was possible to reconstruct the images exactly [30]. However, in this case, the amount of two-point statistics is larger than the 2D or 3D image itself (even if we discard the parts redundant due to symmetry). Such reconstructions do not seem to have much practical application, especially considering the fact that two-point statistics that can be measured experimentally, for example, with the help of small-angle scattering [31] or X-ray tomography [32], is limited to ensemble average over the whole volume or some direction. Ensemble statistics is a kind of compressed two-point statistics and can be effectively utilized to reduce the structural information about the object at hand [33][34][35]. On the other hand, if the aim is to reconstruct 3D information from the 2D slice, using the complete information available from such a slice using full two-point statistics may be beneficial. In this work, we develop a phase-retrieval algorithm of three-dimensional sample reconstruction from its bidimensional cut-section free of the abovementioned disadvantages. The paper is organized as follows. We review the mathematics behind the phase-retrieval algorithm as well as the algorithm itself. Next, we describe rotation as a method for increasing the dimensionality of self-convolution using the fullest possible two-point statistics read from a 2D image. A new phase retrieval algorithm based on additional constraints in Fourier space and static and dynamic adjustments of three-dimensional self-convolution is presented. To demonstrate achieved improvements, we simulate single-phase (or saturated) flow on both original and reconstructed samples of three porous media sam-ples of different genesis. We also evaluate the two-point cluster correlation function known to provide non-trivial connectivity information.
Microstructure m n s is a binary multidimensional array where values 0 and 1 correspond to void and solid voxels, respectively. (s 1 , ..., s N ) : s 1 ∈ (1, 2, ..., S 1 ), s N ∈ (1, 2, ...S N ) is a vector of space coordinates while n enumerates different components present in certain images, and N is the dimension of space to which microstructure appertains. Term correlation is commonly used in signal processing for convolution of the image with inversed self. However, this term in statistics has another meaning. Thus in this article, correlation in the sense of image processing is called convolution, while correlation is obtained as a result of either averaging convolutions by different directions in one image or by image ensemble averaging. Convolution of m n s is defined the following way: where S = (S 1 , ..., S N ), t = (t 1 , ..., t N ) ∈ R n . It contains sufficient data for retrieval of the initial binary image to within a translation and inversion. Fulwood et al. establish the link between Fourier transform of microstructure and self-convolution so that the Gershberg-Saxton algorithm of image retrieval from magnitude of its Fourier transform becomes applicable.
f nn s e 2πisk/(S1·...·Sn) The resulting algorithm 1 consists of main stages. Fourier transform of microstructure on current iteration.

5:
Replacement its magnitude with the magnitude obtained from self-convolution.

6:
Inverse Fourier transform of replacement result.
Satisfying constraints in real space: Using black padding around the input image before calculating its self-convolution and taking into account the position of padding during recovery results in a better convergence and reduction of translation uncertainty [36]. Nevertheless, in the case of approximated convolution reconstructed image with padding tends to have lower porosity in the slices adjacent to padding. Thus usage of padding is appropriate in cases when the role of boundary effects does not seriously influence macroscopic properties like porosity and permeability.
B. The transition from 2d convolution to 3d convolution via rotation One way to obtain a one-dimensional correlation for isotropic microstructure is by averaging its selfconvolution in all directions. Self-convolution has axial symmetry. We propose rotation as a way to translate convolution from bi-dimensional to three-dimensional space. Appropriate transition from 2d to 3d convolution is rather important, as it results in quality of reconstruction. In representative volume of porous media there is linear relationship between expected values of 3d-convolution and 2d-convolution for equal shift values. For binary im- is where p is intersection porosity of original and shifted image. For 3d-convolution for shift (t, 0, 0) expectation is In case of adding black padding and doubling the linear size of an image S = 2S The approximation with properties mentioned above can be used as input for the phase-recovery algorithm. Sample retrieved from this 3d-convolution has a preferential direction which is the convolution rotation axis. Rotation of 2d-convolution around ox, oy, oz axis and averaging the results lead to isotropy of reconstructed image (2). Three-dimensional convolution does not always correspond to existent microstructure, resulting in defects in a reconstructed image such as deformation and noise (figure 1). Adding constraints in Fourier space such as frequency filtering is an efficient way of noise reduction. Further in the article, samples are presented without padding for better presentation and visibility. Though filter adjustment should be performed in each specific case, the width of the corresponding filter in real space should be no less than one voxel. Figure 3 shows the reconstructed samples. The first image is obtained with no filtering. The following images correspond to space filters' widths of 0.5, 1, 1,5, and 2 voxels. Choice of appropriate mask width is the solution to the optimization problem for the cost function. For example, it is possible to minimize the difference between the correlation of the input image and the median correlation calculated of all slices. Stochastic optimization is warranted here because of a significant number of local minima. Thus modified algorithm 2 at step 6 also contains multiplication. Fourier transform of microstructure on current iteration.
Replacement its magnitude with the magnitude obtained from self-convolution.

6:
Frequency filtration with Gaussian filter 7: Satisfying constraints in real space:  Figure 4 shows that transition functions from 2dcorrelations to 3d-correlation for such objects as square, circular disc, and 2d-plurigaussian are non-linear and sufficiently different. We refer to such transition functions for correlations as "Dimensional Transition Functions" or DTF.

Static dimensional transition functions
Thus approximated 3d-convolution should be adjusted later on account of microstructure morphology. Plurigaussian modeling of microstructures is one of the widely used methods [12] for microstructure modeling. Random topography field T is obtained by convolving kernel k(x) with u(x) ∈ U -random field with uniform probability density distribution.
The convolution of the Gaussian kernel with the random field is shown in figure 5.
We implemented an algorithm 3 to generate an approximated microstructure in which 2d-slice corresponds with the input image. Its realization is described thereunder. Calculate the self-convolution of a 2d-image. 2: Average self-convolution by all directions to obtain 1D-correlation. 3: Translate 1D-correlation to 3d-space by rotation. 4: Generate a 3d-random field with uniform distribution. 5: Convolve kernel obtained by rotation of self-convolution and random field with uniform distribution. 6: Truncate the resulting plurigaussian so that the porosity of the middle slice of the resulting microstructure is the closest to the porosity of the input image. 7: Circularly shift truncated plurigaussian for alignment with input slices using a maximum of cross-convolution. 8: Apply morphological erosion or dilation to minimize differences between correlation functions of the input slice and the particular slice of truncated plurigaussian. 9: Calculate correlation function from radial shift for certain 2D-slice of plurigaussian approximation. 10: Calculate correlation function from radial shift for the whole plurigaussian approximation. 11: Plot set of points. Each point corresponds to a particular radial shift where the x coordinate is the value for 2d-correlation, and the y-value is 3d-correlation. 12: Approximate the dependency by power-law to apply it to 3d-autoconvolution obtained by rotation from bidimensional one Such plurigaussian microstructure is used only for building dimensional transition function in each specific case. A synthetic case like circular disk (figure 4) is an example when appropriate convolution transform is essential. Using traditional scaling, we obtain significantly worse results while plurigaussian approximation for circular disk leads to the reconstruction of an object similar to a solid sphere.
In the case of small-scale microstructure, DTF-correction does not play a significant role because the transition function is almost linear. Properties of the resulting sample differ in relation to noise realization because the kernel is convolved with random noise. Thus it is necessary to choose the best realization according to quality metrics such as correlation function. However, in the case of images whose characteristic scale is comparable with its size, plurigaussian approximation does not always give a reasonable guess for dimensional transitional coefficients.

The dynamic dimensional transition function
Quality for building both convolution and DTF influences reconstruction results. The solution is on-the-fly updating of the convolution in such a way that allows adjusting parameters of convolution during phase-retrieval algorithm execution. For example, parameters can be chosen to minimize difference according to a particular metric between a specific slice or slices of reconstructed and input image. The choice of metric is a single issue and depends on the integral characteristic of the 3dimage ( [3]). Approximated 3d-autoconvolution has radial and angular parts. The radial part is relatively smooth in all points except (0, 0, 0), while the angular one somewhat oscillates around zero field with magnitude decreasing with the growth of distance from autocorrelation center (0, 0, 0). Coefficients in the linear combination of these two parts in the resulting approximation for autoconvolution are also optimization parameters. Thus modified algorithm 4 is the following.
where a and b are dimensional transition coefficients.
where rad and ang are radial and angular parts respectively, and c is optimization coefficient. 3: Making random binary noise initial guess for recovered microstructure (m 1 s ) 0 . 4: repeat

5:
Fourier transform of microstructure on current iteration.
Replacement of its magnitude with the magnitude obtained from self-convolution.

7:
Frequency filtration with Gaussian filter which width is optimization parameter d. Calculating cost objective δ j as the difference between the correlation function of the input image and the correlation function of a particular slice of reconstruction.

10:
if δ j is min(δ(1), ..δ(j)) then Satisfying constraints in real space: Updating parameters a, b, c, d according to the optimization algorithm. For example, the correlation function shows good performance for the solid sphere mentioned above, even in the case of the wrong guess of SDTF. Even in case of inappropriate guess for static DTF convolution can be successfully modified for satisfying constraints for 2d-slices ( figure 6). There is no need to optimize parameters a and b due to the linear transition function for porous media. However, angular coefficient c and filter width d still need to be determined. In this case, the algorithm is the following. Replacement its magnitude with the magnitude obtained from self-convolution.

7:
Frequency filtration with Gaussian filter which width is optimization parameter d. Thresholding (m 1 s ) j to conserve porosity of input image.
where th is thresholod and (m 1 s ) 2d is input image.

10:
Calculating cost objective δ j as negative likelihood function (appendix A1) for surface correlation function of the input image and all the slices of reconstruction.

11:
Updating parameters c, d according to the optimization algorithm. While Algorithms 5 and 4 look very similar, they possess some principal differences, the initial guess being the most pronounced one. In Algorithm 5, it is always the same random noise. On the other hand, in Algorithm 4, the random noise is used as the first guess, and each next iteration uses the best previous results as input. Such an implementation allows combining optimization steps consisting of phase recovery with different parameters to equivalent phase-recovery with a more significant number of iterations. The influence of this approach is significant in cases when the number of phase recovery iterations per one optimization step is not enough to reconstruct microstructure from random noise. This seems unnecessary for general porous media (unlike the ball reconstruction problem) due to the less localized distribution of voxels resulting in less iteration number necessary for the convergence. Note that a simple increase in iterations number per optimization step would resolve this problem but is expected to increase overall iterations. Reconstruction of Turing pattern (which is a complex structure often used for algorithm attestation) by algorithm 5 leads to the following results ( figure 7). It can be seen that the slice of the reconstruction, as well as the original, contains both isolated solid and void areas of similar size, but the thickness of patterns is rarely conserved.

The classical simulated annealing reconstruction technique
In order to benchmark the new methodology developed in this work, we compare the reconstruction results against a classical reconstruction based on correlation functions. To make both techniques equal, we applied modified Yeong-Torquato reconstruction [5] and utilized only the two-point probability function S 2 (r) [3] describing the probability that two points separated by a vector displacement r(x 1 , x 2 ) between x 1 and x 2 lie in the pore phase. We calculate S 2 (r) functions in two orthogonal and two diagonal directions, thus, giving four independent CFs for each input 2D image, which are then used separately during reconstruction [29], which are averaged to estimate the other 4 CFs for the unknown third dimension. For any set of correlation functions considered in the Yeong-Torquato technique, matching correlation functions of a given realization with a target CFs set is based on pixel permutations. If a set of two-point correlation functions used in reconstruction is provided in the form of f α 2 (r), where α is a type of CF (different directions in our case) and r is a segment of varying length, the difference between two realizations of the structure can be expressed as the sum of squared differences [5,29]: where f α 2 (r) andf α 2 (r) are the values of the correlation function sets for two realizations (where the former represents a reference set while the latter represents current reconstruction state), w α are weighting parameters chosen based on the input of each CF into energy E for disordered structure [29]. The energy E in eq. 5 is minimized by the simulated annealing optimization. The Metropolis algorithm is used to determine the probability of accepting any permutation p: where T is the temperature of the system, and At initialization, the temperature is chosen so that the probability p for ∆E 0 equals 0.5 [5]. The following cooling schedule based on geometrical progression is used: where k is the time step and λ is a parameter smaller than but close to unity (λ=0.999999 for all reconstructions presented here). An optimized Capek pixel permutation approach based on interface choices [37,38] was used. Periodic boundary conditions were applied for CFs evaluation. The reconstruction procedure was terminated after 10 6 consecutive unsuccessful permutations.

E. Flow simulation
Fluid flow in pore space is described by the Navier-Stokes equation (9) [39].
where v= (v x , v y , v z ), µ, ρ and p are velocity field, viscosity, density and pressure field respectively. In the case of small Reynolds numbers Re = ρvl µ << 1, which is typical for flow in porous media, equation 9 can be transformed into the form (10) where domain Ω is pore space and ∂Ω is its boundary. Permeability K can be determined form Darcy's law (11) where µ is viscosity, L is the distance for which a pressure difference ∆p is applied, and Q is flow rate through the cross-sectional area S. The set of equations (10) has been solved FDMSS solver with 4th order spatial accuracy scheme [40]. To ensure good convergence in all modelling cases either 3.510 3 iterations were performed or we reached error criterion < 0.05 which is based on imbalance for both continuity and motion parts of equation (10).

F. Samples for reconstruction and comparison between techniques
To verify the newly proposed modified phase-retrieval algorithm and compare it against classical technique, we chose three 3D porous media images of different genesis: artificial ceramic [41], sandstone and carbonate rocks [42]. The choice was motivated by a wide range of porosities within these samples and their relative homogeneity and isotropy. Such properties are essential, as we assume the structures to statistically homogeneous so that we can consider r as a scalar distance between pixels while computing correlation functions and isotropy is needed to evaluate the unknown 3rd dimension to perform 2D into 3D reconstruction. We chose a single 2D slice with the porosity value closest to the porosity of the whole original 3D image as input data to both methods. Both simulated annealing and phase retrieval reconstructions are performed with periodic boundary conditions. To compare the reconstruction against the original 3D image, we compute error based on flow simulations using the following definition: where K reconstruction is the permeability of the stochastic reconstruction at hand, and K original is the simulated permeability based on the original 3D XCT image.

III. RESULTS AND DISCUSSIONS
Visual comparison between the original binarized XCT 3D images and their reconstructed replicas revealed no particular information regarding differences between phase-retrieval and simulated annealing methods ( figure  8). To the naked eye, it seems that SA performs somewhat better for ceramic sample, while PR resembles carbonate and sandstone originals in a superior fashion. In contrast to unmodified phase-retrieval from 2D to 3D images [28], our approach produces no apparent artifacts or distortions. Permeabilities of the stochastic replicas provide a much more robust benchmark 9. We immediately observe that both reconstruction methods underestimate permeability -this is to be expected as they are mainly based on the S 2 correlation function, which provides only a limiting information content [43,44] for complex porous media images used for the testing. For high porosity ceramic sample, PR replicas had higher permeabilities than SA, with one of the reconstructions being very similar to the original 3D image. For low porosity sandstone and carbonate samples, SA proved to provide consistently lower errors. To understand the permeability results better and to make the comparison between classical SA and developed PR method fully quantitative morphologically, further analysis in terms of L 2 , SS 2 , and C 2 correlation functions is presented in figures 10, 12. For ceramic sample, PR has higher connectivity (based on C 2 ) and larger pores (based on L 2 ). This results in high permeability values as explained by large and connected pores (visible on 8). For the carbonate sample, PR outperforms SA based on all additional CFs. Correlation functions for sandstone replicas reveal no obvious winner, as phaseretrieval shows better L 2 statistics while diverging more compared to annealing in terms of C 2 . This overall situation is rather interesting, as PR replicas being better in terms of CFs did not guarantee their permeability error K error lower that for SA reconstructions. Note that the S 2 correlation function's match was always perfect by design (both PR and SA minimized the difference).
Putting all the results from 8-12 together, it is impossible to choose either PR or SA as a more favorable approach, and this was not the target of the comparison. All in all, we conclude that both approaches are compa- rable in terms of their accuracy and perform differently for various sample structures. We iteratively improve the structure by voxel permutations with annealing, while phase-retrieval is more of a sequential algorithm. As SA swaps voxels preferentially to the interfaces between binary phases for improved convergence [37,45,46], we introduced surface-surface optimization [47] within PR to counterbalance such permutations. Without this augmentation, i.e., if based on the S 2 correlation function alone, PR performed slightly worse than SA. Thus, the current version of the comparison is deemed fair. Based on the large variety of images used in this work, we are safe to conclude that annealing with S 2 computed in orthogonal and diagonal directions produces consistently better results as compared to full-map S 2 self-convolution and, thus, annealing can be considered to be a polishing step [48]. While classic SA and modified PR method as described here do provide quite similar reconstructions in terms of balanced morphological and permeability assessment, phase-retrieval is computationally much less taxing. Based on our Matlab implementation, phase-retrieval with 40 optimization steps with 20 iterations each took approximately 40 mins for 300 3 voxels replicas. In contrast, simulated annealing took around 2 h for the same problem. While comparison of SA-based reconstruction times is always hampered by stopping criteria, these two 2 h was the average time needed to overcome PR in accuracy based on S 2 correlation statistics used as optimization targets in both algorithms. We believe that our PR implementation has much potential for speeding up and cannot be compared directly against legacy C++ implementation (the SA method). Another interesting aspect of the DDTF approach is its ability to utilize transition functions (4) to create stochastic yet highly deterministic 3D structures from 2D images. For example, this ability could be used to re-create 3D grain morphology from 2D cuts. This way, phase-retrieval can be coupled with process-based reconstructions [25,49] to create three-dimensional shapes during the packing process. Most importantly, in this work, we have explicitly shown that phase-retrieval can utilize more sophisticated structural statistics beyond the S 2 two-point correlation function -we tried SS 2 and C 2 so far, or other morphological metrics such as fractal dimensions [50,51]. This means that fast 3D reconstructions based on modified PR with the incorporation of other additional correlation functions can be used as input data for the SA algorithm [48] or building block for MPS reconstruction [6,20,52], as well as the abovementioned hybridization with process-based methods. Such potential hybridization of different methods and our modification that makes phase-retrieval similar to other main-stream reconstruction techniques paves the way to the overall unification of all approaches, as they solve the very similar optimization problem.

IV. CONCLUSIONS
In this work, we developed a dynamic phase-retrieval stochastic reconstruction algorithm for creating 3D repli-cas from 2D original mages -DDTF. The DDTF is free of artifacts characteristic of previously proposed phaseretrieval techniques. While based on two-point S 2 correlation functions, any correlation function (surface-surface correlation was utilized here) or other morphological metrics can be accounted for during the reconstruction, thus, paving the way to the hybridization of different reconstruction techniques. To test DDTF, we performed reconstructions for three binary porous media samples of different genesis. Based on computed permeability and connectivity (C 2 and L 2 correlation functions), we have shown that the proposed technique in terms of accuracy is comparable to the classic simulated annealing-based reconstruction method but is computationally very effective. This opens the possibility of utilizing DDTF to produce fast or crude replicas further improved by other reconstruction techniques such as SA or process-based methods. Improving the quality of reconstructions based on phase-retrieval by adding additional metrics into the reconstruction procedure is a research activity with a high potential to enhance the methodology significantly.

ACKNOWLEDGMENTS
This work was supported by Russian Science Foundation grant №19-72-10082 (M.V.K. and K.M.G.). The generous cooperation of the authors is within the Flow and Transport in Media with Pores research group (FaT iMP, www.porenetwork.com) and uses some of its software.

Appendix A: Maximal likelihood
The way to estimate probability that certain discredized curve y certatain (x i ) belongs to set of curves y j (x i ), ..., y N (x i ) is the maximal likelihood method. The cost objective δ needed to be minimized in this case is the following.
where mean y mean and standard deviation σ i are defined in A3 and A2 respectively.