Maximum Likelihood based Phase-Retrieval using Fresnel Propagation Forward Models with Optional Constraints

In XPCT, the measured data is sensitive to the 3D variations in the absorption index, $\beta(u,v,w)$, and refractive index decrement, $\delta(u,v,w)$, of the imaged object. Here, $(u,v,w)$ represents the 3D Cartesian coordinate system such that the X-ray propagation is along the $w-$axis and perpendicular to the $u-v$ plane. The absorption index and refractive index decrement are expressed compactly as the complex refractive index $n(u,v,w)=1-\delta(u,v,w)+i\beta(u,v,w)$, where $i$ is the imaginary unit.

As X-rays propagate through an object, the X-ray field undergoes a change in both amplitude and phase. The reduction in X-ray intensity after propagation through an object is a direct measure of the object’s absorption index, $\beta(u,v,w)$. The phase shift of the X-ray field that is induced by the object is a direct measure of its refractive index decrement, $\delta(u,v,w)$. Ignoring constant phase terms, the modulated X-ray field, $f_{O}(u,v)$, is expressed as,

where $f_{I}(u,v)$ is the complex-valued incident X-ray field, $T(u,v)$ is the transmission function such that $T(u,v)=\exp\left\{-A(u,v)-i\phi(u,v)\right\}$, where

Here, $\lambda$ is the wavelength of the monochromatic and coherent X-ray field.

As the X-ray field propagates from the object to the detector, the phase shifts induced by the object manifest as Fresnel diffraction fringes in the X-ray intensity images that are measured at the detector. The X-ray field at the detector plane, $f_{D}(u,v)$, is expressed as a 2D convolution of the X-ray field at the exit plane, $f_{O}(u,v)$, downstream of the sample and the Fresnel impulse response function, $\frac{\exp(ikR)}{i\lambda R}\exp\left\{\frac{i\pi}{\lambda R}\left(u^{2}+v^{2}\right)\right\}$, where $R$ is the object to detector distance and $k=\frac{2\pi}{\lambda}$ is the wavenumber [15, 48, 32]. The X-ray field propagation is depicted in Fig. 1.

The Fresnel transform is more efficiently computed in Fourier frequency space. Let $F_{D}(\mu,\nu)$ and $F_{O}(\mu,\nu)$ denote the 2D Fourier transform of the X-ray fields $f_{D}(u,v)$ and $f_{O}(u,v)$ respectively, where $\left(\mu,\nu\right)$ are the 2D Fourier frequency coordinates. Then,

Information on the phase of the X-ray field is lost since detector measurements are only sensitive to the intensity of the X-ray field. Thus, the detector measurements are modeled as, $|f_{D}(j\Delta,k\Delta)|^{2}$, where $\Delta$ is the width of each detector pixel, $|\cdot|$ denotes the magnitude of a complex number, and $(j,k)$ denote the row and column indices of the detector pixel. The square-root of the normalized detector measurement is then modeled as,

where $f_{I}(u,v)$ is the incident X-ray field from equation (1).

We derive the discretized relation between the incident X-ray field $g_{I}(j,k)$ and the X-ray field $g_{O}(j,k)$ at the exit plane of the object. Given that $x(j,k)$ represents the transmission function in discrete space, we have,

Next, we discretize the relation between the square-root normalized detector measurements $y(j,k)$ in equation (4) and the X-ray field $g_{O}(j,k)$ at the exit plane of the object. To sample the Fresnel transform in equation (3), we substitute $\mu=p\Delta_{\mu}$ and $\nu=q\Delta_{\nu}$, where $(p,q)$ are the discrete frequency coordinates, $\Delta_{\mu}=\frac{1}{N_{u}\Delta}$ is the sampling width along the $\mu$-axis, and $\Delta_{\nu}=\frac{1}{N_{v}\Delta}$ is the sampling width along the $\nu$-axis. Note that $N_{u}$ and $N_{v}$ are the number of discrete coordinates along the $u$-axis and $v$-axis respectively. Thus, a discrete sampling of Fourier space Fresnel transform is given by,

