OSA's Digital Library

Virtual Journal for Biomedical Optics

Virtual Journal for Biomedical Optics

| EXPLORING THE INTERFACE OF LIGHT AND BIOMEDICINE

  • Editors: Andrew Dunn and Anthony Durkin
  • Vol. 6, Iss. 9 — Oct. 3, 2011
« Show journal navigation

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

Sébastien Harasse, Wataru Yashiro, and Atsushi Momose  »View Author Affiliations


Optics Express, Vol. 19, Issue 17, pp. 16560-16573 (2011)
http://dx.doi.org/10.1364/OE.19.016560


View Full Text Article

Acrobat PDF (955 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Phase-contrast X-ray computed laminography is demonstrated for the volume reconstruction of extended flat objects, not suitable to the usual tomographic scan. Using a Talbot interferometer, differential phase measurements are obtained and used to reconstruct the real part of the complex refractive index. The specific geometry of laminography leads to unsampled frequencies in a double cone in the reciprocal space, which degrades the spatial resolution in the direction normal to the object plane. First, the filtered backprojection formula from differential measurements is derived. Then, reconstruction is improved by the use of prior information of compact support and limited range, included in an iterative filtered backprojection algorithm. An implementation on GPU hardware was required to handle the reconstruction of volumes within a reasonable time. A synchrotron radiation experiment on polymer meshes is reported and results of the iterative reconstruction are compared with the simpler filtered backprojection.

© 2011 OSA

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

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

] 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

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

]). Although the frequency response of such scans was well understood, the increasing usage of computers has later allowed digital filtering [3

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]

] 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

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

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]

], using the in-line phase contrast retrieval method.

It is proposed here to use an X-ray Talbot interferometer [6

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]

] 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 – δ + . 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

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]

]). The proposed reconstruction algorithm is based on iterative filtered backprojection [8

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

] (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.

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θ=[cosθ0sinθsinαsinθcosαsinαcosθcosαsinθsinαcosαcosθ]
(1)
A vector r⃗ = [x,y,z]T expressed in object coordinates is transformed into detector coordinates q⃗ = [u,v,w]T as follows :
q=Mθr
(2)

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

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]

] 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)=d2πzTarg[k=1SIk(u,v)exp(2πikS)]
(3)
where Ik 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 zT is the distance between the two gratings, given by
zT=d22λ
(4)
with λ the X-ray wavelength. While other distances are possible, zT is the distance for which the visibility of the fringe is maximum when using the π/2 phase grating [10

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]

]. 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)=λ2πΦ(u,v)u
(5)
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 1, v 2, ⋯ , vn) may be written indifferently as f (v 1, v 2, ⋯ , vn) or f(v⃗) where v⃗ is the vector with same coordinates) :
φu(θ,u,v)=uδ(MθT[uvw])dw
(6)

Taking the Fourier transform and after simplifying, we get the expression of the Fourier slice theorem (cf. for example [1

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

, 11

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

]) in our setup :
φu˜(θ,ωu,ωv)=2πiωuδ˜(MθT[ωuωv0])
(7)
where φu˜ and δ˜ are the Fourier transform of φu and δ respectively.

2.3. Backprojection

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 φ˜ of φu˜ is given by :
{φ˜(θ,ωu,ωv)=12πiωuφu˜(θ,ωu,ωv)ifωu0φ˜(θ,0,ωv)=0
(8)

The backprojection b(x,y,z) is defined here in the object’s local coordinate system as follows :
b(x,y,z)=02πφ(θ,u,v)dθ
(9)
=02π++φ˜(θ,ωu,ωv)exp(2πi(ωuu+ωvv))dωudωvdθ
(10)
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 φu˜:
b(x,y,z)=02π+*φu˜(θ,ωu,ωv)2πiωuexp(2πi(ωuu+ωvv))dωudωvdθ
(11)
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)=02π+*δ˜(MθT[ωuωv0])exp(2πi(ωuu+ωvv))dωudωvdθ
(12)

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 :
[ωxωyωz]=MθT[ωuωv0]
(13)
from which the following condition on ωx, ωy, ωz for ωu and ωv to be defined is derived :
ωx2+ωz2>ωy2tan2α
(14)
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 θ 1, θ 2. Point A is one of the points at the intersection of the two planes, thus acquired two times during the scan.

3. Iterative reconstruction

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

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]

]. 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 3. Let sk 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 Mp different projection angles and each image has Md pixels, then b has dimension M = MpMd.

Updating the estimate sk with the proposed algorithm is done in two consecutive steps as follows :
(sk+1sk+λkhkwithhk=B(bPsk)sk+1𝒞(sk+1)
(20)
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

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

] 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 :
{F(ωu)=cosαsgn(ωu)4πi12(1+cos(2πωua))ifωu0F(0)=0
(21)

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 δ.
{𝒞(s)[i]=s[i]ifδmin<=s[i]<=δmaxandi𝒮𝒞(s)[i]=δminifs[i]<δminandi𝒮𝒞(s)[i]=δmaxifs[i]>δmaxandi𝒮𝒞(s)[i]=0else
(22)
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)=1/2a1/2aF(ωu)exp(2πiωuu)dωu
(23)
where 1/2a is the Nyquist frequency. This can be simplified to the following function defined over the discrete values {na,n ∈ ℕ}:
{f(na)=ncosα4π2(n21)afornevenf(na)=cosα4π2nafornodd
(24)

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 sk by moving along the direction hk with a step size λk, whose value should be carefully computed to guarantee the convergence of the algorithm. Following the idea in [8

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

], hk can be seen as the steepest descent direction of a quadratic functional of sk, so that the IFBP algorithms are minimizing this functional. By computing its integral with respect to sk, hk can be expressed from the gradient of a squared functional R :
hk=skR(sk)=sk[12skTBPskskTBb]
(25)

The choice for the step size λk in Eq. (20) is so that R(sk +1) is minimized in the direction hk, making the IFBP step a steepest descent step. It is computed by analyzing the partial derivative of the squared functional with respect to λk :
R(sk+1)λk=R(sk+λkhk)λk=hkTBPsk+λkhkTBPhkhkTBb
(26)
R is minimum when the derivative is zero, so that the optimal step size λk for steepest descent is obtained by :
λk=hkThkhkTBPhk
(27)

The second step of the algorithm projects sk +1 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

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

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]

].

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

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]

]. 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

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]

] 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 zT of 201 mm from G1. The X-ray detector was composed of a 10 μm P43 phosphor screen (Gd2O2S: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 dSu. 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.

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−6]. 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
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.
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.

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.
Fig. 8 Volume rendering of the polypropylene mesh sample, reconstructed by the iterative method. A video is also available (Media 1).

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

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]

], paintings [21

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]

] or damage in polymer composites [5

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]

]. 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

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]

] 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]

OCIS Codes
(110.6760) Imaging systems : Talbot and self-imaging effects
(110.6955) Imaging systems : Tomographic imaging

ToC Category:
Imaging Systems

History
Original Manuscript: June 16, 2011
Revised Manuscript: August 3, 2011
Manuscript Accepted: August 4, 2011
Published: August 12, 2011

Virtual Issues
Vol. 6, Iss. 9 Virtual Journal for Biomedical Optics

Citation
Sébastien Harasse, Wataru Yashiro, and Atsushi Momose, "Iterative reconstruction in x-ray computed laminography from differential phase measurements," Opt. Express 19, 16560-16573 (2011)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-19-17-16560


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  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]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Supplementary Material


» Media 1: AVI (3664 KB)     

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited