Iterative reconstruction in x-ray computed laminography from differential phase measurements

Sébastien Harasse; Wataru Yashiro; Atsushi Momose

doi:10.1364/OE.19.016560

1. Introduction

Non-destructive imaging of structures inside an opaque object is classically performed by X-ray computed tomography (CT). With synchrotron radiation and the assumption of a parallel beam, it involves scanning the object over 180° around an axis normal to the beam. This however generally requires the object to be smaller than the field of view, which is not always possible depending on the setup. Of particular interest is the case of extended flat samples, for example circuit boards, whose lamellar structures would be visualized. The shape of extended flat samples is so that absorption in directions close to the object plane may be too high to perform a complete CT scan. Alternative scanning geometries are then preferred, two of which are limited angle CT and computed laminography (CL). Limited angle CT [1] simply limits the angular range of CT to angles with lower absorption. Laminography differs by rotating the object around an axis not normal to the beam path, thus avoiding orientations with large absorption. Although laminography requires rotation around 360°, it provides theoretically a more complete sampling of the object than limited angle CT, and is thus studied here.

Also known as circular tomosynthesis, the three dimensional reconstruction of an object slice from a circular scan of the object at an arbitrary angle using X-rays has been performed for a long time (see for example [2]). Although the frequency response of such scans was well understood, the increasing usage of computers has later allowed digital filtering [3] to reduce artifacts from outside of the focused plane. Such X-ray experiments using the phase information are more recent. X-ray phase CL has been studied by Helfen et al. [4, 5], using the in-line phase contrast retrieval method.

It is proposed here to use an X-ray Talbot interferometer [6] to perform the CL scan. Advantages of this imaging method over the in-line phase contrast technique include the possibility to use a detector with a standard pixel size of several microns. In comparison, the in-line phase method requires submicron pixels. Moreover, X-ray Talbot interferometry can deliver three different types of images, namely the absorption, differential phase and visibility reduction. Our purpose is also to extend the applicability of X-ray Talbot interferometry to a wider class of objects compared to CT, namely the flat laterally extended objects.

In this paper, we focus mainly on the differential phase measurements, from which it is possible to estimate by tomographic reconstruction the decrement δ from the real part of the complex refractive index 1 – δ + iβ. First we will derive in section 2 the filtered backprojection formula for CL from differential phase images. The obtained filtered backprojection operator will then be used afterwards in a more effective iterative reconstruction, presented in section 3. Iterative methods generally successfully apply to ill-posed tomographic problems with various constraints like positivity, limited support or gradient sparsity (e.g. [7]). The proposed reconstruction algorithm is based on iterative filtered backprojection [8] (IFBP), with the addition of constraints in the object support and limited range of values. The compact support assumption is a very natural addition to the reconstruction in the case of flat extended objects, as shall be discussed in section 3.1. Because the sampling resolution is limited in directions close to the normal of the object, this prior information on the support helps increase the reconstruction quality. The other constraint limits reconstructed values into a given range. This is of particular importance with our differential phase measurements, which do not contain information about the mean value of the object. Even for non differential measurements like absorption and visibility, limiting the range is still a good idea to prevent algorithms to converge to negative values for example. However, such an iterative algorithm for CL requires the processing of the full volume at once because slices are not independent as in CT. Therefore, to overcome the large dimension of the system, an implementation on a Graphics Processing Unit (GPU) has been performed, allowing the reconstruction of volumes with around 100 million voxels within a few minutes. A synchrotron experiment is presented in section 4, with the imaging of a polymer mesh by CL with a Talbot interferometer. The result obtained by filtered backprojection is compared to the iterative reconstruction.

2. X-ray computed laminography by Talbot interferometry

2.1. Scan geometry

In the present study, we assume the case of a parallel beam illumination. An extension to the cone beam case is however possible as will be discussed in section 3.6. In X-ray computed laminography using synchrotron radiation, the detector and beam direction are fixed, while the flat extended object is rotated about its plane normal direction, which is itself not included in the beam normal plane, unlike CT. The geometric relation between the beam path, object and detector is detailed in Fig. 1. The fixed Cartesian coordinate system of the detector has origin point O and axes (u⃗, v⃗, w⃗), while another coordinate system local to the object is defined with same origin O and axes (x⃗, y⃗, z⃗). The object is rotated about axis y⃗ for 360°, as described by angle θ. The object plane is tilted with a fixed acute angle α from the beam path direction. Angle α should be carefully chosen depending on the application. Smaller angles make the setup closer to the classical CT scan, with the limit case α = 0 being CT. However, using lower angles means that X-rays are almost aligned with the object plane, which is practically infeasible because of the high absorption by the sample. Low angles may also produce phase wrapping problems in differential phase images. Conversely, higher α angles produce a greater loss of resolution in the object plane normal direction.

