Weng-Tai Su, Yi-Chun Hung, Po-Jen Yu, Chia-Wen Lin, Shang-Hua Yang
©SHUTTERSTOCK.COM/PAPAPIG
Visualizing information inside objects is an everlasting need to bridge the world from physics, chemistry, and biology to computation. Among all tomographic techniques, terahertz (THz) computational imaging has demonstrated its unique sensing features to digitalize multidimensional object information in a nondestructive, nonionizing, and noninvasive way. Applying modern signal processing and physics-guided modalities, THz computational imaging systems have now been launched in various application fields in industrial inspection, security screening, chemical inspection, and nondestructive evaluation. In this article, we review recent advances in THz computational imaging modalities in the aspects of system configuration; wave propagation and interaction models; and physics-guided algorithms for digitalizing interior information of imaged objects. Several image restoration and reconstruction issues based on multidimensional THz signals are further discussed, which provides a crosslink among material digitalization, functional property extraction, and multidimensional imager utilization from a signal processing perspective.
In recent years, due to the explosive development of information technology, imaging technology has become more extensive, such as a wide range of applications in consumer electronics, scientific imaging, human–computer interaction, medical imaging, microscopy, and remote sensing. Emphasis is on applied image processing and solving inverse problems using classic algorithms, formal optimization, and modern artificial intelligence techniques. To discover information from invisible to visible, imaging modalities based on different EM spectrum regimes has become an emerging field for establishing digital twins in a virtual world. X-ray imaging is one of the most representative methods to precisely visualize the interior of objects. Despite its capability, high-energy X-ray photons would cause both destructive and ionizing impacts to varieties of object types, preventing further investigations with other object characterization modalities. Functional imaging techniques, such as Raman spectroscopy [1], hyperspectral imaging [2], and temporal imaging [3], were further developed for extracting the physical nature of an object. Recently, THz imaging technology has attracted significant attention due to its noninvasive, nondestructive, nonionizing, material-sensitive, and ultrafast dynamics natures for object visualization and exploration. Through THz–matter interaction, the multifunctional tomographic information of objects can be retrieved even by a noncontact sensing approach.
The THz spectrum, in between microwave and infrared, has often been described as the last frontier of EM wave. As the THz wave can partially penetrate through varieties of optically opaque objects, it carries hidden material information along the traveling path, making this approach a desirable way to visualize the interior of objects without damaging the exterior [4]. Moreover, the rotational, vibrational, torsional frequencies of a great variety of molecules fall in the THz regime. The” self-born spectral barcodes” enrich the retrieved object information with geometrical information and the ingredient information on site. This contact-free, label-free technique has been proven for studying the transport mechanisms of emerging materials with subpicosecond time resolution [5]. Thanks to the great advances of recent high-speed electronics, optoelectronics, and thermal electronics technologies, the heart of the THz imaging systems—THz sources/detectors—becomes much more compact, energy efficient, easy to use, and ready for integrating with commercially available on-chip systems.
In the sub-THz regime, the high-electron-mobility transistor (HMET) and heterojunction field-effect transistor (FET) have been widely used as the critical component of monolithic microwave integrated circuits (MMICs). Those mm-size, high-speed, efficient THz electronics can offer bright THz radiation reaching above 1 THz. The THz photoconductive antenna (PCA) now dominates the market because it can generate/detect coherent EM signals with an ultrabroad covering range from tens of GHz to tens of THz. Notably, these compact THz optoelectronics can be highly integrated with the matured fiber communication, laser, and semiconductor manufacturing industry, making them highly favorable for modular imaging solutions. Considering real-world application schemes, the unmet data acquisition speed drives the need for miniaturizing a massive amount of THz detectors and device integration in a planar platform.
Up-to-date, commercially available THz cameras are mainly based on microbolometer array and FET array technologies. These energy-sensitive THz microbolometer cameras consist of 10,00 to 1 million sensing elements with a sensor size of tens of micrometers operating at higher than 10 Hz, which is now well suitable for real-time THz computational imaging applications. Phase-sensitive or spectral-resolved THz imagers are still gradually reaching a level of maturity for scientific instrument use. Although THz imagers provide a practical way to reveal object interior geometry; specify material information; and reveal ultrafast dynamics further, only limited THz imaging options (e.g., time resolved and frequency resolved) can deliver functional 3D imaging features at the expense of a long data acquisition time for each image voxel. This data acquisition issue currently limits the scope of visualizing 3D objects in size and material variety, so it is worth revisiting the combination of the physical working principles of THz imaging systems; collection of multidimensional THz signals; and computational imaging methods.
Overcoming the signal distortion from different kinds of materials covering the objects is another challenging issue to explore. As the detected signals by those THz cameras and imagers are highly correlated to one or more kinds of object information, such as material identity, ultrafast dynamics, and object geometry, then that raises a question. How can we bridge the THz physical properties and the adopted computational imaging methods to extract invisible objection information?
To answer this question, this work presents a tutorial on state-of-the-art THz computational imaging techniques that combine physics-guided models to achieve high-quality objection information restoration and reconstruction. First, we will focus on the widely used THz computational imaging system configurations that generate the THz dataset in multidimensional domains (e.g., time, frequency). By adequately incorporating wave phenomena into computational imaging methods, we highlight the additional functional imaging and superresolution features that physics-guided models can offer in recent years. From physics-guided models to physics-embedded data-driven methods, we introduce how the data-driven methods can benefit THz image restoration and reconstruction through the assistance of THz multidimensional information, such as spectral characteristics and temporal features. Finally, we demonstrate the effectiveness of physics-embedded data-driven methods for THz computational imaging through experiments on real-world data.
In the past few decades, most THz images have been retrieved by raster scanning an object under test through diffraction-limited spots. Depending on the selection of the THz detector and how THz-radiated waves interact with the illuminated regime, raster-scanned THz imaging systems can reveal plenty of information inside objects, including subsurface structures of culture heritages; contactless conductivity mapping; defect locations in industrial products; and tumor profiles from freshly excised breast tissue. Extended to the raster-scanned configuration, structured-light-based and holographic-based THz imaging systems provide prominent approaches to ease mechanical movement constraints and transfer the data acquisition burden to the computation side. This section will focus on the major THz imaging system configurations to collect object information (Figure 1).
Figure 1. An illustration of THz computational imaging system configurations. (a) THz CS imaging. (b) THz pulse imaging. (c) THz holographic imaging. (d) THz computed tomography (CT).
Instead of raster scanning objects through either moving them or steering the light path, compressive sensing (CS) imaging is a computational method to reconstruct object information with low-rank detection. Depending on the information sparsity of the object and CS modalities, the number of measurements can be greatly reduced compared to the pixel number. This brings three major features to the THz imaging system configuration—high-speed data acquisition; the use of a single THz detector; and no mechanical moving parts—at the expense of increased computational effort and a limited imaging field of view. The basic operation principle of THz CS imaging is to illuminate a series of patterned THz beams on the object. The multiplication of pixelated masks and the vectorized object geometry then record the total transmitted THz signal by a single-pixel (bucket) THz detector as described in the equation ${\bf{s}} = {\bf{Ax}}$, where s, A, and x are the recorded THz signal, the pixelated mask, and the vectorized image to be reconstructed, respectively. Each row of the sensing matrix represents an illumination pattern. To reconstruct the image, we can formulate the following optimization problem: \[\mathop{\min}\limits_{x}{\left\Vert{\bf{s}}{-}{\bf{Ax}}\right\Vert}_{2} + \lambda\phi{(}{x}{)} \tag{1} \] where ${\lambda}$ is a hyperparameter, and the regularizer ${\phi}{(}{x}{)}$ is used to promote the latent physical properties of the solution, such as sparsity and/or low rank. For example, L1-norm is usually utilized to promote the sparsity of the solution. Depending on the formulation of the regularizer in (1), convex and/or nonconvex optimization solvers (e.g., alternating direction method of multipliers; alternative optimization; and gradient descent methods) can be applied to find the global or local optimal solution. A great variety of high-speed, high-modulation-depth THz spatial light modulators have been invented in the past 10 years. With the combination of photoconductive modulation and digital light processing techniques, a real-time THz video (frame rate: 6 fps; pixel number: 32 × 32) has been demonstrated [6].
Pulsed imaging, or time-of-flight imaging, records the traveling THz pulsed signal starting from the THz source (the tested sample) to the THz detector. A typical system configuration is based on a THz time-domain spectroscopy (THz-TDS) system performed in the reflection mode. The time-delayed signals and the echoes reflect inner object structures \[{d} = \frac{{c}{\Delta}{t}}{2n}\] where d is the distance from the sample surface to the sample internal interface; and c, n, and ${\Delta}{t}$ are the speed of light in vacuum, the sample refractive index, and the time delay of the recorded signal. While conventional THz-TDS systems suffered from long data acquisition time due to mechanical components, modern techniques, such as asynchronous optical sampling THz-TDS systems and electronically controlled optical sampling systems, significantly speed up data acquisition time with orders of magnitude—from hours to milliseconds—bringing THz pulsed imaging modality to become a unique but powerful tool to discover multifunctional material properties in wide ranges of spatial, temporal, spectral, and material domains, simultaneously. To this extent, THz pulsed imaging has been allocated for many applications, including drug detection and compressed sensing imaging.
The THz holographic imaging system records both the transmitted THz beam profile and the phase changes through the imaged objects [7]. This recorded THz amplitude/phase mapping represents the thickness, topology, features, and dielectric response of the object, which can further reconstruct the whole object geometry beyond diffraction-limited spatial resolution. A typical THz holographic imaging system is based on a Mach–Zehnder THz interferometer setting—the superposition of the coherent scattered THz beam and the reference THz beam—as described in (2) \[{\bf{H}}{(}{r}{)} = {\bf{U}}_{\bf{r}}{(}{r}{)}^{2} + {\bf{U}}_{\bf{r}}^{*}{\bf{U}}_{\bf{o}}{(}{r}{)} + {\bf{U}}_{\bf{r}}{\bf{U}}_{\bf{o}}^{*}{(}{r}{)} + {\bf{U}}_{\bf{o}}{(}{r}{)}^{2} \tag{2} \] where ${\bf{H}}{(}{r}{),}{\bf{U}}_{\bf{r}}{(}{r}{)}$, and ${\bf{U}}_{\bf{o}}{(}{r}{)}$ are the THz interferogram; the electric field of the reference beam; and the electric field of the coherent scattered THz beam from the object, respectively. Followed by performing Fresnel approximation or the angular spectrum method, the image of the object under test can be reconstructed [7]. Recent trends in time-domain THz holography, one-shot THz holography, and video-rate THz digital holography have brought extensive attention to this field. These modalities extend the system capability to functional 3D tomography at applicable frame rates. With high-fidelity amplitude and phase reconstruction, the depth resolution can reach 2.2 μm (${\lambda} / {284}$ at 0.48 THz), comparable with the conventional optical imaging system.
As the analogy of X-ray computed tomography (CT), a typical THz CT system is composed of a THz source–detector pair and an object-positioning stage that can rotate and move in the local coordinate system. While the THz beam transmits through the objects’ local regime, the loss function of the THz signal, f(x, y) can then be recorded by the THz detector in use. The data content of recorded signals could be field strength, spectral amplitude/phase, or polarization changes. By sampling from different angular and radial positions, one can define the function ${R}{(}{\theta},{\rho}{)}$, the Radon transform, to relate the object information with the recorded signals \[{\bf{R}}{(}{\theta},{\rho}{)} = \mathop{\int}\nolimits_{{-}\infty}\nolimits^{\infty}{\mathop{\int}\nolimits_{{-}\infty}\nolimits^{\infty}{f}}{(}{x},{y}{)}{\delta}{(}{\rho}{-}{x}{\cos}{\theta}{-}{y}{\sin}{\theta}{)}{dxdy}\] where ${\theta}$ and ${\rho}$ are the angular and radial coordinates of the projection line ${(}{\theta},{\rho}{)}$, respectively; and ${\delta}{(}{\cdot}{)}$ is the Dirac impulse. The inverse Radon transform model through filtered backprojection [8] or iterative methods (i.e., simultaneous algebraic reconstruction technique or expectation–maximization) defines the way to reconstruct each object voxel and its physical contents from the projected dataset. In a practical scenario, the static/dynamic measurement misalignments and the error from data discretization significantly impact the reconstructed data, which is a sparse research field waiting to be explored.
Many THz imaging modalities have been developed based on the fundamental properties of THz waves—amplitude, phase, spectrum, polarization, wavefront, chirality, coherence, wave propagation, and light–matter interactions [9]. Physics-guided methods are commonly used in the sequential stages of data preprocessing, image computation, and image postprocessing to translate the object information from the raw data (Figure 2). Although most THz computational imaging studies are based on specific data types, it should be noted that the fusion of multiple data types, such as spatiotemporal–spectral data, can further enhance image quality. With the utilization of prior knowledge (e.g., beam profile, spectral fingerprints, and ultrafast light–matter interaction), recent advances in physics-guided THz imaging have shown extraordinary features, such as noninvasive carrier dynamics mapping and material inspection on a nanoscopic scale [10].
Figure 2. A flowchart of physics-guided THz computational imaging.
Unlike visible light, most THz imaging systems deal with a spatial resolution close to or shorter than the beam wavelength. As the ray optics can no longer explicitly express the diffraction-limited THz beam profile landing on the object, wave optics play essential roles in THz image reconstruction. The mathematical modeling of THz wave propagation based on Gaussian beam theory profiles the spatially varying beam shapes as the point spread function (PSF) along the traveling axis. With the accurate approximate absorption coefficient of the object, beam divergence, and depth of focus, the deconvoluted image’s spatial resolution can reach the subwavelength scale. The phase characteristic of the THz beam is another powerful tool to beat the diffraction limit. The superposition of two coherent THz collimated beams, the planar reference beam and the object beam, forms spatial modulated fringe/speckle patterns that record the optical path differences introduced by the object. Different approaches of scalar diffraction theory (e.g., Fresnel diffraction, Huygens convolution, and Rayleigh–Sommerfeld diffraction), which are applicable to in-line and off-axis holography, have been used for reconstructing object 3D information.
As the image resolution is mainly determined by the detector aperture, hologram size, and off-axis angle, it is possible to achieve a phase resolution of ${\lambda} / {10},{000}$—a nano-scale axial resolution—using a THz beam [7]. Fresnel equations formulate the light–matter interactions on the material interface and inside the object, which is a well-known approach to precisely reconstruct THz images. With the measurable dielectric response prior to the imaged object, it can nicely resolve the thickness variations of multiple stacking layers and the defected regions inside. Along with the combination of the Drude model, which extracts complex conductivity from the material dielectric response, the electrical-property mapping inside the object can be further visualized. By capturing THz light–matter interaction in the near-field regime, research works have demonstrated object information reconstruction based on Fresnel diffraction in the microscopic world.
It is still feasible to record nanoscale THz imaging by introducing a metallic scattering tip that intensifies the local THz electric field in several orders of magnitude near the tip region. Adopting an aperture synthesis technique through a coherent phased-array THz system is another promising trend to achieve beyond diffraction-limited spatial imaging resolution by synthesizing amplitude and phase superpositioned images. Depending on the light–matter interaction properties of the observed object, single-molecule ultrafast motion can be recorded by THz near-field imaging with submolecular resolution. Moreover, data acquisition speed is an everlasting issue for real-time THz imaging. Recent CS imaging modalities quickly adapt to many THz imaging systems, including THz CT, THz pulsed imaging, THz holography, and more system configurations, which maintain the same object visualization features (e.g., hyperspectral, time-resolved, and near-field imaging) with less measurement time [11]. To mitigate the diffraction effect of the reconstructed image caused by the distance between the image plane and object plane, the physics guide of inverse Fresnel diffraction can be utilized with the inverse problem of the THz CS.
Like in other EM bands, the preprocessing and the postprocessing in THz computational imaging can adopt most of the methods developed in signal processing and image processing. The subspace projection on the detected THz spectra achieves the noise removal and, in the meanwhile, preserves the THz spectral features [12]. The statistics-based average over the THz time-resolved signals can be used to suppress the various types of noises [13]. The frequency and wavelet filters are also applied to the THz signal to achieve noise removal [14]. In the image domain, an isotropic diffusion algorithm and a mixture of Gaussian densities are utilized to achieve denoising and image segmentation, respectively [15]. Based on the rapid development of learning-based approaches in visible light imaging, several deep learning models are also demonstrated in THz computational imaging [16]. To date, THz imaging processing techniques are less coupled with physics-related working principles. This could be an exciting playground to extend the further capability of THz computational imaging from the signal processing perspective.
The energy loss encoded in the reflected/transmitted THz signals [Figure 4(a)] can distinguish interior object information based on the Fresnel equation; through THz direct detection, energy-based THz imaging is highly favorable for functional 2D imaging or THz CT. The current works focus on feature extraction methods rather than reconstruction models since different feature extraction approaches can directly reflect corresponding material characteristics. For example, the attenuation spectra of THz spectral amplitude signals directly map to the material conductivity locally. While scanning across the whole field of view, the spatial distribution of conductivity in materials can be easily investigated. In [17], Kawase et al. captured multispectral THz images for identifying illicit drugs. The selected THz frequencies are based on the distinctive spectral fingerprints of each chemical. This way, the spatial distributions and concentration levels of different chemicals can be nicely classified and revealed from THz spectral amplitude signals.
The aforementioned examples demonstrate that one can directly extract various material information through different THz signal features. However, the energy-based methods suffer a great deal from several issues in THz images, such as blurs, noises, low contrast, and low image resolution, resulting in difficulties for effective material analysis. To extend the practical use of THz imaging, research groups have been looking for frequency-resolved and time-resolved THz imaging modalities that exploit additional light–matter interaction information to improve image qualities and functional imaging capabilities.
Time-resolved THz imaging provides meta-information of imaged objects simultaneously, including time of flight; ultrafast phenomena; charged carrier dynamics and distributions; material spectral fingerprints; and object topology. While the THz pulses pass through the object under evaluation, the reflected time-resolved THz electric field signals from incident inner shells arrive at the THz detector at different timestamps [Figures 3 and 4(b)], revealing the geometric information of the object. Through FFT operation, the amplitude/phase spectral responses can be further extracted to reveal the objects’ material information. To date, researchers further extend the capability of time-resolved THz spectroscopy (TRTS) systems to visualize the ultrafast dynamics inside photoconductive objects. The microscopic electrical properties, like conductivity, mobility, and doping mapping, can be further derived through the Drude model.
Figure 3. An example of THz computational imaging based on a time-of-flight system. (a) An illustration of the measured object composed of three paper sheets with letters “Y,” “R,” and “G,” respectively. (b) The THz electric field mapping at different time stamps. (c) The images retrieved by a conventional method and a physics-guided method, respectively. (d) The reconstructed 3D image of the first paper sheet by using the time-of-flight information.
Figure 4. An illustration of THz physics-guided features. (a) The peak-to-peak value ${\Delta}{p} {-} {p}$ of a THz signal is highly correlated to energy loss and can be used for energy-based THz imaging. (b) The time delay ${\Delta}{t}$ of a THz signal can be used to estimate the geometric information of an object. (c) and (d) The (c) phase and (d) amplitude spectra of THz signals associate with the material permittivity and permeability and can be utilized for frequency-resolved THz imaging.
While a TRTS system combines with a scattering-type scanning near-field optical microscopy system, the time-resolved THz hyperspectral nano-imager can noninvasively resolve nanoscale carrier profiling inside advanced materials [10]. To achieve a better data acquisition rate, Zanotto et al. [18] demonstrated a time-domain THz CS imaging system that offers time-resolved, hyperspectral features without the need for mechanical raster scanning. As TRTS is widely used in pure source identification, Lin et al. [12] developed a blind hyperspectral unmixing method—Hyperspectral Penetrating-type Ellipsoidal Reconstruction (HYPERION)—for THz blind source separation. This work shows an unprecedented chemical blind separation mapping capability of time-resolved THz hyperspectral imaging without collecting extensive data or sophisticated model training.
Compared with THz direct detection, THz coherent detection offers not just extremely high dynamic range and low noise equivalent power but the ability to resolve the amplitude/phase spectrum of the THz signals [Figure 4(c) and (d)]. The amplitude and phase dynamics correspond to the optical delay path pass through the imaged object and its material properties (e.g., complex permittivity), which makes THz coherent imaging highly favorable for high-precision 3D imaging. Under a swept-source optical coherence tomography (OCT) configuration, the material spectral response and the 3D object geometry can be easily reconstructed based on the optical frequency-domain reflectometry (OFDR) model [19]. Typically, THz photomixers and THz quantum cascade lasers are the most promising candidates as the swept sources of THz OCT systems due to their excellent phase stability, ultrabroad frequency tuning range, fast frequency-sweeping rate, and compact size. To achieve better-reconstructed imaging quality, the modeling of wave propagation, wavefront profile, and dielectric response at each frequency shall be incorporated in the OFDR model. It should be noted that by synthesizing the broad-frequency-sweeping THz signals, a THz OCT system can form high-quality THz pulses used for time-of-flight THz imaging.
Many imaging methods have been developed based on the light–matter interaction in the THz frequency range in the past decades. Based on THz absorption imaging modalities, the material complex refractive index mapping can be profiled through the Fresnel equation, extracted by the THz power loss while propagating through the tested object boundary. However, the application scopes of those physics-based methods are severely limited since they normally require a sufficient amount of prior knowledge about a measured object to simplify the complex physical models. To break this limitation, data-driven approaches, especially deep neural networks (DNNs), have been attracting intensive attention due to their excellent learning capability.
Data-driven methods are mainly based on deep learning models, which do not resort to any explicit transform model but are learned from representative big data. Deep learning has revolutionized the aforementioned physics-based paradigm in image restoration, for which the physical priors and knowledge can be embedded into data-driven methods. For instance, an appropriate loss function design can help supervise the data-driven model to learn the latent physical priors; the physics-embedded model architecture (i.e., hierarchical learning module) can leverage physical cues in data. Integrating the physical priors into data-driven methods not only can improve their model performance but also enhances noise immunity. However, since a large-scale THz dataset is currently not available, there are only a few works demonstrating the advantages of the physics-embedded data-driven methods. For example, in [20] a wave propagation model is embedded into a deep learning model to reconstruct THz images; an image degradation model with THz point-spread functions is utilized in [21] for data augmentation for THz superresolution image reconstruction.
To explore physics-embedded data-driven THz computational imaging in a broader scope, we specifically choose THz imaging works based on data-driven models with time-resolved THz imaging systems. This is because the retrieved time-resolved signal contains rich curated pixel-wise spectral information, such as the prominent amplitude/phase spectral information corresponding to specific light–matter interaction characteristics of THz waves passing through materials.
Moreover, data-driven models can excellently fuse the different information of THz signals, such as amplitude/phase spectra and the time-resolved THz signals, to achieve superior image restoration [22]. In this section, as an example, we demonstrate the capability of physics-embedded data-driven computational imaging on boosting the performance of a THz-TDS CT system.
To measure the THz tomographic images of an object, as illustrated in Figure 5(b), we established an Asynchronous Optical Sampling (ASOPS) THz-TDS CT system prototype (see also [22]). The system is composed of two asynchronous femtosecond lasers with a central wavelength of 1,560 nm and a power level of 40 mW; a pair of an InGaAs/InAlAs THz PCA source and a low-temperature grown-InGaAs/InAlAs THz PCA detector; a linear and rotation motorized stage; a four plane-convex THz lens with 50-mm focal length; a transimpedance amplifier; and a unit of data acquisition and processing. The repetition rates of the two asynchronous femtosecond lasers are 100 MHz and 100 MHz + 200 Hz, respectively. With the previous configuration, the ASOPS THz-TDS system delivers a 0.1-ps temporal resolution and the THz frequency bandwidth of 5 THz. Additionally, it provides THz pulse signals of 60 µW with a 41.7-dB dynamic range from 0.3 to 3 THz and 516 fs at full width at half-maximum (FWHM) where the FWHM of the THz PSF is 1.25 mm.
Figure 5. The THz dataset preparation process. (a) The ground-truth blueprint utilized for 3D-printed fabrication and the fabricated object hidden in an opaque cup for raster scanning. (b) The raster scanning based on the THz-TDS system. (c) The recorded THz time-resolved signals of two voxels [centered at the brown and blue rectangles in (a)] in the air and Deer body, respectively. (d) The amplitude spectrum of the THz time-resolved signals in (c), where the red points indicate the 12 selected prominent frequencies for generating THz multispectral images. (e) The Time-max image [the maximum values of THz time-resolved signals in (c)]. (f) The THz multispectral images at the 12 selected frequencies in (d).
In the THz-TDS CT system, each object [e.g., the covered 3D-printed deer in Figure 5(a)] is placed on the motorized stage in the THz path between the light source and detector of the THz-TDS system. The seven hand-hold-sized sample objects are fabricated by a fused deposition modeling3D printer with the material of high-impact polystyrene [23] due to its high penetration characteristics in THz bands. Additionally, the sample objects are masked by the paper cup with diameter of 7 cm during the scanning. The physical size of the 3D-printed object (Deer) is approximately 50 mm × 50 mm × 56 mm (width × depth × height), as shown in Figure 5(a). Raster scans are performed on the measured object in multiple view angles, covering a rotational range of 180° (step size: 6° between two neighboring views); a horizontal range of 72 mm (step size: 0.25 mm); and a variable vertical range corresponding to the object height (step size: 0.25 mm). Under this scanning step of 0.25 mm and the 1.25-mm THz PSF, our system delivers spatial resolution in a sub-mm scale. We obtain 30 projections of each object, which are then augmented to 60 projections by horizontal flipping. The ground truths of individual projections are obtained by converting the original 3D printing files into image projections in every view angle, as shown in Figure 5(b).
During the scanning, the THz-TDS CT system profiles each voxel’s THz time-resolved signal with a 0.1-ps temporal resolution, whose amplitude corresponds to the strength of the THz electric field. Based on the dependency between the amplitude of a time-resolved signal and the THz electric field, in conventional THz imaging, the maximum peak of the signal (Time-max [23]) is extracted as the feature for a voxel. The retrieved Time-max images, as illustrated in Figure 5(e), can deliver great signal-to-noise ratio (SNR) and clear object contours. However, the conventional THz imaging based on Time-max features suffers from several drawbacks, such as the undesired contour along the object boundaries; the hollow and holes in the body region; and the blurs in the high spatial-frequency regions. To address this limitation, the spectral information of THz temporal signals can be utilized to supplement the conventional method based on Time-max features since the voxel of the material behaviors are encoded in both the phase and amplitude of different spectral bands, according to the Fresnel equation.
Figure 5(c) illustrates time-domain THz signals measured in the air and in the body and leg of a 3D-printed deer, respectively. While the THz beam is passing through the object, the attenuated THz time-domain signal is encoded with the thickness and material information of the THz-illuminated region. By processing the Time-max images captured at different view angles, the 3D profile of the printed deer can be further reconstructed. Although this conventional way is well fitted for visualizing 3D objects, the inherent diffraction behavior and strong water absorption nature of THz waves induce various kinds of noise sources as well as the loss of object information. This leads to the undesired blurring, distorted, and speckled phenomenon of functional THz images. Existing works have tackled this issue to restore clear images via estimating PSFs [24], image enhancement [25], machine learning [26], and so forth. Their performance is, however, still severely constrained by diffraction-limited THz beams.
To break the limitations mentioned earlier, we show that curated physical guides can help extract the prominent spatiospectral information of hidden objects to facilitate the learning of data-driven models so as to boost the performance of reconstructing deep-subwavelength tomographic images of hidden objects. For example, since THz waves traveling through the object can have distinct light–matter interaction at different frequencies, the supplemental THz multispectral images can provide additional object information (e.g., contours and depth).
As shown in Figure 5, each 2D THz image is composed of an array of THz time-resolved signals with a temporal resolution of 100 fs. By taking Fourier transform on the THz time-resolved signals individually, we can extract their pixel-wise multispectral feature maps [i.e., Figure 5(d) and (f)]. Nevertheless, due to the large number of spectral bands with the measured THz data, it is required to sample a subset of the spectral bands to reduce the number of model parameters. Besides, the THz-TDS system offers relatively better SNR levels in a frequency range of 0.3–1.3 THz. Considering the water absorption lines and transmission windows [27] in 0.3–1.3 THz, we select 12 frequencies at 0.32, 0.57, 0.60, 0.77, 0.80, 1.04, 1.06, 1.09, 1.14, 1.19, 1.21, and 1.24 THz.
Prominent multispectral information, including both amplitude and phase at the selected frequencies, is curated and then employed to restore clear 2D images. Specifically, Figure 6 illustrates multiple 2D THz images of the same object retrieved at the selected frequencies, showing different contrasts and spatial resolutions as these hyperspectral THz image sets have different physical characteristics through the interactions of different THz bands with an object. Specifically, the lower-frequency phase images offer relatively accurate depth information due to their higher SNR level, whereas the higher-frequency phase images offer finer contours and edges because of the shrinking diffraction-limited wavelength sizes (from left to right in Figure 6). The phase also contains, however, a great variety of information of light–matter interaction that could cause learning difficulty for the image restoration task.
Figure 6. An illustration of THz multispectral amplitude and phase images measured from the deer object.
To address this issue, we utilize selected amplitude bands as complementary information. Although the attenuated amplitude bands cannot reflect comparable depth accuracy levels as phase bands, the amplitude bands still present superior SNR and more faithful contours, such as the location information of a measured object.
In summary, fusing the amplitude and phase images from low frequency to high frequency brings the following advantages. Since the low-frequency THz signal provides precise depth (the thickness of an object) and fine edge/contour information in the phase and amplitude, respectively, they together better delineate and restore the object. In contrast, the high-frequency feature maps of amplitude and phase, respectively, provide better edges/contours and precise position information, thereby constituting a better object mask from the complementary features. With these spectral properties of THz images as physical guides, we can extract rich information from a wide spectral range in the frequency domain to simultaneously restore the 2D THz images without any additional computational cost or equipment, which is beneficial for further development of practical THz imaging tasks.
As different EM bands interact with objects differently, THz waves can partially penetrate through various optically opaque materials of objects and carry hidden object tomographic information along the traveling path. This unique feature can guide a new approach to visualize the essence of 3D objects, which other imaging modalities cannot achieve. However, the Time-max images in Figure 5(e) preserve only energy distribution information, which is not enough for learning an effective data-driven model for detailed object information restoration, leading to difficulties in effective material analysis. To address the problem, as mentioned earlier, thanks to the rich object tomographic information carried in the amplitude/phase spectrum of frequency-resolved THz signals corresponding to their optical delays passing through an imaged object with different material properties (e.g., water absorption and object thickness), we can use it to restore the 2D THz images of individual views to facilitate the further development of THz imaging. The key idea is to fuse spatiospectral features with different characteristics on a common ground via data fusion approaches to guide the feature fusion for fine restoration.
Moreover, we reveal in [22] that directly learning from the full spectral information to restore THz images usually leads to an unsatisfactory performance. The main reason is that the full spectral bands of THz signals involve diverse characteristics of materials, noises, and scattered signals, causing difficulties in the model training of DNNs. Instead of learning from the full spectral information of THz signals, this problem can be tackled by extracting prominent multispectral information from the complementary amplitude and phase of THz signals in the curated selection 12 of frequencies bands [22].
To accommodate the curated multispectral amplitude and phase bands for learning the THz image restoration model in a multiscale manner, as depicted in Figure 7, we devise a Subspace-Attention guided Network (SARNet) based on the multiscale U-Net in [28]. The key idea of SARNet is to fuse the multispectral features with different characteristics, selected at the 12 prominent frequencies (i.e., the water absorption lines) of THz signals, on a common ground via deriving the shared latent subspace and discovering the short/long-range dependencies between the amplitude and phase to guide the feature fusion. SARNet can capture such complementary spectral characteristics of materials to restore damaged 2D THz images effectively. In summary, SARNet merges the THz spatiotemporal–spectral data, physics-guided data-driven models, and material properties for high-precision THz tomographic imaging.
Figure 7. (a) The overall network architecture of SARNet consisting of five scale branches, where the finest-scale scale takes the Time-max image as input, and each of the second to fifth takes six images of spectral frequencies (three bands in amplitude and three bands in phase) as inputs. The three gray blocks show the detailed structures of (b) spectral fusion, (c) channel fusion, and (d) Conv-Block. SC: skip connection; ReLU: rectified linear unit.
Specifically, SARNet is composed of an encoder (spectral-fusion module) with five branches of different scales (from the finest to the coarsest) and a decoder (channel-fusion module) with five corresponding scale branches. The spectral fusion block in each scale branch of the encoder involves a subspace-attention-guided fusion module (SAFM), a convolution block (Conv-block), and a downsampler, except for the finest-scale branch that does not employ SAFM. To extract and fuse multispectral features of both amplitude and phase in a multiscale manner, the encoder takes the Time-max image as the input of the finest-scale branch and receives to its second to fifth scale branches 24 images of selected predominant spectral frequencies extracted from the pixel-wise THz signals, where each branch takes six images of different spectral bands (three bands of amplitude and three bands of phase) to extract learnable features from these spectral bands. To reduce the number of model parameters, these 24 amplitude and phase images (from low to high frequencies) are downsampled to four different resolutions and fed into the second to fifth scale branches in a fine-to-coarse manner, as illustrated in Figure 7.
Since the multispectral amplitude and phase feature maps are significantly different in nature, the proposed SAFM fuses these features via first learning a common latent subspace shared between the amplitude and phase features to facilitate associating the short/long-range amplitude-phase dependencies. Projected into the shared latent subspace, the spectral features of amplitude and phase components, along with the downsampled features of the upper layer, can then be properly fused together on a common ground in from the finest to the coarsest scale to derive the final latent code. The detailed design of SAFM can be found in [22].
The Conv-block(L) contains two stacks of L × L convolution, batch normalization, and rectified linear unit (ReLU) operations. Because the properties of the spectral bands of amplitude and phase can be significantly different, we mostly use ${L} = {1}$ to learn the best linear combination of multispectral features to avoid noise confusion and reduce the number of model parameters. The upsampler and downsampler perform ${2}{\times}$ and ${(}{1} / {2}{)}{\times}$ scaling, respectively. The skip connections (SCs) directly pass the feature maps of different spatial scales from the individual encoder branches to the Chanel Fusion blocks of their corresponding branches of the decoder.
In the decoder path, the Channel Fusion block in each scale branch involves an upsampler, a Channel Attention Module (CAM), and a Conv-block. CAM aims to fuse the multiscale features from different spectral bands in the channel dimension using the channel attention mechanism detailed in [22]. The Conv-block has the same functional blocks as that in the encoder. Each decoding-branch receives a “shallower-layer” feature map from the corresponding encoding branch via the SC shortcut and concatenates the feature map with the upsampled version of the decoded “deeper-layer” feature map from its coarser-scale branch. Besides, the concatenated feature map is then processed by CAM to capture the cross-channel attention to complement the local region for image restoration. All the implementation details of SARNet can be found in [22].
Herein, we demonstrate the performance gain brought by the physics guides on THz computational imaging, particularly on THz image restoration and 3D CT. In the experiments, a total of seven sample objects are printed, measured, and aligned for evaluation. (The THz-TDS Image Dataset can be found at https://github.com/wtnthu/THz_Data.)
To the best of our knowledge, there is no method specially designed for restoring THz images besides Time-max. Thus, we compare our method against three representative convolutional neural network (CNN)-based image restoration models, including DnCNN-S [29], NBNet [30], and the baseline U-Net [28], trained on the image pairs of Time-max images and their corresponding ground truths. For objective quality assessment, we adopt two widely used metrics, including the peak SNR (PSNR) and structural similarity index measure, to respectively measure the pixel-level and structure-level similarities between a restored image and its ground truth, as shown in Table 1. We also adopt the mean-squared error (MSE) between the cross-sections of a reconstructed 3D tomography and the corresponding ground truths for assessing the 3D reconstruction accuracy as compared in Table 2.
Table 1. The quantitative comparison (PSNR and SSIM) of THz image restoration performances with different methods on seven test objects.
Table 2. The quantitative comparison of MSE between the cross-sections of a reconstructed 3D tomography and their ground truths with different methods on seven test objects.
Table 1 shows that the four CNN-based data-driven restoration methods all achieve large performance gain over Time-max. Thanks to the physics guides, SARNet outperforms other peer methods in terms of MSE, SSIM, and PSNR as shown in Tables 1 and 2.
For qualitative evaluation, Figure 8 illustrates the 3D reconstructions of Deer, Robot, and Skull (complete 3D CT results on all objects can be found at https://github.com/wtnthu/THz_Tomography_2022), showing that SARNet reconstructs significantly clearer and faithful 3D images with finer details, such as more correct thickness of the body and clear antlers of Deer; accurate facial contour of Skull; and accurate reconstruction of the gun in Robot’s hand, achieving by far the best 3D THz tomography reconstruction quality in the literature. Both the quantitative and qualitative evaluations confirm a significant performance leap with the physics-guided SARNet over the competing methods, demonstrating the high efficacy of curated physical guides on THz CT computational imaging.
Figure 8. An illustration of 3D CT reconstructions on deer, robot, and skull from left to right: (a) Time-max, (b) DnCNN-S [29], (c) NBNet [30], (d) U-Net [28], (e) SARNet, and (f) the ground truth.
THz computational imaging is an emerging research field with great application potential. This article briefly introduced the associations of physical characteristics with major THz computational imaging configurations, including CS imaging, pulse imaging, holographic imaging, and CT. With THz computational imaging, multiple physical characteristics can be modeled and contribute to various applications, such as material exploration, industrial inspection, security screening, chemical inspection, and nondestructive evaluation. With the rich information in multidimensional THz signal behaviors, developing the different THz imaging modalities guided by physical characteristics is the key to the door of THz computational imaging. We have illustrated a THz tomographic imaging example showing how to leverage the prominent spectral information carried in THz time-domain signals, guided by the water absorption profile of THz light–matter interaction, to devise a data-driven CNN model for effectively restoring corrupted THz images. With such well-designed image processing, our experiments confirm a performance leap from the relevant state of the art in THz tomographic imaging. This also sheds light on the strength and opportunities of signal processing in THz computational imaging.
This work was financially supported in part from the Young Scholar Fellowship Program under Grant MOST 110–2636-E-007–017 and in part under Grant 111–2221-E-007–046-MY3 by the Ministry of Science and Technology (MOST) in Taiwan.
Weng-Tai Su (wengtai2008@hotmail.com) received his Ph.D. degree in electrical engineering from National Tsing Hua University, Hsinchu, Taiwan, in 2022. He is currently an intern at Microsoft Taiwan, Taipei 110419, Taiwan, working on face recognition. His research interests include machine learning; image and video processing; and computer vision. He is a Member of IEEE.
Yi-Chun Hung (yichunhung@g.ucla.edu) received his M.S. degree in electrical and computer engineering from the University of California, Los Angeles in 2022. He is currently a research assistant in the Department of Electrical Engineering, National Cheng Kung University, Tainan 701401, Taiwan. His research interests include computational imaging, computer vision, and deep learning.
Po-Jen Yu (jerry321ab@gmail.com) received his B.S. degree in power mechanical engineering from National Tsing Hua University, and he is a research assistant in the research group led by Shang-Hua Yang at National Tsing Hua University, Hsinchu 300044, Taiwan. His research interests include terahertz lens design, manufacturing processes, machine design, and terahertz imaging.
Chia-Wen Lin (cwlin@ee.nthu.edu.tw) received his Ph.D. degree from National Tsing Hua University (NTHU), Taiwan, in 2000. He is currently a professor with the Department of Electrical Engineering and the Institute of Communications Engineering, NTHU, Hsinchu 300044, Taiwan. He served as a Distinguished Lecturer of the IEEE Circuits and Systems Society (2018–2019). He has served on the editorial boards of IEEE Transactions on Multimedia; IEEE Transactions on Image Processing; IEEE Transactions on Circuits and Systems for Video Technology; and IEEE Multimedia. He was the chair of the ICME Steering Committee from 2020 to 2021. He served as the TPC cochair of ICIP 2019 and ICME 2010 and general cochair of VCIP 2018. He received best paper awards from VCIP in 2010 and 2015. He is a Fellow of IEEE.
Shang-Hua Yang (shanghua@ee.nthu.edu.tw) received his Ph.D. degree in electrical engineering from the University of Michigan Ann Arbor in 2016. He is an assistant professor in the Department of Electrical Engineering and the deputy director of the National Tsing Hua University (NTHU) Terahertz Optics and Photonics Center at NTHU, Hsinchu 300044, Taiwan. He has received several prestigious awards, including the IEEE Antennas and Propagation Society Doctoral Research Award (2014); MOST Young Scholar Fellowship (2018); and Human Frontier Science Program Research Grant Award (2020). His research interests include ultrafast electronics, microwave photonics, nanophotonics, computational imaging, artificial intelligence, terahertz optoelectronics, terahertz communication, and integrated terahertz systems. He is a Member of IEEE.
[1] D. A. Long, Raman Spectroscopy. New York, NY, USA: McGraw-Hill, 1977, pp. 1–12.
[2] G. Lu and B. Fei, “Medical hyperspectral imaging: A review,” J. Biomed. Opt., vol. 19, no. 1, p. 010901, 2014, doi: 10.1117/1.JBO.19.1.010901.
[3] B. H. Kolner, “Space-time duality and the theory of temporal imaging,” IEEE J. Quantum Electron., vol. 30, no. 8, pp. 1951–1963, Aug. 1994, doi: 10.1109/3.301659.
[4] C. Jansen et al., “Terahertz imaging: Applications and perspectives,” Appl. Opt., vol. 49, no. 19, pp. E48–E57, 2010, doi: 10.1364/AO.49.000E48.
[5] T. Spies, J. Neu, U. Tayvah, M. Capobianco, B. Pattengale, S. Ostresh, and C. Schmuttenmaer, “Terahertz spectroscopy of emerging materials,” J. Phys. Chem. C, vol. 124, no. 41, pp. 22,335–22,346, 2020, doi: 10.1021/acs.jpcc.0c06344.
[6] R. Stantchev, X. Yu, T. Blu, and E. Pickwell-MacPherson, “Real-time terahertz imaging with a single-pixel detector,” Nature Commun., vol. 11, no. 1, pp. 1–8, 2020, doi: 10.1038/s41467-020-16370-x.
[7] M. S. Heimbeck and H. O. Everitt, “Terahertz digital holographic imaging,” Adv. Opt. Photon., vol. 12, no. 1, pp. 1–59, 2020, doi: 10.1364/AOP.12.000001.
[8] A. C. Kak, “Algorithms for reconstruction with nondiffracting sources,” in Principles of Computerized Tomographic Imaging, A. C. Kak and M. Slaney, Eds. Philadelphia, PA, USA: SIAM, 2001, pp. 49–112.
[9] H. Guerboukha, K. Nallappan, and M. Skorobogatiy, “Toward real-time terahertz imaging,” Adv. Opt. Photon., vol. 10, no. 4, pp. 843–938, 2018, doi: 10.1364/AOP.10.000843.
[10] N. A. Aghamiri, F. Huth, A. J. Huber, A. Fali, R. Hillenbrand, and Y. Abate, “Hyperspectral time-domain terahertz nano-imaging,” Opt. Express, vol. 27, no. 17, pp. 24 231–24 242, 2019, doi: 10.1364/OE.27.024231.
[11] L. Zanotto , R. Piccoli, J. Dong, R. Morandotti, and L. Razzari, “Single-pixel terahertz imaging: A review,” Opto-Electron. Adv., vol. 3, no. 9, p. 09200012, 2020, doi: 10.29026/oea.2020.200012.
[12] C.-H. Lin, Y.-C. Hung, F.-Y. Wang, and S.-H. Yang, “Hyperion: Hyperspectral penetrating-type ellipsoidal reconstruction for terahertz blind source separation,” 2021, arXiv:2109.05425.
[13] M. Skorobogatiy, J. Sadasivan, and H. Guerboukha, “Statistical models for averaging of the pump–probe traces: Example of denoising in terahertz time-domain spectroscopy,” IEEE Trans. THz Sci. Technol., vol. 8, no. 3, pp. 287–298, May 2018, doi: 10.1109/TTHZ.2018.2814820.
[14] Y. Chen, S. Huang, and E. Pickwell-MacPherson, “Frequency-wavelet domain deconvolution for terahertz reflection imaging and spectroscopy,” Opt. Express, vol. 18, no. 2, pp. 1177–1190, 2010, doi: 10.1364/OE.18.001177.
[15] X. Shen, C. R. Dietlein, E. Grossman, Z. Popovic, and F. G. Meyer, “Detection and segmentation of concealed objects in terahertz images,” IEEE Trans. Image Process., vol. 17, no. 12, pp. 2465–2475, Dec. 2008, doi: 10.1109/TIP.2008.2006662.
[16] H. Park and J.-H. Son, “Machine learning techniques for THz imaging and time-domain spectroscopy,” Sensors, vol. 21, no. 4, p. 1186, Feb. 2021, doi: 10.3390/s21041186.
[17] K. Kawase, Y. Ogawa, Y. Watanabe, and H. Inoue, “Non-destructive terahertz imaging of illicit drugs using spectral fingerprints,” Opt. Express, vol. 11, no. 20, pp. 2549–2554, 2003, doi: 10.1364/OE.11.002549.
[18] L. Zanotto, R. Piccoli, J. Dong, D. Caraffini, R. Morandotti, and L. Razzari, “Time-domain terahertz compressive imaging,” Opt. Express, vol. 28, no. 3, pp. 3795–3802, 2020, doi: 10.1364/OE.384134.
[19] T. Nagatsuma, H. Nishii, and T. Ikeo, “Terahertz imaging based on optical coherence tomography,” Photon. Res., vol. 2, no. 4, pp. B64–B69, 2014, doi: 10.1364/PRJ.2.000B64.
[20] J. Li et al., “Spectrally encoded single-pixel machine vision using diffractive networks,” Sci. Adv., vol. 7, no. 13, p. eabd7690, 2021, doi: 10.1126/sciadv.abd7690.
[21] X. Yang, D. Zhang, Z. Wang, Y. Zhang, J. Wu, B. Wu, and X. Wu, “Super-resolution reconstruction of terahertz images based on a deep-learning network with a residual channel attention mechanism,” Appl. Opt., vol. 61, no. 12, pp. 3363–3370, 2022, doi: 10.1364/AO.452511.
[22] W.-T. Su, T.-H. Chao, S.-H. Yang, and C.-W. Lin, “Seeing through a black box: Toward high-quality terahertz tomographic imaging via multi-scale spatio-spectral image fusion,” in Proc. Eur. Conf. Comput. Vis., Oct. 2022.
[23] Y.-C. Hung, T.-H. Chao, P. Yu, and S.-H. Yang, “Terahertz spatio-temporal deep learning computed tomography,” Opt. Express, vol. 30, no. 13, pp. 22,523–22,537, 2022, doi: 10.1364/OE.461439.
[24] D. Popescu and A. Hellicar, “Point spread function estimation for a terahertz imaging system,” EURASIP J. Adv. Signal Process., vol. 2010, no. 1, p. 575,817, 2010, doi: 10.1155/2010/575817.
[25] T. Wong, M. Kahl, P. Bolívar, and A. Kolb, “Computational image enhancement for frequency modulated continuous wave (FMCW) THz image,” J. Infrared, Millimeter, THz Waves, vol. 40, no. 7, pp. 775–800, 2019, doi: 10.1007/s10762-019-00609-w.
[26] M. Ljubenovic, S. Sazrafkan, J. D. Beenhouwer, and J. Sijbers, “CNN-based deblurring of terahertz images,” in Proc. 15th Int. Conf. Comp. Vis. Theory Appl., 2020, pp. 323–330.
[27] D. M. Slocum, E. J. Slingerland, R. H. Giles, and T. M. Goyette, “Atmospheric absorption of terahertz radiation and water vapor continuum effects,” J. Quantitative Spectrosc. Radiative Transfer, vol. 127, pp. 49–63, Sep. 2013, doi: 10.1016/j.jqsrt.2013.04.022.
[28] O. Ronneberger, P. Fischer, and T. Brox, “U-Net: Convolutional networks for biomedical image segmentation,” in Proc. Int. Conf. Med. Image Comput. Comput.-Assisted Intervention, 2015, pp. 234–241, doi: 10.1007/978-3-319-24574-4_28.
[29] K. Zhang, W. Zuo, Y. Chen, D. Meng, and L. Zhang, “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising,” IEEE Trans. Image Process., vol. 26, no. 7, pp. 3142–3155, Jul. 2017, doi: 10.1109/TIP.2017.2662206.
[30] S. Cheng, Y. Wang, H. Huang, D. Liu, H. Fan, and S. Liu, “NBNet: Noise basis learning for image denoising with subspace projection,” in Proc. IEEE/CVF Int. Conf. Comput. Vis. Pattern Recognit., 2021, pp. 4896–4906.
Digital Object Identifier 10.1109/MSP.2022.3198807