Open Access
Add to BibSonomy
Abstract
X-ray diffusive dark-field imaging, which allows spatially unresolved microstructure to be mapped across a sample, is an increasingly popular tool in an array of settings. Here, we present a new algorithm for phase and dark-field computed tomography based on the x-ray Fokker-Planck equation. Needing only a coherent x-ray source, sample, and detector, our propagation-based algorithm can map the sample density and dark-field/diffusion properties of the sample in 3D. Importantly, incorporating dark-field information in the density reconstruction process enables a higher spatial resolution reconstruction than possible with previous propagation-based approaches. Two sample exposures at each projection angle are sufficient for the successful reconstruction of both the sample density and dark-field Fokker-Planck diffusion coefficients. We anticipate that the proposed algorithm may be of benefit in biomedical imaging and industrial settings.
© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
X-ray imaging, both in projection and in computed tomography (CT) [1], is an indispensable tool for visualizing samples in a non-invasive manner. The technique has been adopted across a vast array of fields, such as medicine, security, and manufacturing. Phase-contrast x-ray imaging [2], and its extension to computed tomography, has been a major focus for several decades within the x-ray imaging research community. In recent times, the emergence of x-ray diffusive dark-field imaging and x-ray dark-field computed tomography [3–8] has attracted considerable interest. Diffusive dark-field imaging (‘dark-field imaging’ henceforth) measures a signal that arises from diffuse scattering of the illuminating x-ray wavefield from sub-pixel features within the sample, and hence allows the influence of such sub-pixel features to be detected. This approach is already finding application in clinical settings [9–11]. This dark-field signal has primarily been accessed by using imaging methods that place additional optical elements in the experimental configuration, such as analyzer-based imaging [12], grating interferometry [3,13,14], speckle-based/structured illumination imaging [15–21], beam tracking [22] and edge illumination [23,24]. Recent work has shown that the dark-field signal may also be captured using propagation-based x-ray imaging [25–29].
In propagation-based imaging [30–32], x-rays of sufficient spatial coherence are transmitted through a sample and undergo free-space propagation through a distance $\Delta$ from the sample to the detector plane. This free-space propagation converts sample-induced phase variations in the x-ray beam into intensity variations at the detector plane, which can then be measured. Additionally, unresolved microstructure contained within the volume of the sample will cause a fraction of the x-rays to be diffusely scattered, resulting in a local visibility reduction in the intensity pattern recorded at the detector, which is interpreted as dark-field signal [28]. We assume this unresolved microstructure to be spatially random, meaning that the unresolved microstructure is randomly distributed through the volume of the sample. We exclude the effects of partially ordered or fully ordered microstructure (e.g. crystals), which are beyond the scope of the present paper. Traditionally, such diffuse scattering has either not been considered or has been assumed to be negligible, in which case the image formation process in propagation-based phase-contrast imaging can be modeled by the transport-of-intensity equation (TIE) [33]. Several software packages for retrieving phase information from propagation-based x-ray images based on the TIE have been developed [34–36], which along with the adaption of TIE-based phase-retrieval algorithms to multi-material objects [37] and the use of partially spatially coherent radiation [38], have made propagation-based imaging a widely accessible technique.
In this paper, we extend the dark-field retrieval method presented in our previous work [28] to CT. This approach is based on the following model for image contrast seen when a coherent x-ray wavefield propagates, namely the Fokker-Planck equation for paraxial wave optics [25,39]:
Fig. 1. Experimental schematic to perform computed tomography using propagation-based imaging. The blurry region of the snippet of the image shown at the detector plane arises as a result of the sample-induced diffusion due to SAXS. This sample-induced diffusion can be extracted, along with the sample thickness, by using two intensity measurements at different propagation distances. By collecting intensity measurements at different sample rotation angles, a tomographic reconstruction of the sample density and dark-field signal can be performed.
The x-ray Fokker-Planck equation describes how the optical energy carried by the illuminating wavefield is transversely redistributed as it propagates downstream of the sample. A fraction of the optical energy carried by the illuminating wavefield is converted to SAXS, while the remaining fraction is coherently transported, resulting in a bifurcation into coherent and diffusive energy-flow channels. The coherent energy-flow channel is associated with phase effects (propagation-based phase-contrast), while the diffusive energy-flow channel is associated with x-ray dark-field effects (position-dependent visibility reduction). Returning to a point made earlier, phase effects are linked to the refraction of the emergent x-rays due to the sample, which results in bright/dark fringes in the recorded intensity pattern, which increase in width and visibility with propagation [30,31]. Dark-field effects arise due to the diffuse scattering of the emergent x-rays by the unresolved spatially-random microstructure contained in the sample, with this diffuse scattering blurring the local contrast [41]. These phase and dark-field effects are seen differently and evolve differently with propagation distance, as evidenced by the different powers of $\Delta$ that appear in the final two terms of Eq. (1). This observation is central to the reconstruction method presented in this paper. Such phase and dark-field effects can be seen in Fig. 1 below, which shows the propagation-based imaging setup that corresponds to Eq. (1). It should be noted that we use the terms ‘sample density’ and ‘linear attenuation coefficient’ interchangeably, where ‘sample density’ is taken to mean ‘sample number density’. This equivalence is justified via the relationship (cf. Eq. (6) in Ref. [42]):
where $\mu (x,y,z)$ is the linear attenuation coefficient, $\sigma$ is the total x-ray cross section for one scatterer and $\rho (x,y,z)$ is the number density of scatterers in the sample.As additional context for this work, we note that the concept of spatial resolution and how it can be improved is a key factor in determining the usefulness of a CT reconstruction technique. In order to improve spatial resolution, CT reconstructions that utilize the TIE phase retrieval method of Paganin et al. [34] have sometimes incorporated deconvolution filters [35] or taken source-size blurring into account [38], and recently have shown an improvement in spatial resolution by incorporating the discrete mathematics used in processing digital images in the phase retrieval algorithms [43,44]. The present paper extends the method presented in Leatham et al. [28] to the realm of CT, and also aims to show that 1) a dark-field CT can be captured without additional optics, using only propagation-based imaging, through the introduction of a tomographic Fokker-Planck linear diffusion coefficient which quantifies the diffusion of the x-ray beam by a given voxel of the sample, and 2) by separating out phase and dark-field effects, an improvement in spatial resolution can be seen in the recovered sample density image when compared to the TIE-based phase retrieval method of Paganin et al. [34]. Such an improvement in spatial resolution may be anticipated for the following reason. Since the projected intensity pattern is locally blurred by dark-field effects, the spatial resolution of both the sample thickness projections and CT sample density slices is reduced. If dark-field effects are not taken into account, such as in the TIE phase retrieval method [34], the spatial resolution of the recovered images suffers.
The paper is structured as follows. Section 2 details our Fokker-Planck reconstruction method for x-ray phase and dark-field CT without optical elements, and section 3 shows the results of the experiment we performed to demonstrate the ability of our method to tomographically reconstruct phase and dark-field signals in three dimensions. Section 4 discusses the broader implications of our work and section 5 outlines avenues for future research, as well as providing some concluding remarks.
2. Method
The CT reconstruction method we propose here builds directly on the phase and dark-field retrieval method presented in Leatham et al. [28]. In this approach, outlined in Fig. 2, a single-material approximation is used to write the sample-exit-surface intensity and phase terms in the Fokker-Planck equation as a function of sample thickness only. Two unknowns, the sample projected thickness $T$ and sample diffusion coefficient $D$, are then determined by solving this single-material Fokker-Planck equation using data obtained at two different propagation distances $\Delta _{1}$ and $\Delta _{2}$. Given such propagation-based intensity images captured at two different sample-to-detector propagation distances ($\Delta _{1}\neq \Delta _{2}$), at a fixed energy, over a set of sample projection angles, denoted by $\Theta$, and assuming quasi-monochromatic $z$-directed plane waves of incident intensity $I_{0}$ illuminate the sample, our method for dark-field CT constitutes the following. Firstly, using the two propagation-based intensity images captured at each projection angle, we calculate the sample transmission at each angle, $t(x,y;\Theta )$, according to (cf. Ref. [28]):
The final step of our algorithm is to perform CT reconstruction from the calculated sample transmission and dark-field signals. In particular, we input the calculated sample transmissions into the CT reconstruction code of Ref. [36] to obtain the three-dimensional distribution of the linear attenuation coefficient, $\mu (x,y,z)$, for instance by using filtered-back-projection [46]. The parameter $\gamma$ determines the sharpness of boundaries between adjacent materials in the reconstructed volume. In order to perform CT reconstruction from the projected dark-field signals, we define the quantity $\mu _{D}$, which we refer to as the ‘Fokker-Planck linear diffusion coefficient’, via
The quantity $\mu _{D}$ quantifies how much a given voxel of the sample diffuses the x-ray beam, independent of the set-up used to capture this measurement. From the projected dark-field signals given by Eq. (5), we create projections of $\exp [-D(x,y;\Theta )]$, which we input into the CT reconstruction code of Ref. [36] in order to obtain the three-dimensional distribution of the Fokker-Planck linear diffusion coefficient, $\mu _{D}$, on account of Eq. (8). Note that the CT code used in this paper, and indeed most CT codes, assume that the input is ‘transmitted intensity’, and hence take the negative logarithm of the inputted projections. Hence, having exponentiated images, such as sample transmission given by Eq. (4) and projections of $\exp [-D(x,y;\Theta )]$, is crucial to obtaining the desired distribution of $\mu (x,y,z)$ and $\mu _{D}(x,y,z)$ from such CT codes. A summary of our algorithm workflow can be found in Fig. 2.Fig. 2. Flowchart summarizing the propagation-based method described in this section for dark-field and phase CT.
In the case of a single-material sample with position-dependent projected thickness $T$, whose statistical properties are independent of depth $z$ for any particular transverse location $(x,y)$, Eq. (8) reduces to
where functional dependencies have been omitted. The result given by Eq. (9) could be anticipated based on the following logic. The root-mean-square (RMS) apex half-angle, $\theta _{s}$ is proportional to the square root of the sample thickness [40]. Given the relationship in Eq. (2), this implies that the dark-field signal is proportional to the sample thickness. Taking the Fokker-Planck linear diffusion coefficient, $\mu _{D}$, to be the factor of proportionality, one arrives at Eq. (9).To close this section, we link the Fokker-Planck linear diffusion coefficient to the dark-field extinction coefficient, as employed in grating interferometry. Imaging methods that make use of reference patterns, for example, grating-based imaging or single-grid imaging, often speak in terms of a visibility reduction that arises due to unresolved sample microstructure contained in the sample [5,47]. This visibility reduction, $V$, can be expressed in terms of the dark-field signal via (see Eq. (132) of Ref. [40], also mentioned in Ref. [39])
where $p$ is the period of the analyzer grating. This visibility reduction is used to tomographically quantify grating-interferometry-based dark-field via the broadly-adopted sample linear diffusion coefficient / dark-field extinction coefficient, $\mu _{d}$, for a single-material sample [5,47]: Equating Eq. (10) and Eq. (11), and solving for $\mu _d$, one obtains: where Eq. (9) has been used in the second step. As the propagation distance, $\Delta$, and the period of the reference pattern, $p$, are properties of the experimental setup, Eq. (12) provides a way to convert the Fokker-Planck linear diffusion coefficient $\mu _{D}$, which is independent of experimental setup parameters, to the more broadly-adopted set-up-dependent linear diffusion coefficient, $\mu _{d}$.3. Experimental results
In order to demonstrate the applicability of our CT reconstruction method outlined in the previous section, we collected propagation-based x-ray phase-contrast images at the Australian Synchrotron, on the Imaging and Medical Beamline (IMBL) in Hutch 3B. The sample consisted of a 1.3 cm diameter polymethyl methacrylate (PMMA) tube with nylon wire wrapped around the outside of the tube, creating some background texture that was part of the sample, as well as three PMMA microtubes sitting inside the larger tube. With reference to the $0^\circ$ projections (see Fig. 3), the bottom-left microtube contained water, the bottom-right microtube contained PMMA microspheres and the upper microtube contained agarose powder. The PMMA microspheres were supplied by the company Polysciences Inc. (catalogue number 26305-500). The diameters of the microspheres were distributed within the range 1-10 µm, with a mean diameter of 6 µm, and the density of the microspheres was approximately ${1.19}\;\textrm{g/cm}^{3}$. The agarose powder contained particles spanning a relatively large range of diameters, with these particles being heterogeneous in size and shape. The diameters of the particles contained in the agarose powder were distributed between 20-160 µm, with an average diameter of approximately 80 µm. The sample was placed on a dedicated table located approximately 130 m from the source, where x-ray photons of energy 25 keV were produced from a 2 T dipole bending magnet. At this energy, the wavenumber is $k = 1.27{\times}10^{-11}\;\textrm{m}^{-1}$ and the complex refractive index of PMMA has $\delta = 4.26{\times}10^{-7}$ and $\beta = 1.81{\times}10^{-10}$, so that $\gamma =2353.96$ and $\mu_{s} = {45.9}\;\textrm{m}^{-1}$. This value of $\gamma$ was used to control the sharpness of the boundaries that appear in our sample when calculating the sample transmission, and subsequently the sample projected thickness, as well as the projected dark-field signals.
Fig. 3. Experimental images of the sample, which includes three PMMA tubes, containing water, agarose powder and microspheres, from left to right, shown here in projection, taken at $\Theta =0^\circ$. All images are shown on a linear grayscale, where $a$ is the minimum value and $b$ is the maximum value. Top row: propagation-based intensity images taken at (a) 0.5 m propagation distance and (b) 2 m propagation distance (for the grayscale $a=0.6$, $b=1.1$), and sample thickness reconstructions using (c) TIE phase retrieval from 0.5 m intensity data and (d) Fokker-Planck (FP) phase retrieval via Eq. (4) ($a=$ 2 mm, $b=$ 8 mm). Bottom row: (e) dark-field ($D$) reconstruction using the projection image captured at a propagation distance of 2 m ($a=0$, $b=7.07\times 10^{-11}$), (f) close-up regions of the green and blue boxes in (a) and (b), respectively and (g) sample thickness profiles taken over the blue and red lines in panels (c) and (d).
The IMBL’s ‘Ruby’ detector [48], used to image the sample, is equipped with a commercial macro lens with a variable focal distance, allowing for the pixel size to be set between 6.3 µm and 22.3 µm. For our experiment, the pixel size used was 12.3 µm. We moved the detector to distances of 0.5 m and 2 m downstream of the sample, and performed a scan of the sample at each propagation distance, with each scan containing $1810$ projections collected over $181^\circ$, with a rotation step of $0.1^\circ$ between each projection. These propagation distances provided a balanced trade-off between visualizing phase and dark-field effects, with the phase effects visually apparent at both distances and subtle dark-field blurring effects seen at the larger distance. We also collected $30$ images of both the flat-field and the dark-current during each scan. Each scan took approximately $15$ minutes to complete, noting that we made no particular effort to minimize radiation dose in this first proof-of-principle experiment.
Prior to analysis, the projections captured at each distance were flat-field and dark-current corrected, then resized to account for slight magnification differences using the source-to-sample distance and finally registered to sub-pixel accuracy using the ‘phase_cross_correlation’ function of the open-source python library scikit-image [49], with an upsampling factor of $1000$ used as the argument of this function. All image processing and data analysis was performed using Python3 code and X-tract software [36] on the Australian Synchrotron Compute Infrastructure (ASCI) desktop. Taking each pair of propagation-based images (Figs. 3(a)-(b)), Eqs. (4) and (5) were used to reconstruct two-dimensional projection images of the sample transmission (and hence single-material thickness) (Fig. 3(d)) and dark-field signal (Fig. 3(e)) at each projection angle. The sample thickness reconstruction calculated using the TIE phase retrieval [34] is also provided in Fig. 3(c) for comparison.
Panel (f) of Fig. 3 demonstrates the increase in strength of dark-field effects seen with increasing propagation distance, which causes a local visibility reduction in the intensity image captured at 2 m, as compared to the intensity image captured at 0.5 m [28]. This local visibility reduction is seen as regions of increased blurriness in the measured intensity data, with the edges of the wrapped fibres becoming less visible, in particular where the projected thickness of microstructures increases (e.g. bottom-left of these regions relative to top-right). Figure 3(c)/(d) shows that the sample thickness reconstructions constructed using TIE phase retrieval [34] and our present method appear to be identical by eye. However, as discussed in Leatham et al. [28], a Fokker-Planck approach to thickness reconstruction provides increased spatial resolution relative to TIE phase retrieval, since TIE phase retrieval misinterprets diffusive dark-field contrast as slowly changing sample thickness. This can be seen from panel (g), where the blue density profile (TIE) and the red density profile (Fokker-Planck) are very similar, but the blue/TIE profile is more slowly-varying on the left of the line profile where dark-field signal coming from the PMMA microspheres locally reduces contrast/resolution. The Fokker-Planck-based reconstruction of the sample transmission using Eq. (4) at each rotation angle serves a dual purpose. Firstly, the computed set of sample transmissions can be input into a computed tomography reconstruction to create a map of the linear attenuation coefficient. Secondly, the sample transmission at each rotation angle is used as a part of the calculation of the dark-field signal at each rotation angle according to Eq. (5).
Using the reconstructed sample transmission projections and the measured intensity data at 2 m, we computed the dark-field signal projections arising from SAXS, induced by the unresolved sample microstructure according to Eq. (5). It is evident from Fig. 3(e) that the PMMA microspheres generate the strongest dark-field signal, with a weaker diffuse scattering signal coming from the agarose powder and the water. We suspect that the total volume of agarose powder is not sufficient to provide a very strong dark-field signal, based on the weak contrast observed in the projected thickness image shown in Fig. 3(d), corresponding to a weak attenuation signal. Due to the heterogeneity in the shape of the particles that made up the agarose powder, the agarose powder likely did not pack as efficiently as the PMMA microspheres, meaning there were not as many air/powder interfaces as seen for the PMMA microspheres, resulting in this weaker dark-field signal. Additionally, with particle diameter ranging from 20-160 µm, fewer particles would be seen in a given projection than for the 1-10 µm PMMA microspheres. The weak dark-field signal emanating from the agarose powder may be better retrieved by using a larger propagation distance in the dark-field reconstruction. By comparison, it can be seen that the projected sample thickness is largest in the microtube containing water, from Fig. 3(d). This demonstrates the complementary nature of the sample thickness and the dark-field. Note that some fully-resolved features appear in the projected dark-field signal shown in Fig. 3. The appearance of such fully-resolved features may be due to slight mismatches in the fringes of the measured intensity images (e.g. Fig. 3(a) and (b)) and the ‘no dark-field’ estimated intensity image given by the $(\frac {1}{\Delta _{1}^{2}}-\frac {\gamma }{2k \Delta _{1}}\nabla ^{2}_{\perp })t(x,y;\Theta )$ term in Eq. (5). Additionally, such fully-resolved features may provide a weak edge-scattering signal, which can contribute to the projected dark-field signal here, since edge-scattering is not explicitly taken into account in our present analysis. Fully-resolved features, and any background signal more generally, could potentially be removed from the projected dark-field signal by including in the field of view a region which is known to not scatter x-rays, or regions where the sample is absent. Such regions can act as a reference point of zero dark-field when reconstructing the dark-field in projection, which could be utilized by subtracting the mean value of the reconstructed dark-field signal in such regions [40]. Note that since the dark-field signal is retrieved by inverting the Laplacian operator, the point of zero dark-field is ambiguous without such a potential reference point.
We then used a Fast Filtered-Back-Projection based on a grid reconstruction tomography scheme [50] with a Hamming window to create maps of the linear attenuation coefficient and Fokker-Planck linear diffusion coefficient from the full set of projections, as outlined in the method section of this paper, via Eq. (8). Slices of the linear attenuation coefficient and Fokker-Planck linear diffusion coefficient through the bottom and middle of our sample can be found in Fig. 4.
Fig. 4. CT reconstructions of our sample. Top row: slices of the linear attenuation coefficient reconstructed using (a) TIE phase-retrieval (bottom slice), (b) Fokker-Planck (FP) phase-retrieval Eq. (4) (bottom slice), (c) TIE phase-retrieval (middle slice) and (d) FP phase-retrieval Eq. (4) (middle slice). The grayscales are shown with $a = 0\;\textrm{cm}^{-1}$, $b = 0.48\;\textrm{cm}^{-1}$ for (a)-(d). Bottom row: (e) bottom slice of the Fokker-Planck linear diffusion coefficient, (f) middle slice of the Fokker-Planck linear diffusion coefficient ($a = 0\;\textrm{cm}^{-1}$, $b = 7.53{\times}10^{-11}\;\textrm{cm}^{-1}$) and (g) profiles taken over the blue and red lines in (c) and (d), showing a spatial resolution improvement in the reconstruction of the linear attenuation coefficient when using the FP method compared to TIE phase retrieval.
To the eye, the CT slices of the linear attenuation coefficient shown in Fig. 4(a) and (b), as well as in (c) and (d), appear identical. However, as can be seen from the intensity profile given in panel (g) of Fig. 4, there is a spatial resolution improvement in the reconstruction of the linear attenuation coefficient when using the Fokker-Planck method described here compared to TIE-based phase retrieval, which comes directly from the increase in the spatial resolution of the projected sample thickness. This can also be seen by noting that one can get from the Fokker-Planck slice (Fig. 4(d)) to the TIE slice (Fig. 4(c)) by blurring with a two-dimensional symmetric Gaussian function with a standard deviation of 1-2 pixels. This increase in spatial resolution is ultimately a consequence of the Fokker-Planck approach taken to reconstructing the sample transmission (as in our method), and hence sample thickness, which unlike TIE phase retrieval, distinguishes between phase and dark-field effects. It is worth reiterating that the increase in spatial resolution found here is only due to separating out phase and dark-field effects using two intensity images. This is in contrast to how spatial resolution improvement in CT reconstructions from TIE phase retrieval has been achieved in the past, primarily through using deconvolution filters to qualitatively boost spatial resolution [35]. Such an approach could be added to further increase spatial resolution, as appropriate. The use of discrete mathematics [43,44], through the discrete Fourier transform, may also further enhance the spatial resolution increase seen in our Fokker-Planck approach. The other key result that arises from explicitly taking dark-field effects into account is the ability of our method to tomographically reconstruct a dark-field signal in three dimensions, through the Fokker-Planck linear diffusion coefficient. The reconstructed slices of the Fokker-Planck linear diffusion coefficient, Fig. 4(e) and Fig. 4(f), highlight a strong diffusive signal from the PMMA microspheres through both the top and the middle of the sample, while the water and the agarose powder provide a much weaker signal. The diffusive signal seen across the tube containing the PMMA microspheres is relatively uniform, which reflects the relatively uniform packing of the microspheres. The Fokker-Planck linear diffusion coefficient provides information about the sample which is complementary to the linear attenuation coefficient, much in the same way that the projected dark-field signal is complementary to the projected sample thickness. The non-uniform background seen within the sample holder may be due to slight mismatches in the high-contrast fringes within the measured intensity images and the ‘no dark-field’ estimated intensity image in projection, which are then back-projected. Additionally, image noise may also contribute to this non-uniform background, as image noise can induce a dark-field signal [40]. The apparent dark-field signal seen from the PMMA microtubes is likely an edge signal. This apparent dark-field signal may be due to the absence of a reference point of zero dark-field, as well as slight mismatches in the intensity fringes (described above), across the PMMA microtubes.
4. Discussion
This paper presents a new CT phase and dark-field reconstruction algorithm for two-distance propagation-based imaging, which builds upon the Fokker-Planck-based retrieval method of Leatham et al. [28]. To our knowledge, the use of propagation-based imaging for diffusive dark-field computed tomography has only been demonstrated in papers by Gureyev et al. [26] and Aminzadeh et al. [27]. A comparison between the method presented in this manuscript and that presented in Refs. [26] and [27] would be an interesting pursuit for further investigation. Of interest in future work would also be a comparison of the method presented here to dark-field CT methods which make use of optical elements, such as those methods utilizing grating interferometry [51] or analyzer-based imaging [52]. As some preliminary remarks on this front, we note that dark-field CT methods employing additional optical elements in the experimental setup have traditionally required multiple sample exposures, in the vicinity of seven exposures, associated with a higher sample dose (noting that a sliding window approach [53] can help with this). Our method, on the other hand, only needs two sample exposures for each tomographic projection, and hence may provide a dose-saving advantage over such dark-field CT methods. The impact on spatial resolution would have to be investigated. Additionally, with regards to spatial coherence requirements, our method requires that the characteristic length scale associated with source-size blurring be less than the characteristic length scale associated with diffuse scattering due to sample-induced SAXS, so that the diffusive dark-field signal associated with position-dependent SAXS can be detected [28].
In order to transition the method presented in this manuscript to clinical or industrial settings, potential optimization conditions would need to be addressed. Perhaps the biggest limiting factor of our method is the need for two sample exposures at each rotation angle, as this affects the total scan time. Ideally in the clinical and industrial settings, one would like to minimize sample dose, and for the scan to be as quick as possible, so that a large number of samples can be examined in a given amount of time. One way in which the need for two sample exposures could be removed in our method is by utilizing an imaging setup with either a beam-splitting geometry or a semi-transparent detector [54]. Two detectors placed at the desired propagation distances would be required in both geometries, with an object capable of reflecting and transmitting the x-rays being inserted after the first detector in the beam-splitting geometry. The first detector would need to be semi-transparent in the semi-transparent detector geometry, so as to allow ideally 50% of the x-rays to pass through to the second detector. In each geometry, the two detectors would be fixed at the chosen propagation distances and so in principle there should be no image artifacts arising from misalignment between the two propagation distances. Each of these imaging geometries would allow for propagation-based intensity images to be captured at two propagation distances at the same time, i.e. with a single sample exposure. These imaging geometries may also allow our method to be integrated with previous work demonstrating how a hundred or even a thousand CT projections can be captured per second in phase-contrast computed tomography [55,56], thereby enabling time-resolved tomography studies. Fast dark-field tomography could alternatively be achieved by applying the $\mu _D$ approach described here to a recently-published variation on Fokker-Planck propagation-based dark-field imaging, where two energies are used instead of two propagation distances [29].
Propagation-based imaging and in particular, the TIE phase retrieval method of Paganin et al. [34], has been applied not only in the context of x-ray imaging, but also to visible-light microscopy [57], electron microscopy [58] and neutron imaging [42]. Given that the phase and dark-field retrieval method of Leatham et al. [28], which underpins the method presented in this paper, is a diffusive generalization of the method of Paganin et al. and makes use of the same imaging setup, it stands to reason that the method of Leatham et al. and by extension the method contained in this paper, may also have a domain of validity that extends beyond x-rays and into visible-light, electron and neutron tomography.
5. Directions for future research
Here we outline several possible directions for future research based on the work presented in this paper. The first of these potential directions is to generalize our CT reconstruction method to samples composed of arbitrarily many materials. For such samples, the key quantities to be reconstructed are the attenuation of the illuminating wavefield, the phase shift of the illuminating wavefield induced by the sample, and the dark-field signal. This could be achieved by using three different propagation distances to disentangle the attenuation, phase and dark-field signals in Eq. (1) in projection. One could then use standard CT techniques to create maps of the sample linear attenuation coefficient, $\mu (x,y,z)$, the real decrement of the complex refractive index of the sample, $\delta (x,y,z)$, and the Fokker-Planck linear diffusion coefficient, $\mu _D(x,y,z)$.
A second possible direction would be to extend the CT reconstruction method presented in this paper to the case of diffusion tensor tomography by considering samples that contain sub-pixel features that scatter in a preferred direction, i.e. samples that produce a ‘directional dark-field’ signal [59,60]. This could be achieved by replacing the single diffusion coefficient $D(x,y)$ with a dimensionless rank-two symmetric diffusion tensor [25,40]:
Additionally, the method presented in this paper might be applied to lung imaging. Recent studies, such as that by Willer et al. [9], have demonstrated that dark-field chest x-rays have the potential to differentiate between healthy and diseased lung tissue in humans, while conventional chest x-rays cannot. Propagation-based methods are well-suited to smaller samples, with a number of small-animal biomedical research studies using a propagation-based set-up to better understand lung health [61–63]. For this reason, it would be an interesting avenue for future research to investigate the ability of our CT method to (i) differentiate between healthy and diseased lung tissue using animal subjects for a variety of different lung disease models, and (ii) to quantify the extent of the disease. Our method could also be used in biomedical research, for example to measure the effectiveness of treatments for lung diseases [64], or to assess the severity of disease [65].
A further potential interesting avenue for future research would be to investigate how the three-dimensional map of the Fokker-Planck linear diffusion coefficient could be used to quantify properties of the sample microstructure. Such properties could include the size or surface-area-to-volume ratio of the microstructure, $a(x,y,z)$ [40]. A reconstruction of such properties can be seen as a ‘removal’ of the influence of the x-ray beam, as the Fokker-Planck linear diffusion coefficient encodes both the influence of the x-ray beam used to illuminate the sample and the properties of the sample, while the surface-area-to-volume ratio, for example, is purely a sample property. Since microstructure size may also be determined from the sample linear diffusion coefficient ($\mu _{d}$) in grating-based imaging [47], it is natural to wonder whether there is a correspondence between propagation-based CT as described in this paper and imaging methods which use reference patterns. Indeed, since dark-field signal in the former context and visibility reduction in the latter context are related via Eq. (10), and $\mu _{d}$ can be obtained via [5]
we have the correspondence as shown in Fig. 5 below. Note that converting a dark-field signal to a visibility reduction can be seen as ‘adding in’ the properties of the reference pattern/mask to the x-ray beam and sample properties, while the process of inferring sample properties from $\mu _{d}$ can be seen as a ‘removal’ of the influence of the x-ray beam and reference pattern/mask properties, leaving only the sample properties.Fig. 5. Correspondence between dark-field signal in the context of propagation-based imaging and visibility reduction in imaging methods that make use of reference patterns [5,47]. This figure also indicates the relationship between dark-field signal/visibility reduction and sample microstructure properties, such as surface-area-to-volume ratio.
To summarize, this paper has presented a computed tomography extension of the phase and dark-field retrieval method of Ref. [28], based on the Fokker-Planck generalization of the transport-of-intensity equation, allowing for the reconstruction of the sample linear attenuation coefficient and Fokker-Planck linear diffusion coefficient using the experimental technique of propagation-based imaging. In particular, we introduced a new measure, $\mu _{D}$, to quantify how much a given voxel of the sample will diffuse the x-ray beam, independent of the experimental set-up used to capture this measurement. The method presented provides an improvement in spatial resolution for the reconstruction of the sample linear attenuation coefficient CT slices compared to TIE-based phase retrieval. Additionally, our method provides complementary sample information about unresolved microstructure through the Fokker-Planck linear diffusion coefficient, providing one of the first demonstrations of x-ray diffusive dark-field tomography without optics. The resulting map of the Fokker-Planck linear diffusion coefficient could potentially be used in the future to infer sample measurements, such as the surface-area-to-volume ratio of the spatially unresolved sample microstructure.
Funding
Australian Government (Research Training Program (RTP) Scholarship); Australian Research Council (DP230101327, FT180100374).
Acknowledgments
This research was undertaken on the Imaging and Medical Beamline at the Australian Synchrotron, part of ANSTO, under proposal 18642. In particular, we thank Dr. Daniel Hausermann, Dr. Chris Hall and Dr. Anton Maksimenko for their help with the experiment.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. F. Natterer, The Mathematics of Computerized Tomography (Society for Industrial and Applied Mathematics, 2001).
2. M. Endrizzi, “X-ray phase-contrast imaging,” Nucl. Instrum. Methods Phys. Res., Sect. A 878, 88–98 (2018). [CrossRef]
3. F. Pfeiffer, M. Bech, O. Bunk, et al., “Hard-x-ray dark-field imaging using a grating interferometer,” Nat. Mater. 7(2), 134–137 (2008). [CrossRef]
4. Z. Wang, K. Kang, Z. Huang, et al., “Quantitative grating-based x-ray dark-field computed tomography,” Appl. Phys. Lett. 95(9), 094105 (2009). [CrossRef]
5. M. Bech, O. Bunk, T. Donath, et al., “Quantitative x-ray dark-field computed tomography,” Phys. Med. Biol. 55(18), 5529–5539 (2010). [CrossRef]
6. F. Prade, K. Fischer, D. Heinz, et al., “Time resolved x-ray dark-field tomography revealing water transport in a fresh cement sample,” Sci. Rep. 6(1), 29108 (2016). [CrossRef]
7. R. Burkhardt, T. Gora, A. A. Fingerle, et al., “In-vivo X-ray dark-field computed tomography for the detection of radiation-induced lung damage in mice,” Phys. Imaging Radiat. Oncol. 20, 11–16 (2021). [CrossRef]
8. B. K. Blykers, C. Organista, M. N. Boone, et al., “Tunable x-ray dark-field imaging for sub-resolution feature size quantification in porous media,” Sci. Rep. 11(1), 18446 (2021). [CrossRef]
9. K. Willer, A. A. Fingerle, W. Noichl, et al., “X-ray dark-field chest imaging for detection and quantification of emphysema in patients with chronic obstructive pulmonary disease: a diagnostic accuracy study,” Lancet Digit. 3(11), E733–E744 (2021). [CrossRef]
10. F. T. Gassert, T. Urban, A. Kufner, et al., “Dark-field X-ray imaging for the assessment of osteoporosis in human lumbar spine specimens,” Front. Physiol. 14, 1217007 (2023). [CrossRef]
11. M. Viermetz, N. Gustschin, C. Schmid, et al., “Initial characterization of dark-field CT on a clinical gantry,” IEEE Trans. Med. Imaging 42(4), 1035–1045 (2023). [CrossRef]
12. E. Pagot, P. Cloetens, S. Fiedler, et al., “A method to extract quantitative information in analyzer-based x-ray phase contrast imaging,” Appl. Phys. Lett. 82(20), 3421–3423 (2003). [CrossRef]
13. W. Yashiro, Y. Terui, K. Kawabata, et al., “On the origin of visibility contrast in x-ray Talbot interferometry,” Opt. Express 18(16), 16890–16901 (2010). [CrossRef]
14. F. Schaff, F. Prade, Y. Sharma, et al., “Non-iterative directional dark-field tomography,” Sci. Rep. 7(1), 3307 (2017). [CrossRef]
15. K. S. Morgan, D. M. Paganin, and K. W. Siu, “X-ray phase imaging with a paper analyzer,” Appl. Phys. Lett. 100(12), 124102 (2012). [CrossRef]
16. S. Berujon, H. Wang, and K. Sawhney, “X-ray multimodal imaging using a random-phase object,” Phys. Rev. A 86(6), 063813 (2012). [CrossRef]
17. H. Wang, Y. Kashyap, and K. Sawhney, “Quantitative X-ray dark-field and phase tomography using single directional speckle scanning technique,” Appl. Phys. Lett. 108(12), 124102 (2016). [CrossRef]
18. M. C. Zdora, “State of the art of x-ray speckle-based phase-contrast and dark-field imaging,” J. Imaging 4(5), 60 (2018). [CrossRef]
19. M.-C. Zdora, P. Thibault, H. Deyhle, et al., “Tunable X-ray speckle-based phase-contrast and dark-field imaging using the unified modulated pattern analysis approach,” J. Instrum. 13(05), C05005 (2018). [CrossRef]
20. S. J. Alloo, D. M. Paganin, K. S. Morgan, et al., “Dark-field tomography of an attenuating object using intrinsic x-ray speckle-tracking,” J. Med. Imag. 9(3), 031502 (2022). [CrossRef]
21. R. Smith, F. De Marco, L. Broche, et al., “X-ray directional dark-field imaging using Unified Modulated Pattern Analysis,” PLoS One 17(8), e0273315 (2022). [CrossRef]
22. F. A. Vittoria, M. Endrizzi, P. C. Diemoz, et al., “X-ray absorption, phase and dark-field tomography through a beam tracking approach,” Sci. Rep. 5(1), 16318 (2015). [CrossRef]
23. M. Endrizzi, P. C. Diemoz, T. P. Millard, et al., “Hard X-ray dark-field imaging with incoherent sample illumination,” Appl. Phys. Lett. 104(2), 024106 (2014). [CrossRef]
24. A. Doherty, S. Savvidis, C. Navarrete-León, et al., “Edge-illumination x-ray dark-field tomography,” Phys. Rev. Appl. 19(5), 054042 (2023). [CrossRef]
25. D. M. Paganin and K. S. Morgan, “X-ray Fokker–Planck equation for paraxial imaging,” Sci. Rep. 9(1), 17537 (2019). [CrossRef]
26. T. E. Gureyev, D. M. Paganin, B. Arhatari, et al., “Dark-field signal extraction in propagation-based phase-contrast imaging,” Phys. Med. Biol. 65(21), 215029 (2020). [CrossRef]
27. A. Aminzadeh, B. D. Arhatari, A. Maksimenko, et al., “Imaging breast microcalcifications using dark-field signal in propagation-based phase-contrast tomography,” IEEE Trans. Med. Imaging 41(11), 2980–2990 (2022). [CrossRef]
28. T. A. Leatham, D. M. Paganin, and K. S. Morgan, “X-ray dark-field and phase retrieval without optics, via the Fokker-Planck equation,” IEEE Trans. Med. Imaging 42(6), 1681–1695 (2023). [CrossRef]
29. J. N. Ahlers, K. M. Pavlov, M. J. Kitchen, et al., “X-ray dark-field via spectral propagation-based imaging,” arXiv, arXiv:2309.15874 (2023). [CrossRef]
30. A. Snigirev, I. Snigireva, V. Kohn, et al., “On the possibilities of x-ray phase contrast microimaging by coherent high-energy synchrotron radiation,” Rev. Sci. Instrum. 66(12), 5486–5492 (1995). [CrossRef]
31. P. Cloetens, R. Barrett, J. Baruchel, et al., “Phase objects in synchrotron radiation hard x-ray imaging,” J. Phys. D: Appl. Phys. 29(1), 133–146 (1996). [CrossRef]
32. S. W. Wilkins, T. E. Gureyev, D. Gao, et al., “Phase-contrast imaging using polychromatic hard x-rays,” Nature 384(6607), 335–338 (1996). [CrossRef]
33. M. R. Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am. 73(11), 1434–1441 (1983). [CrossRef]
34. D. Paganin, S. C. Mayo, T. E. Gureyev, et al., “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc. 206(1), 33–40 (2002). [CrossRef]
35. T. Weitkamp, D. Haas, D. Wegrzynek, et al., “ANKAphase: software for single-distance phase retrieval from inline X-ray phase-contrast radiographs,” J. Synchrotron Rad. 18(4), 617–629 (2011). [CrossRef]
36. T. E. Gureyev, Y. Nesterets, D. Ternovski, et al., “Toolbox for advanced x-ray image processing,” Proc. SPIE 8141, 81410B (2011). [CrossRef]
37. M. A. Beltran, D. M. Paganin, K. Uesugi, et al., “2D and 3D X-ray phase retrieval of multi-material objects using a single defocus distance,” Opt. Express 18(7), 6423–6436 (2010). [CrossRef]
38. M. A. Beltran, D. M. Paganin, and D. Pelliccia, “Phase-and-amplitude recovery from a single phase contrast image using partially spatially coherent radiation,” J. Opt. 20(5), 055605 (2018). [CrossRef]
39. K. S. Morgan and D. M. Paganin, “Applying the Fokker–Planck equation to grating-based x-ray phase and dark-field imaging,” Sci. Rep. 9(1), 17465 (2019). [CrossRef]
40. D. M. Paganin, D. Pelliccia, and K. S. Morgan, “Paraxial diffusion-field retrieval,” Phys. Rev. A 108(1), 013517 (2023). [CrossRef]
41. Ya. I. Nesterets, “On the origins of decoherence and extinction contrast in phase-contrast imaging,” Opt. Commun. 281(4), 533–542 (2008). [CrossRef]
42. D. M. Paganin, M. Sales, P. M. Kadletz, et al., “Revisiting neutron propagation-based phase-contrast imaging and tomography: use of phase retrieval to amplify the effective degree of brilliance,” Phys. Rev. Appl. 19(3), 034005 (2023). [CrossRef]
43. D. M. Paganin, V. Favre-Nicolin, A. Mirone, et al., “Boosting spatial resolution by incorporating periodic boundary conditions into single-distance hard-x-ray phase retrieval,” J. Opt. 22(11), 115607 (2020). [CrossRef]
44. J. A. Pollock, K. S. Morgan, L. C. P. Croton, et al., “Precise phase retrieval for propagation-based images using discrete mathematics,” Sci. Rep. 12(1), 18469 (2022). [CrossRef]
45. D. M. Paganin, Coherent X-ray Optics (Oxford University Press, 2006).
46. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE Press, 1988).
47. S. K. Lynch, V. Pai, J. Auxier, et al., “Interpretation of dark-field contrast and particle size selectivity in grating interferometers,” Appl. Opt. 50(22), 4310–4319 (2011). [CrossRef]
48. C. Hall, D. Hausermann, A. Maksimenko, et al., “Detectors for the Imaging and Medical Beam Line at the Australian Synchrotron,” J. Instrum. 8(06), C06011 (2013). [CrossRef]
49. S. van der Walt, J. L. Schönberger, J. Nunez-Iglesias, et al., and the scikit-image contributors, “scikit-image: image processing in Python,” (2023). https://scikit-image.org/docs/stable/api/skimage.registration.html.
50. B. A. Dowd, G. H. Campbell, R. B. Marr, et al., “Developments in synchrotron x-ray computed microtomography at the National Synchrotron Light Source,” Proc. SPIE 3772, 224–236 (1999). [CrossRef]
51. N. Bevins, J. Zambelli, Z. Qi, et al., “X-ray dark-field computed tomography using a grating interferometer setup,” Proc. SPIE 7622, 76220P (2010). [CrossRef]
52. J. G. Brankov, M. N. Wernick, Y. Yang, et al., “A computed tomography implementation of multiple-image radiography,” Med. Phys. 33(2), 278–289 (2006). [CrossRef]
53. I. Zanette, T. Weitkamp, G. Le Duc, et al., “X-ray grating-based phase tomography for 3D histology,” RSC Adv. 3(43), 19816–19819 (2013). [CrossRef]
54. R. P. Carnibella, A. Fouras, and M. J. Kitchen, “Single-exposure dual-energy-subtraction X-ray imaging using a synchrotron source,” J. Synchrotron Radiat. 19(6), 954–959 (2012). [CrossRef]
55. F. García-Moreno, P. H. Kamm, T. R. Neu, et al., “Using X-ray tomoscopy to explore the dynamics of foaming metal,” Nat. Commun. 10(1), 3762 (2019). [CrossRef]
56. F. García-Moreno, P. H. Kamm, T. R. Neu, et al., “Tomoscopy: time-resolved tomography for dynamic processes in materials,” Adv. Mater. 33(45), 2104659 (2021). [CrossRef]
57. P. K. Poola and R. John, “Label-free nanoscale characterization of red blood cell structure and dynamics using single-shot transport of intensity equation,” J. Biomed. Opt. 22(10), 106001 (2017). [CrossRef]
58. A. C. Y. Liu, D. M. Paganin, L. Bourgeois, et al., “Projected thickness reconstruction from a single defocused transmission electron microscope image of an amorphous object,” Ultramicroscopy 111(8), 959–968 (2011). [CrossRef]
59. T. H. Jensen, M. Bech, O. Bunk, et al., “Directional x-ray dark-field imaging,” Phys. Med. Biol. 55(12), 3317–3323 (2010). [CrossRef]
60. T. H. Jensen, M. Bech, I. Zanette, et al., “Directional x-ray dark-field imaging of strongly ordered systems,” Phys. Rev. B 82(21), 214103 (2010). [CrossRef]
61. M. L. Siew, A. B. te Pas, M. J. Wallace, et al., “Positive end-expiratory pressure enhances development of a functional residual capacity in preterm rabbits ventilated from birth,” J. Appl. Physiol. 106(5), 1487–1493 (2009). [CrossRef]
62. C. S. Stahr, C. R. Samarage, M. Donnelley, et al., “Quantification of heterogeneity in lung disease with image-based pulmonary function testing,” Sci. Rep. 6(1), 29438 (2016). [CrossRef]
63. R. Gradl, M. Dierolf, B. Günther, et al., “In vivo dynamic phase-contrast x-ray imaging using a compact light source,” Sci. Rep. 8(1), 6788 (2018). [CrossRef]
64. K. S. Morgan, M. Donnelley, N. Farrow, et al., “In vivo x-ray imaging reveals improved airway surface hydration after a therapy designed for cystic fibrosis,” Am. J. Respir. Crit. Care Med. 190(4), 469–471 (2014). [CrossRef]
65. M. Frank, F. T. Gassert, T. Urban, et al., “Dark-field chest X-ray imaging for the assessment of COVID-19-pneumonia,” Commun. Med. 2(1), 147 (2022). [CrossRef]
