Layer-Resolving Terahertz Light-Field Imaging Based on Angular Intensity Filtering Method

Abstract

Terahertz focal plane array imaging methods, direct camera imaging and conventional light field imaging methods are incapable of resolving and separating layers of multilayer objects. In this paper, for the purpose of fast, high-resolution and layer-resolving imaging of multilayer structures with different reflection characteristics, a novel angular intensity filtering (AIF) method based on terahertz light-field imaging is purposed. The method utilizes the extra dimensional information from the 4D light field and the reflection characteristics of the imaging object, and the method is capable to resolve and reconstruct layers individually. The feasibility of the method is validated by experiment on both “idealized” and “practical” multilayer samples, and the advantages in performance of the method are proven by quantitative comparison with conventional methods.

1. Introduction

A terahertz wave is an electromagnetic wave with wavelengths between 30–3000 μm or a frequency between 0.1–10 THz. Thanks to the unique characteristics in high penetration against dielectric materials and low electron energy [1,2,3], the applications of terahertz imaging techniques in the fields in which see-through imaging and inspection are demanded, such as industrial non-destructive testing (NDT) [3,4,5], heritage conservation [6,7,8], security [9,10,11] and aerospace [12], have been widely concerned and studied. For the NDT requirements of multilayer structures with complex optical characteristics and potential defect forms, terahertz imaging methods have advantages in efficiency and flexibility over conventional NDT methods such as thermography and ultrasound imaging [13,14]. The improvements in strong terahertz source and compact terahertz detectors with high sensitivity techniques enabled the development of high-performance terahertz imaging methods, namely, the terahertz pulsed imaging (TPI) [15] and terahertz focal plane array (FPA) imaging [16,17,18].

TPI is a point-by-point imaging method based on terahertz time domain spectroscopy (TDS) system, using a pulsed terahertz beam. Terahertz pulses generated with ultrafast optic methods contain an extremely wide frequency spectrum and rich information in both time and frequency domain. The transversal structures and optical characteristics of the objects could be resolved by applying signal processing methods such as deconvolution and wavelet transformation on the terahertz waveform. Thus, TPI is capable of resolving and reconstructing the transversal structures of the object with a high resolution, or support fine analysis on potential defects [19,20,21]. However, The TPI method relies on point-to-point raster scanning to acquire the lateral information, which limited the effectiveness and further application potential of the TPI method.

The improvement in high-power terahertz beams and integrated terahertz detector technologies enabled fast speed, large field of view and high lateral resolution imaging with terahertz FPA technology. In the application fields demanding inspection for large scale objects, such as marine antifoul and anticorrosion coating, terahertz FPA imaging is a perfect complement for TPI technology. However, for multi-layered objects with complex structures and optical characteristics, on the one hand, conventional FPA methods have no depth-resolving capabilities, and are incapable of resolving features of different layers; on the other hand, the interested features could be submerged by reflections of other layers and unable to be extracted and reconstructed, which limited the capability of conventional FPA imaging methods.

By taking advantages of the polarization characteristics of reflected or transmitted waves, polarization imaging methods are widely studied in the terahertz waveband for revealing hidden features of the object and improving the signal-to-noise ratio (SNR) performances [22,23,24]. However, when the optical configurations are restricted, polarization imaging methods are unable to solve the problems of layer-resolving and submerged features.

For problems mentioned above, the light field imaging (LFI) technique is a promising solution. As a novel computation imaging method based on incoherent FPA detectors, light field imaging (LFI) method resolve depth, radiometry and spatial information by acquiring and utilizing 4D light field information of the object [25]. In visible wavebands, the improvements in integrated array detectors enabled practical 4D light field acquisition methods, which allowed LFI technique applying [26] in multiple application fields, such as dynamic refocusing, depth of field (DoF) extension [27], 3D reconstruction [28,29,30], super-resolution reconstruction [31,32], synthetic aperture (SA) imaging and high-speed video capture [33]. In infrared and terahertz waveband, several studies on the LFI technique have also been reported [34,35].

In our previous works, by introducing the LFI Synthetic Aperture (LFI-SA) technique, the distortion and resolving power issues of conventional active illumination FPA imaging systems has been solved. By utilizing the full 4D light field information, it is promising to achieve layer resolving, separation and reconstruction as well.

In this paper, a novel angular intensity filtering (AIF) method based on terahertz light-field imaging for fast, high-resolution and layer-resolving imaging of multilayer structures with different reflection characteristics is purposed. The AIF method resolves and reconstructs the layers beneath the surface by utilizing the extra dimensional information from the 4D light field and the reflection characteristics of the imaging object, which are incapable to distinguish with direct FPA imaging and conventional LFI-SA methods. A reflection-type active illumination terahertz LFI system based on FPA detectors has been established and imaging experiments of both “idealized” multilayer sample and “practical” marine antifoul coating sample with a rusted layer have been initiated, and the feasibility and performance of our AIF method is validated.

2. Related Works

The depth resolving capability makes LFI a promising tomography method when working together with the penetration and non-destruction features of terahertz wave. In visible waveband, light field tomography based on light field microscopy system is proven feasible [36,37]. However, for typical light field imaging systems, such as hand-held light field cameras or camera arrays, the depth resolution of conventional LFI depth estimation algorithms is limited by the angular resolution of the LFI system, which makes it unfeasible to resolve thin multilayer materials directly with a typical light field imaging system [33].

Moreover, conventional LFI depth estimation algorithms based on defocusing and correspondence largely based on the assumption of ideal Lambertian model and the simplification of wave optic effects [38]. For terahertz imaging systems based on active illumination, the characteristics of light source, optics and detectors made the robustness of conventional methods unacceptable for terahertz LFI results [39]. For this reason, the resolving and imaging of such structures relies on the differences on radiometry characteristics of layers.

