Open Access
Add to BibSonomyAbstract
This paper proposes a method to enhance terahertz reflection tomographic imaging by interference cancellation between layers. When the gap between layers is small, the signal reflected on the upper layer interferes with that on the lower layer, which degrades the quality of the reconstructed tomographic image in the lower layer. The proposed method estimates the upper-layer reflection signal by system modeling, which is then eliminated from the acquired signal. In this way, it can provide the correct lower-layer reflection signal, thereby improving the quality of the lower-layer tomographic image. The performance of the proposed method was confirmed using computer simulation data and real terahertz reflection data.
© 2016 Optical Society of America
1. Introduction
Tomography imaging systems using terahertz (THz) are widely employed in the fields of medical imaging, defect detection, and security [1–7]. The THz reflection tomography system inputs a THz pulse to an object and acquires the reflected signal. The acquired signal contains the results of reflection from all layers of the object and each reflection appears at a different time region due to different layer depths. After analyzing the acquired signal in the time region corresponding to the layer of interest, the tomographic image of the layer is reconstructed [8,9].
Figure 1 shows the operation of a reflection tomography system that inputs a THz pulse p(t) to an object with two layers, L1 and L2, and acquires a reflected signal x(t), where t is a sampled time index. The acquired signal x(t) is the sum of two reflected signals, u1(t) and u2(t), on L1 and L2, respectively. To capture the true characteristics of each layer, u1(t) and u2(t) need to be separated from their sum x(t) = u1(t) + u2(t). u1(t), u2(t) and x(t) vary with pixel position, but the position index is omitted for the sake of simple expression in this paper.
Fig. 1 Acquisition of reflected signal in THz reflection tomography system for object with two layers.
In order to investigate the nature of the THz reflected signal, x(t) was obtained from a THz reflection tomography system [8,9]. If strong reflection occurs on L1, u2(t) ≈0 and x(t) ≈u1(t) because of the limited penetration through L1. Figure 2 shows x(t) ≈u1(t) obtained in this manner. u1(t) has a decaying region after its positive peak position. This is a common property of THz reflected signals, regardless of the reflection coefficient, when generated by the same input THz pulse. For a two-layer object, therefore, if u1(t) does not decay completely to zero before the reflection on L2 begins, which happens when a gap between layers is small, x(t) in a time region corresponding to L2 is not equal to u2(t) because the residue of u1(t) is added to u2(t). That is, due to the dispersive nature of the THz reflected signal, u1(t) may interfere with u2(t) when acquiring x(t). In this case, the true characteristics of L2 cannot be extracted from x(t), and the quality of the L2 tomographic image decreases.
To solve this problem, in this paper, a method for canceling the interference between layers in x(t) is proposed. The key operations of the proposed method consist of three steps: modeling of reflection on L1, estimation of u1(t), and estimation of u2(t) from x(t) as well as the estimated u1(t). The tomographic image of L2 is reconstructed from the estimated u2(t), and the quality will be better than that obtained through the conventional method using the original x(t). The reflection modeling is based on system modeling and adaptive filter techniques. The proposed method can be applied to an existing THz reflection tomography system without additional hardware or modification, because it only changes the procedure to reconstruct the image from the given acquired reflection data.
2. Proposed interference canceling method
2.1 Overall procedure
When x(t) in all pixel positions of the object is given, a time region for the layer of interest is defined and the tomographic image for the layer is reconstructed from x(t) in this time region using the conventional method. In this paper, the image pixel values are defined by the positive peak values of x(t) in the given time region [8,9]. In this way, the meaningful layers, Li, and the time region for each of those layers are defined.
In order to apply the proposed method to x(t), the input THz pulse p(t) is required. In general, however, p(t) is not measured separately. Hence, p(t) is estimated from the given x(t). In the proposed method, x(t) from the position with the strongest reflection on L1 is chosen as the estimated p(t). Since this x(t) coincides with a signal reflected on the first layer, it receives no interference from the upper layers. In addition, it contains little modification because of the strength of the reflection. Based on this argument, the estimated p(t) will be very close to the true input.
The proposed method is applied to x(t) in each pixel position and independently determines each pixel value for the image. Using x(t) and the estimated p(t), a system that models the reflection on L1 is determined through system modeling and adaptive filter techniques. Then, u1(t) is estimated by simulating the reflection on L1 using the determined system. u2(t) is estimated with u2′(t) = x(t) − u1′(t), where u1′(t) is the estimated u1(t). Finally, the positive peak value of u2′(t) in the time region for L2 is determined and the tomographic image for L2 is reconstructed. The same operation can be applied to objects with more than two layers in a repetitive manner, which will be explained in Section 3.2.
2.2 System modeling for reflection
In this paper, reflections are modeled by a LTI (linear time-invariant) system. For simple computation, an FIR (finite impulse response) structure is used, and the reflected signal is expressed by a linear combination of time-shifted inputs p(t − k), k = 0, …, K, where K is the order of the system. When the impulse response of the system is h(t), the reflected signal becomes x(t) = p(t) * h(t) = , where * is a convolution operation.
In general, the modeling of unknown system is implemented as in Fig. 3(a). An input p(t) is applied to the unknown system and its output x(t) is computed. The same p(t) is also applied to a varying system g(t) and the output y(t) is computed. Then, the optimal gopt(t) is determined such that the mean-square-error between x(t) and y(t) is minimized. Since the two systems conduct equivalent operations as far as the relation between input and output is concerned, gopt(t) is now modeling the unknown system.
Fig. 3 Block diagram for system modeling. (a) Conventional block diagram for general system modeling. (b) Block diagram for modeling reflection on L1 in the proposed method.
In order to estimate u1(t), a system that models the reflection on L1 is required. If the system representing the signal modification by the reflection on L1 and L2 is h1(t) and h2(t), respectively, then the reflected signal becomes u1(t) = p(t) * h1(t) and u2(t) = p(t − Δ) * h2(t) = p(t) * h2(t − Δ), where Δ is the two-way time delay between L1 and L2. Here, the signal modification by the penetration through layers is ignored. Hence, a block diagram can be designed for modeling the reflection as shown in Fig. 3(b). However, this cannot correctly model the reflection on L1, because the resulting gopt(t) models [h1(t) + h2(t − Δ)], not h1(t) alone.
In order to model h1(t) alone using Fig. 3(b), the gap between two layers is utilized. In x(t), the reflection on L2 appears after the reflection on L1 by Δ. The start time t = Ri for ui(t), i = 1, 2, is first determined so that ui(t) = 0 for 0 ≤ t < Ri. To do so, the shape of p(t) is investigated and a time gap Q between the start position and the positive peak position of p(t) is determined as shown in Fig. 4(a), where p(t) was obtained by the proposed method from x(t). If only moderate signal change occurs by reflection, it can be assumed that the reflected signal also has the same time gap of Q between the start and the positive peak positions. Therefore, when x(t) is given, Ri can be set to (positive peak position of x(t) in Li region − Q) as shown in Fig. 4(b), where the Li region represents a time region on the t-axis corresponding to layer Li. Consequently, in R1 ≤ t < R2, u2(t) = 0 and x(t) is generated using only the h1(t), resulting in x(t) = u1(t) = p(t) * h1(t). Based on this property, if gopt(t) is determined such that E(g) in Eq. (1) is minimized with a search region R1 ≤ t < R2, then gopt(t) models the reflection on L1.
Fig. 4 Determination of Ri for reflection modeling. (a) Input THz pulse p(t) and the time gap Q between the start and the positive peak positions. (b) R1 and R2 determined from x(t).
2.3 Estimation of reflected signal
Since gopt(t) models the reflection on L1, u1′(t) can be estimated by simulating the reflection using p(t) and gopt(t), resulting in u1′(t) = p(t) * gopt(t). However, this is valid only for R1 ≤ t < R2, because such gopt(t) was optimized in R1 ≤ t < R2 as seen in Eq. (1). Therefore, it is necessary to determine u1′(t) for t ≥ R2 in addition.
It is apparent that the actual u1(t) is generated by the fixed operation of h1(t) during all t. Therefore, gopt(t), despite being optimized for R1 ≤ t < R2, still models the reflection on L1 for t ≥ R2; the only limitation in this process is that the property of u1(t) in t ≥ R2 is not considered when determining gopt(t). Consequently, for all t, even when u2(t) ≠ 0, u1(t) can be estimated as u1′(t) = p(t) * gopt(t). The key concept for the proposed estimation method is to model the system in a possible region (R1 ≤ t < R2) and to copy the model to the other region (t ≥ R2) where system modeling is impossible due to nonzero u2(t) in x(t). Finally, after eliminating u1′(t) from x(t), the estimated u2(t) is determined by u2′(t) = x(t) − u1′(t) = x(t) − p(t) * gopt(t). The tomographic image of L2 is then reconstructed from u2′(t), not from x(t).
2.4 Gradient-based optimization
In order to implement the proposed method, it is necessary to determine gopt(t), which minimizes E(g) in Eq. (1). When arbitrary x(t) and p(t) are given, however, this problem cannot be analytically solved, and a gradient-based adaptive algorithm is used for this purpose.
Since E(g) is a quadratic function of g(t), it has one global minimum point. Hence, starting from an arbitrary initial g(t), gopt(t) can be reached if g(t) is repeatedly adjusted such that E(g) decreases. To do this, a direction of adjustment is necessary. According to the adaptive algorithm theory, the negative gradient of E(g) is the direction to gopt(t) [10]. The gradient of E(g) is given in Eq. (2). Hence, the update equation for g(t) in Eq. (3) continuously decreases E(g), as long as an appropriate update constant μ is used [10].
Therefore, in the proposed method, gopt(t) is determined by running the update of g(t) iteratively until convergence when a decrease of E(g) by the new g(t) is very small. μ is selected empirically after investigating the behavior of E(g) reduction as the update continues. The computational load of the proposed method depends on the number of iterations until convergence. Many techniques for fast convergence are reported such as signal whitening [10], and they can be utilized for fast operation of the proposed method.
3. Performance analysis
3.1 Computer simulation data
The performance of the proposed method is analytically measured using computer simulation data. The input signal p(t) is obtained from the THz reflection system as before, which is shown in Fig. 4(a), where Q = 20. The computer simulation data are generated assuming a virtual two-layer-object as shown in Fig. 5(a). The gap between L1 and L2 is set to 15 samples on the t-axis, and the number of pixels is 60 × 65 = 3900.
Fig. 5 Computer simulation data used in the performance analysis. (a) Test object used to generate computer simulation data. (b) System impulse response h0(t) used to model the reflection. (c) Simulated reflection signal on each layer, u1(t) and u2(t). (d) Acquired signal, x(t) = u1(t) + u2(t). (e) Estimated reflection signal, u1′(t) and u2′(t), by the proposed method.
The reflection is modeled by a system impulse response h0(t) given in Fig. 5(b), where the system order is K = 50. It was verified that h0(t) outputs a signal that resembles the actual reflected signal. In addition, in order to consider the effect of minute variation in the reflection property for different positions and the measurement noise, a random noise n(t) is added to h0(t), where the power of n(t) is set to 1% and 10% of h0(t) power. The reflection on L1 is then modeled by h1(t) = h0(t) + n(t). The center square area, S, in L2 has a strong reflection and is modeled by h2(t) = 2[h0(t) + n(t)], and the background B has a weak reflection and is modeled by h2(t) = h0(t) + n(t). Figure 5(c) shows an example of simulated reflection signal in the S area of L2, where strong reflection on L2 occurs. In the L2 region on the t-axis depicted in the plots, u1(t) ≠ 0, which confirms that u1(t) interferes with u2(t). So x(t) = u1(t) + u2(t) in the L2 region has a peak value that is different from the true peak for u2(t). The proposed method is applied to x(t) = u1(t) + u2(t) in Fig. 5(d), and determines gopt(t) over the search region R1 ≤ t < R2. After 500 iterations in the g(t) update process in Eq. (3), the performance of system modeling almost saturates. The resulting u1′(t) = p(t) * gopt(t) and u2′(t) = x(t) − u1′(t) are shown in Fig. 5(e), and both are very close to their true signals shown in Fig. 5(c). The tomographic image of L2 is reconstructed from the peak value of u2′(t) in the L2 region.
Figure 6 shows the tomographic image of L2 reconstructed by the conventional method using x(t) and by the proposed method using u2′(t), along with the error image between the true and the reconstructed images, where the true image is determined from u2(t). It can be seen that the proposed method yields an image of higher quality.
Fig. 6 Tomographic image of L2 reconstructed by the conventional method and by the proposed method when the layer gap is 15 samples. (a) 1% noise addition to h0(t). (b) 10% noise addition to h0(t).
The analytical performance of the proposed method is measured in terms of PSNR (peak signal-to-noise ratio) between the true image and the reconstructed image for the layer gaps between 13 and 22 samples. Figure 7 shows PSNR for the conventional method and the proposed method when 1% and 10% noise is added to h0(t), along with the PSNR increase by the proposed method. As the layer gap increases, the PSNR increase by the proposed method becomes smaller because the effect of layer interference is decreasing.
Fig. 7 Performance of the proposed method measured in terms of PSNR(dB) for computer simulation data. (a) 1% noise addition to h0(t). (b) 10% noise addition to h0(t).
From these evaluation results using simulation data that resembles the actual reflection signal, it is verified that the interference of L1 is effectively cancelled from x(t) and that the resulting tomographic image of L2 has a higher quality than the image without interference cancellation.
3.2 Real THz reflection data
Experiments on THz reflection tomography are conducted and real THz reflection data are collected in real environments. In order to analytically measure the performance of the proposed method, a test object with four layers consisting of flat materials with uniform characteristics is used, as shown in Fig. 8(a). Each layer has a specific shape to be imaged. A THz pulse p(t) is transmitted to the test object and the reflected signal x(t) is acquired.
Fig. 8 Performance analysis for real THz reflection data. (a) Test object used to acquire the reflected signal. (b) x(t) in position with strong reflection on L2 . (c) x(t) in position with strong reflection on L3. (d) Estimated reflection signal, u1′(t) and u2′(t), for x(t) in (b) by the proposed method.
As the first step in the proposed method, p(t) is determined from x(t). On L1, the strongest reflection occurs in the square-shape area, so x(t) in the lower right position is chosen as p(t), and this p(t) was shown in Fig. 4(a). The interference cancellation for L2 is performed in the following procedure. Figure 8(b) shows an example of x(t) in a position with strong reflection on L2, where two distinct positive peaks represent the reflections on L1 and L2. Since Q = 20 for p(t), Ri, i = 1, 2, is set to (positive peak position of x(t) in Li region − 20), as shown in Fig. 8(b). Inserting these Ri's into Eq. (3), gopt(t) for each pixel position is determined and the image is reconstructed from u2′(t) = x(t) − p(t) * gopt(t).
Interference cancellation for L3 and L4 is different from that for L2 in the search region of gopt(t). Figure 8(c) shows an example of x(t) in a position with strong reflection on L3, where x(t) = u1(t) + u2(t) + u3(t). u1(t) and u2(t), the reflected signal on L1 and L2, may interfere with u3(t), the reflected signal on L3. But, the shape of p(t) indicates that the layer interference does not proceed as far as two layers. For L3, therefore, it is necessary to cancel only the interference from L2, which requires system modeling for only the reflection on L2. Subsequently, the search region for gopt(t) in Eq. (3) is changed to R2 ≤ t < R3 as shown in Fig. 8(c), where R3 is set to (positive peak position of x(t) in L3 region − 20). With this search region, gopt(t) is determined and u3′(t) is set to x(t) − p(t) * gopt(t) as before. Finally, the tomographic image of L3 is reconstructed from u3′(t). Similarly, for L4, gopt(t) is determined using the search region R3 ≤ t < R4, where R4 is set to (positive peak position of x(t) in L4 region − 20), and the tomographic image of L4 is reconstructed from u4′(t) = x(t) − p(t) * gopt(t).
Figure 8(d) shows the results of u1′(t) and u2′(t) for x(t) in Fig. 8(b). The estimated u1′(t) is not zero in the L2 region, which confirms that interference actually occurs between L1 and L2. Figure 9 shows the tomographic image of each layer reconstructed by the conventional method and by the proposed method. In this experiment, however, the true u1(t) and u2(t) and the true image are not provided, hence the analytical quality cannot be measured either. Therefore, in this paper, another method is used to verify the effect of interference cancellation and the enhancement of image quality.
Fig. 9 Tomographic image reconstructed by the conventional method and by the proposed method for real THz reflection data. (a) L2. (b) L3. (c) L4.
All materials used in the experiments are flat with uniform characteristics. Therefore, in an ideal case, all positions in each shape area and the background for each layer should have the same x(t) and the same pixel values. Based on this argument, a virtual true image is defined for the performance analysis, which has area-wise constant pixel values. This constant pixel value, however, is not known in the experiments. Hence, the best guess for the constant pixel value is the area-wise average for all acquired pixel values. Therefore, from the tomographic image reconstructed using the conventional method, the virtual true image is defined by averaging the pixel values for each area. The PSNR between the virtual true image and the reconstructed image is then measured.
Table 1 shows the PSNR for the conventional method and the proposed method. For all layers, the proposed method has a higher PSNR than the conventional method, which implies that the proposed method cancels the interference between layers and outputs an image closer to the virtual true image. From these results, the performance of the proposed method for real THz reflection data is confirmed.
Table 1. Performance of the Proposed Method Measured in Terms of PSNR(dB) for Real THz Reflection Data
As Table 1 shows, the increase in the PSNR for real THz reflection data is not large. The layer gap of 0.6mm corresponds approximately to 22 samples in the sampled time domain. Therefore, the reason for small PSNR increase is that the layer gap in the test object is large and the effect of layer interference is small. Hence, even when the layer interference is cancelled out correctly by the proposed method, the PSNR increase may not be large, which can be confirmed by the results in Fig. 7. Another experiment using a test object with a smaller layer gap and severe interference is under investigation.
4. Conclusion
In this paper, a method to enhance THz reflection tomographic imaging by interference cancellation between layers is proposed. From the acquired reflection data, the proposed method first determines a system that models the reflection on the upper layer, then estimates the upper-layer reflected signal by simulating the reflection using the modeled system. Afterward, it eliminates the estimated upper-layer interference signal from the acquired signal, resulting in an interference-canceled signal. By reconstructing the tomographic image using the interference-canceled signal, instead of the original reflection signal, an improvement in image quality is achieved. The proposed method was applied to computer simulation data and real THz reflection data, and enhancement of the image quality by the proposed method, compared to a method without interference cancellation, was confirmed. The proposed method can be applied to an existing THz reflection tomography system without additional hardware or modification, because it only changes the procedure to reconstruct the image from the acquired reflection data
Acknowledgments
This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (No. 2009-0083512 and 2015R1A5A1A95022841), and by the Research Grant of Kwangwoon University in 2016.
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, R. H. Jacobsen, and M. C. Nuss, “T-ray imaging,” IEEE J. Sel. Top. Quantum Electron. 2(3), 679–692 (1996). [CrossRef]
3. D. M. Mittleman, S. Hunsche, L. Boivin, and M. C. Nuss, “T-ray tomography,” Opt. Lett. 22(12), 904–906 (1997). [CrossRef] [PubMed]
4. J.-H. Son, “Terahertz electromagnetic interactions with biological matter and their applications,” J. Appl. Phys. 105(10), 102033 (2009). [CrossRef]
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. J. Y. Park, H. J. Choi, G.-E. Nam, K.-S. Cho, and J.-H. Son, “In vivo dual-modality terahertz/magnetic resonance imaging using superparamagnetic iron oxide nanoparticles as a dual contrast agent,” IEEE Trans. Terahertz Sci. Technol. 2(1), 93–98 (2012). [CrossRef]
8. S.-H. Cho, S.-H. Lee, C. Nam-Gung, S. J. Oh, J.-H. Son, H. Park, and C.-B. Ahn, “Fast terahertz reflection tomography using block-based compressed sensing,” Opt. Express 19(17), 16401–16409 (2011). [CrossRef] [PubMed]
9. K. Kim, D.-G. Lee, W.-G. Ham, J. Ku, S.-H. Lee, C.-B. Ahn, J.-H. Son, and H. Park, “Adaptive compressed sensing for the fast terahertz reflection tomography,” IEEE Trans. Terahertz Sci. Technol. 3(4), 395–401 (2013). [CrossRef] [PubMed]
10. B. Widrow and S. D. Stearns, Adaptive Signal Processing (Prentice-Hall, 1985).
