OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 19, Iss. 6 — Mar. 14, 2011
  • pp: 5105–5117
« Show journal navigation

Investigation on reconstruction methods applied to 3D terahertz computed tomography

B. Recur, A. Younus, S. Salort, P. Mounaix, B. Chassagne, P. Desbarats, J-P. Caumes, and E. Abraham  »View Author Affiliations


Optics Express, Vol. 19, Issue 6, pp. 5105-5117 (2011)
http://dx.doi.org/10.1364/OE.19.005105


View Full Text Article

Acrobat PDF (1090 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

3D terahertz computed tomography has been performed using a monochromatic millimeter wave imaging system coupled with an infrared temperature sensor. Three different reconstruction methods (standard back-projection algorithm and two iterative analysis) have been compared in order to reconstruct large size 3D objects. The quality (intensity, contrast and geometric preservation) of reconstructed cross-sectional images has been discussed together with the optimization of the number of projections. Final demonstration to real-life 3D objects has been processed to illustrate the potential of the reconstruction methods for applied terahertz tomography.

© 2011 Optical Society of America

1. Introduction

In the field of 3D imaging, X-Ray Computed Tomography (CT) is an ubiquitous technique which provides cross-sectional images of an object by analyzing the radiation transmitted by the sample through different incidence angles. This technique can visualize 3D dense materials such as wood, bones and biological tissues but cannot be easily applied to soft materials such as plastics, papers or paintings owing to the low absorption of the X-Ray radiation. Alternatively, terahertz (THz) radiation proposes attractive features such as non-invasive and non-destructive analysis, transparency and good penetration depth through various materials [1

1. W. L. Chan, J. Deibel, and D. M. Mittleman, “Imaging with terahertz radiation,” Rep. Prog. Phys. 70, 1325–1379 (2007). [CrossRef]

]. All these properties make THz radiation very promising for direct applications in non-destructive inspection, detection of illicit drugs [2

2. K. Kawase, Y. Ogawa, Y. Watanabe, and H. Inoue, “Non-destructive terahertz imaging of illicit drugs using spectral fingerprints,” Opt. Express 11, 2549–2554 (2003). [CrossRef] [PubMed]

] and explosives [3

3. Y. C. Shen, T. Lo, P. F. Taday, B. E. Cole, W. R. Tribe, and M. C. Kemp, “Detection and identification of explosives using terahertz pulsed spectroscopic imaging,” Appl. Phys. Lett. 86, 241116 (2005). [CrossRef]

], art conservation [4

4. K. Fukunaga and M. Picollo, “Terahertz spectroscopy applied to the analysis of artists materials,” Appl. Phys. A 100, 591–597 (2010). [CrossRef]

, 5

5. E. Abraham, A. Younus, J.-C. Delagnes, and P. Mounaix, “Non-invasive investigation of art paintings by terahertz imaging,” Appl. Phys. A 100, 585–590 (2010). [CrossRef]

] and biomedical applications [6

6. S. Y. Huang, Y. X. J. Wang, D. K. W. Yeung, A. T. Ahuja, Y. T. Zhang, and E. Pickwell-MacPherson, “Tissue characterization using terahertz pulsed imaging in reflection geometry,” Phys. Med. Biol. 54, 149–160 (2009). [CrossRef]

]. For 3D THz imaging, reflection tomography and cross section measurements have been proposed using for instance the time-of-flight of reflected THz pulses [7

7. J. Takayanagi, H. Jinno, S. Ichino, K. Suizu, M. Yamashita, T. Ouchi, S. Kasai, H. Ohtake, H. Uchida, N. Nishizawa, and K. Kawase, “High-resolution time-of-flight terahertz tomography using a femtosecond fiber laser,” Opt. Express 17, 7549–7555 (2009). [CrossRef] [PubMed]

, 8

8. K. Iwaszczuk, H. Heiselberg, and P. U. Jepsen, “Terahertz radar cross section measurements,” Opt. Express 18, 26399–26408 (2010). [CrossRef] [PubMed]

]. Especially in 2002, THz CT has been proposed for 3D imaging of opaque materials [9

9. B. Ferguson, S. Wang, D. Gray, D. Abbot, and X. C. Zhang, “T-ray computed tomography,” Opt. Lett. 27, 1312–1314 (2002). [CrossRef]

]. Ferguson et al. demonstrated that cross-sectional images can be obtained by measuring the transmitted amplitude and phase of broadband THz pulses at multiple projection angles. However, it has been emphasized that several peaks in the THz waveform obtained with a time-domain spectrometer strongly complicate the signal analysis [10

10. E. Abraham, A. Younus, C. Aguerre, P. Desbarats, and P. Mounaix, “Refraction losses in terahertz computed tomography,” Opt. Commun. 283, 2050–2055 (2010). [CrossRef]

]. Although THz CT seems powerful, few papers have been published since the technique suffers from strong limitations [9

9. B. Ferguson, S. Wang, D. Gray, D. Abbot, and X. C. Zhang, “T-ray computed tomography,” Opt. Lett. 27, 1312–1314 (2002). [CrossRef]

, 11

11. S. Wang, B. Ferguson, D. Abbott, and X. C. Zhang, “T-ray imaging and tomography,” J. Biol. Phys. 29, 247–256 (2003). [CrossRef]

, 12

12. S. Wang and X. C. Zhang, “Pulsed terahertz tomography,” J. Phys. D: Appl. Phys. 37, R1–R36 (2004). [CrossRef]

, 13

13. M. M. Awad and R. A. Cheville, “Transmission terahertz waveguide-based imaging below the diffraction limit,” Appl. Phys. Lett. 86, 221107 (2005). [CrossRef]

, 14

14. X. Yin, B. W. H. Ng, B. Ferguson, and D. Abbott, “Wavelet based local tomographic image using terahertz techniques,” Digit. Signal Process. 19, 750–763 (2009). [CrossRef]

, 15

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

].

The first limitation concerns the diffraction effects and Fresnel losses experienced by the propagation of the THz wave through the sample [9

9. B. Ferguson, S. Wang, D. Gray, D. Abbot, and X. C. Zhang, “T-ray computed tomography,” Opt. Lett. 27, 1312–1314 (2002). [CrossRef]

]. Until now, most samples that have been imaged using THz CT were made of polystyrene or similar low refractive index materials. However, as soon as the refractive index of the sample is in the order or greater than 1.5, which represents the majority of realistic samples, the THz beam is strongly refracted by the sample and the transmitted signal is very difficult to detect. Therefore, in THz CT, it is always difficult to discriminate if the imaging contrast arises from sample absorption or deviation of the THz beam by refraction. One possibility of differentiation is to perform pulsed THz CT and measure the phase of the transmitted signal [9

9. B. Ferguson, S. Wang, D. Gray, D. Abbot, and X. C. Zhang, “T-ray computed tomography,” Opt. Lett. 27, 1312–1314 (2002). [CrossRef]

].

The second limitation is the long acquisition time required by pulsed THz system based on time-domain spectrometer, since it used a point-to-point measurement associated with a temporal sampling and the rotation of the sample. This limitation can be reduced by the utilization of a continuous wave THz source even if in this case the phase information of the object can not be extracted [16

16. K. L. Nguyen, M. L. Johns, L. F. Gladden, C. H. Worral, P. Alexander, H. E. Beere, M. Pepper, D. A. Ritchie, J. Alton, S. Barbieri, and E. H. Linfield, “Three-dimensional imaging with a terahertz quantum cascade laser,” Opt. Express 14, 2123–2129 (2006). [CrossRef] [PubMed]

, 17

17. 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. SPIE 7837, 783709 (2010).

]. The limitation is also directly connected to the number of projection data. To reduce this parameter, a recent study proposed the development of depth-resolving THz imaging with tomosynthesis, which is similar to CT except that the number of projections is much smaller [18

18. N. Sunaguchi, Y. Sasaki, N. Maikusa, M. Kawai, T. Yuasa, and C. Otani, “Depth-resolving terahertz imaging with tomosynthesis,” Opt. Express 17, 9558–9570 (2009). [CrossRef] [PubMed]

]. The authors used only five projections instead of generally 15 to 20 in THz CT. However, the efficiency of their system was mainly limited to thin and wide samples (50 sheets of post-it notes). Another method to perform fast THz imaging consists in illuminating the sample by a focused THz line and detecting the THz radiation by non-collinear electro-optical time-to-space conversion in a non-linear crystal [19

19. T. Yasuda, T. Yasui, T. Araki, and E. Abraham, “Real-time two-dimensional terahertz tomography of moving objects,” Opt. Commun. 267, 128–136 (2006). [CrossRef]

, 20

20. T. Yasui, K. Sawanaka, A. Ihara, E. Abraham, M. Hashimoto, and T. Araki, “Real-time terahertz color scanner for moving objects,” Opt. Express 16, 1208–1221 (2008). [CrossRef] [PubMed]

].

Finally, it is obvious that 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. Usually, a back-projection of the filtered projections (BFP) is employed as the standard reconstruction method [21

21. G. T. Herman, Image Reconstruction From Projections : The Fundamentals of Computerized Tomography (Academic Press Inc., 1980).

]. The analysis is based on Radon inverse transformations from the projection data. In the case of THz CT, to our knowledge, only this BFP method has been employed since it is proposed in most conventional mathematical softwares. However, in X-Ray CT, it is well-known that BFP suffers from several disadvantages such as beam hardening, noise sensitivity and geometric degradation in case of insufficient number of projections. Consequently, alternative iterative reconstruction methods have been proposed for X-Ray CT such as the Simultaneous Algebraic Reconstruction Technique (SART) [22

22. A. H. Andersen and A. C. Kak, “Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm,” Ultrason. Imaging 6, 81–94 (1984). [CrossRef] [PubMed]

] and the Ordered Subsets Expectation Maximization (OSEM) method [23

23. L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Trans. Med. Imaging1, 113–122 (1982). [CrossRef] [PubMed]

, 24

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

]. These alternative iterative reconstruction methods have never been applied to THz CT. In this paper, to the first time to our knowledge, we propose to test and compare these three reconstruction methods (BFP, SART, OSEM) applied to THz CT. Especially, a particular interest will concern the optimization of the number of projections associated with the preservation of the image quality. For this study, we employed a flexible and easy-to-align system based on a continuous millimeter wave source and a commercial pyroelectric sensor.

The paper is organized as follows. In the first part, the experimental setup and the spatial resolution of the system will be presented. Then, the reconstruction methods will be exposed and evaluated taking into account the number of projections and the quality (intensity, contrast, noise, geometric preservation) of the final reconstructed cross-sections. Finally, challenges for precise and efficient tomographic reconstruction of complex refraction index materials will be underlined in order to demonstrate the potential of the proposed 3D millimeter wave tomo-graphic scanner.

2. Experimental setup

The experimental setup of the 3D millimeter wave tomographic scanner is shown in Fig. 1(a). The output beam of a compact millimeter wave Gunn diode coupled with a horn antenna is collimated using an off-axis parabolic mirror M (f′ = 150 mm). Two different Gunn diodes have been used: a 80 GHz diode coupled to a frequency tripler delivering 0.3 mW at 240 GHz and a 110 GHz diode delivering 20 mW. The advantage of the 240 GHz source is essentially the better spatial resolution (wavelength 1.25 mm), whereas the second source has a much larger output power. The THz beam is then focused with a Teflon lens L (f′ = 60 mm) on the sample S which is positioned on three-axes XY-θ motorized stages, the angle θ corresponding to a rotation around the vertical Y-axis. For the 110 GHz source, Fig.1(b) shows the 2D transversal profile of the THz beam at the beam waist visualized using a photothermal THz convertor (Teracam-Alphanov, [25

25. C. Pradere, J.-P. Caumes, D. Balageas, S. Salort, E. Abraham, B. Chassagne, and J.-C. Batsale, “Photothermal converters for quantitative 2D and 3D real-time terahertz imaging,” Quant. Infrared Thermog. 7, 217–235 (2010). [CrossRef]

]). This result shows that, at the sample position S, the beam profile is homogeneous with a Gaussian circular shape (4 mm beam diameter, measured at FWHM) in agreement with the theoretical values obtained from the propagation of Gaussian beams. Similarly (data not shown), at the same position, a FWHM beam diameter of 2 mm has been measured with the 240 GHz source. These results indicate 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. Consequently, the system is more adapted for the visualization of sub-centimeter structures within large size object, typically more than (100 × 100 × 100) mm3. Moreover, to properly apply THz CT reconstruction algorithms and avoid depth resolution degradation, the target’s lateral extent in the direction of the THz wave vector has to be in the order of the confocal parameter of the THz beam [11

11. S. Wang, B. Ferguson, D. Abbott, and X. C. Zhang, “T-ray imaging and tomography,” J. Biol. Phys. 29, 247–256 (2003). [CrossRef]

]. Here, we can estimate that this confocal parameter is about 20 mm, which is compatible with the size of the tested samples. Symmetrically, the focal volume is imaged using a similar arrangement on a commercial low-cost pyroelectric sensor (Spectrum Detector Inc.). For data acquisition with the pyroelectric sensor, the THz beam is modulated at 20 Hz by an optical chopper and the amplitude of the transmitted THz is acquired with a lock-in amplifier.

Fig. 1 (a) Experimental setup. C: optical chopper, L: Teflon lens (f′ = 60 mm), M: off-axis parabolic mirror (f′ = 150 mm), S: sample, D: pyroelectric detector. (b) 2D spatial profile of the THz beam waist at the sample position (110 GHz source) visualized with a photothermal THz convertor (Teracam).

In order to perform a 3D reconstruction of volumetric object, the experimental procedure is organized as follows [17

17. 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. SPIE 7837, 783709 (2010).

]. First, we record a 2D transmission image of the sample by moving the object in the X and Y directions with a scan step of 0.5 mm or 1 mm in both directions. With a scan speed of 5 pixels/sec, the acquisition time for a (100 × 100) pixels image size is about 30 minutes. Then, the sample is rotated by a rotation step in order to provide a different visualization of the object. The operation is repeated Nθ times from θ = 0° to θ = 180° and we finally get a set of Nθ projections of the sample corresponding to the different angles of visualization. For instance, if Nθ = 18, = 10° and the total acquisition time is about 9 hours for the complete 3D visualization of the sample. However, to visualize only a single cross-sectional image of the sample, corresponding to a 100 pixels horizontal scan associated to 18 different projections, this acquisition time is reduced to approximately 6 minutes. From these projection data, we are able to construct the sinogram of the object which represents, for a given horizontal slice, the evolution of the transmitted THz amplitude as a function of the rotation angle. Finally we apply the BFP, SART or OSEM algorithms to reconstruct the final 3D volume of the sample.

3. Tomographic reconstruction methods

Tomography is used to reconstruct the volume of an object from the set of projections done from the exterior of the object [21

21. G. T. Herman, Image Reconstruction From Projections : The Fundamentals of Computerized Tomography (Academic Press Inc., 1980).

, 26

26. 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).

]. This technique, widely developed in X-Ray CT scan imaging, is modeled by the Radon theorem [27

27. J. Radon, “Über die Bestimmung von Funktionen durch ihre Integralwerte langs 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.

]. The direct transform ℛ describes the projection line acquisition. It maps a 2D function defined by f(x,y) into a 1D projection along an angle θ and a module ρ. It is defined by the following formula:
θ(ρ)=f(x,y)δ(ρxcosθysinθ)dxdy
(1)
where θ and ρ are respectively the angular and radial coordinates of the projection line (θ,ρ) and δ is the Dirac impulse. Then, ℛθ(ρ) represents the absorption sum undergone by the ray along the line. As an illustration, let us consider a black foam parallelepiped (absorbance α = 0.02 cm−1, refractive index 1.1) with a (41 × 49)mm2 cross section. The parallelepiped is drilled by two holes (15 mm diameter): one hole with air and another hole filled by a Teflon cylinder (absorbance α < 1 cm−1, refractive index 1.45) having a small cylindrical air hole inside (6 mm diameter) (Fig. 2(a)). For one horizontal cross-section, the acquisition along several angles gives a sinogram composed with Nθ lines, corresponding to the number of projections, and Nρ columns corresponding to the number of pixels in the horizontal direction. For instance, Fig. 2(b) represents the sinogram acquired with the 110 GHz source with Nθ = 72 projections ( = 2.5°) and Nρ = 128 horizontal pixels (0.5 mm scan step). The black color corresponds to a maximum transmitted signal whereas the white color indicates a decrease of the transmitted signal. Similar color coding has been applied throughout the paper. Low signal intensity can be due to either absorption either deviation of the THz beam due to refraction. Here, owing to the low absorbance of Teflon, mainly refraction losses have to be considered since the Teflon cylinder acts as a lens. Before the application of any reconstruction method, this sinogram already reflects the global geometric structure of the sample cross-section with especially the presence of the refracting Teflon material. A more precise description of the THz imaging of this sample will be presented in section 4.

