The terahertz (THz, 0.1–10 THz) wave, which has garnered increasing attention, lies between the microwave and infrared bands, and possesses several unique physical properties, including a nonionizing nature, low photon energy, good penetrability, and a distinctive fingerprint spectrum. These attributes make THz waves highly promising for applications in detection and have proven viable in fields such as medical diagnostics^[1], noncontact and nondestructive detection^[2], and water content measurement^[3]. Extensive research has focused on developing THz detection systems. For instance, Song et al. proposed a two-dimensional (2D) scanning THz system at 325 GHz using two-tone square-law detection^[4], achieving vector imaging and an amplitude standard deviation of 0.23 mV. Similarly, Cao et al. demonstrated a comparable THz imaging system at 300 GHz^[5] that reached an amplitude standard deviation of 0.13 mV. In contrast, Dülme et al. reported a THz imaging system at 300 GHz that utilized a second reference photodiode to eliminate additional phase-noise terms^[6], achieving an amplitude standard deviation of approximately 0.02 mV and an ultralow phase-noise of 0.034°. Our current work^[7], which was based on an optical frequency comb (OFC) and two-stage mixing, accomplished an amplitude standard deviation of 0.032 mV. Nevertheless, most of these systems rely on 2D mechanical scanning schemes, which limit their ability to accommodate various sample shapes and sizes, particularly sections of elongated samples. Alternatively, the THz computed tomography (CT) scheme is highly suitable for such samples, as it enables visualization of internal structures, showcasing the three-dimensional (3D) configurations of samples^[8].

However, the image quality and resolution of THz-CT cannot match those of X-ray CT due to its larger beam waist radius^[9]. Consequently, additional devices are typically required to improve detection performance^[10–12]. Moreover, sources in most previous THz-CT systems are predominantly based on continuous-wave (CW) gas lasers, frequency multipliers of Gunn tubes, or pulsed THz time-domain spectroscopy (THz-TDS)^[13–17]. These methods suffer from restricted electronic bandwidth and limited frequency tunability in CW systems, and low signal-to-noise ratio and complex signal processing in pulsed systems. In this context, photonics-enabled CW THz sources are recommended to address these challenges, demonstrating advantageous performance in imaging systems^[18].

In this paper, we propose a high-performance CT-based nondestructive detection method using the photonics-enabled THz CW scheme that incorporates custom digital algorithms. In the proof-of-concept experiment, an OFC and a unitraveling carrier photodiode are employed for THz signal generation, with a Schottky barrier diode serving as the receiver for CT detection. Detection information is extracted by a lock-in amplifier, and the original CT sinograms are processed using our specialized THz imaging algorithms. The experimental nondestructive CT demonstration exhibits stable performance over 24 h, achieving a 0.5 mm error margin and enhanced image quality without additional devices.

The experimental setup of the proposed scheme is illustrated in Fig. 1. A laser (NKT) emits a 16 dBm optical signal at 193.414 THz. A phase modulator (PM, EOSPACE) driven by a 33 GHz amplified radio-frequency (RF) signal with a power of 25 dBm generates an OFC. Subsequently, a wave shaper (WS, Finisar) filters out two optical frequency components at 193.249 and 193.579 THz, serving as the optical local oscillator (LO) and carrier, respectively. The spectrum of the OFC and the WS selection is depicted in the inset of Fig. 1. By adjusting the line spacing of the OFC and selecting different orders of comb lines, different frequencies can be obtained. Please note that the photonics-enabled THz source offers an extensive tuning range, while the eventual frequency should take into account other devices’ operation bandwidths in the system. A 3 MHz RF signal with a power of 12 dBm is divided into two branches. One branch modulates the optical carrier via an intensity modulator (IM, EOSPACE) operating at its linear point to generate $\pm 1$st harmonics, while the other serves as a reference signal for subsequent phase locking. The modulated optical signal and the optical LO are power-balanced by erbium-doped fiber amplifiers (EDFAs), coupled, and launched into a unitraveling carrier photodiode (UTC-PD, PMAN-13001, NTT) to generate THz waves at 329.997, 330, and 330.003 GHz, which radiate into the air via a horn antenna. Since the UTC-PD is a polarization-sensitive device, different polarization states of the input optical signal can lead to varying optoelectronic conversion efficiencies. To maximize the utilization of the input optical signal’s power and increase the conversion efficiency, it is necessary to use a polarization controller (PC) to adjust the polarization state of the optical signal before it enters the UTC-PD^[19].

