1. Introduction

Phase contrast imaging techniques can enhance the contrast of an image, over that produced by attenuation, by rendering the changes imparted by an object on the phase of a wavefield visible. This is particularly beneficial when the object displays weak attenuation contrast, as is often the case for biological specimens. Phase contrast imaging techniques have been developed for many imaging modalities including visible light microscopy, X-ray and neutron imaging, ultrasonography and magnetic resonance imaging. Conventional X-ray imaging remains the best approach for producing high spatial resolution images with the ability to penetrate through thick, optically opaque materials, and is thus widely employed for clinical and materials science applications. X-ray phase contrast modalities include crystal-based interferometry, grating interferometry, propagation-based imaging (PBI) and analyser-based imaging (ABI). For detailed reviews on these techniques see [1–3] and references therein.

ABI affords very high spatial and angular resolution imaging, provides high sensitivity to weak phase shifts and ultra-small angle X-ray scattering (USAXS) induced by the sample, and yields images essentially free from incoherently (i.e., with modification in the wavelength) scattered radiation. ABI requires a post-sample crystalline optic (the ‘analyser’) to produce phase contrast and requires highly temporally coherent radiation, but its requirement for high spatial coherence is relaxed (cf. PBI) [1–3]. In ABI, transverse phase shifts, with phase gradient components parallel to the analyser’s diffraction plane (Fig. 1 ), are rendered visible upon Bragg reflection from the analyser that specularly reflects the incident beam with the reflectivity depending on the angle of incidence of the X-ray beam. Thus, phase gradients in the wavefield just after the object are converted into intensity variations recorded by the detector placed after the analyser crystal. The angularly dependent reflectivity, called the ‘rocking curve’, can be measured by recording the reflected intensity upon rotating the crystal through the Bragg condition. The shape of experimental rocking curves depends on the choice of the crystal type and orientation, the crystal thickness, and the divergence and energy spectrum of the incident illumination. Rocking curve thickness oscillations can be measured when the incident wavefield is highly coherent [4]; these oscillations are readily smoothed by imperfections in the imaging system. In this setting the entire rocking curve may well be approximated by a bell-shaped curve [5–7]. This broadens the angular range over which image reconstruction can be accurately performed. Figure 1 shows the ABI setup used in this research and Fig. 2 shows examples of measured rocking curves.

 

Fig. 1 Schematic of the experimental analyser-based X-ray imaging setup used in this research (not to scale).

Download Full Size | PDF

 

Fig. 2 (a) Measured rocking curves (RCs) for the transmitted and diffracted beams through the Si(1 1 1) Laue analyser using 26 keV X-rays. This plot shows the intrinsic rocking curves measured with no sample in the beam and with the beam transmitted through the lung of a mouse. The rocking curves measured through the lung show reduced intensity due to absorption; a shift in peak position due to refraction; and broadening by USAXS within the sample. Each curve is measured using a single pixel from 260 angular measurements in 0.1 arc second angular steps. The vertical (dashed) orange line shows the angle at which images in Fig. 3 were recorded. (b) The diffracted (intrinsic) rocking curve from (a) is here symmetrised by mirroring the intensity on the positive side of the curve to the negative side. Both a symmetric Gaussian and Pearson type VII function are fit to the experimental data.

Download Full Size | PDF

Contrast in ABI is produced by the combined effects of absorption, refraction and scattering. If the sample contains features smaller than the detector resolution, they may scatter the beam such that a given pixel may register X-rays reflected by the analyser at numerous angles, with a reflectivity determined by the rocking curve [8]. X-rays scattered to angles within the angular acceptance of the rocking curve, from either USAXS or refraction, can result in a significant change in contrast. Scattering of photons can be evidenced by measuring the rocking curve with the sample placed in the beam, the width of which broadens as a result of scattering (see Fig. 2). This long-known effect [9] results from the convolution of the sample’s scattering distribution with the intrinsic rocking curve of the analyser [6–14].

Numerous research teams have sought to improve upon existing methods of extracting quantitative information from analyser-based X-ray phase contrast images. For non-homogeneous samples that do not produce USAXS, attenuation and phase information can be separated using two (or more) analyser-based images acquired at distinct orientations of the analyser. Many phase retrieval algorithms have been developed to address this problem using various approximations. Nesterets et al. [15] showed that a quantitative reconstruction can be achieved under the ‘weak-object’ approximation, such that phase and intensity variations across the sample are very small. However, this approximation will be violated by many macroscopic objects of interest in the biomedical or materials science communities. Paganin et al. [16] showed that quantitative phase retrieval is possible from two or four images if the coherent transfer function of the system can be linearized as a function of spatial frequency. However, measuring the complex transfer function of the system can be difficult, even for a monochromatic plane wave [17], which imposes a practical challenge on the application of their method.

Most ABI phase retrieval methods invoke the geometrical optics approximation (GOA), which is valid when the phase of the wave is slowly varying over the length scale of the extinction length of the crystal or N_T>>1, where N_T is the Takagi number [12, 18]. Therefore, the GOA may fail for instance at boundaries between materials of highly variable refractive index. However, aligning the crystal away from the steep slopes of the rocking curve can significantly extend the validity of the GOA across the entire image space [5].

