In this study, we propose a THz computed tomography (CT) method based on phase contrast, which retrieves the phase shift information at each data point through a phase modulation technique using a Mach-Zehnder interferometer with a continuous wave (CW) source. The THz CT is based on first-generation CT, which acquires a set of projections by translational and rotational scans using a thin beam. From the phase-shift projections, we reconstruct a spatial distribution of refractive indices in a cross section of interest. We constructed a preliminary system using a highly coherent CW THz source with a frequency of 0.54 THz to prove the concept and performed an imaging experiment using phantoms to investigate its imaging features such as artifact-immune imaging, quantitative measurement, and selective detection.
© 2013 Optical Society of America
The terahertz (THz) wave is an electromagnetic wave with a frequency that lies between radio and infrared frequencies, ranging from 0.1 to 10 THz, corresponding to a wavelength of 3 mm to 30 μm. The recent developments in THz optical devices have accelerated the research on THz imaging [1–8]. A variety of imaging methods have been proposed to take advantage of the characteristics of THz waves. For instance, transillumination imaging, which is similar to X-ray imaging, can be realized based on the high transmissivity of THz waves for certain materials such as paper, wood, powder, semiconductors, and plastic. Consequently, this imaging can also be extended to the tomographic mode, analogous to X-ray computed tomography (CT) . Early THz tomographic imaging was based on the measurement of the time-of-flight of the reflected pulses . This technique reproduces the 3D refractive index profiles of objects consisting of well-separated layers of different refractive indices, which provides an extremely high depth resolution on the order of 1 μm. Subsequently, a multi-color THz-CT system using time-domain spectroscopy was developed [11–13]. In the experiments leading to this development, the tomographic information was reproduced based on the same CT reconstruction algorithm as that used in X-ray CT. This imaging system simultaneously forms cross-sectional images by absorption and phase contrasts in a sample with a complicated structure and a spatial resolution on the order of sub millimeters.
In spite of the excellent imaging properties of these systems, the lack of tractable sources and detectors hampered their practical application. For THz-wave generation, femtosecond lasers were employed in many cases; these lasers generate THz waves using a photoconductive antenna or nonlinear optical crystals. However, femtosecond lasers operate stably only under the conditions in which the temperature and humidity are controlled with precision. In addition, they are sensitive to mechanical or air vibrations. The recent advent of continuous wave (CW) sources has provided a suitable alternative to circumvent the issues related to femtosecond lasers. Although backward wave oscillator (BWO) has long been used as a CW source in THz regions, it has no satisfactory specifications with regard to life-time, stability, and output intensity. Recently, CW sources based on millimeter-wave frequency-multiplier circuits using semiconductor diodes have been developed for radio-astronomical use , which are immune to variations in temperature and vibrations and have stability in output intensity as well as a long life-time. Furthermore, until recently, bolometers cooled by liquid nitrogen or liquid helium have been used as THz-wave detectors; however, this technology was recently replaced by detectors based on Schottky-barrier diodes (SBD), which have been developed for commercial use and which operate stably even at ordinary temperatures without high voltage . A combination of such a CW source and SBD detector is expected to be used for the realization of a practically compact and portable THz imaging system .
On the other hand, from the viewpoint of CT imaging, the effects of reflection and refraction are substantial as compared to X-rays when the sample is bulky, because the refractive indices in THz regions are relatively large. Therefore, the mismatch in refractive indices causes unintended intensity dissipation in the incident direction, owing to refraction, reflection, and scattering, especially at the boundaries between air and material. As the unintended intensity dissipation is apparently regarded as attenuation at the boundary, the adoption of an imaging protocol similar to that of X-ray CT based on attenuation contrast results in the appearance of remarkable artifacts in the reconstructed image at the boundary, leading to impaired quantitative observations . In contrast, the phase-shift is not influenced by the intensity dissipation at the boundary, and hence, artifact-free reconstruction is expected as long as the transmitted wave can be detected.
In this study, we propose a THz-CT method based on phase contrast with significantly reduced artifacts using a highly coherent CW source and a semiconductor detector. The THz-CT is based on first-generation CT, which acquires a set of projections by translational and rotational scans with the use of a thin beam . At each data-point, the phase-shift is retrieved through a phase modulation technique using a Mach-Zehnder interferometer with a CW source. From the projections of the phase-shift, we reconstruct a spatial distribution of refractive indices in a cross section of interest. We constructed a preliminarily system to prove the concept and performed an imaging experiment using a polystyrene foam phantom to demonstrate the effectiveness of the proposed method.
2. Phase shift retrieval
Figure 1 shows a conceptual diagram of the proposed imaging system. In order to retrieve the phase shift from the transmitted intensity obtained from a single detector, we used a Mach-Zehnder interferometer. The thin parallel THz beam with high coherence from a CW source is split into the signal and local oscillator beams by the beam splitter BS1. The signal beam impinges on a sample, which is mounted on positioning devices so as to enable translation and rotation, and it is subjected to a variety of optical phenomena, such as refraction, reflection, scattering, and absorption, governed by inhomogeneous distribution of refractive indices in the sample. The transmitted beam retaining the propagation direction is mixed with the local oscillator beam, whose path length is adjusted using the reference mirror, at the beam splitter BS3. The mixed beam is detected by the SBD detector.
We provide below a quantitative discussion on how the phase-shift is retrieved from the measured interference intensity. First, we consider the case in which there is no sample. We represent the signal wave by , where k = 2π / λ, and AS and zS are the amplitude and path length of the signal arm, respectively. On the other hand, we denote the local oscillator wave by , where AL and zL are the amplitude and path length of the local oscillator arm, respectively. The output signal is proportional to the following value:19], must be added to the phase term in Eq. (2), we omit it because it is not related to the following discussion with respect to phase retrieval. On rearranging this expression, we obtainEq. (13) involves difference and division, it is apparently sensitive to measurement noise. We discuss the stability of estimation in Eq. (13). Generally, when the indirectly measured quantity ϕ derived from directly measured quantities Xi (i = 1,2,3,…,n), which are independent stochastic variables with means or true values mi and variances σi2, is represented as , the variance of ϕ, , is from the formula of error propagation given asEq. (13), where , , , and . Thus, we obtain the following expression after some calculations:Eqs. (9) to (12), , where IA0 and ϕ0 are true values of IA and ϕ, respectively. Therefore,Eq. (16) that σ12, σ22, σ32, and σ42 evenly contribute to . Therefore, Eq. (13) is a good-natured estimator of ϕ.
On the other hand, with the assumption that the refractive index n is unity in air, n – 1 is zero in the region outside the sample. Thus, from Eq. (8),
On the other hand, from Eq. (4), we find the value to be proportional to . Similar to the above calculation, assuming that the absorption term κ is zero in the region outside the sample, we obtain
3. Experimental setup
Figure 2 shows a schematic of the THz-CT imaging system based on the Mach-Zehnder interferometer. A frequency-multiplier CW source with a frequency of 540 GHz and a wavelength of 555 µm (0.7 mW, Virginia Diodes, Inc.) is used as the light source. The amplitude of the THz wave from the source is modulated at 30 kHz with a rectangular waveform generated by a function generator for lock-in detection. The beam splitters (a), (b), and (c), shown in Fig. 2, are made from high-resistivity silicon wafer. The THz wave is collimated to a thin parallel beam with convex lenses and parabolic mirrors. The beam size of the signal beam before convex lens (b) and the reference beam are 6 mm in diameter. Here, we need multiple mirrors and lenses for producing a thin parallel beam because the output wave from the source is divergent, leading to reduced incident intensity in front of a sample.
The beam after interference is detected by the SBD detector (WR-1.2, Virginia Diodes, Inc.). The detected signal was fed into a preamplifier (with a gain of 5000), then into a lock-in amplifier (time constant 30 ms), and finally to a personal computer through a data acquisition card. The dynamic range of this system is 65 dB.
The CT data-acquisition protocol is based on first-generation CT. Translating a sample from one edge to the other edge of the sample using a translational positioning device at a predefined step yields a single projection data at a projection angle. Then, the projection data acquisition procedure is repeated over 360° while rotating the sample using a rotational positioning device at a predefined angular step. The cross-sectional image is reconstructed from the set of projections using the filtered back projection (FBP) method . From Eq. (17), it can be seen that the reconstructed image corresponds to a distribution of refractive indices in the sample, because the projection data is a line integral of n – 1.
In the first-generation CT, the spot size of the beam should be held constant in the measurement range because the quality of the projection images depends on the uniformity of the beam cross section. Accordingly, we evaluated the beam profile using the knife-edge method. Figures 3(a) and 3(b) show the beam profiles across and along the beam, as a function of the knife position. In Fig. 3(a), the red dotted curve is that for measured data, and the blue solid curve is that for differential data obtained from the raw measured data. A beam diameter of less than 2.5 mm was realized in the 40 mm depth around the focus, where the focal length of the convex lens positioned just upstream of the object was 101.6 mm. The sample was measured within this region of the beam. Thus, the spatial resolution of the CT system is broadly 2.5 mm.
4. Imaging experiment
In order to validate the imaging characteristics of the system, we imaged a polystyrene foam sphere of 30 mm in diameter and having two channels of 5 mm in diameter. Figures 4(a) and 4(b) show the bird-view and top-view photographs of the phantom, respectively. Figure 4(c) is a schematic of the cross-section of the phantom at the level where the diameter of the cross-section is maximum. The sample was fixed in the translational and rotational stages such that the channels were perpendicular to the incident beam. After a series of 0.1-mm-step translational scanning procedures were completed, the sample was rotated in steps of 1°. A series of data acquisition procedures were performed over 360°, resulting in 360 projection data. Then, we performed a series of measurements by moving the reference mirror by λ/8 using the translational stage. We finally obtained four sets of projection data in a similar manner.
The raw sinograms obtained for the difference in path lengths between the signal and local oscillator arms of 0, λ/4, λ/2, and 3λ/4 are shown in Figs. 5(a)–5(d), respectively, where each image was 400 × 360 pixels in size. The brightness varies with the difference in path length. Figure 6(a) shows the sinogram image of the phase-shift estimated by Eq. (13) from the four sets of raw sinogram images shown in Fig. 5. Figure 6(b) shows a line profile of the first row of the sinogram, which is indicated by a red line in Fig. 6(a). We observe that a jump of 2π occurs at the boundaries between air and the sample. This is owing to the fact that the phase-shift is wrapped between –π and π because of the arctangent in Eq. (13). Therefore, the wrapped sinogram presents discontinuities at the boundaries. It is necessary to unwrap the sinogram to obtain projections from which the cross-sectional images are reconstructed. We applied the phase-unwrap algorithm devised by Cusack et al . Figure 7(a) is an unwrapped sinogram obtained from Fig. 6(a), and Fig. 7(b) is a line profile of the first row of Fig. 7(a). Through this algorithm, the wrapped phase-shift is satisfactorily unwrapped. Next, we reconstructed the cross section from the sinogram using the FBP algorithm with a Shepp-Logan filter . Figure 8(a) shows the phase-contrast THz-CT image of the phantom at the level, at which the cross-sectional diameter is maximum. For comparison, we reconstructed an absorption-contrast image (Fig. 8(b)) from the same four sets of the raw sinograms using Eq. (18). In both the reconstructed images, the pixel values were normalized in 8-bit gradation. Figure 8(c) compares the profiles of both these images along the red lines shown in Figs. 8(a) and 8(b), by plotting the pixel value as a function of the position. The blue and red curves in Fig. 8(c) correspond to the profiles of the phase- and attenuation-contrast CT images, respectively. In the attenuation-contrast CT, remarkable artifacts are observed at the boundary between air and polystyrene foam regions. In addition, the pixel value in some regions is much lower than that in the air region, although the pixel value in air is zero. These facts indicate that the attenuation-contrast CT can enable neither artifact-free reconstruction nor quantitative measurement. On the other hand, no remarkable artifacts are observed in the phase-contrast CT image.
As discussed in Section 2, pixel values in the reconstructed phase-contrast CT image represent the refractive indices. Using Eq. (17), we recalculated an average pixel value denoted by the yellow rectangle (120 × 60 pixels) indicated in Fig. 8(a). The estimated refractive index in polystyrene foam ranges from 1.02 to 1.035, and the average and standard deviation are 1.026 and 0.003, respectively. In reference , it was reported that the refractive index of polystyrene foam ranges from 1.016 to 1.022. The values are similar to our estimated average value.
Although the estimated refractive indices in air regions should be almost zero, those in the channel regions take finite values as shown in Fig. 8(c). This is because the diameter of the channels is relatively large as compared to the diameter of the beam cross section. If the beam diameter is sufficiently less than the channel diameter, the nonconformity is dissolved.
5. Selective detection from scattered wave
The proposed imaging system is based on the Mach-Zehnder interferometer. The reasons behind the use of the interferometer are the ability to measure phase-shift information using a CW source with a high coherence and the ability to selectively detect the forward-scattered directional component from multiply and divergently scattered waves emerging from the sample through the mixing of the signal beam with the local oscillator beam. The ability of selective detection of the forward-scattered waves was experimentally investigated using a sample to induce multiple scattering. Figures 9(a) and 9(b) show a photograph of the sample and its schematic, respectively. The sample consists of 0.5-mm-diameter ceramic balls made of Zirconia, which are randomly arranged on a sheet of hard paper and enable Mie scattering. In order to monitor the beams, the detector system with the lens L4 (Fig. 1) set on the positioning devices was horizontally and vertically raster-scanned in the plane perpendicular to the optical axis. The scanning range was 60 × 60 mm2 with a center corresponding with the cross point between the optical axis and the scanning plane, and both the horizontal and vertical scanning steps were of 0.5 mm.
First, we monitored the signal beam without the sample by stopping the local oscillator with a beam stopper located between the beam splitters BS2 and BS3 (Fig. 1). Figure 10(a) shows a distribution of the beam intensities in 3D representation viewed from top with colored contour regions. Figure 10(b) shows a line profile along the red line indicated in Fig. 10(a).
Next, we monitored the signal beam by inserting the scattering sample between the lenses L2 and L3. Figures 10(c) and 10(d) show the distribution of beam intensity and the line profile, respectively. From Fig. 10(c), we observe that the distribution appears rough, owing to multiple scattering, and the shape of distribution remarkably changes from that in Fig. 10(a). Further, two peaks are generated, as shown in Fig. 10(d).
Subsequently, we monitored the local oscillator beam by stopping the signal beam with a beam stopper located between the mirror and BS3. Figures 10(e) and 10(f) show the distribution of the beam intensity and the line profile, respectively. Here, in order to clearly demonstrate the effect of selective detection, we intentionally deformed the local oscillator beam profile. As shown in Fig. 10(f), the smaller peaks are observed in the right tail of the main peak.
Finally, we monitored the mixed beam with the sample. Figures 10(g) and 10(h) show the distribution of beam intensity and the line profile, respectively. Comparing Figs. 10(e) and 10(g), we can observe that the distribution of the local oscillator beam resembles that of the mixed beam, while the distribution of the mixed beam appears relatively smooth. Similarly, from Figs. 10(f) and 10(h), the profile of the local oscillator beam resembles that of the mixed beam, while the intensity of the mixed beam is slightly lower than that of the local oscillator beam. The normalized cross-correlation between Figs. 10(e) and 10(g) is 0.821. Although the shape of profile slightly varies by changing the difference in path length between signal and local arms, the normalized cross-correlation between profiles of reference beam and mixed beam falls within the range from 0.75 to 0.85, and the range in the profile of mixed beam where the intensity is not zero is invariant. Thus, the shape of mixed beam greatly receives restriction of the shape of reference beam.
Next, we measured the mixed signals while gradually changing the difference in path length between signal and local arms by translating the reference mirror at a step of 2 μm at the five detector positions, P1 to P5, which are represented as red x-marks on the horizontal lines in Figs. 10(d), 10(f), and 10(h). Figure 11 shows the mixed signals, where the horizontal and vertical lines are the deference in path length and detected intensity, respectively. We can see that signals oscillate in a sinusoidal manner at P2, P3, and P4, where the intensities of mixed beam are not zero, while no signals are observed at P1 and P5, where the intensities of mixed beam are almost zero. We quantitatively consider the signal at P2. From Figs. 10(d) and 10(f), intensities of reference and signal beams are 0.350 mV and 0.002 mV, respectively. Using these values, ID = 0.352 mV and IA = 0.053 mV from Eqs. (6) and (7), respectively. On the other hand, from Fig. 11, the maximum (M) and minimum (m) of signal at P2 are 0.392 mV and 0.306 mV, respectively. Assuming that the signal is a sinusoidal function while it is distorted because the sample has inhomogeneous regions in refractive index, the center of oscillation and amplitude are obtained as 0.349 mV and 0.0043 mV from (M + m) / 2 and (M – m) / 2, which are regarded as the estimated values of ID and IA, respectively. In spite of the rough assumption that the signal is a sinusoidal function, these values are close to the above values obtained from Eqs. (6) and (7) using experimental data from Fig. 10. Also for signals at P3 and P4, similar results are derived as shown in Table 1. From the results, we can see that interference is observed in the region where the intensity of mixed beam is not zero. It addition, the amplitude of mixed signal depends on the intensity of signal beam, and the center of oscillation of mixed signal depends on the intensity of reference beam. Therefore, we can conclude that the interferometer selectively detects only the forward-scattered component, propagating in the same direction as that of the local oscillator beam and preserving the phase information, from among multiply and divergently scattered signal waves emerging from the sample when the signal beam is mixed with the local oscillator beam. This is because the imaging system could quantitatively reproduce the cross section with no remarkable artifacts in spite of the mismatch in the refractive index, as described in Section 4.
Furthermore, the results verify the significant potential as follows: In this research, the imaging system is based on the first-generation CT. This leads to a long data acquisition time, which hampers its practical application. On the other hand, if volumetric parallel beams are available for use as an incident beam and a local oscillator beam, we can obtain a single projection at a time using a 2D detector. Therefore, in order to collect a set of projections required for CT reconstruction, we simply rotate the sample without translation, which results in a drastic reduction in the total data acquisition time. In such an imaging system with a 2D detector, the problem of cross talk is expected to occur. However, the results described above indicate that, if each flux in the volumetric local oscillator beam is parallel to other fluxes, the flux can selectively detect only the forward-scattered component propagating in the same direction as that of the flux from a divergent scattered signal beam. Therefore, we will be possible to obtain a single projection at a time without cross-talk.
We proposed a novel THz-CT imaging method based on phase-contrast using a CW source. We constructed a preliminary system for validating the concept. The system acquires the projections of phase-shift using a phase modulation technique with the Mach-Zehnder interferometer and reconstructs the distribution of refractive indices from the measured projections. The experiments using a physical phantom showed effectiveness in that artifact-free quantitative reconstruction was feasible.
The present system has two major problems: One is that the incident intensity in front of a sample is insufficient owing to the use of multiple mirrors and lenses for forming a thin parallel beam. The other is that the total data-acquisition time is very long because the system is based on first-generation CT, in which a set of projection data must be sequentially collected point by point using a thin beam. However, the excellent property of selective detection, which was proved in this study, will surpasses these difficulties. The use of volumetric parallel beams as both signal and local oscillator beams enable us to obtain a single projection at a time by using a 2D detector. Thereby, the data-acquisition time can be drastically reduced. Further, only a few mirrors and lenses are required to form volumetric parallel beams. Our future work will involve construction of a practical system using a volumetric parallel beam and a 2D detector.
References and links
2. R. M. Woodward, V. P. Wallace, D. D. Arnone, E. H. Linfield, and M. Pepper, “Terahertz pulsed imaging of skin cancer in the time and frequency domain,” J. Biol. Phys. 29(2/3), 257–259 (2003). [CrossRef] [PubMed]
4. T. Yasuda, T. Iwata, T. Araki, and T. Yasui, “Improvement of minimum paint film thickness for THz paint meters by multiple-regression analysis,” Appl. Opt. 46(30), 7518–7526 (2007). [CrossRef] [PubMed]
5. T. Kiwa, J. Kondo, S. Oka, I. Kawayama, H. Yamada, M. Tonouchi, and K. Tsukada, “Chemical sensing plate with a laser-terahertz monitoring system,” Appl. Opt. 47(18), 3324–3327 (2008). [CrossRef] [PubMed]
6. S. R. Murrill, E. L. Jacobs, S. K. Moyer, C. E. Halford, S. T. Griffin, F. C. De Lucia, D. T. Petkie, and C. C. Franck, “Terahertz imaging system performance model for concealed-weapon identification,” Appl. Opt. 47(9), 1286–1297 (2008). [CrossRef] [PubMed]
7. Y. Kawada, T. Yasuda, H. Takahashi, and S.-i. Aoshima, “Real-time measurement of temporal waveforms of a terahertz pulse using a probe pulse with a tilted pulse front,” Opt. Lett. 33(2), 180–182 (2008). [CrossRef] [PubMed]
8. Y. Kawada, T. Yasuda, H. Takahashi, and S. Aoshima, “Real-time measurement of temporal waveforms of a terahertz pulse using a probe pulse with a tilted pulse front,” Opt. Lett. 33(2), 180–182 (2008). [CrossRef] [PubMed]
9. C. Kak and M. Slanery, “Principles of Computerized Tomographic Imaging,” New York: IEEE Press (1987).
12. S. Wang, B. Ferguson, and X.-C. Zhang, “Pulsed terahertz tomography,” J. Phys. D Appl. Phys. 37(4), R1–R36 (2004). [CrossRef]
13. E. Abraham, A. Younus, C. Aguerre, P. Desbarats, and P. Mounaix, “Refraction losses in terahertz computed tomography,” Opt. Commun. 283(10), 2050–2055 (2010). [CrossRef]
14. D. Porterfield, J. Hesler, T. Crowe, W. Bishop, and D. Woolard, “Integrated terahertz transmit / receive modules,” Proc. of 33rd European Microwave Conference, 1319–1322 (2003).
15. A. Dobroiu, M. Yamashita, Y. N. Ohshima, Y. Morita, C. Otani, and K. Kawase, “Terahertz imaging system based on a backward-wave oscillator,” Appl. Opt. 43(30), 5637–5646 (2004). [CrossRef] [PubMed]
16. 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(6), 5105–5117 (2011). [CrossRef] [PubMed]
18. J. Hsieh, Computed Tomography Principles, Design, Artifacts, and Recent Advances, Second Edition (John Wiley & Sons, Inc. & SPIE, 2009).
20. A. Rosenfeld and C. Kak, Digital Picture Processing, 2nd Ed., Vol. I (Academic Press, 1982).
22. G. Zhao, M. Mors, T. Wenckebach, and P. C. M. Planken, “Terahertz dielectric properties of polystyrene foam,” J. Opt. Soc. Am. B 19(6), 1476–1479 (2002). [CrossRef]