Open Access
Add to BibSonomyAbstract
Here we present a method for the removal of multi-material artifacts which occur during the application of a single material phase retrieval procedure to X-ray tomographic data sets. For the phase retrieval we chose the most common method which is the single material filter. The correction method which we describe in the following has been designed for samples consisting of three distinct materials, hence effectively two different material interfaces. Furthermore the material phase with the strongest X-ray interaction needs to show sufficient absorption in order to allow for segmenting this phase through application of a grey value threshold. If these conditions are fulfilled the method is easy to apply through post processing as is shown for the volume images of two sample types.
© 2015 Optical Society of America
1. Introduction
Computed tomography is a method for producing X-ray volume images of a given sample. Be it for the medical diagnosis of a patient or for the non-destructive testing of materials, the underlying method is the same: A series of projection views onto the sample over an angular range from 0° to 180° (or 360°) is converted to a three-dimensional volume image through a dedicated back-projection algorithm. If, hereby, two material phases show too little difference in X-ray absorption, or if the whole sample is largely transparent to the X-rays, volume information can only be obtained through phase contrast techniques [1, 2]. The latter highlights the 3D image of the sample at its materials interfaces, thus enhancing the visibility of structural details which would remain hidden otherwise. Phase contrast imaging in combination with a dedicated phase retrieval filter [3–6] yields volume images (called holo-tomography) which resemble X-ray absorption CT, yet with a stronger materials contrast, which has been retrieved from the phase contrast.
At synchrotron facilities the user can choose at which propagation (sample-detector) distance the inline phase contrast data is recorded, and also whether or not several such distances are used for the phase retrieval. Here we restrict our view to a single distance phase contrast measurement for which the propagation length d fulfills the Weitkamp criterion, hence d ≤ zc = (2Δx)2/λ, with Δx the detector pixel size and λ the X-ray wavelength [7]. This single distance near field scenario is the most common at synchrotrons and in the laboratory for which D. Paganin developed a simple linear phase retrieval method [8], which is the most robust of several similar algorithms [9].
This method (like others) has however the drawback of retrieving the phase correctly only for one material interface, which is why it is often referred to as single material phase retrieval. In the case of several distinct material interfaces one has to choose which will be treated correctly, while imaging artifacts are to be expected at the remaining others (e.g. ternary metallic alloys, biological samples with bone, soft tissue and air). Figure 1(b) and 1(c) show such imaging artifacts in the CT slice of a rose stem which may arise in a sample with three materials if the single material phase retrieval is adapted either to the softer 1(b) or to the hard tissue 1(c).
Fig. 1 Axial CT-slices of a dry rose stem. (a) Overview over the upper right quadrant (1.6 mm wide), (b, c) close up onto the structure showing typical multi material artifacts, examples for artifacts indicated by arrows. If the phase retrieval is applied to the softer tissue (cells), it causes a blurring around the harder calcium oxalate (CaOx, black) crystals, some of which lie above or below the displayed slice (b). If it is applied to the CaOx, phase contrast artifacts (fringe contrast) will remain for the softer material (gray) (c).
Another approach to phase retrieval for multi material samples has been presented in [10,11]. In that work several phase retrieval filters are applied to reconstruct a series of volume images, each optimized for one material interface. Then in each volume the corresponding materials phase is masked and the final image is composed of the sum of these masked volumes. In principle, this approach is more complete and also more complex than ours, but in practice the results are of the same quality.
For the sake of completeness we would like to point out that phase retrieval from multi distance CT measurements are also prone to multi material artifacts. Since multi distance holotomograms are generally retrieved by a contrast transfer function method (e.g. the so called mixed approach [5]) dedicated multi material corrections have to be used [12].
The method which we will introduce allows for the phase retrieval of multiple materials interfaces, which leads to the title “multi material phase retrieval”. We chose an approach which corrects the multi material artifact by post-processing the volume image with a 3D Fourier filter. The method is simple, quick and effective although it relies on a relatively crude assumption.
2. Theory and methods
2.1. Holo-tomography
In X-ray holo-tomography the volume image is built from the following steps: (1) A series of single distance inline phase contrast projection images is recorded for sample rotations from 0° to 180°, (2) the single material phase retrieval filter is applied in the same manner to all projections, and (3) the volume image is built from the filtered projections by numerical back projection methods. Our proposed multi material correction applies after step (3), i.e. it is a volume image post processing step.
A phase retrieved projection image is calculated from the measured intensities I(x, y) at the pixel coordinates x and y, as well as from the flat field image I0(x, y) according to [8]:
with the (radially symmetric) spatial frequency u (reciprocal variable to the sample coordinate r2 = |r⃗|2), ℱ2D the two-dimensional Fourier transform, μ(r⃗) the absorption coefficient at point r⃗ in the sample and ∫ dl the line integral along a straight line from the source to the detector pixel. The phase retrieval parameter p2 = 4π2dδ/μ has to be chosen for a given material and X-ray wavelength. Here, δ is the real increment of the refractive index n = 1 − δ + iβ and d is the propagation distance. We assume a material-air interface here, for a material-material interface retrieval see [10]. In that case we obtain p2 = 4π2dΔδ/Δμ with Δδ and Δμ the differences between the two materials, but the filter also differs in other aspects.We abstain from using a single material assumption (δ ∝ μ) to rescale the resulting images to δ, because this would set the values of δ for all materials except one incorrectly. Furthermore if a single material Fourier filter is used, the material interfaces of a multi material sample will still be distorted, but the gray values in homogeneous areas will be in absorption units, and correct for all materials.
Both single and two materials phase retrievals include a (Fourier) filter kernel of the form:
which applies through multiplication in reciprocal space. Since Kp(u) is radially symmetric it is sufficient to use one reciprocal coordinate u2 = ux2 +uy2 in two and u2 = ux2 +uy2 +uz2 in three dimensions. This kernel has the properties of a low-pass filter, hence its effect is to smooth or blur images. In this case, the blurring inverts the prior oversharpening from the phase contrast effect.The phase retrieval parameter p determines the strength of the retrieval procedure, higher values yield a stronger effect. In general, material interfaces with a larger difference in X-ray absorption require smaller values of p, e.g. the value of p for calcium-air is smaller than the value for carbon-air. In practice, the value of p is chosen empirically according to the optimal visual quality of the result.
From the inverse Fourier transform of Eq. (2) we can calculate the filter kernel in real-space (convolution kernel):
From this equation we can deduce the range of the filter from the value of p, which scales the exponential term. Accordingly the necessary padding length for the application of a fast Fourier transform (FFT) can be estimated through p as well.
2.2. Preconditions
The aim is to deduce a method which — by applying a volume filter — removes the multi material image artifacts which are created during the application of the single material phase retrieval. By combining the single material phase retrieval with this correction method we obtain multi material phase retrieved volume images.
First we assume that the sample contains three different materials and effectively two types of materials interfaces for each of which a phase retrieval parameter exists. For one materials interface the difference in X-ray absorption shall be weak, the third material shall display a higher absorption. Both transitions from this strongly absorbing material to the two materials of weaker absorption can be treated well by a single material phase retrieval approach with the same parameter p.
As the initial motivation for the application of phase contrast is the retrieval of the interface between the two less absorbing materials in the sample these assumptions cover the majority of the real cases where multi materials artifacts occur. If the sample consisted only of materials which have strongly different X-ray attenuation, they would be sufficiently resolved by a mere absorption CT.
2.3. Correction method
The two materials interfaces shall be labeled A and B, whereby A designates the interface with the stronger attenuation contrast. Before back-projecting, the single material phase retrieval filter is applied with a relatively small value of p (for the interface A).
After back-projecting we now attempt to apply a volume filter in order to alter the strength of the phase retrieval, so that it matches the material interface B. In practice this filter shall later solely apply to those parts of the volume which contain materials interfaces of type B in order to obtain the correct phase retrieval for both interfaces.
The reason for choosing a volume filter is that in the projection images both materials interfaces A and B are superimposed and thus cannot be masked. The justification for applying phase retrieval to a volume image stems from the exchangeability of the inverse Radon transform and a Fourier filter. Yet, normally phase retrieval in Eq. (1) is applied to the projection images and before taking the logarithm of the intensities. By applying phase retrieval as a volume filter we basically linearize the logarithm in Eq. (1) for the sample thickness which is relatively crude but works well in practice [13].
In order to alter the strength of the phase retrieval, we deduce a new filter kernel from Eq. (2). First we invert the phase retrieval for interface A by dividing by KA(u) in Fourier space then we multiply with KB(u) to obtain the phase retrieval shift filter. The new filter kernel is
The corresponding inverse Fourier transform gives the new (real space) convolution kernel:
where δ(r) is the Dirac-δ -distribution. The exponential term in this equation should hence give the functional form of the blurring artifacts (see Fig. 1(b)) created by applying too much phase retrieval for a certain material interface.In order to apply this filter selectively to one materials interface (B) only, the volume image has to be masked for the other interface (A). Generally the necessary image segmentation can pose a problem yet in many examples simple threshold segmentation suffices. The binary mask MA(r⃗) is generated from the volume image V(r⃗) by applying the threshold t in the following manner:
In order to remove noise from the segmentation a binary opening algorithm is applied to MA(r⃗) which introduces also minor smoothing to the mask (structuring element: 3×3×3 cube). Additionally one could further smooth out the mask (e.g. by convolving with a Gaussian kernel corresponding to the resolution of the imaging system) in order to remove unrealistic sharp edges at the materials interface A. In this work we abstain from this further smoothing. The remaining volume segment (“inverse mask”), shall be referred to as MB(r⃗) = 1 − MA(r⃗).
From this idea stem two ways by which multi material artifacts can be removed, depending on whether the interface A or B is retrieved prior to application of the method. In both cases the Fourier filter in Eq. (4) is applied locally. We shall call them forward and backward correction.
2.4. Forward correction
Here, the single material holo-tomography for materials interface A (using a smaller value for p) is built first: VA(r⃗) is optimized for the material interface at the material with the strongest X-ray attenuation. We thereby evade the blurring artifacts which could be seen in Fig. 1(b). Then we apply segmentation to create MA and MB (inverse mask only containing material interfaces of type B). Phase contrast edges included by the threshold are filtered out by the binary opening. Then we amplify the phase retrieval in MB by applying KA→B(u). The remaining fringe contrast in MB is thus removed, the result is added to VAMA. The method summarizes to:
2.5. Backward correction
The starting point here is the holo-tomography VB(r⃗) for which stronger phase retrieval has been applied. Here, the segmentation of material interfaces of type A poses a problem because of the additional blurring artifacts (cf. Fig. 1(b) which extend into regions of materials interface B. The correct filter would need to be applied to the blurring artifacts but exclude the superimposed structure. As this is impossible, a somewhat smaller portion of the volume is segmented instead, only including the voxels of the material with the higher absorption. The backward correction summarizes to:
For this variant, the filter steps are shown in Fig. 2.
Fig. 2 The filter steps (of the backward correction) shown for the rose stem sample from Fig. 1. The gray value ranges are (0, 5.5) for (c), (−1.5, 0.5) for (d) and (0, 2.2) for the rest. The first step consists of separating the single material phase retrieved image (a) into the masked area (c) and the rest (b). Next, the correction filter KA→B(u) is applied to (b), resulting in (d). The final result (e) is achieved by summing (b) and (d). The blurring artefacts can be seen in (c) and (d) with opposite sign.
For the backward correction, the simple mask generation by threshold does not work well, because the strong blurring artifacts smeared out the interface A and some voxels are falsely assigned to MB. The mask can be improved by iterating the correction in the following way: The first correction is applied for a relatively low threshold t (removing most of the blurring). The result of this first correction is the used to generate a mask. The next correction is generated from this mask, then the threshold is raised and a new mask is calculated. This is iterated until the mask only contains the material with the highest absorption.
3. Results
The measurements shown here were recorded on the ID19 beamline at the European synchrotron radiation facility (ESRF) in Grenoble, France. For improving the image sharpness a Wiener deconvolution was used in addition to the phase retrieval.
3.1. Rose stem sample
The first sample is a dry rose stem (already shown in Fig. 1) which contains air, dry wood (mostly carbon, hence weak absorption) and embedded crystals of calcium oxalate [14] (CaOx, strong absorption). The effective pixel sampling was 1.1 μm at 20 keV photon energy.
First we tested the validity of the approximation (linearization of the log in Eq. (1)) that allows us to exchange the inverse Radon transform (= volume reconstruction) with a Fourier filter (Eq. (4)). There are no visible differences between the application of the phase retrieval prior or after the volume reconstruction. To be precise, the magnitude of the difference in terms of grey values was one order of magnitude smaller than the image noise level of the data. We conclude that the linear approximation is valid.
Figure 3(a) shows an axial slice of the single material holo-tomography which has been adjusted to the air-wood interface (pB = 63 μm): A strong blurring around the CaOx crystals (black) can be observed. Figure 3(b) shows the same holo-tomography, but optimized for the CaOx-air/wood materials interface (pA = 24 μm): No blurring occurs but the phase retrieval for the weaker materials interface (wood-air) is visibly incomplete. The result of the backward correction is shown in Fig. 3(c): Despite the uncertainties which arise from the segmentation of the mask MA the method yields a good result when the phase retrieval parameter pA is raised to 39 μm. Three iterations for the refinement of the mask were applied to yield Fig. 3(c). The result of the forward correction is shown in Fig. 3(d): Both materials interfaces are retrieved accurately within one application of Eq. (7).
Fig. 3 Axial CT-slices from the holo-tomographies of the rose stem sample. Without (top) and with (bottom) multi material correction applied. The visible gray scale range is [0.0, 2.4] for (a), (c) and (d) and [−1.5, 3.9] for (b), the gray value for the CaOx (black) is ouside of the visible gray value range. Both materials have varying densities. The image area shown is 0.72 mm wide and 0.54 mm high.
Figure 4 displays histograms of the volume grey values as well as line plots through the slices shown in Fig. 3, for the backward correction of the single material phase retrieval of VB (wood-air) (left) as well as for the forward correction of VA (CaOx-air/wood) (right). This allows for a more quantitative analysis of the filter results. The histograms show that the gray values that should be unchanged by the filter remain the same (left: lower gray values, right: higher gray values). In comparison, the backward correction reconstructs smaller gray values for the CaOx parts. The line plots should ideally be part-wise constant with identical gray values for one material. The filter results from the multi material correction are close to that and much closer than the original single material holo-tomographies.
Fig. 4 Volume histograms (top) and line plots (bottom) for the rose stem sample. The position of the line plot in the sample is indicated in Fig. 3(b). Single material holotomographies (dashed lines) for wood-air and CaOx-wood/air, results of the multi material corrections (solid lines).
3.2. Metal alloy sample
The second sample is a metallic alloy composed of aluminum and silicon (weak materials interface) as well as strongly absorbing inter-metallic phases (e.g. copper). The effective pixel sampling was 0.7 μm at 19 keV photon energy.
Figure 5 shows the two results of single material phase retrieval for the strong Fig. 5 (a) and the weak Fig. 5(b, c) materials interfaces as well. The result of the forward correction is shown in (d), based on VA (Fig. 5 (a)). For the Cu-AlSi/Si materials interface a value of pA = 7.1 μm was used for the phase retrieval, while pB = 20 μm was optimal for the AlSi-Si interface. Applying the stronger phase retrieval pB to the projections data set posed difficulties: in addition to the usual blurring artifacts, strong image artifacts (streaks) resulted from the violation of the single material criterion (b). Therefore (c) was obtained from applying KA→B (Eq. (4)) to the whole volume VA (a), effectively avoiding these artifacts.
Fig. 5 Axial CT-slices from the holo-tomographies of the metal alloy sample. Single material phase retrieved ((a), (b) and (c)), forward corrected result (d), based on (a). In addition to the blurring artifacts, there are horizontal streaks in (b), which appear only if the phase retrieval is applied on projections. These straks originate from the rotational centre of the CT, which is to the left of the image area. They are not visible in (c), where the phase retrieval is applied as a volume filter starting from the weak projection-based phase retrieval (a). Identical gray scales, Cu (black) is outside of the visible gray value range. The image area shown is 0.29 mm wide and 0.22 mm high.
Backward correction did not yield good results for the metal alloy sample and is therefore not shown here.
4. Discussion
For the rose stem sample, visual impression of the result of the backward correction is as good as the results from the forward correction: Both materials interfaces appear to be retrieved correctly without remaining blurring artifacts. For the backward correction the initial segmentation of MA was incomplete, so this result is somewhat surprising. Also, the gray values for the CaOx in the backward corrected result seem to be reconstructed too low. For the metal alloy sample, the results of the backward correction did not yield good results.
The forward correction method applies well to both samples, no remaining image artifacts are visible. This had to be expected since our approximation of linearizing the logarithm in Eq. (1) proved to be fully valid and segmentation worked well for both samples. This leads us to favor the forward correction as a more universal correction method.
For the metal alloy, the phase retrieval applied to the volume reconstruction (Eq. (4)) works better than the projections-based phase retrieval (Eq. (1)). The volume phase retrieval applies locally with a known range, therefore it cannot produce long-range artifacts as observed with the projections-based phase retrieval. Thus the volume phase retrieval may actually be the more robust algorithm in the case of a single material interface.
5. Conclusions
With the method we presented here, a correct multi material holo-tomography can be performed for a variety of samples which would yield multi material artifacts, when single material phase retrieval is applied in the common way. Our method is simple, quick and yields good visual results. It should in principle apply to even more than two materials interfaces, provided that segmentation of the latter is feasible, but this remains to be tested.
Method implementation
An implementation of the method in the python programming language can be requested from the authors per e-mail. It includes a tool with a graphical user interface in PyQt and a fast preview. To avoid memory shortage, it applies the method to overlapping segments of the full volume data, thus allowing the tool to work on a normal computer.
Acknowledgments
We would like to thank Stefan Gerth for the measurements, Alexander Rack for supervising the beamtime and Jonas Dittmann for discussion and help with the method.
References and Links
1. P. Cloetens, W. Ludwig, J. Baruchel, D. Van Dyck, J. Van Landuyt, J.P. Guigay, and M. Schlenker, “Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron radiation x rays,” Appl. Phys. Lett. 75(19), 2912–2914 (1999). [CrossRef]
2. A. Snigirev, I. Snigireva, V. Kohn, S. Kuznetsov, and I. Schelokov, “On the possibilities of x-ray phase contrast microimaging by coherent high-energy synchrotron radiation,” Rev. Sci. Instrum. 66(12) 5486–5492 (1995). [CrossRef]
3. M. Boone, Y. De Witte, M. Dierick, J. Van den Bulcke, J. Vlassenbroeck, and L. Van Hoorebeke, “Practical use of the modified Bronnikov algorithm in micro-CT,” Nucl. Instr. Meth. Phys. Res. 267(7), 1182–1186 (2009). [CrossRef]
4. A. Groso, R. Abela, and M. Stampanoni, “Implementation of a fast method for high resolution phase contrast tomography,” Opt. Express 14(18), 8103–8110 (2006). [CrossRef] [PubMed]
5. J.P. Guigay, M. Langer, R. Boistel, and P. Cloetens, “Mixed transfer function and transport of intensity approach for phase retrieval in the Fresnel region,” Opt. Lett. 32(12), 1617–1619 (2007). [CrossRef] [PubMed]
6. T. E. Gureyev, S. Mayo, S.W. Wilkins, D. Paganin, and A.W. Stevenson, “Quantitative in-line phase-contrast imaging with multienergy x rays,” Phys. Rev. Lett. 86(25), 5827 (2001). [CrossRef] [PubMed]
7. T. Weitkamp, D. Haas, D. Wegrzynek, and A. Rack, “ANKAphase: software for single-distance phase retrieval from inline X-ray phase-contrast radiographs,” J. Synchrotron Radiation 18(4), 617–629 (2011). [CrossRef]
8. D. Paganin, S.C. Mayo, T.E. Gureyev, P.R. Miller, and S.W. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” Journal of Microscopy 206(1), 33–40 (2002). [CrossRef] [PubMed]
9. A. Burvall, U. Lundström, P.A.C. Takman, D. H. Larsson, and H.M. Hertz, “Phase retrieval in X-ray phase-contrast imaging suitable for tomography,” Opt. Express 19(11), 10359–10376 (2011). [CrossRef] [PubMed]
10. M. A. Beltran, D. Paganin, K. Uesugi, and M.J. Kitchen, “2D and 3D X-ray phase retrieval of multi-material objects using a single defocus distance,” Opt. Express 18(7), 6423–6436 (2010). [CrossRef] [PubMed]
11. M. A. Beltran, D. Paganin, K.K.W. Siu, A. Fouras, S.B. Hooper, D.H. Reser, and M.J. Kitchen, “Interface-specific x-ray phase retrieval tomography of complex biological organs,” Phys. Med. Biol. 56(23), 7353 (2011). [CrossRef] [PubMed]
12. M. Langer, P. Cloetens, A. Pacureanu, and F. Peyrin., ”X-ray in-line phase tomography of multimaterial objects,” Opt. Lett. 37:(11),2151–2153 (2012). [CrossRef] [PubMed]
13. S. Zabler, H. Riesemeier, P. Fratzl, and P. Zaslansky, “Fresnel-propagated imaging for the study of human tooth dentin by partially coherent x-ray tomography,” Opt. Express 14(19), 8584–8597 (2006). [CrossRef] [PubMed]
14. U. Matsushima, A. Hilger, W. Graf, S. Zabler, I. Manke, M. Dawson, G. Choinka, and W.B. Herppich, “Calcium oxalate crystal distribution in rose peduncles: Noninvasive analysis by synchrotron X-ray micro-tomography,” Postharvest Biol. Tec. 72, 27–34 (2012). [CrossRef]