An early ABI phase retrieval algorithm, based on the GOA, employed a linear approximation to the rocking curve slopes and utilized an image acquired on each slope [19]. Whilst this provides an elegantly simple reconstruction its accuracy is limited to very small changes in beam direction within the linear approximation. That work has been extended by numerous authors, with some of the methods quantitatively compared in articles by Huang et al. [20], Hu et al. [21] and Diemoz et al. [22]. The valid range of refraction angles has been extended by taking into account non-linear terms of the rocking curve [23], improving the fit to the rocking curve using polynomial curve fitting [24] and Gaussian functions [12, 20, 21, 24, 25]. A few articles compared these (and other) ABI phase retrieval methods and found Gaussian curve fitting to provide the best reconstruction [20, 22, 24]. However, Fig. 2(b) clearly shows that a Gaussian function can be a poor fit to experimental rocking curves, particularly at the peak and tails, a problem long known for fitting X-ray diffraction peaks [26]. This is due to the ‘long-slit geometry’ of ABI, whereby scattering is integrated in a direction perpendicular to the diffraction plane, broadening the tails of the rocking curve [7]. Significantly improved curve fitting, and hence phase retrieval, has been achieved using pseudo-Voigtian [7] and Pearson type VII functions [5] (see Fig. 2).

In the presence of a scattering object, phase information can be accurately recovered by accounting for the broadening of the rocking curve. This generally requires multiple images of the object at different analyser orientations to measure the rocking curve with the sample in the beam. The various reconstruction methods use statistical measures of the change of the rocking curve’s peak position, width and integrated reflectivity to respectively separate refraction, scattering and absorption information [6, 7, 10–14]. At least three images are required for each technique, but more accurate information is generally obtained with larger data sets [6]. Independent tests of competing algorithms reveal that the best reconstructions are obtained when rocking curves are measured with and without the sample at every pixel in the image and each is fit with a Gaussian [20, 22]. However, Suhonen et al. [7] showed that highly quantitative reconstructions can be obtained upon fitting the rocking curves with pseudo-Voigtian functions, that in principle should be more accurate than Gaussian curve fitting. Furthermore, it was since shown that Pearson VII functions provide better fitting to scattering curves measured using ABI than both Gaussian and pseudo-Voigtian functions [27].

Herein we demonstrate that: (i) Pearson VII functions can be fit to rocking curves to provide highly quantitative reconstructions of phase and amplitude information, and that a sample’s USAXS distribution can be quantitatively reconstructed via a deconvolution process (Sections 2.2 and 4.1); (ii) in the absence of USAXS, quantitative phase retrieval can be performed using two simultaneously acquired analyser-based images with the Laue geometry (Sections 2.3 and 4.2); and (iii) quantitative phase, absorption and scatter retrieval can be performed using two simultaneously acquired images via an iterative reconstruction technique (Sections 2.4 and 4.2). We finish with concluding remarks in Section 5.

2. Analyser-based phase retrieval under the geometrical optics approximation (GOA)

2.1 Rocking curve parameterisation using the Pearson type VII distribution

Experimentally measured rocking curves inevitably contain noise that can detrimentally affect image reconstruction. To overcome this we fit each rocking curve with a smooth function that is continuous over all angular space. Hall et al. [26] showed that the Pearson type VII distribution of the form

2.2 Multiple image radiography (MIR)

The process of separating the effects of refraction, absorption and USAXS using multiple images, obtained at different orientations, is often called multiple image radiography (MIR) [11]. Since fitting the experimental rocking curves with Gaussian functions can improve MIR reconstructions [20, 22], and given that Pearson VII functions can provide significantly better rocking curve fitting than Gaussians (Fig. 2) [5, 27], we have exploited this function to extract refraction, absorption and USAXS information from multiple ABI projections.

Using only the diffracted beam from a near-perfect Laue crystal we measured the rocking curve for every pixel in the image space with and without the object in the beam. Upon fitting with Pearson VII functions (Eq. (1)) we use the angular location of the rocking curve peak to define the angle of the beam incident on the crystal with $\left({\theta }_s\right)$ and without $({\theta }_0)$ the sample. The difference gives the change in beam direction (at each pixel) due to refraction:

The intensity ratio of images acquired at ${\theta }_s$ and ${\theta }_0$ provides an almost scatter-free attenuation map, here called the peak attenuation image. The image appears scatter-free since most scattered photons are rejected by the analyser. Taking the ratio of the integral of the rocking curve with the sample to that without is a measure of the intensity transmitted by the sample and provides a map of X-ray absorption. Since Eq. (1) does not lend itself to integration by analytic methods for arbitrary values of a, c and m, we integrate each curve numerically over a range of a few milliradians. We call the resultant reconstruction the average absorption image, as denoted in ref. [12]. Finally, the USAXS distribution can be determined by deconvolving the intrinsic rocking curve from the object rocking curve. Deconvolution is inherently unstable, particularly when using signals containing noise. We overcome this by using the noise-free Pearson VII approximations to the rocking curves and perform numerical deconvolution by Fourier space division using Wiener deconvolution. Wiener deconvolution uses a small regularizing parameter to prevent instability of the division-by-zero type. To map the object scattering function into a two-dimensional image we fit the resultant scattering curve $P\left({\theta }_{Scat}\right)$ with another Pearson VII function and calculate the half-width at half-maximum (HWHM) of the curve for each pixel as

2.3 Laue dual-image phase retrieval in the absence of USAXS