Given equation (6), we can now express the relation between the DFT $G_{D}(p,q)$ of $g_{D}(i,j)$ and DFT $G_{O}(p,q)$ of $g_{O}(i,j)$ as,

We use edge padding in the space domain for $g_{O}(j,k)$ to avoid circular convolution artifacts. In edge padding, the pixel values along an edge are copied to the neighboring padded regions. Since the detector only measures the intensity of the X-ray field, the square root of the normalized measurement in the absence of noise is expressed as,

where $g_{D}(j,k)$ is the discrete X-ray field in the detector plane given by the inverse discrete Fourier transform (IDFT) of $G_{D}(p,q)$ and $|\cdot|$ denotes the magnitude of a complex number.

In practice, the noise in detector measurements is modeled using Poisson statistics. Due to the variance stabilizing property of the square root transformation of Poisson random variables [49], we can model the noise statistics of the square root normalized measurement $y(j,k)$ (experimentally determined using equation (25) in the supplementary document) as a Gaussian distribution with constant variance. Thus, the forward model for the square root normalized measurement is given by,

where $n(j,k)$ is additive Gaussian noise with a constant variance for all $j$ and $k$.

The forward model that relates $y(j,k)$ to $x(j,k)$ is obtained by combining equations (5), (7), and (9). For convenience of notation, we will use matrix-vector notation for mathematical formulation of the forward model in discrete space. In section III-B1, we present a forward model where $x(j,k)$ is unconstrained, which does not impose restrictions on the material composition of the imaged object. In this case, the $x(j,k)$ is a complex valued number that uniquely encodes the X-ray phase shift and total X-ray absorption by the object. In section III-B2, we present a forward model where $x(j,k)$ is constrained, which restricts the material composition of the imaged object. In this case, the feasible solution space of the complex valued $x(j,k)$ is restricted such that it can only be expressed as the function of a real valued $z(j,k)$. This constraint is useful to restrict the feasible solution space for the phase shift and X-ray absorption.

To reconstruct x without any constraints, we utilize X-ray images at several propagation distances. Let $\textbf{y}_{l}$ denote the square root normalized measurements at a propagation distance of $R_{l}$. Then, the negative log-likelihood function is $l(\textbf{y};\textbf{x})=\frac{1}{\sigma^{2}}\sum_{l=1}^{L}\left|\left|\textbf{y}_{l}-\left|\textbf{H}_{l}\textbf{x}\right|\right|\right|_{2}^{2}$, where $L$ is the total number of distances. Here, the noise in each element of $\textbf{y}_{l}$ is approximated to be additive Gaussian with variance $\sigma^{2}$ and $\left|\left|\cdot\right|\right|^{2}_{2}$ denotes the squared $l^{2}$ norm of a vector. Thus, the unconstrained x is estimated by solving the following optimization problem,

We call this approach as unconstrained NLPR (U-NLPR).

Using measurements at a single propagation distance and in the absence of sparsity constraints, it is difficult to independently reconstruct both the phase and absorption of the X-ray field after propagation through the imaged objects. Information on the phase and absorption are entangled in the X-ray intensity measurements. Hence, we impose restrictions on the composition of the object to halve the number of unknowns during phase-retrieval using equation (13).

Let $\textbf{y}_{l}$ denote the square root normalized measurements at a propagation distance of $R_{l}$. Then, the negative log-likelihood function is $l(\textbf{y};\textbf{z})=\frac{1}{\sigma^{2}}\sum_{l=1}^{L}\left|\left|\textbf{y}_{l}-\left|\textbf{H}_{l}\textbf{z}^{\alpha+i\gamma}\right|\right|\right|_{2}^{2}$. We estimate z by solving the following optimization problem,

The complex valued transmission function is then estimated as $\hat{\textbf{x}}=\hat{\textbf{z}}^{\alpha+i\gamma}$. This method is called constrained NLPR (C-NLPR).