In visible wavebands, the LFI-based specular removal algorithms are applied to separate light wave components with different radiometry characteristics from the objects which are capable to resolve semi-transparent, overlapped and layered objects, and demonstrated advantages in performance and flexibility compared with conventional methods.

Tao et al. purposed two different specular separation and removal algorithms based on LFI. The line consistency-based algorithm distinguishes the reflection characteristics of the object according to the line consistency in the RGB space, which seems not to be compatible with terahertz images without multi-color channel information [40]. The k-means clustering specular removal algorithm distinguishes the specular and diffuse properties of light with the bidirectional reflection distribution function (BRDF) model and k-means clustering specular removal algorithm and extracts the specular properties while re-evaluating and re-clustering the refocused pixels iteratively [41]. The latter algorithm is applied to terahertz LFI images and compared with our results in Section 4.2 and Section 4.3.

In the terahertz waveband, several studies on the feasibility and application of LFI have been started. Jain et al. verified the feasibility of LFI in the terahertz waveband with an experimental virtual camera array system [35]. In previous works, we established a practical active illumination terahertz LFI system based on gas-pumped terahertz sources and terahertz FPA detectors, and verified the performance improvements of terahertz LFI in resolving power and SNR performance compared with conventional methods with transparent-type LFI experiments [39].

3. Theory and Methods

3.1. The BRDF Model and Reflection Properties of Surfaces

The differences in reflection characteristics can be described with the dichromatic bidirectional reflectance distribution function (BRDF) model [40,41]. The dichromatic BRDF model is used for describing the angular-independent diffuse and angular-dependent specular properties of reflection. In the dichromatic BRDF model, the electromagnetic wave reflected from a certain object surface can be expressed as:

(1)$f\left(\lambda ,\mathsf{\Theta }\right)=g_d\left(\lambda \right)f_d+g_s{f}_s\left(\mathsf{\Theta }\right)$

where the $f\left(\lambda ,\mathsf{\Theta }\right)$ represents the reflection multipliers of the object at wavelength $\lambda$ and viewing angle $\mathsf{\Theta }$, $g_d$ and $g_s$ represent the spectral reflectance of diffuse and specular reflection at wavelength $\lambda$, $f_d$ and $f_s$ represents the reflection multipliers of diffuse and specular reflection at viewing angle $\mathsf{\Theta }$. Considering the $g_s$ is wavelength-independent for the dielectric materials and the $f_d$ is angle-independent for the diffuse properties, the formula can be simplified as:

(2)$f_{\lambda }\left(\mathsf{\Theta }\right)=g_{\lambda }{}_d{f}_d+g_{\lambda }{}_s{f}_s\left(\mathsf{\Theta }\right)$

Additionally, the reflection intensity at viewing angle $\mathsf{\Theta }$ can be expressed as:

(3)$I\left(\mathsf{\Theta }\right)={\displaystyle \sum }{I}_{L_k}\left(g_{\lambda }{}_d{f}_d+g_{\lambda }{}_s{f}_s\left(\mathsf{\Theta }\right)\right)·\stackrel{\rightharpoonup }{n}·{\stackrel{\rightharpoonup }{l}}_k$

where the $I_{L_k}$ represents the intensity of the light source $k$, the $\stackrel{\rightharpoonup }{n}$ represents the direction of the surface normal, and the ${\stackrel{\rightharpoonup }{l}}_k$ represents the incident direction of the light source k. When there is only one single light source, the formula can be simplified as:

(4)$I\left(\mathsf{\Theta }\right)=I_L\left(g_{\lambda }{}_d{f}_d+g_{\lambda }{}_s{f}_s\left(\mathsf{\Theta }\right)\right)·\stackrel{\rightharpoonup }{n}·\stackrel{\rightharpoonup }{l}$

When a narrow terahertz beam reflects on a certain surface, if the surface is relatively specular, it has a smaller $g_{\lambda }{}_d$, and a larger $g_{\lambda }{}_s$, which leads to a sharper $f_s\left(\mathsf{\Theta }\right)$ distribution and indicates a narrower and more concentrated reflection, as shown in Figure 1a. If the surface is relatively less specular, it has a smaller $g_{\lambda }{}_d$ and a larger $g_{\lambda }{}_s$, which leads to a flatter $f_s\left(\mathsf{\Theta }\right)$ distribution and indicates a more even diffusion, as shown in Figure 1b. Such differences in BRDF characteristics can be used to distinguish the reflected beam from surfaces of different types.

Active illumination terahertz imaging relies on collimated and concentrated beams as illumination. When a collimated beam interacts with materials with different characteristics, the differences in intensity and scattering of the reflected beam would be significant.

Figure 2 shows the BRDF distribution measurement results of typical glossy and rough surfaces in a terahertz waveband. Because of the wavelength characteristics, the coherence of the light source, and the performances of the detectors, the measured BRDF distribution shows different characteristics compared with the BRDF distribution of same types of materials in visible waveband. For example, the $f_d$ of diffuse and scattered surfaces are not exactly flat and the $f_s$ of glossy surfaces are not unimodal. However, the BRDF distribution of different materials still shows significant differences, which are capable of being distinguished by algorithms.

3.2. The Interaction between Terahertz Beam and Multilayer Materials

We now consider the case of multilayer structures consisting of layers with different reflection characteristics. The multilayer structure consists of two layers: the surface layer and the base layer. The surface layer is specular, uniform and transparent. It reflects and transmits the incident terahertz wave, and it has a high $g_s$, an ignorable $g_d$ and a sharper $f_s\left(\mathsf{\Theta }\right)$ distribution. The base layer is relatively less specular. It scatters and absorbs the incident terahertz wave, and it has a higher $g_d$, a much lower $g_s$ and a flatter $f_s\left(\mathsf{\Theta }\right)$ distribution.

