1. 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 $\mu$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 $\mu$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}/(\lambda z_\textrm {eff})$ is comparable to or less than 1, diffraction effects are significant. Here, $a$ is the feature size, $\lambda$ is the wavelength and $z_\textrm {eff}$ is the effective distance given by $z_\textrm {eff}^{-1}=z_1^{-1}+z_2^{-1}$ where $z_1$ is the source-sample distance and $z_2$ is the sample-detector distance [7]. In the case of parallel illumination, $z_\textrm {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 [9] for small samples, so, in practice, $\lambda$ will tend to increase as $a$ decreases, also leading to a smaller Fresnel number.

Because nanoscale resolution is increasingly common, it is timely to reconsider diffraction effects in tomography. The tomography of integrated circuit interconnects, for example, has grown from a synchrotron demonstration [10] to vastly larger demonstrations using synchrotron sources and ptychography, developed in the two decades since diffraction patterns of non-periodic objects were first inverted from experimental data [11]. (By “inverted”, we mean that the real-space images of the samples generating the diffraction were found.) The application of x-ray nanotomography and other imaging techniques to the problem of integrated circuit interconnects has been reviewed recently [12]. In addition to x-ray tomography, Fresnel diffraction has been identified to be a problem in electron tomography of thick samples [13].

Alternative methods of phase retrieval for propagation-based imaging are also used [14]. Recently, the conditions for propagation-based imaging in a laboratory x-ray nanotomography instrument were characterized [7]. Similarly, there are recent studies which incorporate diffraction effects into a ray-tracing-based Monte Carlo simulation of x-ray images [15,16].

Aside from ptychography, which operates in the far field, one common approach to diffraction-based tomography is based on the Transport of Intensity Equation (TIE) [17–19] and is known as Paganin’s algorithm. This has been implemented in the open source code Ankaphase [20]. Paganin’s algorithm is based on applying a single step of Euler’s method for solving differential equations to the TIE. Its validity requires $F>\!>1$ and a monochromatic source, which are far different from the conditions considered in this paper. The condition $F>\!>1$ is equivalent to requiring the propagation distance to be much less than the Rayleigh length in the case of Gaussian beams. Another technique used at synchrotrons is to acquire phase maps for each projection by imaging at several defocus distances, followed by tomographic reconstruction of the reconstructed phase. An initial report of a polystyrene sample reconstructed with 0.95 $\mu$m voxels found phase shifts in excess of 2$\pi$ [21].

Accordingly, there is a need to develop tomographic methods that are based on more robust assumptions for application to broadband tube sources in instruments with spatial resolution near or below 1 $\mu$m. To encourage practical use, these methods should not be much slower than existing projection-based methods, they should allow for a source spectrum in a principled way, and allow operation for intermediate Fresnel numbers, say $F=0.1$ to $F=1$, i.e., between the validity conditions for Paganin’s algorithm and ptychography. The finite bandwidth of the source spectrum necessarily results in partial coherence of the beam. Partially coherent propagation has been considered in accelerator design over the last decade [22,23].

Some of us observed diffraction effects in a tomographic reconstruction of a graded-index optical fiber [24]. An earlier study observed similar features [25]. In both studies, the diffracting regions were avoided in the analysis of the images.

2. Theory and algorithm

In this work, an algorithm for partially coherent x-ray computed tomography, including Fresnel diffraction, is developed and applied to tilt series images of an optical fiber. 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]. Our approach is similar to a likelihood-based method, which appeared recently [27]; however, our treatment of the Fresnel propagator goes beyond the TIE and Paganin’s algorithm [17]. The inclusion of Fresnel diffraction removes some reconstruction artifacts. The use of the Bayesian prior probability distribution removes others.

Diffraction can be described by scalar wave theory as well as by classical electromagnetism. Polarization effects in non-magnetic materials are important when the wavelength is comparable to the feature size [28]. Although we consider non-magnetic systems, vector x-ray nanotomography 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 $\mu$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.

2.1 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\psi }$ is observed for each source spectrum $j$ and viewing condition $\psi$ is given by

The projection integral [33]

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