In the free-space stage, two parabolic mirrors (GCC-501101, DHC) focus the THz waves onto a sample, and two additional parabolic mirrors collimate the THz waves toward a receiving antenna.

Regarding reception, the THz waves are converted into electrical signals at 3 and 6 MHz in a Schottky barrier diode (SBD, VDI) through self-mixing. The sample information is embedded within the output amplitude of the SBD. These signals are then amplified by intermediate frequency amplifiers (AMPs, SHF804B) with a gain of 44 dB. A lock-in amplifier (LIA, Zurich Instruments) extracts the sample information from the 3 MHz signal.

In the experiment, continuous line scans are conducted by moving the sample through the THz focus at a velocity of 28.13 mm/s with steps of 5 mm. After each line scan, the sample is rotated by 1 deg until a full 360-deg scan is completed. The total acquisition time is approximately 30 min, primarily constrained by the motor speeds and the number of transceivers.

We assume that the output of the laser is $$E_l(t)=A_l\exp \{j\{[{\omega }_l+{\omega }_{\text{nl}}(t)]t+{\phi }_{\text{nl}}(t)\}\},$$(1)where $A_l$ and ${\omega }_l$ are the amplitude and angular frequency of the laser, and ${\omega }_{\text{nl}}(t)$ and ${\phi }_{\text{nl}}(t)$ denote the angular frequency fluctuation and phase noise of the laser, respectively. After the modulation, we choose two components of the OFC as the optical carrier and LO, which can be given as $$E_c(t)=A_c\exp \{j\{[{\omega }_c+{\omega }_{\text{nl}}(t)]t+{\phi }_{\text{nl}}(t)\}\},$$(2)$$E_{\mathrm{LO}}(t)=A_{\mathrm{LO}}\exp \{j\{[{\omega }_{\mathrm{LO}}+{\omega }_{\text{nl}}(t)]t+{\phi }_{\text{nl}}(t)\}\},$$(3)where $A_c$, $A_{\mathrm{LO}}$, ${\omega }_c$, and ${\omega }_{\mathrm{LO}}$ are the amplitudes and angular frequencies of the optical carrier and LO, respectively. Then, the optical carrier is modulated by the RF signal through the IM to generate harmonics, which can be expressed as $$\begin{gathered} E_{IM}(t)=A_c\{J_0\left(\frac{\pi V_m}{2V_{\pi }}\right)·\exp \{j[{\omega }_c+{\omega }_{\text{nl}}(t)]t+{\phi }_{\text{nl}}(t)\} \\ +J_1\left(\frac{\pi V_m}{2V_{\pi }}\right)·\exp \{j\{[{\omega }_c+{\omega }_{\text{nl}}(t)+\mathrm{\Delta }\omega ]t+{\phi }_{\text{nl}}(t)\}\} \\ +J_{-1}\left(\frac{\pi V_m}{2V_{\pi }}\right)·\exp \{j\{[{\omega }_c+{\omega }_{\text{nl}}(t)-\mathrm{\Delta }\omega ]t+{\phi }_{\text{nl}}(t)\}\}\}, \end{gathered}$$(4)where $J_x(\pi V_m/2V_{\pi })$ is the Bessel function of the order $x$, $V_{\pi }$ is the half-wave voltage of the IM, and $V_m$ and $\mathrm{\Delta }\omega$ are the amplitude and angular frequency of the RF source, respectively. $\mathrm{\Delta }\omega$ can be considered as constant, since the fluctuation of the commercial RF source is negligible compared to the optical source. Higher-order frequency components are also neglected here because of very small amplitudes. Then the signals are combined with the optical LO and injected into the UTC-PD. In this case, the angular frequency fluctuation ${\omega }_{\text{nl}}(t)$ and the phase noise ${\phi }_{\text{nl}}(t)$ of the laser in signals can be eliminated through heterodyning, and the THz signals passing through the sample can be described as $$\begin{gathered} E_{\mathrm{THz}}(t)=A_{T1}\exp \{j[({\omega }_c-{\omega }_{\mathrm{lo}})t+{\mathrm{\Phi }}_1]\} \\ +A_{T2}\exp \{j[({\omega }_c+\mathrm{\Delta }\omega -{\omega }_{\mathrm{lo}})t+{\mathrm{\Phi }}_2]\} \\ +A_{T3}\exp \{j[({\omega }_c-\mathrm{\Delta }\omega -{\omega }_{\mathrm{lo}})t+{\mathrm{\Phi }}_3]\}, \end{gathered}$$(5)where $A_{Tx}$ and ${\mathrm{\Phi }}_x$ represent frequency-dependent attenuated amplitudes and phase shifts of the three THz components, respectively. Since $\mathrm{\Delta }\omega$ is quite small and the $\pm 1$st harmonics are symmetric, we can get ${\mathrm{\Phi }}_1\approx {\mathrm{\Phi }}_2\approx {\mathrm{\Phi }}_3$ and $A_{T2}\approx A_{T3}$. Subsequently, the signal is detected and self-mixed by the SBD, and the output can be expressed as $$\begin{gathered} E_{\mathrm{SBD}}(t)=A_{\mathrm{DC}}+A_{T1}{A}_{T2}\exp [j(\mathrm{\Delta }\omega ·t+{\mathrm{\Phi }}_2-{\mathrm{\Phi }}_1)] \\ +A_{T1}{A}_{T3}\exp [j(\mathrm{\Delta }\omega ·t+{\mathrm{\Phi }}_3-{\mathrm{\Phi }}_1)] \\ +A_{T2}{A}_{T3}\exp [j(2\mathrm{\Delta }\omega ·t+{\mathrm{\Phi }}_3-{\mathrm{\Phi }}_2)] \\ \approx A_{\mathrm{DC}}+2A_{T1}{A}_{T2}\exp [j(\mathrm{\Delta }\omega ·t)] \\ +A_{T2}{A}_{T3}\exp [j(2\mathrm{\Delta }\omega ·t)], \end{gathered}$$(6)where $A_{\mathrm{DC}}$ denotes the direct current (DC) component. As shown in Eq. (6), the information of the sample is within the amplitude of the SBD output, which can be extracted using the LIA.

