Open Access
Add to BibSonomyAbstract
In this paper, a new fast terahertz reflection tomography is proposed using block-based compressed sensing. Since measuring the time-domain signal on two-dimensional grid requires excessive time, reducing measurement time is highly demanding in terahertz tomography. The proposed technique directly reduces the number of sampling points in the spatial domain without modulation or transformation of the signal. Compressed sensing in spatial domain suggests a block-based reconstruction, which substantially reduces computational time without degrading the image quality. An overlap-average method is proposed to remove the block artifact in the block-based compressed sensing. Fast terahertz reflection tomography using the block-based compressed sensing is demonstrated with an integrated circuit and parched anchovy as examples.
©2011 Optical Society of America
1. Introduction
Terahertz imaging and spectroscopy systems with terahertz electromagnetic wave (T-ray) have been successfully applied to various fields including security, inspection, and biomedical imaging [1–6]. T-ray biomedical imaging has great potential due to the interaction of T-ray with water molecules in the biological tissue noninvasively. Terahertz tomography is an emerging area for investigating small biological samples nondestructively [7,8].
T-ray imaging by raster scanning is the most widely used system. To reduce the long acquisition time by raster scanning, various techniques have been proposed. Array detectors have been developed to mitigate speed limitations, which, require expensive and complex equipment [9,10]. A fast T-ray imaging method based on compressed sensing (CS) has been proposed [11–14]; it acquires data in a much lower rate than the Nyquist sampling criterion by utilizing the image sparseness in general. As an extreme case, a T-ray imaging using a single detector without raster scanning was proposed in [12]. The T-ray beam passing the object is spatially modulated by a random pattern mask, and the resulting beam is focused by a lens and recorded by a single fixed detector. It is equivalent to recording a linear combination of the object beam multiplied by the random pattern mask. If a sufficient number of linear combinations of the object beams with different masks are acquired, a perfect reconstruction is possible using CS. Although the number of random pattern masks is much smaller than that of the pixels to be reconstructed, it is still quite a large number. Additional time and mechanism to change the number of masks are needed, which limit applications of the technique, especially for high resolution imaging.
T-ray scanning in spatial domain to measure the time domain spectroscopy (TDS) can be used for terahertz tomography as well as imaging. The block-based compressed sensing techniques have been proposed to reduce the computational burden of the reconstruction as well as measurement time [15,16]. If the object to be imaged is nonstationary, then the block-based CS is advantageous over conventional CS in minimizing L1 norm [11]. A computer simulation of the block-based CS has been done with two dimensional natural images, where the computational gain was addressed [15]. In this paper, we extend the technique to T-ray tomography. The block-based CS is similar to the conventional CS for the acquisition end, i.e., random sampling in a reduced sampling rate. Reconstruction is done in a block-by-block manner. Discontinuities at the boundaries of the block may be observed by the block-based CS. Wiener filtering with 3x3 window was performed in [15] to reduce such block artifact. We propose an overlap-average method, where the reconstruction blocks are overlapped and average value is taken for the pixels multiply reconstructed [16]. Note that no extra acquisition is needed. Wiener filtering is a simple approach to reduce the block artifact; however, some resolution degradation is unavoidable. On the other hand, the proposed ensemble average over the overlapped blocks needs not compromise the resolution, though it has larger computational load. Since the block-based CS can make reconstruction in parallel with the acquisition, computational load is not an issue if the reconstruction for a block is done faster than the acquisition of the next block.
The objective of this paper is to develop a fast tomography method that reduces measurement time by applying the compressed sensing technique to the spatial domain. The block-based reconstruction substantially reduces processing time without degrading the image quality. By applying the proposed method to the conventional T-ray time-domain spectroscopy system without any additional hardware, it was verified that tomographic images with qualities equivalent to those by the conventional full scan are obtained more quickly.
2. The block-based compressed sensing
Spectral density of an ordinary image including terahertz biomedical image and tomography tends to be sparse since most of energy is distributed in low spatial frequency band [17], which makes it possible to apply CS to reduce measurement time. Sparse sampling in the spatial domain for terahertz tomography may be represented as
where m is the measured vector of time signal, μ is the image to be reconstructed, and A is the sampling matrix having 1’s for measured indices. In Eq. (1) sparse sampling is applied in a transverse plane, and t denotes the time delay which is equivalent to depth at a given spatial location. Since compressed sensing is applied to the spatial domain, not to the time domain, further description of the algorithm is made for a given t without loss of generality. If the number of measured TDS is L and that of pixels for the image of size n x n is N (= n2), then m and μ are represented as column vectors of size Lx1 and Nx1, and A as a matrix of size LxN. Since L is usually much smaller than N, reconstruction finds a solution to the underdetermined Eq. (1). By utilizing sparseness of the image in a transformed domain, reconstruction may be done by minimizing ||Fμ||1 subject to ||m - Aμ||2 < ε, where F denotes a transform whose transformed coefficients appear sparse, e.g., Fourier transform or discrete cosine transform (DCT), || · ||p denotes Lp norm, and ε is a constant denoting tolerance or error bound for the measured data. Since the solution to minimize L1 norm is not given in a closed form, an iterative solution by FOCUSS algorithm is used [18,19].An analytical solution to Eq. (1) with a weighted L2 norm minimization is known. The FOCUSS algorithm is a recursive linear estimation based on a weighted pseudo-inverse solution. The weights at each step are derived from the transformed coefficients of the solution (estimated image) of the previous iterative step. Thus, a large transformed coefficient of the solution becomes larger as the number of iterations increases, which enforces the sparseness in the transformed domain. For a fast convergence, the weighting values of small coefficients are truncated with a threshold. The choice of initial image affects the convergence speed of the algorithm. Linear interpolation of the measured data is used as the initial estimate.
A simple example to show the rationale of the block-based compressed sensing is given in Fig. 1 . The signal in Fig. 1(a) has four distinct frequencies, 1 KHz, 3 KHz, 2 KHz, and 4 KHz in 4 intervals, respectively. The four frequencies are solely observed if a short-time Fourier transform is applied to each segmented signal as shown in (b), while many more frequencies are observed if Fourier transform is applied to the entire signal as shown in (c). As is seen in this example, if a signal is nonstationary, then the segmented (block-based) reconstruction is desirable in utilizing the signal sparseness. From the view point of sampled signal, the segmented signal is more easily approximated with fewer basis functions. In more realistic applications such as multidimensional imagery including tomography, the spatial correlation between pixels decreases as the distance increases. This property is fully utilized in image and video compression methods, such as JPEG and MPEG, where the image is divided into small blocks and each block is independently encoded [17].
Fig. 1 Nonstationary signal having four different frequencies, 1 KHz, 3 KHz, 2 KHz, and 4 KHz in 4 intervals (a). Four distinct spectra are observed if a segmented Fourier transform is applied to each segmented signal (b). Spectrum of the signal if Fourier transform is applied to the entire signal (c).
The main computational load of the reconstruction is the matrix inversion of size M x M if the number of sampling points with the compressed sensing is M. Thus, the complexity of the reconstruction is O(M 3) [20]. Let Q be the ratio of the number of pixels in the entire image to that in a block. Since the number of sampling points in a block is reduced by the same factor of Q, the complexity of the reconstruction for a block is reduced by the factor Q3. The number of blocks to be reconstructed is increased by the factor of Q if the blocks are laid without overlap. By considering an overlap factor K which will be described below, the computational gain using the block-based CS is Q3/(Q · K) = Q2 / K .
The optimal block size of the block-based CS depends on the local image characteristics. In general, a stationary region is better processed in a single block in order to fully utilize the correlation information, whereas an uncorrelated region is better divided into blocks for better reconstruction as well as lower computational load. In this paper, a fixed block size is used over the entire image, thus some discontinuities may be observed at the boundaries of the block if a region is not properly divided into blocks.
The overlap-average method is proposed to reduce such block artifact, where blocks are overlapped, and each block is reconstructed independently. Then, the average value is taken if a pixel is involved in multiple blocks due to the overlap. The overlap-average method does not increase the number of data measurements because the overlap is only for the reconstruction end defined on the existing acquired data. Therefore, the conventional CS and proposed block-based CS have the same measurement time. The overlap-average method alleviates the effect of inadequate division of regions by the fixed-sized blocks.
The overlap factor K is defined as the number of different reconstruction blocks in which a pixel belongs. If the overlap blocks are generated by shifting a fixed number of pixels (S) sequentially in both horizontal and vertical directions, the overlap factor is given by (B/S)2 for a block size of BxB except some pixels near the boundaries of the image. For instance if a block size of 16x16 is used for the reconstruction, and the overlap blocks are generated by shifting 2 pixels sequentially, the overlap factor is given by (16/2)2=64.
Sampling procedures are as follows. Once the block size and compression factor are chosen, then the scanning region is regularly divided into quarter-sized blocks without overlap. The number of sampling points in each quarter-sized block is kept constant over the entire image, and the locations of the sampling points are chosen randomly in each quarter-sized block. For instance, if the block size is 16x16 with a compression factor of 4, then the object is divided into 8x8 quarter blocks, and (8x8)/4=16 sampling locations are chosen randomly in each quarter block. The reason we set a quarter-sized block for random sampling is to make the sampling more uniformly to prevent some blocks, including those generated by the overlap, from having few sampling points by chance, where the block-based CS may not work properly.
The overlap-average method is tested by a computer simulation. A phantom composed of vertical, horizontal, and diagonal lines with variable widths and a simulation phantom used in computed tomography (CT) is devised. The size of the image is 256x256 (=65536). Reconstructed images by the block-based CS before and after applying the overlap-average method are shown in Fig. 2 in the first and second rows, respectively. Interpolated images by the cubic interpolation technique using the same data set are also shown in the third row for comparison. The compression factors are 8, 4, and 2 in the first, second, and third columns, respectively. The block size is 16x16, and blocks are overlapped by sequentially shifting 2 pixels in both horizontal and vertical directions (K=64). The block artifacts are observed in (a-c), especially in the CT phantom, while they are mostly removed by the overlap-average method as seen in (d-f). In this example, the computational gain by the block-based CS is 2563/(256*64) = 1024. Compared to the interpolated images (g-i), the reconstructed images by the block-based CS show higher resolution. For example, the vertical and horizontal lines are seriously degraded in the cubic interpolated images (g, h) when the compression factors are high, while they are mostly resolvable in the block-based CS images (d, e). A quantitative evaluation by the peak-signal-to-noise-ratio is summarized in Table 1 .
Fig. 2 Reconstructed images by the block-based CS before and after applying the overlap-average method are shown in the first and second rows, respectively. Cubic interpolated images are shown in the third row for comparison. The compression factors are 8, 4, and 2 in the first, second, and third columns, respectively.
Table 1. Evaluation of PSNR [dB] for the Images in Fig. 2 for Various Compression Factors by the Block-Based CS and Cubic Interpolation. Figure Numbers are Shown in Parenthesis.
3. T-ray tomography by block-based CS
A typical reflection-type terahertz time-domain spectroscopy system shown in Fig. 3 is used for reflection tomography. The T-ray is emitted from the terahertz transmitter of InAs wafer to the object; the reflected time domain signal of the T-ray is acquired by a single detector of photoconductive switch fabricated on low-temperature-grown GaAs (LT-GaAs). Translational and rotational scanning mechanisms are equipped with the system.
The time resolution of the T-ray is set to 70.3fs by the time delay stage, which is equivalent to the distance of 21.1μm in vacuum along the transmission (depth) direction, while the spatial resolution in the transverse plane is set to 250μm by the scanning mechanism. The number of sampling points in time domain is 512, which is used to resolve the depth information. The pulsed T-ray driven by a femtosecond Ti:sapphire laser has a full bandwidth of 0.36 THz centered at 0.55THz. Due to the different spectral response and chromatic aberration effects of T-rays with finite bandwidth, higher resolution may be obtained if a T-ray with narrow bandwidth is used. Instead of using a hardware filter, we applied a time-frequency analysis to obtain both temporal and spectral information for higher resolution tomography. Since the transmission distance of T-ray by the unit time delay (21μm divided by the refractive index inside the object) is 10~20 times higher than the spatial resolution (250μm), a number of planes may be used for time-frequency analysis without loss of resolution in general for isotropic resolution. Short-time Fourier transform is used for time-frequency analysis [21]. Since the depth information is directly obtained from the delay [7], the tomographic image is a two dimensional T-ray response of a particular frequency obtained from a series of two dimensional raster scanned images for a given range of time delays.
Multi-slice reflection tomographic images of an integrated circuit (IC) are shown in Fig. 4 . TDS is measured by conventional raster scanning. The number of sampling points is 120x240. The time taken for measuring one TDS (pixel) is 0.2s. Thus the total data acquisition time is 120×240×0.2=5760s, neglecting the time taken for moving the object in the raster scanning. Spatial resolution is 250μm in both horizontal and vertical directions. A two dimensional peak image is shown in (b). Four reconstructed planes are shown in (c-f), where different structures are clearly observed, while they are superimposed in the 2-D peak image (b). A rectangular time widow of 2.25ps (32 time delays) equivalent to spectral resolution of 0.444 THz is used for the short-time Fourier transform. The time resolution of 2.25ps is equivalent to a depth resolution of 375μm (slice thickness) inside IC. The spectral response at 0.67 THz appears most appropriate.
Fig. 4 Photograph of an integrated circuit (a), two dimensional terahertz peak image (b), and four tomographic images at the depths of 0.14mm, 2.06mm, 2.44mm, and 3.19mm are shown from 4 different time delays (c-f).
Applying the block-based compressed sensing to the terahertz tomography is shown in Fig. 5 , where only one sectional image is shown. Implementation of the block-based CS from the conventional TDS system is straightforward by programming the scanner to move predefined sampling positions instead of a raster scan. In this experiment, the TDS data previously measured by the raster scan are used according to the sampling positions of the block-based CS instead of separate acquisition for each compression factor, by which experimental variations or drifts during the measurements are excluded in the comparison. The block size is 8x8, and the overlap factor is 64; thus, the computational gain is about 3160. The block size is chosen by considering PSNR and computational time. PSNR appears not sensitive to the block size if the block size is equal to or larger than 8x8. The reconstructed image via full scan is shown in (a), while those obtained by the block-based CS are shown in (b)-(d) with the compression factors of 8, 4, and 2, respectively. As seen in Fig. 5, the reconstructed images with the compression factors of 4 and 2 (c, d) appear similar to that by full scan (a) with PSNRs above 30dB. Some degradation is observed in the reconstructed image with the compression factor of 8 (b).
Fig. 5 Terahertz tomographic images obtained with full scan (a), and by the block-based CS with compression factors of 8 (b), 4 (c), and 2 (d), respectively.
Computational load of the CS is much larger than that of raster scan. For the block-based CS, such computations can be done in parallel with the acquisition, thus little additional time is needed if the reconstruction of a block is done faster than the acquisition of the next block. This is true mostly due to a small block size with less computation and the technology advances in multi-core processors. This is different from CS where complete data is acquired before reconstruction. Thus imaging time can be reduced by the same factor of compression.
The PSNRs of the reconstructed images by the block-based CS and cubic interpolation are summarized in Table 2 . In the evaluation, the image by full scan is used as a reference image. As seen in Table 2, the block-based CS provides higher PSNRs compared to those by cubic interpolation, although the differences are not as high as those in Table 1. This is partially due to the experimental setup of the reflection terahertz tomography with an incident angle of 32°, which results in image blur as pointed out in reference [22].
Table 2. PSNR of the Reconstructed Images by the Cubic Interpolation and Block-Based CS for a Plane Shown in Fig. 5.
The preliminary results of terahertz tomography applied to bio-sample are shown in Fig. 6 . They are sectional images of parched anchovy. The number of sampling points in the spatial domain is 56x168. A photograph of the sample is shown in (a) with dissection image (b). Four sectional images are shown in (c)-(f) with the transverse resolution of 250μm and the slice thickness of 105μm. The short-time Fourier transform is applied for a rectangular window of size 8. Some internal organs are distinctively observed in the terahertz images (c-f) similar to those in the dissection image.
Fig. 6 Terahertz tomographic images of the parched anchovy. Photograph of the sample and dissection image are shown in (a) and (b), respectively. Terahertz sectional images are shown in (c)-(f) with transverse resolution of 250μm. Slice thickness is 105μm.
Applying the block-based CS to the terahertz tomography of parched anchovy is shown in Fig. 7 , where only two sectional planes are exemplary shown. The block size is 8x8, and overlap factor is 64; thus, the computational gain is about 330. Images reconstructed by full scan are shown in (a, e), and those by the block-based CS with compression factors of 8 (b, f), 4 (c, g), and 2 (d, h) are shown. As seen in Fig. 7, the reconstructed images with compression factors of 4 and 2 appear similar to those by full scan with PSNRs above 35dB. Some degradation is observed in the reconstructed images with the compression factor of 8. The PSNRs of the reconstructed images by the block-based CS and cubic interpolation are evaluated in Table 3 . In the evaluation, the images by full scan are used as reference. As seen in Table 3, the block-based CS provides higher PSNRs compared to the cubic interpolation. The PSNR differences become larger with higher compression factors.
Fig. 7 Terahertz tomographic images obtained by the block-based CS. Reconstructed images are shown by full scan (a, e), by the block-based CS with compression factors of 8 (b, f), 4 (c, g), and 2 (d, h).
Table 3. Average PSNR of the images by the cubic interpolation and block-based CS.
4. Conclusion
In this paper, a new fast tomography method with conventional T-ray imaging system is proposed using under-sampled data in the spatial domain by the block-based CS method. The proposed method does not require any additional hardware and achieves fast acquisition and reconstruction without image degradation. Therefore, it provides a general solution to the fast spatial-domain T-ray tomography.
By applying the proposed method to the conventional T-ray imaging system, tomographic images were obtained and their quality was analyzed. It was shown that 25% measurements can deliver acceptable quality images for the test objects. It was demonstrated that the block-based CS provides superior performances over the conventional cubic interpolation from both simulation and experiments. Tomographic images of integrated circuit and parched anchovy were shown to demonstrate the potential applications of terahertz tomography.
Acknowledgments
This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (No. 2009-0083512) and by a grant from the Korean Health Technology R&D Project of the Ministry for Health, Welfare & Family Affairs, Republic of Korea (A101954).
References and links
1. B. B. Hu and M. C. Nuss, “Imaging with terahertz waves,” Opt. Lett. 20(16), 1716–1718 (1995). [CrossRef] [PubMed]
2. D. M. Mittleman, M. Gupta, R. Neelamani, R. G. Baraniuk, J. V. Rudd, and M. Koch, “Recent advances in terahertz imaging,” Appl. Phys. B 68(6), 1085–1094 (1999). [CrossRef]
3. J.-H. Son, “Terahertz electromagnetic interactions with biological matter and their applications,” J. Appl. Phys. 105(10), 102033 (2009). [CrossRef]
4. A. J. Fitzgerald, V. P. Wallace, M. Jimenez-Linan, L. Bobrow, R. J. Pye, A. D. Purushotham, and D. D. Arnone, “Terahertz pulsed imaging of human breast tumors,” Radiology 239(2), 533–540 (2006). [CrossRef] [PubMed]
5. S. J. Oh, J. Kang, I. Maeng, J.-S. Suh, Y.-M. Huh, S. Haam, and J.-H. Son, “Nanoparticle-enabled terahertz imaging for cancer diagnosis,” Opt. Express 17(5), 3469–3475 (2009). [CrossRef] [PubMed]
6. S. J. Oh, J. Choi, I. Maeng, J. Y. Park, K. Lee, Y.-M. Huh, J.-S. Suh, S. Haam, and J.-H. Son, “Molecular imaging with terahertz waves,” Opt. Express 19(5), 4009–4016 (2011). [CrossRef] [PubMed]
7. D. M. Mittleman, S. Hunsche, L. Boivin, and M. C. Nuss, “T-ray tomography,” Opt. Lett. 22(12), 904–906 (1997). [CrossRef] [PubMed]
8. S. Wang and X.-C. Zhang, “Pulsed terahertz tomography,” J. Phys. D Appl. Phys. 37(4), R1–R36 (2004). [CrossRef]
9. Z. Jiang and X. C. Zhang, “Terahertz imaging via electrooptic effect,” IEEE Trans. Microw. Theory Tech. 47(12), 2644–2650 (1999). [CrossRef]
10. J. Xu and X. C. Zhang, “Terahertz wave reciprocal imaging,” Appl. Phys. Lett. 88(15), 151107 (2006). [CrossRef]
11. D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory 52(4), 1289–1306 (2006). [CrossRef]
12. W. L. Chan, K. Charan, D. Takhar, K. F. Kelly, R. G. Baraniuk, and D. M. Mittleman, “A single-pixel terahertz imaging system based on compressed sensing,” Appl. Phys. Lett. 93(12), 121105 (2008). [CrossRef]
13. Y. C. Shen, L. Gan, M. Stringer, A. Burnett, K. Tych, H. Shen, J. E. Cunningham, E. P. J. Parrott, J. A. Zeitler, L. F. Gladden, E. H. Linfield, and A. G. Davies, “Terahertz pulsed spectroscopic imaging using optimized binary masks,” Appl. Phys. Lett. 95(23), 231112 (2009). [CrossRef]
14. K. H. Jin, Y. Kim, D. S. Yee, O. K. Lee, and J. C. Ye, “Compressed sensing pulse-echo mode terahertz reflectance tomography,” Opt. Lett. 34(24), 3863–3865 (2009). [CrossRef] [PubMed]
15. L. Gan, “Block compressed sensing of natural images,” in Proc. Int. Conf. Digital Signal Processing, pp.403–406, Cardiff, UK (2007).
16. S. H. Cho, J.-H. Park, S.-H. Lee, H. Park, J.-H. Son, and C. B. Ahn, “Direct block-based compressed sensing for THz reflection tomography,” in Proc. Int. 2nd THz-Bio Workshop, pp.92–93, Seoul, Korea (2011).
17. T. Acharya and P.-S. Tsai, JPEG2000 standard for image compression: concepts, algorithms and VLSI architectures (Wiley-Interscience, 2005).
18. I. F. Gorodnitsky, J. S. George, and B. D. Rao, “Neuromagnetic source imaging with FOCUSS: a recursive weighted minimum norm algorithm,” Electroencephalogr. Clin. Neurophysiol. 95(4), 231–251 (1995). [CrossRef] [PubMed]
19. H. Jung, J. C. Ye, and E. Y. Kim, “Improved k-t BLAST and k-t SENSE using FOCUSS,” Phys. Med. Biol. 52(11), 3201–3226 (2007). [CrossRef] [PubMed]
20. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical recipes in C: the art of scientific computing 2nd ed. (Cambridge University Press, 1992).
21. S. Qian and D. Chen, Joint time-frequency analysis: methods and applications (Prentice Hall, 1996).
22. S.-H. Ding, Q. Li, R. Yao, and Q. Wang, “High-resolution terahertz reflective imaging and image restoration,” Appl. Opt. 49(36), 6834–6839 (2010). [CrossRef] [PubMed]
