Maximum Likelihood Based Phase-Retrieval Using Fresnel Propagation Forward Models With Optional Constraints

X-ray phase-contrast tomography (XPCT) is widely used for high contrast 3D imaging using either synchrotron or laboratory microfocus X-ray sources. XPCT enables an order of magnitude improvement in image contrast of the reconstructed material interfaces with low X-ray absorption contrast. The dominant approaches to 3D reconstruction using XPCT relies on the use of phase-retrieval algorithms that make one or more limiting approximations for the experimental configuration and material properties. Since many experimental scenarios violate such approximations, the resulting reconstructions contain blur, artifacts, or other quantitative inaccuracies. Our solution to this problem is to formulate new iterative non-linear phase-retrieval (NLPR) algorithms that avoid such limiting approximations. Compared to the widely used state-of-the-art approaches, we show that our proposed algorithms result in sharp and quantitatively accurate reconstruction with reduced artifacts. Unlike existing NLPR algorithms, our approaches avoid the laborious manual tuning of regularization hyper-parameters while still achieving the stated goals. As an alternative to regularization, we propose explicit constraints on the material properties to constrain the solution space and solve the phase-retrieval problem. These constraints are easily user-configurable since they follow directly from the imaged object's dimensions and material properties.


I. INTRODUCTION
Propagation-based X-ray phase-contrast tomography (XPCT) at synchrotron beamlines is widely used for 3D imaging due to the high-intensity, monochromatic, spatially coherent, and parallel-beam properties of synchrotron X-rays.XPCT is also a popular tool for obtaining higher edge contrast in laboratory micro-focus X-ray CT systems.XPCT is used for 3D reconstruction of a wide variety of objects in biology [1]- [3], material science [4]- [8], medical imaging [9], [10], and paleontology [11], [12].In synchrotron XPCT 1 , the object is exposed to a parallel beam of X-rays and the X-ray intensity is recorded by a 2D detector at several rotation angles of the object (Fig. 1).Phase-contrast is a function of the Fresnel number and the phase shift induced by the object on the X-ray field.A detailed discussion of the parameters governing phase-contrast is provided in [13].
In our application (Fig. 1), we adjust the Fresnel number by modifying the propagation distance to achieve adequate phase-contrast in the X-ray images.XPCT can either be performed at a single object-to-detector distance or be repeated at several object-to-detector distances.
Fig. 1.Parallel-beam X-ray phase-contrast tomography (XPCT) experiment at a synchrotron user-facility.The 3D object is rotated along an axis and 2D detector measurements are periodically acquired at several rotation angles.The object-to-detector propagation distance, R, is adjusted to produce phasecontrast in the measured X-ray images.The labels for the X-ray fields are defined in section II.
The solution to the inverse problem of object reconstruction is a combination of two sequential steps.First, we perform 2D phase-retrieval at each tomographic view to reconstruct the complex-valued images of the transmission X-ray field, f O (u, v) in Fig. 1, or its phase shift from the detector images |f D (u, v)| 2 [14], [15].Henceforth, we shall refer to the phase component of the recovered field f O (u, v) as the phase image.Next, we use a tomographic reconstruction algorithm to reconstruct the 3D refractive index decrement of the object from the phase images [16].In our phase-retrieval algorithms, we reconstruct the 2D phase image at each tomographic view angle from X-ray images at one or more propagation distances.Phase-retrieval is a non-linear inverse problem since it relies on the inversion of a non-linear analytical forward model that relates the phase to the X-ray images.The tomographic reconstruction step is used to reconstruct the 3D refractive index decrement from the retrieved phase images at all the views.Reconstruction of the refractive index decrement from the phase images is a linear inverse problem that is solved using any analytical or iterative CT reconstruction method, including the widely used filtered back projection (FBP) algorithm [16].
The popular choice for phase-retrieval algorithms in XPCT relies on the inversion of forward models that express the mea-sured X-ray images 2 as a linear transformation of the phase image.For XPCT using X-ray images at a single propagation distance, phase-retrieval reduces to the application of a digital linear filter on the individual X-ray images [14], [17]- [24].These methods also constrain the material composition by approximating the absorption index to be zero or enforcing a proportionality relation between the refractive index decrement and the absorption index.In the case of multi-distance XPCT, phase-retrieval is performed by the inversion of approximate linear relations between the multi-distance X-ray images and the phase [15], [25]- [29].Some of these methods avoid the approximations of zero-absorption and phase-absorption proportionality [15], [28], [29].Gureyev et al. [30], [31] derives linear analytic formulas for phase-retrieval in the Fresnel region for coherent and partially-coherent radiation.While such linear phase-retrieval (LPR) algorithms are computationally fast, the forward modeling approximations severely restrict the material composition and experimental design.If LPR algorithms are used outside their range of validity, they produce reconstructions with quantitative inaccuracies, false (unwanted) artifacts, and/or image blur.
To address the limitations of linear phase-retrieval (LPR), several non-linear phase-retrieval (NLPR) algorithms have been proposed [32]- [40].The algorithms in [33]- [36], [38]- [40] reconstruct the phase images with regularization to promote sparsity in projection space.Alternatively, the algorithms in [32], [37] directly reconstruct the 3D object using tomography consistency conditions to improve phase-retrieval by leveraging information from across view angles.While the latter approach is arguably better to constrain the solution space, it is also highly compute intensive while also being difficult to parallel compute.Hence, we pursue the former disjoint approach of phase-retrieval followed by tomographic reconstruction of the object.
NLPR algorithms [32]- [40] are generally more accurate than LPR algorithms since they avoid the linear approximations made in the forward measurement model.NLPR estimates the phase images by minimizing a distance measure between the X-ray images and the output of a non-linear forward measurement model.They also use certain sparsity constraints in the form of prior models or regularization functions to limit the solution space for the phase-retrieval problem.The use of such explicit regularization functions requires the user to manually tune one or more parameters to achieve the claimed performance improvements of artifact reduction, improved sharpness, and quantitative accuracy.In contrast, our approach achieves these benefits while avoiding parameter tuning due to the absence of regularization functions.Optionally, our approach also supports implicit constraints in our non-linear forward model that follow directly from the material properties.
In recent years, several deep learning (DL) approaches [41]- [47] have also been proposed for phase-retrieval in XPCT.Using DL, the authors demonstrated improved performance without the excessive computational cost of NLPR approaches. 2 Or, simple functional forms such as affine and/or logarithm transforms of the X-ray images are expressed as a linear transform of the phase.
However, the limitation of DL approaches is the need for representative training data to train the neural network models.In contrast, both LPR and NLPR algorithms do not require any training data and are broadly applicable to a wide range of objects.
In this paper, we present new non-linear phase-retrieval (NLPR) algorithms for both single-distance and multi-distance XPCT.A limited version of this manuscript with preliminary results was published in the form of a conference proceedings paper [48].Our NLPR algorithms also support constraints on the material composition.In particular, we demonstrate the use of the phase-absorption proportionality constraint (mathematically equivalent to the single material constraint) for single-distance XPCT.First, we formulate discrete nonlinear forward models using the Fresnel transform in Fourier frequency space.This model expresses the measured data as an analytical non-linear function of the phase images.Next, we formulate an objective function that is a measure of the distance between the measured data and the forward model output.Unlike existing non-linear approaches, we do not use prior-models or regularization functions in the objective function to achieve the desired goals of artifact-free, sharp, and accurate reconstructions.Finally, we iteratively solve for the phase images that minimize the objective function while also satisfying any constraint on the material composition.
We initialize our NLPR algorithms with the phase and absorption images that are estimated using LPR algorithms.We perform a comprehensive investigation of initialization strategies by varying the LPR method used for initialization and its regularization parameter.The characteristics and contributions of our proposed NLPR algorithms are - • Maximum likelihood formulation of phase-retrieval by inversion of non-linear forward measurement models.• Fresnel propagation based non-linear forward models with optional constraints.• Phase-retrieval by minimization of non-convex objective functions using the LBFGS algorithm.• Quantitatively accurate and sharp reconstructions with reduced artifacts compared to LPR approaches.• Avoids manual tuning of regularization hyper-parameters due to the absence of regularizing prior models.• Our recommended initialization strategy also avoids hyper-parameter tuning.
• Open-source software with documentation.We implemented our algorithms using the PyTorch framework and python programming language.For faster computation, our algorithms can also be run on multiple GPUs.We have released our phase-retrieval software under an opensource license along with documentation at the link https: //github.com/phasetorch/phasetorch.

