Posterior regularization method for phase removal of shale nano-structure imaging in space domain

X-Ray computed tomography is a non-destructive method that is used, among many applications, to study the size, shape, 3D structures and interconnections of pores in shale. We use phase retrieval methods to deal with the “edge enhancement” effect caused by phase shift. The process of phase retrieval can be described by the transport-of-intensity equation (TIE). But this is an ill-posed problem. The existing methods focus on phase retrieval in the frequency domain. To tackle the ill-posedness, we propose a new method whose main idea is to solve this problem in space domain with a regularization technique. We study a synthetic shale model and simulate the projection data. Then we apply three methods to retrieve the phase: conventional method in frequency domain, direct solving method and iterative Tikhonov regularization method in space domain. Finally, we use the standard filtered back-projection (FBP) method to present the outcome. By analyzing the results, we find advantages of the new method: more stability and fewer artifacts under noise perturbations. The study shows that relative errors of the new method are nearly 1% of that of the traditional method based on frequency domain, and hence the new method is promising for the practical data processing.

Introduction

X-ray computed tomography (CT) is an important diagnostic approach in healthcare and is a non-destructive detecting technology. Due to its particular advantages, it is widely used in industry and medical imaging. CT relies on X-ray flux measurements from different angles to form an image. If the input X-ray photons are mono-chromatic, the X-ray intensities follow the Beer-Lambert law from the microscopic view (Natterer, 2001): $\Delta I/I=f\left(x\right)\Delta x.$ In this equation, $f\left(x\right)$ is the attenuation coefficient (function) of the object, $\Delta x$ is the small distance that the ray travels and $\Delta I$ is the intensity loss. Let $I_0$ be the initial intensity of the beam; $L$ be the straight line, and let $I_1$ be intensity of the beam after having passed the object. It follows from above equation that $I_1/I_0=\exp \left\{-\underset{L}{{\displaystyle \int }}f\left(x\right)dx\right\}.$ We need to find the attenuation coefficient (function) $f$ from the above integral equation, which can be used to describe the object’s structure. In practice, the integral can be regarded as a function on $R^2$ being mapped into the set of its line integrals, and it can be expressed as the Radon transform (Natterer, 2001). Therefore the reconstruction problem of CT can be performed via inversion of the Radon transform in $R^2$.

The common CT reconstruction algorithms can be divided into two categories: analytical methods and iterative methods (Natterer, 2001). In practice, we use methods like filtered back-projection algorithm (FBP) and the algebraic reconstruction technique (ART) correspondingly (Shepp and Logan, 1983; Gordon et al., 1970). The former one utilizes explicit formula which may cause ring-artifacts, whereas the latter can suppress some artifacts but requires longer computing time than the former one. Some commercially available CT scanners use fan-beam projections, but they are not applicable for the reconstruction of shale structure by using the synchrotron radiation-based X-ray CT. Synchrotron radiation allows us to select monochromatic X-ray beams from the full spectrum and conserve enough intensity for imaging. Because of the small size of the shale sample, parallel geometry can be used to conduct the experiment.

Conventional X-ray (CT) is based on the contrast of the absorption. However, a wide range of samples used in biology and material exhibit very weak absorption contrast, nevertheless can produce significant phase shifts in the X-ray beam (Bronnikov, 2006). When the distance between the sample and the CCD (CCD stands for “charged couple device”) becomes larger, the projected image becomes a mixture of the absorption-contrast image and the phase-contrast image. As for shale, it contains a lot of weak absorption organic materials which are similar with pores in the absorption-contrast image. Although the conventional X-ray CT algorithm may not distinguish the difference between the pore and weak absorption materials, the phase shift that X-ray undergoes when passing through the material may provide better image quality. Experiments have demonstrated that the intensity distribution changes a lot with different distances between the sample and the CCD (Bronnikov, 2006). This leads to a fact that the Radon transform cannot be used to simulate the practical process of the CT imaging exactly. Therefore, it is necessary to use the information about the phase. The “edge enhancement” effect caused by the phase shift will affect the final imaging (Chen et al., 2009). So we should use the phase information to reduce these artifacts.

X-ray phase contrast imaging has received much more attention in recent years, either in biomedical sciences (Bravin et al., 2013) or in material sciences (Mayo et al., 2012a; 2012b). The object can be depicted by its refractive index distribution, $n\left(r\right)=1-p\left(r\right)+i\beta \left(r\right)$, where $p$ denotes the phase shift caused by the object, the imaginary part $\beta$ describes the absorption, and the vector $r$ is the spatial coordinate. The parameter $\beta$ has a relationship with the linear attenuation coefficient $\mu$: $\mu =4\pi /\left(\lambda \cdot \beta \right)$. Based on different kinds of optical devices, there are many commonly used phase retrieval methods like: in-line phase contrast (Lee, 2015), grating interferometer (Nesterets and Wilkins, 2008), analyzer-based method (Bravin, 2003), and coded aperture (Olivo et al., 2012). The advantages of the in-line phase contrast imaging are its simple arrangement of optical components and insensitiveness to misalignment. And because there is no mask like interferometer or analyzer in this method, the flux will not be lost (Lee, 2015).

The computational procedures of the X-ray phase contrast imaging can be outlined as follows (Lee, 2015).

• Phase retrieval

- Obtain the measured intensity function

- Calculate the Fourier transform of the measurement

- Apply a user-defined filter in frequency domain

- Perform the inverse Fourier transform of the filtered result

- Obtain the phase image through FBP or inversion

• Filtered back-projection (FBP)

- Calculate the Fourier transform of the phase

- Frequency filtering

- Calculate the inverse Fourier transform of the filtered result

- Back-projection

In above steps, the phase retrieval is used for dealing with the edge effect, while the FBP is used to present the image results.

