Computational Image Enhancement for Frequency Modulated Continuous Wave (FMCW) THz Image

Since its original inception in the early nineties [1], THz imaging has demonstrated a very large potential for contact-free analysis, nondestructive testing and stand-off detection in a wide variety of application fields, such as semiconductor industry, biology, medicine, material analysis, quality control, and security [2, 3, 4]. In many of these application fields, THz imaging is competing with established imaging methodologies, such as optical inspection or X-ray imaging. Compared with imaging in the optical or X-ray parts of the electromagnetic spectrum, THz imaging is significantly limited in its spatial resolution due to the substantially longer wavelength of the associated frequencies.

A wide range of technological approaches to realize THz imaging systems have been demonstrated, including femtosecond laser based scanning systems [1, 5], synthetic aperture systems [6, 7], and hybrid systems [8]. Several of these approaches using pulses or using frequency modulated THz signals, allow to sense temporal or phase shifts to the object’s surface, making 3D THz imaging possible.

Despite the fact that in most of these approaches THz imaging is performed close to the diffraction limit, the comparatively low spatial resolution associated with the THz radiation wavelengths significantly hampers the application range. There is a huge interest to increase the spatial resolution of this approach beyond the diffraction limit, in order to make this technique competitive in comparison to established methods such as X-ray imaging [9].

THz imaging below the diffraction limit is an emerging area [3], which can roughly be classified into two alternatives: by system enhancements or by computational approaches. System enhancements include for example, interferometric sensing [10] to increase the depth sensitivity in THz time-of-flight imaging, or near-field sensing approaches [11] which demonstrate a nanometer scale lateral resolution by circumventing the diffraction limit. Computational image enhancement techniques aim at improving the resolution by utilizing numerical procedures and additional signal or system information, e.g. from overlapping oversampled THz imaging signals, without introducing additional equipment effort and cost.

Depending on the THz image acquisition mode, there are several alternative approaches for computational image enhancement. THz imaging superresolution (also referred to as high-resolution or image restoration) is often associated to spatial resolution enhancement in the xy-direction [12, 13, 14, 15, 16]. In contrast, depth resolution enhancement is associated to improvement in azimuth direction (z-direction) [17, 18, 19].

In this paper, we propose a novel method to enhance THz image resolution beyond diffraction limit, attaining a lateral (xy) resolution increase and a depth (z) accuracy increase. The concept is demonstrated for a Frequency Modulated Continuous Wave (FMCW) THz scanning system operating at $514-640$GHz, but can certainly be adapted to other THz imaging techniques. Our approach comprises the following technical contributions:

• •

  Parameter extraction by complex signal fitting model in z-direction for each pixel which allows for the acquisition of non-planar targets, incorporating

  • –

    an accurate estimation of the per-pixel distance to the object surface, and

  • –

    a proper reconstruction of the reflection intensity as a measure for the object’s material-properties.

• •

  Based on the accurately reconstructed reflection intensity, we are able to apply more complex, state-of-art 2D blind deconvolution techniques in order to improve the spatial xy-resolution beyond what is achievable with traditional (e.g. gaussian kernel) deconvolution procedures.

Sections 3 and 4 give a brief overview of the THz 3D imaging system and the proposed method, respectively. The details of the curve fitting procedure are described in Sec. 5. In Sec. 7, the evaluation of the computational result of the proposed method are depicted.

The majority of prior research focuses on the lateral resolution of 2D THz images, where the Lucy-Richardson deconvolution algorithm [20, 21] is one of the most frequently used methods. Knobloch et al. [22] firstly applied Lucy-Richardson deconvolution on THz images. Xu et al. [12] proposed a THz time-domain spectroscopy (THz-TDS) image high-resolution reconstruction model incorporating a 2D wavelet decomposition and a Lucy-Richardson deconvolution to reconstruct a high-resolution THz image from four low-resolution images and to reconstruct a high-resolution image from single degraded 2D low-resolution image. Li et al. [13] proposed to use the Lucy-Richardson deconvolution algorithm for a coherent THz 2D imaging system. Ding et al. [14] used the Lucy-Richardson deconvolution for a THz reflective 2D imaging system.

