x-ray, phase imaging, computational imaging

Conventional x-ray imaging relies on differences in attenuation in a material to produce image contrast. While useful for differentiating structures with large density variations, the subtle differences in attenuation for low-density materials can be difficult to detect. X-ray phase imaging, on the other hand, relies on differences in phase delay which is typically several orders of magnitude larger than attenuation.

However, most methods of producing x-ray phase images rely on specialized synchrotron sources, small and low power microfocus sources or the careful alignment of several precision gratings. We demonstrate that focusing polycapillary optics can produce small focal spots from conventional x-ray sources to enable phase imaging. Moreover, in conjunction with focusing optics, the use of a simple, low cost wire mesh to structure the beam can significantly improve phase reconstructions.

When an x ray beam traverses an object, the interaction can be modeled through both attenuation, a loss of power due to absorption and scattering, and phase delay, an apparent advance of the waves making up the beam. Conventional x-ray imaging utilizes only differences in attenuation to produce image contrast, which results in poor contrast for low density materials. The phase delay accrued by x rays in traversing an object, on the other hand, can produce a signature that (for low-density materials) can be roughly 1000 times larger than attenuation at energies from 10-30 keV.

The impact of phase delay on x rays can be quantified through the complex index of refraction

Equation 1

The photoelectric attenuation of the incident beam’s intensity is then given by multiplication by a transmission factor

Equation 2

Where the attenuation 𝜇(𝐱) = ∫ 𝛽(𝐱, 𝑧)𝑑z is the integral of the complex part of the refractive index taken along ray trajectories through the object. The phase difference of x rays after exiting the sample compared with those which traverse the air outside the sample is

Equation 3

with the integral again being taken over ray paths through the sample.

Under a geometric approximation that the rays represent the direction of intensity flow, rays are deflected due to phase differences at the exit face of the sample such that 𝛉 = (2𝜋/𝜆 )∇𝜙 is a vector with units of radians representing the direction and angle of deflection due to phase. This leads to a method of phase retrieval know as propagation-based x-ray phase imaging (PBXPI), where the contact image is compared with the image after propagating some distance from the object.

Differences in the images can be attributed to ray deflection/phase gradient and then phase can be reconstructed by solving a differential equation. A schematic of a PBXPI measurement is illustrated in Figure 1.

Fig. 1 Sketch of propagation-based imaging. A small source (microfocus or the focal spot of a polycapillary optic) illuminates an object. Rays are attenuated and deflected an amount proportional to the phase gradient in the contact plane. The total deflection depends on the phase gradient multiplied by the propagation distance.

The mathematical expression describing PBXPI is known at the transport of intensity equation and is expressed as

Equation 4

Where
𝛻𝒙 is a gradient taken with respect to the contact plane coordinates,
𝑅₂ is the object-to-detector distance and
𝑀 = (𝑅₁+𝑅₂)/ 𝑅₁ is the magnification of the system
where
𝑅₁ is the source-to-object distance.

If both the contact plane image and detector-plane images are measured, and the distances and wavelength of the x rays are known, this reduces to a second order, linear partial differential equation that can be solved for the object phase. However, there are several practical considerations that make comparing these two images extremely difficult.

Deflection angles tend to be sub-milliradian, requiring that 𝑅₂ be on the order of tens of centimeters in order to achieve intensity variations that span the hundreds of microns necessary for measurement with high-resolution x ray cameras. To solve above equation without significant artifacts requires aligning the two images to the accuracy of a single pixel while accounting for magnification. This process is simply impractical in most cases. Therefore, several other methods have been proposed that utilize only a single image at a plane displaced from the object.

The weak attenuation model assumes that attenuation is negligible so that the contact image is uniform in intensity, so that 𝐼₀(𝐱) ≈ 𝐼₀ = constant. The phase then satisfies an equation of the form

Equation 5

where g denotes a combination of the measured, propagated intensity, constants, and the uniform contact image estimated from the background or mean value of the propagated image. The phase can be recovered by solving this Poisson’s equation using a variety of techniques. Here we apply a Fourier-transform-based, Poisson solver due to its speed and the ability to easily analyze the model’s shortcomings in the Fourier domain.

The Fourier domain solver takes advantage of the fact that the negative Laplacian operator is equivalent to multiplication by 𝐻(𝑢) = 4𝜋²𝑢² in the Fourier domain. This leads to a phase solution of

Equation 6

Artifacts in the reconstruction arise from the transfer function 𝐻⁻¹(𝑢) diverging as 𝑢 → 0 which means that low spatial frequency components of phase (those that correspond to larger spatial structures) must be very strong to produce similar components in the propagated intensity. Therefore, the bulk of the object attenuation will be reconstructed as extremely strong phase. This is commonly treated with regularization which removes low spatial frequencies below some user defined parameter 𝛼 from the reconstruction by using 𝐻−1(𝑢) = 𝑢²/(4𝜋²(𝑢⁴+𝛼⁴)). Figure 2 shows a comparison of 𝐻⁻¹ for WA and WA with Tikhonov regularization.