There are two basic assumptions when we use the information of phase to retrieve the structure of the object: the absorption is negligible ($\mu \approx 0$) (Bronnikov, 1999), and the phase and absorption are proportional to each other ($p\propto \beta$) (Beltran et al., 2002; Paganin et al., 2002; Wu and Liu, 2005). When the object is shale, the absorption of the light is inevitable, so we should use the latter assumption. We can use the phase information to rectify the final image through the proper choice of the ratio between parameter values $p$ and $\beta$. Even though there exit some exceptions like Shack-Hartman sensor (Lane and Tallon, 1992), most of the existing methods for the analysis of X-ray phase contrast data requires application of a filter in the Fourier domain just like the procedures listed above (Bronnikov, 1999; Beltran et al., 2002; Paganin et al., 2002; Wu and Liu, 2005). Some experiments have shown that the X-ray in-line phase contrast CT technology based on the modified Bronnikov method can be used in the experimental investigation of high polymer mixed materials and medical samples (Liu et al., 2012). Synchrotron radiation X-ray propagation-based phase-contrast computed tomography has been used to reconstruct the structure of plant organs (Ye et al., 2013). The essence of these problems is solving the transport of intensity equation (TIE) (Teague et al., 1983). In previous researches, people have studied the convergence of the regularized solutions of this question in frequency domain (Sixou, 2015). The iterative Tikhonov regularization method was developed for solving the TIE in spatial domain (Tang and Wang, 2017). In Wang and Tang (2018), method and device for nano-scale imaging with spatial domain technique was also addressed.

In this paper, we propose a new method based on space domain. In this new method, the problem can be expressed as $u_e=\mathcal{A}f+\text{error}$, where $\mathcal{A}$ is an operator with second-order differential term, $u_e$ is the observation data and $f$ is the projection data which only contains the information of the absorption. The model can be expressed as a least squares variational problem ${\min }_{f}T\left[f\right]:=\frac{1}{2}{‖u_e-\mathcal{A}f‖}_{L_2}^2,$ where“$:=$” means “defined by”. The minimal norm least squares solution can be expressed as $f={\mathcal{A}}^{+}{u}_e$, where ${\mathcal{A}}^{+}$ denotes the pseudo-inverse of $\mathcal{A}$. However, we find that with the improvement of the model’s resolution, the condition number of the operator $\mathcal{A}$ becomes huge. And the spectrum of the self-adjoint operator ${\mathcal{A}}^{*}\mathcal{A}$ is much worse than that of the operator $\mathcal{A}$; so when noise exists, it is hard to form a correct image. This paper discusses the ill-posedness of the phase retrieval problem and proposes a posterior regularization method for phase retrieval in space domain.

We will compare three methods: traditional transformed-based frequency-domain method, direct normal equation solving method (equaling to the least-squares method), and a posterior regularization method. To be simplicity in context, the above three methods are named as TFDM, LSM and PRM, respectively. These methods will be addressed in details in the following sections. Finally, some conclusions are given.

Theory

Transport-of-intensity equation

Let $\left(x_1,x_2,x_3\right)$ be a Cartesian coordinate system of the object, the detector is located with the $\left(x,y\right)$ coordinate system, and the direction of the $y$ axis coincides with the direction of the $x_3$ axis (see Figure 1).

FIGURE 1

FIGURE 1. Coordinate system for the object and the plane of observation.

We assume that the source $U_i$ is the incident monochromatic plane wave (considering the synchrotron radiation generates a highly collimated beam, this assumption is easy to achieve). The wave field downstream of the object can be described as

$U_{\theta }\left(x,y\right)=T_{\theta }\left(x,y\right)U_i,(1)$

where $T_{\theta }\left(x,y\right)$ is the complex transmission function of the object for the scanning angle $\theta$, which contains the absorption of the light and the phase shift and is of the form

$T_{\theta }\left(x,y\right)=\exp \left[-\frac{1}{2}{\mu }_{\theta }\left(x,y\right)+i{\varphi }_{\theta }\left(x,y\right)\right].(2)$

The absorption function ${\mu }_{\theta }$ and the phase shift ${\varphi }_{\theta }$ follow the Beer-Lambert’s law and have the following integral formulas:

${\mu }_{\theta }\left(x,y\right)=\int \mu \left(x_1,x_2,y\right)\delta \left(x-x_1\cos \theta -x_2\sin \theta \right)d{x}_{1}d{x}_2,(3)$

${\varphi }_{\theta }\left(x,y\right)=\frac{2\pi }{\lambda }\int g\left(x_1,x_2,y\right)\delta \left(x-x_1\cos \theta -x_2\sin \theta \right)d{x}_{1}d{x}_2,(4)$

where $\mu$ is the linear attenuation coefficient of the object, $g\left(x_1,x_2,x_3\right)=\text{real}\left[n\left(x_1,x_2,x_3\right)\right]-1$, where $\text{real}\left(.\right)$ means the real part and $\delta \left(.\right)$ is the Dirac delta function. The two equations describe the 2D Radon transforms, and the intensity of the plane which has a distance $d$ from the object can be expressed as (Bronnikov, 2015)

$I_{\theta }^d\left(x,y\right)={\left|h_d**U_{\theta }\right|}^2,(5)$

where the double asterisk “**” denotes two-dimensional convolution and $h_d$ is the Fresnel propagator which can be defined as:

$h_d\left(x,y\right)=\frac{\exp \left(ikd\right)}{i\lambda d}\exp \left[i\frac{\pi }{\lambda d}\left(x^2+y^2\right)\right],(6)$

where $\lambda$ is the wavelength. If we consider the near field condition $\lambda d\ll D^2$ ($D$ is the size of the object) and assume that ${\mu }_{\theta }$ varies insignificantly, we obtain the intensity distribution at a small distance as (Brenner, 1952; Bronnikov, 1999):