II. X-RAY MEASUREMENT PHYSICS
In this section, we present a mathematical formulation for the physics of measurement in XPCT.

A. Modulation of the X-ray Field
In XPCT, the measured data is sensitive to the 3D variations in the absorption index, β(u, v, w), and refractive index decrement, δ(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 − δ(u, v, w) + iβ(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, β(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, δ(u, v, w).Ignoring constant phase terms, the modulated Xray field, f O (u, v), is expressed as, where Here, λ is the wavelength of the monochromatic and coherent X-ray field.

B. Measurement of Fresnel Propagated 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, , where R is the object to detector distance and k = 2π λ is the wavenumber [15], [32], [48].The X-ray field propagation is depicted in Fig. 1.
The Fresnel transform is more efficiently computed in Fourier frequency space.Let F D (µ, ν) and F O (µ, ν) denote the 2D Fourier transform of the X-ray fields f D (u, v) and f O (u, v) respectively, where (µ, ν) are the 2D Fourier frequency coordinates.Then, (3) 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∆, k∆)| 2 , where ∆ is the width of each detector pixel, | • | 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).

III. FORWARD MODEL
The forward models formulated in this section will be used in the phase-retrieval algorithms that estimate the transmission function from the measurements y(j, k).To formulate a forward model, we first translate the continuous space expressions in section II to discrete space.Let x(j, k) denote a sampling of the transmission function T (u, v) in equation (1).Let g D (j, k), g O (j, k), and g I (j, k) denote the sampled discrete representation of the continuous space X-ray fields f D (u, v), f O (u, v), and f I (u, v) respectively 3 .Let G D (p, q), G O (p, q), and G I (p, q) represent the discrete Fourier transform (DFT) coefficients of g D (j, k), g O (j, k), and g I (j, k) respectively.

A. Discretization of Measurement Model
We derive the discretized relation between the incident Xray 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 µ = p∆ µ and ν = q∆ ν , where (p, q) are the discrete frequency coordinates, ∆ µ = 1 Nu∆ is the sampling width along the µ-axis, and ∆ ν = 1 Nv∆ is the sampling width along the ν-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 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 |•| 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.

B. Non-Linear Forward Models for XPCT
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 Xray 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.
1) Unconstrained x: Let y, x, and n be column vectors whose elements include a raster ordering of the discrete samples y(j, k), x(j, k), and n(j, k) respectively.The matrix H is an operator that when left-multiplied to x returns the Fresnel transform of x.The forward model term |Hx| computes the magnitude of the IDFT of the Fresnel transform (equation ( 7)) of the DFT of x(j, k) 4 .The vector form of the unconstrained forward model is, where |•| denotes element-wise magnitude.The DFT and IDFT operations are implemented using fast Fourier transform algorithms.Equation (10) expresses the dependence of the real valued measurement vector y on the complex valued transmission vector x.
2) Constrained x: In the absence of data y at multiple propagation distances, unique reconstruction of x is achieved by imposing constraints on the imaged sample.For singledistance XPCT, it is typical to use constraints such as single material [17], phase-absorption proportionality [21], or zero absorption [20].However, such constraints may introduce inaccuracies and artifacts for objects that do not satisfy these constraints.
The single-material and phase-absorption proportionality constraints are mathematically expressed as, Alternatively, under the zero-absorption constraint, we have, The two constraints in equations ( 11) and ( 12) are equivalent to expressing the complex valued x in terms of a real valued z such that, where α and γ are real valued constants that are dependent on the material composition.For the phase-absorption proportionality constraint (equation ( 11)), we set the ratio of γ over α equal to the ratio of the refractive index decrement δ to the absorption index β.For the zero-absorption constraint in equation ( 12), we set α = 0. We discuss the possible choices for α and γ and the associated trade-offs in section IV-B.Equation ( 13) halves the dimensionality of the phaseretrieval problem by expressing the complex valued vector x in equation ( 10) as a function of a real valued vector z.Given the constraint in equation ( 13), the forward model that expresses y in terms of z is given by,

