Constrained Non-Linear Phase Retrieval for Single Distance X-ray Phase Contrast Tomography
CConstrained Non-Linear Phase Retrieval for Single Distance X-ray Phase Contrast Tomography
K. Aditya Mohan ? , Dilworth Y. Parkinson † , Jefferson A. Cuadra ? ; ? Lawrence Livermore National Laboratory, Livermore, CA USA; † Lawrence Berkeley National Laboratory, Berkeley, CA USA;
Abstract
X-ray phase contrast tomography (XPCT) is widely used for3D imaging of objects with weak contrast in X-ray absorptionindex but strong contrast in refractive index decrement. To recon-struct an object imaged using XPCT, phase retrieval algorithmsare first used to estimate the X-ray phase projections, which isthe 2D projection of the refractive index decrement, at each view.Phase retrieval is followed by refractive index decrement recon-struction from the phase projections using an algorithm such asfiltered back projection (FBP). In practice, phase retrieval is mostcommonly solved by approximating it as a linear inverse problem.However, this linear approximation often results in artifacts andblurring when the conditions for the approximation are violated.In this paper, we formulate phase retrieval as a non-linearinverse problem, where we solve for the transmission function,which is the negative exponential of the projections, from XPCTmeasurements. We use a constraint to enforce proportionalitybetween phase and absorption projections. We do not use con-straints such as large Fresnel number, slowly varying phase, orBorn/Rytov approximations. Our approach also does not requireany regularization parameter tuning since there is no explicitsparsity enforcing regularization function. We validate the per-formance of our non-linear phase retrieval (NLPR) method usingboth simulated and real synchrotron datasets. We compare NLPRwith a popular linear phase retrieval (LPR) approach and showthat NLPR achieves sharper reconstructions with higher quanti-tative accuracy.
Introduction
Traditional X-ray computed tomography (XCT) relies on X-ray absorption contrast for imaging. However, many samples,e.g. animal soft tissue in biological imaging and material inter-faces with similar atomic number elements, have low absorptioncontrast. To solve this problem, phase contrast, which is typi-cally many orders of magnitude larger than absorption contrast,is induced by increasing the propagation distance from object todetector [1, 2] (Fig. 1). This technique where tomographic mea-surements with both absorption and phase contrast are acquiredat large object-to-detector propagation distances is called X-rayphase contrast imaging (XPCT).In XPCT, phase contrast in the measurements is related tothe X-ray phase, which is the 2D projection of the refractive in-dex decrement values of the object. In contrast, X-ray absorptionis quantified as the projection of the object’s absorption index.However, the non-linear dependence of X-ray measurements onthe absorption and phase projections makes it challenging to solvethe inverse problem of object reconstruction. Typically, a phaseretrieval algorithm is first used to estimate the phase projections
Figure 1.
System diagram of XPCT. In XPCT, the propagation distancefrom object to detector is increased to allow the phase variations in the X-ray field exiting the object to develop into intensity variations called phasecontrast at the detector. at each tomographic view. Then, the refractive index decrementvalues in 3D are reconstructed using an algorithm such as filteredback projection (FBP) [16].Several phase retrieval algorithms have been proposed forboth single distance and multi-distance XPCT [2–11]. Thesemethods either estimate the transmission function, which is thenegative exponential of the projections, or directly solve for theprojections. Among these methods, the most widely used ap-proaches for single distance phase retrieval use linear analyticalclosed form solutions to either recover the transmission functionor the projections [3–6]. A comparison between the various lin-ear phase retrieval approaches for single distance XPCT is in [2].Recently, non-linear approaches to single distance phase retrievalhave been presented in [1,10]. These non-linear approaches claimsubstantial improvements over linear approaches using non-linearforward models and sparsity enforcing regularization functions.However, various reasons including the need for regularizationparameter tuning and computational complexity have preventedthe widespread use of these methods.In this paper, we present a novel phase retrieval method thatuses a non-linear forward model and a constraint to enforce a lin-ear relation between the absorption and phase projections. Thelinear constraint can either be used to enforce zero values for theabsorption projection or enforce proportionality between absorp-tion and phase projections. The non-linear forward model usesFresnel transform to express the X-ray measurements at the detec-tor plane as a function of the transmission function. Our forwardmodel does not make any approximations made by linear phaseretrieval approaches such as large Fresnel number, slowly varyingphase, or Born/Rytov approximations [2].Our non-linear phase retrieval (NLPR) approach uses nu-merical optimization to estimate the transmission function suchthat the XPCT measurements best match its predictions from thenon-linear forward model. Since absorption and phase are pro-jection space quantities that are related non-linearly to the trans-
IS&T International Symposium on Electronic Imaging 2020Computational Imaging https://doi.org/10.2352/ISSN.2470-1173.2020.14.COIMG-146This is a work of the U.S. Government. ission function, we propose a novel formulation in transmissionspace that automatically enforces a linear relation between ab-sorption and phase projections during optimization. Since our ap-proach does not use any sparsity enforcing regularization functionto model prior information, significant manual effort and time of-ten spent in tuning the regularization parameter is saved. Further-more, since phase retrieval at each view can be run independently,our approach is easily parallelized on a multi-core processor orhigh performance computing (HPC) systems.
Measurement Physics
The physics of X-ray interaction with an object is defined interms of the object’s 3D variation in absorption index, b ( u , v , w ) ,and refractive index decrement, d ( u , v , w ) , where ( u , v , w ) rep-resent the 3D Cartesian coordinates. The absorption index, b ( u , v , w ) , modulates the amplitude (intensity) of the incomingX-ray field and the refractive index decrement, d ( u , v , w ) , modu-lates the X-ray phase. After propagation through the object, theeffective change in X-ray intensity is quantified by the absorptionprojection, A ( u , v ) = pl Z w b ( u , v , w ) dw , (1)and the change in X-ray phase is given by the phase projection, f ( u , v ) = pl Z w d ( u , v , w ) dw , (2)where X-ray propagation is assumed to be along the w axis and l is the X-ray wavelength.Next, we will express the X-ray field exiting the object, f O ( u , v ) , as a function of the X-ray field incident on the object, f I ( u , v ) (Fig. 1). The transformation from f I ( u , v ) to f O ( u , v ) isdependent on the transmission function T ( u , v ) defined as, T ( u , v ) = exp ( A ( u , v ) i f ( u , v )) , (3)where i = p T ( u , v ) C is complex valued, and A ( u , v ) R and f ( u , v ) R are real valued. Then, the X-ray field f O ( u , v ) exiting the object is given by, f O ( u , v ) = f I ( u , v ) T ( u , v ) . (4)The X-ray field f O ( u , v ) undergoes Fresnel diffraction as itpropagates towards the detector. If F D ( µ , n ) and F O ( µ , n ) denotethe Fourier transforms of f D ( u , v ) and f O ( u , v ) , respectively, then, F D ( µ , n ) = F O ( µ , n ) exp ⇣ i pl R ⇣ µ + n ⌘⌘ , (5)where we have ignored constant phase terms, R is the object-to-detector distance, and ( µ , n ) are the 2D Fourier frequency coor-dinates. However, each detector pixel is only able to measure theintensity of the incoming X-ray field. Hence, the measurement atdetector pixel ( j , k ) is modeled as,˜ y ( j , k ) = | f D ( j D , k D ) | , (6)where D is the pixel width of the detector and | · | denotes magni-tude. In our application, the phase and magnitude of the incidentX-ray field f I ( u , v ) is not known. The magnitude of f I ( u , v ) couldbe estimated if the bright-field (a.k.a. flat-field) measurementswere made at the object’s plane of X-ray incidence. However, thebright-field measurements are made when the object-to-detectordistance is R , which is the distance optimized for phase contrast.Hence, we make some simplifying assumptions on f I ( u , v ) thatalso simplifies the subsequent non-linear phase retrieval method.First, we assume that the incident X-ray field f I ( u , v ) is a planewave with constant phase. Hence, without loss of generality, weassume f I ( u , v ) has zero phase since any information on constantphase terms are lost when the detector makes intensity measure-ments. We also make an approximation where the effect of themagnitude of f I ( u , v ) is ignored by normalizing ˜ y ( j , k ) with itsbright and dark fields. Let b ( j , k ) denote the bright-field mea-surements acquired by the detector in the absence of the objectat a object-to-detector distance of R and d ( j , k ) denote the darkfield measurements acquired in the absence of X-rays. Let N u and N v denote the number of detector pixels along the u axis and v axis, respectively. Then, let the square root of the normalizeddetector image be given by, y ( j , k ) = s ˜ y ( j , k ) d ( j , k ) b ( j , k ) d ( j , k ) , (7)where 1 j N u and 1 k N v . Forward Model
In this section, we derive a forward model in discrete spacethat expresses the square root measurements y ( j , k ) in terms ofthe transmission function. Let the discrete space equivalent ofthe transmission function T ( u , v ) be denoted by ˜ x ( j , k ) . How-ever, there is no unique solution to the problem of estimating thecomplex valued ˜ x ( j , k ) from the real valued y ( j , k ) . To solve theuniqueness problem, we constrain the feasible solution space forthe transmission function ˜ x ( j , k ) . In particular, we assume that thethe complex valued ˜ x ( j , k ) can be represented in terms of a realvalued x ( j , k ) such that,˜ x ( j , k ) = x a + i g ( j , k ) , (8)where a and g are real valued constants used to constrain the do-main space of ˜ x ( j , k ) . Thus, the problem reduces to estimation ofthe real valued x ( j , k ) from y ( j , k ) .For the forward model, we use the discrete frequency equiv-alent of equation (5). Let z ( j , k ) be the discrete space equivalentof the continuous space X-ray field f D ( u , v ) at the detector plane.Let the discrete Fourier transforms (DFT) of z ( j , k ) and ˜ x ( j , k ) be denoted by Z ( p , q ) and ˜ X ( p , q ) , respectively, where ( p , q ) areDFT coordinates. Then, Z ( p , q ) = ˜ X ( p , q ) H ( p , q ) , (9)where H ( p , q ) is the Fresnel transform in discrete Fourier spacegiven by, H ( p , q ) = exp ⇣ i pl R ⇣ p D µ + q D n ⌘⌘ . (10)Here, D µ = N u D and D n = N v D are the widths of the Fourier fre-quency bins along the µ and n frequency axes (equation (5)), re-spectively. Then, z ( j , k ) is obtained by inverse discrete Fourier IS&T International Symposium on Electronic Imaging 2020Computational Imaging ransform (IDFT) of Z ( p , q ) . Thus, the forward model in discretespace is, y ( j , k ) = | z ( j , k ) | + w ( j , k ) , (11)where w ( j , k ) is additive Gaussian noise.For mathematical convenience, we will express the forwardmodel in vector form. Let y , z , x , and ˜ x be vectors in raster or-der corresponding to y ( j , k ) , z ( j , k ) , x ( j , k ) , and ˜ x ( j , k ) , respec-tively. The Fresnel transform in equation (9) including the DFTand IDFT operations are coded into a matrix H such that z = H ˜ x .Then, the forward model that expresses the measurement y interms of ˜ x = x a + i g is given by, y = Hx a + i g + w . (12)where | · | denotes element-wise magnitude and w is noise vectorconsisting of all w ( j , k ) in raster order. Phase Retrieval
Linear Phase Retrieval (LPR)
As a baseline, we use the linear phase retrieval (LPR) methodpresented in Paganin et al. [4]. This method assumes a homoge-neous object and solves for the transmission function as a linearanalytical transformation of the normalized measurements. It as-sumes data acquisition at large Fresnel numbers [2] which causesexcessive smoothing and artifacts at large object-to-detector dis-tances. The homogeneous assumption for the object is equivalentto assuming that the ratio of the refractive index decrement andthe absorption index is a constant given by db . In essence, thisphase retrieval method is a low pass filter where the amount ofsmoothing increases with the object-to-detector distance, R , andthe ratio db . Non-Linear Phase Retrieval (NLPR)
Our approach to non-linear phase retrieval (NLPR) is tosolve for x from square root measurements y using the non-linearforward model in equation (12). Thus, NLPR is formulated as,ˆ x = argmin x y Hx a + i g , (13)where x j R , H j , k C , y j R , | · | denotes element-wise mag-nitude, and ||·|| computes vector magnitude. Note that x j and y j denote the j th element of the vectors x and y , respectively. Simi-larly, H j , k denotes the element along the j th row and k th columnof the matrix H .To solve the uniqueness problem during phase retrieval, weimpose constraints to restrict the feasible solution space when op-timizing for the transmission function. In particular, we use aconstraint that forces the phase projection to be proportional tothe absorption projection [2]. The phase proportional to absorp-tion constraint equivalently assumes that the value at every pointinside the object has the same ratio of the refractive index decre-ment, d , and the absorption index, b . The phase proportional toabsorption constraint is obtained by setting a = g = db . Al-ternatively, the zero absorption constraint can also be used by asuitable setting of a and g . To solve equation (13), we use the limited memory Broy-den–Fletcher–Goldfarb–Shanno (LBFGS) algorithm with positiv-ity constraints [18, 19]. Since LBFGS uses the gradient infor-mation for optimization, we derive the gradients of the objectivefunction in equation (13) using algorithmic differentiation and thederivative tables in [12]. We use the LBFGS implementation inthe open source software library NLopt [17]. The estimated resultfrom LPR is used as an initial estimate when solving equation(13) using LBFGS.Once ˆ x is computed, the next step is to perform tomographicreconstruction of the refractive index decrement of the object.Note that ˆ x is a vector containing the estimates ˆ x ( j , k ) in rasterorder. From equation (2), we see that the phase projections are re-lated to the line integral of the refractive index decrement. Basedon equations (3) and (8), the phase projections in discrete spaceare computed as g ( log ( ˆ x ( j , k ))) . Once the phase projections arecomputed at each view, the refractive index decrement is recon-structed using filtered back projection (FBP). Results
Simulated Data (a) Single-material spheres (b) Multi-material spheres
Figure 2.
Simulated normalized measurements at first view angle. (a)and (b) show measurements for single material spheres phantom and multi-material spheres phantom, respectively.
In this section, we use simulated data to validate the per-formance of NLPR. We simulate a single material phantom anda multi-material phantom each containing three spheres of radii4 µ m , 6 µ m , and 5 µ m , respectively . The spheres in the singlematerial phantom have the same refractive index decrement, d , of1 . ⇥ and absorption index, b , of 4 . ⇥ . These d and b values correspond to the compound Silicon Carbide (SiC) at aX-ray energy of 20 keV . The three spheres of the multi-materialphantom have absorption index values of 4 . ⇥ (10 ⇥ larger), 4 . ⇥ , and 4 . ⇥ , respectively, while the re-fractive index decrement values are 1 . ⇥ ,1 . ⇥ , and3 . ⇥ (2 ⇥ larger), respectively.The shape of the ground-truth volume containing the sphereswas 96 ⇥ ⇥
128 with a voxel width of 0 . µ m . From theground-truth volume, we used equation (12) to simulate X-ray im-ages (radiographs) of size 48 ⇥
64 at a pixel width of 0 . µ m .Radiographs were simulated with Poisson noise at 64 differentviews equally spaced over an angular range of 180 . The object-to-detector distance was 100 mm and X-ray energy was 20 keV .Before any comparisons with the proposed method, the ground-truth volume was downsampled to a size of 48 ⇥ ⇥
64. The sim-ulated measurements from the single material and multi-material From left to right in each image of Fig. 2 and Fig. 3
IS&T International Symposium on Electronic Imaging 2020Computational Imaging etrieved phase and refractive index decrement reconstruction for single material spheres (a) Ground-Truth (b) LPR (c) LPR-Sharp (d) NLPR(e) Ground-Truth (f) LPR+FBP (g) LPR-Sharp+FBP (h) NLPR+FBP
Retrieved phase and refractive index decrement reconstruction for multi-material spheres (i) Ground-Truth (j) LPR (k) LPR-Sharp (l) NLPR(m) Ground-Truth (n) LPR+FBP (o) LPR-Sharp+FBP (p) NLPR+FBP
Figure 3.
Simulated data reconstruction comparison. (a-d) and (i-l) show the ground-truth and retrieved phases using various methods. (e-h) and (m-p) showthe ground-truth and reconstructed refractive index decrement using various methods. (a-h) and (i-p) compares results for single material spheres data andmulti-material spheres data, respectively. LPR refers to the method in Paganin et al. [4] that uses the true value of R = mm . LPR-Sharp is the method inPaganin et al. [4] but assuming an incorrect lower value of R = mm for improved sharpness. LPR results in excessive smoothing and streak artifacts as shownin (f,g,n,o) while NLPR produces an accurate artifact-free reconstruction of the object with sharp edges. Figure 4.
Line profile comparison of the reconstructed refractive indexdecrement through the center of the bottom two circles (spherical cross-sections) in Fig. 3 (e,f,h). phantoms at the first view angle are shown in Fig. 2. Beforerunning LPR and NLPR, each X-ray image was padded to 1 . ⇥ the original size in each dimension using edge padding. For the LBFGS optimization in NLPR, we used a value of 10 for therelative threshold on the optimization parameters, x , in NLopt[17]. In equation (13), we used a = g = R = mm . Alternatively, LPR-Sharp achievessharper images by using the approach in Paganin et al. [4] butwith the lower and incorrect object-to-detector distance of R = mm . From Fig. 3 (e-h) and Fig. 3 (m-p), we can see that LPRresults in blurred reconstructions with streak artifacts while NLPRproduces accurate reconstructions without any artifacts. Also, adrastic lowering of R from its true value of 100 mm still resulted in IS&T International Symposium on Electronic Imaging 2020Computational Imaging uantifying performance for single material spheres
GT LPR+FBP LPR-Sharp+FBP NLPR+FBP . µ m . µ m . µ m . µ m . µ m . µ m . µ m . µ m . µ m . µ m . µ m . µ m (a) Areas of circles (spherical cross-sections) in Fig. 3 (e-h) LPR+FBP LPR-Sharp+FBP NLPR+FBPPhase . ⇥ .
36 4 . ⇥ Ref. In. Dec. . ⇥ . ⇥ . ⇥ (b) RMSE computed over whole reconstructed volume Quantifying performance for multi-material spheres
GT LPR+FBP LPR-Sharp+FBP NLPR+FBP . µ m . µ m . µ m . µ m . µ m . µ m . µ m . µ m . µ m . µ m . µ m . µ m (c) Areas of circles (spherical cross-sections) in Fig. 3 (m-p) LPR+FBP LPR-Sharp+FBP NLPR+FBP . ⇥ . ⇥ . ⇥ . ⇥ . ⇥ . ⇥ . ⇥ . ⇥ . ⇥ (d) RMSE within the region of each circle in Fig. 3 (n-p) Figure 5. (a,b) and (c,d) quantifies the accuracy and sharpness of the re-fractive index decrement reconstruction for the single material spheres dataand multi-material spheres data, respectively. (a) and (c) gives the areas forthe three circles in Fig. 3 (e-h) and Fig. 3 (m-p), respectively. The areaswith NLPR closely matches the ground-truth (GT) while LPR and LPR-Sharpoverestimate the area. (b) computes the root mean squared error (RMSE)over the whole volume for single material spheres reconstruction. (d) com-putes the RMSE within each circle for the multi-material reconstructed im-ages in Fig. 3 (n-p). reconstructions with significant amount of blur that also containphase contrast fringes as shown in Fig. 3 (g,o). From Fig. 3 (m,p),we can see that NLPR accurately reconstructs the shape in spiteof the presence of materials with different and varying values for db . In Fig. 3 (n-p), while it was assumed that db = db for the three spheres were 35, 350, and 700.We also observed that the smoothing nature of LPR overesti-mates the area of region occupied by the spheres. To demonstratethis effect, we segment each circle in Fig. 3 (e-h,m-p), computethe area of each circle, and compare the computed areas in Fig.5 (a,c). The region within each red bounding box as shown inFig. 3 (e) is extracted and segmented using Otsu thresholding[14]. After thresholding, each circular area is computed by sum-ming the number of positive pixels and multiplying by the squareof the pixel width. From Fig. 5 (a,c), we can see that the ar-eas with NLPR almost exactly matches the ground-truth whilethe areas are severely overestimated when using LPR. For the sin-gle material data, the root mean squared error (RMSE) in Fig. 5(b) between the reconstructed and ground-truth values over thewhole volume shows that NLPR has higher quantitative accuracythan LPR. We also note that LPR-Sharp results in severely de-graded quantitative accuracy due to the incorrect assumption forthe object-to-detector distance R . In Fig. 5 (d), RMSE estimates (a) Radiograph (b) Sinogram(c) LPR (d) NLPR Figure 6.
Experimental data measurements and reconstructed phases.(a) shows the measurements at the first view angle. (b) shows the sinogramover all the views along the center row of (a). (c) and (d) show the retrievedphase using LPR and NLPR, respectively. computed within the region of each circle in Fig. 3 (n-p) is pre-sented for multi-material data. We can see that LPR has a lowerRMSE than NLPR for the first sphere with db =
35, while NLPRhas lower RMSE for the other two spheres.
Experimental Data
In order to validate NLPR with experimental data, we used asynchrotron to perform X-ray phase contrast imaging of an objectconsisting of SiC fibers. Radiographs were acquired at 256 dif-ferent views equally spaced over an angular range of 180 degrees.The shape of each radiograph is 320 ⇥
324 with a pixel widthof 0 . µ m . The X-ray energy was 20 keV and the object-to-detector distance was fixed at R = mm which ensured signif-icant amount of phase contrast in the measured images as shownin Fig. 6 (a,b). Fig. 6 (a) shows the normalized measurements atthe first view angle and Fig. 6 (b) shows the sinogram along thecenter row of Fig. 6 (a) over all the view angles. Before runningLPR and NLPR, each X-ray image was padded to 1 . ⇥ the orig-inal size in each dimension using edge padding. For the LBFGSoptimization in NLPR, we used a value of 10 for the relativethreshold on the optimization parameters, x , in NLopt [17]. Inequation (13), we used a = g = IS&T International Symposium on Electronic Imaging 2020Computational Imaging a) LPR+FBP (b) NLPR+FBP(c) LPR+FBP (d) NLPR+FBP (e) MTF from (c) and (d)
LPR+FBP NLPR+FBP67.40 µ m µ m µ m µ m (f) Areas of circles within the red boxes in (a) and (b). Figure 7.
Comparison between FBP reconstructed images from LPR andNLPR phase projections. (a) and (b) show reconstructed image slices corre-sponding to the center row of Fig. 6 (a) and the sinogram in Fig. 6 (b). (c)and (d) zooms and compares the circle bounded by a red rectangle on theleft side of (a) and (b), respectively. (e) compares the modulation transferfunction (MTF) computed from the circles in (c) and (d). (f) compares theareas of the circles bounded by the two red rectangles in (a) and (b).
The modulation transfer function (MTF) of the circles in Fig. 7(c,d) are computed using the method in [15] and shown in Fig.7 (e). We can see that the MTF with NLPR lies above the MTFwith LPR, which indicates sharper images with higher resolutionwhen using NLPR. We employ the same approach used to pro-duce Fig. 5 (a,c) and use it to compute the areas of the circularregions within the red bounding boxes in Fig. 7 (a,b). We cansee from the computed areas in Fig. 7 (f) that the areas when us-ing NLPR with FBP is lower than the areas when using LPR withFBP. Furthermore, the areas in Fig. 7 (f) can be compared withthe areas computed for simulated data in Fig. 5 (a). From Fig. 5(a), we have reason to conclude that LPR has overestimated theareas of the circular regions in Fig. 7 (f).While our phase retrieval algorithms assumed only one ma-terial with a constant d / b , we see that there is a second materialthat is used to hold the fibers in place as shown within the redrectangular box in Fig. 8 (a,b). Fig. 8 (c,d) zooms into the redrectangular region in Fig. 8 (a,b) and clearly shows the sharpnessimprovement when using NLPR. Thus, the desirable properties of (a) LPR+FBP (b) NLPR+FBP(c) LPR+FBP (d) NLPR+FBP Figure 8.
Comparison between FBP reconstructed images from LPR andNLPR phase projections. (a) and (b) compares the th cross-section slicefrom the refractive index decrement reconstructions. (c) and (d) zooms intothe area within the red rectangular region of (a) and (b). NLPR also apply to materials that do not have the same specifiedvalue of d / b within the sample. Conclusions
We presented a non-linear phase retrieval (NLPR) methodand validated its performance using both simulated and experi-mental data. We compared NLPR with a popular linear phaseretrieval (LPR) approach presented in Paganin et al. [4]. In ourexperiments, we showed that NLPR resulted in more quantita-tively accurate reconstructions than LPR. We also showed thatLPR significantly overestimated the object’s area due to blurringwhile NLPR accurately estimated the object’s area by reducingthe blur.
Acknowledgments
LLNL-PROC-803597. This work was performed under theauspices of the U.S. Department of Energy by Lawrence Liver-more National Laboratory under contract DEAC52- 07NA27344.Lawrence Livermore National Security, LLC. LDRD Fundingwith tracking number 19-ERD-022 was used on this project. Thisresearch used resources of the Advanced Light Source, a DOEOffice of Science User Facility under contract no. DE-AC02-05CH11231. We thank Jeff Kallman, Jean-Baptiste Forien, andKyle Champley from LLNL for useful discussions.This document was prepared as an account of work spon-sored by an agency of the United States government. Neitherthe United States government nor Lawrence Livermore NationalSecurity, LLC, nor any of their employees makes any warranty,expressed or implied, or assumes any legal liability or respon-sibility for the accuracy, completeness, or usefulness of any in-formation, apparatus, product, or process disclosed, or representsthat its use would not infringe privately owned rights. Referenceherein to any specific commercial product, process, or serviceby trade name, trademark, manufacturer, or otherwise does notnecessarily constitute or imply its endorsement, recommendation,
IS&T International Symposium on Electronic Imaging 2020Computational Imaging r favoring by the United States government or Lawrence Liver-more National Security, LLC. The views and opinions of authorsexpressed herein do not necessarily state or reflect those of theUnited States government or Lawrence Livermore National Se-curity, LLC, and shall not be used for advertising or product en-dorsement purposes. The United States Government retains andthe publisher, by accepting the article for publication, acknowl-edges that the United States Government retains a non-exclusive,paid-up, irrevocable, world-wide license to publish or reproducethe published form of this article or allow others to do so, forUnited States Government purposes.
References [1] K. A. Mohan, X. Xiao and C. A. Bouman, Direct model-based to-mographic reconstruction of the complex refractive index, IEEE Int.Conf. on Im. Proc. (ICIP), pg. 1754-1758 (2016).[2] A. Burvall, U. Lundstr¨om, P. A. C. Takman, D. H. Larsson, and H.M. Hertz, Phase retrieval in X-ray phase-contrast imaging suitable fortomography, Opt. Express, 19, pg. 10359-10376 (2011).[3] A. V. Bronnikov, Reconstruction formulas in phase contrast tomogra-phy, Opt. Comm., 171, 46, pg. 239 – 244 (1999).[4] D. Paganin, S. C. Mayo, T. E. Gureyev, P. R. Miller, and S. W.Wilkins, Simultaneous phase and amplitude extraction from a singledefocused image of a homogeneous object, J. of Microscopy, 206, pg.33-40 (2002).[5] T. E. Gureyev, T. J. Davis, A. Pogany, S. C. Mayo, and S. W. Wilkins,Optical phase retrieval by use of first Born and Rytov-type approxi-mations, Appl. Opt., 43, 12, pg. 2418–2430 (2004).[6] X. Wu, H. Liu, and A. Yan, X-ray phase-attenuation duality and phaseretrieval, Opt. Lett., 30, 4, pg. 379–381 (2005).[7] M. Langer, P. Cloetens, J.-P. Guigay, and F. Peyrin, Quantitative com-parison of direct phase retrieval algorithms in in-line phase tomogra-phy, Med. Phys., 35, pg. 4556-4566 (2008).[8] J. Guigay, M. Langer, R. Boistel, and P. Cloetens, Mixed transferfunction and transport of intensity approach for phase retrieval in theFresnel region, Opt. Lett., 32, pg. 1617-1619 (2007).[9] J. Petruccelli, L. Tian, and G. Barbastathis, The transport of inten-sity equation for optical path length recovery using partially coherentillumination, Opt. Express, 21, pg. 14430-14441 (2013).[10] V. Davidoiu, B. Sixou, M. Langer, and F. Peyrin, Nonlinear ap-proaches for the single-distance phase retrieval problem involvingregularizations with sparsity constraints, Appl. Opt., 52, pg. 3977-3986 (2013)[11] V. Davidoiu, B. Sixou, M. Langer and F. Peyrin, Non-linear iterativephase retrieval based on Frechet derivative and projection operators,2012 9th IEEE Int. Symp. on Bio. Im. (ISBI), pg. 106-109 (2012).[12] A. S. Jurling and J. R. Fienup, Applications of algorithmic differ-entiation to phase retrieval algorithms, J. Opt. Soc. Am. A, 31, pg.1348-1359 (2014).[13] M. A. Herraez, D. R. Burton, M. J. Lalor, and M. A. Gdeisat, Fasttwo-dimensional phase-unwrapping algorithm based on sorting by re-liability following a noncontinuous path, J. Appl. Opt., 41, 35, pg.7437 (2002).[14] N. Otsu, A Threshold Selection Method from Gray-Level His-tograms, IEEE Tran. on Sys., Man, and Cybernetics, 9, 1, pg. 62-66(1979).[15] E07 Committee, ASTM E1695-95, Standard Test Method forMeasurement of Computed Tomography (CT) System Performance,ASTM International. [16] A. C. Kak and Malcolm Slaney, Principles of Computerized Tomo-graphic Imaging, Society of Ind. and Appl. Math., 2001[17] S. G. Johnson, The NLopt nonlinear-optimization package,http://github.com/stevengj/nlopt[18] J. Nocedal, Updating quasi-Newton matrices with limited storage,Math. Comput., 35, pg. 773-782 (1980).[19] D. C. Liu and J. Nocedal, On the limited memory BFGS method forlarge scale optimization, Math. Prog., 45, pg. 503-528 (1989).
Author Biography
K. Aditya Mohan (ORCID: https://orcid.org/0000-0002-0921-6559;Email: [email protected]) received his B.Tech. degree in electronics andcommunication engineering from National Institute of Technology Kar-nataka in 2010. He received his M.S. and Ph.D. degrees in electricaland computer engineering from Purdue University in 2014 and 2017, re-spectively. He is currently a Signal and Image Processing Engineer inthe Computational Engineering Division at Lawrence Livermore NationalLaboratory. His research interests include computational imaging, in-verse problems, and machine learning.Dula Parkinson (Email: [email protected]) received his BS inchemistry from Brigham Young University (2001), his PhD in physicalchemistry from UC Berkeley (2006), and completed a postdoctoral fel-lowship at UC San Francisco (2010). Since then he has worked as astaff scientist at the Advanced Light Source at Lawrence Berkeley Na-tional Laboratory, focusing on instrumentation and image processing forsynchrotron microCT.Jefferson Cuadra (Email: [email protected]) received his BS andPhD in mechanical engineering from the New Jersey Institute of Technol-ogy (2010) and from Drexel University (2015), respectively. He is part ofthe Materials Engineering Division at Lawrence Livermore National Lab-oratory in Livermore, CA. His work has been on metrology applicationsusing X-ray CT, in addition to fatigue and fracture analysis of materialsusing FEM and nondestructive testing and evaluation.