TOMOGRAPHIC RECONSTRUCTION APPARATUS AND REMOVING DIFFRACTION EFFECTS IN A TOMOGRAPHIC IMAGE

Abstract

Removing diffraction effects in a tomographic image includes: obtaining an empirical image of a sample; producing an initial wave at a radiation source; forward propagating the initial wave from the radiation source toward a detector; receiving by a sample the initial wave; forward propagating the initial wave through the sample and accumulating phase and amplitude information to produce a phase accumulated wave; back propagating the phase accumulated wave; and forward propagating the phase accumulated wave 208 while treating Fresnel diffraction, such that the empirical image is reconstructed by projections and diffraction via maximum likelihood, a Bayesian prior probability distribution, and a Fresnel propagator.

Description

BRIEF DESCRIPTION

Disclosed is a process for removing diffraction effects in a tomographic image, the process comprising: obtaining an empirical image of a sample, the empirical image comprising diffraction effects; producing an initial wave at a radiation source; forward propagating, in a forward propagation direction from the radiation source toward a detector; receiving, at a radiation source proximate surface of the sample that is interposed between the radiation source and the detector, the initial wave from the radiation source; forward propagating, in the forward propagation direction, the initial wave through the sample and accumulating phase and amplitude information of the initial wave as the initial wave propagates through the sample toward the detector to produce a phase accumulated wave at a detector proximate surface of the sample according to

U(X,Y,Z_f)≈exp{ikZ_f−k∫_Zi^ZfdZ[iδ_E(X,Y,Z)+β_E(X,Y,Z)]}U(X,Y,Z_i);

back propagating, in a reverse propagation direction in absence of the sample, the phase accumulated wave from the detector proximate surface to the radiation source; and forward propagating, in the forward propagation direction from the radiation source toward the detector in absence of the sample, the phase accumulated wave while treating Fresnel diffraction according to

$U\left(X_b,Y_b,Z_b\right)=\frac{e^{{\mathrm{ikZ}}_{\mathrm{ba}}}}{i\lambda Z_{\mathrm{ba}}}\int {\mathrm{dX}}_a{\mathrm{dY}}_a\exp \left\{i\frac{k}{2Z_{\mathrm{ba}}}\left[{\left(X_b-X_a\right)}^2+{\left(Y_b-Y_a\right)}^2\right]\right\}\times U\left(X_a,Y_a,Z_a\right),,$

such that the empirical image of the sample is reconstructed by projections and diffraction via maximum likelihood, a Bayesian prior probability distribution, and a Fresnel propagator.

Disclosed is a process for removing diffraction effects in a tomographic image, the process comprising: determining b(X₁, Y₁) from a plurality of projections and an estimate of the material parameters δ and β; applying the Fourier transform to b(X₁, Y₁) to find B(X₀, Y₀); determining Ũ(X₂, Y₂, Z₂) by multiplying B(X₀, Y₀) by a quadratic phasor and prefactors and applying an inverse Fourier transform according to

$\tilde{U}\left(X_2,Y_2,Z_2\right)=\frac{1}{{\lambda }^3{Z}_{10}^2{Z}_{20}}\int \int {\mathrm{dX}}_0{\mathrm{dY}}_0\exp \left\{i\frac{k}{2}\left(\frac{1}{Z_{20}}-\frac{1}{Z_{10}}\right)\left(X_0^2+Y_0^2\right)\right\}\times B\left(X_0,Y_0\right)\exp \left\{-i\frac{k}{Z_{20}}\left(X_2{X}_0+Y_2{Y}_0\right)\right\}$

BRIEF DESCRIPTION OF THE DRAWINGS

The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawings will be provided by the Office upon request and payment of the necessary fee.

The following description cannot be considered limiting in any way. Various objectives, features, and advantages of the disclosed subject matter can be more fully appreciated with reference to the following detailed description of the disclosed subject matter when considered in connection with the following drawings, in which like reference numerals identify like elements.

FIG. 1 shows a tomographic reconstruction apparatus, according to some embodiments.

FIG. 2 shows a tomographic reconstruction apparatus, according to some embodiments.

FIG. 3 shows 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, according to some embodiments, according to the Example.

FIG. 4 shows (A) a visible light image and (B) a scanning electron microscope (SEM) image of a cross section of the optical fiber, according to the Example. The regions shown are: (I) Cladding, (II) Inner Coating, and (III) Outer Coating. The scale bars are each 20 μm.

FIG. 5 shows (a) the motion of the optical fiber in the source-detector frame x{circumflex over ( )}, y{circumflex over ( )}, {circumflex over ( )}z undergoing rotation, a, around the system rotation axis y{circumflex over ( )} can be described by a vector r″ with a maximum tilt y with respect to the rotational axis, and an initial rotation position with respect to the projection angle of maximum tilt, α0, according to the Example. (b) The detected slope in the projection image is given by the angle between the projected fiber vector r{circumflex over ( )}′(α) and the detector y axis, y{circumflex over ( )}′. If the detector frame x{circumflex over ( )}′, y{circumflex over ( )}′ is not perfectly aligned with the system frame, the observed projected slope y′(α) has an additional offset β and is given by Eq. (30).

FIG. 6 shows, according to the Example, 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°. The central 100 slices used in the reconstruction are shown.

FIG. 7 shows, according to the Example, 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.

FIG. 8 shows, according to the Example, magnified images show how the Bayesian term eliminates the stipple in the center of the image. (left) Center of FIG. 7(c), maximum likelihood diffraction. (right) Center of FIG. 7(d), Bayesian diffraction.

FIG. 9 shows, according to the Example, 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.

FIG. 10 shows, according to the Example, 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 the message passing interface (MPI).

DETAILED DESCRIPTION

A detailed description of one or more embodiments is presented herein by way of exemplification and not limitation.