When collimated and concentrated terahertz waves is incident on the multilayer structure with a small incident angle, most of the energy reflects at the surface sharply and the remaining energy is transmitted through the surface layer, then scattered and reflected by the base layer, as shown in Figure 3. When the loss of the surface layer is ignored, the BRDF distribution of the multilayer structure can be expressed as:

(5)$f_{\mathrm{total}}\left(\mathsf{\Theta }\right)=f_{\mathrm{surface}}\left(\mathsf{\Theta }\right)+f_{\mathrm{base}}\left(\mathsf{\Theta }\right)=g_{{\mathrm{surface}}_{\mathrm{d}}}{f}_{{\mathrm{surface}}_{\mathrm{d}}}+g_{{\mathrm{surface}}_{\mathrm{s}}}{f}_{{\mathrm{surface}}_{\mathrm{s}}}\left(\mathsf{\Theta }\right)+g_{{\mathrm{base}}_{\mathrm{d}}}{f}_{{\mathrm{base}}_{\mathrm{d}}}+g_{{\mathrm{base}}_{\mathrm{s}}}{f}_{{\mathrm{base}}_{\mathrm{s}}}\left(\mathsf{\Theta }\right)$

where $f_{\mathrm{surface}}\left(\mathsf{\Theta }\right)$ and $f_{\mathrm{base}}\left(\mathsf{\Theta }\right)$ represent the BRDF distribution of the surface layer and the base layer at certain viewing angle $\mathsf{\Theta }$, respectively. When the $g_{{\mathrm{surface}}_{\mathrm{d}}}$ can be ignored, and the $f_{{\mathrm{base}}_{\mathrm{s}}}\left(\mathsf{\Theta }\right)$ is flat, the above formula can be simplified as:

(6)$f_{\mathrm{total}}\left(\mathsf{\Theta }\right)\approx g_{{\mathrm{surface}}_{\mathrm{s}}}{f}_{{\mathrm{surface}}_{\mathrm{s}}}\left(\mathsf{\Theta }\right)+f_{\mathrm{base}}$

which means, among the total reflected wave, only the $f_{\mathrm{surface}}$ varies with the viewing angle, whereas the $f_{\mathrm{base}}$ does not.

Thus, the reflected terahertz wave from multilayer structure consists of $f_{\mathrm{surface}}$ property from the surface layer, and the $f_{\mathrm{base}}$ property from the base layers. In the case that $f_{\mathrm{base}}$ is not higher than $f_{\mathrm{surface}}$, on the one hand, the flatter $f_{{\mathrm{base}}_{\mathrm{s}}}\left(\mathsf{\Theta }\right)$-made $f_{\mathrm{base}}$ is lower than the peak of $f_{{\mathrm{base}}_{\mathrm{s}}}$ on any direction; on the other hand, the absorption of the base layer further attenuated the $f_{\mathrm{base}}$ on all directions, which resulted in the reflection intensity of the base and base layer being significantly weaker than the reflection intensity of the surface layer.

In practical situations, the reflection from the base layer comes from the remaining component of the incident beam after the reflection and absorption of the surface layer. When considering the loss of the surface layer, the reflection from the base layer can be expressed as:

(7)$f_{\mathrm{total}}\left(\mathsf{\Theta }\right)=f_{\mathrm{surface}}\left(\mathsf{\Theta }\right)+\left(1-A_{\mathrm{reflection}}\right){\left(1-A_{\mathrm{absorption}}\right)}^2{f}_{\mathrm{base}}$

in which the $A_{\mathrm{reflection}}$ is the reflection loss of the surface layer and the $A_{\mathrm{absorption}}$ is the absorption loss of the surface layer. Correspondingly, the $f_{\mathrm{base}}$ could be obtained by:

(8)$f_{\mathrm{base}}=\frac{f_{\mathrm{total}}\left(\mathsf{\Theta }\right)-f_{\mathrm{surface}}\left(\mathsf{\Theta }\right)}{\left(1-A_{\mathrm{reflection}}\right){\left(1-A_{\mathrm{absorption}}\right)}^2}=\frac{f_{\mathrm{total}}\left(\mathsf{\Theta }\right)-f_{\mathrm{surface}}\left(\mathsf{\Theta }\right)}{\left(1-\int f_{\mathrm{surface}}\left(\mathsf{\Theta }\right)d\mathsf{\Theta }\right){\left(1-A_{\mathrm{absorption}}\right)}^2}$

When using the conventional imaging method with FPA detectors, the reflections from different layers are received indistinguishably and the characteristics of the base layer are submerged in the reflection of surface layer, and make it impossible to distinguish the image result from either of the layers.

Changing the angle of the acquisition system and light source to avoid receiving from the direction of the strongest specular reflection seems to be a reasonable way to image from the beneath layers. However, considering that the angular distribution of real scattering of rough materials in terahertz waveband is still concentrated in a relatively small range, and the performance of terahertz FPA detectors, it is unable to achieve an acceptable SNR with this method.

3.3. The Angular Intensity Filtering Method

The angular intensity filtering (AIF) method is a method for resolving and reconstructing the images of individual layers beneath the surface from a multilayer structure, by utilizing the angular intensity-resolving capabilities of 4D LFI and the differences in BRDF characteristics of different layers. Specifically, on the basis of conventional LFI-SA methods, the AIF method generates the saturation map and specular map of the 4D light field information, distinguishes and filters the pixels from different layers and reconstructs the image of individual layers from the filtered 4D light field information.