The scalar constraint parameters of $\alpha$ and $\gamma$ are dependent on the material composition. They also influence the speed of convergence of C-NLPR. Hence, it is important to intelligently set the parameters of $\alpha$ and $\gamma$. For a pure-phase object with zero absorption, we can choose $\alpha=0$ and $\gamma=1$. However, reconstruction of pure-phase objects will not be investigated in this paper. If a sample is homogeneous, i.e., consists of a single material [17] or satisfies the phase-absorption proportionality constraint [21], then we have several choices for setting $\alpha$ and $\gamma$. Let $\delta$ and $\beta$ denote the scalar values of the refractive index decrement and absorption index under these constraints, i.e., we assume $\delta(u,v,w)/\beta(u,v,w)=\delta/\beta\,\,\forall u,v,w$.

We implemented our NLPR algorithms using python programming language and PyTorch framework [51]. PyTorch uses algorithmic differentiation (or automatic differentiation) to compute gradients of the objective functions that are used for minimization. In equation (15), algorithmic differentiation is used to compute the gradient of $l(\textbf{y};\textbf{x})$ with respect to x. Similarly, algorithmic differentiation is used to compute the gradient of $l(\textbf{y};\textbf{z})$ with respect to z in equation (16).

We use the LBFGS algorithm [52, 53] to solve the optimization problems in equations (15) and (16). The gradients computed using algorithmic differentiation are used by the LBFGS optimization algorithm to reconstruct x in (15) and z in (16). In this paper, we use the LBFGS implementation in [54] with a history size of $64$. The maximum number of iterations is capped at $10^{4}$ and we use Wolfe line search for optimal selection of step-size. We use a convergence criteria that automatically stops the LBFGS iterations based on the convergence of reconstruction and objective function values. We stop LBFGS when the following conditions are met for $M$ consecutive iterations,

• •

  Average of the absolute differences in the reconstruction (z in C-NLPR or x in U-NLPR), expressed as a percentage, is less than $L_{r}$. The difference is computed between the reconstructions at consecutive two iterations.

• •

  The absolute difference in the objective function ($l\left(\textbf{y};\textbf{z}\right)$ in C-NLPR or $l\left(\textbf{y};\textbf{x}\right)$ in U-NLPR), expressed as a percentage, is less than $L_{c}$.

In this paper, we choose $M=5$, $L_{c}=1\%$, and $L_{r}=0.5\%$ for all the simulated and experimental results in section VI. We recognize that this particular setting for the convergence criteria may appear overly stringent. Our particular choice was designed to ensure sufficient convergence for all data presented in this paper. However, for any given data set, the convergence criteria may be relaxed for reduced run-time.

To perform tomographic reconstruction, we need to estimate discrete sampled representations of $A(u,v)$ and $\phi(u,v)$ (from equation (2)). Let the vectors A and $\boldsymbol{\phi}$ be the discrete representations of $A(u,v)$ and $\phi(u,v)$ in raster order.

The reconstruction of refractive index decrement produced by U-NLPR and FBP may contain low frequency artifacts. These artifacts have been well-documented in the research literature [15, 59, 60, 31] on multi-distance phase-retrieval algorithms. Low frequency artifacts refer to slowly varying spurious artifacts in the reconstructions due to the loss of low frequency information in the measurements. Multi-distance phase-retrieval does not use material constraints such as the phase-absorption proportionality. Instead, they rely on measurements at a wide range of propagation distances for inversion.

We will investigate the loss of low frequency information using the transfer function in equation (3). Equation (3) expresses the X-ray field $F_{D}(\mu,\nu)$ at the detector as a function of the X-ray field at the exit-plane of the object $F_{O}(\mu,\nu)$ in Fourier space. In the limit as the frequency components $\mu$ and $\nu$ approach zero, we have,