IV. NON-LINEAR PHASE-RETRIEVAL (NLPR)
In this section, we will formulate phase-retrieval algorithms using the maximum likelihood (ML) estimation framework.The reconstruction x is such that it minimizes the negative log-likelihood function l(y; x), which is a measure of the statistical likelihood of the measurement data y given the transmission function x.Under the ML framework, we perform phase-retrieval by solving the optimization problem of x = arg min x l(y; x).First, we formulate an approach to estimate the unconstrained transmission function x.Next, we present an approach to estimate x such that it also satisfies the constraint in equation (13).Finally, we solve for the X-ray absorption and phase shift at the exit plane of the object from the estimated x.We do not use a prior likelihood model to enforce sparsity in reconstruction of x.Instead, we initialize x using conventional linear phase-retrieval algorithms, which mimics the role of regularization and leads to improved solutions [50].

A. Unconstrained NLPR (U-NLPR)
To reconstruct x without any constraints, we utilize Xray images at several propagation distances.Let y l denote the square root normalized measurements at a propagation distance of R l .Then, the negative log-likelihood function is l(y; x) = 1 (15) We call this approach as unconstrained NLPR (U-NLPR).

B. Constrained NLPR (C-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 y l denote the square root normalized measurements at a propagation distance of R l .Then, the negative loglikelihood function is l(y 2 .We estimate z by solving the following optimization problem, The complex valued transmission function is then estimated as x = ẑα+iγ .This method is called constrained NLPR (C-NLPR).
The scalar constraint parameters of α and γ 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 α and γ.For a pure-phase object with zero absorption, we can choose α = 0 and γ = 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 α and γ.Let δ and β denote the scalar values of the refractive index decrement and absorption index under these constraints, i.e., we assume δ(u, v, w)/β(u, v, w) = δ/β ∀u, v, w.
1) C-NLPR / One-α: For this constraint, we set α = 1 such that x = z 1+iγ from equation (13).From equations ( 2) and ( 13), we see that |x| = z is a discretization of the absolute value for the transmission function The dynamic range of z is determined by the corresponding dynamic ranges for A(u, v) that is in-turn determined by β(u, v, w).For low X-ray absorption materials, if z n denotes the n th element of z, then z n ≈ 1 since β ≈ 0. The precision of 32-bit floating point numbers becomes progressively inadequate to accurately represent the dynamic range of z n as it approaches 1.Hence, α = 1 may lead to very slow convergence or numerical instabilities at very low values of β.Since z is the discretized representation of exp {−A(u, v)}, z γ is a discretization of exp {−ϕ(u, v)} only when γ = δ/β.Hence, α = 1 and γ = δ/β leads to the parameterization of x = z 1+iδ/β , which only requires knowledge of the ratio of δ/β.This method will be referred to as C-NLPR/One-α.
2) C-NLPR / One-γ: For this constraint, we set γ = 1 such that x = z α+i from equation (13).From equations ( 2) and ( 13), we see that z i is a discretization of the phase component of the transmission function exp {−iϕ(u, v)}.The dynamic range of z is determined by the corresponding dynamic range of ϕ(u, v) or δ(u, v, w).For objects with a high refractive index decrement, if z n denotes the n th element of z, then z n ≈ 0 if δ is very large.Hence, γ = 1 may lead to very slow convergence or numerical instabilities for highly refractive materials.Since z is a discretization of exp {−ϕ(u, v)}, z α will represent exp {−A(u, v)} only when α = β/δ.Hence, γ = 1 and α = β/δ leads to the parameterization of x = z β/δ+i , which only requires knowledge of the ratio of δ/β.This method will be referred to as C-NLPR/One-γ.
3) C-NLPR / TrOpt-α, γ: For this constraint, we control the dynamic range of z depending on the individual known values of δ and β.We set α and γ such that z n approximately varies between a preset low value of T l and a high value of 1.For T l , we simply choose a value that is greater than 0 but also an order of magnitude less than 1.In this paper, we choose T l = 0.01.Then, the values for α and γ are such that, • z n = 1 in the absence of any material along the ray path for the n th pixel.• z n = T l when a material with refractive index decrement δ and absorption index β lies along the entirety of the n th ray.Ideally, z n should only vary between T l and 1 in the absence of measurement non-idealities and a perfect choice for δ and β.In practice, z n varies approximately within this chosen transmission range, which reduces the risk of numerical instabilities.Thus, our choice for α and γ are, where ∆ is the pixel width.The number of pixels in the X-ray images along the u−axis and v−axis are N u and N v respectively.However, this particular choice of α and γ requires approximate knowledge of both δ and β.In contrast, the previous two choices that are described in sections IV-B1 and IV-B2 only require knowledge of the ratio δ/β.The constraint obtained using equations ( 17) and ( 18) will be referred to as C-NLPR/TrOpt-α, γ (TrOpt indicates optimized transmission).
In section VI, we demonstrate that α = 1, γ = δ β provides the best trade-off between stable optimization and good performance while only requiring knowledge of the ratio δ/β.