In addition to a 2D deconvolution algorithm, Hou et al. [15] proposed a method to enhance THz image quality by applying a Finite Impulse Response Filter in time domain. Ahi and Anwar [16] proposed a deconvolution method based on their THz image generation model for a THz far-field 2D imaging system. Xie et al. [23] proposed to use a Markov Random Field (MRF) model for THz superresolution.

For the THz blur kernel estimation, Ding et al. [7] proposed to use a Range Migration Algorithm (RMA). Gu et al. [24] proposed to use a physical model for THz point spread function (PSF) as a Gaussian kernel based on the optical parameters of a THz-TDS system.

To enhance the depth accuracy, Walker et al. [17] proposed a time domain deconvolution method to increase the time domain zero padding factor for a THz reflection imaging. Chen et al. [18] proposed a hybrid Frequency-Wavelet Domain Deconvolution (FWDD) method to improve calculation of impulse response functions. Takayanagi et al. [19] proposed a deconvolution method based on Wiener filtering in the time domain. Dong et al. [25] proposed a sparse deconvolution method of impulse response function on multi-layered structure.

In the context of THz 3D imaging reconstruction method Shift Migration (PSM) algorithm to reconstruct a 3D image, Sun et al. [26] proposed to extend PSM as Enhanced Phase Shift Migration (EPSM) to improve computational efficiency. Liu et al. [27] proposed a 3D Wavenumber Scaling Algorithm (3D WSA) to reconstruct the entire 3D holographic image.

To our best knowledge, the method in this paper is the first attempt that combines depth accuracy with lateral resolution enhancement in order to achieve a jointly improved spatial resolution and accuracy in both, $xy$- and $z$-direction.

Our electronic THz imaging system is based on hollow-waveguide multipliers and mixers, operating in a frequency modulated continuous wave (FMCW) mode for measuring depth information. The components are operating around a frequency of $577$GHz with a bandwidth of $126$GHz. More details on the experimental approach are described in[7].

Fig. 1 shows the imaging geometry. Both transmitter (Tx) and receiver (Rx) are mounted on the same platform. The imaging unit, consisting of Tx, Rx and optical components, are moved along the x and y direction using stepper motors and linear stages. This imaging unit takes a depth profile of the object at each lateral position, in order to acquire a full 3D image. The data is acquired with a lateral step size of $262.5\mu m$ in xy-direction. During measurement, the motor controller and the data acquisition are synchronized to enable on the fly measurements. An adequate integration time and velocity is chosen in order to provide enough time for the acquisition of 1400 samples per depth profile and 36 averages per sample. The total acquisition time for such an averaged depth profile is $5$ms.

Transmitter (Tx) and receiver (Rx) are mounted in a monostatic geometry. A beam splitter and two hyperbolic lenses focus the beam radiated from the Tx and reflected from the sample into the Rx. By measuring minimum dimension which obtains more than 3dB intensity difference, the resolution of the setup is measured as $793.7\mu m$ using a metallic USAF 1951 Resolving Power Test Target scaled to the THz frequency range, which is close to the theoretical ideal expectation of $622\mu m$ (numerical aperture NA$=0.508$ @ $578$GHz). Due to the monostatic scanning approach no optical magnification occurs. The system’s depth resolution $\Delta d$ is defined by

where $c$ is the speed of light in air, and $B$ is the system bandwidth. For our system, we have $\Delta d=1210\mu m$

For the FMCW operation, a voltage controlled oscillator (VCO) is tuned from $14.28-17.78$GHz (see Fig. 2). The signal at the output of the VCO is distributed to the Tx and the Rx using a power-splitter. For transmission the signal is then upconverted to $514-640$GHz with a chain of multipliers. After downconversion to the intermediate frequency (IF) range using a sub-harmonic mixer, which is fed with the 6th harmonic of the VCO signal, the signal is digitized with 10 MS (mega-samples) per second sampling rate and transferred to the host PC. An additional delay between the Tx and Rx path creates a frequency offset in the intermediate frequency signal for proper data acquisition.

In this part, we apply per-pixel parameter extraction in z-direction in order to represent the measured complex signal by a complex reflection model. As each pixel is treated independently, we simplify notation by dropping the pixel-location using, e.g., $\hat{u}[z_{i}]$ for $\hat{u}(x,y)[z_{i}]$.