The spatial resolution achieved by x-ray computed tomography has improved by over five orders of magnitude from 2 mm in early medical imaging. Groups using synchrotron sources have pioneered the reduction of spatial resolution from the 15 μm resolution achieved in 1993 on rock samples to about 10 nm achieved in recent years using ptychography, a lensless technique which depends intimately on the analysis of diffraction patterns. Hard x-ray projective nanotomography at synchrotrons is also practiced. 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. Such sources are nearly universally used in laboratory-based x-ray tomography. Although diffraction effects in tomography are a well-researched field, geometrical projection remains the dominant paradigm. Typically, when the Fresnel number F=a²/(λ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_eff⁻¹=z₁⁻¹+z₂⁻¹ where z₁ is the source-sample distance and z₂ is the sample-detector distance. In the case of parallel illumination, z_eff=z₂. 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. As the feature size a decreases, the Fresnel number decreases. Additionally, less energetic x rays typically provide for efficient interactions for small samples, so, in practice, λ 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 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. Inverted refers to 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. In addition to x-ray tomography, Fresnel diffraction has been identified to be a problem in electron tomography of thick samples.

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

Aside from ptychography, which operates in the far field, one common approach to diffractionbased tomography is based on the Transport of Intensity Equation (TIE) and is known as Paganin's algorithm. This has been implemented in the open source code Ankaphase. 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 μm voxels found phase shifts in excess of 2π.

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

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

It has been discovered that an algorithm for partially coherent x-ray computed tomography, including Fresnel diffraction reconstructs images of a sample with projections and diffraction, using maximum likelihood and a Bayesian process. The Fresnel propagator used herein goes beyond TIE and Paganin's algorithm. The inclusion of Fresnel diffraction removes reconstruction artifacts, and the use of the Bayesian prior probability distribution removes others.

Tomographic reconstruction apparatus 200 acquires empirical image 211 and can remove diffraction effects in the tomographic image 211. In an embodiment, with reference to FIG. 1 and FIG. 2, removing diffraction effects in a tomographic image includes: obtaining an empirical image 211 of a sample 203, the empirical image 211 comprising diffraction effects; producing an initial wave 202 at a radiation source 201; forward propagating, in a forward propagation direction 207 from the radiation source 201 toward a detector 205; receiving, at a radiation source proximate surface 209 of the sample 203 that is interposed between the radiation source 201 and the detector 205, the initial wave 202 from the radiation source 201; forward propagating, in the forward propagation direction 207, the initial wave 202 through the sample 203 and accumulating phase 204 and amplitude information of the initial wave 202 as the initial wave 202 propagates through the sample 203 toward the detector 205 to produce a phase accumulated wave 208 at a detector proximate surface 206 of the sample 203 according to

U(X,Y,Z_f)≈exp{ikZ_f−k∫_Zi^ZfdZ[iδ_E(X,Y,Z)+β_E(X,Y,Z)]}U(X,Y,Z_i);

back propagating, in a reverse propagation direction 210 in absence of the sample 203, the phase accumulated wave 208 from the detector proximate surface 206 to the radiation source 201; and forward propagating, in the forward propagation direction 207 from the radiation source 201 toward the detector 205 in absence of the sample 203, the phase accumulated wave 208 while treating Fresnel diffraction according to

$U\left(X_b,Y_b,Z_b\right)=\frac{e^{{\mathrm{ikZ}}_{\mathrm{ba}}}}{i\lambda Z_{\mathrm{ba}}}\int {\mathrm{dX}}_a{\mathrm{dY}}_a\exp \left\{i\frac{k}{2Z_{\mathrm{ba}}}\left[{\left(X_b-X_a\right)}^2+{\left(Y_b-Y_a\right)}^2\right]\right\}\times U\left(X_a,Y_a,Z_a\right),,$

such that the empirical image 211 of the sample 203 is reconstructed by projections and diffraction via maximum likelihood, a Bayesian prior probability distribution, and a Fresnel propagator.

In an embodiment, forward propagating, in the forward propagation direction 207, the initial wave 202 through the sample 203 includes determining h_i from I_n, according to

$I_{j\psi }=\int \mathrm{dE}D\left(E\right)I_j^{\left(0\right)}\left(E\right)\exp \left(-\sum _i{\alpha }_i\left(E\right)\sum _{\overrightarrow{r}}{f}_{\overrightarrow{r}i}{A}_{\overrightarrow{r}\psi }\right).$

In an embodiment, forward propagation, in the forward propagation direction 207, the initial wave 202 through the sample 203 includes minimizing an objective function the difference between the scalar diffraction equation within the Fresnel approximation and the empirical image 211. In an embodiment, the objective function is)

L _MAP({right arrow over (n)}|{right arrow over (f)})=L_ML({right arrow over (n)}|{right arrow over (f)})+g({right arrow over (f)}).

In an embodiment, the prior probability distribution is

$g\left(\overrightarrow{f}\right)=\sum _{\langle \overrightarrow{r}{\overrightarrow{r}}^{\prime })}{c}_{\overrightarrow{r}{\overrightarrow{r}}^{\prime }}{❘f_{\overrightarrow{r}}-f_{{\overrightarrow{r}}^{\prime }}❘}^p.$

In an embodiment, forward propagation further includes maximizing the objective function; determining the gradient of the objective function according to the equation

$\frac{\partial I_{j\psi }}{\partial f_{\overrightarrow{r}i}}=-\int \mathrm{dE}D\left(E\right)I_j^{\left(0\right)}\left(E\right){\alpha }_i\left(E\right)A_{\overrightarrow{r}\psi }\exp \left(-\sum _{\overrightarrow{r}\overrightarrow{r}}{\alpha }_{\overrightarrow{r}}\left(E\right)f_{\overrightarrow{r}i}{A}_{\overrightarrow{r}\psi }\right);$

and subjecting values of the objective function and the gradient of the objective function to a BFGS algorithm.

In an embodiment, forward propagating, in the forward propagation direction 207 from the radiation source 201 toward the detector 205 in absence of the sample 203, further includes determining the scalar wave according to

$\tilde{U}\left(X_2,Y_2,Z_2\right)=\frac{1}{{\lambda }^3{Z}_{10}^2{Z}_{20}}\int \int {\mathrm{dX}}_0{\mathrm{dY}}_0\exp \left\{i\frac{k}{2}\left(\frac{1}{Z_{20}}-\frac{1}{Z_{10}}\right)\left(X_0^2+Y_0^2\right)\right\}\times B\left(X_0,Y_0\right)\exp \left\{-i\frac{k}{Z_{20}}\left(X_2{X}_0+Y_2{Y}_0\right)\right\}.$

In an embodiment, forward propagating, in the forward propagation direction 207 from the radiation source 201 toward the detector 205 in absence of the sample 203, further includes determining the intensity according to

I(X₂,Y₂,Z₂)=|U(X₂,Y₂,Z₂)|²=|Ũ(X₂,Y₂,Z₂)|².