Minimizing $L_\textrm {MAP}$ gives the maximum a posteriori (MAP) estimate whereas minimizing the negative of the log likelihood $L_\textrm {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

Bouman and Sauer [26] give a convergence proof for the case of $1<p\le 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\vec 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_{\vec r i}\ge 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(\vec 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.

2.2 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_{f\!i}=Z_f-Z_i$, and $\delta$ and $\beta$ are related to the complex index of refraction of the material by

The wavevector $k=\frac {2\pi }{\lambda }$, is related to the photon energy by $E = \hbar c k$ where $\hbar$ 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.

2.2.1 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

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

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 by

The only difference between $U(X_2,Y_2,Z_2)$ and $\tilde U(X_2,Y_2,Z_2)$ is a phasor which does not appear in the intensity $I(X_2,Y_2,Z_2)$.

Algorithmically, we do the following:

The pixels sizes $\Delta _i$ on the source, sample, and detector planes, for $i=0,1,2$ respectively, obey $N \Delta _0 \Delta _1 = \lambda Z_{10}$, $N \Delta _0 \Delta _2 = \lambda Z_{20}$, $\Delta _2 = M \Delta _1$, and $M = {Z_{20}}/{Z_{10}}$. The geometric magnification $M$ of the pixels is consistent with the Fourier scaling theorem [42].

2.3 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.3.1 makes the assumption plausible. Hence, there are $O(N^3)$ operations per view, just as for the projective case. In the projective case, the $O(N^3)$ operations arise from making $O(N^2)$ projections each summing $O(N)$ voxels.

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 $\psi$. Here, we further decompose $\psi =(\vec s_2, v)$ where $\vec 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

The projection is given by Eq. (2), with the complex ${\tilde \alpha _i(E)}$ being analogous to the real $\alpha _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

Equation (25) implies a large computational savings since there are $O(N)$ more instances of $\frac {\partial L_{ML}} {\partial f_{\vec r i}}$ than $\frac {\partial P_{\vec s_1 v i} } {\partial f_{\vec r i}}$ and $A_{\vec r s_1 v}$ is very easy to compute.

In the case of partial coherence, we may approximate $G(\vec s_2, \vec 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(\vec s_2, \vec 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 $\vec s_2$, we loop over a set of virtual pixels on the sample denoted by $\vec s_1$. In practice the points $\vec s_1$ and $\vec 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 on angles {

 Find projections through sample

 Use Fresnel propagator to find $U_2$

 Loop on projections {

 Find $\partial U_1 / \partial P_{\vec s1 i}$

 Find $\partial U_2 / \partial P_{\vec s1 i}$ over the non-negligible region of the Green’s function

 Find $\partial I_{j \vec s_2 v} / \partial P_{\vec s1 i}$

 Find $\partial L / \partial P_{\vec s1 i}$

 Multiply that by the system matrix and accumulate terms in the gradient along one projection

 }

 }

The algorithm is very similar to the projective case. In terms of scaling, the time increases by a constant factor if $G(\vec s_2, \vec s_1, E, v)$ has order $N$ non-zero entries where there are $N^2$ pixels. This is because we need $O(N)$ steps to find each projection, so we can afford to do $O(N)$ operations thereafter. Overall scaling is as $O(N^4)$, taking into account $O(N^2)$ projections per view, $O(N)$ operations per projection, and $O(N)$ view points. 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.

2.3.1 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

The first-order interference pattern on the plane $z=z_2$ for a particular wavevector $k$ is given by

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.

 

Fig. 1. Fresnel diffraction (left) with full coherence and (right) with partial coherence. The fully coherent case is the interference of a point source and a plane wave at 10 keV. The partially coherent case uses the spectral distribution of the source centered on 10 keV with a standard deviation of 2 keV. Both use the experimental geometry.

Download Full Size | PDF

3. 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 $\mu$m $\pm$ 2.5 $\mu$m, a silica cladding layer with a diameter of 125 $\mu$m$~\pm$ 1 $\mu$m, and an acrylic coating with a diameter of 245 $\mu$m$~\pm ~$10 $\mu$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.

3.1 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.

 

Fig. 2. (A) Visible light image and (B) SEM image of a cross section of the optical fiber. The regions shown are: (I) Cladding, (II) Inner Coating, and (III) Outer Coating. The scale bars are each 20 $\mu$m.

Download Full Size | PDF

3.2 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_\textrm {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 $\times$ 495 spaced over $180^\circ$ with an angular increment of $0.225^\circ$. The tilt series images were trimmed to 266 $\times$ 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 $\gamma '(\alpha )$ of the fiber to that of a projection of a rigid cylinder and correcting for global tilt using

 

Fig. 3. (a) The motion of the optical fiber in the source-detector frame $\hat x$,$\hat y$,$\hat z$ undergoing rotation, $\alpha$, around the system rotation axis $\hat y$ can be described by a vector $\hat r$ with a maximum tilt $\gamma$ with respect to the rotational axis, and an initial rotation position with respect to the projection angle of maximum tilt, $\alpha _0$. (b) The detected slope in the projection image is given by the angle between the projected fiber vector $\hat r'(\alpha )$ and the detector $y$ axis, $\hat y'$. If the detector frame $\hat x'$,$\hat y'$ is not perfectly aligned with the system frame, the observed projected slope $\gamma '(\alpha )$ has an additional offset $\beta$ and is given by Eq. (30).

Download Full Size | PDF

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 $\beta$ is 0.14(1)$^\circ$, indicating that the tilt axis and the detector are nearly aligned, and the angle $\gamma$ to be 3.95(1)$^\circ$ 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)$^\circ$, a value which should be and is consistent with $\gamma$.

 

Fig. 4. The image of the optical fiber with the largest angle between the fiber direction and the rotation axis after alignment by Arec3d. The image has an observation angle of 121.275$^\circ$. The central 100 slices used in the reconstruction are shown.

Download Full Size | PDF

3.3 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 $\times$ 100 detector pixels after alignment for each of 801 views. The detector pixel size is 2.691 $\mu$m. Half of the length of the diagonal of the detector after alignment is 382 $\mu$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 $\times$ 186 $\times$ 70 cubic voxels, each of which is 1.5 $\mu$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,\ldots ,65$, and there are $801\times 266\times 100=21\,306\,660$ values of $\psi$.

At an x-ray photon energy of 10 keV, hence a wavelength $\lambda$ 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\propto 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.

 

Fig. 5. Each image of the graded-index optical fiber is the average of 14 reconstructed slices from the center of the image. The reconstructed values were scaled from 0 to 1 in each panel and then a cube root was applied to emphasize the smaller values. The gray scale (right) applies to each subfigure: (a) maximum likelihood projection (b) Bayesian projection (c) maximum likelihood diffraction, and (d) Bayesian diffraction.

Download Full Size | PDF

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.

 

Fig. 6. Magnified images show how the Bayesian term eliminates the stipple in the center of the image. (left) Center of Fig. 5(c), maximum likelihood diffraction. (right) Center of Fig. 5(d), Bayesian diffraction.

Download Full Size | PDF

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 $\mu$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.

 

Fig. 7. Averages of the reconstructed values over 8 central slices and rings described in the text. The units of the reconstruction are arbitrary, but the same for both subfigures. The green lines refer to projective tomography and the blue lines to diffractive tomography. The dashed lines use the maximum likelihood approximation, whereas the solid lines include a Bayesian prior. (a) Averages are taken centered about the cladding-inner coating boundary. The dashed vertical lines mark the manufacturer’s specification for the core-cladding boundary and the cladding-coating boundary. (b) Averages described in the text are taken with the inner-outer coating boundary as defined by the projective, maximum likelihood reconstruction. Positions are relative to that boundary.

Download Full Size | PDF

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 $\mu$m offset from the core-cladding center by 4.5 $\mu$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 $\mu$m, vs. 6 $\mu$m for the projective Bayesian reconstruction and 5 $\mu$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 $\mu$m. Hence, the smaller value is a better description of the sample.

3.4 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.

Table 1. Computer time per iteration for the four methods for the problem described in the text running on an 8 processor computer with a clock speed of 3.7 GHz. The number of iterations is given at right. The program was written in Fortran 90 with MPI.

View Table

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.

4. 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.