To eliminate internal blur and noise while preserving the edges of the original sinograms, a THz-CT system requires additional image processing to address the issues arising from the larger focus.

1. 1)First, preprocessing is crucial for the removal of salt and pepper noise; otherwise, the noise will have a devastating impact on the results of CT detection. We employ a median filter to carry out the initial filtering^[20].
2. 2)Second, the sinograms obtained by the THz-CT have blurred edges, which need to be sharpened to enhance the edge sharpness. We utilize the Sobel operator for edge sharpening, as it offers precise edge localization^[21]. The sharpening can be described as $$\begin{gathered} G_x=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]\ast A,G_y=\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right]\ast A, \\ G=\sqrt{G_x^2+G_y^2}\approx |G_x|+|G_y|, \\ {A}^{\prime }=A+\alpha G, \end{gathered}$$(7)where “$\ast$” denotes convolution, $G_x$ and $G_y$, respectively, represent the edges of the original sinogram $A$ extracted from the horizontal and vertical directions, and $\alpha$ ($\alpha \ge 0$) is the sharpening coefficient.
3. 3)After sharpening, a bilateral filter is used to filter noise while preserving the edges^[22]. We assume the weights of pixel value ($r$) and spatial distance ($s$) are $$G_r(p)=\exp (-\frac{{\Vert I_p-I_q\Vert }^2}{2{\sigma }_r^2}),$$(8)$$G_s(p)=\exp (-\frac{{\Vert p-q\Vert }^2}{2{\sigma }_s^2}),$$(9)where $I_x$ denotes the value of pixel $x$, $q$ is the center of the filter window, and $\sigma$ represents the standard deviation. The bilateral filter can be expressed as $$\begin{gathered} \mathrm{BF}=\frac{1}{W_q}\sum _{p\in S}{G}_s(p)G_r(p)\ast I_p \\ =\frac{1}{W_q}\sum _{p\in S}\exp (-\frac{{\Vert p-q\Vert }^2}{2{\sigma }_s^2})\exp (-\frac{{\Vert I_p-I_q\Vert }^2}{2{\sigma }_r^2})\ast I_p, \end{gathered}$$(10)where $W_q$ is the weight sum of each pixel value in the filter window for normalizing the weight. $W_q$ can be expressed as $$\begin{gathered} W_q=\sum _{p\in S}{G}_s(p)G_r(p) \\ =\sum _{p\in S}\exp (-\frac{{\Vert p-q\Vert }^2}{2{\sigma }_s^2})\exp (-\frac{{\Vert I_p-I_q\Vert }^2}{2{\sigma }_r^2}). \end{gathered}$$(11)In flat regions, the pixel value weight $G_r$ in the filter window is similar, making the spatial distance weight $G_s$ the dominant factor influencing the filtering effect. In edge regions, the pixel value weight $G_r$ on the same side of the edge is similar and substantially larger (or smaller) than the $G_r$ on the opposite side of the edge. In this scenario, the $G_r$ on the opposite side has minimal impact on the filtering result, and the edge information is protected effectively, demonstrating good adaptability.
4. 4)Finally, the processed sinograms can be reconstructed into images using the inverse Radon transform-based CT backprojection methods, such as direct backprojection (DBP) and filtered backprojection (FBP)^[23].