C. Optimization Algorithm
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(y; x) with respect to x.Similarly, algorithmic differentiation is used to compute the gradient of l(y; 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 (y; z) in C-NLPR or l (y; x) 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.

D. Initialization
1) Multi-distance linear phase-retrieval: The performance of U-NLPR is dependent on the initial estimate for x that is used to initialize the optimization in equation (15).If the phase and absorption in equation ( 2) are set to zero, then each element of the vector x is 1.
For multi-distance phase-retrieval using U-NLPR, we show in section VI that zero-initialization for the phase/absorption leads to improved results compared to the state-of-the-art conventional approaches to multi-distance phase-retrieval.However, we also demonstrate that initialization using conventional phase-retrieval methods can further improve the performance of U-NLPR.In this paper, we explore the use of the following phase-retrieval (PR) methods for initialization of U-NLPR.
1) Contrast Transfer Function (CTF) PR: Derived by Taylor expansion of the X-ray transmission field under the assumptions of weak absorption and slowly varying phase shift [15], [26].

2) Transport of Intensity Equation (TIE) PR: Derived by
Taylor expansion of the X-ray transmission field under the assumption of small propagation distances [15], [55].3) Mixed PR: Extends the validity of the approximations in CTF and TIE by forming a new hybrid PR that combines CTF and TIE [15], [29].Unlike U-NLPR, these conventional phase-retrieval methods use a regularization parameter α ′ that must be fine-tuned for acceptable performance [15].The value of this regularization parameter α ′ will impact both the performance of the conventional phase-retrieval method used for initialization of U-NLPR and the U-NLPR algorithm.To achieve the best per-formance for U-NLPR without the need for parameter tuning, we use the Contrast Transfer Function (CTF) phase-retrieval [15] with a sufficiently low value for the regularization given by, where ν is a very small value and the term (BC − A 2 ) is defined 5 in the reference [15].In this paper, we set ν = 10 −8 .As a pre-processing step before running multi-distance phase-retrieval algorithms including U-NLPR, the X-ray images must be registered such that the object appears at the same location in the X-ray images at all the propagation distances.
2) Single-distance linear phase-retrieval: For singledistance phase-retrieval using C-NLPR, we investigate the use of zero initialization and Paganin PR [17] for the phase images.Paganin PR solves the transport of intensity equation while assuming a single-material object (phase-absorption proportionality).