In an embodiment, a process for removing diffraction effects in a tomographic image includes: determining b(X₁, Y₁) from a plurality of projections and an estimate of the material parameters δ and β; applying the Fourier transform to b(X₁, Y₁) to find B(X₀, Y₀); determining Ũ(X₂, Y₂, Z₂) by multiplying B(X₀, Y₀) by a quadratic phasor and prefactors and applying an inverse Fourier transform according to

$\tilde{U}\left(X_2,Y_2,Z_2\right)=\frac{1}{{\lambda }^3{Z}_{10}^2{Z}_{20}}\int \int {\mathrm{dX}}_0{\mathrm{dY}}_0\exp \left\{i\frac{k}{2}\left(\frac{1}{Z_{20}}-\frac{1}{Z_{10}}\right)\left(X_0^2+Y_0^2\right)\right\}\times B\left(X_0,Y_0\right)\exp \left\{-i\frac{k}{Z_{20}}\left(X_2{X}_0+Y_2{Y}_0\right)\right\}$

In an embodiment, the process further includes performing the previous steps for each photon energy in a spectrum of energies of radiation that produces the tomographic image. In an embodiment, process further includes looping over viewing angles of the detector 205. In an embodiment, the process further includes, for each of the viewing angles: determining projections of the initial wave 202 through the sample 203; and determining the scalar wave U₂ of the initial wave 202 on the detector 205. In an embodiment, the process further includes looping over the projections of the initial wave 202. In an embodiment, the process further includes, for each of the projections: determining ∂U₁/∂P_{{right arrow over (s)}1i}; determining ∂U₂/∂P_{{right arrow over (s)}1i}; determining ∂I_{j{right arrow over (s)}2v}/∂P_{{right arrow over (s)}1i}; and determining ∂L/∂P_{{right arrow over (s)}1i}. In an embodiment, the process further includes multiplying ∂L/∂P_{{right arrow over (s)}1i} by the system matrix and accumulating terms in a gradient along a projection.

In an embodiment, determining the scalar wave U₂ of the initial wave 202 on the detector 205 occurs according to

$U_2\left({\overrightarrow{s}}_2,E,v\right)=\sum _{{\overrightarrow{s}}_1}G\left({\overrightarrow{s}}_2,{\overrightarrow{s}}_1,E\right)U_1\left({\overrightarrow{s}}_1,E,v\right).$

In an embodiment, the radiation source 201 produces the initial wave 202 such that the initial wave 202 is partially coherent.

In an embodiment, the radiation source 201 is an x-ray tube although other sources such as electron sources can be used.

The process for removing diffraction effects in a tomographic image is illustrated further by the following Example, which is non-limiting.

EXAMPLE

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

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. Although herein considers non-magnetic systems, vector x-ray nanotomography has been used to image magnetic systems on the micrometer scale where polarization dependence is important. Here, the x-ray wavelength, typically 100 μm, is small compared to the minimum feature size, which is near 1 μm. Hence, scalar diffraction theory is sufficient. Moreover, considered is the case where the width of the detector is small compared to the effective sample-detector distance, so the Fresnel propagator can be used.

Projective Tomography.

Code can perform projective tomography using a multi-spectral source. The process herein builds on projective tomography as presented herein. The multi-spectral, multi-material notation 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, a single spectrum and single basis material are used although physically more than one material can be 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

$\begin{array}{cc}{f}_{j\psi }=\int \mathrm{dE}D\left(E\right)I_j^{\left(0\right)}\left(E\right)\exp \left(-\sum _i{\alpha }_i\left(E\right)\sum _{\overrightarrow{r}}{f}_{\overrightarrow{r}i}{A}_{j\psi }\right)& \left(1\right)\end{array}$