Fig. 2 Comparison of transfer functions for WA and WA with Tikhonov regularization.

Phase-attenuation duality (PAD) assumes that Compton scattering is the dominant process contributing to attenuation.

Under this assumption, both the phase and attenuation are proportional to electron density and therefore proportional to each other with proportionality constant 𝛾 that depends on the Compton scattering cross section.

This means that the contact image intensity can be written as a function of phase,

Equation 7

This allows Eq. (4) to be written as a differential equation for the contact image in terms of the (measured) propagated image

Equaation 8

where 𝐶 = (𝑅₂𝜆𝛾 )/ (4𝜋𝑀). This may be solved using Fourier transform methods

Equation 9

with 𝐻⁻¹(𝑢) = 1/(1+4𝜋²𝐶𝑢²) illustrated in Fig. 3.

Fig. 3 Shape of the transfer function for both the WA and PAD models.

Unlike the WA model, PAD’s transfer function does not diverge as 𝑢 → 0 which mitigates the need for regularization. The limitation of PAD is that the assumption that Compton scattering is dominant is only true at sufficiently high x-ray energies. The phase reconstruction can suffer significantly if these assumptions are violated.

The single material (SM) approximation assumes that a single, known material is being imaged. In this case, the values of 𝛿 and 𝛽 may be computed (typically at the mean energy of the x-ray beam).

As in the PAD model, attenuation is proportional to phase so the contact image can be written as

Equation 10

As in PAD, this leads to an equation for the propagated intensity of the form given in Eq. (8) with 𝐶 = (𝜆𝑅₂𝛿) /(4𝜋𝛽𝑀). The solution has the form of Eq. (9). The transfer function has the same shape as PAD, as illustrated in Fig. 3, and therefore the SM model does not typically require additional regularization. The drawback is an obvious one: if there is more than one material present, this approximation will not produce accurate phase.

Mesh-based imaging utilizes a metallic, wire mesh with a period on the order of 100 microns. The mesh is coarse enough to be directly imaged on an x-ray camera. A reference image is acquired with no object in place. With the object in place, the mesh lines will be distorted due to x-ray deflection. By measuring the distortion of the mesh, the deflection can be quantified and used to produce measurements of the phase gradient perpendicular to the mesh lines.

Since the mesh is periodic, by inserting the Fourier series expansion of the mesh into Eq. (4), one can obtain an expression that describes attenuation and phase information contained at each mesh harmonic. Phase reconstruction using this method is illustrated in Fig. 4.

Fig. 4 Diagram of the mesh-based imaging system

First, images are acquired and Fourier transformed both with and without the object in place. Each harmonic in the Fourier domains is windowed with a window of diameter 2𝑝𝜋/p, where p is the mesh period.

These windowed results are inverse Fourier transformed back to the space domain. Each harmonic image with the object in place is divided by its corresponding image without the object in place. The resulting image satisfies an expression of the form

Equation 11

central (0,0) harmonic contains only the attenuation image. The higher order harmonics contain phase derivatives in directions perpendicular to the mesh wires which can be isolated from attenuation by division with the (0,0) harmonic.

The preceding models all assume a pure point source of x rays, though this can be readily extended to plane-wave illumination (for example in a synchrotron) by moving the point source infinitely far from the object. Extended laboratory sources are partially coherent which results in a degradation of the phase reconstruction. This typically limits source sizes to roughly 50 microns or smaller, though the mesh-based imaging relaxes this requirement somewhat.

Instead of using small, low-power microfocus sources, we instead use focusing polycapillary optics with a large spot source to create a smaller focal spot with a more restricted divergence.

Polycapillary optics are arrays of hundreds of thousands of hollow glass tubes that guide the x rays by total internal reflection along their length. The physics and applications of focusing polycapillary optics are extensively reviewed. In the results presented here, an optic was used to produce a 100 um focal spot coupled with a 35W Rh source.

Fig. 5 Harmonic windowing for mesh-based imaging. (a) The four central harmonics. (b) Windowing 2 harmonics

Fig. 6 Diagram of a focusing polycapillary optic