E. Estimation of Phase and Absorption
To perform tomographic reconstruction, we need to estimate discrete sampled representations of A(u, v) and ϕ(u, v) (from equation ( 2)).Let the vectors A and ϕ be the discrete representations of A(u, v) and ϕ(u, v) in raster order.
1) U-NLPR: The X-ray absorption, Â, and phase shift, φ, that is induced by the object on the incident X-ray field is, where x(I) and x(R) are the imaginary and real parts of x respectively.Here, log (•) and tan −1 (•)6 are element-wise vector operators.The phase estimated using U-NLPR is wrapped if the dynamic range for the phase exceeds 2π.Hence, we use phase unwrapping [57], [58] to unwrap the phase images φ prior to tomographic reconstruction.
2) C-NLPR: The X-ray absorption and phase images for C-NLPR are computed as, Since the phase φ obtained using equation ( 21) is already unwrapped, we do not need to apply an explicit phase unwrapping procedure.Note that equations ( 20) and ( 21) are for a single angular view of the CT scan.Hence, we need to repeatedly apply equations (20) and (21) to compute the absorption and phase images for each view independently.

V. TOMOGRAPHIC RECONSTRUCTION
In XPCT, measurement data is acquired at several rotation angles of the object.Let ϕ (n) and A (n) denote the phase and absorption images at view index n.The φ(n) and Â(n) from equations ( 20) and ( 21) are proportional to the projections of the refractive index decrement and absorption index at view n.
From equation (2), we see that the phase shift and absorption terms divided by the wavenumber are linear projections of the refractive index decrement and absorption index respectively.Filtered back projection (FBP) [16] is a popular algorithm that is widely used for reconstruction of X-ray absorption index from its linear projections.Thus, we use FBP to also reconstruct the refractive index decrement from its linear projections, ϕ (n) divided by the wavenumber, at all the views.In this paper, we do not investigate reconstruction of the absorption index.

A. Low Frequency Information Loss
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], [31], [59], [60] 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.Multidistance 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 (µ, ν) at the detector as a function of the X-ray field at the exit-plane of the object F O (µ, ν) in Fourier space.In the limit as the frequency components µ and ν approach zero, we have, Thus, the value for F D (µ, ν), 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 (µ, ν) 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 ϕ(u, v) as the sum of a constant ϕ 0 and a zero-mean phase term φ(u, v), i.e., The constant phase term ϕ 0 factors out as the scalar multiple exp {−iϕ 0 } 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 {−iϕ 0 } is lost since exp {−iϕ 0 } is a multiplying factor for f D (j∆, k∆).