As discussed in Section 3.1, considering the case of simplest dual-layered structure, for an input 4D light field $L_F$, the saturation map $L_{\mathrm{saturation}}$ can be expressed as:

(9)$L_{\mathrm{saturation}}\left(x,y,u,v\right)=L_F\left(x,y,u,v\right)·H\left(L_F\left(x,y,u,v\right)-T\right)\ast K_0$

where the $H$ is the Heaviside step function, the $T$ is the intensity threshold to distinguish the $f_{\mathrm{surface}}$ and $f_{\mathrm{base}}$ and the $K_0$ is the convolution kernel for pre-processing the image. The saturation map indicates the pixels directly infected by the specular properties of the surface layer. The specular map $L_{\mathrm{spec}}$ calculated from the $L_{\mathrm{saturation}}$ can be expressed as:

(10)$L_{\mathrm{spec}}\left(x,y,u,v\right)=H\left(\left(L_{\mathrm{saturation}}\ast K\right)\left(x,y,u,v\right)\right)$

where $K$ is the convolution kernel for filtering the pixels infected by the specular properties. The specular map indicates all the pixels directly and indirectly affected by the specular properties of the surface layer. Thus, the filtered reconstruction result of the 4D light field can be expressed as:

(11)$E_{\mathrm{surface}}\left(x_0,y_0\right)={\displaystyle \iint }{L}_{z\left(x_0,y_0\right)}\left(x^{\prime },y^{\prime },u,v\right)\cdot L_{\mathrm{spec}}\left(x’,y’,u,v\right)dudv$

(12)$E_{\mathrm{base}}\left(x_0,y_0\right)={\displaystyle \iint }{L}_{z\left(x_0,y_0\right)}\left(x^{\prime },y^{\prime },u,v\right)\cdot \left(1-L_{\mathrm{spec}}\left(x^{\prime },y^{\prime },u,v\right)\right)dudv$

where the $E_{\mathrm{surface}}$ and $E_{\mathrm{base}}$ are the reconstruction result of the surface layer and the base layer, respectively, and the $L_{z\left(x_0,y_0\right)}$ is the pixels contributed to the reconstruction of the position $\left(x_0,y_0\right)$. Thus, the resolving and imaging of individual layers discussed in Section 3.1 has been realized. The framework of the AIF algorithm implementation is shown in Figure 4.

4. Experiment and Results

4.1. Experiment Setup

A reflection-type active illumination terahertz 4D LFI system is established to acquire the 4D light field information. Specifically, the terahertz beam emitted from the source is collimated and illuminated on the imaging object, then the static 4D light field including both the specular and the diffuse properties is acquired by the virtual camera array system on the proper angle and position, as shown in Figure 5. The incident angle of terahertz beam onto the object is 20°. The camera plane of the virtual camera array is vertical to the reflected beam from the object. The distance between the source and the object is 1600 mm, and the distance between the object and the camera plane is 600 mm.

In the system, a coherent SIFIR-50 CO2-pumped continuous-wave terahertz source is used as a terahertz source, and an INO IRXCAM-THz-384 terahertz camera module is uses for acquisition. The SIFIR-50 source works on the frequency of 2.52 THz with an output power of 50 mW. The IRXCAM-THz-384 uses uncooled microbolometer FPA array detectors. The pixel size, resolution and frame rate of the detectors are 35 μm, 384 × 288 px and 50 fps, respectively. The camera module has a silicon optic lens of 44 mm focal length. To filtering the infrared wave in the environment, an additional Teflon filter is installed to the optics. The virtual camera array is simulated by using the camera module capturing light field slices at different positions of the camera plane. A 2D motorized translation stage is used for driving the camera module on the camera plane.

Two samples are designed and used in the experiments below. An “idealized” sample is used for validating the feasibility and effectiveness of the AIF method, and a “practical” marine antifoul coating sample is used for evaluating the performance of the AIF method.

4.2. The Feasibility Evaluation of AIF Method on Multilayer Resolving

In this section, the feasibility of the AIF method in correctness and image quality of reconstruction is validated by acquiring and reconstructing the light field data from the “idealized” sample, and compare the reconstruction results quantitatively.

An “idealized” multilayer sample with different BRDF characteristics is used as a sample for the experiments. The multilayer sample is composed of tightly stacked thin silicon wafer slices as the surface layer, and a rusted metal plate as the base layer, as shown in Figure 6. The silicon wafer slice is 200 μm in thickness, polished on both sides, and transflective to terahertz waves with ignorable scatter and diffuse. The rusted metal plate has a rusted layer on the surface, which scatters and diffuses the terahertz beam. The thickness of the rusted layer is about 40 μm in average. The base layer of the sample is individually used as reference sample to generate the specular-free reference image. The number of sub-images acquired in the experiment is 11 × 31, and the spatial interval is 5 mm.

The reconstruction methods are used for reconstructing the base layer of the “idealized” sample, and the reconstruction result is compared with the reference image acquired and reconstructed from the reference sample to evaluate the performances of the methods. Ideally, the reconstruction results should match the reference image.

Here, the conventional LFI-SA method and k-means clustering specular removal algorithm proposed by Tao et al. are used as comparison of the AIF method. Tao’s algorithm is optimized for light fields of visible waveband and RGB color space, but it is still available in a terahertz waveband.

For the reference image, a Gaussian noise addition of σ = 0.001 causes barely appreciable degradation of image and starts affecting the quality and recognition of image. For this reason, the reference with a Gaussian noise addition of σ = 0.001 is used as comparison with evaluate the recognizability and quality of the reconstruction image as well.