The requirement for multiple images in MIR may be problematic for medical applications where dose minimization is important. Moreover, the requirement to change the angle of the analyser between exposures in most ABI phase retrieval approaches (including MIR) prevents quantitative phase retrieval when imaging dynamic processes. Simultaneously recording the X-rays diffracted by and transmitted through a sufficiently thin analyser provides two images that yield differing phase contrast [28, 29]. Chapman et al. [29] provided a semi-quantitative demonstration that phase retrieval can be performed using these images alone, under the GOA, using a linear approximation to diffracted and transmitted rocking curves. Tomographic reconstructions using simulated images for the Laue geometry were also performed under the GOA by Bushuev and Guskova [4]. We have previously demonstrated experimentally that it is possible to perform quantitative phase retrieval using a Laue analyser [30]. In that work we employed nearest-neighbor linear interpolation to fit the intrinsic rocking curve for phase retrieval purposes. Herein we demonstrate a simplified GOA approach using Pearson VII fitting to diffracted and transmitted rocking curves for performing phase retrieval.

Under the GOA the intensity in the image plane is the product of the intensity transmitted through the object $(I_R)$ and the angularly dependent rocking curve [4, 29]:

Unfortunately, no simple analytic solution can be found for θ given the quotient $R\left(\theta \right)/T\left(\theta \right)$ when both $R\left(\theta \right)$ and $T\left(\theta \right)$ are expressed as Pearson VII functions. Instead, we can divide the measured diffracted rocking curve by the transmitted curve and fit the resultant curve with a Pearson VII function, i.e.:$\frac{I_D}{I_T}=\frac{I_{R}R\left(\theta \right)}{I_{R}T\left(\theta \right)}=c{\left[1+{\theta }^2/\left(m{a}^2\right)\right]}^{-m}.$ (5)

The centroid of this curve has been defined here as zero. Rearranging Eq. (5) allows one to solve for the angle of the beam incident on the crystal:

Since the ratio in Eq. (5) is a two-to-one function of theta, it provides a non-unique solution that cannot differentiate to which side of the Bragg peak the beam is aligned, as evidenced by the two branches of the square root in Eq. (6). Hence we must assume that the post-sample beam is aligned to the same side of the curve as the incident beam. Since the rocking curve is never truly symmetric, due to the Borrmann effect (see e.g. [31]), the curves can be reflected about their center to improve parameterization of the least squares regression (see Fig. 2).

In the absence of any sample, Eq. (6) yields ${\theta }_0$. With a sample, Eq. (6) yields ${\theta }_s$. The angular deflection by the sample can then be found using Eq. (2). To solve for the apparent absorption image we fit the transmitted rocking curve with an inverted Pearson VII function:

New parameters have been introduced to avoid confusion with the function used to fit the rocking curve ratio in Eq. (5). The apparent absorption is then found by rearranging Eq. (7):

2.4 Laue dual-image phase retrieval in the presence of USAXS: an iterative approach

Here we seek a phase retrieval algorithm capable of extracting the refraction, absorption and scattering information from just two simultaneously acquired images. We desire an analytic solution that accounts for convolution of the scattering distribution with the rocking curves. A simple analytic solution can be found if both functions are approximated using Gaussians. However, we have seen that Pearson VII functions can provide a much better fit to rocking curves. Analytical convolutions involving Pearson VII functions may not be possible if the decay constant m (see Eq. (1)) is neither an integer nor a half-integer, or may only be expressed as a power series that is not compact enough for easy use [26,32]. An analytic solution is therefore not practical using Pearson VII distributions. Convolutions can instead be performed numerically, with fast implementation possible via Fast Fourier Transforms.

To perform phase retrieval that accounts for three unknowns from just two images requires a priori knowledge about the sample. By choice, we require that the scattering distribution of the sample be known as a function of the projected thickness of the scattering volume. This requires that MIR be performed on a test sample of variable thickness, containing the same type of scattering medium found in the object of interest. For biological objects this could involve using excised tissue samples. The projected thickness of the scattering medium must either be known a priori or it may be measured from the MIR reconstructed absorption image or refraction angle map if either the linear attenuation coefficient (μ_l) or refractive index decrement $\left(\delta \right)$ is respectively known. If the scattering medium is held within some container the ${\mu }_l$ and δ values for the containers must also be known. For example, if the attenuation map is used then the projected thickness of the scattering medium may be extracted using the Beer-Lambert attenuation law. To calculate the projected thickness of a single material from the refraction angle map we first determine the phase shift of the wavefield relative to vacuum under the GOA as $\phi \left(y\right)=k\int \mathrm{\Delta }\theta \left(y\right)dy$ [3], where $k=2\pi /\lambda$ and λ is the X-ray wavelength. Here the integral is parallel to the y-z diffraction plane (Fig. 1). Under the projection approximation, the phase of the wavefield beyond a single homogeneous material (including voids) is directly proportional to the projected thickness as $\phi =-k\delta t$ [3]. Combining these equations, the projected thickness is