B. Quantitative Evaluation
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 δ (bg) .Next, we compute the average value of the reconstruction inside the material-of-interest, denoted as δ (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, δ (m:dif ) is the quantitatively accurate refractive index decrement that is comparable to the theoretical values.

A. Simulated Data
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 Xray 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 × 80 × 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 phasecontrast are shown in Fig. 2. First, we simulate a singlematerial 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 (δ) of 1.67×10 −6 and absorption index (β) of 4.77 × 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 (δ = 1.67×10 −6 , β = 4.77 × 10 −9 ), Teflon (δ = 1.1 × 10 −6 , β = 9.09 × 10 −10 ), Alumina (δ = 2.03 × 10 −6 , β = 3.97 × 10 −9 ), and Polyimide (δ = 7.61 × 10 −7 , β = 3.21 × 10 −10 ) as the materials.The Xray 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 multimaterial objects, the radii of the spheres were chosen to be  19).Irrespective of the regularization, U-NLPR produces lower NRMSE and higher SSIM when compared to the phase-retrieval method used for initialization.The point of intersection of the vertical orange dotted line with the solid red line indicates the performance for U-NLPR with CTF initialization using the regularization from equation (19), which achieves the best performance while avoiding manual tuning of the regularization hyper-parameter.
12 µm, 16 µm, 20 µm, and 24 µm respectively.We simulate X-ray images of size 80 × 128 with pixel width of 1.29 µ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 Multi-distance Phase-Retrieval (PR) followed by FBP Reconstruction of the Multi-Material Object  tomographic reconstruction of the refractive index decrement from the phase images.Fig. 3 (a) is the normalized root mean squared error (NRMSE) 7 between the reconstructions and the ground-truth.Fig. 3 (b) is the structural similarity index measure (SSIM) 8 .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 phaseretrieval 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 7 NRMSE is normalized using the averaged l 2 norm of the ground-truth. 8For 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.
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  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 hyperparameters.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 singlematerial 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-α method from section IV-B1 in Fig. 6. Surprisingly, we also observe that the qualitative performance does not degrade for the multimaterial 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 suboptimal reconstructions with large amounts of noise as shown in Fig. 6 (b).In Fig. 6 and Fig. 7, we used the C-NLPR/One-α method described in section IV-B1.
The convergence speed and reconstruction quality of C-NLPR are strongly influenced by the choice of α and γ as evidenced in Fig. 8. Fig. 8    The maximum number of iterations was fixed at 10 4 for the LBFGS optimization.While C-NLPR/TrOpt-α, γ has the best convergence for all cases, it requires knowledge of both δ and β.Alternatively, C-NLPR/One-α only uses knowledge of the ratio δ/β while also matching C-NLPR/TrOpt-α, γ in convergence speed for "True-δ, β" and "High-δ".Our recommendation is to use C-NLPR/One-α since it achieves the best overall performance while only requiring knowledge of the ratio δ/β.Here, our objective is primarily to select feasible α and γ values that impact the convergence speed of the algorithm.Thus, small errors in the ratio of δ/β or the individual values of δ, β are inconsequential.Only an orderof-magnitude change in α, γ will have a noticeable impact on convergence.

B. Experimental Data
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 µm, 28 µm, 125 µm, and 20 µ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× lens was used to collect X-ray images, resulting in an effective pixel size of approximately 0.65 µm.A total number of 1312 X-ray images were recorded over 180 • 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 × 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 multidistance 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 phasecontrast from the Advanced Light Source (ALS) Beamline 8.3.2.CT data of SiC fibers was acquired at an object-todetector distance of R = 98 mm and X-ray energy of 20 keV .The size of each X-ray image was 324 × 320 and the pixel width was 0.645 µ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 singledistance phase-retrieval algorithms are shown in Fig. 10.Fig. 10 (a, e) and (b, f) show the reconstructions using Paganin and C-NLPR/One-α 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-α 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-α even in the presence of a second material that is not modeled by the constraints of C-NLPR.

VII. CONCLUSION
For propagation-based X-ray phase-contrast tomography (XPCT), we presented new non-linear phase-retrieval algorithms (NLPR) to reconstruct the phase shift induced by the object on the X-ray field.Then, we demonstrated reconstruction of the refractive index decrement in 3D from the phase shift images using tomographic reconstruction algorithms.Our approaches do not require any manual tuning of image quality related hyper-parameters such as regularization.Our NLPR algorithms are suitable for both single-distance and multidistance XPCT while also supporting constraints on the material composition.For single-distance XPCT, we demonstrated phase-retrieval (PR) under the constraint of single-material or phase-absorption proportionality.Our NLPR algorithms produced the best reconstructions based on both quantitative metrics of accuracy as well as qualitative evaluation of artifact and noise reduction.The superior performance of NLPR is a result of employing non-linear measurement models that are more accurate than the linear approximate models used by existing linear PR approaches.For multi-distance XPCT, we show that zero-initialization of the phase for NLPR may be sufficient, but the best performance is achieved by initializing with the Contrast Transfer Function (CTF) PR.For single-distance XPCT, we show that NLPR produces the best reconstruction when initialized with Paganin PR.

σ 2 L
l=1 ||y l − |H l x||| the unconstrained x is estimated by solving the following optimization problem,

Fig. 2 .Fig. 3 .
Fig. 2. Simulation of phase-contrast CT data.(a, b) and (c, d) show the simulated single-material homogeneous object and the multi-material heterogeneous object respectively.(a, c) show a planar slice (u − v axes) perpendicular to the rotation axis and (b, d) show a planar slice (u − w axes) parallel to the rotation axis.The slices in (a-d) pass through the center of the object.For the single-material object, the normalized X-ray image at a propagation distance, R, of 200 mm is shown in (e).For the multi-material object, the normalized X-ray images at propagation distances of 10 mm (F N = 2.68), 200 mm (F N = 0.13), and 400 mm (F N = 0.07) are shown in (f), (g), and (h) respectively.F N is the Fresnel number that is defined as F N = ∆ 2λR .All X-ray images in (e-h) are at the first tomographic view.Phase-contrast fringes are visible in (e), (g), and (h).From (f-h), we observe that increasing R also increases the strength of phase-contrast.