The reconstructed result of the samples is shown in Figure 7, in which (a) and (b) are the result of the reference sample, and (c), (d), and (e) are the result of the multilayer sample. The Figure 7a is the original reference image. Figure 7b is degraded image of Figure 7a by adding Gaussian noise of σ = 0.001, used as contrast. Figure 7c–e are the reconstruction result with conventional LFI-SA method, k-means clustering specular removal algorithm by Tao et al., and AIF method, respectively.

Compared with the conventional LFI-SA method and Tao’s algorithms, the AIF method showed advantages in quality and efficiency of reconstruction of the base layer. In the result of the AIF method, the specular reflection of the surface layer is efficiently filtered, and the scatter and diffuse of the base layer is correctly displayed.

Theoretically, conventional LFI-SA method is capable to resolve depth with any depth resolution when the synthetic aperture size and acquisition spacing is allowed. However, it is infeasible to resolve layers of sub-millimeter spacing with conventional LFI-SA methods in practice, especially for multilayer structures. As shown in Figure 7c, the conventional LFI-SA method failed to resolve the characteristics of the base layer. This may be caused by two reasons. First, as explained in Section 3.1, the reflection of the base layer is submerged by the reflection of the surface layer. Second, in the practical depth estimation algorithms based on defocusing and correspondence, the depth resolution relies on the resolution of the light field in all four $xyuv$ dimensions and the optical resolving power of the acquisition devices, which is limited for terahertz LFI system as well.

Besides the conventional LFI-SA method, Tao’s method showed incompatibility to terahertz light fields as well. Specifically, the remaining of specular reflection from the surface layer is still significant, and the characteristics of the base layer are still submerged and indistinguishable. This may be caused by two reasons. First, terahertz images have no multi-color channel, which reduced the available dimensions of the k-means clustering method and affected the sensitivity of the algorithm. Second, the specular reflection affects the adjacent pixels in all four dimensions of the light field, which the k-means clustering algorithm is unable to handle. In contrast, the result of our method showed a better quality and similarity to the reference image than any other reconstruction methods.

In the case of quantitatively analysis, the performance of the methods and the correctness of the results are measured by the structural similarity (SSIM) [42] between the reconstructed images and the reference image.

The SSIM measurement of the results of ours, Tao’s, conventional LFI-SA methods and the noise-induced reference image are compared for quantitative evaluation of performance as well. As shown in Table 1, in terms of SSIM value, the AIF method is significantly superior to Tao’s method, which cannot open the gap with the conventional LFI-SA method, and is ineffective to separate reflections of different layers in the terahertz waveband.

4.3. Layer Separation and Reconstruction of Practical Multilayer Sample

In this section, the performance and robustness of the AIF method is further validated by applying the AIF method on “practical” marine antifoul coating samples with non-uniform coatings and a rusted base layer.

The “practical” sample is a multilayer marine antifoul coating on a metal basement and rusted layer, as shown in Figure 8. The sample has two layers of coating, and the total thickness of the coatings is about 300 μm. The base layer of the sample is metal and has a rusted layer of 50 μm on average between the coatings and the base layer. The number of sub-images acquired in the experiment is 11 × 31, and the spatial interval is 5 mm. In this experiment, we concentrate on the separation between the base layer and the total coating layers, and the extraction of the rusted features on the base layer.

Figure 9 is the measured BRDF distribution of the coating layers and the rusted base layer, respectively. As shown in the figure, because of the BRDF characteristics of the coatings, the reflection of the base layer is submerged by the reflection of the coating layers in acquired 4D light field as well, such as the idealized sample described in Section 4.2.

Moreover, compared with the silicon plate of the “idealized” sample, the marine coating showed a more significant non-uniformity in BRDF characteristics. Specifically, the BRDF distribution of the marine coating is less concentrated and regular. This is caused by the reflection characteristics and quality of the coatings, and made it more difficult to separate the layers with AIF method. Here, the k-means-clustering-based specular removal algorithm of Tao et al. is used as comparison of the AIF method. Both methods are used for separating and reconstructing the surface and the base layer of the “practical” sample, respectively.

For the reconstruction result of the base layer, considering there are no available “ideal” references for quantitatively evaluation, the reconstruction result of the base layer is not evaluated quantitatively, as in Section 4.2. Instead, because both of the compared methods separate and reconstruct layers from the full 4D light field data, the coating layer reconstruction result is used as a supplement to evaluate the correctness and quality of the base layer reconstruction results.

Figure 10a–c is the reconstruction results of the coating layer with the conventional LFI-SA method, Tao’s algorithm and the AIF method, respectively. For the coating layer, considering the BRDF characteristics of the marine coatings and the base layer, the reconstruction result of conventional LFI-SA method is used as reference. Ideally, the reconstruction result of the coating layer should match the result of the conventional LFI-SA method. As shown in Figure 10b, the result of Tao’s algorithm showed a poor quality in separating and reconstructing the coating layer. Specifically, the reflection from the base layer is enhanced instead of suppressed in the reconstruction result, which made the features of the coating layer indistinguishable. On the other hand, as shown in Figure 10c, in the result of our AIF method, the reflection from the base layer is significantly filtered and suppressed, and the reconstruction result matches the result of the conventional LFI-SA method.

Combined with the analysis in Section 4.2, the result of LFI-SA could be used as a reference image of the coating layer reconstruction with other methods for quantitatively analysis, considering the differences in reflection intensity between the coating layer and the base layer. As shown in Table 2, when using the result of LFI-SA method as reference, the AIF method result showed a significant better performance than that of Tao’s method as well.

Figure 10d,e is the reconstruction result of the base layer with k-means-clustering-based specular removal algorithm and AIF method, respectively. For the base layer, as shown in Figure 10d, the k-means clustering method failed to reveal and reconstruct the corrosion condition of the base layer.

Besides the inherent limitation of the k-means clustering method in terahertz waveband, the non-uniformity in BRDF distribution of the coating layer is another main issue. On the one hand, at the position where the reflection of the coating layer is weak, the reflection of the base layer is over-filtered. On the other hand, at the position where the reflection of the coating layer is strong, the reflection of the base layer is still submerged after the filtering. This limited the performance of the k-means clustering methods in multilayer resolving. Moreover, since the k-means clustering method is unable to obtain a reliable reconstruction result of the coating layer, the reconstruction of the base layer is also considered to not be reliable.

Meanwhile, our AIF method showed a better performance in the reconstruction of the base layer. As shown in Figure 10e, the specular reflection of the coating layers is correctly filtered. At the parts with both strong and weak reflection of the coating layer, the rusted characteristics of the base layer are clearly reconstructed.

As a contrast, Figure 10f showed the imaging result with single camera and the same optical configurations. As shown in the figure, in the majority of the image, features of the base layer are submerged and indistinguishable because of the saturation caused by specular reflection and strong attenuation caused by non-uniformity of the coating layer. Although some of the features on the base layer are still visible, they are significantly attenuated and distorted because of the light source characteristics and limited physical aperture.

5. Conclusions

In this paper, for the problem of fast, layer-resolving imaging on thin multilayer structures, a novel AIF method based on active illumination LFI-SA is proposed. An active illumination terahertz LFI-SA system is established, and both “idealized” and “practical” multilayer samples are used in the experiment for the validation and evaluation of the method. In both of the experiment on “idealized” and “practical” multilayer samples, the AIF method is proven effective on layer-resolving imaging and reconstruction.

Currently, only two-layer separation and reconstruction on layers with significant BRDF distribution divergence and clear interface is initiated. For object with more layers or more complex BRDF characteristics, the AIF method is promising in layer-resolving imaging and interested characteristics extraction as well. For this reason, the AIF method is also promising in the application fields which require fast, see-through and layer-resolving inspection which conventional methods are unable to handle.

Author Contributions

Conceptualization, C.Z., N.L., J.Z. and Y.Z.; methodology, N.L.; software, N.L.; validation, C.Z., N.L., J.Z. and Y.Z.; formal analysis, N.L.; investigation, N.L.; resources, C.Z., J.Z. and Y.Z.; data curation, N.L.; writing—original draft preparation, N.L.; writing—review and editing, C.Z.; visualization, N.L.; supervision, C.Z.; project administration, C.Z.; funding acquisition, C.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Key Scientific Instrument and Equipment Development Projects of China, grant number DH–2012YQ14005; and National Basic Research Program of China (973 Program), grant number 2014CB339806.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. The reflection of terahertz beam on surfaces of different BRDF-characteristics, (a) the specular reflection and (b) the diffuse reflection.

Figure 2. The BRDF-distribution measurement result of surfaces with different reflection characteristics.

Figure 3. The reflection of terahertz beam on “idealized” multilayer structure model.

Figure 4. The framework of the AIF algorithm implementation.

Figure 5. System setup of the reflection-type active illumination terahertz LFI system.

Figure 6. The structure of the “idealized” multilayer sample used in the experiment.

Figure 7. The reconstruction result of the reference sample and the base layer of the “idealized” multilayer sample with different methods. (a) The reference image. (b) The reference image with Gaussian noise of σ = 0.001. (c) The conventional LFI-SA method. (d) The Tao’s method based on k-means clustering. (e) The AIF method.

Figure 8. The structure of the “practical” sample of marine antifoul coating used in the experiment.

Figure 9. The BRDF-distribution measurement result of the coating layers and the rusted base layer of the marine antifoul coating sample.

Figure 10. The reconstruction results of the “practical” marine antifoul coating sample with different methods. (a) The conventional LFI-SA method. (b) Coating layer with the Tao’s method based on k-means clustering. (c) Coating layer with the AIF method. (d) Base layer with the Tao’s method. (e) Base layer with the AIF method. (f) Imaging result with single camera, as contrast.

Table 1

SSIM results of the “idealized” sample with different reconstruction methods.

Methods and Results	SSIM Value
AIF method (Ours)	0.6419
Conventional LFI	0.3629
K-means clustering (Tao et al.)	0.4204
Reference with Gaussian noise, σ = 0.001	0.6668

Table 2

SSIM results of the coating layer with different reconstruction methods.

Methods and Results	SSIM Value
AIF method (Ours)	0.9239
K-means clustering (Tao et al.)	0.5945

References