Next we exploit the knowledge that $P\left(\theta \right)$ must have a mean scattering angle of zero. Furthermore, under plane wave illumination, the net phase excursion across any object must be zero so that the net refraction angle through that object is also zero. To express this mathematically, define the boundaries of a projection in, say, the y-direction to begin at position A and terminate at B. Outside the object the phase of the plane wave remains unchanged, hence the phase excursion $\mathrm{\Delta }\phi ={\phi }_B-{\phi }_A=0.$ Note that by first subtracting ${\theta }_0$ from the reconstruction, any curvature of the incident illumination is removed. Under the GOA, the local beam direction is proportional to the phase gradient; i.e., $\mathrm{\Delta }\theta \left(y\right)=\left(1/k\right)\left[\partial \phi \left(y\right)/\partial y\right]$ [3]. The net refraction angle is found by integrating $\mathrm{\Delta }\theta \left(y\right)$ from points A toB:${\displaystyle {\int }_A^B\mathrm{\Delta }\theta (y)dy=\frac{1}{k}}{\displaystyle {\int }_A^B\frac{\partial \phi (y)}{\partial y}}dy=0,$ hence the net refraction angle across any object is zero. If the location of the scattering volume is known the remaining area can be masked out of the remaining calculations. The edge enhancement provided by phase contrast can greatly assist in the demarcation of such boundaries.

Armed with this a priori knowledge, we construct an iterative phase retrieval procedure that begins by ignoring the effects of USAXS and uses the phase retrieval methodology outlined in Section 2.3. Where there is appreciable USAXS the reconstructed refraction angle map will not be quantitative. Since the average refraction angle across the object should be zero, errors in the reconstruction can be reduced by subtracting the average refraction angle along each column of the data parallel to the diffraction plane. We consider this the 0th order reconstruction for $\mathrm{\Delta }\theta$ since no scattering information has been reconstructed. Integrating the refraction angle map parallel to the diffraction plane (e.g., from point A to the particular value of y) estimates the scattering volume’s projected thickness at a given y-coordinate, as per Eq. (9). Using the previously modeled scattering distribution $P(\theta ;t)$ (see e.g., Eq. (3)), an estimate of the scattering distribution can be constructed. For each pixel in the region-of-interest, $P(\theta ;t)$ can be numerically convolved with the intrinsic rocking curve to obtain an estimate of the sample rocking curve. The phase retrieval procedure defined in Section 2.3 can then be followed to find the new estimate of $\mathrm{\Delta }\theta (x,y)$. We denote this the 1st order reconstruction that provides an estimate for both $\mathrm{\Delta }\theta (x,y)$and $P(\theta ;t;x,y)$. This procedure, including subtracting the average refraction angle, can be repeated until the optimal solution, with respect to either the known object thickness or the measured intensity (see below), is found. In practice we have found that averaging the$P(\theta ;t)$parameters between the current and the previous iteration leads, on average, to better convergence toward the optimal solution.

This iterative procedure requires that any ambiguity regarding the sign of $\mathrm{\Delta }\theta$ be resolved. Therefore,${\theta }_0$ must be set sufficiently far from the Bragg peak that minimal X-rays are deflected over the Bragg peak. If the projected thickness of the scattering volume is known a priori then the optimal scattering parameters can be immediately used. Alternatively one can use the reconstructed data until the iterative process converges to the known thickness. With no a priori knowledge of we instead start with the 0th order reconstruction of$\mathrm{\Delta }\theta$and then reconstruct the diffracted (or transmitted) intensity using, for example, Eqs. (2) and (5). When the reconstructed intensity best matches the experimentally measured intensity the optimal solution is obtained. Figure 4 provides a detailed flow-chart for this iterative procedure. Having successfully tested our iterative phase retrieval algorithm on noise-free simulated data (unpublished), we herein report our findings on experimental data in Section 4.2.

 

Fig. 4 Flow-chart for the dual-image iterative phase retrieval algorithm. * represents an optional step.

Download Full Size | PDF

3. Materials and methods

3.1. Image acquisition

All experiments were performed using synchrotron radiation at SPring-8, Japan, with images acquired in Hutch 3 of beamline 20B2 in the Biomedical Imaging Centre (proposal 2009A1882). For the Laue analyser we employed a nominally 100 μm thick silicon wafer that was part of a monolithic silicon slab, connected at the base and sides, to minimize warping of atomic planes. The thickness was designed to provide a strong diffracted beam at the Bragg peak with minimal transmittance at 25 keV, which is a suitable energy for phase contrast imaging of small animals [33], for the Si $(111)$ Bragg reflection. In practice we found this condition was best met using 26 keV X-rays, hence we employed a monochromatic 26 keV beam for all experiments reported herein. We used a Si $(111)$ double-bounce monochromator in a non-dispersive setup with the analyser. The source-to-object distance was around 210 m, giving an almost parallel beam incident on the analyser. The $(1\overline{1}0)$ face of the analyser wafer had an area of 50 mm × 40 mm, which is sufficiently large to image the chest of a mouse in a single exposure. Figure 1 illustrates the experimental setup.

Both diffracted and transmitted beams from the Laue analyser were recorded simultaneously using a single X-ray sensitive area detector. We employed a 4000 × 2672 pixel Hamamatsu CCD camera (C9300-124) with a tapered fiber optic bonded between the CCD chip and the 20 µm thick gadolinium oxysulfide (Gd₂O₂S; P43) phosphor. The taper ratio was 1.8:1, converting the native pixel size of 9 µm to an effective pixel size of 16.2 µm.

Rocking curve measurements were made by recording images of the beam as the crystal was rotated through the Bragg peak in 260 steps of 0.1 arc seconds. The exposure time for each image was 300 ms, resulting in a net exposure time of 78 s. This many images are not typically required for MIR, but it ensured the long rocking curve tails were well sampled.

3.2. Image processing