Fig. 4 .
Fig. 4. Tomographic reconstructions of the refractive index decrement from the phase images produced by the CTF, TIE, Mixed, and U-NLPR (proposed) phase-retrieval (PR) algorithms.(a-e) show planar slices along the u − v axes that pass through the center of the reconstruction volume.(a), (b), and (c) are using the TIE, CTF, and Mixed PR algorithms respectively.For the conventional PR methods of TIE, CTF, and Mixed, we present the best performing reconstruction at the regularization with the highest SSIM (from Fig. 3).(d) shows the reconstruction using U-NLPR that is initialized with zeros for the phase and absorption.(e) shows the reconstruction using U-NLPR that is initialized with CTF at the predetermined fixed regularization in equation (19).The gray values in (a-e) are scaled between −2.53 × 10 −7 and 2.03 × 10 −6 .Compared to the conventional PR reconstructions in (a-c), U-NLPR reduces streak artifacts and noise as shown in (d, e).Reconstruction slices along the u − w axes are shown in Fig S1 of the supplementary document.

Fig. 5 .
Fig. 5. (a) is the tomographic reconstruction (u − v axes) using Gerchberg-Saxton PR (GSPR) that is initialized with TIE PR at the best regularization parameter.(b) is the SSIM as a function of the regularization parameter for the PR that is used as initialization.GSPR produces artifacts that resemble Fresnel diffraction fringes in (a).U-NLPR with CTF initialization at the fixed regularization of equation (19) (intersection of the orange dotted line and the red solid line) has higher SSIM than GSPR at any regularization value.

Fig. 6 .
Fig. 6.Tomographic reconstructions of the refractive index decrement from the phase images produced by the single-distance phase-retrieval (PR) algorithms of Paganin and C-NLPR/One-α (proposed).(a-c) show reconstructions of the single material object.(d, e) show reconstructions of the multi-material object.(a-e) show planar slices along the u − v axes that pass through the center of the reconstruction volume.(a) and (d) show reconstructions using Paganin PR.(b) shows the reconstruction using C-NLPR that is initialized with zeros for the phase image (label C-NLPR/0-Initial).(c) and (e) show reconstructions using C-NLPR that is initialized using Paganin PR (label C-NLPR/Pag-Initial).The gray values in (a-e) are scaled between −2.09 × 10 −7 and 1.67 × 10 −6 .C-NLPR with Paganin initialization produces the best reconstruction that minimize noise and artifacts.

Fig. 7 .
Fig. 7. NRMSE and SSIM between reconstruction and ground-truth for the refractive index decrement as a function of the (inaccurate) propagation distance R that is input to Paganin phase-retrieval (PR).The true propagation distance R used for simulation is 200 mm.The NRMSE is lowest and SSIM is highest when R is set equal to its true value of 200 mm.Initializing C-NLPR with zeros for the phase results in sub-optimal performance that is worse than Paganin PR.C-NLPR initialized with Paganin PR (at the correct R) has the best performance.Here, C-NLPR refers to C-NLPR/One-α.

TrueFig. 8 .
Fig. 8. Quantification of performance and speed of C-NLPR for various choices of α and γ as described in sections IV-B1, IV-B2, and IV-B3.Our analysis is repeated for various δ, β values.True-δ, β uses the ground-truth δ, β of SiC.High-δ refers to a δ that is 10× the δ of SiC.Low-β refers to a β that is 0.02× the β of SiC.(a) shows the performance metrics of RMSE/SSIM within each square block.(a) uses a linear colormap where yellow indicates the best, green is better, blue is bad, and dark blue is the worst performance.(b, c, d) show the number of LBFGS iterations as a function of the view index.Corrupted reconstructions caused by numerical instabilities is indicated by "nan" (not-a-number) in (a) and shaded with a translucent color in (c).From (a, c), we see that C-NLPR/One-γ results in "nan" due to numerical instabilities for High-δ.Hence, we recommend avoiding C-NLPR/One-γ for large ratios of δ/β.With C-NLPR/One-α, while we achieve good overall performance in (a), the number of iterations may reach large numbers for Low-β as shown in (d).With C-NLPR/TrOpt-α, γ, we achieve good overall performance (see (a)) without the need for a large number of iterations (see (b, c, d)).Legend for (c, d) is same as in (b).