Thus, the value for $F_{D}(\mu,\nu)$, in the limit of zero frequencies, does not change with the propagation distance $R$ since $R$ is a parameter of only the Fresnel transfer function in equation (3). The detector only measures X-ray intensities (equation (4)). To compensate for this loss in phase information, we acquire measurements at varying propagation distances $R$ such that the phase-contrast fringes vary in magnitude and thickness. The phase information that is encoded in the phase-contrast fringes can be retrieved by acquiring data at multiple distances. However, since $F_{D}(\mu,\nu)$ does not vary sufficiently with distance $R$ at the low frequencies, we obtain low frequency artifacts in the reconstruction due to insufficient information.

During phase-retrieval, we cannot reconstruct the average value of the phase irrespective of the number of multi-distance measurements. Let us represent the phase $\phi(u,v)$ as the sum of a constant $\phi_{0}$ and a zero-mean phase term $\tilde{\phi}(u,v)$, i.e.,

The constant phase term $\phi_{0}$ factors out as the scalar multiple $\exp\left\{-i\phi_{0}\right\}$ in equations (1) and (3). Since the square root normalized measurements $y(j,k)$ in equation (4) measure only the X-ray intensity, information on $\exp\left\{-i\phi_{0}\right\}$ is lost since $\exp\left\{-i\phi_{0}\right\}$ is a multiplying factor for $f_{D}(j\Delta,k\Delta)$.

Our approach to extracting quantitative information from the refractive index decrement reconstructions is to compare the reconstructed values for the material-of-interest with the reconstructed values of the background air. We assume that the theoretical refractive index decrement of the background air is $0$. Unfortunately, this process of quantitative evaluation cannot be automated since determination of the background air region is non-trivial and application dependent. First, we compute the average value of the reconstruction in the background, denoted as $\delta^{(bg)}$. Next, we compute the average value of the reconstruction inside the material-of-interest, denoted as $\delta^{(m)}$. Finally, we compute the background subtracted refractive index decrement of the material-of-interest as,

Background subtraction is necessary due to the loss of the average value for the refractive index in the reconstruction. After background subtraction, $\delta^{(m:dif)}$ is the quantitatively accurate refractive index decrement that is comparable to the theoretical values.

In this section, we compare the performance of our new NLPR algorithms against existing approaches using simulated phase-contrast CT data of a 3D object at a monochromatic X-ray energy of $20\,keV$. To avoid using the same forward model for both simulation and inversion, we simulate X-ray images from finely sampled objects at very high resolutions using equation (8). The simulated phase-contrast CT data, the stack of X-ray images over all views, is then sub-sampled to a size of $128\times 80\times 128$ by block-averaging the simulated images over non-overlapping square-shaped windows of pixels. Finally, we simulate Poisson-like noise in the measurements by adding zero-mean Gaussian noise with a standard deviation that is $0.1\%$ of the simulated data (equation (9)).

Cross-sectional slices through the simulated object along with the corresponding simulated X-ray images with phase-contrast are shown in Fig. 2. First, we simulate a single-material object (Fig. 2 (a, b)) that consists of several spheres with differing diameters but made using the same material. The material used for the spheres is SiC with a refractive index decrement ($\delta$) of $1.67\times 10^{-6}$ and absorption index ($\beta$) of $4.77\times 10^{-9}$ at an X-ray energy of $20\,keV$. For this object, the normalized X-ray image with phase-contrast at a propagation distance of $R=200\,mm$ is shown in Fig. 2 (e). Next, we simulate a multi-material object (Fig. 2 (c, d)) that consists of four spheres with varying diameters and differing atomic composition. For the spheres, we chose SiC ($\delta=1.67\times 10^{-6},\beta=4.77\times 10^{-9}$), Teflon ($\delta=1.1\times 10^{-6},\beta=9.09\times 10^{-10}$), Alumina ($\delta=2.03\times 10^{-6},\beta=3.97\times 10^{-9}$), and Polyimide ($\delta=7.61\times 10^{-7},\beta=3.21\times 10^{-10}$) as the materials. The X-ray images with phase-contrast for this heterogeneous object are shown in Fig. 2 (f-h) at propagation distances, $R$, of $10\,mm$, $200\,mm$, and $400\,mm$. For both the single and multi-material objects, the radii of the spheres were chosen to be $12\,\mu m$, $16\,\mu m$, $20\,\mu m$, and $24\,\mu m$ respectively. We simulate X-ray images of size $80\times 128$ with pixel width of $1.29\,\mu m$ at $128$ tomographic views equally spaced over an angular range of $180^{0}$. For phase-retrieval (PR), we use the square root of the normalized X-ray images (equation (25) in supplementary document).