Images were acquired using an Oxford Ultrabright Rh source operated at 42.5 kVp and 0.84 mA, which produced an x ray spectrum with a mean energy of 21 , near the Kα energy of 20 keV. Focusing optic XOS, Inc. serial number 3204 produced a roughly Gaussian spot of 100 µm FWHM at an output focal length of f=49 mm. In all cases, the sample used was a nylon (C₁₂H₂₂N₂O₂) rod of diameter 1.6 mm with density of 1.15 g/cm³. All images were acquired on a Radeye HR detector with 22 µm pixels. The source-to-object distance was 𝑅1 = 35 cm . For the propagation based imaging, the object-to-detector distance was 𝑅₂ = 55 cm. For the mesh-based imaging, , the object-to-mesh distance was was 𝑅_om = 14.8 cm and the mesh-to-camera distance was was 𝑅_mc = 20.8 cm. The mesh period was 198 µm, and the exposure was 300 mAs.

Quantitative phase at 21 keV was calculated from theory to be approximately 100 radians. Using the weak attenuation (WA) model, the reconstructed phase depth of the rod and even its spatial profile varied widely with regularization parameter as shown in Fig. 7. For low regularization, α =10⁻⁵, attenuation in the bulk of the rod was interpreted as phase and the reconstructed phase was too large by two orders of magnitude. For moderate regularization that confined the reconstruction to feature sizes roughly equal to the rod’s diameter, 𝛼 = 10^−3.5, the reconstructed phase is within a factor of 3 the true phase, but significant artifacts remain at the edge of the object due to attenuation being misinterpreted as phase. For strong regularization, 𝛼 = 10^−3.0 , spatial frequency and therefore feature size that is allowed in the reconstruction is so significantly limited that only the fringes at the edge of the rod are visible and the reconstructed profile is flat in the center of the rod.

Fig. 7 Phase reconstruction (radians) using weak attenuation for (from left to right) low regularization, 𝛼 = 10^−5, moderate regularization, 𝛼 = 10^−3.5, and strong regularization, 𝛼 = 10^−3.0. At 21 keV, the theoretical value of phase should be 100 radians. The upper row is the processed phase map and the bottom row is the profile at the marked location.

The reconstructed phase depth profiles produced from the pahse attenuation duality (PAD) model are shown in Fig. 8. As expected, because PAD incorrectly attributes attenuation mostly to Compton scattering, the phase estimates are off by a factor of between 2-3 from the theoretically computed value of 100 radians at 21 keV. Because the beam is polychromatic, the model was also testing assuming a beam energy of 15 keV, which gives a theoretical value of 139 radians, and a reconstruction which is also off by a similar factor.

Fig. 8 Phase reconstruction (in radians) using phase-attenuation duality with an assumed beam energy of (left) 15 keV and (right) 21 keV. The theoretical value of phase at 15 keV is 139 radians while the reconstruction is approximately 340. At 21 keV, the theoretical value is 100 radians with a reconstructed value of roughly 260 radians.

For the single material (SM) assumption, the reconstructed phase depth profiles are shown in Fig. 9. The reconstruction is fairly accurate at 150 measured radians (compared to 139 predicted) at 15 keV, but at the beam mean energy of 21 keV, the reconstruction is approximately 220 radians compared with 100 predicted.

Fig. 9 Phase reconstruction (in radians) using a single material assumption with an assumed beam energy of (left) 15 keV and (right) 21 keV. The theoretical value of phase at 15 keV is 139 radians while the reconstruction is approximately 170. At 21 keV, the theoretical value is 100 radians with a reconstructed value of roughly 260 radians.

The mesh-based imaging yields differential phase contrast (DPC) images which are directional derivatives of the phase.

By computing the vertical derivative of the phase using the mesh-based method and integrating it along a vertical line through the DPC image, as shown in Fig. 10, the quantitative phase of the rod can be estimated. Utilizing 60 different lines through the rod yields results of 107.6 radians ± 11.6 radians at 21 keV. We suspect this improvement is due to the fact that a 100 µm spot is larger than those typically used for propagation-based phase imaging which likely impacted the quality of our phase reconstructions. By utilizing a mesh-based method, the phase signature was substantially enhanced, enabling reasonable quantitative phase even with this larger focal spot.

Fig. 10 Vertical DPC image of the phase and integral. The integrated phase was compared along 60 lines through the rod and yielded a mean value of 107.6 radians ± 11.6 radians.

We have demonstrated the ability of focusing polycapillary optics to focus the x rays from large spot conventional sources to spots of sufficient size and quality for propagation and mesh-based phase contrast imaging. While the quantitative phase results for the WA, PAD and SM models using an optic of 100 µm were not accurate, likely because of the size of the focal spot, the use of a wire mesh enabled quantitative phase imaging with this same focal spot by boosting the phase signature. Future work will include incorporating source size in the reconstructions, which can be done with deconvolution or modifying the constant C in SM and PAD reconstructions. Since the reconstruction and theoretical phase vary significantly over the spectrum of the source, the source spectrum will also be incorporated into the models.