X-ray computed tomography using partially coherent Fresnel diffraction with application to an optical fiber.

A reconstruction algorithm for partially coherent x-ray computed tomography (XCT) including Fresnel diffraction is developed and applied to an optical fiber. The algorithm is applicable to a high-resolution tube-based laboratory-scale x-ray tomography instrument. The computing time is only a few times longer than the projective counterpart. The algorithm is used to reconstruct, with projections and diffraction, a tilt series acquired at the micrometer scale of a graded-index optical fiber using maximum likelihood and a Bayesian method based on the work of Bouman and Sauer. The inclusion of Fresnel diffraction removes some reconstruction artifacts and use of a Bayesian prior probability distribution removes others, resulting in a substantially more accurate reconstruction.


Introduction
The spatial resolution achieved by x-ray computed tomography has improved by over five orders of magnitude from 2 mm in early medical imaging [1]. Groups using synchrotron sources have pioneered the reduction of spatial resolution from the 15 µm resolution achieved in 1993 on rock samples [2] to about 10 nm achieved in recent years using ptychography [3], a lensless technique which depends intimately on the analysis of diffraction patterns. Hard x-ray projective nanotomography at synchrotrons is also practiced [4]. The capabilities of laboratory sources have also increased markedly, with widely available commercial instruments offering a resolution of 3 µm for large industrial parts and a few vendors offering resolutions below 100 nm for smaller samples.
A quarter century back, it was discovered that diffraction effects are also visible in x-ray images from tube sources despite their broad-band nature [5]. Such sources are nearly universally used in laboratory-based x-ray tomography. Although diffraction effects in tomography are a well-researched field [6], geometrical projection remains the dominant paradigm. Typically, when the Fresnel number F = a 2 /(λz eff ) is comparable to or less than 1, diffraction effects are significant. Here, a is the feature size, λ is the wavelength and z eff is the effective distance given by z −1 eff = z −1 1 + z −1 2 where z 1 is the source-sample distance and z 2 is the sample-detector distance [7]. In the case of parallel illumination, z eff = z 2 . The rule for addition of inverse distances is at the heart of the thin lens law, and has been used for many years in x-ray phase-coherent imaging, including a beamline with variable magnification [8]. As the feature size a decreases, the Fresnel number decreases. Additionally, less energetic x rays typically provide for efficient interactions has been used to image magnetic systems on the micrometer scale [29] where polarization dependence is important [30]. Here, the x-ray wavelength, typically 100 pm, is small compared to the minimum feature size, which is near 1 µm. Hence, scalar diffraction theory is sufficient. Moreover, we will consider the case where the width of the detector is small compared to the effective sample-detector distance, so we can use the Fresnel propagator.