A quantitative comparison of the refractive index decrement reconstructions using various multi-distance phase-retrieval algorithms is presented in Fig. 3. At each tomographic view, we use phase-retrieval to estimate the phase images that are a measure of the phase shift induced by the object on the X-ray field. Then, we use filtered back projection (FBP) to perform a tomographic reconstruction of the refractive index decrement from the phase images. Fig. 3 (a) is the normalized root mean squared error (NRMSE)⁷⁷ 7 NRMSE is normalized using the averaged $l^{2}$ norm of the ground-truth. between the reconstructions and the ground-truth. Fig. 3 (b) is the structural similarity index measure (SSIM)⁸⁸ 8 For SSIM, the inputs are linearly scaled such that the minimum and maximum values of ground-truth are mapped to $-1$ and $1$ respectively. We use Gaussian weights with standard deviation of $8$ pixels.. We use the NRMSE and SSIM implementations from scikit-image [58]. We compute NRMSE and SSIM for the entire volume of the spheres after background subtraction as described in the section V-B. For the conventional phase-retrieval methods [15, 28, 29] of CTF, TIE, and Mixed, we plot (dashed lines) the NRMSE and SSIM measures as a function of its regularization parameter. To adequately reduce NRMSE and increase SSIM using TIE phase-retrieval, we must carefully tune its regularization hyper-parameter to lie within a narrow range of parameter values. An arbitrary choice of either a very low or high parameter value for regularization will lead to sub-optimal results with TIE. For CTF and Mixed phase-retrieval, it is possible to achieve low NRMSE and high SSIM by choosing a sufficiently small value for the regularization. As the regularization parameter is reduced, the performance of CTF does not degrade while Mixed leads to a marginal reduction in performance. For CTF, we can minimize the NRMSE and maximize SSIM by choosing a sufficiently small fixed value for the regularization parameter. In particular, we fix the regularization parameter value using equation (19) to avoid the need for parameter tuning.

We use U-NLPR from section IV-A for multi-distance phase-retrieval without any material constraints. The performance of U-NLPR is influenced by the estimates of the phase and absorption that is used for initialization of the optimization in equation (15). In Fig. 3, we investigate the use of phase images from CTF, TIE, or Mixed phase-retrieval as initial estimates for initialization of U-NLPR. The performance of these conventional methods is dependent on the chosen value of the regularization parameter. In Fig. 3, we plot (solid lines) the NRMSE and SSIM for U-NLPR as a function of the regularization parameter of the conventional method (i.e., CTF, TIE, and Mixed) used for initialization. We observe that U-NLPR always produces lower NRMSE and higher SSIM when compared to the estimates from the method used for initialization. Importantly, we observe that the best performance of U-NLPR is achieved by initialization using the CTF algorithm with a sufficiently low value for the regularization as specified in equation (19). While the Mixed approach may also be suitable for initialization, we prefer the CTF approach due to the ease of choosing the regularization using equation (19). Note that the SSIM with Mixed PR decreases slightly at very low regularization values.

The reconstructions at the optimal regularization with the highest SSIM using the conventional phase-retrieval algorithms of TIE, CTF, and Mixed approaches are shown in Fig. 4 (a-c). While such an optimal selection is perhaps an unfair advantage to the conventional algorithms, we did this to compare against the best possible reconstructions. In Fig. 4 (d), we show the reconstruction using U-NLPR that is initialized with zero values for the phase and absorption. In Fig. 4 (e), we show the U-NLPR reconstruction that is initialized using CTF phase-retrieval that used the regularization of equation (19). Since U-NLPR does not use any regularization, we avoid the need for the tedious manual tuning of regularization hyper-parameters. From Fig. 4 (a-c) and Fig. 2 (c), we observe that conventional phase-retrieval methods produce spurious streak artifacts. While the artifacts are faint in the case of TIE, both CTF and Mixed produce strong streak artifacts. When initialized with zeros for the phase and absorption, U-NLPR in Fig. 4 (d) produce a better reconstruction with significantly fewer streak artifacts than the conventional methods. When initialized using CTF with the fixed regularization of equation (19), U-NLPR produces the best reconstruction with reduced artifacts as shown in Fig. 4 (e).