The diffracted and transmitted images must first be separated and aligned. With respect to the transmitted image the diffracted image may be distorted by detector distortions and/or poor alignment or miscutting of the crystal causing beam magnification. Detector distortions should be carefully measured and corrected. In this case we have removed non-linear distortions imparted by the tapered fiber optic. Distortions due to the analyser can be removed by appropriate geometric transforms; namely translation, rotation and scaling. To that end we placed three steel ball-bearing fiducial markers around the perimeter of the field-of-view (FOV) to define three reference points in each image for distortion correction.

For all imaging techniques presented here the detector dark current was first subtracted from each image for normalization of all images. For the dual-image phase retrieval algorithm it was necessary to measure the intrinsic diffracted and transmitted rocking curves (and their ratio) for every pixel in the image due to the local fluctuations in pixel sensitivity, resulting primarily from the fiber optic taper (see Fig. 1 and the hexagonal background features in Fig. 3 ). This notably increased the reconstruction time. To symmetrize the rocking curves before applying the Pearson VII fitting, a Gaussian function, summed with a quadratic function, was fit to the diffracted rocking curve to locate the peak center about which to reflect the function. For MIR these peak centroids were used to measure the refraction angle map.

 

Fig. 3 (a) Diffracted and (b) transmitted images of a mouse thorax recorded using a Si(1 1 1) Laue analyser crystal with 26 keV X-rays. Field-of-view (FOV): 21×21 mm². The mouse is submerged in a tube of water, hence the attenuation is similar everywhere but appears darker through the more attenuating bones and brighter through the air filled (low density) lung tissue. The bright curved (horizontal) line through the lung reveals the ventral edge of the lung following the shape of the diaphragm. The centers of the black circles on the ventral edge indicate the position at which the rocking curves in Fig. 2 were recorded.

Download Full Size | PDF

For the iterative phase retrieval process it is only necessary to apply the algorithm to parts of an image where there is a known scattering volume. In this work we have chosen to manually select portions of the image. We note that when objects of different refractive index overlap in projection, integration of the refraction angle map will only be valid for estimating the thickness of the scatterer when the refraction signal from the scatterer is the dominant signal. If this is not the case, further a priori knowledge of the scatterer’s projected thickness is required. In principle, it may be possible to automatically locate the region containing the scattering material using just two images, as discussed below.

3.3. Sample preparation

Two objects were employed to test the phase retrieval methods outlined above. The first consisted of a 40 × 40 × 6 mm³ polymethyl methacrylate (PMMA; C₅H₈0₂) block. Two polished cylindrical cavities of 4 mm diameter and 30 mm length were drilled inside the block. The cavities were positioned in the beam with their lengths aligned horizontally and perpendicular to the beam. In this orientation the cavities produced strong phase gradients within the diffraction plane. The upper cavity was left void whilst the lower cavity was filled with hollow glass microspheres (CEL-STAR^®, Z-27; Tokai Kogyo, Japan) with diameters ranging from 10 to 120 µm (63 µm average) and inner diameters of 8 to 12 μm. The small size and large refractive index gradients of the microspheres were expected to cause significant USAXS.

The second object was an adolescent Balb/c nude male mouse killed prior to imaging by an overdose of nembutal in accordance with regulations set by the Monash University animal ethics committee and the SPring-8 Animal Care and Use Committee. The mouse was placed in a cylinder filled with water to remove strong phase gradients at the skin/air boundaries. An endotracheal tube was inserted via a tracheotomy into the mid-cervical trachea and connected to a custom designed small animal ventilator [34]. The lungs were inflated to a fixed pressure of 25 cmH₂O that was kept constant throughout the duration of the imaging sequences.

4. Results and discussion

4.1. Multiple image radiography (MIR) results

Figure 5 shows the scattering distribution $P(\theta )$ of a section of mouse lung tissue measured by deconvolving the Pearson VII functions fitted to the diffracted rocking curves in Fig. 2. A regularization parameter of 1 × 10⁻¹³ was found to be the smallest value to give consistent results for the Wiener deconvolution for all samples. Weak oscillations seen at the tails of the curve in Fig. 5 reveal acceptably small artifacts from the deconvolution. The integrated probability of this curve is very close to unity (>0.99), which validates the use of the deconvolution process as it accounts for virtually all X-rays scattered at ultra-small angles. Overlaid with this plot is a Pearson VII fit to the scattering curve. The decay parameter $m<3$ indicates that the tails of the scattering curve are very broad, as expected in the ‘long-slit geometry’ [7].

 

Fig. 5 Angular X-ray scattering distribution$P(\theta )$ (thick black curve) of a section of mouse lung tissue highlighted in Fig. 3. The thin red curve shows the Pearson type VII function fit to the scattering distribution.

Download Full Size | PDF

Figure 6 shows the reconstructed MIR images of the mouse thorax using the procedure outlined in Section 2.2. Figure 6(a) shows the refraction angle map, from which we see that the strongest phase gradients in the mouse thorax arise at the outer boundaries of the lung tissue, particularly adjacent to the diaphragm (dark band) where the projected thickness of lung tissue rapidly changes. We also see similar phase gradients arising from bubbles located in the upper abdomen, most likely indicating gas bubbles in the stomach (lower right part of image indicated by hollow black arrow). Weaker phase gradients are visible at the bone/tissue interfaces. We note that the refraction at the air/lung interfaces is so much stronger than at the bone/tissue interfaces that the bones are barely visible where the lungs and bone structures overlap in projection. This justifies the use of the single-material assumption for the iterative phase retrieval of the lung in Section 4.2. Larger airways (bronchi; black arrow) can be seen branching into the lung lobes, but the trachea (white arrow) is less apparent since the strong phase gradients produced by the vertical trachea are perpendicular to the diffraction plane.

 