where f_{{right arrow over (r)}i} is the concentration of basis material i at voxel {right arrow over (r)}, D(E) is the detection efficiency at photon energy E, and A_{{right arrow over (r)}ψ} is the system matrix which consists of the lengths of line segments travelling through voxel {right arrow over (r)} for viewing condition ψ. The viewing condition ψ is a joint set of viewing angles and projections. There is one projection per detector pixel. Continuing, I_j⁽⁰⁾ (E) 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_{{right arrow over (r)}i} for each voxel at {right arrow over (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, and used are those of the Center for X-ray Optics (CXRO) that also offers the complex index of refraction, which is used for Fresnel diffraction.

The projection integral

$\begin{array}{cc}{P}_{i\psi }=\sum _{\overrightarrow{r}}{f}_{\overrightarrow{r}i}{A}_{\overrightarrow{r}\psi }& \left(2\right)\end{array}$

appears in Eq. (1). If there are N×N×N voxels in the reconstruction, then A_{{right arrow over (r)}ψ} can have no more than 3N non-zero values for a given ψ. In principle, the dimensions of A are N³ by the product of the number of detector pixels times the number of viewing angles. The number of detector pixels is expected to be O(N²), and the number of viewing angles is expected to be O(N). For example, for Nyquist-limited sampling in 2D projection tomography, the number of viewing angles should be π/2N or greater. Hence A has O(N⁶) elements of which O(N⁴) are non-zero. This matrix is too large to store, so it is computed as needed.

The best solution is obtained 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, the function is

L _MAP({right arrow over (n)}|{right arrow over (f)})=L_ML({right arrow over (n)}|{right arrow over (f)})+g({right arrow over (f)}),  (3)

where {right arrow over (n)} has elements n_J, which are the counts observed at the joint index of observation J=(jψ), {right arrow over (f)} is the proposed reconstruction whose components are the f_{{right arrow over (r)}i} introduced above, and g({right arrow over (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:

$\begin{array}{cc}{L}_{\mathrm{ML}}\left(\overrightarrow{n}❘\overrightarrow{f}\right)=\sum _j\left[\ln n_J!-n_J\ln I_J\left(\overrightarrow{f}\right)+I_J\left(\overrightarrow{f}\right)\right].& \left(4\right)\end{array}$

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. The prior probability distribution is

$\begin{array}{cc}g\left(\overrightarrow{f}\right)=\sum _{\langle \overrightarrow{r}\overrightarrow{r}\rangle }{c}_{\overrightarrow{r}\overrightarrow{r}}{❘f_{\overrightarrow{r}}-f_{{\overrightarrow{r}}^{\prime }}❘}^p& \left(5\right)\end{array}$

where the sum is taken over neighboring pairs. The process specializes to isotropic cubic voxels where the neighboring pairs include six faces with C_{{right arrow over (r)}{right arrow over (r)}′}=1, twelve edges with c_{{right arrow over (r)}{right arrow over (r)}′}=1/√{square root over (2)}, and eight vertices with c_{{right arrow over (r)}{right arrow over (r)}}′=1/√{square root over (3)}.

For the case of 1<p≤2, one can use 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. Herein is adopted p=1.1. 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.

The objective function is maximized by initializing the f_{{right arrow over (r)}i} with random numbers and applying the L-BFGS-B code of Ref. (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:

$\begin{array}{cc}\frac{\partial I_{j\psi }}{\partial f_{\overrightarrow{r}i}}=-\int \mathrm{dE}D\left(E\right)I_j^{\left(0\right)}\left(E\right){\alpha }_i\left(E\right)A_{\overrightarrow{r}\psi }\exp \left(-\sum _{\overrightarrow{r}\overrightarrow{r}}{\alpha }_{i^{\prime }}\left(E\right)f_{\overrightarrow{r}\overrightarrow{r}}{A}_{\overrightarrow{r}\psi }\right).& \left(6\right)\end{array}$

A constraint is that each material in each voxel is non-negative, i.e., f_{{right arrow over (r)}i}≥0. The BFGS algorithm is used with a 128 dimensional subspace, a value which was near optimal in a study with the program on a different example. This figure is more than the range of 3 to 20 recommended by Ref. If the Bayesian term g({right arrow over (f)}) is included, 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, treat Fresnel diffraction as follows: accumulate changes in both phase and amplitude on projections through the sample. Neglect diffraction within the sample itself but include it when considering propagation from the sample to the detector. The equation is

U(X,Y,Z_f)≈exp{ikZ_fi−k∫_Zi^ZfdZ[iδ_E(X,Y,Z)+β_E(X,Y,Z)]}U(X,Y,Z_i),  (7)

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

n _E=1−δ_E+iβ_E.  (8)

The wavevector

$k=\frac{2\pi }{\lambda },$

is related to the photon energy by E=ℏck where ℏ is the reduced Planck constant and c is the speed of light. Below, tabulated values of the complex index of refraction for silica is 2.2 g/cm³. The distinction between the index of refraction n(E) and the observations m should be clear from context.

Fresnel propagator.

Assume that there is a point source for the x rays and 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. The spherical wave is 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. 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. However, the instant implementation does not invoke the theorem explicitly or make use of transformed variables.

The Fresnel propagator that takes the scalar amplitude from a plane Z=Z_a to a plane Z=Z_b is given by

$\begin{array}{cc}U\left(X_b,Y_b,Z_b\right)=\frac{e^{{\mathrm{ikZ}}_{\mathrm{ba}}}}{i\lambda Z_{\mathrm{ba}}}\int \int {\mathrm{dX}}_a{\mathrm{dY}}_a\exp \left\{i{\frac{k}{2Z_{\mathrm{ba}}}\left[{\left(X_b-X_a\right)}^2+\left(Y_b-Y_a\right)\right]}^2\right\}\times U\left(X_a,Y_a,Z_a\right),& \left(9\right)\end{array}$

where (X_a, Y_a) are Cartesian coordinates in the Z=Z_a plane with a similar relation for (X_b, Y_b), Z_ba=Z_b−Z_a, and the domain of integration is the whole Z=Z_a plane. U is energy-dependent, but the variable E is suppressed to keep the formulas readable.

The diffraction calculation proceeds with these planes: Z=Z₀ the plane including the point source, Z=Z₁ the mid-point of the sample, and Z=Z₂ 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₁. In the Fresnel propagation portion of the calculation, the projected phase is assigned to the Z=Z₁ plane.

Starting from a point source of unit integrated strength, taken to be at (X₀=0, Y₀=0, Z₀), the scalar wave on the Z=Z₁ plane, without considering the effect of the sample, is

$\begin{array}{cc}{U}^{\left(0\right)}\left(X_1,Y_1,Z_1\right)=\frac{e^{{\mathrm{ikZ}}_{10}}}{i\lambda Z_{10}}\exp \left\{i\frac{k}{2Z_{10}}\left(X_1^2+Y_1^2\right)\right\}.& \left(10\right)\end{array}$

The effect of the sample treated in the projection approximation is to introduce a phasor via

U(X₁,Y₁,Z₁)=b(X₁,Y₁)U⁽⁰⁾(X₁,Y₁,Z₁)  (11)

with

b(X,Y)=exp{−∫_Zi^ZfdZ[ikδ_E(X,Y,Z)+kβ_E(X,Y,Z)]}.  (12)

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, this approximation is arrived at by compressing the sample into a small strip while preserving the projected mass. In this implementation, 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, the scalar wave is

$\begin{array}{cc}\begin{array}{c}U\left(X_0,Y_0,Z_0\right)=\frac{e^{{\mathrm{ikZ}}_{01}}}{i\lambda Z_{01}}\int \int {\mathrm{dX}}_1{\mathrm{dY}}_1\\ \exp \left\{i\frac{k}{2Z_{01}}\left[{\left(X_0-X_1\right)}^2+{\left(Y_0-Y_1\right)}^2\right]\right\}U\left(X_1,Y_1,Z_1\right)\\ =\frac{1}{{\lambda }^2{Z}_{10}^2}\exp \left\{-i\frac{k}{2Z_{10}}\left(X_0^2+Y_0^2\right)\right\}\times \\ \int \int {\mathrm{dX}}_1{\mathrm{dY}}_{1}b\left(X_1,Y_1\right)\exp \left\{i\frac{k}{Z_{10}}\left(X_0{X}_1+Y_0{Y}_1\right)\right\}.\end{array}& \left(13\right)\end{array}$

In going from line 1 to line 2 in Eq. (13), the quadratic phase factor in (X₁, Y₁) cancels. The bandwidth of b(X₁, Y₁) is comparable to functions that arise in projection tomography

For the wave on the detector, propagate forward from the source plane Z=Z₀ to the detector plane Z=Z₂ resulting in

$\begin{array}{cc}U\left(X_2,Y_2,Z_2\right)=\frac{1}{{\lambda }^2{Z}_{10}^2}\frac{1}{i\lambda Z_{10}}\exp \left\{{\mathrm{ikZ}}_{20}+i\frac{k}{2Z_{20}}\left(X_2^2+Y_2^2\right)\right\}\times \int \int {\mathrm{dX}}_0{\mathrm{dY}}_0\exp \left\{i\frac{k}{2}\left(\frac{1}{Z_{20}}-\frac{1}{Z_{10}}\right)\left(X_0^2+Y_0^2\right)\right\}\times B\left(X_0,Y_0\right)\exp \left\{-i\frac{k}{Z_{20}}\left(X_2{X}_0+Y_2{Y}_0\right)\right\}& \left(14\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}B\left(X_0,Y_0\right)=\int \int {\mathrm{dX}}_1{\mathrm{dY}}_{1}b\left(X_1,Y_1\right)\exp \left\{i\frac{k}{Z_{10}}\left(X_0{X}_1+Y_0{Y}_1\right)\right\}.& \left(15\right)\end{array}$

Equation (14) describes free space propagation performed as if the sample is not present. Physically, the sample is already accounted through b(X₁, Y₁). Computationally, the scalar wave is determined

$\begin{array}{cc}\tilde{U}\left(X_2,Y_2,Z_2\right)=\frac{1}{{\lambda }^3{Z}_{10}^2{Z}_{20}}\int \int {\mathrm{dX}}_0{\mathrm{dY}}_0\exp \left\{i\frac{k}{2}\left(\frac{1}{Z_{20}}-\frac{1}{Z_{10}}\right)\left(X_0^2+Y_0^2\right)\right\}\times B\left(X_0,Y_0\right)\exp \left\{-i\frac{k}{Z_{20}}\left(X_2{X}_0+Y_2{Y}_0\right)\right\}& \left(16\right)\end{array}$

and then find the intensity by squaring the amplitude, i.e.,

I(X₂,Y₂,Z₂)=|U(X₂,Y₂,Z₂)|²=|Ũ(X₂,Y₂,Z₂)|².  (17)

The only difference between U(X₂, Y₂, Z₂) and U (X₂, Y₂, Z₂) is a phasor which does not appear in the intensity I(X₂, Y₂, Z₂).

Algorithmically, perform the following:

• • 1. find b(X₁, Y₁) from the projections and an estimate of the material parameters δ and β at any point,
  • 2. apply the Fourier transform to find B(X₀, Y₀),
  • 3. multiply by the quadratic phasor in Eq. (16), then
  • 4. apply the inverse Fourier transform to find Ũ (X₂, Y₂, Z₂) after including the prefactors in Eq. (16);
  • 5. these steps are performed for each photon energy in the spectrum.

The pixels sizes Δ_i on the source, sample, and detector planes, for i=0, 1, 2 respectively, obey NΔ₀Δ₁=λZ₁₀, NΔ₀Δ₂=λZ₂₀, Δ₂=MΔ₁, and M=Z₂₀/Z₁₀. The geometric magnification M of the pixels is consistent with the Fourier scaling theorem.

Gradient of likelihood for the Fresnel propagator.

The code relies on the BFGS algorithm 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. This has been achieved previously for projective tomography. In the Fresnel case, the calculation is similar but assume that the influence of a given projection on the detector affects O(N) of the N² detector pixels. The analytic example makes the assumption plausible. Hence, there are O(N³) operations per view, just as for the projective case. In the projective case, the O(N³) operations arise from making O(N²) 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 ψ. Here, further decompose ψ=({right arrow over (s)}₂, v) where {right arrow over (s)}₂ is a detector pixel and v is a viewing angle. For the Fresnel case,

$\begin{array}{cc}{I}_{j{\overrightarrow{s}}_{2}v}=\sum _E{I}_{\mathrm{jE}}^{\left(0\right)}{❘U_s\left({\overrightarrow{s}}_2,E,v\right)❘}^2.& \left(18\right)\end{array}$

This equation is similar to Eq. (1) in the projective case. In turn

$\begin{array}{cc}{U}_2\left({\overrightarrow{s}}_2,E,v\right)=\sum _{{\overrightarrow{s}}_1}G\left({\overrightarrow{s}}_2,{\overrightarrow{s}}_1,E\right)U_1\left({\overrightarrow{s}}_1,E,v\right).& \left(19\right)\end{array}$

Here, U₂ and U₁ 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

$\begin{array}{cc}{U}_1\left({\overrightarrow{s}}_1,E,v\right)=\exp \left\{\sum _i{\hat{\alpha }}_i\left(E\right)P_{{\overrightarrow{s}}_1\mathrm{vi}}\right\}{U}_1^{\left(0\right)}\left({\overrightarrow{s}}_1,E,v\right)& \left(20\right)\end{array}$

where P_{{right arrow over (s)}1vi} is a projection and U₁⁽⁰⁾ (s1, E, v) is the scalar wave at the sample plane if the sample were absent and

α_i(E)=−k[β_i(E)+iδ_i(E)].  (21)

The projection is given by Eq. (2), with the complex {tilde over (α)}_i(E) being analogous to the real {tilde over (α)}_i(E).

BFGS requires the derivatives of Eq. (3), which in turn requires the derivatives of Eq. (18).

$\begin{array}{cc}\begin{array}{c}\frac{\partial I_{j{\overrightarrow{s}}_{2}v}}{\partial f_{\overrightarrow{r}i}}=\sum _E{I}_{\mathrm{jE}}^{\left(0\right)}\frac{\partial }{\partial f_{\overrightarrow{r}i}}{❘U_2\left({\overrightarrow{s}}_2,E,v\right)❘}^2\\ =2\sum _E{I}_{\mathrm{jE}}^{\left(0\right)}\mathrm{Re}\left[U_2^{*}\left({\overrightarrow{s}}_2,E,v\right)\frac{\partial }{\partial f_{\overrightarrow{r}i}}{U}_2\left({\overrightarrow{s}}_2,E,v\right)\right].\end{array}& \left(22\right)\end{array}$

Equation (22) requires the following expression

$\begin{array}{cc}\frac{\partial }{\partial f_{\overrightarrow{r}i}}{U}_2\left({\overrightarrow{s}}_2,E,v\right)=\sum _{{\overrightarrow{s}}_1}G\left({\overrightarrow{s}}_2,{\overrightarrow{s}}_1,E\right)\frac{\partial }{\partial f_{\overrightarrow{r}i}}{U}_1\left({\overrightarrow{s}}_1,E,v\right).& \left(23\right)\end{array}$

Equation (23) requires the derivative of Eq. (20) namely

$\begin{array}{cc}\begin{array}{c}\frac{\partial }{\partial f_{\overrightarrow{r}i}}{U}_1\left({\overrightarrow{s}}_2,E,v\right)=\left[{\tilde{\alpha }}_i\left(E\right)\frac{\partial P_{{\overrightarrow{s}}_1{r}_i}}{\partial f_{\overrightarrow{r}i}}\right]\exp \left\{\sum _i{\tilde{\alpha }}_i\left(E\right)\text{?}\right\}\\ U_1^{\left(0\right)}\left({\overrightarrow{s}}_1,E,v\right)\\ ={\tilde{\alpha }}_i\left(E\right)A_{\overrightarrow{r}{s}_{1}v}{U}_1\left({\overrightarrow{s}}_1,E,v\right)\end{array}\text{?}\text{indicates text missing or illegible when filed}& \left(24\right)\end{array}$

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 as

$\begin{array}{cc}\frac{\partial L_{\mathrm{ML}}}{\partial f_{\overrightarrow{r}i}}=\frac{\partial L_{\mathrm{ML}}}{\partial \text{?}}{A}_{\overrightarrow{r}{s}_{1}v}.\text{?}\text{indicates text missing or illegible when filed}& \left(25\right)\end{array}$

Equation (25) implies a large computational savings since there are O(N) more instances of

$\frac{\partial L_{\mathrm{ML}}}{\partial f_{\overrightarrow{r}i}}$

than

$\frac{\partial P_{{\overrightarrow{r}}_1\mathrm{vi}}}{\partial f_{\overrightarrow{r}i}}$

and A{right arrow over ( )}_rs1v is very easy to compute.

In the case of partial coherence, approximate G({right arrow over (s)}₂, {right arrow over (s)}₁, 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. Refer to the region of G({right arrow over (s)}₂, {right arrow over (s)}₁, E), which is represented 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 {right arrow over (s)}₂, loop over a set of virtual pixels on the sample denoted by si. In practice the points {right arrow over (s)}₁ and {right arrow over (s)}₂ will be in 1:1 correspondence, i.e., use one projection per detector pixel. Use the fast Fourier transforms of FFTW 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₂
  • Loop on projections {
    • Find ∂U₁/∂P_{{right arrow over (s)}1i}
    • Find ∂U₂/∂P_{{right arrow over (s)}1i over the non-negligible region of the Green's function}
    • Find ∂I_{j{right arrow over (s)}2v}/∂P_{{right arrow over (s)}1i}
    • Find ∂L/∂P_{{right arrow over (s)}1i}
    • 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({right arrow over (s)}₂, {right arrow over (s)}₁, E, v) has order N non-zero entries where there are N² pixels. This is because one needs O(N) steps to find each projection, so one can afford to do O(N) operations thereafter. Overall scaling is as O(N⁴), taking into account O(N²) projections per view, O(N) operations per projection, and O(N) view points. The projective case also scales as O(N⁴) albeit with a smaller coefficient. The projective case can 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, start the diffractive solution from the projective solution. In testing using random dot patterns, it was found that the diffractive program could has a finite basin of attraction when starting from random numbers.

Analytic Example of Partial Coherence.

For tomography, there is a need to find the complex index of refraction at each voxel in some defined portion of space. In order to use BFGS, one needs the gradient of the diffraction pattern with respect to changes in the index of refraction at a given voxel. Here, an analytic example is used 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 there is a single voxel that differs from the vacuum value by a differential amount. 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. The voxel is located at (x₁, y₁, z₁). The detector plane is located in the plane z=z₂. The unscattered wave is denoted by U and the scattered wave by U⁽¹⁾ and is given by

$\begin{array}{cc}{U}^{\left(1\right)}\left(x_2,y_2,z_2;x_1,y_1,z_1\right)=\frac{U_0{U}_1}{z_1{z}_{21}}\exp \left[i\frac{k}{2z_1}\left(x_1^2+y_1^2\right)\right]\times \exp \left\{i\frac{k}{2z_{21}}\left[{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2\right]\right\}& \left(26\right)\end{array}$

which may be found by considering Eq. (9) with a point source.

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

$\begin{array}{cc}\begin{array}{c}I\left(k,x_2,y_2;x_1,y_1\right)=2\mathrm{Re}{U}^{*}\left(x_2,y_2,z_2\right)U^{\left(1\right)}\left(x_2,y_2,z_2;x_1,y_1,z_1\right)\\ =C_1\mathrm{Re}{U}_1\times \\ \exp \left\{i\frac{k}{2}\left[\begin{array}{c}\frac{{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2}{z_2-z_1}-\\ \frac{x_2^2+y_2^2}{z_2}+\frac{x_1^2+y_1^2}{z_1}\end{array}\right]\right\},\end{array}& \left(27\right)\end{array}$

which is a first order expansion of Eq. (17). Next, assume the source intensity is normally distributed in k, hence also normally distributed in the photon energy. One wants to find

$\begin{array}{cc}\begin{array}{c}I\left(x_2,y_2;x_1,y_1\right)={\int }_{-\infty }^{\infty }\mathrm{dkI}\left(k,x_2,y_2;x_1,y_1\right)\\ =\frac{C_1}{\sqrt{2\pi }{\sigma }_k}\mathrm{Re}{U}_1{\int }_{-\infty }^{\infty }\mathrm{dk}\exp \left[-\frac{{\left(k-k_0\right)}^2}{2{\sigma }_k^2}+\mathrm{ik}\xi \right]\\ =C_1\exp \left(-\frac{{\sigma }_k^2{\xi }^2}{2}\right)\mathrm{Re}{U}_1\exp \left({\mathrm{ik}}_0\xi \right)\end{array}& \left(28\right)\end{array}$

where ξ is half of the term in brackets in Eq. (27).

To understand the Gaussian envelope, rearrange as

$\begin{array}{cc}\begin{array}{c}\xi =\frac{1}{2}\left[\begin{array}{c}\frac{{\left(x_2-x_1\right)}^2}{z_{21}}-\frac{x_2^2}{z_2}+\frac{x_1^2}{z_1}+\\ \frac{{\left(y_2-y_1\right)}^2}{z_{21}}-\frac{y_2^2}{z_2}+\frac{y_1^2}{z_1}\end{array}\right]\\ =\frac{{\left(x_2-M{x}_1\right)}^2+{\left(y_2-M{y}_1\right)}^2}{2{\mathrm{Mz}}_{21}}.\end{array}& \left(29\right)\end{array}$

The result shows the diffraction pattern is centered on the geometric projection of (x₁, y₁) onto the z₂ 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. 3.

Empirical Results.

A graded index optical fiber is studied. 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 showed that 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, visible light microscopy and scanning electron microscopy (SEM) were performed 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 50× 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. 4. 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, Calif., 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₁=16 mm and the sample to detector distance was Z₂₁=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 20× was used. A tilt series was acquired 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. 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 y′(α) of the fiber to that of a projection of a rigid cylinder and correcting for global tilt using

γ′(α)=tan⁻¹[cos(α−α₀)tan γ]+β,  (30)

where α is the angle of the tomographic projection, α₀ 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. 5. Both α and α₀ denote angles in the same space as the angles of the tilt series. The angles α and α₀ are unrelated to a above; similarly, {circumflex over (r)} and {circumflex over (r)}′(α) have no relation to the symbols r and denoting the voxels.

The tilt series images of the optical fiber were then further aligned using the program Arec3d. 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. 6. The digits in parentheses represent one standard error in the fit. The angle between the fiber and the reconstruction axis from an analysis of the reconstructed stack was 3.87(14°), a value which should be and is consistent with y.

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; 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 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. Take D(E)=1, i.e., detector efficiency is assumed to be ideal. 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.

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 the assumption of Poisson statistics in Eq. (4). To make contact with the formalism, reconstruct with a single basis material, so i=1, and use a single spectrum, so j=1, with 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 μm, 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, express the range of Fresnel numbers for the Gaussian model spectrum using uncertainty notation as F=0.91(18). Since F is near 1, diffraction effects can be a moderate correction to geometric optics.

Results are shown in FIG. 7. 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. 7(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. FIG. 7(c) is a reconstruction using maximum likelihood, but including the effect of Fresnel diffraction. Compared to FIG. 7(a), feature 1, at the cladding-coating boundary and feature 2 at the coating-air boundary are edge diffraction artifacts (in FIG. 7(a)) which are greatly reduced when diffraction is included in the physical model of reconstruction as shown in FIG. 7(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. 7(c). When a Bayesian prior is used, the stipple is eliminated in the projective maximum likelihood reconstruction of FIG. 7(b) as well as the Bayesian/diffraction reconstruction, shown in FIG. 7(d). The stipple is better seen in FIG. 8, a blow-up of the central region.

Cylindrical averages are shown in FIG. 9. In FIG. 9(a), the cylinders were defined as 1 pixel wide rings about the center of the cladding on the 8 central slices. In FIG. 9(a), the reconstructions are seen to give about the same continuous decrease in the core region, followed by a gentler decrease in the cladding. The core region met expectations from the manufacturer's specification, but the gentler 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. 9(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. 9(b) occurs just past the end of FIG. 9(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. 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. 4(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 FIG. 10. 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 FIG. 10. 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. The results are not strictly comparable because of the structure of the starting points. The Fresnel calculations required more iterations. Overall, estimate a one order of magnitude penalty for the Fresnel calculation compared to its projective counterpart.

The algorithm treats 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. 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. The boundary with both visible light microscopy and scanning electron microscopy was observed, confirming a conjecture made earlier.

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. The algorithm presented here can account for such features with a computational time which scales like projective tomography.

The processes described herein may be embodied in, and fully automated via, software code modules executed by a computing system that includes one or more general purpose computers or processors. The code modules may be stored in any type of non-transitory computer-readable medium or other computer storage device. Some or all the methods may alternatively be embodied in specialized computer hardware. In addition, the components referred to herein may be implemented in hardware, software, firmware, or a combination thereof.

Many other variations than those described herein will be apparent from this disclosure. For example, depending on the embodiment, certain acts, events, or functions of any of the algorithms described herein can be performed in a different sequence, can be added, merged, or left out altogether (e.g., not all described acts or events are necessary for the practice of the algorithms). Moreover, in certain embodiments, acts or events can be performed concurrently, e.g., through multi-threaded processing, interrupt processing, or multiple processors or processor cores or on other parallel architectures, rather than sequentially. In addition, different tasks or processes can be performed by different machines and/or computing systems that can function together.

Any logical blocks, modules, and algorithm elements described or used in connection with the embodiments disclosed herein can be implemented as electronic hardware, computer software, or combinations of both. To clearly illustrate this interchangeability of hardware and software, various illustrative components, blocks, modules, and elements have been described above generally in terms of their functionality. Whether such functionality is implemented as hardware or software depends upon the particular application and design constraints imposed on the overall system. The described functionality can be implemented in varying ways for each particular application, but such implementation decisions should not be interpreted as causing a departure from the scope of the disclosure.

The various illustrative logical blocks and modules described or used in connection with the embodiments disclosed herein can be implemented or performed by a machine, such as a processing unit or processor, a digital signal processor (DSP), an application specific integrated circuit (ASIC), a field programmable gate array (FPGA) or other programmable logic device, discrete gate or transistor logic, discrete hardware components, or any combination thereof designed to perform the functions described herein. A processor can be a microprocessor, but in the alternative, the processor can be a controller, microcontroller, or state machine, combinations of the same, or the like. A processor can include electrical circuitry configured to process computer-executable instructions. In another embodiment, a processor includes an FPGA or other programmable device that performs logic operations without processing computer-executable instructions. A processor can also be implemented as a combination of computing devices, e.g., a combination of a DSP and a microprocessor, a plurality of microprocessors, one or more microprocessors in conjunction with a DSP core, or any other such configuration. Although described herein primarily with respect to digital technology, a processor may also include primarily analog components. For example, some or all of the signal processing algorithms described herein may be implemented in analog circuitry or mixed analog and digital circuitry. A computing environment can include any type of computer system, including, but not limited to, a computer system based on a microprocessor, a mainframe computer, a digital signal processor, a portable computing device, a device controller, or a computational engine within an appliance, to name a few.

The elements of a method, process, or algorithm described in connection with the embodiments disclosed herein can be embodied directly in hardware, in a software module stored in one or more memory devices and executed by one or more processors, or in a combination of the two. A software module can reside in RAM memory, flash memory, ROM memory, EPROM memory, EEPROM memory, registers, hard disk, a removable disk, a CD-ROM, or any other form of non-transitory computer-readable storage medium, media, or physical computer storage known in the art. An example storage medium can be coupled to the processor such that the processor can read information from, and write information to, the storage medium. In the alternative, the storage medium can be integral to the processor. The storage medium can be volatile or nonvolatile.

While one or more embodiments have been shown and described, modifications and substitutions may be made thereto without departing from the spirit and scope of the invention. Accordingly, it is to be understood that the present invention has been described by way of illustrations and not limitation. Embodiments herein can be used independently or can be combined.

All ranges disclosed herein are inclusive of the endpoints, and the endpoints are independently combinable with each other. The ranges are continuous and thus contain every value and subset thereof in the range. Unless otherwise stated or contextually inapplicable, all percentages, when expressing a quantity, are weight percentages. The suffix (s) as used herein is intended to include both the singular and the plural of the term that it modifies, thereby including at least one of that term (e.g., the colorant(s) includes at least one colorants). Option, optional, or optionally means that the subsequently described event or circumstance can or cannot occur, and that the description includes instances where the event occurs and instances where it does not. As used herein, combination is inclusive of blends, mixtures, alloys, reaction products, collection of elements, and the like.

As used herein, a combination thereof refers to a combination comprising at least one of the named constituents, components, compounds, or elements, optionally together with one or more of the same class of constituents, components, compounds, or elements.

All references are incorporated herein by reference.

The use of the terms “a,” “an,” and “the” and similar referents in the context of describing the invention (especially in the context of the following claims) are to be construed to cover both the singular and the plural, unless otherwise indicated herein or clearly contradicted by context. It can further be noted that the terms first, second, primary, secondary, and the like herein do not denote any order, quantity, or importance, but rather are used to distinguish one element from another. It will also be understood that, although the terms first, second, etc. are, in some instances, used herein to describe various elements, these elements should not be limited by these terms. For example, a first current could be termed a second current, and, similarly, a second current could be termed a first current, without departing from the scope of the various described embodiments. The first current and the second current are both currents, but they are not the same condition unless explicitly stated as such.

The modifier about used in connection with a quantity is inclusive of the stated value and has the meaning dictated by the context (e.g., it includes the degree of error associated with measurement of the particular quantity). The conjunction or is used to link objects of a list or alternatives and is not disjunctive; rather the elements can be used separately or can be combined together under appropriate circumstances.

Claims

1. A process for removing diffraction effects in a tomographic image, the process comprising: obtaining an empirical image of a sample, the empirical image comprising diffraction effects;producing an initial wave at a radiation source;forward propagating, in a forward propagation direction from the radiation source toward a detector;receiving, at a radiation source proximate surface of the sample that is interposed between the radiation source and the detector, the initial wave from the radiation source;forward propagating, in the forward propagation direction, the initial wave through the sample and accumulating phase and amplitude information of the initial wave as the initial wave propagates through the sample toward the detector to produce a phase accumulated wave at a detector proximate surface of the sample according to U(X,Y,Zf)≈exp{ikZf−k∫ZiZfdZ[iδE(X,Y,Z)+βE(X,Y,Z)]}U(X,Y,Zi);back propagating, in a reverse propagation direction in absence of the sample, the phase accumulated wave from the detector proximate surface to the radiation source; andforward propagating, in the forward propagation direction from the radiation source toward the detector in absence of the sample, the phase accumulated wave while treating Fresnel diffraction according to

2. The process of 1, wherein forward propagating, in the forward propagation direction, the initial wave through the sample comprises determining f{right arrow over (r)}i from Ijψ according to the equation

3. The process of claim 2, wherein forward propagation, in the forward propagation direction, the initial wave through the sample further comprises minimizing an objective function the difference between the scalar diffraction equation within the Fresnel approximation and the empirical image.

4. The process of claim 3, wherein the objective function is LMAP({right arrow over (n)}|{right arrow over (f)})=LML({right arrow over (n)}|{right arrow over (f)})+g({right arrow over (f)}).

5. The process of claim 4, wherein the prior probability distribution is

6. The process of claim 5, wherein forward propagation further comprises maximizing the objective function; determining the gradient of the objective function according to the equation

7. The process of claim 6, wherein forward propagating, in the forward propagation direction from the radiation source toward the detector in absence of the sample, further comprises determining the scalar wave according to

8. The process of claim 7, wherein forward propagating, in the forward propagation direction from the radiation source toward the detector in absence of the sample, further comprises determining the intensity according to I(X2,Y2,Z2)=|U(X2,Y2,Z2)|2=|Ũ(X2,Y2,Z2)|2.

9. The process of claim 1, further comprising looping over viewing angles of the detector.

10. The process of claim 9, further comprising, for each of the viewing angles: determining projections of the initial wave through the sample; anddetermining the scalar wave U2 of the initial wave on the detector.

11. The process of claim 10, further comprising looping over the projections of the initial wave.

12. The process of claim 11, further comprising, for each of the projections: determining ∂U1/∂P{right arrow over (s)}1i;determining ∂U2/∂P{right arrow over (s)}1i;determining ∝Ij{right arrow over (s)}2v/∂P{right arrow over (s)}1i; anddetermining ∂L/∂P{right arrow over (s)}1i.

13. The process of claim 12, further comprising multiplying ∂L/∂P{right arrow over (s)}1i by the system matrix and accumulating terms in a gradient along a projection.

14. The process of claim 9, wherein determining the scalar wave U2 of the initial wave on the detector occurs according to

15. The process of claim 1, wherein the radiation source produces the initial wave such that the initial wave is partially coherent.

16. The process of claim 1, wherein the radiation source comprises an x-ray tube.

17. A process for removing diffraction effects in a tomographic image, the process comprising: determining b(X1, Y1) from a plurality of projections and an estimate of the material parameters δ and β;applying the Fourier transform to b(X1, Y1) to find B(X0, Y0);determining Ũ(X2, Y2, Z2) by multiplying B(X0, Y0) by a quadratic phasor and prefactors and applying an inverse Fourier transform according to

18. The process of claim 17, further comprising performing the steps in claim 17 for each photon energy in a spectrum of energies of radiation that produces the tomographic image.