We also compared our U-NLPR algorithm with a composite approach that uses TIE phase-retrieval followed by a Gerchberg-Saxton phase-retrieval (GSPR) algorithm in Fig. 5. While GSPR [50, 61] is primarily used for far-field diffraction imaging, it has also been successfully applied for imaging in the Fresnel region. Fig. 5 (a) is the refractive index reconstruction using GSPR and TIE initialization. We observe artifacts that are similar to Fresnel diffraction fringes in Fig. 5 (a). U-NLPR with CTF initialization at the regularization from equation (19) has higher SSIM than GSPR at any regularization.

We visually compare the tomographic reconstruction performance of single-distance phase-retrieval algorithms in Fig. 6. Qualitative comparisons of phase-retrieval methods for single-material object and multi-material objects are shown in Fig. 6 (a-c) and Fig. 6 (d, e) respectively. Paganin phase-retrieval [17] produces streak artifacts as shown in Fig. 6 (a, d). We investigate the performance of C-NLPR from section IV-B for non-linear phase-retrieval using the single-material constraint. If C-NLPR is initialized with zeros for the phase images, C-NLPR reduces artifacts but significantly enhances the noise in Fig. 6 (b). C-NLPR using Paganin phase-retrieval for initialization provides the best quality reconstruction as shown in Fig. 6 (c, e). We use the C-NLPR/One-$\alpha$ method from section IV-B1 in Fig. 6. Surprisingly, we also observe that the qualitative performance does not degrade for the multi-material object.

Fig. 7 shows a quantitative analysis of the reconstruction performance for various single-distance phase-retrieval algorithms. With Paganin phase-retrieval, it is common to achieve sharper reconstructions by artificially lowering the propagation distance, $R$, that is input to the algorithm. However, from Fig. 7, we see that an inaccurate setting for $R$ compared to its true value leads to sub-optimal NRMSE and SSIM values. The best performance for Paganin is achieved when $R$ is set equal to the true propagation distance of $200\,mm$. Initializing C-NLPR with Paganin phase-retrieved images results in the lowest NRMSE and highest SSIM among all approaches. However, initializing C-NLPR with zero phase images produces sub-optimal reconstructions with large amounts of noise as shown in Fig. 6 (b). In Fig. 6 and Fig. 7, we used the C-NLPR/One-$\alpha$ method described in section IV-B1.

The convergence speed and reconstruction quality of C-NLPR are strongly influenced by the choice of $\alpha$ and $\gamma$ as evidenced in Fig. 8. Fig. 8 (a) compares the NRMSE and SSIM for different constraint choices of C-NLPR/TrOpt-$\alpha,\gamma$ (section IV-B3), C-NLPR/One-$\alpha$ (section IV-B1), and C-NLPR/One-$\gamma$ (section IV-B2). The rows of Fig. 8 (a) correspond to different choices for the material properties of $\delta,\beta$ while the columns cycle through the various $\alpha,\gamma$ settings used for imposing the single-material constraint. The row label “True-$\delta,\beta$” indicates that the simulated object used the true $\delta,\beta$ values for SiC at $20\,keV$. The label “High-$\delta$” indicates a very large value for $\delta$ that is $10\times$ higher than the $\delta$ of SiC. The label “Low-$\beta$” indicates a very low value for $\beta$ that is $1/50\times$ the $\beta$ of SiC. For the various materials characterized by its $\delta,\beta$ values, we see that C-NLPR/TrOpt-$\alpha,\gamma$ consistently produces the lowest NRMSE. Using the SSIM measure, the results are a tie between C-NLPR/TrOpt-$\alpha,\gamma$ and C-NLPR/One-$\alpha$. However, the latter only uses knowledge of the ratio $\delta/\beta$.