Our FMCW-THz 3D imaging system is calibrated by amplitude normalization with respect to an ideal metallic reflector. Thus, we achieve an ideal rectangular frequency amplitude signal response, which, after being Fourier transformed, results in an ideal $\operatorname*{sinc}$ function $v(z)$ as continuous spatial signal amplitude. In this model, we assume single layer reflection, and the extension to multi-layer reflection is part of future work. As we assume a single layer reflection, the complex signal model $v(z)$ is modeled as a modulated $\operatorname*{sinc}$ function:

In Eq. (4), $A$ is the electric field amplitude, $\mu$ and $\sigma$ are the mean (i.e. the depth) and the width of the $\operatorname*{sinc}$ function, respectively, $\omega$ is the angular frequency of the sinusoidal carrier and $\phi$ is the depth-related phase shift. We formulate the parameter extraction as a complex curve fitting process with respect to optimizing and minimizing the energy function $\varepsilon_{c}$

where the subscripts ${re}$ and ${im}$ denote the real and the imaginary part of a complex number, respectively, and $Z_{f}$ is the fitting window (see Sec. 5.1).

Because of the highly non-linear optimization involved in the curve fitting process, a direct application of Eq. (5) to a non-linear solver potentially results in local minima and does not lead to robust results. Therefore, we apply the following optimization steps which achieve a robust complex curve fitting:

1. 1.

   Estimate the signal’s maximum magnitude z-position $z_{m}$ in order to localize the curve fitting window (see Sec. 5.1).

2. 2.

   Apply a curve fitting to the magnitude signal leading to initial values for $A_{m},\mu_{m},\sigma_{m}$ (see Sec. 5.2).

3. 3.

   Estimate the initial phase value $\phi_{m}$ using a phase matching with respect to the angle of complex signal $\angle\hat{u}$ (see Sec. 5.3).

4. 4.

   Based on the initial values $A_{m},\mu_{m},\sigma_{m},\phi_{m}$ the optimization is performed by minimizing the energy $\varepsilon_{c}$ (see Eq. (5)) using the Trust Region Algorithm [30].

5. 5.

   Reconstruct an intensity image $I_{v}$ and an depth image $d_{v}$ based on the curve fitting result (sec Sec. 5.4).

Our curve fitting approach significantly reduces the per-pixel intensity inhomogeneity (see Sec. 7.2.2). The subsequent lateral deconvolution algorithm discussed in Sec. 6 involves the numerical solution of an ill-posed inverse problem of finding the blur kernel and enhancing the image’s sharpness at the same time, which is very sensitive to noise. Therefore, correcting the intensity yields two advantages, i.e., it stabilizes the numerical deconvolution process and it prevents wrong interpretations of intensity variation as structural or material transitions. We will provide comparison in Fig. 13 and discuss the lateral resolution enhancement with and without the per-pixel parameter extraction in Sec. 7.3.1.

After the curve fitting, the reconstructed intensity $I_{v}$ is a 2D image with real and positive values. We can now apply a state-of-art 2D deconvolution algorithm to enhance the xy-domain resolution. In contrast to prior work we use a total variation (TV) blind deconvolution algorithm to improve the spatial resolution [31, 32].

In deconvolution image model, our reconstructed intensity image $I_{v}$ can be expressed as a blurred observation of a sharp image $I_{d}$

where $h$ is the spatially invariant point spread function (PSF) (also known as blur kernel) and $\eta$ is the noise; $\ast$ denotes convolution. Blind deconvolution methods allow for estimating the blur kernel directly from the data, which is, however, an ill-posed inverse problem that requires a prior knowledge in order to deduce a robust result [33]. In this paper, we utilize the sparse nature of intensity gradients [34] and choose a state-of-art TV blind deconvolution algorithm that minimizes

Here, $\|I_{d}\ast h-I_{v}\|_{1}$ is commonly referred as a data term, $\lambda$ is a regularization parameter and $\|\nabla I_{d}\|_{1}$ is the TV-regularization (or prior) that enforces the gradient of the resulting deblurred image $I_{d}$ to be sparse. As by-product, the blind deconvolution yields the estimated PSF $h$. We obtain our final results using the implementation from Xu et al. [31, 32].

In this section we evaluate our image enhancement approach using the following two datasets