Fig. 2 (a) Photograph of original object: parallelepiped black foam (41 × 49) mm2 with 2 holes, diameter 15 mm (1 hole with air and 1 hole containing a Teflon cylinder with a 6 mm cylindrical air hole inside). (b) Sinogram with Nθ = 72 projections (lines) and Nρ = 128 samples per projection (columns). Acquisition with the 110 GHz source.

The reconstruction process is given by the inverse Radon transform which recovers the original domain from the projections. Given a sinogram 𝒮 containing several ℛθ(ρ) values (θ ∈ [0,π[ and ρ ∈ ℝ), the discrete inverse Radon transform computes each image pixel as follows:
I(i,j)=iθ=0Nθ1ρ=Nρ2Nρ2Wθ(iθ)(ρ)A(θ,ρ),(i,j)
(2)
where:
  • θ(iθ)=iθπNθ,,
  • A(θ,ρ),(i,j) is the weight-matrix defining the weight value between each pixel and each projection line,
  • Wθ(ρ)=ν=Nρ2Nρ2|ν|(ρs=Nρ2Nρ2θ(ρs)ei2πρsν)ei2πρν.

First, this inversion filters each projection ℛθ in the Fourier domain with the ramp filter |ν| to increase geometric details. Second, it computes the pixel values from the filtered projections Wθ(iθ). This method is thus denoted back-projection of filtered projections (BFP). As already pointed out in the introduction, this technique is widely used in THz CT imaging since it is proposed in most of CT software tools.

However, it is known that BFP suffers from several drawbacks such as noise sensitivity and beam hardening. Especially the beam hardening phenomenon induces artifacts such as cupping, streaks and blurring because rays from some projection angles are hardened to a differing extent than rays from other angles, confusing the reconstruction algorithm. Moreover BFP is very sensitive to the projection number Nθ. If this number is too low, beam hardening artifacts will be enhanced. This number Nθ is also very critical because it acts on the global acquisition time, which is particularly important in THz CT imaging. Consequently, in this paper, we also investigated other reconstruction methods allowing to reduce the projection number (i.e. the acquisition time) while preserving reconstruction quality.

Particularly, several CT imaging systems are based on iterative reconstructions such as SART [22

22. A. H. Andersen and A. C. Kak, “Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm,” Ultrason. Imaging 6, 81–94 (1984). [CrossRef] [PubMed]

] or OSEM [24

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

] algorithms. The SART algebraic method is based on the Karcz-marz algorithm used to approach the solution of the linear equation system I = ATR, where I is the image, R is the sinogram and A is the weight-matrix [28

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

]. SART is 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 Rθsk (computed from Ik,s−1 by using direct Radon transform). A super-iteration k is over when all the projections have been used. Consequently, pixel update using SART is computed as follows:
Ik,s(i,j)=Ik,s1(i,j)+λρ=0Nρ1A(θs,p),(i,j)[θs(ρ)Rθsk(ρ)Dθs(ρ)]ρ=0Nρ1A(θs,ρ),(i,j)
(3)
where Dθs(ρ)=i=0W1j=0H1A(θs,ρ),(i,j) is the norm of the segment (θs, ρ) crossing the image and (W × H) is the image size. Iterations are performed until the convergence of the solution. The initial I0,0 image is usually an uniform image.

OSEM algorithm [23

23. L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Trans. Med. Imaging1, 113–122 (1982). [CrossRef] [PubMed]

] is another iterative process which slightly differs from the SART. The update is done from a subset of several projections at once and the error correction is multiplicative:
Ik+1(i,j)=I(i,j)iθ=0Nθ1ρ=0Nρ1A(θ,ρ),(i,j)θ(ρ)Rθk(ρ)iθ=0Nθ1ρ=0Nρ1A(θ,ρ),(i,j)
(4)

With the OSEM method, the convergence of the solution is longer than with the SART method because OSEM directly uses all the projections at once. Moreover, after convergence to the solution, OSEM is very sensitive to any solution divergence so that it needs non-trivial regularizations to reduce induced artifacts. However, OSEM reconstruction quality is commonly higher than SART and BFP ones. If we assume that the complete computation time is 1 second for BFP, we can estimate that this time increases to 4 seconds for SART and 7 seconds for OSEM. In all cases, even if the OSEM computation time is longer than other methods, it is still negligible compared to the long acquisition time required for THz CT.

4. Quality and accuracy according to the projection number and the reconstruction method

Fig. 3 Sinograms of two metallic bars (12 mm diameter) separated by 50 mm, with a projection number Nθ = 72. (a) Ideal theoretical sinogram calculated by direct Radon transformation of a synthetic model. (b) Acquisition using the millimeter wave tomographic scanner with the 110 GHz source. Same scale and range as in (a).
Fig. 4 Cross sections of two metallic bars (12 mm diameter) separated by 50 mm. (a) Ideal synthetic cross section of the sample. (b) BFP reconstruction from 3(b). (c) SART reconstruction from 3(b). (d) OSEM reconstruction from 3(b). (b) to (d): same scale as in (a).
Fig. 5 Intensity profiles along the horizontal line intercepting the center of both metallic bars (from Fig. 4).

Fig. 6 Manufactured sample presented in Fig. 2. Reconstructions using sinograms with 12, 18, 24, 36, 72 projections and the BFP (a), SART (b) and OSEM (c) methods. Same scale has been used for all cross sections.

Afterwards, we analyzed the image quality on a quantitative point of view. Numerically, a complete comparison of BFP, SART and OSEM has already been performed from synthetic digital images [29

29. B. Recur, “Qualité et Précision en Reconstruction Tomographique : Algorithmes et Applications,” Ph.D. thesis, LaBRI, Bordeaux 1 University (2010).

]. It clearly indicates that, if the number of projections is sufficient (to be defined depending on the experimental technique), OSEM image quality is worse (especially contour accuracy) than BFP and SART images. However, in case of THz CT, we are mostly interested in the evolution of the image quality as a function of the projection number. Especially, the key-point concerns the preservation of the image quality for a limited number of projections (typically less than 20 projections). Quantitatively, the image quality can be determined using the structural similarity (SSIM) parameter measured between a reference image I and a transformed image It. The SSIM parameter is given by [30

30. 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. 13600–612 (2004). [CrossRef] [PubMed]

]:
SSIM(I,It)=l(I,It)c(I,It)r(I,It)
(5)
where l(I,It), c(I,It) and r(I,It) are the intensity, contrast and geometric equivalence rates between both images, respectively. The image I is the reference corresponding to the 72 projection sinogram, whereas the images It are reconstructed from Nθ = 12 to Nθ = 36. Table 1 details the three equivalence rates obtained by BFP (Table 1(a)), SART (Table 1(b)) and OSEM (Table 1(c)) methods according to the projection number. Here, it is important to notice that these tables are independent and cannot be compared each others since, for each table, the reconstructed image obtained with Nθ = 72 is considered as the reference with all metric values set to unity. However, for each reconstruction method and depending on the projection number, the tables indicate the image degradation with respect to the intensity, contrast and geometric preservation. We can notice that all metric values slightly decrease whatever the method, but r(I,It) especially decreases significantly with BFP (indicated in red color in Table 1(a)). Then, the SSIM parameter of BFP is mainly deteriorated because the geometrical aspect of the image is not well preserved. For SART, the global image quality is already optimized with Nθ = 18. For SART, it appears that the method is almost insensitive to the projection number. However, as pointed out previously, we have to keep in mind that this method also provides a general vagueness in the final reconstructed image.

Table 1. Details of the SSIM parameter for BFP (a), SART (b) and OSEM (c). l(I,It): intensity, c(I,It): contrast, r(I,It): geometric equivalence rates. Red color indicates significant results.

table-icon
View This Table

For each reconstruction method, the evolution of the SSIM parameter is also represented as a function of the projection number Nθ in Fig. 7. Independently from each others, the curves illustrate the reconstruction quality losses depending on the number of projections. As previously noticed in Table 1, we can remark that the SART and OSEM algorithms compute almost without quality loss whatever the projection number. Inversely, the BFP quality decreases if the projection number is less than 25. This observation is essential in case of THz CT since the projection number is limited owing to the long acquisition time. Here, we clearly point out that, for a limited number of projection data (typically less than 20), SART and OSEM methods are more appropriated for efficient THz CT.

Fig. 7 Manufactured sample presented in Fig. 2. SSIM parameter as a function of the projection number and the reconstruction method.

5. 3D reconstructions

Fig. 8 White foam parallelepiped (30 mm cube size) drilled by two oblique metallic bars (6 mm diameter). (a) Photograph of the 3D sample. (b) BFP reconstruction ( Media 1). (c) SART reconstruction ( Media 2). (d) OSEM reconstruction ( Media 3). Acquisition with the 240 GHz source.

The last sample is a wooden Russian doll Matriochka (total height 160 mm, Fig. 9(a)). The THz 3D video recordings of the sample volume reconstruction (Fig. 9(b) for BFP ( Media 4), Fig. 9(c) for SART ( Media 5), Fig. 9(d) for OSEM ( Media 6)) clearly reveal the shape of the outer doll and the presence of the smaller doll inside the main doll (measured size 95 mm). This last example of 3D reconstructions of complex objects emphasize the potential of the SART method which provides high quality contour and reasonable artifacts despite the insufficient number of projections.

Fig. 9 Wooden Russian doll Matriochka (total height 160 mm). (a) Photograph of the 3D sample. (b) BFP reconstruction ( Media 4). (c) SART reconstruction ( Media 5). (d) OSEM reconstruction ( Media 6). Acquisition with the 110 GHz source.

6. Conclusion

With a millimeter wave source associated with an infrared temperature sensor, we investigated the reconstruction methods applied to THz CT. The system is easy-to-align, low-cost and portable which represents a major advantage compared to standard time-domain spectroscopic imaging. After comparing with the commonly used BFP reconstruction algorithm, we pointed out the ability of the SART and OSEM algorithms to properly reconstruct cross-sectional images. The SART method applied to THz CT obtains equivalent accuracy and quality than BFP. However, for a limited number of projections (less than 25), we noticed a quantitative degradation of the BFP reconstruction whereas SART and OSEM methods can already offer an optimized reconstruction quality. Consequently, iterative methods allow to reduce the acquisition time (i.e. the projection number) without noticeable quality and accuracy losses. These quantitative results, well-known in standard tomography such as X-Ray CT, are here pointed out for THz CT imaging, which provides a new insight for the improvement of 3D THz imaging.

References and links

1.

W. L. Chan, J. Deibel, and D. M. Mittleman, “Imaging with terahertz radiation,” Rep. Prog. Phys. 70, 1325–1379 (2007). [CrossRef]

2.

K. Kawase, Y. Ogawa, Y. Watanabe, and H. Inoue, “Non-destructive terahertz imaging of illicit drugs using spectral fingerprints,” Opt. Express 11, 2549–2554 (2003). [CrossRef] [PubMed]

3.

Y. C. Shen, T. Lo, P. F. Taday, B. E. Cole, W. R. Tribe, and M. C. Kemp, “Detection and identification of explosives using terahertz pulsed spectroscopic imaging,” Appl. Phys. Lett. 86, 241116 (2005). [CrossRef]

4.

K. Fukunaga and M. Picollo, “Terahertz spectroscopy applied to the analysis of artists materials,” Appl. Phys. A 100, 591–597 (2010). [CrossRef]

5.

E. Abraham, A. Younus, J.-C. Delagnes, and P. Mounaix, “Non-invasive investigation of art paintings by terahertz imaging,” Appl. Phys. A 100, 585–590 (2010). [CrossRef]

6.

S. Y. Huang, Y. X. J. Wang, D. K. W. Yeung, A. T. Ahuja, Y. T. Zhang, and E. Pickwell-MacPherson, “Tissue characterization using terahertz pulsed imaging in reflection geometry,” Phys. Med. Biol. 54, 149–160 (2009). [CrossRef]

7.

J. Takayanagi, H. Jinno, S. Ichino, K. Suizu, M. Yamashita, T. Ouchi, S. Kasai, H. Ohtake, H. Uchida, N. Nishizawa, and K. Kawase, “High-resolution time-of-flight terahertz tomography using a femtosecond fiber laser,” Opt. Express 17, 7549–7555 (2009). [CrossRef] [PubMed]

8.

K. Iwaszczuk, H. Heiselberg, and P. U. Jepsen, “Terahertz radar cross section measurements,” Opt. Express 18, 26399–26408 (2010). [CrossRef] [PubMed]

9.

B. Ferguson, S. Wang, D. Gray, D. Abbot, and X. C. Zhang, “T-ray computed tomography,” Opt. Lett. 27, 1312–1314 (2002). [CrossRef]

10.

E. Abraham, A. Younus, C. Aguerre, P. Desbarats, and P. Mounaix, “Refraction losses in terahertz computed tomography,” Opt. Commun. 283, 2050–2055 (2010). [CrossRef]

11.

S. Wang, B. Ferguson, D. Abbott, and X. C. Zhang, “T-ray imaging and tomography,” J. Biol. Phys. 29, 247–256 (2003). [CrossRef]

12.

S. Wang and X. C. Zhang, “Pulsed terahertz tomography,” J. Phys. D: Appl. Phys. 37, R1–R36 (2004). [CrossRef]

13.

M. M. Awad and R. A. Cheville, “Transmission terahertz waveguide-based imaging below the diffraction limit,” Appl. Phys. Lett. 86, 221107 (2005). [CrossRef]

14.

X. Yin, B. W. H. Ng, B. Ferguson, and D. Abbott, “Wavelet based local tomographic image using terahertz techniques,” Digit. Signal Process. 19, 750–763 (2009). [CrossRef]

15.

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]

16.

K. L. Nguyen, M. L. Johns, L. F. Gladden, C. H. Worral, P. Alexander, H. E. Beere, M. Pepper, D. A. Ritchie, J. Alton, S. Barbieri, and E. H. Linfield, “Three-dimensional imaging with a terahertz quantum cascade laser,” Opt. Express 14, 2123–2129 (2006). [CrossRef] [PubMed]

17.

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. SPIE 7837, 783709 (2010).

18.

N. Sunaguchi, Y. Sasaki, N. Maikusa, M. Kawai, T. Yuasa, and C. Otani, “Depth-resolving terahertz imaging with tomosynthesis,” Opt. Express 17, 9558–9570 (2009). [CrossRef] [PubMed]

19.

T. Yasuda, T. Yasui, T. Araki, and E. Abraham, “Real-time two-dimensional terahertz tomography of moving objects,” Opt. Commun. 267, 128–136 (2006). [CrossRef]

20.

T. Yasui, K. Sawanaka, A. Ihara, E. Abraham, M. Hashimoto, and T. Araki, “Real-time terahertz color scanner for moving objects,” Opt. Express 16, 1208–1221 (2008). [CrossRef] [PubMed]

21.

G. T. Herman, Image Reconstruction From Projections : The Fundamentals of Computerized Tomography (Academic Press Inc., 1980).

22.

A. H. Andersen and A. C. Kak, “Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm,” Ultrason. Imaging 6, 81–94 (1984). [CrossRef] [PubMed]

23.

L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Trans. Med. Imaging1, 113–122 (1982). [CrossRef] [PubMed]

24.

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]