Projective tomography
Earlier, some of us presented a code to perform projective tomography using a multi-spectral source [31]. Our method builds on projective tomography so the formulation is presented again to make this paper self-contained. We retain the multi-spectral, multi-material notation, which is incorporated into the code. In the context of tomography, the use of two spectra is often called "dual energy." In the application in this paper, we only use a single spectrum and single basis material, although physically, we have more than one material present in the sample. Materials in the sample such as acrylic and silica are differentiated by a single real number denoting the density of the basis material in a given voxel.
The intensity I jψ is observed for each source spectrum j and viewing condition ψ is given by where f ⃗ ri is the concentration of basis material i at voxel ⃗ r, D(E) is the detection efficiency at photon energy E, and A ⃗ rψ is the system matrix which consists of the lengths of line segments travelling through voxel ⃗ r for viewing condition ψ. The viewing condition ψ is a joint set of viewing angles and projections. We have one projection per detector pixel. Continuing, is a source spectrum and j is the spectrum index. The interaction of material i and the beam is represented by an absorption coefficient α i (E). The process of "reconstruction" is the determination of the f ⃗ ri for each voxel at ⃗ r and each basis material i. In Eq. (1), all indices are given explicitly. The voxel size will be large compared to a wavelength but small compared to the system dimensions. All other variables are treated as known, either because they are measured, are parameters from the experiment, or, in the case of the spectrum given by an assumed model. The absorption coefficients are tabulated by several sources. We use those of the Center for X-ray Optics (CXRO) [32] which also offers the complex index of refraction, which we need for Fresnel diffraction. The projection integral [33] appears in Eq. (1). If there are N × N × N voxels in the reconstruction, then A ⃗ rψ can have no more than 3N non-zero values for a given ψ [34,35]. In principle, the dimensions of A are N 3 by the product of the number of detector pixels times the number of viewing angles. We expect the number of detector pixels to be O(N 2 ) and the number of viewing angles to be O(N). For example, for Nyquist-limited sampling in 2D projection tomography, the number of viewing angles should be π 2 N or greater [36]. Hence A has O(N 6 ) elements of which O(N 4 ) are non-zero. This matrix is too large to store, so it is computed as needed.
We find the best solution by minimizing an objective function which considers both the differences of the predictions of Eq. (1) from the observations as well as prior information about the reconstruction. Specifically, we will use the function suggested by Bouman and Sauer [26,37], which is given by where ⃗ n has elements n J , which are the counts observed at the joint index of observation J = (jψ), ⃗ f is the proposed reconstruction whose components are the f ⃗ ri introduced above, and g( ⃗ f ) is the prior distribution. The first term in the log-likelihood function is derived assuming that each observation n J obeys the Poisson distribution with mean I J : Minimizing L MAP gives the maximum a posteriori (MAP) estimate whereas minimizing the negative of the log likelihood L ML alone is the Maximum Likelihood (ML) method. In this work, we follow the suggestion of Bouman and Sauer [26] and take the prior probability distribution to be where the sum is taken over neighboring pairs. We specialize to isotropic cubic voxels where the neighboring pairs include six faces with c ⃗ r⃗ r ′ = 1, twelve edges with c ⃗ r⃗ r ′ = 1/ √ 2, and eight vertices with c ⃗ r⃗ r ′ = 1/ √ 3. Bouman and Sauer [26] give a convergence proof for the case of 1<p ≤ 2 and recommend p = 1.1 for the material inspection problem, corresponding to distinct levels with abrupt edges. The optical fiber conforms to this with the exception of the graded index in the core. We adopt p = 1.1 in the present study. The prior distribution is restricted to the case of a single basis material. An appropriate prior distribution for the multi-material case is a research topic [38].
The objective function is maximized by initializing the f i⃗ r with random numbers and applying the L-BFGS-B code of Ref. [39] (hereafter BFGS). The BFGS algorithm uses values of the function and its gradient. It is possible to find these together with considerable reuse of intermediate results. Optimization using the BFGS algorithm requires the calculation of the gradient of the objective function: We impose the constraint that each material in each voxel is non-negative, i.e., f ⃗ ri ≥ 0. We apply the BFGS algorithm with a 128 dimensional subspace, a value which was near optimal in a study with the program on a different example [31]. This figure is more than the range of 3 to 20 recommended by Ref. [39]. If the Bayesian term g( ⃗ f ) is included, we proceed in two stages: first, the ML solution is found, then it is used as a starting point for the Bayesian reconstruction.

Tomography with the Fresnel propagator
Following Paganin [40], we treat Fresnel diffraction as follows: we accumulate changes in both phase and amplitude on projections through the sample. We neglect diffraction within the sample itself, but we include it when considering propagation from the sample to the detector. The key equation is Eq. (2.39) of Ref. [40], which, with a small change of notation, is Here, X, Y, and Z are coordinates in the source-detector frame, U is the complex-valued scalar wave, Z i and Z f are two planes surrounding the sample, Z fi = Z f − Z i , and δ and β are related to the complex index of refraction of the material by The wavevector k = 2π λ , is related to the photon energy by E = ℏck where ℏ is the reduced Planck constant and c is the speed of light. In our example below, we use tabulated values of the complex index of refraction [32,41] for silica at 2.2 g/cm 3 . We hope that the distinction between the index of refraction n(E) and the observations n J will be clear from context.

Fresnel propagator
We assume that there is a point source for the x rays: we neglect the finite extent of the source. Straightforward application of the Fresnel propagator (i.e., from the sample plane to the detector plane) is challenging numerically because of the need to treat the outgoing spherical wave, which has rapid oscillations far off-axis, particularly at large propagation distances. Here, we treat the spherical wave as an analytic function with the sample imposing a slowly-varying multiplicative modification in phase and amplitude given by Eq. (7). The sampling requirements are the same as that of projective tomography. We implement the forward model with three steps: (1) propagating a point source to the sample plane and modulate it using the sample transmission function; (2) back propagate the scalar wave to the source plane; and (3) propagate the scalar wave to the detector plane. The propagation steps are done through free space (i.e., without the sample) after the phase of the scalar wave in the sample plane has been determined by the projection integrals of Eq. (7). The method is equivalent to the application of the Fourier magnification theorem [40,42]. However, our implementation does not invoke the theorem explicitly or make use of transformed variables.
The Fresnel propagator [43] that takes the scalar amplitude from a plane Z = Z a to a plane Z = Z b is given by where (X a , Y a ) are Cartesian coordinates in the Z = Z a plane with a similar relation for and the domain of integration is the whole Z = Z a plane. U is energy-dependent, but we suppress the variable E to keep the formulas readable.
The diffraction calculation proceeds with these planes: Z = Z 0 the plane including the point source, Z = Z 1 the mid-point of the sample, and Z = Z 2 the detector plane. Since the outgoing spherical wave is locally a plane wave, the direction for the projection is the radial direction away from the source. The initial and final points in Eq. (7), Z i and Z f , are near Z = Z 1 . In the Fresnel propagation portion of the calculation, we assign the projected phase to the Z = Z 1 plane.
Starting from a point source of unit integrated strength, taken to be at (X 0 = 0, Y 0 = 0, Z 0 ), the scalar wave on the Z = Z 1 plane, without considering the effect of the sample, is The effect of the sample treated in the projection approximation is to introduce a phasor via The projection integral of Eq. (12) is similar to the one found in Eq. (1), although one is an integral over the complex index of refraction and the other is an integral over the real absorption coefficient. Physically, we arrive at this approximation by compressing the sample into a small strip while preserving the projected mass. In our implementation, we neglect the small difference between projections parallel to the Z axis and those directed away from the source. Like U, b is energy dependent.
Equation (10) is not easily computed numerically because of the quadratic phase factor. However, propagating back through free space to the plane of the origin [43], the scalar wave is In going from line 1 to line 2 in Eq. (13), the quadratic phase factor in (X 1 , Y 1 ) cancels. The bandwidth of b(X 1 , Y 1 ) is comparable to functions that arise in projection tomography. We are interested in the wave on the detector. Hence, we propagate forward from the source plane Z = Z 0 to the detector plane Z = Z 2 resulting in where Equation (14) describes free space propagation performed as if the sample is not present. Physically, we have already accounted for the sample through b(X 1 , Y 1 ). Computationally, we find the scalar wave bỹ and then we find the intensity by squaring the amplitude, i.e., The only difference between U(X 2 , Algorithmically, we do the following: 5. these steps are performed for each photon energy in the spectrum.

Gradient of likelihood for the Fresnel propagator
The code relies on the BFGS algorithm [39] to maximize the log likelihood (or the weighted log likelihood in the Bayesian case). The algorithm requires the calculation of the gradient. In the case of Fresnel diffraction, the challenge is to organize the algorithm so that the time to calculate the gradient scales like the time to calculate the projections. We have achieved this previously for projective tomography [31]. In the Fresnel case, the calculation is similar but we will assume that the influence of a given projection on the detector affects O(N) of the N 2 detector pixels. The analytic example in subsection 2. The change in log-likelihood L with respect to a change in a voxel value is given by Eqs. (3) and (6). The observed intensities I J are jointly indexed by J, which was previously decomposed into a spectrum index j and a joint viewing-angle and detector pixel variable ψ. Here, we further decompose ψ = (⃗ s 2 , v) where ⃗ s 2 is a detector pixel and v is a viewing angle. For the Fresnel case, we may write This equation is similar to Eq. (1) in the projective case. In turn Here, U 2 and U 1 are the scalar waves on the detector and sample planes, respectively, and G is the Fresnel propagator, a Green's function. The scalar wave function is given by where P ⃗ s 1 vi is a projection and U (0) 1 (⃗ s 1 , E, v) is the scalar wave at the sample plane if the sample were absent andα The projection is given by Eq. (2), with the complexα i (E) being analogous to the real α i (E). BFGS requires the derivatives of Eq. (3), which in turn requires the derivatives of Eq. (18).
Equation (22) requires the following expression Equation (23) requires the derivative of Eq. (20) namely using Eq. (2), the definition of the projection, in the final line. These equations need to be reconsidered for efficiency. Equations (3) and (6) can be recast [44] as Equation (25)  and A ⃗ rs 1 v is very easy to compute. In the case of partial coherence, we may approximate G(⃗ s 2 , ⃗ s 1 , E) by a short-range function, anticipating a cancellation at long range after the average over energies. An analytic example of such cancellation is given in the next section. We refer to the region of G(⃗ s 2 , ⃗ s 1 , E) which we represent numerically as the "non-negligible region." The organization of the Fresnel gradient is similar to the projective case. Instead of looping over the detector pixels ⃗ s 2 , we loop over a set of virtual pixels on the sample denoted by ⃗ s 1 . In practice the points ⃗ s 1 and ⃗ s 2 will be in 1:1 correspondence, i.e., we use one projection per detector pixel. We use the fast Fourier transforms of FFTW [45] to implement the Fresnel propagator. The gradient calculation, simplified for one energy E, one material i, and one view v, proceeds as follows: Loop The projective case also scales as O(N 4 ) albeit with a smaller coefficient. We anticipate the projective case will converge more readily because there are fewer approximations in the calculation of its gradient (specifically, the approximation of short-range influence of the partially coherent Fresnel propagator).
As a practical matter, we start the diffractive solution from the projective solution. In testing using random dot patterns, we found that the diffractive program could has a finite basin of attraction when starting from random numbers.

Analytic example of partial coherence
For tomography, we need to find the complex index of refraction at each voxel in some defined portion of space. In order to use BFGS, we need the gradient of the diffraction pattern with respect to changes in the index of refraction at a given voxel. Here, we give an analytic example of a point in the plane interacting with an unscattered wave that motivates the short-range approximation made numerically above. The gradient of the diffraction pattern is essentially the interference pattern of a point source in the sample and the rest of the wave.
Suppose we have a single voxel that differs from the vacuum value by a differential amount. We want to calculate the change in the diffraction pattern to first order. Also, suppose that there is a multi-wavelength point source located at the origin; the source plane is z 0 = 0. The voxel is located at (x 1 , y 1 , z 1 ). The detector plane is located in the plane z = z 2 . The unscattered wave is denoted by U and the scattered wave by U (1) and is given by which may be found by considering Eq. (9) with a point source.
The first-order interference pattern on the plane z = z 2 for a particular wavevector k is given by which is a first order expansion of Eq. (17). Next, we assume the source intensity is normally distributed in k, hence also normally distributed in the photon energy. We want to find where ξ is half of the term in brackets in Eq. (27).
To understand the Gaussian envelope, we can rearrange The result shows the diffraction pattern is centered on the geometric projection of (x 1 , y 1 ) onto the z 2 plane. The first-order partially coherent diffraction pattern of a displaced point is centered on the projection of the scattering point onto the detector plane under geometric magnification. It has the functional form of the coherent diffraction pattern multiplied by a Gaussian. In contrast to the coherent diffraction pattern, the partially coherent diffraction pattern is short ranged. An example is shown in Fig. 1.

Experiment and results
We study a graded index optical fiber, Thorlabs GIF625 (Newtown, NJ, USA). The optical fiber has a core composed of germanium-doped silica with a diameter of 62.5 µm ± 2.5 µm, a silica cladding layer with a diameter of 125 µm ± 1 µm, and an acrylic coating with a diameter of 245 µm ± 10 µm, where the uncertainties are the tolerances stated by the manufacturer. The core region has a graded index of refraction with the density of germanium increasing toward the center.

Microscopy
An XCT study of the fiber by some of us appeared earlier [24]. There appeared to be two layers in the acrylic coating, which was a surprise as it did not appear in the manufacturer's specifications. Accordingly, we performed visible light microscopy and scanning electron microscopy (SEM) to check for the presence of such a double layer. For these results, the optical fiber was cut using a ceramic scribe, and a small sample piece was attached to an aluminum sample mounting stub using conductive carbon tape. A visible light microscope equipped with a 50x objective and a Zeiss Gemini 300 Variable Pressure SEM were used to image the sample. The visible light and SEM images are shown in Fig. 2. Each image shows that the coating comprises two layers separated by a very thin boundary. Higher magnification SEM images did not reveal evidence of a meaningful gap between the layers.

X-ray image acquisition and alignment
The sample was imaged using a Zeiss Versa XRM-500 (Concord, CA, USA) using an x-ray tube voltage of 60 kV. The sample was attached to a 3 mm diameter aluminum nail, which was inserted into a pin-vise holder. The source to sample distance was z 1 = 16 mm and the sample to detector distance was z 21 = 25 mm, leading to z eff = 9.76 mm and a geometric magnification M = 2.563 for the x rays. In addition, an optical magnification of 20x was used. We obtained a tilt series about a single rotation axis of 801 images of size 495 × 495 spaced over 180 • with an angular increment of 0.225 • . The tilt series images were trimmed to 266 × 465 pixels during alignment, which included both translations and rotations of the images. Reconstructions were performed on 100 central slices of the aligned images.
The tilt series images of the optical fiber were aligned in a two-step process. First, edge detection was used to extract the left and right external edges of the core followed by fits to the line using the RANSAC algorithm [46]. From these edges, shifts of integer pixels were made in horizontal dimension to give a preliminary alignment. The rotational alignment was done by fitting the extracted slopes γ ′ (α) of the fiber to that of a projection of a rigid cylinder and correcting for global tilt using where α is the angle of the tomographic projection, α 0 gives the initial rotational position of the sample, γ is the tilt of the cylinder with respect to the rotation axis, and β is the misalignment angle between the projection axis, i.e., the detector y axis, and the rotation axis, as shown in Fig. 3. Both α and α 0 denote angles in the same space as the angles of the tilt series. The angles α and α 0 are unrelated to α i above; similarly,r andr ′ (α) have no relation to the symbols ⃗ r and ⃗ r ′ denoting the voxels.  3. (a) The motion of the optical fiber in the source-detector framex,ŷ,ẑ undergoing rotation, α, around the system rotation axisŷ can be described by a vectorr with a maximum tilt γ with respect to the rotational axis, and an initial rotation position with respect to the projection angle of maximum tilt, α 0 . (b) The detected slope in the projection image is given by the angle between the projected fiber vectorr ′ (α) and the detector y axis,ŷ ′ . If the detector framex ′ ,ŷ ′ is not perfectly aligned with the system frame, the observed projected slope γ ′ (α) has an additional offset β and is given by Eq. (30).
The tilt series images of the optical fiber were then further aligned using the program Arec3d [47]. It is a model-based alignment method, where a least-squares formulation of the tomographic inversion is solved while simultaneously refining the alignment by aligning the experimental projection images to computational projections of the reconstruction. Due to sufficient accuracy during pre-alignment, the rotational alignment in Arec3d was disabled. Additionally, owing to the relatively small amount of features in the image, a high-pass filter was used in the alignment step of the algorithm.
Arec3d reports angle β is 0.14(1) • , indicating that the tilt axis and the detector are nearly aligned, and the angle γ to be 3.95(1) • as shown in Fig. 4. The digits in parentheses represent one standard error in the fit. We also find the angle between the fiber and the reconstruction axis from an analysis of the reconstructed stack to be 3.87 (14) • , a value which should be and is consistent with γ.

Reconstructions
Reconstructions were performed on the aligned images with a code developed in house. Results from a version including projections and scatter corrections have appeared [31]; however, scatter corrections were not used in the present project. Instead, the new development is the Fresnel tomography applied above. Additionally, a Bayesian prior probability distribution that was developed for the material inspection problem [26] was used optionally.
The main parameters used in the calculation are given next. The x-ray spectrum is taken to have 65 photon energies with a Gaussian intensity centered at 10 keV, with a standard deviation of 2 keV running from 2 keV to 18 keV. We take D(E) = 1, i.e., detector efficiency is assumed to be ideal. We retain 266 × 100 detector pixels after alignment for each of 801 views. The detector pixel size is 2.691 µm. Half of the length of the diagonal of the detector after alignment is 382 µm, whereas the source to sample distance is 41 mm. The ratio is 0.009<<1, so a validity condition for the Fresnel approximation is obeyed.
We reconstruct into 186 × 186 × 70 cubic voxels, each of which is 1.5 µm on a side. The background intensity was taken to be 1318.30 counts, which is the square of the signal-to-noise of the background. This scaling is required due to our assumption of Poisson statistics in Eq. (4). To make contact with the formalism, we reconstruct with a single basis material, so i = 1, we have a single spectrum, so j = 1, we have 65 photon energies, so the subscript E = 1, . . . , 65, and there are 801 × 266 × 100 = 21 306 660 values of ψ.
At an x-ray photon energy of 10 keV, hence a wavelength λ of 124 pm, using the formula in the introduction, the Fresnel number F = 0.91, taking a to be the detector pixel size referred to the sample. Noting that F ∝ E, where E is the photon energy, we may express the range of Fresnel numbers for our Gaussian model spectrum using uncertainty notation as F = 0.91 (18). Since F is near 1, we expect diffraction effects to be a moderate correction to geometric optics.
Results are shown in Fig. 5. The optical fiber has three components, the core (the center to feature 4), the cladding, from feature 4 to feature 1, and the coating, from feature 1 to feature 2. The core region has a graded index of refraction. The cladding consists of silica without germanium. Surrounding the optically active area is an acrylic coating, provided for mechanical strength. The reconstruction of Fig. 5(a) was performed with projective tomography with the maximum likelihood method. It is very similar to one obtained with vendor-supplied software, which was presented earlier [24]. Figure 5(c) is a reconstruction using maximum likelihood, but including the effect of Fresnel diffraction. Compared to Fig. 5(a), feature 1, at the cladding-coating boundary and feature 2 at the coating-air boundary are edge diffraction artifacts (in Fig. 5(a)) which are greatly reduced when diffraction is included in the physical model of reconstruction as shown in Fig. 5(c).
There is a stipple artifact in the core region, between the two labels "4". The stippling is also present in the maximum likelihood diffraction reconstruction of Fig. 5(c). When a Bayesian prior is used, the stipple is eliminated in the projective maximum likelihood reconstruction of Fig. 5(b) as well as the Bayesian/diffraction reconstruction, shown in Fig. 5(d). The stipple is better seen in Fig. 6, a blow-up of the central region.
Cylindrical averages are shown in Fig. 7. In Fig. 7(a), the cylinders were defined as 1 pixel wide rings about the center of the cladding on the 8 central slices. In Fig. 7(a), the reconstructions are seen to give about the same continuous decrease in the core region, followed by a more gentle decrease in the cladding. The core region meets our expectations from the manufacturer's specification, but the more gentle decrease in the cladding region was not expected from the manufacturer's specifications; instead, a flat region was expected. The biggest differences between the approximations comes at the boundary of the cladding and the inner coating, feature 1. Here, the projective reconstructions show a sharp decrease, limited in many cases by the value 0. The diffraction approximation shows a fall-off in density at the boundary within 3 µm or 2 pixels. In contrast, the inclusion of diffraction yields a more physical boundary. There are moderate differences due to the inclusion of the Bayesian prior. In Fig. 7(b), averages were taken to highlight the boundary between the inner coating and outer coating, feature 3. The boundary is approximately circular with a radius of 94.5 µm offset from the core-cladding center by 4.5 µm. Thus, the zero position in Fig. 7(b) occurs just past the end of Fig. 7(a), and the inner coating regions overlap. The diffraction pattern is sufficiently fine that it was necessary to track the individual voxels composing it. This was done by locating zero values on the projective maximum likelihood reconstruction on each layer. The zero values are the center of the diffraction feature. These values did not form complete rings, so the values were augmented by hand interpolation using ImageJ [48]. Dilations of the images on each slice were performed successively to locate voxels a given distance from the zero values. The voxels identified with a given distance were used for all 4 reconstructions, which are created from the same aligned experimental data.
The main result is that the inclusion of Fresnel diffraction leads to a much smaller central dip, with 2 to 3 times less amplitude and a full-width at half-maximum of 3 µm, vs. 6 µm for the projective Bayesian reconstruction and 5 µm for the projective maximum likelihood reconstruction. The SEM image in Fig. 2(b) shows that the layer between the inner and outer coating is much smaller than 1 µm. Hence, the smaller value is a better description of the sample.

Computing times
Computer times per iteration are given in Table 1. These show that the Fresnel diffraction leads to a penalty factor of 3.2 per iteration compared to the projective method, whereas the Bayesian term leads to a penalty factor of only 1.07. The number of iterations is also given in Table 1. Only the projective, maximum likelihood calculation started from random numbers. Its result was used to start the Fresnel, maximum likelihood calculation and the projective, Bayesian calculation. The latter was used to start the Fresnel, Bayesian calculation. The stopping point is given by the "medium convergence" parameter suggested by the BFGS code authors [39]. The results are not strictly comparable because of the structure of the starting points. As anticipated just above Sec. 2.3.1, we are not surprised that the Fresnel calculations required more iterations. Overall, we estimate a one order of magnitude penalty for the Fresnel calculation compared to its projective counterpart.

Conclusions
We have presented an algorithm to treat partial coherence with a computational time that scales like its projective counterpart. A code implementing the algorithm is used to reconstruct a graded-index optical fiber with projections and diffraction, and using maximum likelihood and a Bayesian method based on Bouman and Sauer [26,37]. The two key mathematical circumstances which make this possible are the use of projections to find the effect of the sample on a scalar wave and the short-range influence of a partially coherent sum of scalar waves.
The use of the Bayesian prior removes a stippling effect in the core region. The inclusion of Fresnel diffraction allows a better description of the material close to internal material boundaries than projective tomography. On the other hand, the projective reconstruction was more sensitive to an unanticipated boundary interior to the coating of two nominally identical materials, which could be a useful diagnostic feature if properly interpreted. We observed the boundary with both visible light microscopy and scanning electron microscopy, confirming a conjecture made earlier [24]. X-ray tubes are the most common x-ray source for tomography and have a broadband spectra. This leads to partial coherence. For x-ray nanotomography, it is easy for the partial coherence to lead to diffraction features in what otherwise appears to be a normal image [24]. The algorithm which we have presented here can account for such features with a computational time which scales like projective tomography. We hope that it finds widespread application.