Fig. 6 Multiple image radiography of an adolescent mouse thorax. (a) Refraction angle image. (b) Scatter HWHM image. (c) Integrated absorption image. (d) Peak absorption image. FOV: 21 × 21 mm². X-ray energy: 26 keV. Net exposure time: 78 s.

Download Full Size | PDF

The scatter distribution map of the mouse thorax is seen in Fig. 6(b). This image, which is shown in color to adequately display the large dynamic range, reveals strong USAXS produced by the aerated lung tissue. Weaker scattering signals arise from the bones and from the bubbles in the abdomen. Figure 6(c) is the integrated absorption map of the thorax that clearly reveals the high density of the bones (black) against the soft tissue and water (grey). The low density aerated lung tissue is visible as a lighter shade of grey. With no phase contrast it is more difficult to observe the boundaries of the lung tissue. Figure 6(d) is the peak absorption image, which rejects most of the scattered radiation and thus causes the highly scattering lung tissue (cf. Figure 6(b)) to appear dark against surrounding soft tissues.

Figure 7 shows the reconstructed MIR images for the PMMA block with horizontally aligned cylindrical cavities. The hollow upper cavity provides strong phase gradients parallel to the diffraction plane (Fig. 7(a)). Interestingly, the edges of the hollow cavity provide a clear scattering signal (Fig. 7(b)), which results from large phase gradients producing a measureable dispersion of the X-ray beam within the width of a pixel. This phenomenon is related to Young’s boundary wave and the related geometrical theory of diffraction [35]. The same cavity provides only a small amount of attenuation contrast, as revealed in Fig. 7(c) and (d).

 

Fig. 7 Multiple image radiography of two cylindrical cavities in a PMMA block. The lower cavity is filled with hollow glass microspheres. (a) Refraction angle image. (b) Scatter HWHM image. (c) Integrated absorption image. (d) Peak absorption image. FOV: 25 × 16 mm². X-ray energy: 26 keV. Net exposure time: 78 s.

Download Full Size | PDF

The presence of the hollow glass microspheres in the lower cavity is revealed primarily in the scattering map (Fig. 7(b)) and the peak absorption image (Fig. 7(d)), which indicate the strong scattering effect of overlapping microspheres; akin to the microstructure of the minor airways of the mouse (Fig. 6). The microspheres are only evident in the refraction angle map (Fig. 7(a)) as fine speckles and, in the middle of the cavity, as a more obvious structure where the spheres appear to be inhomogeneously packed (hollow black arrow). Further evidence of inhomogeneous packing can be seen in the scattering map (Fig. 7(b)), particularly on the left hand edge where the cap appears to have compressed the spheres. In Fig. 7(a), (c) and (d) a faint horizontal band is also evident (black arrows in Fig. 7(a)). This resulted from a spurious reflection from the analyser altering the shape of the rocking curve at that position, thereby affecting the curve fitting procedure. Numerous horizontal bands are also visible in these images that arise from intensity variations within the beam caused by carbon build-up on the monochromator crystals. These streaks are normally removed by accounting for the intensity of the incident beam in the reconstruction; note their absence in Fig. 6(c) and (d). In this case the beam position clearly drifted significantly between acquiring images of the sample and the direct beam. Such artifacts can be avoided by regular cleaning of the monochromator, improving the stability of the optics, or by reducing data acquisition times.

4.2 Laue dual-image phase retrieval results

In this section we demonstrate dual-image phase retrieval for the same objects shown in section 4.1. We begin by calculating the refraction angle map using Eqs. (2) and (6) and the apparent absorption image using Eq. (8) by neglecting any USAXS produced by the sample. We then employ the iterative method to account for USAXS in areas of strong scattering.

Figure 8 shows the refraction angle maps of the same PMMA block shown in Fig. 7. In this instance the reflectivity of the diffracted beam was set to about 10% of the peak Bragg reflectivity. Figure 8(a) reveals the refraction angle map calculated neglecting USAXS. The result is expected to look similar to that provided in Fig. 7(a). Neglecting the effects of scatter clearly results in a poor reconstruction of the lower, microsphere-filled cavity. We note that the dark appearance of this cavity could be exploited for automated detection of high scattering ROIs using morphological operations such as thresholding or region-growing, etc., which remains the subject of future research.

 

Fig. 8 Refraction angle $(\mathrm{\Delta }\theta )$ images of the PMMA phantom calculated using two simultaneously recorded images. Here the diffracted intensity was approximately 10% of the Bragg peak reflectivity. FOV: 25 × 16 mm². X-ray energy: 26 keV. Net exposure time: 300 ms.

Download Full Size | PDF

To account for scattering in the iterative reconstruction procedure, we first measure the average parameters of the scattering distribution as a function of the projected scatterer thickness from the MIR data (Fig. 7). In this instance the projected thickness can be readily calculated using the mathematical projection of a cylinder, perpendicular to its axis. For each value of equal cavity thickness we take the average value for each scattering parameter (see Fig. 7(b)). Ideally a single polynomial function could be used to describe each such parameter as a function of thickness. However, we see in Fig. 7(b) that the scattering distribution has a sharp change at the edge of the cavity that is also seen in the cavity containing no microspheres, as discussed above. We therefore employed an ‘edge’ function to describe scattering within the first seven pixels (~113 µm) of the cavity boundary together with polynomial functions to map the average scattering coefficients within the remainder of the cavity. Figure 9(a) shows the polynomial fit to the average width parameterof the scattering distribution as a function of the projected cavity thickness away from its edges.

 