25.

C. Pradere, J.-P. Caumes, D. Balageas, S. Salort, E. Abraham, B. Chassagne, and J.-C. Batsale, “Photothermal converters for quantitative 2D and 3D real-time terahertz imaging,” Quant. Infrared Thermog. 7, 217–235 (2010). [CrossRef]

26.

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

27.

J. Radon, “Über die Bestimmung von Funktionen durch ihre Integralwerte langs 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.

28.

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]

29.

B. Recur, “Qualité et Précision en Reconstruction Tomographique : Algorithmes et Applications,” Ph.D. thesis, LaBRI, Bordeaux 1 University (2010).

30.

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. 13600–612 (2004). [CrossRef] [PubMed]

OCIS Codes
(100.6890) Image processing : Three-dimensional image processing
(120.5800) Instrumentation, measurement, and metrology : Scanners
(170.3010) Medical optics and biotechnology : Image reconstruction techniques
(110.6795) Imaging systems : Terahertz imaging
(110.6955) Imaging systems : Tomographic imaging

ToC Category:
Imaging Systems

History
Original Manuscript: January 19, 2011
Revised Manuscript: February 16, 2011
Manuscript Accepted: February 16, 2011
Published: March 2, 2011

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

Citation
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)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-6-5105


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. W. L. Chan, J. Deibel, and D. M. Mittleman, “Imaging with terahertz radiation,” Rep. Prog. Phys. 70, 1325–1379 (2007). [CrossRef]
  2. K. Kawase, Y. Ogawa, Y. Watanabe, and H. Inoue, “Non-destructive terahertz imaging of illicit drugs using spectral fingerprints,” Opt. Express 11, 2549–2554 (2003). [CrossRef] [PubMed]
  3. Y. C. Shen, T. Lo, P. F. Taday, B. E. Cole, W. R. Tribe, and M. C. Kemp, “Detection and identification of explosives using terahertz pulsed spectroscopic imaging,” Appl. Phys. Lett. 86, 241116 (2005). [CrossRef]
  4. K. Fukunaga and M. Picollo, “Terahertz spectroscopy applied to the analysis of artists materials,” Appl. Phys., A Mater. Sci. Process. 100, 591–597 (2010). [CrossRef]
  5. E. Abraham, A. Younus, J.-C. Delagnes, and P. Mounaix, “Non-invasive investigation of art paintings by terahertz imaging,” Appl. Phys., A Mater. Sci. Process. 100, 585–590 (2010). [CrossRef]
  6. S. Y. Huang, Y. X. J. Wang, D. K. W. Yeung, A. T. Ahuja, Y. T. Zhang, and E. Pickwell-MacPherson, “Tissue characterization using terahertz pulsed imaging in reflection geometry,” Phys. Med. Biol. 54, 149–160 (2009). [CrossRef]
  7. J. Takayanagi, H. Jinno, S. Ichino, K. Suizu, M. Yamashita, T. Ouchi, S. Kasai, H. Ohtake, H. Uchida, N. Nishizawa, and K. Kawase, “High-resolution time-of-flight terahertz tomography using a femtosecond fiber laser,” Opt. Express 17, 7549–7555 (2009). [CrossRef] [PubMed]
  8. K. Iwaszczuk, H. Heiselberg, and P. U. Jepsen, “Terahertz radar cross section measurements,” Opt. Express 18, 26399–26408 (2010). [CrossRef] [PubMed]
  9. B. Ferguson, S. Wang, D. Gray, D. Abbot, and X. C. Zhang, “T-ray computed tomography,” Opt. Lett. 27, 1312–1314 (2002). [CrossRef]
  10. E. Abraham, A. Younus, C. Aguerre, P. Desbarats, and P. Mounaix, “Refraction losses in terahertz computed tomography,” Opt. Commun. 283, 2050–2055 (2010). [CrossRef]
  11. S. Wang, B. Ferguson, D. Abbott, and X. C. Zhang, “T-ray imaging and tomography,” J. Biol. Phys. 29, 247–256 (2003). [CrossRef]
  12. S. Wang and X. C. Zhang, “Pulsed terahertz tomography,” J. Phys. D Appl. Phys. 37, R1–R36 (2004). [CrossRef]
  13. M. M. Awad and R. A. Cheville, “Transmission terahertz waveguide-based imaging below the diffraction limit,” Appl. Phys. Lett. 86, 221107 (2005). [CrossRef]
  14. X. Yin, B. W. H. Ng, B. Ferguson, and D. Abbott, “Wavelet based local tomographic image using terahertz techniques,” Digit. Signal Process. 19, 750–763 (2009). [CrossRef]
  15. 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]
  16. K. L. Nguyen, M. L. Johns, L. F. Gladden, C. H. Worral, P. Alexander, H. E. Beere, M. Pepper, D. A. Ritchie, J. Alton, S. Barbieri, and E. H. Linfield, “Three-dimensional imaging with a terahertz quantum cascade laser,” Opt. Express 14, 2123–2129 (2006). [CrossRef] [PubMed]
  17. 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. SPIE 7837, 783709 (2010).
  18. N. Sunaguchi, Y. Sasaki, N. Maikusa, M. Kawai, T. Yuasa, and C. Otani, “Depth-resolving terahertz imaging with tomosynthesis,” Opt. Express 17, 9558–9570 (2009). [CrossRef] [PubMed]
  19. T. Yasuda, T. Yasui, T. Araki, and E. Abraham, “Real-time two-dimensional terahertz tomography of moving objects,” Opt. Commun. 267, 128–136 (2006). [CrossRef]
  20. T. Yasui, K. Sawanaka, A. Ihara, E. Abraham, M. Hashimoto, and T. Araki, “Real-time terahertz color scanner for moving objects,” Opt. Express 16, 1208–1221 (2008). [CrossRef] [PubMed]
  21. G. T. Herman, Image Reconstruction From Projections: The Fundamentals of Computerized Tomography (Academic Press Inc., 1980).
  22. A. H. Andersen, and A. C. Kak, “Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm,” Ultrason. Imaging 6, 81–94 (1984). [CrossRef] [PubMed]
  23. L. A. Shepp, and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Trans. Med. Imaging 1, 113–122 (1982). [CrossRef] [PubMed]
  24. 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]
  25. C. Pradere, J.-P. Caumes, D. Balageas, S. Salort, E. Abraham, B. Chassagne, and J.-C. Batsale, “Photothermal converters for quantitative 2D and 3D real-time terahertz imaging,” Quant. Infrared Thermog. 7, 217–235 (2010). [CrossRef]
  26. 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).
  27. J. Radon, “¨Uber die Bestimmung von Funktionen durch ihre Integralwerte langs 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.
  28. 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]
  29. B. Recur, “Qualité et Précision en Reconstruction Tomographique: Algorithmes et Applications,” Ph.D. thesis, LaBRI, Bordeaux 1 University (2010).
  30. 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]

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: MOV (322 KB)     
» Media 2: MOV (285 KB)     
» Media 3: MOV (215 KB)     
» Media 4: MOV (631 KB)     
» Media 5: MOV (371 KB)     
» Media 6: MOV (354 KB)     

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited