Open Access
Add to BibSonomyAbstract
In this paper, a model of the beam propagation is developed according to the physical properties of THz waves used in THz computed tomography (CT) scan imaging. This model is first included in an acquisition simulator to observe and estimate the impact of the Gaussian beam intensity profile on the projection sets. Second, the model is introduced in several inversion methods as a convolution filter to perform efficient tomographic reconstructions of simulated and real acquired objects. Results obtained with three reconstruction methods (BFP, SART and OSEM) are compared to the techniques proposed in this paper. We will demonstrate an increase of the overall quality and accuracy of the 3D reconstructions.
© 2011 Optical Society of America
1. Introduction
Terahertz (THz) technology, and especially THz spectro-imaging, is now a well-established tool in the field of contact-free and non-destructive testing (NDT). While commercial devices become available, the fundamental issues of the imaging process [1] often remain overlooked. In the field of 3D imaging, X-Ray Computed Tomography (CT) is an omnipresent technique which provides cross-sectional images of an object by analyzing the radiation transmitted by the sample through different incidence angles. This technique provides 3D visualization of dense materials such as human body, biological tissues and also starts to enter in the industrial field. Nevertheless, this powerful technique cannot be easily applied to soft materials such as plastics, papers or paintings owing to the low absorption of the X-Ray radiation for example. Alternatively, THz radiation proposes attractive features such as non-invasive and non-destructive analysis, transparency and good penetration depth through various materials specially below 1 THz. All these remarkable properties make THz radiation very efficient for direct applications in non-destructive inspection and the security field. Most of the experimental demonstrations were performed in 2D whereas THz CT, which is a powerful technique, the literature remains scattered [1–6].
The first reason concerns the diffraction effects and Fresnel losses experienced by the propagation of the THz wave through the sample [7]. Secondly, although low NEP ( at room temperature) - and consequently fast detectors [8] - have been developed and successfully applied in imaging [9], they are not yet integrated in an array configuration. Detectors are thus mono- to few-pixels large, resulting in a long scanning process in order to acquire 2D images. Reconstruction algorithms allowing accurate reconstruction from a limited number of slices are crucial in order to circumvent these current technical limitations, but also for future development. The limitation is also directly connected to the number of projection data. At least, an important consideration in CT concerns the choice of the reconstruction method to be able to visualize the different cross-sectional images and the final 3D volume of the sample. In X-Ray, iterative reconstruction methods have been proposed for X-Ray CT such as the Simultaneous Algebraic Reconstruction Technique (SART) [10] and the Ordered Subsets Expectation Maximization (OSEM) [11,12]. These well-established techniques inherited from X-Ray know-how were successfully applied in a first approximation to THz CT.
In this paper, we propose to improve the physical modeling by adding an heuristic propagation of a gaussian beam. We aim to take into account a more realistic physical behavior of the electromagnetic waves used in THz CT. This approach is justified by the fact that the reconstruction methods assume a ray tracing approach while the THz beam profile is very far from this basic hypothesis. There is indeed a trade-off range in the focusing condition of the diffracting THz wave. While 2λ-waist beams can be qualitatively treated as a ray of light, tightly focused ones (down the ultimate resolution of λ/2) can badly accommodate with the ray tracing approach of all the reconstruction methods used in X-Ray. This model is first included in an acquisition simulator to observe and estimate the impact of the gaussian beam distribution on the acquisition images. Second, the model is tested in several reconstruction methods as a convolution/deconvolution filter to perform a more efficient tomographic reconstruction of simulated and real acquired objects. Results obtained with three usual methods (BFP, SART and OSEM) are compared to the new methods proposed in the paper. The discussion focuses on the efficiency of optimized algorithms to increase the overall quality and accuracy of the reconstructions.
2. THz computed tomography
2.1. Experimental setup properties
The experimental setup of the 3D millimeter wave tomographic scanner is shown in Fig. 1(a). The output beam of a compact millimeter wave (mmw) Gunn diode coupled with a horn antenna is collimated using an off-axis parabolic mirror. The Gunn diode used is a 80 GHz diode coupled to a frequency tripler delivering 0.3 mW at 240 GHz (wavelength 1.25 mm). The THz beam is then focused with a Teflon® lens (f = 75 mm focal length and D = 50.8 mm diameter) on the sample which is positioned on a three-axes motorized stage comprising the X,Y and θ movements, respectively. Figure 1(b) shows the 2D transversal profile of the THz beam at the beam waist visualized using a photothermal THz convertor. This result shows that, at the sample position, the beam profile is homogeneous with a Gaussian circular shape (2 mm beam diameter, measured at FWHM) in agreement with the theoretical values obtained from the propagation of Gaussian beam models. This result indicates that the spatial resolution of the 3D millimeter wave tomographic scanner is limited to a few millimeters owing to the long wavelength of the emitting sources.
Fig. 1 (a) Experimental setup. L1 - L2 : HDPE lenses, PM1 - PM2 : Parabolic Mirrors. (b) 2D spatial profile of the THz beam waist at the sample position (240 GHz source) visualized with a pyroelectric detector.
In order to perform a 3D reconstruction of volumetric objects, the experimental procedure is organized as follows [13]. First, we record a 2D transmission image of the sample in the X and Y directions (1 × 1 mm2). This is repeated for Nθ = 18 projections from θ = 0 to 180°. From these projection data, we are able to construct the sinograms of the object which represent, for a given horizontal slice, the evolution of the transmitted THz intensity as a function of the rotation angle. Finally, we apply tomographic algorithms to reconstruct the final 3D volume of the sample.
2.2. Usual tomographic reconstruction methods
Tomography is a volume reconstruction process widely used in X-Ray scan imaging. The volume imaging acquired object is obtained from a set of projections measured from the exterior of the objet [14, 15]. Each projection value, denoted ℛθ(ρ), represents the absorption sum undergone by the ray along the line (θ,ρ). It is modeled by the Radon transform [16] as follows :
where θ and ρ are respectively the angular and radial coordinates of the projection line (θ,ρ), δ is the Dirac impulse and f(x,y) is the absorption function. The acquisition along several angles of an horizontal cross-section represents a sinogram composed with Nθ lines, corresponding to the number of projections, and Nρ columns, corresponding to the number of acquired values per projection. For instance, the sinogram Fig. 2(b) corresponds to an acquisition, at the indicated cross-section, of the object represented Fig. 2(a) with Nθ = 18 and Nρ = 46.Fig. 2 (a) Medicine box acquired at the indicated cross-section. (b) Corresponding sinogram with Nθ = 18 and Nρ = 46. (c) BFP, (d) SART, (e) OSEM results.
The inverse Radon transform, denoted backprojection of filtered projections (BFP) [14, 17], is usually performed to recover the original domain from the projections. First, this inversion filters each projection in the Fourier domain to increase geometric details. Second, it computes each f(x,y) value as a sum of filtered projections crossing (x,y).
Apart from this direct method, iterative techniques originally developed for X-Ray scan imaging, have been introduced for THz CT imaging in [17]. For instance, the simultaneous algebraic reconstruction technique (SART) [10, 18] approaches the solution of the linear equation system I = AT R, where I is the image, R is the sinogram and A is the weight-matrix using an iterative process following k ∈ [0···Niter[. Each sub-iteration s, 0 ≤ s < Nθ, updates each pixel of the image Ik,s by comparing the original projection ℛθs with (computed from Ik,s−1) as follows :
where :A super-iteration k is completed when all the projections have been used. Iterations in k are performed until the convergence of the solution. The initial I0,0 image is usually an uniform image. The scheme on Fig. 3 summarizes the iterative method process. The update step of the algorithm depends on the iterative technique used. In the following, we consider the SART and the OSEM (ordered subset expectation maximization) [12] techniques, which only differ on the number of projection subsets and the error update estimation.
The images on Figs. 2(c)–2(e) represent the reconstructions from the sinogram on Fig. 2(b) using the three methods. In our example, these results reveal a part of a tablet inside the medicine box.
3. From propagation beam modeling to tomographic reconstruction optimizations
3.1. Propagation beam observation and modeling
In X-ray CT transmission processes, the beam shape can be considered constant because of the smallness of the wavelength with respect to the object size so that there is no diffraction. For the same reason, the intensity distribution is uniform over the beam cross-section. This property allows Radon model to be used directly in all reconstruction methods.
In THz CT imaging, the propagation beam is close to a Gaussian distribution which is both determined by the THz wave properties and, eventually, the lens used to enforce the beam focus. The radius of the beam (from the beam axis) has its minimum value w0 at the so-called beam waist. According to the wavelength λ, the radius at the position z from the beam waist is :
where is the Rayleigh range. Moreover, the intensity distribution over the cross-section is given by : where :As an illustration, Fig. 4(a) shows beam cross-sections obtained with a 240 GHz source. These 2D images are obtained using a photothermal THz sensor positioned along the Z-axis with a displacement of 5 mm between each orthogonal image profile. Note that the energy decreases from the center to the edge of the beam following a Gaussian distribution. A profile at the central plane of the beam is given in Fig. 4(b). It shows the variation of the radius (horizontally) and the Gaussian distribution of the intensity (vertically).
Fig. 4 (a) Propagation beam acquired along the Z-axis with a 5 mm step. (b) Beam propagation along the Z-axis.
Using the Eq. (3) and Eq. (4), the observed Gaussian beam can be computed to simulate a given source. For instance, Fig. 5 shows the 3D profile and energy distribution of a beam simulating the 240 GHz propagation. When the beam is focalized by a lens, THz energy profile is influenced by the lens which the energy distribution is treated by the Gaussian model. If the propagation beam distribution is unknown, the model could be estimated from each observation of the intensity profiles given in Fig. 4(a) using, for example, a non trivial 3D Newton-Gauss algorithm.
Fig. 5 Simulated propagation of the beam for a 240 GHz source (Intensity (a.u), position in meters), FWHM = 2 mm at beam waist according to previously measured beam.
3.2. Simulated acquisition
Using the Gaussian beam modeling, we develop now the simulated acquisition from such a beam. Figure 6 illustrates the θ-projection acquisition. An object is modeled in a N2 pixel image and the Gaussian beam is given at the same size and resolution. The object is rotated by θ from the original angular position to perform easily the column-by-column convolutions with the Gaussian beam. These convolutions simulate at each depth of the projection point of view (POV), the manner that the beam acquires the object. Then, they lead to a modified object corresponding to what it is viewed by the sensors. Finally, the acquisition of N samples leads to the projection acquired through the Gaussian beam. The process is easy to extend in 3D : each 1D convolution is replaced by a 2D convolution between each depth-image and the corresponding 2D cross-section of the beam.
Mathematically, this process is obtained by the following optimized version of the Eq. (1) :
where I(x,y) is the Gaussian propagation model Eq. (4), * is the convolution operator and (x0, y0) is the center of the beam waist. The f(x, y) domain is rotated by θ to assume ϕ = 0. It leads to an horizontal acquisition allowing an easier computation of the column-by-column convolutions. Note that the integral along y is simplified by fixing the value y = ρ.As an illustration, let us consider the image Fig. 7(b) representing the four metallic bars of the object in Fig. 7(a). An usual Radon acquisition is shown in Fig. 8(a). Simulated acquisition using Gaussian beam convolutions is given in Fig. 8(b). In such an acquisition, the edges are not well defined anymore. Then, the simulated sinogram characteristics are closer to the real acquisition given in Fig. 8(c). It puts in evidence the better reliability of our model to simulate the THz physical acquisition phenomenon and moreover, the discrepancy between Gaussian beam and ray-tracing approaches.
Fig. 7 (a) Object composed of four metallic bars (two bars diameter 10mm on the top and on the left, one bar diameter 12mm on the right and one bar diameter 8mm on the bottom), (b) Synthetic model of the object.
Fig. 8 (a) Usual Radon acquisition from synthetic model, (b) Simulated acquisition using Gaussian beam convolutions of the synthetic model, (c) Acquisition of the real object. Grey scale indicates absorption rate (full absorption in white, and background air absorption tends to 0).
3.3. Tomographic reconstruction optimizations
Up to here, the modeling and the influence of the Gaussian beam have been introduced. Now, we deal with its consideration in the reconstruction methods in order to try to reduce induced artifacts.
Using the BFP, the idea consists in inverting the convolution proposed in Fig. 5 as described on the scheme Fig. 9. The projection acquired through a Gaussian beam, denoted , is retro-projected into a temporary image sized N2. The image is deconvolued column-by-column with the beam model to obtain Iθ. This one is rotated by the angle −θ and added to the result image. The Wiener deconvolution method [19] is used to avoid well-known problems of direct deconvolution in Fourier space.
Using iterative methods, the projection calculated through a simulated Gaussian beam at an iteration k is given by the discrete version of the Eq. (5). Then, the update step comparison between computed and measured acquisitions becomes consistent because the projections are both obtained through a simulated/real Gaussian beam. Moreover, similarly than BFP, the obtained error projection is deconvolued before it is added in the Ik+1 image. Everything else in the algorithm Fig. 3 remains the same.
4. Results and discussion
4.1. Discussion from simulated acquisition
Images in Figs. 10(a)–10(c) correspond to the results obtained with BFP, SART and OSEM from the theoretical sinogram Fig. 8(a). Figures 10(d)–10(f) represent the images obtained with the usual methods from the simulated sinogram Fig. 8(b). The distortions of each metallic bar circular shape are a consequence of the Gaussian distribution encountered during the acquisition. Inversely, images on Figs. 10(g)–10(i), which are obtained from the simulated acquisition using the optimized methods, provide a closer circular shape of the metallic bars.
Fig. 10 (a)(b)(c) BFP, SART, OSEM results from the ideal sinogram Fig. 8(a). (d)(e)(f), BFP, SART, OSEM from simulated Gaussian beam acquisition 8(b). (g)(h)(i) BFP, SART and OSEM reconstructions from the same sinogram using optimized algorithms.
To show the new method efficiency, quality and accuracy of the results are now compared to each other using image-based comparison metrics. First, the quality preservation is estimated using the equivalence rate between two images. The image properties are compared to each other using the Structural SIMilarity (SSIM) criterion [20] :
where 0 ≤ l(I,J) ≤ 1 (resp. 0 ≤ c(I,J) ≤ 1) is the global intensity (resp. contrast) comparison between two images I and J, and 0 ≤ r(I,J) ≤ 1 is the correlation coefficient. 0 ≤ SSIM(I,J) ≤ 1 gives the quality rate of J compared to the reference I. The closer the SSIM value is to 1, the better the reconstruction quality is.Table 1 gives quality results for each image on Fig. 10 according to the theoretical signal Fig. 7(b) chosen as reference. Quality results from the images Figs. 10(a)–10(c) show that iterative reconstructions are better than BFP, especially the OSEM method. Results from images Figs. 10(d)–10(f) highlight the quality loss induced by the Gaussian beam acquisition : SSIM decreases at least by 6%. The detail of SSIM metrics shows that the geometrical properties of the reconstructed images are degraded. The optimized techniques proposed in this paper allow a better reconstruction of the geometry whatever the method. For instance, r increases by 5% with OSEM, and SSIM is globally more accurate with optimized BFP and OSEM. However, it generates luminance and contrast losses with the SART method.
Table 1. Comparison between reconstruction from sinogram acquired through a Gaussian beam and the original domain. It relates on the quality rate of each method to reconstruct the original image.
To complete the quality analysis, an accuracy study is proposed by a comparison of the horizontal profiles of the top bar. Profiles given by the BFP, SART and OSEM reconstructions are shown in Figs. 11(a)–11(c). Whatever the method, the theoretical profile (red curves) is not well recovered even if it is reconstructed from a theoretical sinogram (green curves). However, it is close to the great profile whereas reconstructions from Gaussian beam simulation lead to a spreading of the profile (blue and violet curves). Using the optimized techniques, the spreading is significantly reduced (violet curves), especially using the optimized SART.
Fig. 11 Top bar profiles reconstructed with BFP (a), SART (b) and OSEM (c) : from the theoretical sinogram (Green) and from the simulated Gaussian acquisition with usual (blue) and optimized (Violet) methods. (d) : Top bar profiles reconstructed with the optimized methods. (Green) : BFP. (Blue) : SART. (Violet) : OSEM. (Red) : Theorical signal.
Figure 11(d) compares the profiles obtained by the optimized techniques and the theoretical signal. These results confirm that the optimized SART gives the minimal spreading, followed by the OSEM and finally by the BFP technique. Then, for a given acquisition, SART provides a better accuracy than BFP and OSEM.
4.2. 2D and 3D reconstructions from real acquisitions
Reconstructions computed from the real acquisition Fig. 8(c) are proposed in Fig. 12. The lesser spreading of the metallic bars observed with the optimized methods confirm the theoretical analysis results discussed above.
Fig. 12 (a)(b)(c) Usual BFP, SART and OSEM results from the real acquisition Fig. 8(d). (d)(e)(f) Optimized BFP, SART and OSEM results.
Images given in Fig. 13 show 3D reconstructions of the medicine box. These visualizations are obtained from the stacks of reconstructed cross-sections, computed one-by-one from the 3D acquisition (stack of sinograms). Using the standard methods, each cross-sectional image is directly computed from its corresponding sinogram. Inversely, the optimized methods need all the sinograms at once to take into account the 3D propagation beam model. From the results, we can remark a better quality and accuracy of details in general. For example, the contact area (defined by the high-intensity pixels) between the tablet and the box is easier to estimate. Indeed, as we can see in the cross-sectional images in Fig. 14, this area is less thick on the results obtained with optimized methods. Such qualitative analysis improvement highlights the need to take into account the propagation beam into the tomographic reconstructions for THz CT.
Fig. 13 3D reconstructions of the medicine box using usual BFP (a), SART (b) and OSEM (c) compared to the optimized BFP (d), SART (e) and OSEM (f) methods.
Fig. 14 Estimation of the contact area between the tablet and the medicine box (red pixels) on the images obtained with standard BFP (a), SART (b) and OSEM (c) in one hand and optimized BFP (d), SART (e) and OSEM (f) on the other hand.
5. Conclusion
We demonstrated a development of an improved beam propagation model taking into account the physical properties of THz waves used in THz computed tomography (CT) scan imaging. This more realistic model is clearly justified by the fact that previous models were simply relying on well-known and commonly applied X-Ray reconstruction methods, which assume a ray tracing approach, whereas the THz beam profile is very far from this basic hypothesis. Therefore, we included in an acquisition simulator the Gaussian beam intensity profile in order to observe and to estimate its impact on the projection sets. The simulated sinogram characteristics are found closer to the real acquisition and demonstrate the better reliability of our model to simulate the THz physical acquisition process. Moreover, the discrepancy between Gaussian beam and ray-tracing approaches has been studied. At least, the model is introduced in several inversion methods as a convolution filter to perform efficient tomographic reconstructions of simulated and real acquired objects. Results obtained with three different reconstruction methods (BFP, SART and OSEM) have been compared on the basis of the efficiency of optimized algorithms. We demonstrate an increase of the overall quality and accuracy of the reconstructions. The extra time computation due to these optimized algorithms (≈ 50%, 9 s to 15 s with SART) is negligeable with respect to the current acquisition time (≈ 5 h). This rigorous quantitative approach, associated with a detector array sensor (THz camera), will provide a new interest of 3D THz CT imaging.
Acknowledgments
Development and Optimization of THz NDT of Aeronautics Composite multilayered structures) - Contract N 266320 / FP7-AAT-2010-RTD-1 - http://www.dotnac-project.eu/.
References and links
1. B. Ferguson, S. Wang, D. Gray, D. Abbot, and X. C. Zhang, “T-ray computed tomography,” Opt. Lett. 27, 1312–1314 (2002). [CrossRef]
2. S. Wang, B. Ferguson, D. Abbott, and X. C. Zhang, “T-ray imaging and tomography,” J. Biol. Phys. 29, 247–256 (2003). [CrossRef]
3. S. Wang and X. C. Zhang, “Pulsed terahertz tomography,” J. Phys. D: Appl. Phys. 37, R1–R36 (2004). [CrossRef]
4. M. M. Awad and R. A. Cheville, “Transmission terahertz waveguide-based imaging below the diffraction limit,” Appl. Phys. Lett. 86, 221107 (2005). [CrossRef]
5. X. Yin, B. W. H. Ng, B. Ferguson, and D. Abbott, “Wavelet based local tomographic image using terahertz techniques,” Digital Signal Process. 19, 750–763 (2009). [CrossRef]
6. A. Brahm, M. Kunz, S. Riehemann, G. Notni, and A. Tünnermann, “Volumetric spectral analysis of materials using terahertz-tomography techniques,” Appl. Phys. B 100, 151–158 (2010). [CrossRef]
7. E. Abraham, A. Younus, C. Aguerre, P. Desbarats, and P. Mounaix, “Refraction losses in terahertz computed tomography,” Opt. Commun. 283, 2050–2055 (2010). [CrossRef]
8. S. Nadar, H. Videlier, D. Coquillat, F. Teppe, and M. Sakowicz, “Room temperature imaging at 1.63 and 2.54 THz with field effect transistor detectors,” J. Appl. Phys. 108, 054508 (2010). [CrossRef]
9. A. El Fatimy, J.C. Delagnes, A. Younus, E. Nguema, F. Teppe, W. Knap, E. Abraham, and P. Mounaix, “Plasma wave field effect transistors as a resonant detector for imaging applications up to one terahertz for terahertz imaging”, Opt. Commun. 282(15), 3055–3058 (2009). [CrossRef]
10. A. H. Andersen and A. C. Kak, “Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm,” Ultrasonic Imaging 6, 81–94 (1984). [CrossRef] [PubMed]
11. L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Trans. Med. Imaging 1, 113–122 (1982). [CrossRef] [PubMed]
12. H. M. Hudson and R. S. Larkin, “Accelerated image reconstruction using ordered subsets of projection data,” IEEE Trans. Med. Imaging 13, 601–609 (1994). [CrossRef] [PubMed]
13. A. Younus, S. Salort, B. Recur, P. Desbarats, P. Mounaix, J-P. Caumes, and E. Abraham, “3D millimeter wave tomographic scanner for large size opaque object inspection with different refractive index contrasts,” in Millimetre Wave and Terahertz Sensors and Technology III, K.A. Krapels and N.A. Salmon, eds., Proc. SPIE7837, 783709 (2010). [CrossRef]
14. G. T. Herman, Image Reconstruction From Projections : The Fundamentals of Computerized Tomography (Academic Press, 1980).
15. P. Toft, “The Radon Transform : Theory and Implementation,” Ph.D. thesis, Department of Mathematical Modelling, Section for Digital Signal Processing, Technical University of Denmark (1996).
16. J. Radon, “Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten,” Ber. Ver. Sachs. Akad. Wiss. Leipzig, Math-Phys. Kl 69, 262–277 (1917). In German. An english translation can be found in S. R. Deans: The Radon Transform and Some of Its Applications.
17. B. Recur, A. Younus, S. Salort, P. Mounaix, B. Chassagne, P. Desbarats, J-P. Caumes, and E. Abraham, “Investigation on reconstruction methods applied to 3D terahertz computed tomography,” Opt. Express 19, 5105–5117 (2011). [CrossRef] [PubMed]
18. R. Gordon, R. Bender, and G. T. Herman, “Algebraic Reconstruction Techniques (ART) for Three-dimensional Electron Microscopy and X-ray Photography,” J. Theor. Biol. 29, 471–481 (1970). [CrossRef] [PubMed]
19. A. P. Dhawan, R. M. Rangayyan, and R. Gordon, “Image restoration by Wiener deconvolution in limited-view computed tomography,” Appl. Opt. 24, 4013–4020 (1985). [CrossRef] [PubMed]
20. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment : From error visibility to structural similarity,” IEEE Trans. Image Process. 13, 600–612 (2004). [CrossRef] [PubMed]