Fig. 9 .
Fig. 9. Experimental data reconstruction comparison of the refractive index decrement for various multi-distance phase-retrieval (PR) algorithms.(a-e) show planar reconstruction slices along the u − v axes passing through the center of the volume.The reconstruction in (a-e) is cropped to show the region within the interior of the sample holder.(f-j) and (k-o) zooms into two different regions of the reconstructions in (a-e).Since the quantitative values vary substantially between PR methods (see section V-A), we scale the gray values of each image individually between the value percentiles of 5% and 95%.Both CTF and TIE produce substantial low frequency artifacts in (a), (b), (g), and (l).Mixed reduces the artifacts but has increased noise as shown in (h, m).U-NLPR, irrespective of the initialization, produces the best reconstructions that minimize noise and artifacts.In particular, U-NLPR with CTF initialization (using the fixed regularization in equation (19)) produce less intense artifacts at the center of the images in (d, e) when compared to zero-initialization.Unlike U-NLPR, the CTF, TIE, and Mixed PR results are at optimal regularization values that were manually chosen to achieve the best visual quality of reconstructions.

Fig. 10 .
Fig. 10.Tomographic reconstruction of the refractive index decrement from phase images retrieved using single-distance phase-retrieval (PR) methods.(a) and (b) show the center slice (u − v axes) of the reconstruction using Paganin PR and C-NLPR/One-α respectively.(c) and (d) show line profile comparisons between Paganin and C-NLPR along the yellow marked lines in (a) and (b).In (c, d), the large dynamic range for Paganin and C-NLPR/One-γ indicate the presence of streak artifacts since these line-profiles are in the background region of (a, b).Both C-NLPR/TrOpt-α, γ and C-NLPR/One-α reduce streaking artifacts as illustrated in (c, d).(e, f) zooms into a region along a different slice in the presence of a second material between the SiC fibers, which demonstrates the sharper reconstruction using C-NLPR/One-α when compared to Paganin PR.Images in (a, b, e, f) are scaled between −1.79 × 10 −7 and 1.43 × 10 −6 .(c, d) use a common legend that is indicated at the top of the plots.(g) is the modulation transfer function (MTF) for the sharpness of the disc inside the red square in (a, b).(g) indicates sharper reconstructions using C-NLPR/One-α and C-NLPR/TrOpt-α, γ since the corresponding curves are above the curve for Paganin PR at the higher frequencies.Our recommendation is to use C-NLPR/One-α that only requires knowledge of δ/β.

Fig. 10 (
Fig.S7in 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 phasecontrast from the Advanced Light Source (ALS) Beamline 8.3.2.CT data of SiC fibers was acquired at an object-todetector distance of R = 98 mm and X-ray energy of 20 keV .The size of each X-ray image was 324 × 320 and the pixel width was 0.645 µ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.S8of the supplementary document.Tomographic reconstructions of the refractive index decrement from phase images reconstructed using various singledistance phase-retrieval algorithms are shown in Fig.10.Fig.10 (a, e) and (b, f) show the reconstructions using Paganin and C-NLPR/One-α 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-α 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-α 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α and C-NLPR/TrOpt-α, γ.A convergence analysis for the reconstructions in Fig.10is 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.
Fig. S4.Quantitative comparison of the refractive index reconstructions for the simulation experiment described in Fig. S2 and S3.(a) and (b) shows the NRMSE and SSIM as a function of the regularization parameter for the conventional PR.Since U-NLPR does not use regularization, the parameter along the horizontal axis is for the conventional PR methods used as initialization.U-NLPR consistently out-performs the conventional PR method used as initialization.U-NLPR with CTF initialization (U-NLPR/CTF-Initial) at the pre-determined regularization in equation (19) is the best result that also avoids parameter tuning (intersection of the vertical orange dotted line and the solid red line).

Fig. S6 .
Fig. S6.Normalized X-ray images at the first view angle from synchrotron phase-contrast CT of Al, Al 2 O 3 , P P , and P ET fibers.(a), (b), and (c) show the X-ray images at propagation distances of 30 mm, 100 mm, and 250 mm respectively.The intensity range for the gray-values is between 0.4 and 1.2.Phase-contrast fringes are stronger in (c) than in (a) due to the larger propagation distance of (c).

Fig. S7 .
Fig. S7.NRMSE for Al and PET fibers vs. regularization of conventional PR methods for the multi-distance experimental data.The dashed lines show the NRMSE for the conventional PR methods of TIE, CTF, and Mixed.The solid lines are the NRMSE plots for U-NPLR that is initialized using one of the conventional PR methods.While U-NLPR does not use regularization, its performance nevertheless varies with the regularization parameter of the conventional PR used for initialization.U-NLPR out-performs all the conventional PR methods irrespective of the initialization.U-NLPR with CTF initialization (regularization from equation (19) indicated by the dotted vertical orange line) has the best performance without the need for parameter tuning.

Fig
Fig. S8.(a) is the X-ray image of SiC fibers at the first tomographic view and propagation distance of R = 98 mm.(b) zooms into the image in (a) to better visualize the phase-contrast fringes.The intensity range of gray-values is between 0.53 and 1.97.