According to the configuration in Fig. 1, we experimentally established and conducted a photonic THz-CT system. To characterize the stability and estimate the noise of our system, we calculated the mean and standard deviation of the amplitude captured by the LIA over 24 h, as shown in Fig. 2. Please note that the system gradually stabilizes as thermal and environmental fluctuations decrease over time, leading to data stability improvement after 24 h. Previously, the smallest standard deviation of amplitude reported by us based on OFC was 0.032 mV^[13]. In this work, we enhanced the amplitude stability and achieved a standard deviation of 0.016 mV by removing the effects of secondary mixers. Moreover, the stability performance remained excellent after 24 h. Our system represented the lowest reported amplitude standard deviation above the 300 GHz band, and more detailed analysis on the THz OFC transmitter and its development have also been discussed in our previous work^[7].

To validate our proposed system and analyze the error, we first used two hex wrenches concealed within a shelter for detection. As depicted in Figs. 3(a) and 3(b), the spacings between the two hex wrenches were 15 mm in the vertical direction and 31 mm in the horizontal direction, respectively. We detected the targets for extensive experiments, and one representative set of the images using DBP and FBP is displayed in Figs. 3(c) and 3(d). Both hex wrenches were clearly visible in the images, with the FBP-based image exhibiting higher definition but larger background noise. The measured distances between the two targets were 15.5 mm in the vertical direction and 31.5 mm in the horizontal direction, slightly exceeding the actual measurements. This discrepancy was mainly attributed to the effects of the shelter, whose refractive index was higher than that of air, causing the THz beam to refract. Additionally, the vibration of the motors in motion may also contribute to the error.

To further evaluate the performance of our imaging algorithms, we used a sample of the letter “J” made of polymer material concealed within a section of pipe, fixed on a platform, as depicted in Fig. 4(a). The sample is made by Sculpteo’s VeroWhite Resin material, with a permittivity of 2.72, and the pipes are standard plastic water pipes made of synthetic resin. We obtained images based on DBP and FBP directly from the original CT sinogram [Fig. 4(b)], as shown in Figs. 4(c) and 4(d). The letter appeared blurred, especially in the DBP scheme, and the FBP-based image had many messy lines, which seriously deteriorated the image quality.

Next, we applied our imaging algorithms to the original CT sinogram. After preprocessing, salt and pepper noise was effectively diminished, leading to a noticeable reduction in noise levels in the subsequent images. Applying sharpening enhanced the clarity of edges, as displayed in Fig. 5. While there was an improvement in the imaging results, the extent of enhancement was constrained as the sharpening coefficient increased. Although the edges of the letter became clearer, the main body remained blurry. Subsequently, applying a bilateral filter filtered out noise except for the edges, as shown in Fig. 6. Different sizes of the filter window yielded different filtering effects. A small filter window made it challenging to distinguish between edge and flat regions, resulting in overly preserved edges; conversely, a large filter window tended to overlook subtle edge variations, leading to distortions in the imaging results.

To demonstrate the efficacy of our imaging algorithms, we detected and imaged samples featuring the letters “Z,” “J,” and “U,” with the results displayed in Fig. 7. The letters Z, J, and U were readily recognizable in the images based on DBP and FBP after processing with our specialized THz imaging algorithms. Compared to the images without processing, our imaging algorithms exhibited significant adaptability to both CT imaging schemes. However, there were still distortions at the head and tail of the letters, primarily caused by the scattering of the THz beam due to the corner angles of the samples. Addressing these distortions was crucial for further improving image quality. One potential solution to mitigate this issue was the application of axicon lenses, which can help tackle the scattering effects more effectively^[12].

In conclusion, we propose a photonic CW THz-CT system operating above 300 GHz. Utilizing an OFC and specialized THz imaging algorithms, our system achieves flexible frequency tunability across a wide THz frequency range, which effectively reduces blur and distortion in detection. The experimental system demonstrates stable performance over 24 h, high resolution with a 0.5 mm error margin, as well as nondestructive CT imaging capability. We believe our system holds significant potential for various THz-CT applications, such as medical examination, food evaluation, and defect detection.