In Fig. 8 (b-d), we plot the number of iterations for the LBFGS optimization algorithm as a function of the tomographic view index. Fig. 8 (b), Fig. 8 (c), and Fig. 8 (d) are plots of the number of iterations (to meet our convergence criteria) for “True-$\delta,\beta$”, “High-$\delta$”, and “Low-$\beta$” respectively. For the true values of $\delta,\beta$ of SiC, C-NLPR/One-$\gamma$ has the fastest convergence in Fig. 8 (b) but increased NRMSE in the first row of Fig. 8 (a). In contrast, both C-NLPR/TrOpt-$\alpha,\gamma$ and C-NLPR/One-$\alpha$ converge slower but achieve lower NRMSE. From Fig. 8 (c) and second row of Fig. 8 (a), we see that numerical instabilities associated with C-NLPR/One-$\gamma$ result in slow convergence and corrupted reconstructions due to the high $\delta$ value (see section IV-B2). From Fig. 8 (c), C-NLPR/One-$\alpha$ is slow to converge and has a high NRMSE as indicated in the third row of Fig. 8 (a) (see section IV-B1). The maximum number of iterations was fixed at $10^{4}$ for the LBFGS optimization. While C-NLPR/TrOpt-$\alpha,\gamma$ has the best convergence for all cases, it requires knowledge of both $\delta$ and $\beta$. Alternatively, C-NLPR/One-$\alpha$ only uses knowledge of the ratio $\delta/\beta$ while also matching C-NLPR/TrOpt-$\alpha,\gamma$ in convergence speed for “True-$\delta,\beta$” and “High-$\delta$”. Our recommendation is to use C-NLPR/One-$\alpha$ since it achieves the best overall performance while only requiring knowledge of the ratio $\delta/\beta$. Here, our objective is primarily to select feasible $\alpha$ and $\gamma$ values that impact the convergence speed of the algorithm. Thus, small errors in the ratio of $\delta/\beta$ or the individual values of $\delta,\beta$ are inconsequential. Only an order-of-magnitude change in $\alpha,\gamma$ will have a noticeable impact on convergence.

For comparison of phase-retrieval methods using experimental data, a bundle of fibers comprising PolyEthylene Terephthalate (PET), PolyPropylene (PP), Aluminum ($Al$), and Aluminum Oxide ($Al_{2}O_{3}$) with respective diameters of $200\,\mu m$, $28\,\mu m$, $125\,\mu m$, and $20\,\mu m$ were enclosed in a borosilicate capillary glass and scanned for tomography at the X-ray imaging beamline 8.3.2 of the Advanced Light Source, Berkeley, California. Monochromator was set to deliver X-ray beam with energy of $22\,keV$. A PCO edge camera combined with $10\times$ lens was used to collect X-ray images, resulting in an effective pixel size of approximately $0.65\,\mu m$. A total number of $1312$ X-ray images were recorded over $180^{\circ}$ range and using a $500\,ms$ exposure time. The object was placed at approximately $21\,m$ from the X-ray source and the detector was placed at object-to-detector propagation distances of $30\,mm$, $100\,mm$, and $250\,mm$ for each scan. The size of each X-ray image was $1686\times 2532$. The normalized X-ray images are shown in Fig. S6 of the supplementary document.