Fig. 9 Scattering width parameter a (see Eq. (1)) used for iterative phase retrieval of (a) the PMMA cavity containing hollow glass microspheres and (b) mouse lung tissue.

Download Full Size | PDF

Setting the mean refraction angle within the lower cavity to zero for each column yields Fig. 8(b), which qualitatively and quantitatively more closely resembles the upper cavity. We label this first approximation the 0th order iterative reconstruction. We then numerically integrate the refraction angle map to estimate the projected thickness of the cavity containing the scatterers (Eq. (9)) to estimate the scattering distribution within the cavity. The rocking curve parameters were then iteratively updated and the refraction angle map recalculated to find the map that yielded the projected thickness closest to the true projected thickness.

To achieve a highly quantitative reconstruction we allowed each column within the scattering region-of-interest (ROI) to be optimized independently. Figure 10 shows the root-mean-square (RMS) difference between the known and reconstructed projected cavity thicknesses for two different columns across the microsphere filled cavity. The iterative procedure was limited to 20 iterations since each iteration over the whole cavity took almost an hour on a personal computer with a 2.8 GHz dual core processor. This figure reveals that iteration does not monotonically improve the reconstruction, but instead displays oscillatory behavior. However, each shows a notable improvement in the reconstruction after a few iterations. We use the global minimum of this RMS difference function as the optimal reconstruction for each column. Using a different optimization for each column can, however, result in streak artifacts in the resultant image. We therefore selectively smooth the final image in Fourier space to remove only the high frequency vertical artifacts. The final, smoothed iterative reconstruction is shown in Fig. 8(c).

 

Fig. 10 RMS differences between iteratively reconstructed and known projected glass microsphere thickness for the reconstruction shown in Fig. 8(c). Black dashed lines there show the corresponding positions L1 and L2.

Download Full Size | PDF

Figure 11 displays the results of the refraction angle map of the same object shown in Fig. 8 using a different pair of images. In this case the reflectivity of the diffracted beam was around 1% of the maximum reflectivity. These images are noticeably noisier than the reconstructions of Fig. 8 due to the weaker signal in the diffracted image. Figure 11(a) shows a poor reconstruction of the microsphere-filled cavity when USAXS is ignored. The 0th order reconstruction (Fig. 11(b)) is this time qualitatively incorrect. A bright band, showing positive refraction angles, is evident at the bottom of the cavity where the refraction angle should be negative. The expected black band is only restored upon iteratively altering the shape of the rocking curve parameters used for the reconstruction, as evidenced by Fig. 11(c). Finally, Fig. 12 shows the reconstructed apparent absorption images and scattering distributions that coincide with the reconstructions shown in Fig. 8 and Fig. 11.

 

Fig. 11 Refraction angle $(\mathrm{\Delta }\theta )$ images of the PMMA phantom calculated using two simultaneously recorded images. Here the diffracted intensity was approximately 1% of the Bragg peak reflectivity. FOV: 25 × 16 mm². X-ray energy: 26 keV. Net exposure time: 300 ms.

Download Full Size | PDF

 

Fig. 12 Apparent absorption images of the PMMA phantom reconstructed at (a) 10% and (b) 1% Bragg reflectivity, respectively. (c) and (d) are scatter HWHM images of the microsphere-filled cavity at the same respective reflectivities as (a) and (b). X-ray energy: 26 keV. Net exposure time: 300 ms.

Download Full Size | PDF

A quantitative comparison of the reconstruction approaches is provided in Fig. 13 . Here we numerically integrate the reconstructions of the refraction angle of the cavities to show the projected thickness of the cavities via Eq. (9). Figure 13(a) shows the reconstructions for a single column across the hollow cavity compared to the known (modeled) projected thickness of the cavity. The MIR reconstruction provides a very close reconstruction of the projected cavity thickness. Using just two simultaneous acquired images with the diffracted intensity at 10% of the peak reflectivity provides a noisier but highly quantitative reconstruction of the cavity. A noisier but still quantitative reconstruction is obtained using two images with the diffracted intensity at 1% of the peak reflectivity. Figure 13(b) shows the same sequence of reconstructions for the cavity filled with glass microspheres. Not surprisingly, each reconstruction is slightly less quantitative than the reconstruction of the hollow cavity due to the added complication of the highly scattering microspheres.

 

Fig. 13 Projected thickness profiles of the (a) hollow cavity and (b) microsphere filled cavity shown reconstructed using MIR (Fig. 7), and dual-image ABI at 10% reflectivity (Fig. 8) and 1% reflectivity (Fig. 11) (see white dashed lines in figures for profile locations).

Download Full Size | PDF