1. 1.

   MetalPCB: A nearly planar “USAF” target etched on a circuit board (Fig. 3). The dataset has been acquired using the setup described in Sec. 3 and has the resolution $N_{x}=446$, $N_{y}=446$, $N_{z}=1400$. The lateral per-pixel distance (i.e. the sensor mechanical movement distance) is $262.5{\mu}m$. The target is etched in the standard size scale of USAF target MIL-STD-150A.

2. 2.

   StepChart: A metallic step chart with steps varying from $4000\mu m$ to $50\mu m$, and a reference plane to locate the reference zero position (Fig. 5). The dataset has also been acquired using the setup in Sec. 3 with the resolution $N_{x}=575$, $N_{y}=113$, $N_{z}=1400$. The lateral per-pixel distance is $262.5{\mu}m$.

It should be noted, that our system has a lateral resolution of $622\mu m$ as its ideal diffraction limit, an experimentally measured lateral resolution $793.7\mu m$ and a depth accuracy of $\Delta d=1210\mu m$ as its diffraction limits (see Sec. 3).

We apply and compare the following images, reconstruction and image enhancement methods. The images are named as $\langle image\rangle\_\langle method\rangle\mbox{-}\langle kernel\rangle$, where $image$ is the intensity image applying deconvolution methods, $method$ is the method applied to the intensity image, and $kernel$ is the deblur kernel adopted in the deconvolution method.

1. 1.

   Refer: The reference intensity image $I_{u}$; no reconstruction using curve fitting and no image enhancement applied.

2. 2.

   Reconst: The reconstructed intensity image $I_{v}$ using curve fitting; no image enhancement applied.

3. 3.

   Refer_Xu: image enhancement using Xu et al. [31, 32] applied to reference intensity Refer.

4. 4.

   Reconst_Xu: image enhancement using Xu et al. [31, 32] applied to reconstructed intensity Reconst. This is our proposed method.

5. 5.

   Refer_LR-G: using Lucy-Richardson [20, 21] with gaussian kernel applied to reference intensity Refer.

6. 6.

   Reconst_LR-G: using Lucy-Richardson [20, 21] with gaussian kernel applied to reconstructed intensity Reconst.

7. 7.

   Refer_LR-Xu: using Lucy-Richardson [20, 21] with sparse kernel extracted from Xu et al. [31, 32] applied to reference intensity Refer.

8. 8.

   Reconst_LR-Xu: using Lucy-Richardson [20, 21] with sparse kernel extracted from Xu et al. [31, 32] applied to reconstructed intensity Reconst.

In the following, we will first deduce the optimal window size for the quality control of the fitting using the MetalPCB dataset (Sec. 7.1). In Sec. 7.2, we will discuss the depth accuracy using the StepChart dataset. In Sec. 7.3, we will discuss the lateral resolution on MetalPCB dataset. All intensity images are normalized to a perfect metal reflection, and displayed using MATLAB’s perceptionally uniform colormap parula [35].

In this paper, we propose a THz computational image enhancement method to enhance the lateral resolution and depth accuracy beyond the diffraction limit. The method is based on a complex curve fitting on THz 3D image azimuth axis, and a blind deconvolution method on the lateral domain. The experiment results show that our curve fitting method enhances the depth accuracy to $91\mu m$. Because of this enhanced depth accuracy, the experiment also shows that the proposed reconstruction method achieves improved intensities on homogeneous, non-planar material regions without introducing additional distortions or noise.

Based on the reconstructed intensities, we apply several lateral blind deconvolution methods. Evidently, all three examined THz deconvolution methods benefit from the enhanced intensity image using our method. In comparison to the classical Lucy-Richardson deconvolution algorithm, the experiments show that the proposed blind deconvolution method achieves the best horizontal resolution enhancement from $794.3\mu m$ to $346.2\mu m$, yielding an improvement factor of $2.29$. In terms of intensity contrast, the proposed method clearly outperforms earlier approaches. Moreover, the experiments show that the proposed method is able to emphasize fine silk texture embedded within a polymer material.

This research was funded by the German Research Foundation (DFG) as part of the research training group GRK 1564 Imaging New Modalities.

This is a pre-print of an article published in Journal of Infrared, Millimeter, and Terahertz Waves. The final authenticated version is available online at: https://doi.org/10.1007/s10762-019-00609-w