Tomographic reconstructions of the refractive index decrement from phase images reconstructed using various multi-distance phase-retrieval algorithms are shown in Fig. 9. The regularization parameters for the conventional phase-retrieval methods of CTF, TIE, and Mixed are manually tuned to obtain the best visual quality of reconstructions. Substantial low frequency reconstruction artifacts are observed in Fig. 9 (a, b, g, l) for CTF and TIE phase-retrieval. Mixed phase-retrieval produces large amounts of reconstruction noise as shown in the zoomed images of Fig. 9 (h, m). We note that the noise with Mixed can be reduced by adjusting the regularization but it also severely reduces the quality of reconstruction. U-NLPR produces the best reconstructions in Fig. 9 (d, e, i, j, n, o) that substantially reduces both low frequency artifacts and noise when compared to the conventional methods. Surprisingly, U-NLPR with zero-initialization results in a similar level of reconstruction quality as U-NLPR that is initialized using CTF with regularization from equation (19). Note that some of the artifacts in the centers of Fig. 9 (a-e) are ring artifacts [62] and we do not explore the correction of ring artifacts in this paper. The complicated interaction between ring-artifact removal and phase-retrieval algorithms necessitate further investigation. Fig. S7 in the supplementary document presents a quantitative comparison of the phase-retrieval performance.

For comparison of single-distance phase-retrieval algorithms, we obtained experimental X-ray CT data with phase-contrast from the Advanced Light Source (ALS) Beamline 8.3.2. CT data of SiC fibers was acquired at an object-to-detector distance of $R=98\,mm$ and X-ray energy of $20\,keV$. The size of each X-ray image was $324\times 320$ and the pixel width was $0.645\,\mu m$. X-ray images were acquired at $256$ different views equally spaced over an angular range of $180$ degrees. A normalized X-ray image is shown in Fig. S8 of the supplementary document.

Tomographic reconstructions of the refractive index decrement from phase images reconstructed using various single-distance phase-retrieval algorithms are shown in Fig. 10. Fig. 10 (a, e) and (b, f) show the reconstructions using Paganin and C-NLPR/One-$\alpha$ respectively. The line profiles in Fig. 10 (c, d) are along the yellow lines in Fig. 10 (a, b). They highlight the reduction in streak artifacts using C-NLPR/One-$\alpha$ when compared to Paganin phase-retrieval. The zoomed images in Fig. 10 (e, f) also demonstrate the significant enhancement in sharpness using C-NLPR/One-$\alpha$ even in the presence of a second material that is not modeled by the constraints of C-NLPR. Fig. 10 (g) is a plot of the modulation transfer function (MTF) [63] that indicates sharper reconstructions from C-NLPR/One-$\alpha$ and C-NLPR/TrOpt-$\alpha,\gamma$. A convergence analysis for the reconstructions in Fig. 10 is shown in Fig. 11. Fig. 11 (a) is a plot of the time elapsed and the objective function as a function of the LBFGS iterations. Fig. 11 (b) is a plot of the percentage change in the objective function and the estimated values as a function of the iteration number. The total time for phase-retrieval at one view was $6.3$ seconds on a NVIDIA Tesla V100 GPU.

Bright-field (a.k.a. flat-field) and dark-field measurements are acquired to appropriately normalize the CT scans. Bright-field refers to measurements made with the X-ray beam but without the object and dark-field refers to measurements made in the absence of the X-ray beam. Typically, bright-field measurements are made at the same propagation distances as used for the phase-contrast CT scans.

For normalization, we make certain simplifying approximations on $f_{I}(u,v)$ that also serve to simplify the subsequent problem of phase-retrieval. First, we approximate the incident field $f_{I}(u,v)$ as a plane wave with constant phase. Without loss of generality, we assume that this constant phase is zero since any information on constant phase terms is lost (equation (4)). Next, we normalize the detector measurements using the bright- and dark-fields. Lastly, we compute the square root of the normalized detector image before running the algorithms presented in section IV. Let $b(j,k)$ and $d(j,k)$ denote the bright-field and dark-field measurements respectively. If $\tilde{y}(j,k)$ denotes the detector image with the object, the square root of the normalized detector image is given by,

Here, $(\tilde{y}(j,k)-d(j,k))$ is a measure of $\left|f_{D}(j\Delta,k\Delta)\right|^{2}$ and $(b(j,k)-d(j,k))$ is a measure of $\left|f_{I}(j\Delta,k\Delta)\right|^{2}$.