A quantitative comparison of the different approaches, averaged across the full length of the cavities, is provided in Table 1 . MIR always provides a quantitative reconstruction of the cavity with the least variability, which is expected due to the large number of images ensuring a high quality reconstruction. However, MIR consistently underestimates the thickness of the cavity as a result of the refraction at the edge of the sample being underestimated – a typical result for GOA based reconstructions [20, 22]. Our dual-image iterative reconstructions can result in the projected thickness being under- or over-estimated, as seen in Fig. 13 and Table 1. We note from Table 1 that the significantly reduced number of images (2 vs. 260) increases the amount of variability in the data appreciably. However, the variability is worst for the images obtained with the diffracted image recorded at the 1% reflectivity level. Although the intensity in the transmitted image $(I_T)$ increases with the decrease in the diffracted intensity $(I_D)$, the noise in reconstructed images increases as$I_D$tends to zero, due to the division by$I_D$in Eq. (6). Thus, for dual-image reconstructions, although the crystal should be detuned sufficiently far from the Bragg peak to avoid ambiguity of the sign of the refraction angle (Eq. (6)), it is important to have sufficient intensity in both images to suppress noise in the final reconstructions.

Table 1. Volume fraction (%) and FOM (in parentheses) of the cylindrical cavities measured using MIR and dual-image reconstructions.

View Table

Table 1 also provides a figure-of-merit (FOM) for the three reconstructions. The FOM is the (contrast-to-noise ratio)²/(entrance skin dose) [36], where the skin dose was measured to be approximately 23 mGy/s. The contrast was measured from average absorption through the cavities compared with that through the PMMA block. We note that the FOM is very low for the MIR data due to the long exposure time and noise introduced by the motion of the beam. This could be significantly improved by recording fewer images over the rocking curve, so long as sufficient points (at least three) are used to allow adequate fitting of the curves (ideally six or more [6]). Even so, the dual-image method will always yield a high FOM due to the comparatively low dose.

Figure 14 shows the reconstructed absorption, scattering and refraction angle maps of the mouse lung images shown in Fig. 6, here calculated using the dual-image iterative algorithm. Average scattering parameters were first measured as a function of the projected lung thickness from the MIR data. We use Eq. (9) to estimate the lung thickness by approximating the chest as a single material, neglecting the small refraction effect from the bone/tissue interfaces. Figure 9(b) shows the polynomial fitting applied to the scattering width parameter. Since no averaging could be performed the scattering parameters were highly variable. This variability, together with the comparatively smooth boundary between aerated lung and surrounding soft tissue, meant that no ‘edge function’ was required for this parameterization.

 

Fig. 14 Refraction angle images of the mouse thorax calculated using two simultaneously recorded images. Here the diffracted intensity was approximately 15% of the Bragg peak intensity. (a) Initial refraction angle reconstruction. (b) Best iterative refraction angle reconstruction. (c) Scatter HWHM image. (d) Apparent absorption image. FOV: 21 × 21 mm². X-ray energy: 26 keV. Net exposure time: 300 ms.

Download Full Size | PDF

Figure 14(a) shows the initial refraction angle map calculated before accounting for USAXS. Figure 14(b) shows the optimized refraction angle map after 20 iterations applied to a manually selected ROI containing the lung tissue. The optimal iteration number used for each column across the highly scattering lung tissue is plotted in Fig. 15(a) . Qualitatively, the refraction angle map of Fig. 14(b) is in good agreement with the MIR result of Fig. 6(a). We note that whilst the refraction angle measured within the lung is in good agreement with the expected values, bones outside the lung appear slightly dark compared with those in Fig. 6(a) as we have not parameterized the scattering by the bones for this reconstruction. Figure 14(c) and (d) reveal the optimized scattering distribution map and the apparent absorption image of the mouse thorax. Both are in excellent agreement with the MIR results of Fig. 6(b) and (c).

 

Fig. 15 (a) Optimal number of iterations used to generate the refraction angle map for the mouse lung shown in Fig. 14(b). (b) Projected lung volume calculated from Fig. 6(b) (MIR – thick black line) and Fig. 14(b) (dual-image iterative reconstruction – thin red line) summed along vertical columns.

Download Full Size | PDF

A quantitative comparison between the refraction angle map of lung tissue measured using MIR (Fig. 6(a)) and dual-image phase retrieval (Fig. 14(b)) is provided in Fig. 15(b). The refraction angle has been integrated parallel to the (vertical) diffraction axis to estimate the projected lung tissue thickness (Eq. (9)). Figure 15(b) shows the sum of the projected thickness along each column to reveal the lung air volume measured within each column. The dual-image phase retrieval yields an average airway volume 97.5% ± 17.8% of the volume measured using the 260-image MIR reconstruction. Furthermore, from Fig. 15 we note that the number of iterations required for an optimal solution is approximately proportional to the projected thickness of the scattering body.

5. Conclusions

We have developed two phase retrieval techniques for analyser-based phase contrast X-ray imaging that provide information about an object’s X-ray absorption, refraction and scattering properties. The first measures the complete rocking curves of the post-analyser diffracted beam with and without the sample. Fitting these curves with Pearson type VII functions enables highly quantitative reconstructions to be obtained. The second method requires only two images to be acquired and uses an iterative approach to reconstruct the object properties based on a priori knowledge of how scattering by the object varies as a function of the scattering object’s projected thickness. Whilst the reconstructions are expectedly noisier than the multiple image approach, the dose to the sample is substantially reduced. Moreover, simultaneously detecting the X-rays transmitted through a thin analyser crystal in the Laue geometry allows quantitative phase retrieval to be performed without rotating the analyser. This will enable quantitative reconstructions to be performed on images of dynamically recorded processes such as moving biological organisms or mechanical systems.