Fig. 1 Geometry of the computed laminography scan. Angle α is fixed between 0 and 90°, whereas angle θ ranges from 0 to 360°. Base vectors (u⃗, v⃗, w⃗) are represented at the detector position for clarity, although the coordinate system origin is point O, at the object position.

Download Full Size | PDF

The relation between object coordinates and detector coordinates is defined by matrix M_θ, composition of rotation by angle α around axis x⃗ and rotation by angle θ around axis y⃗ :

M_{θ} = [\begin{matrix} cos θ & 0 & - sin θ \\ sin α sin θ & cos α & sin α cos θ \\ cos α sin θ & - sin α & cos α cos θ \end{matrix}]

A vector r⃗ = [x,y,z]^T expressed in object coordinates is transformed into detector coordinates q⃗ = [u,v,w]^T as follows :

\vec{q} = M_{θ} \vec{r}

2.2. Differential phase measurements

In the considered Talbot interferometry setup, coherent X-rays are diffracted by a grating G1 placed just after the object. The material and thickness of G1 are designed so as to apply a π/2 phase shift to the incident rays passing through it. This produces the Talbot effect, where self images of the grating are repeatedly generated at specific distances downstream. Phase shift Φ of the rays produced by the object distorts these self images. The analysis of a self image by fringe scanning [9] using the absorption grating G2 allows the retrieval of the beam deflection angle φ_u(θ, u,v) produced by the phase shift by the object, in the direction of u⃗, normal to the grating lines and the beam path :

φ_{u} (θ, u, v) = \frac{d}{2 π z_{T}} arg [\sum_{k = 1}^{S} I_{k} (u, v) exp (- \frac{2 π i k}{S})]

where I_k are the acquired images for the S steps of the fringe scanning process, i is the imaginary unit, d is the pitch of the gratings and z_T is the distance between the two gratings, given by

z_{T} = \frac{d^{2}}{2 λ}

with λ the X-ray wavelength. While other distances are possible, z_T is the distance for which the visibility of the fringe is maximum when using the π/2 phase grating [10]. The beam deflection angles φ_u constitute the measurements of the object used for the volume reconstruction. Because the beam deflection is proportional to the derivative of the phase shift, our measurements are called differential phase measurements :

φ_{u} (θ, u, v) = \frac{λ}{2 π} \frac{\partial Φ (u, v)}{\partial u}

Also, φ_u can be expressed with the integral of δ along the beam path, differentiated in the u⃗ direction (To simplify notations, the value of any multivariate function f at point (v ₁, v ₂, ⋯ , v_n) may be written indifferently as f (v ₁, v ₂, ⋯ , v_n) or f(v⃗) where v⃗ is the vector with same coordinates) :

φ_{u} (θ, u, v) = \frac{\partial}{\partial u} \int_{- \infty}^{\infty} δ (M_{θ}^{T} [\begin{matrix} u \\ v \\ w \end{matrix}]) d w

Taking the Fourier transform and after simplifying, we get the expression of the Fourier slice theorem (cf. for example [1, 11]) in our setup :

\tilde{φ_{u}} (θ, ω_{u}, ω_{v}) = 2 π i ω_{u} \tilde{δ} (M_{θ}^{T} [\begin{matrix} ω_{u} \\ ω_{v} \\ 0 \end{matrix}])

where

\tilde{φ_{u}}

and $\tilde{δ}$ are the Fourier transform of φ_u and δ respectively.

2.3. Backprojection

Inverting relation (6) with backprojection would allow the retrieval of the refractive index function δ. However, despite similarities with the Radon transform, the presence of rotation matrix M_θ and the derivative operation complicates the problem. Still, computing the backprojection of φ_u is needed because it is a basic step for reconstruction algorithms, and in particular leads to the filtered backprojection (FBP) method. In the literature, the reconstruction of an object slice by filtered backprojection of differential measurements has first been proposed in [12] in the case of optical tomography, showing how the usual ramp filter of classical tomography is replaced by a sign function filter. This technique has then been used in X-ray tomography by Talbot interferometry [10]. The other important related result in the literature concerns the transfer function of the circular tomosynthesis scan. It has been studied in the past [2] in the case of X-ray absorption measurements. The digital filtering of the backprojected measurements has consequently been performed in [3]. A framework for filtered backprojection reconstruction has then been proposed in [13]. In this section, we deal with the case of the laminography scan from differential measurements, which is a combination of these two previous problems. The frequency response and inverse filter are derived in our case.

To deal with the differential measurements, an integral φ(θ, u, v) of projections φ_u(θ, u, v) along u⃗ is considered. For each value of v, the unknown mean of the result is arbitrarily set to 0. φ(θ, u, v) is then proportional to the phase shift produced by the object along the given ray, plus an unknown additional term depending on v. In reciprocal space, the integral $\tilde{φ}$ of $\tilde{φ_{u}}$ is given by :