$I_{\theta }^d\left(x,y\right)=I_{\theta }^0\left[1-\frac{\lambda d}{2\pi }{\nabla }^2{\varphi }_{\theta }\right],(7)$

where $I_{\theta }^0$ denotes the intensity at the object plane. From Eq. 7, we find that the projected image consists of the absorption information $I_{\theta }^0$ and the phase information ${\varphi }_{\theta }$, which may generate the “edge enhancement” effect. With the assumption that ${\mu }_{\theta }$ varies slightly, and let $T\left(r\right)$ be the projected thickness of the object, we obtain the TIE equation (Paganin et al., 2002)

$\left(-\frac{dp}{\mu }{\nabla }^2+1\right)e^{-\mu T\left(r\right)}=\frac{I_{\theta }^d}{I^{\text{in}}}.(8)$

The right-hand side of Eq. 8 is the observation, the left represents the effect caused by the phase. The conventional method to deal with this problem is to take the Fourier transform to both sides of the above equation to get (Paganin et al., 2002)

$\mathcal{F}\left\{e^{-\mu T\left(r\right)}\right\}=\mu \frac{\mathcal{F}\left\{I_{\theta }^d\right\}/I^{\text{in}}}{dp{\left|k_{\perp }\right|}^2+\mu },(9)$

where $k_{\perp }$ is the coordinate in space domain, the symbol “$\perp$” stands for the plane perpendicular to the direction of the X-ray travels. This method is equivalent to using a frequency filtering to get the projected thickness $T\left(r\right)$. Finally, we can use the filtered back projection to get the parameter of the model. We call the transformed-based frequency-domain method as “TFDM”.

Instead of this final step, here we solve Eq. 8 directly in the space domain. To do so, we first build an operator equation to reformulate the Eq. 8 to get $e^{-\mu T\left(r\right)}$, then use the filtered back projection to get the projected thickness of the object. Equation 8 can be expressed in the form of the first kind operator equation as:

$u=\mathcal{A}f,(10)$

where $\mathcal{A}=\left(1-\frac{dp}{\mu }{\nabla }^2\right)$, $f=e^{-\mu T\left(r\right)}$ and $u=I_{\theta }^d/I^{\text{in}}$.

General solutions of the operator equation of the first kind

All data are noisy, so Eq. 10 should be rewritten as:

$u_e=\mathcal{A}f+\text{error}.(11)$

Therefore we cannot guarantee that $u_e$ belongs to the range of the operator. So we need to find the approximate solution that minimizes the residual norm of the simulated data to the observed data on some metric space:

$f:=\underset{f}{\mathrm{argmin}}{‖\mathcal{A}f-u_e‖}^2,(12)$

where “$\mathrm{argmin}$” stands for argument of the minimum. Because $\mathcal{A}$ includes the differential operator (see Eq. 10) and its spectrum is unbounded, the discrete equations derived by the difference scheme or the finite element method are always ill-conditioned. This ill-conditioning will be increasingly severe with the shrinking of the discrete scaling, therefore the phase retrieval problem is an ill-posed problem (Xiao et al., 2003). The strategy to solve the ill-posed problem is to find the solution of a series of well-posed problems.

This method for solving ill-posed problems is called a Tikhonov regularizing method (Tikhonov and Arsenin, 1977; Wang et al., 2010). Tikhonov regularization means a variational approach based on minimization of a smoothing functional as follows (Xiao et al., 2003)

$M^{\alpha }\left[f,u_e\right]={\rho }_U^2\left(\mathcal{A}f-u_e\right)+\alpha \mathrm{\Omega }\left[f\right],(13)$

where ${\rho }_U$ is a measure in the observation space $U$ and $\mathrm{\Omega }\left[f\right]$ is a non-negative and continuous functional defined in the solution space $F$.

Numerical methods for minimizing $M^{\alpha }\left[f,u\right]$ are based in general on three approaches (Wang, 2007; Freeden et al., 2010; Tang and Wang, 2017; Wang et al., 2020):

(1) Direct methods, i.e., solving an Euler equation: ${\left(M^{\alpha }\left[f,u_e\right]\right)}^{\prime }=0$.

(2) Iterative Tikhonov regularization methods: i.e., solving the above Euler equation iteratively with higher-order convergence rate.

(3) Optimization methods, i.e., solving an optimization problem $M^{\alpha }\left[f,u_e\right]\to \min$ via gradient or Newton’s method.

We will consider a posterior regularization method with optimal regularization parameter selection. Details will be given in the following section.

Solving methodology

Discrete Tikhonov regularization

To simplify the notation, we denote the discrete form of the operator equation as

$Af=u_e,(14)$

where $A$ is the discretization of the operator $\mathcal{A}=\left(1-\frac{dp}{\mu }{\nabla }^2\right)$. For the second-order differential operator in $\mathcal{A}$, we apply the five-point finite difference scheme