1. Mittleman, M.D.; Jacobsen, R.H.; Nuss, M.C. T-ray imaging. IEEE J. Sel. Top. Quantum Electron.; 1996; 2, pp. 679-692. [DOI: https://dx.doi.org/10.1109/2944.571768]

2. Mittleman, M.D.; Jacobsen, R.H.; Nuss, M.C. Twenty years of terahertz imaging. Opt. Express; 2018; 26, 9417. [DOI: https://dx.doi.org/10.1364/OE.26.009417] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29715894]

3. Guillet, J.-P.; Recur, B.; Frederique, L.; Bousquet, B.; Canioni, L.; Manek-Hönninger, I.; Desbarats, P.; Mounaix, P. Review of Terahertz Tomography Techniques. J. Infrared Millim. Terahertz Waves; 2014; 35, pp. 382-411. [DOI: https://dx.doi.org/10.1007/s10762-014-0057-0]

4. Zhong, H.; Xu, J.; Xie, X.; Yuan, T.; Reightler, R.; Madaras, E.; Zhang, X.-C. Nondestructive defect identification with terahertz time-of-flight tomography. IEEE Sens. J.; 2005; 5, pp. 203-208. [DOI: https://dx.doi.org/10.1109/JSEN.2004.841341]

5. Tanabe, T.; Watanabe, K.; Oyama, Y.; Seo, K. Polarization sensitive THz absorption spectroscopy for the evaluation of uniaxially deformed ultra-high molecular weight polyethylene. NDT E Int.; 2010; 43, pp. 329-333. [DOI: https://dx.doi.org/10.1016/j.ndteint.2010.03.001]

6. Groves, R.M.; Pradarutti, B.; Kouloumpi, E.; Osten, W.; Notni, G. 2D and 3D non-destructive evaluation of a wooden panel painting using shearography and terahertz imaging. NDT E Int.; 2009; 42, pp. 543-549. [DOI: https://dx.doi.org/10.1016/j.ndteint.2009.04.002]

7. Jackson, J.; Mourou, M.; Whitaker, J.; Duling, I.; Williamson, S.; Menu, M.; Mourou, G. Terahertz imaging for non-destructive evaluation of mural paintings. Opt. Commun.; 2008; 281, pp. 527-532. [DOI: https://dx.doi.org/10.1016/j.optcom.2007.10.049]

8. Caumes, J.P.; Younus, A.; Salort, S.; Chassagne, B. Terahertz tomographic imaging of XVIIIth Dynasty Egyptian sealed pottery. Appl. Opt.; 2011; 50, pp. 3604-3608. [DOI: https://dx.doi.org/10.1364/AO.50.003604]

9. Kapilevich, B.; Pinhasi, Y.; Arusi, R.; Anisimov, M.; Hardon, D.; Litvak, B.; Wool, Y. 330 GHz FMCW Image Sensor for Homeland Security Applications. J. Infrared Millim. Terahertz Waves; 2010; 31, pp. 1370-1381. [DOI: https://dx.doi.org/10.1007/s10762-010-9702-4]

10. Marchese, L.E.; Terroux, M.; Dufour, D.; Bolduc, M.; Chevalier, C.; Genereux, F.; Jerominek, H.; Bergeron, A. Case study of concealed weapons detection at stand-off distances using a compact, large field-of-view THz camera. Micro- and Nanotechnology Sensors, Systems, and Applications VI; SPIE: Baltimore, MD, USA, 2014.

11. Kemp, M.C.; Taday, P.F.; Cole, B.E.; Cluff, J.A.; Fitzgerald, A.J.; Tribe, W.R. Security applications of terahertz technology. Terahertz for Military and Security Applications; SPIE: Orlando, FL, USA, 2003; pp. 44-52.

12. Stoik, C.; Bohn, M.; Blackshire, J. Nondestructive evaluation of aircraft composites using reflective terahertz time domain spectroscopy. NDT E Int.; 2010; 43, pp. 106-115. [DOI: https://dx.doi.org/10.1016/j.ndteint.2009.09.005]

13. Tu, W.; Zhong, S.; Shen, Y.; Incecik, A. Nondestructive testing of marine protective coatings using terahertz waves with stationary wavelet transform. Ocean Eng.; 2016; 111, pp. 582-592. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2015.11.028]

14. Tu, W.; Zhong, S.; Incecik, A.; Fu, X. Defect feature extraction of marine protective coatings by terahertz pulsed imaging. Ocean Eng.; 2018; 155, pp. 382-391. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2018.01.033]

15. Takayanagi, J.; Jinno, H.; Ichino, S.; Suizu, K.; Yamashita, M.; Ouchi, T.; Kasai, S.; Ohtake, H.; Uchida, H.; Nishizawa, N. et al. High-resolution time-of-flight terahertz tomography using a femtosecond fiber laser. Opt. Express; 2009; 17, pp. 7533-7539. [DOI: https://dx.doi.org/10.1364/OE.17.007533]

16. Marchese, L.; Bolduc, M.; Tremblay, B.; Doucet, M.; Oulachgar, H.; Noc, L.L.; Williamson, F.; Alain, C.; Jerominek, H.; Bergeron, A. A microbolometer-based THz imager. Proceedings of the SPIE Orlando; Orlando, FL, USA, 5–9 April 2010; 7671.

17. Dufour, D.; Marchese, L.; Terroux, M.; Oulachgar, H.; Généreux, F.; Doucet, M.; Mercier, L.; Tremblay, B.; Alain, C.; Beaupré, P. et al. Review of terahertz technology development at INO. J. Infrared Millim. Terahertz Waves; 2015; 36, pp. 922-946. [DOI: https://dx.doi.org/10.1007/s10762-015-0181-5]

18. Al Hadi, R.; Sherry, H.; Grzyb, J.; Zhao, Y.; Forster, W.; Keller, H.M.; Cathelin, A.; Kaiser, A.; Pfeiffer, U.R. A 1 k-Pixel Video Camera for 0.7–1.1 Terahertz Imaging Applications in 65-nm CMOS. IEEE J. Solid-State Circuits; 2012; 47, pp. 2999-3012. [DOI: https://dx.doi.org/10.1109/JSSC.2012.2217851]

19. Singh, S.P.; Jha, A.K.; Akhtar, M.J. A contactless thickness measurement of multilayer structure using terahertz time domain spectroscopy. Proceedings of the 2015 IEEE Conference on Antenna Measurements & Applications (CAMA); Chiang Mai, Thailand, 30 November–2 December 2015; pp. 1-4.

20. Krimi, S.; Klier, J.; Jonuscheit, J. Highly accurate thickness measurement of multi-layered automotive paints using terahertz technology. Appl. Phys. Lett.; 2016; 109, pp. 7518-7526. [DOI: https://dx.doi.org/10.1063/1.4955407]

21. Krimi, S.; Torosyan, G.; Beigang, R. Advanced GPU-Based Terahertz Approach for In-Line Multilayer Thickness Measurements. IEEE J. Sel. Top. Quantum Electron.; 2016; 23, pp. 1-12. [DOI: https://dx.doi.org/10.1109/JSTQE.2016.2646520]

22. Van der Valk, N.C.J.; van der Marel, W.A.M.; Planken, P.C.M. Terahertz polarization imaging. Opt. Lett.; 2005; 30, pp. 2802-2804. [DOI: https://dx.doi.org/10.1364/OL.30.002802]

23. Doradla, P.; Alavi, K.; Joseph, C.S.; Giles, R.H. Detection of colon cancer by continuous-wave terahertz polarization imaging technique. J. Biomed. Opt.; 2013; 18, 090504. [DOI: https://dx.doi.org/10.1117/1.JBO.18.9.090504]

24. Ying, L.; Shen, Y.; Zhao, G. Terahertz polarization imaging based on the continuous wave terahertz radiations. Proceedings of the Spie the International Society for Optical Engineering; Beijing, China, 17–19 May 2015.

25. Levoy, M. Light Fields and Computational Imaging. Computer; 2006; 39, pp. 46-55. [DOI: https://dx.doi.org/10.1109/MC.2006.270]

26. Ng, R. Digital Light Field Photography. Ph.D. Thesis; Stanford University: Stanford, CA, USA, 2006.

27. Georgiev, T.; Lumsdaine, A. Depth of Field in Plenoptic Cameras; Eurographics: Munich, Germany, 2009.

28. Adelson, E.H.; Wang, J.Y.A. Single lens stereo with a plenoptic camera. IEEE Trans. Pattern Anal. Mach. Intell.; 1992; 2, pp. 99-106. [DOI: https://dx.doi.org/10.1109/34.121783]

29. Zhang, Y.; Lv, H.; Liu, Y.; Wang, H.; Wang, X.; Huang, Q.; Xiang, X.; Dai, Q. Light-Field Depth Estimation via Epipolar Plane Image Analysis and Locally Linear Embedding. IEEE Trans. Circuits Syst. Video Technol.; 2016; 27, pp. 739-747. [DOI: https://dx.doi.org/10.1109/TCSVT.2016.2555778]

30. Lüke, J.P.; Marichal-Hernández, J.G.; Rosa, F.; Rodríguez-Ramos, J.M. A prototype of a real-time single lens 3D camera. Proceedings of the 3DSA’10; Tokyo, Japan, 19–21 May 2010; pp. 19-106.

31. Bishop, T.E.; Zanetti, S.; Favaro, P. Light field superresolution. Proceedings of the Computational Photography (ICCP), 2009 IEEE International Conference on IEEE; San Francisco, CA, USA, 16–17 April 2009; pp. 1-9.

32. Carles, G.; Downing, J.; Harvey, A.R. Super-resolution imaging using a camera array. Opt. Lett.; 2014; 39, pp. 1889-1892. [DOI: https://dx.doi.org/10.1364/OL.39.001889] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24686631]

33. Wiburn, B. High Performance Imaging Using Arrays of Inexpensive Cameras. Ph.D. Thesis; Stanford University: Stanford, CA, USA, 2004.

34. Coelho, P.A.; Tapia, J.; Pérez, F.; Torres, S.N.; Saavedra, C. Infrared light field imaging system free of fixed-pattern noise. Sci. Rep.; 2017; 7, 13040. [DOI: https://dx.doi.org/10.1038/s41598-017-13595-7]

35. Jain, R.; Grzyb, J.; Pfeiffer, U.R. Terahertz light-field imaging. IEEE Transa. Terahertz Sci. Technol.; 2016; 6, pp. 649-657. [DOI: https://dx.doi.org/10.1109/TTHZ.2016.2584861]

36. Levoy, M.; Ng, R.; Adams, A.; Footer, M.; Horowitz, M. Light Field Microscopy; ACM SIGGRAPH 2006 Papers ACM: New York, NY, USA, 2006; pp. 924-934.

37. Lin, X.; Wu, J.; Zheng, G.; Dai, Q. Camera array based light field microscopy. Biomed. Opt. Express; 2015; 6, pp. 3179-3189. [DOI: https://dx.doi.org/10.1364/BOE.6.003179]

38. Tao, M.W.; Hadap, S.; Malik, J.; Ramamoorthi, R. Depth from Combining Defocus and Correspondence Using Light-Field Cameras. Proceedings of the 2013 IEEE International Conference on Computer Vision IEEE; Sydney, Australia, 3–6 December 2013.

39. Lyu, N.; Zuo, J.; Zhao, Y.; Zhang, C. Terahertz Synthetic Aperture Imaging with a Light Field Imaging System. Electronics; 2020; 9, 830. [DOI: https://dx.doi.org/10.3390/electronics9050830]

40. Tao, M.W.; Su, J.; Wang, T.C.; Malik, J.; Ramamoorthi, R. Depth Estimation and Specular Removal for Glossy Surfaces Using Point and Line Consistency with Light-Field Cameras. IEEE Trans. Pattern Anal. Mach. Intell.; 2016; 38, pp. 1155-1169. [DOI: https://dx.doi.org/10.1109/TPAMI.2015.2477811]

41. Tao, M.W.; Wang, T.C.; Malik, J.; Ramamoorthi, R. Depth Estimation for Glossy Surfaces with Light-Field Cameras. Proceedings of the European Conference on Computer Vision; Zurich, Switzerland, 6–12 September 2014.

42. Horé, A.; Ziou, D. Image quality metrics: PSNR vs. SSIM. Proceedings of the 20th International Conference on Pattern Recognition, ICPR 2010; Istanbul, Turkey, 23–26 August 2010.