{\begin{array}{l} \tilde{φ} (θ, ω_{u}, ω_{v}) & = & \frac{1}{2 π i ω_{u}} \tilde{φ_{u}} (θ, ω_{u}, ω_{v}) if ω_{u} \neq 0 \\ \tilde{φ} (θ, 0, ω_{v}) & = & 0 \end{array}

The backprojection b(x,y,z) is defined here in the object’s local coordinate system as follows :

b (x, y, z) = \int_{0}^{2 π} φ (θ, u, v) d θ

= \int_{0}^{2 π} \int_{- \infty}^{+ \infty} \int_{- \infty}^{+ \infty} \tilde{φ} (θ, ω_{u}, ω_{v}) exp (2 π i (ω_{u} u + ω_{v} v)) d ω_{u} d ω_{v} d θ

with u and v obtained by relation (2). Using Eq. (8), we get the expression of the backprojection directly from the Fourier transform of the beam deflection measurements

\tilde{φ_{u}}

:

b (x, y, z) = \int_{0}^{2 π} \int_{- \infty}^{+ \infty} \int_{ℝ *} \frac{\tilde{φ_{u}} (θ, ω_{u}, ω_{v})}{2 π i ω_{u}} exp (2 π i (ω_{u} u + ω_{v} v)) d ω_{u} d ω_{v} d θ

Note that the inner integral is over ℝ^* = ℝ\{0}, to take into account the special case ω_u = 0. Using the Fourier slice theorem, Eq. (7), we obtain :

b (x, y, z) = \int_{0}^{2 π} \int_{- \infty}^{+ \infty} \int_{ℝ *} \tilde{δ} (M_{θ}^{T} [\begin{matrix} ω_{u} \\ ω_{v} \\ 0 \end{matrix}]) exp (2 π i (ω_{u} u + ω_{v} v)) d ω_{u} d ω_{v} d θ

From this last expression of the backprojection, changing the coordinates of the integral from ω_u, ω_v, θ to Cartesian coordinates ω_x ω_y ω_z in the object reciprocal space will give the relation between functions b and δ. The relation between old and new variables is :

[\begin{matrix} ω_{x} \\ ω_{y} \\ ω_{z} \end{matrix}] = M_{θ}^{T} [\begin{matrix} ω_{u} \\ ω_{v} \\ 0 \end{matrix}]

from which the following condition on ω_x, ω_y, ω_z for ω_u and ω_v to be defined is derived :

ω_{x}^{2} + ω_{z}^{2} > ω_{y}^{2} {tan}^{2} α

The strict nature of this inequality comes from the fact that ω_u = 0 is excluded from the integral in Eq. (12). The frequency region not satisfying condition (14) is not acquired by the measurements and is therefore lost in the reconstruction. This is a double cone with aperture 2α and axis ω_y, as illustrated in Fig. 2.

Fig. 2 Reciprocal space of the object, showing the unsampled region as a double cone of aperture 2α and axis ω_y. Two projection planes (yellow and green) are also displayed, for different angles θ ₁, θ ₂. Point A is one of the points at the intersection of the two planes, thus acquired two times during the scan.

Download Full Size | PDF

The backprojection written in the Cartesian reciprocal space finally reads :

b (x, y, z) = \underset{C}{∭} \tilde{δ} ([\begin{matrix} ω_{x} \\ ω_{y} \\ ω_{z} \end{matrix}]) \frac{2 exp (2 π i (ω_{x} x + ω_{y} y + ω_{z} z))}{cos α \sqrt{ω_{x}^{2} + ω_{z}^{2} - ω_{y}^{2} {tan}^{2} α}} d ω_{x} d ω_{y} d ω_{z}

where C is the subset of ℝ³ satisfying Eq. (14). This can then be interpreted as the filtering of the object function by H, which is the frequency response of the projection backprojection process, defined in the reciprocal space as :

{\begin{array}{l} H (ω_{x}, ω_{y}, ω_{z}) & = & \frac{2}{cos α \sqrt{ω_{x}^{2} + ω_{z}^{2} - ω_{y}^{2} {tan}^{2} α}} & if & ω_{x}^{2} + ω_{z}^{2} > ω_{y}^{2} {tan}^{2} α \\ H (ω_{x}, ω_{y}, ω_{z}) & = & 0 & else \end{array}

The obtained frequency response differs from the previous work in [3] and [13] only in the inequality condition, which is strict in the case of our differential measurements. This is however a crucial difference, because it shows how the null spatial frequency, corresponding to the mean value of the object function, is lost during the imaging process.

Similarly to previous work on laminography in the absorption case [13,14], the filtered back-projection algorithm can now be derived in the case of differential measurements, filtering projection data in order to cancel the frequency response. Noting that $ω_{u}^{2} = ω_{x}^{2} + ω_{z}^{2} - ω_{y}^{2} {tan}^{2} α$ , the frequency response H, expressed in the projection reciprocal space is the following function G, independent of θ and ω_v :

G (θ, ω_{u}, ω_{v}) = H (M_{θ} [\begin{matrix} ω_{u} \\ ω_{v} \\ 0 \end{matrix}]) = \frac{2}{cos α | ω_{u} |}

Multiplying the expression inside the integral of (11) by the inverse filter 1/G gives the FBP formula from the beam deflection measurements :

b_{fbp} (x, y, z) = \int_{0}^{2 π} \int_{- \infty}^{+ \infty} \int_{ℝ^{*}} \frac{cos α sgn (ω_{u})}{4 π i} \tilde{φ_{u}} (θ, ω_{u}, ω_{v}) exp (2 π i (ω_{u} u + ω_{v} v)) d ω_{u} d ω_{v} d θ

This means that the beam deflection data are multiplied by the following inverse filter F _inv in reciprocal space before backprojection :

F_{inv} (ω_{u}) = \frac{cos α sgn (ω_{u})}{4 π i}

3. Iterative reconstruction

Reconstruction for CL is usually performed by a FBP algorithm [4, 5], similar to the one presented for differential phase measurements in the previous section. However, because this reconstruction implicitly assumes that the unsampled double cone region in the reciprocal space is zero, it leads to typical artifacts in the form of directional blur in directions corresponding to missing frequencies. Attempts to reduce the effect of the missing spatial frequencies have been proposed in [13], by a further damping of low frequencies that are cropped by the double cone. However this assumes that low frequencies (and thus quantitative values) are not to be reconstructed. In [15], a FBP solution is also proposed, in which the filter is designed so that the blurred out-of-plane objects appear less structured, compared to the typically observed ring-shaped noise. A different solution to the problem of missing frequencies involves the introduction of prior information about the imaged object. However, the FBP algorithm cannot deal with such information, so that iterative estimation methods should be preferred.

3.1. Object constraints

Because CL aims at imaging flat extended objects, a prior knowledge on the general shape of the object comes naturally. A limited support of the object refractive index function δ is assumed in the vertical, y⃗ direction. This is interestingly in opposition to the classical limited support constraint used in CT (the limit case α = 0), which is defined on axes x⃗ and z⃗ so that the object is fully contained in the field of view. In our case, It is not required that the prior support perfectly matches the actual support of the sample. An approximative rectangular box shaped prior support is empirically obtained from the FBP reconstruction.

A second constraint that can be used is the limited range of the values taken by function δ. The expected values for δ should always be positive. An upper limit is also defined depending on the object. Non negativity and support constraints have been classically used in the field of optical microscope tomography [16]. Using these two classical constraints, it is proposed here to compensate for the missing frequencies of the CL scan, by estimating the object in agreement with both the measurements and priors.

3.2. Iterative filtered backprojection

The proposed reconstruction algorithm is based on iterative filtered backprojection (IFBP), with the addition of support and range constraints. The detector space is naturally discretized in square pixels whose size a × a depends on the equipment. Similarly, the object space is discretized in a finite grid of voxels with equivalent size a ³. Let s_k denote the estimate of the object at iteration k. It is a vector with dimension N, the number of voxels in the grid. b is a vector representing the measured beam deflection images φ_u. If there are M_p different projection angles and each image has M_d pixels, then b has dimension M = M_pM_d.

Updating the estimate s_k with the proposed algorithm is done in two consecutive steps as follows :

(\begin{array}{l} s_{k + 1}^{'} & \leftarrow & s_{k} + λ_{k} h_{k} & with & h_{k} = B (b - P s_{k}) \\ s_{k + 1} & \leftarrow & 𝒞 (s_{k + 1}^{'}) \end{array}

where P is the projection operator, modeling the imaging process leading to the measurement of a differential phase image, given an vector in object space. It is conceptually a matrix with dimension M × N. B is a general backprojection operator with dimension N × M, meaning that filtering is involved in the operator. 𝒞 is a nonlinear function which projects the object to the prior constraint. λ_k is a scalar factor whose value is chosen to guarantee convergence.

Matrices P and B are of course not explicitly expressed as this is infeasible with useable data sizes because of the limited memory. Instead, projections onto pixels and backprojections onto voxels are directly calculated. In the proposed CL reconstruction from differential phase measurements, projection and backprojection operators are defined as follows :

The operator P takes the integral of the input along each ray, then differentiates the result along axis u⃗.
The operator B performs a filtering of its input in detector space, whose action is to integrate along axis u⃗, compensate for the frequency response of backprojection, and damp high frequencies. The operator then backprojects values to the object space according to the geometry.

It has been demonstrated [8] how such an additional filtering step in the backprojection can improve the initial convergence speed of iterative algorithms. The filter that is used in this study is a combination of the FBP filter and a Hann window damping high frequencies. This damping is necessary to reduce reconstruction artifacts appearing after several iterations. Referring to the FBP filter, Eq. (19), the filter is defined as follows in the detector reciprocal space :

{\begin{array}{l} F (ω_{u}) & = & \frac{cos α sgn (ω_{u})}{4 π i} \cdot \frac{1}{2} (1 + cos (2 π ω_{u} a)) & if & ω_{u} \neq 0 \\ F (0) & = & 0 \end{array}

Function 𝒞 in Eq. (20) constrains its input to the prior information of compact support and limited range of values. More specifically, 𝒞 (s) is the closest vector to s satisfying the constraints. Let 𝒮 be the set of voxel indices in [0,N – 1], belonging to the prior support, and let δ _min and δ _max be the prior bounds for the values of function δ.

{\begin{array}{l} 𝒞 (s) [i] & = & s [i] & if δ_{min} < = s [i] < = δ_{max} and i \in 𝒮 \\ 𝒞 (s) [i] & = & δ_{min} & if s [i] < δ_{min} and i \in 𝒮 \\ 𝒞 (s) [i] & = & δ_{max} & if s [i] > δ_{max} and i \in 𝒮 \\ 𝒞 (s) [i] & = & 0 & else \end{array}

where s[i] is the i ^th element of vector s.

3.3. Filtering in detector space

To avoid numerical errors due to the discretization, the filtered backprojector works exclusively in real space. A straightforward implementation of filter F, Eq. (21), would be by computing the fast Fourier transform (FFT) of the input, then multiplying by the filter, and finally take the inverse Fourier transform. This is a suboptimal way of proceeding because it introduces aliasing when the number of samples is low. The explicit expression of the filter kernel f for F in real space is a more precise approach :

f (u) = \int_{- 1 / 2 a}^{1 / 2 a} F (ω_{u}) exp (2 π i ω_{u} u) d ω_{u}

where 1/2a is the Nyquist frequency. This can be simplified to the following function defined over the discrete values {na,n ∈ ℕ}:

{\begin{array}{l} f (n a) & = & \frac{n cos α}{4 π^{2} (n^{2} - 1) a} & for n even \\ f (n a) & = & \frac{cos α}{4 π^{2} n a} & for n odd \end{array}

The one-dimensional convolution of the backprojector’s input with kernel f, in direction u⃗, produces the filtered image that is then backprojected onto the voxels in object space.

3.4. Steepest descent

The IFBP step of the reconstruction, in the first step of (20), updates the current estimate s^k by moving along the direction h_k with a step size λ_k, whose value should be carefully computed to guarantee the convergence of the algorithm. Following the idea in [8], h_k can be seen as the steepest descent direction of a quadratic functional of s_k, so that the IFBP algorithms are minimizing this functional. By computing its integral with respect to s_k, h_k can be expressed from the gradient of a squared functional R :

h_{k} = - \frac{\partial}{\partial s_{k}} R (s_{k}) = - \frac{\partial}{\partial s_{k}} [\frac{1}{2} s_{k}^{T} B P s_{k} - s_{k}^{T} B b]

The choice for the step size λ_k in Eq. (20) is so that R(s_k ₊₁) is minimized in the direction h_k, making the IFBP step a steepest descent step. It is computed by analyzing the partial derivative of the squared functional with respect to λ_k :

\frac{\partial R (s_{k + 1})}{\partial λ_{k}} = \frac{\partial R (s_{k} + λ_{k} h_{k})}{\partial λ_{k}} = h_{k}^{T} B P s_{k} + λ_{k} h_{k}^{T} B P h_{k} - h_{k}^{T} B b

R is minimum when the derivative is zero, so that the optimal step size λ_k for steepest descent is obtained by :

λ_{k} = \frac{h_{k}^{T} h_{k}}{h_{k}^{T} B P h_{k}}

The second step of the algorithm projects s′_k ₊₁ onto the constraints with function 𝒞.

3.5. Implementation on GPU hardware

CL reconstruction using an iterative algorithm is a computationally intensive task, for which memory amount and calculation time must be considered. Unlike parallel-beam CT where each horizontal object slice has independent measurements, the full volume estimate is needed to compute the projections along inclined rays. The implementation of CL reconstruction thus requires a lot of computer memory and processing power for any usable data size. A solution using a Graphics Processing Unit (GPU) card has been developed, allowing the parallelization of projection and backprojection tasks. The card is a Tesla C2070 from NVIDIA, based on the CUDA architecture. It has enabled the reconstruction of volumes with a reasonable speed. For example, the reconstruction of 98 million voxels (700 × 200 × 700) from 500 projections with 1344 × 495 pixels is performed in 20 minutes in the current implementation, after 10 iterations.

Our implementation of the projection operator is ray-driven, which means that each ray (one per detector’s pixel) is processed in a separate thread. Conversely, the backprojection operator is voxel-driven, in that each voxel is processed by one thread. In this implementation choice, the two operators are then unmatched in the sense that the rays considered for projection and backprojection are slightly different. This however favors parallelism of the algorithm, and has been shown to reduce ring artifacts [17, 18].

3D texture memory is used to speed up the bilinear and trilinear interpolation during projection and backprojection. Because it is read-only memory, another writeable memory space is used, so that the full volume and projection space are actually stored twice in the graphics unit. This implementation decision has also influenced the choice of steepest descent as our optimization method, compared to a conjugate gradient algorithm that would have been more memory-consuming.

3.6. Cone beam geometry

As stated in section 2.1, the proposed imaging method assumes parallel beam illumination, as delivered by synchrotron facilities. Although not detailed in this paper, we think that the extension to cone beam illumination is feasible with the following modifications. Concerning the experimental setup, the grating interferometer is readily applicable to cone beam illumination, by adapting the pitch of the gratings according to the magnification ratio, as presented by [19]. The volume reconstruction is more challenging because the Fourier slice theorem is not valid in cone beam geometry. In addition to the geometric differences due to the perspective transform of the cone beam setup, it is also necessary to adapt the filter applied to projection data. This filter must compensate for the frequency response of the cone beam projection-backprojection. An approach weighting the projection data, similar to the FDK [20] method should be necessary.

4. Experiment

Experiments of X-ray computed laminography using a Talbot interferometer have been performed at the beamline BL20XU at SPring-8 synchrotron facility in Hyogo, Japan. The interferometer has been set up optimally for the selected wavelength λ = 0.7A. Å 5.3 μm pitch gold grating was used for G1, with a nominal thickness of 1.9 μm, so that it produces a π/2 phase shift at λ. The absorption gold grating G2 had the same pitch, with a higher thickness of 30 μm. It was placed at a distance z_T of 201 mm from G1. The X-ray detector was composed of a 10 μm P43 phosphor screen (Gd₂O₂S:Tb⁺ powder) coupled to a cooled CCD camera (Hamamatsu Photonics C4742-98-24A), with an effective pixel size a = 4.34μm. The number of pixels in acquired images was set to 1344 × 495.

The sample was a polypropylene mesh with a wire diameter of 215 μm and opening size of 250 μm (cf. Fig. 3(a)). The beam width was around 5 mm, so that the field of view was much smaller than the sample. This and the small thickness of the mesh are properties that make CL a good candidate for imaging. Moreover, the structures of the sample are well understood, cover the entire field of view and have sharp edges to facilitate the measurement of spatial resolution. The mesh is supported by a 0.5 mm thick acrylic plate, and placed on a cylindrical holder in acrylic, as pictured in Fig. 3(b). The diameter of the holder was 40 mm and the effect of its wall across the X-ray beam was assumed to be negligible. The angle α between the rotation axis and the beam normal plane was set to 20°. Smaller angles α were found to cause strong phase wrapping in the differential phase images at some angles θ. 20° was a good compromise between the quality of the measurements and the need for a small angle to reduce the region of missing frequencies. The object was rotated for 360°, with an angle step of 0.72°, so that projections along 500 angles were acquired. For each angle, S = 5 steps of fringe scanning were performed, by successively acquiring a transmission image then translating the absorption grating by $\frac{d}{S} \vec{u}$ . The exposure time for one transmission image was 0.4 s.

Fig. 3 (a) An optical microscope image of the polypropylene mesh, showing the weave pattern. The wire diameter is 215 μm and the opening is 250 μm. (b) Setup of the mesh sample during CL scan. The X-ray beam transmits through the mesh, acrylic plate support and cylindrical acrylic holder.

Download Full Size | PDF

An example of differential phase image obtained from fringe scanning is shown in Fig. 4, computed as in Eq. (3). The differential direction u⃗ is horizontal, from right to left. A comparison between the FBP reconstruction and the iterative method is presented in Fig. 5. The FBP algorithm, as described by Eq. (18), has the advantage of being a fast reconstruction method, requiring only one filtering step and one backprojection step. It gives an reliable first estimation of the object. However its implicit assumption that unsampled spatial frequencies are zero produces strong artifacts. Most notably, strong linear blur artifacts appears in directions with angle α from the vertical axis, as can be seen in the vertical slice. Moreover, the unsampled mean value of the refractive index decrement δ is also assumed to be zero, producing unrealistic negative values. The introduction of range and support constraints in the iterative algorithm helps reducing the drawbacks of FBP, at the price of a more computationally expensive method. The reduction of blur stripes can be observed, and values for δ have been constrained in the range [0 10⁻⁶]. The presented result was obtained after a fixed number of 10 iterations. In Fig. 6, line profiles in horizontal and vertical directions passing through the mesh wire are presented. The improvement of the iterative reconstruction over the filtered backprojection is mainly seen in the vertical direction, where the low spatial frequencies are better estimated. From such profiles it is also possible to get a measure of the spatial resolution. Assuming that the mesh wire has a constant δ, line profiles actually show the response of the imaging system to sharp edges. A simple measurement of the resolution is the 10% to 90% distance. In the case of the 215 μm wide structures of the presented mesh sample, it has been estimated to 30 μm in the horizontal direction and 40 μm in the vertical direction.

Fig. 4 Beam deflection angle image of the mesh sample obtained by fringe scanning

Download Full Size | PDF

Fig. 5 Comparison of FBP and iterative reconstruction methods. (a) Horizontal slice normal to y⃗ obtained by FBP. (b) Same slice obtained by iterative reconstruction. (c) Vertical slice normal to z⃗ obtained by FBP. (d) Same vertical slice obtained by iterative reconstruction, showing fewer artifacts. Red dotted lines show the position of horizontal h and vertical v profiles in Fig. 6. The vertical extent of image (d) matches the prior support of the object.

Download Full Size | PDF

Fig. 6 Estimated refractive index decrement δ along the horizontal and vertical profiles defined in Fig. 5(d). The iterative reconstruction (red solid line) has no negative values and square shaped profiles as we can expect from the mesh sample. The filtered backprojection (black dotted line) shows unrealistic negative values as well as strong distortion of the ideal square shape in the vertical direction.

Download Full Size | PDF

The effectiveness of the Hann filtering before backprojection is illustrated in Fig. 7, showing how it prevents high frequency artefacts from appearing after several iterations. A volume render of the estimated refractive index decrement δ is finally presented in Fig. 8, showing details in the weave pattern.

Fig. 7 Magnified region of a reconstructed horizontal slice after 30 iterations, with and without Hann filtering before backprojection. (a) Without Hann filtering, undesirable high frequency lines are present in the image. (b) With Hann filtering, these artefacts are successfully reduced.

Download Full Size | PDF

Fig. 8 Volume rendering of the polypropylene mesh sample, reconstructed by the iterative method. A video is also available (Media 1).

Download Full Size | PDF

5. Conclusion

The imaging of flat extended objects by X-ray phase computed laminography using a Talbot interferometer has been presented. After deriving the filtered backprojection algorithm in this particular setup, an iterative reconstruction method has been proposed to deal with the typical directional artifacts. The iterative scheme is based on the iterative filtered backprojection algorithm, with a Hann windowed sign function filter, and additional constraints on the range of values and support of the estimated object. The obtained solution helps improving the reconstruction of refractive index decrement δ through a better estimation of the unsampled frequencies. The proposed method has been demonstrated by the imaging of a polypropylene mesh sample at a synchrotron radiation facility. By an implementation of the algorithm on GPU hardware, it is now possible to obtain the reconstruction result in a few minutes. Absorption or phase laminographic imaging has already found applications in several fields, allowing for example the inspection of microsystem devices [14], paintings [21] or damage in polymer composites [5]. We believe that the proposed constrained iterative reconstruction can improve results in such applications compared to the classical filtered backprojection. Another advantage is that using the same X-ray Talbot interferometry setup, the image of visibility reduction [22] can be obtained and used for the estimation of properties in microstructures producing small-angle scattering. Finally, although we have assumed in this paper a parallel illumination from synchrotron facilities, the extension to the cone beam geometry of conventional laboratory sources is certainly feasible, and will be dealt with in future work.

Acknowledgments

The synchrotron radiation experiment was performed with the approval of the SPring-8 advisory committee (proposal 2010A1164). This work was supported by the SENTAN project of the Japan Science and Technology Agency (JST).

References and links

1. F. Natterer, The Mathematics of Computerized Tomography (Society for Industrial and Applied Mathematics, 2001). [CrossRef]

2. D. G. Grant, “Tomosynthesis: a three-dimensional radiographic imaging technique,” IEEE Trans. Biomed. Eng. 19(1), 20–28 (1972). [CrossRef] [PubMed]

3. H. Matsuo, A. Iwata, I. Horiba, and N. Suzumura, “Three-dimensional image reconstruction by digital tomosynthesis using inverse filtering,” IEEE Trans. Med. Imaging 12(2), 307–313 (1993). [CrossRef] [PubMed]

4. L. Helfen, T. Baumbach, P. Cloetens, and J. Baruchel, “Phase-contrast and holographic computed laminography,” Appl. Phys. Lett. 94(10), 104103 (2009). [CrossRef]

5. F. Xu, L. Helfen, A. J. Moffat, G. Johnson, I. Sinclair, and T. Baumbach, “Synchrotron radiation computed laminography for polymer composite failure studies,” J. Synchrotron Radiat. 17(2), 222–226 (2010). [CrossRef] [PubMed]

6. A. Momose, S. Kawamoto, I. Koyama, Y. Hamaishi, K. Takai, and Y. Suzuki, “Demonstration of X-ray Talbot interferometry,” Jpn. J. Appl. Phys. 42(7B), 866–868 (2003). [CrossRef]

7. E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” Phys. Med. Biol. 53, 4777 (2008). [CrossRef] [PubMed]

8. D. Lalush and B. Tsui, “Improving the convergence of iterative filtered backprojection algorithms,” Med. Phys. 21(8), 1283–1286 (1994). [CrossRef] [PubMed]

9. J. H. Bruning, D. R. Herriott, J. E. Gallagher, D. P. Rosenfeld, A. D. White, and D. J. Brangaccio, “Digital wavefront measuring interferometer for testing optical surfaces and lenses,” Appl. Opt. 13(11), 2693–2703 (1974). [CrossRef] [PubMed]

10. A. Momose, W. Yashiro, Y. Takeda, Y. Suzuki, and T. Hattori, “Phase tomography by X-ray Talbot interferometry for biological imaging,” Jpn. J. Appl. Phys. 45, 5254–5262 (2006). [CrossRef]

11. J. Als-Nielsen and D. McMorrow, Elements of Modern X-ray Physics, 2nd ed. (John Wiley and Sons, 2011). [CrossRef]

12. G. W. Faris and R. L. Byer, “Three-dimensional beam-deflection optical tomography of a supersonic jet,” Appl. Opt. 27(24), 5202–5212 (1988). [CrossRef] [PubMed]

13. G. Lauritsch and W. H. Härer, “A theoretical framework for filtered backprojection in tomosynthesis,” Proc. SPIE 3338, 1127–1137 (1998). [CrossRef]

14. L. Helfen, A. Myagotin, P. Mikulík, P. Pernot, A. Voropaev, M. Elyyan, M. Di Michiel, J. Baruchel, and T. Baumbach, “On the implementation of computed laminography using synchrotron radiation,” Rev. Sci. Instrum. 82, 063702 (2011). [CrossRef] [PubMed]

15. G. M. Stevens, R. Fahrig, and N. J. Pelc, “Filtered backprojection for modifying the impulse response of circular tomosynthesis,” Med. Phys. 28, 372–380 (2001). [CrossRef] [PubMed]

16. O. Nakamura, S. Kawata, and S. Minami, “Optical microscopy tomography. II. Nonnegative constraint by a gradient-projection method,” J. Opt. Soc. Am. A 5(4), 554–561 (1988). [CrossRef]

17. G. L. Zeng and G. T. Gullberg, “Unmatched projector/backprojector pairs in an iterative reconstruction algorithm,” IEEE Trans. Med. Imaging 19, 548–555 (2000). [CrossRef] [PubMed]

18. R. Guedouar and B. Zarrad, “A comparative study between matched and mis-matched projection/back projection pairs used with ASIRT reconstruction method,” Nucl. Instrum. Methods Phys. Res. A 619(1–3), 225–229 (2010). [CrossRef]

19. M. Engelhardt, J. Baumann, M. Schuster, C. Kottler, F. Pfeiffer, O. Bunk, and C. David, “High-resolution differential phase contrast imaging using a magnifying projection geometry with a microfocus x-ray source,” Appl. Phys. Lett. 90, 224101 (2007). [CrossRef]

20. L. A. Feldkamp, L. C. Davis, and J. W. Kress, “Practical cone-beam algorithm,” J. Opt. Soc. Am. A 1, 612–619 (1984). [CrossRef]

21. K. Krug, L. Porra, P. Coan, A. Wallert, J. Dik, A. Coerdt, A. Bravin, M. Elyyan, P. Reischig, L. Helfen, and T. Baumbach, “Relics in medieval altarpieces? Combining X-ray tomographic, laminographic and phase-contrast imaging to visualize thin organic objects in paintings,” J. Synchrotron Radiat. 15, 55–61 (2008). [CrossRef]

22. F. Pfeiffer, M. Bech, O. Bunk, P. Kraft, E. F. Eikenberry, Ch. Brönnimann, C. Grünzweig, and C. David, “Hard-X-ray dark-field imaging using a grating interferometer,” Nat. Mater. 7, 134–137 (2008). [CrossRef] [PubMed]

Iterative reconstruction in x-ray computed laminography from differential phase measurements

Abstract

1. Introduction

2. X-ray computed laminography by Talbot interferometry

2.1. Scan geometry

2.2. Differential phase measurements

2.3. Backprojection

3. Iterative reconstruction

3.1. Object constraints

3.2. Iterative filtered backprojection

3.3. Filtering in detector space

3.4. Steepest descent

3.5. Implementation on GPU hardware

3.6. Cone beam geometry

4. Experiment

5. Conclusion

Acknowledgments

References and links

Supplementary Material (1)

Cited By

Figures (8)

Equations (27)

Optics Express