$\left\{\begin{array}{l}\frac{{\partial }^2{f}_{i,j}}{\partial x^2}\approx \frac{{\delta }^{2}f}{\delta x^2}=\frac{1}{{\left(\delta x\right)}^2}\left(a_1{f}_{i,j}+a_2{f}_{i+1,j}+a_3{f}_{i+2,j}+a_4{f}_{i-1,j}+a_5{f}_{i-2,j}\right),\\ \frac{{\partial }^2{f}_{i,j}}{\partial y^2}\approx \frac{{\delta }^{2}f}{\delta y^2}=\frac{1}{{\left(\delta y\right)}^2}\left(a_1{f}_{i,j}+a_2{f}_{i,j+1}+a_3{f}_{i,j+2}+a_4{f}_{i,j-1}+a_5{f}_{i,j-2}\right).\end{array}\right.(15)$

The coefficients $a_i$ ($i=1,\cdots ,5$) can be obtained through the Taylor expansion and solving the corresponding equations as: $a_1=-1/12$; $a_2=4/3$; $a_3=-5/2$; $a_4=4/3$ and $a_5=-1/12$.

When the light travels through the air, we assume that the intensity at the detector equals the intensity at the object plane. This is the situation that happens at the boundary of the projection data. This condition means, the operator $A$ that we want, does not change the input data at the boundary, hence it satisfies the Dirichlet boundary condition. When we obtain the coefficients $a_i$ ($i=1,\cdots ,5$) derived from Taylor expansion, the finite difference matrix is asymmetric. The condition number of the matrix $A$ will be very large, which leads to unstable solutions of the corresponding linear systems. To tackle this problem, we need to find asymmetric matrix which can express the second-order differential. Considering the boundary condition we have mentioned above, we can simply adjust the first two rows and the last two rows of the matrix to achieve our purpose. This new matrix $A$ does not change the boundary value, even though it becomes a symmetric matrix. And the linear system with the symmetric matrix reads as

$\frac{dp}{\mu }\frac{1}{{\left(\mathrm{d}x\right)}^2}\left[\begin{array}{l}1/120-1/1200\cdots 0\\ 0-5/44/3-1/120\cdots 0\\ -1/124/3-5/24/3-1/12\cdots 0\\ \ddots \ddots \ddots \ddots \ddots ⋮\\ -1/124/3-5/24/3-1/12\\ 0\cdots 0-1/124/3-5/40\\ 0\cdots 00-1/1201/12\end{array}\right]\left[\begin{array}{c}f\left(r1\right)\\ f\left(r2\right)\\ f\left(⋮\right)\\ f\left(⋮\right)\\ f\left(\text{end}-1\right)\\ f\left(\text{end}\right)\end{array}\right]=\left[\begin{array}{c}b\left(r1\right)\\ b\left(r2\right)\\ b\left(⋮\right)\\ b\left(⋮\right)\\ b\left(\text{end}-1\right)\\ b\left(\text{end}\right)\end{array}\right],(16)$

where $\mathrm{d}x=\delta x=\delta y$, ri refers to the $i$-th position and the dimension of the matrix $A$ depends on the resolution of the CCD, e.g., 512-by-512, or 2048-by-2048. In our model, the size of the matrix $A$ is 512-by-512.

Using this finite difference scheme, we can construct discrete operator $A$, and the problem is transformed to the solution of a system of linear algebraic Eq. 14.

With the improvement of the spatial resolution, the condition number of the operator $A$ becomes large caused by the finite difference scheme. If we want to find the minimum norm of the residual in $l_2$ space (i.e., least-square method), we need to solve a normal equation $A^{T}Af=A^T{u}_e$. But the normal equation is still ill-conditioned, as the distribution of its spectra is heavily decentralized. We call the method of solving the above normal equation as “LSM”.

The regularization method of solving Eq. 14 refers to minimize the following objective functional,

$\underset{f}{\min }{J}^{\alpha }\left[f,u\right]:={‖Af-u_e‖}_{l_2}^2+\alpha {‖f‖}_{l_2}^2,(17)$

where $\alpha$ is a regularization parameter which is larger than 0, then a key issue is how to select a optimal value of the regularization parameter $\alpha$ in a posteriori way. We call the method of solving the above regularized problem as “PRM”. Details about a posterior selection of the regularization parameter $\alpha$ are outlined in Appendix A.

Numerical experiments

Modeling parameters

In order to verify the validity and advantages of the new algorithm, we carried out numerical simulations in the case of parallel-beam CT. The model we used is shown in Figure 2. In which, the grey values of the phantom represents the $\beta$ values of the material, which range from 0 to $20.5\times {10}^{-10}$. The parameters we used refer to the shale’s basal components (Wang et al., 2015). The size of the model is set to be $512\times 512$ pixels. The scanning angular range is $\left[0°,180°\right]$ and the step is $1°$. The energy of the X-ray is 8Kev and the wave length is 154.06pm (Zschornack, 2007). The distance between the detector and the sample is set to be 30 cm. The ratio between the phase shift and the absorption is set as1200. Let ${\text{data}}_{\text{true}}$ denote the model parameter, ${\text{data}}_{\text{res}}$ be its reconstruction. We define the relative error as follows:

$\text{Err}=\frac{{‖{\text{data}}_{\text{true}}-{\text{data}}_{\text{res}}‖}_2}{{‖{\text{data}}_{\text{true}}‖}_2}.(18)$

FIGURE 2

FIGURE 2. Linear absorption model.

The model’s main component is quartz, and the component percentages are listed as in Table 1.

TABLE 1

TABLE 1. Theoretical $\beta$ values for materials in the model at E=30 Kev.

Simulation results

The line integral projection is displayed as in Figure 3, where we discretize the model by decomposing it into $512\times 512$ pixels and assume the model parameter to be constant in each small square. We use point-by-point scanning instead of the exact mathematical description of the model because of the fringe effect. For each straight line $L\left(r_i,\theta \right)$, we add up the length of $\left(L\left(r_i,\theta \right)\cap d_j\right)$ times the absorption value in $d_j$ to get the $u\left(r_i\right)$ (see Figure 3). After this step, we apply Eq. 8 to obtain the simulation data.

FIGURE 3

FIGURE 3. Acquisition geometry.

In our experiments, as mentioned before, we compare the three methods respectively: TFDM—the frequency domain method with phase retrieval, LSM—the direct method in space domain (without regularization), and PRM-the proposed regularization method in space domain.

Noiseless data

First we consider the ideal case, i.e., the data is noiseless. As a comparison, we first use the FBP method to deal with the projection without considering phase property. The reconstruction result is illustrated in Figure 4. The Err will be 227.2 and we will only get some blurred information at the edge of the structure. Therefore, the reconstruction result is unsatisfactory. This indicates that the effect caused by the phase information is inevitable.

FIGURE 4

FIGURE 4. FBP result without phase retrieval.

We apply the aforementioned three methods to the phase retrieval problem. The reconstruction results are shown in Figures 5–7, respectively. The Errs for the three methods are 7.99, 0.018 and 0.018 correspondingly. The inversion results indicate that for noiseless data, the reconstructions from the direct method and the regularization method are similar, and are better than that of the frequency domain method.

FIGURE 5

FIGURE 5. The inversion result using TFDM for noiseless data.

FIGURE 6

FIGURE 6. The inversion result using LSM for noiseless data.

FIGURE 7

FIGURE 7. The inversion result using PRM for noiseless data.

Noisy cases

A random noise with Gaussian distribution was added to the phantom projection. Two levels of noise were used, so that two new data sets were generated. First, we try to simulate the effect caused by small noise, where the noise-level equaling 0.1% is added to the true data. Figures 8–10 display the simulation results using the TFDM, LSM and PRM with 0.1% Gaussian noise in each case. The Errs for the three methods are 4.38, 1.51 and 0.0482, correspondingly.

FIGURE 8

FIGURE 8. The inversion result using TFDM for noisy data with noise-level 0.1%.

FIGURE 9

FIGURE 9. The inversion result using LSM for noisy data with noise-level 0.1%.

FIGURE 10

FIGURE 10. The inversion result using PRM for noisy data with noise-level 0.1%.

If we add large noise (noise-level = 1%), the results with the TFDM, LSM and PRM, are shown in Figures 11–13, respectively. The Errs for the three methods are 10.85, 14.02 and 0.0627 correspondingly. To quantitatively illustrate differences of the three methods, we draw a profile in Figure 14, which shows the results calculated by above mentioned three methods, in which the black line is the model parameter (true value), the red line is the result calculated by the regularization method (PRM) proposed in this paper, the green line is the result calculated by the LSM and the blue line is the result calculated by the TFDM in the frequency domain.

FIGURE 11

FIGURE 11. The inversion result using TFDM for noisy data with noise-level 1%.

FIGURE 12

FIGURE 12. The inversion result using LSM for noisy data with noise-level 1%.

FIGURE 13

FIGURE 13. The inversion result using PRM for noisy data with noise-level 1%.

FIGURE 14

FIGURE 14. Comparison of three methods on their profiles for noisy data with noise-level 1%.

Due to the ill-posedness of the nano-scale CT imaging problem, the direct solution is not stable (LSM) when noise exists. When the noise level is 1%, it is even hard to get anything from the result as is shown in Figure 12. Under the same experiment condition, the signal-noise ratio will be low especially when the sample is in micro size. Therefore, even if the direct method has a fast computational speed, it is still hard to image well for the practical data, as noise always exists for practical measurements.

Even though the results given by the frequency domain method (TFDM) can show the structure information under some noise levels, if we observe a particular position where the pore located, the results, calculated by the frequency domain method, show significant difference. When the value of the model parameter is high, there exists serious artifact around them. So the results from the frequency method may be affected by the values of pixels around them. The worst thing for the frequency domain method is that its relative error is much higher than that of the proposed method in space domain. We also performed some numerical experiments under different situations (e.g., change the resolution, change the distance between sample and detector), the method we proposed behaved well.

Practical shale sample results

An organic-rich shale sample is selected for this study. This sample is a Qiliao shale from Chongqing Shizhu-Qiliao profile belonging to the five peak group-Longmaxi (LMX) group, China. It is composed of quartz (71%), feldspar (8%), pyrite (3%), dolomite (4%), and illite-smectite-chlorite (14%). LMX formation shale at Shizhu-Qiliao outcrop mainly consisted of dark carbonaceous shale with pyrite and quartz veins. To carry out the experiment, the Nano-TXM samples were prepared by an FEI Helios Nano-Lab 600 DualBeam FIB-SEM. Nano-TXM experiments were carried out at beamline BL01B1, the National Synchrotron Radiation Research Center (NSRRC) in Hsinchu, Taiwan (Wang et al., 2016). Equipped with Zernike-phase-contrast capability, the instrument can take images of light materials such as biological specimens. The spatial resolution of the microscope is $0.9\Delta r/m$, where $\Delta r$ is the outermost width of the zone-plate, with the diffraction order being an odd number (Song et al., 2007). They provided 2D micrograph and 3D tomography at spatial resolutions of 50 nm, with a first-order diffraction of a Fresnel zone plate at an X-ray energy of 8keV; the image field of view was 15 × 15 μm² for the first-order diffraction of the zone plate. The exposure time is 60s for a 2D image. Meanwhile, the phase term is retrieved by the Zernike’s phase contrast method. The gold phase ring is positioned at the back focal plane of the objective zone-plate retards or advances the phase of the zero-th order diffraction by π/2 resulting a recording of the phase contrast images at the detector. The angle scanning is 180°, and there are 181 snapshots that are recorded in our system. These data will be used to perform an inversion.

As an illustration, the projection data of the first slice is shown in Figure 15. If we do not consider removal of phase effect, the inversion result is shown in Figure 16. With removal phase effect, we make a comparison of three methods: TFDM, LSM and PRM. As mentioned in theoretical simulations, the proposed regularization method (PRM) performs best. The inversion results using TFDM, LSM and PRM are shown in Figures 17–19, respectively. It is clear from Figure 19 that the high resolution result is obtained using PRM. To be focus, we zoomed the same part of each figure of the reconstruction in Figure 20, it clearly shows that the proposed regularization method (PRM) yields the better reconstruction result with high resolution. We also plot the histograms of the reconstruction results in Figure 21, where each output grayscale image has 64 bins. Finally, we patch all slices using PRM together to generate a three-dimensional in Figure 22. It illustrates that honeycomb coal-like pyrite and pores are observed clearly.

FIGURE 15

FIGURE 15. Projection data of the first slice.

FIGURE 16

FIGURE 16. One slice of the reconstruction results: without phase removal.

FIGURE 17

FIGURE 17. One slice of the reconstruction results: TFDM.

FIGURE 18

FIGURE 18. One slice of the reconstruction results: LSM.

FIGURE 19

FIGURE 19. One slice of the reconstruction results: PRM.

FIGURE 20

FIGURE 20. Comparison of reconstruction results with zoomed parts: (A) without phase removal, (B) TFDM, (C) LSM and (D) PRM.

FIGURE 21

FIGURE 21. Histograms of reconstruction results: (A) without phase removal, (B) TFDM, (C) LSM and (D) PRM.

FIGURE 22

FIGURE 22. Three-dimensional display of the reconstruction results of PRM.

Discussion and conclusion

From methodological view of point, our method works for complex, finely layers systems prevalent in most shales, this is because, for shale samples, although the absorption of light cannot be ignored like biological samples, when the object leaves the receiving screen for a certain distance, the phase factor will cause certain interference to the projection data, forming an “edge enhancement” effect, which is necessary to combine the basic theory of phase contrast imaging to remove the interference of phase factors on the projection. Therefore, a phase retrieval model based on TIE equation in space domain is established. During solution of the TIE equation, since the discretization of differential operators leads to the instability of the problem, the regularization technique is employed, and a posterior regularization method in space domain (PRM) is proposed for solving the phase retrieval problem, meanwhile selection of an optimal posterior regularization parameter is adopted, which accelerates the convergence speed of the iterations. Due to an optimal regularization parameter is achieved, the solution reaches optimal. Comparison with the standard method in frequency domain (TFDM) and the direct method in space domain (LSM) shows that the regularization method yields an accurate reconstruction and is robust to noise. Though this method is theoretically sound, there is still a restriction: we assume that the ratio between the real and imaginary parts of the refraction coefficient is homogeneous which may not be practical. So when we want to deal with the real data, we should choose the suitable energy and ratio (between values $p$ and $\beta$) based on some priori information like the absorption edge of the target compositions. In addition, the above theoretical model achieves the phase retrieval in 2D, which means that we assume the model does not change at the $y$ direction in Figure 1. So when we deal with the actual data, we may apply finite difference equations in two dimensions, which leaves for future work.

We remark that other imaging technique such as scanning electron microscopy equipped with a focused ion beam (FIB-SEM) can be used for nano-scale imaging of pore structures of shale. Though the resolution is high, it belongs to the destructive method. We believe that a combination of FIB-SEM (supplying high quality actual pore volume fraction) with nano-scale CT imaging will be a good choice for large amount of data analysis, e.g., employing the artificial intelligence technique, which will be also an interesting topic.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

The research is supported by Doctoral Scientific Research Foundation of Liaocheng University (Grant No. 318051911), the Innovation Research Program of the Chinese Academy of Sciences (Grant No. ZDBS-LY-DQC003) and the key program IGGCAS-2019031.

Acknowledgments

We are grateful for constructive comments from four reviewers and the editor, which improve the quality of the paper. The Center of Big Data and AI in Earth Sciences of IGGCAS in supplying facilities for simulation and practical data processing is also acknowledged.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Beltran, M. A., Paganin, D. M., Uesugi, K., and Kitchen, M. J. (2010). 2D and 3D X-ray phase retrieval of multi-material objects using a single defocus distance. Opt. Express. 18, 6423–6436. doi:10.1364/OE.18.006423

PubMed Abstract | CrossRef Full Text | Google Scholar

Bravin, A., Coan, P., and Suortti, P. (2013). X-Ray phase-contrast imaging: From pre-clinical applications towards clinics. Phys. Med. Biol. 58, 1–35. doi:10.1088/0031-9155/58/1/R1

PubMed Abstract | CrossRef Full Text | Google Scholar

Bravin, A. (2003). Exploiting the X-ray refraction contrast with an analyzer: The state of the art. J. Phys. 36, A24–A29. doi:10.1088/0022-3727/36/10A/306

CrossRef Full Text | Google Scholar

Brenner, H. (1952). On the asymptotic evaluation of diffraction integrals with a special view to the theory of defocusing and optical contrast. Physica 18, 469–485. doi:10.1016/S0031-8914(52)80079-5

CrossRef Full Text | Google Scholar

Bronnikov, A. V. (2006). Phase-contrast CT: Fundamental theorem and fast image reconstruction algorithms. Proc. Spie. 6318, 63180Q. doi:10.1117/12.679389

CrossRef Full Text | Google Scholar

Bronnikov, A. V. (1999). Reconstruction formulas in phase-contrast tomography. Opt. Commun. 171, 239–244. doi:10.1016/S0030-4018(99)00575-1

CrossRef Full Text | Google Scholar

Bronnikov, A. V. (2002). Theory of quantitative phase-contrast computed tomography. J. Opt. Soc. Am. A 19, 472–480. doi:10.1364/JOSAA.19.000472

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, R. C., Yangle, W., Tiqiao, X., Huitong, X., Zhou, G., and Deng, B. (2009). Preliminary results for X-ray phase contrast micro-tomography on the biomedical imaging beamline at SSRF. Nucl. Phys. 32, 241–245. doi:10.1088/1674-1137/33/8/010

CrossRef Full Text | Google Scholar

Fan, S. F., Ma, Q. H., and Wang, Y. F. (2006). A new method for choosing regularization parameter with perturbed operator and noisy data. J. Beijing Normal Univ. Nat. Sci. 42, 25–31. doi:10.1007/s10444-011-9203-6

CrossRef Full Text | Google Scholar

Freeden, W., Nashed, M. Z., Sonar, T., and Heine, C. (2010). Handbook of geomathematics. Heidelberg, Germany: Springer.

Google Scholar

Gordon, R., Bender, R., and Herman, G. T. (1970). Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography. J. Theor. Biol. 29, 471–481. doi:10.1016/0022-5193(70)90109-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Lane, R. G., and Tallon, M. (1992). Wave-front reconstruction using a Shack-Hartmann sensor. Appl. Opt. 31, 6902–6908. doi:10.1364/AO.31.006902

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, P. C. (2015). Phase retrieval method for in-line phase contrast X-ray imaging and denoising by regularization. Opt. Express. 23, 100668–110679. doi:10.1364/OE.23.010668

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H. Q., Wang, Y., Ren, Y., Xue, Y., He, Y., and Guo, H. (2012). Investigation on X-ray micro-computed tomography suitable for organic compound samples based on modified Bronnikov algorithm. Acta. Opt. Sin. 32, 320–326. doi:10.3788/aos201232.0434001

CrossRef Full Text | Google Scholar

Mayo, S. C., Stevenson, A. W., and Wilkings, S. W. (2012a). In-line phase-contrast x-ray imaging and tomography for materials science. Mater 5, 937–965. doi:10.3390/ma5050937

CrossRef Full Text | Google Scholar

Mayo, S. C., Tulloh, A. M., Trinchi, A., and Yang, Y. S. (2012b). Data-constrained microstructure characterisation with multi-spectrum X-ray micro-CT. Microsc. Microanal. 18, 524–530. doi:10.1017/S1431927612000323

PubMed Abstract | CrossRef Full Text | Google Scholar

Natterer, F. (2001). The mathematics of computerized tomography. New York, USA: Society for Industrial and Applied Mathematics.

Google Scholar

Nesterets, Y. I., and Wilkins, S. W. (2008). Phase-contrast imaging using a scanning-double-grating configuration. Opt. Express. 16, 5849–5867. doi:10.1364/OE.16.005849

PubMed Abstract | CrossRef Full Text | Google Scholar

Olivo, A., Diemoz, P. C., and Bravin, A. (2012). Amplification of the phase contrast signal at very high x-ray energies. Opt. Lett. 37, 915–917. doi:10.1364/OL.37.000915

PubMed Abstract | CrossRef Full Text | Google Scholar

Paganin, D., Mayo, S. C., Gureyev, T. E., Miller, P. R., and Wilkins, S. W. (2002). Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object. J. Microsc. 206, 33–40. doi:10.1046/j.1365-2818.2002.01010.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Shepp, L. A., and Logan, B. F. (1983). The Fourier reconstruction of a head section. J. Opt. Soc. Am. 73, 1434–1441. doi:10.1109/TNS.1974.6499235

CrossRef Full Text | Google Scholar

Sixou, B. (2015). Regularization methods for phase retrieval and phase contrast tomography. Abstr. Appl. 1–11. doi:10.1155/2015/943501

CrossRef Full Text | Google Scholar

Song, Y. F., Chang, C. H., Liu, C. Y., Chang, S. H., Jeng, U. S., Lai, Y. H., et al. (2007). X-ray beamlines for structural studies at the NSRRC superconducting wavelength shifter. J. Synchrotron Rad. 14, 320–325. doi:10.1107/S0909049507021516

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, W., and Wang, Y. F. (2017). Iterative regularization methods for phase retrieval TIE equation in space domain. Chin. J. Geophys. 60 (5), 1851–1860. doi:10.6038/cjg20170520

CrossRef Full Text | Google Scholar

Teague, M. R. (1983). Deterministic phase retrieval: A green’s function solution. J. Opt. Soc. Am. 73, 1434–1441. doi:10.1364/JOSA.73.001434

CrossRef Full Text | Google Scholar

Tikhonov, A. N., and Arsenin, V. Y. (1977). Solutions of ill-posed problems. New York, USA: John Wiley & Sons.

Google Scholar

Wang, Y. D., Yang, Y., Liu, K., Ren, Y., Tan, H., and Deng, B. (2015). Quantitative and multi-scale characterization of pore connections in tight reservoirs with micro-CT and DCM. Bull. Mineral. Pet. Geochem. 34, 86–92. doi:10.3969/j.issn.1007-2802.2015.01.010

CrossRef Full Text | Google Scholar

Wang, Y. F. (2007). Computational methods for inverse problems and their applications. Beijing, China: Higher Education Press.

Google Scholar

Wang, Y. F., Fan, S. F., Leonov, A. S., Lukyanenko, D. V., Yagola, A. G., Wang, L. H., et al. (2020). Non-smooth regularization and fast optimization algorithm for micropore reconstruction of shale. Chin. J. Geophys. 63 (5), 2036–2042. doi:10.6038/cjg2020M0684

CrossRef Full Text | Google Scholar

Wang, Y. F., and Tang, W. (2018). Method and device for nano-scale imaging based on spatial phase retrieval technique. Patent 2, ZL201610566803. doi:10.1016/j.ultramic.2011.10.012

CrossRef Full Text | Google Scholar

Wang, Y. F., and Xiao, T. Y. (2001). Fast realization algorithms for determining regularization parameters in linear inverse problems. Inverse. Probl. 17, 281–291. doi:10.1088/0266-5611/17/2/308

CrossRef Full Text | Google Scholar

Y. F. Wang, A. G. Yagola, and C. C. Yang (Editors) (2010). (Heidelberg, Germany: Springer).Optimization and regularization for computational inverse problems and applications.

Google Scholar

Wang, Y., Pu, J., Wang, L. H., Wang, J. Q., Jiang, Z., Song, Y. F., et al. (2016). Characterization of typical 3D pore networks of Jiulaodong formation shale using nano-transmission X-ray microscopy. Fuel 170, 84–91. doi:10.1016/j.fuel.2015.11.086

CrossRef Full Text | Google Scholar

Wu, X. Z., and Liu, H. (2005). X-ray cone-beam phase tomography formulas based on phase-attenuation duality. Opt. Express. 13, 6000–6014. doi:10.1364/OPEX.13.006000

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, T. Y., Yu, S. G., and Wang, Y. F. (2003). Numerical methods for the solution of inverse problems. Beijing, China: Science Press.

Google Scholar

Ye, L. L. (2013). X-ray phase contrast micro-tomography and its application in quantitative 3D imaging study of wild ginseng characteristic microstructures. Acta. Opt. Sin. 33, 365–370. doi:10.3788/AOS201333.1234002

CrossRef Full Text | Google Scholar

Zschornack, G. (2007). Handbook of X-ray data. New York, USA: Springer.

Google Scholar

Appendix A

A posterior selection of the regularization parameter.

The minimization problem (17) is equivalent to solving the linear equations:

$\left(A^{T}A+\alpha I\right)f=A^T{u}_e,(\mathrm{A1})$

where $\alpha >0$. Clearly $f_e^{\alpha }={\left(\alpha I+A^{T}A\right)}^{-1}{A}^T{u}_e$ can be an approximation of the exact solution $f_T=A^{+}u$ for a particular choice of the regularization parameter $\alpha$, where $A^{+}$ denotes the Moore-Penrose generalized inverse. Usually there are two ways to choose the regularization parameter (Tikhonov and Arsenin, 1977; Wang, 2007):

(1) A prior way: it requires knowledge of some source condition about the solution, e.g., $f_T=A^{+}u\in \text{Range}{\left(A^{T}A\right)}^{\nu }$ for some $\nu >0$. Clearly this is hard to satisfy.

(2) A posterior way: we can choose an optimal parameter ${\alpha }^{\text{opt}}$ to balance the trade-off between the least squares and the regularized term. The main process will be described below.

We suppose the noise level $e$ is bounded and satisfies $‖u-u_e‖\le \mathrm{e}\le ‖u_e‖$. Let $f\left(\alpha \right)=f_e^{\alpha }$, then the regularization parameter $\alpha$ satisfies the non-linear equation (Wang and Xiao, 2001; Xiao et al., 2003; Fan et al., 2006):

$\phi \left(\alpha \right)={‖Af\left(\alpha \right)-u_e‖}^2-e^2=0.(\mathrm{A2})$

Since $\phi \left(\alpha \right)$ is an infinitely differentiable function and is monotone increasing, there must be a root such that the above non-linear equation holds. To accelerate convergence of finding a root, we consider the Newton’s root-finding method, which can be written as

${\alpha }_{new}=\alpha -\frac{\phi \left(\alpha \right)}{{\phi }^{\prime }\left(\alpha \right)}.(\mathrm{A3})$

The derivative of the non-linear function $\phi \left(\alpha \right)$ can be calculated through

${\phi }^{\prime }\left(\alpha \right)=-2\alpha \left(\frac{\mathrm{d}}{\mathrm{d}\alpha }\left[f\left(\alpha \right)\right],f\left(\alpha \right)\right),(\mathrm{A4})$

where $\left(,.,\right)$ stands for the inner product.

To solve for the parameter $\alpha$, we need to solve the following linear equations:

$\left\{\begin{array}{l}\left(A^{*}A+\alpha I\right)f\left(\alpha \right)=A^{*}{u}_e,\\ \left(A^{*}A+\alpha I\right)f^{\prime }\left(\alpha \right)=-f\left(\alpha \right).\end{array}\right.(\mathrm{A5})$

Thus we use the Newton iterative formula for the parameter $\alpha$:

${\alpha }_{k+1}={\alpha }_k-\frac{{‖Af\left({\alpha }_k\right)-u_e‖}^2-e^2}{-2{\alpha }_k\left(f^{\prime }\left({\alpha }_k\right),f\left({\alpha }_k\right)\right)}.(\mathrm{A6})$

Having specified the system matrix $A$ and the data $u_e$, we can outline the implementation of the algorithm as follows:

(1) Input initial values of the regularization parameter ${\alpha }_0$, error-level $e$, tolerance $\epsilon >0$, maximum iteration number $k_{\max }$ and the maximum scanning angle ${\theta }_{\max }$; Input the matrix $A$ and the data $u_{e,\theta }$; Set $\theta :=0$ and $k:=0$;

(2) Solve the Eq. (A5) using the Gaussian eliminations;

(3) Compute the function values $\phi \left({\alpha }_k\right)$ and ${\phi }^{\prime }\left({\alpha }_k\right)$:

$\left\{\begin{array}{l}\phi \left({\alpha }_k\right)=\left(Af\left({\alpha }_k\right)-u_{e,\theta },Af\left({\alpha }_k\right)-u_{e,\theta }\right)-e^2,\\ {\phi }^{\prime }\left({\alpha }_k\right)=-2{\alpha }_k\left(f^{\prime }\left({\alpha }_k\right),f\left({\alpha }_k\right)\right);\end{array}\right.(\mathrm{A7})$

(4) Compute the parameter ${\alpha }_{k+1}$:

${\alpha }_{k+1}={\alpha }_k-\frac{\phi \left({\alpha }_k\right)}{{\phi }^{\prime }\left({\alpha }_k\right)};(\mathrm{A8})$

(5) If $‖{\alpha }_{k+1}-{\alpha }_k‖<\epsilon$ or $k=k_{\max }$, GOTO Step 6; Otherwise, let $k:=k+1$, and GOTO Step 2;

(6) Output $f_{\theta }^{*}=f_{{\alpha }_k}$. If $\theta ={\theta }_{\max }$, stop; Otherwise, $\theta :=\theta +1$ and GOTO Step 2.