1. Introduction

Turbid media, composed of wavelength-scale inhomogeneous particles, can strongly scatter the light, resulting in randomly-distributed rays [1]. Multiple scattering imposes a fundamental limit on the focusing of light due to the optical diffusion [2]. This in turn confines the effective depth of focus to only few multiples of the scattering mean free path [3]. Recent advances in optical imaging, such as wavefront shaping [4] or utilizing the memory effect by angular scanning of the scattered light’s speckle patterns [5], have addressed some of these limitations. However, they require either the measurement of the scattered beam on both sides of the medium [4], or a set of angular-resolved scans of a coherent laser beam [5], thus necessitating long acquisition times and complicated setups. Similarly, optical techniques utilizing only the multiply-scattered photons hold an exciting promise in deep tissue diagnosis [6]. Nonetheless, broadband spectroscopy using multiply-scattered diffuse waves for extraction of intra- or inter-molecular resonance absorption modes has not been demonstrated.

In many crystalline solids and polar gases, the timescales of the weak collective motions and the relaxation phonon modes correspond to the oscillation periods of THz waves [7]. Therefore, THz spectroscopy can provide complementary information to intra-molecular interactions probed by techniques such as infrared spectroscopy [8–12]. However, when a sample is immersed in a turbid medium, although the multiply-scattered waves can still carry the spectral fingerprints, the scattering artifacts distort or obscure the spectroscopic features [13]. The cross-correlations between diffuse THz waves have been proved to reveal the dynamics of the scatterers [14]. Nevertheless, spectroscopic analysis of dielectric resonant modes using diffuse THz waves has not been possible. Different analytical and numerical approaches, such as dense medium theory [15] or effective medium theory incorporating the Mie scattering amplitudes [16], have been proposed to estimate the THz scattering loss in granular materials and improving the resolution in terahertz time-domain spectroscopy (THz-TDS) imaging [17]. The scattering-induced attenuation in ballistic or quasi-ballistic THz beams in the pressed pellets of $\alpha$-lactose monohydrate and $\alpha$-D-glucose has been investigated [18]. A blind demodulation algorithm based on low-rank matrix recovery and alternating minimization has been proposed to remove the sweep distortion originated from densely-layered structures [19]. Signal processing techniques such as wavelet [20] and cepstral analysis [21] have been used to mitigate the scattering artifacts in diffuse THz waves in transmission geometry. Wavelet-based spectral decomposition has been utilized to identify the resonant frequencies of THz reflection coefficients in the presence of severe rough surface scattering [22–24].

Despite many previous attempts, accurate reflection chemical imaging and characterization of materials cloaked with a scattering medium has not been demonstrated. In this work, we present a novel computational technique for creating spatially- and spectrally-resolved chemical maps through a diffusive cloak using THz-TDS. Our approach is based on the bimodality coefficient spectrum, which we define using the skewness and kurtosis of spectral images of a sample. The bimodality coefficient has the unique property of quantifying the modal characteristics of an image histogram distribution. We first decompose the measured reflection spectra using the maximal overlap discrete wavelet transform (MODWT). We reconstruct the spectral images using only a subset of MODWT coefficients such that the bimodality coefficient is maximized at resonant frequencies. To demonstrate the robustness of our technique in reflection imaging of complex samples, we created several spectroscopic Boehler star targets composed of different chemicals with overlapping spectral features. In all the test samples, the scattering effects are designed to be so severe that the raw spectra could not be used to resolve the shape of the Boehler star target or identify the chemicals. However, the experimental results prove that our approach can resolve diffraction-limited spectroscopic images of the sample buried beneath the scattering layer. Importantly, our approach does not rely on hardware modifications to a standard THz spectroscopy setup.

2. Methods

2.1 THz PHASR Scanner

The THz-TDS measurements are obtained using our PHASR (Portable HAndheld Spectral Reflection) Scanner, described in detail previously [25]. In this design, all the optical components shown schematically in Fig. 1(a) are placed inside a 3D-printed fiber-coupled housing (dimensions: $37.3\times 14\times 25.1$ cm$^3$), shown in Fig. 1(b). In the PHASR Scanner, the TERA ASOPS (Asynchronous OPtical Sampling) dual-fiber-laser THz spectrometer (Menlo Systems, Inc., Newton, NJ, USA) is incorporated into a handheld, collocated, telecentric imaging system. A THz beam generated by the photoconductive antenna (PCA) in the emitter (E) is collimated using a TPX lens (L1) with 50 mm focal length. The collimated beam is directed towards a gimballed mirror (GM) using a high-resistivity silicon beam splitter (BS). The mirror gimbal is mounted on a two-axis motorized system composed of a goniometer and a rotational stage. It raster scans the collimated beam over the aperture of a custom-made telecentric f-$\theta$ lens [26]. Therefore, the focused beam is always in parallel to the optic axis and has a constant focal spot-size of approximately 3.1$\times$2.4 mm$^2$. The reflected beam retraces the path of the incident beam back to the BS, where it is directed towards the second collimating lens (L2) and is focused on the PCA inside the detector (D). Previously, we have conducted multiple experiments on the utility of the PHASR Scanner for biomedical sensing applications [27–29].

 

Fig. 1. (a),(b) Schematic and photograph of the optical components inside the PHASR Scanner housing. (c) Boehler star covered by a diffusive scattering cloak with a fixed thickness of about 4.63 mm. Example incident (black) and reflected (blue) THz-TDS waveforms through the scattering layer are shown. (d) Image of the target formed using the extinction coefficients, $\varepsilon (f)$, at 0.53 THz without a scattering layer in the beam path (the ground truth). The radius of the red circle in the center is 2.1 mm, corresponding to an 1.11-mm spatial resolution. (e) Wavelength-normalized transport mean free path, $l_{\text {tr}}$, of the scattering layers of 0.3 and 0.5 g/cm$^3$ volume density. The black and magenta lines show the wavelength-normalized propagation distance of one-way (toward the target) and round-trip (back to the detector) travels through the scattering layer, respectively.

Download Full Size | PDF

2.2 Imaging target

To assess the imaging performance in spatially resolving the THz resonant modes, the sample chemicals are placed inside the petals of a Boehler star resolution target. Figure 1(c) shows a 3D-printed Boehler star whose petals are filled with $\alpha$-lactose monohydrate with particle size smaller than 10 $\mu$m. The dashed magenta square in Fig. 1(c) shows a 12$\times$12-mm$^2$ field-of-view. Figure 1(d) illustrates the image of the Boehler star without the presence of any diffusive medium in the beam path. The color axis in Fig. 1(d) is given by the extinction coefficients, $\varepsilon (f)$, at 0.53 THz. A THz-TDS measurement is obtained at each 0.25$\times$0.25-mm$^2$ pixel area, resulting in a 48$\times$48-pixel image of the region enclosed by the magenta square in Fig. 1(c). The petals of $\alpha$-lactose monohydrate are clearly discernible in Fig. 1(d). Moreover, the radius of the unresolved circle in the center corresponds to the diffraction-limited imaging resolution at 0.53 THz (focal length = 40 mm, numerical aperture (NA) = 0.40). Figure 1(d) serves as the ground truth for comparison to imaging in the presence of diffusive scattering.

2.3 Turbid scattering medium

The imaging target is cloaked by a scattering layer as shown in Fig. 1(c). The scattering layer is made from low-density polyethylene (LDPE) powders of 180-$\mu$m mean particle size. In order to create a scattering layer with known spatial dimensions, we created a 3D-printed housing to be placed on the top of the sample and filled with LDPE. The LDPE powders result in significant Mie scattering over the measurement bandwidth of $f = 0.2\text {--}1.1$ THz. LDPE particles are chosen because they are nonabsorbent in the THz region, and their refractive index is independent of the wavelength. Therefore, the pulse attenuation can only be associated with the scattering. To investigate the robustness of our technique, we study two different cloaks with 0.3 and 0.5 $\text {g}/\text {cm}^3$ particle volume density. The thickness of the cloaks, i.e., the beam propagation distance in the turbid medium, is kept the same. In electromagnetic scattering theory, traveling a distance longer than the Boltzmann transport mean free path, $l_{\text {tr}}$, randomizes the beam into ballistic and diffusive fields [30]. Given that the density of LDPE powder is approximately equal to 0.91, we can determine the volume density of the scattering layer by measuring the weight of the particles beforehand. The scattering mean free path, $l_\text {s}$, modifies a coherent propagation according to [31],

Figure 2(a) shows the image of the Boehler star in the presence of the cloak with 0.3 $\text {g}/\text {cm}^3$ volume density. The color axis represents the $\varepsilon (f)$ at 0.53 THz. It is evident that the petals of $\alpha$-lactose are not resolved in Fig. 2(a). Here, because a measured $\varepsilon$ is composed of the absorption coefficient of the sample, the scattering-induced attenuation, and the Mie scattering-induced spectral artifacts, the petals of $\alpha$-lactose are obscured. Spatial averaging over disjointed measurements is considered the standard approach to mitigate the scattering effects in THz spectroscopy [21,33]. However, it cannot be used in imaging applications and for characterizing samples made from multiple chemicals when the area covered by each material is not known in advance. The multiresolution spectral analysis approach presented here can mitigate the scattering effects without utilizing spatial averaging. This method relies on decomposition of the measured extinction coefficients using the wavelet transform and reconstruction from wavelet coefficients at specifically-chosen frequency scales. Moreover, unlike the cepstral filtering technique for mitigating the THz scattering artifacts, which requires a priori knowledge on the characteristics of the scattering particles to design a low- or band-pass filter [21,34], our technique is robust in that its implementation is independent of the scattering medium or the properties of the characteristic spectral resonance features of the sample materials.

 

Fig. 2. (a) Image of the $\alpha$-lactose Boehler star covered by a scattering layer having a 0.3 $\text {g}/\text {cm}^3$ volume density. The color axis is given by the measured $\varepsilon (f)$ at 0.53 THz. (b) Image of the same sample formed using $\sum _{j=3}^4D_j$(0.53 THz) as the color map. (c),(d) Similar to (a),(b) for the scattering layer with 0.5 $\text {g}/\text {cm}^3$ volume density. (e),(f) Normalized histograms of the images in (a) and (b), yielding the bimodality coefficient values of $\beta =0.46 < 5/9$ and $\beta =0.56 > 5/9$, respectively. A probability distribution function (red line) is fitted to each histogram.

Download Full Size | PDF

2.4 Wavelet multiresolution analysis of extinction coefficients

The discrete wavelet transform (DWT) decomposes a signal into sets of wavelet and scaling coefficients. Each set describes the localized changes in the signal at a specific variation scale. The scale is an interval of signal over which localized, weighted averages of the signal (i.e., the scaling coefficients) and the differences of those averages (i.e., the wavelet coefficients) are calculated. Application of DWT directly to the extinction coefficients results in wavelet scales corresponding to dyadic increments of the frequency sampling interval given by $\sigma _j=2^{j-1}\delta f$ [35,36]. Here, $\sigma _j$ is the wavelet scale of the $j^{\text {th}}$ level of DWT. The $\delta f$ is the frequency-domain sampling interval. Its value depends on the sampling rate of measured THz pulses, $f_s$, in addition to their size, $N$, and is given by $\delta f=f_s/N$. For example, given a $\delta f\approx 0.02$ THz, the 4$^\text {th}$ and 5$^\text {th}$-level wavelet coefficients correspond to $\sigma _4=2^3\times 0.02=0.16$ THz and $\sigma _5=2^4\times 0.02=0.32$ THz intervals in the extinction coefficients, respectively. The wavelet multiresolution analysis (MRA) of an extinction coefficient measured at a pixel with spatial coordinates of $(x,y)$, $\varepsilon _{(x,y)}$, using the maximal overlap DWT (MODWT) at $J$ levels of decomposition is given by [35,36],

Figure 2(b) shows the image of the Boehler star formed using MODWT MRA in the presence of the cloak with 0.3 g/cm$^3$ volume density. The color axis in Fig. 2(b) is $\sum _{j=3}^4\tilde {D}_j$(0.53 THz), which represents the reconstruction of $\varepsilon$ from only the third- and fourth-level wavelet coefficients. The coefficients of $\sum _{j=3}^4\tilde {D}_j(f)$ are associated with the cepstral content of $\varepsilon$ in the quefrency range of $[1/(2^{5}\delta f):1/(2^3\delta f)]$. It can be seen that the petals of $\alpha$-lactose are clearly resolved in Fig. 2(b). In addition, the spatial resolution is comparable to the image of the sample obtained without a cloak in the beam path, shown in Fig. 1(d). In particular, the diameter of the red circle in the center of Fig. 2(b) is found to be 4.2 mm, corresponding to a 1.1-mm spatial resolution using $\frac {\text {4.2}\times \pi }{2*\text {N}_{\text {petals}}}=\text {1.1}$, where $\text {N}_{\text {petals}}$ is the number of petals in the target [38]. This resolution is larger than the Rayleigh diffraction limit at 0.53 TH because, due to the use of a solid Boehler star mold, the sample powders are not densely-packed under a hydraulic press as is customary when making pellets. Instead, powders of $\alpha$-lactose are loosely-packed in each petal, which results in additional light scattering by the extra air voids between the particles. Figure 2(c) shows the image of the Boehler star obtained in the presence of the cloak with 0.5 g/cm$^3$ volume density. Similar to Fig. 2(a), the color axis in Fig. 2(c) represents $\varepsilon$(0.53 THz). Figure 2(d) illustrates the result of reconstruction of $\varepsilon$ from $\sum _{j=3}^4\tilde {D}_j$. It can be noticed that the petals of $\alpha$-lactose are resolved in the image formed using $\sum _{j=3}^4\tilde {D}_j$(0.53 THz). In this paper, MODWT is calculated using the least asymmetric mother wavelet filter with four vanishing moments, i.e., sym4 or LA(8) mother wavelet [36,39], at $J=6$ levels of decomposition.

2.5 Bimodality coefficient spectrum

In broadband THz spectral imaging of a mixture sample, images formed at resonance absorption frequencies provide a higher visual contrast [10]. We employ this basic concept to find the MODWT details coefficients that can separate the resonant spectral features from the scattering artifacts. In the examples of Figs. 2(a)-(d), we use the higher-order image statistics to determine a specific combination of MODWT details vectors. This combination is utilized in the MODWT MRA given by Eq. (2) [20]. In particular, we use the skewness and kurtosis to define a new spectroscopic concept named the bimodality coefficient spectrum. Bimodality coefficient can measure the modal characteristics of an image histogram distribution. Skewness (a measure of histogram asymmetry) and kurtosis (a measure of histogram tailedness) represent the third- and fourth-order standardized moments of an image around its mean, respectively [40]. The skewness of an image with $n$ pixels is given by [40],

To find a unique combination of MODWT details coefficients, we calculate the $\beta$ of the images formed at a resonant frequency using all possible combinations of $\tilde {D}_j$ in the MRA reconstruction of Eq. (2). For example, $\sum _{j=2}^4\tilde {D}_j$ and $\sum _{j=3}^4\tilde {D}_j$ are two forms of such combinations. At $J=6$ levels of decomposition, there are $2^{J+1}-1=127$ unique combinations. Subsequently, we select the combination that yields the image with highest $\beta$ value to reconstruct the $\varepsilon$.

2.6 Principal component analysis

In hyperspectral imaging, the measured spectra are composed of contributions from disparate substances [44]. Therefore, several spectral unmixing techniques have been proposed to distinguish the constituent materials and their proportions in each spectrum measurement [45]. Principal component analysis (PCA) is an unsupervised spectral unmixing approach and a mathematically optimal way to treat similar problems [45]. We can also treat the problem of extracting the resonance absorption signatures from the $\varepsilon$ in the presence of scattering as a spectral unmixing case.

PCA identifies a set of orthogonal axes, onto which the projection of the data conveys the highest variance [46]. These orthogonal axes are known as principal components (PC). The PCs along a desired dimension, here the frequency axis, can be obtained using the eigen-decomposition of the sample covariance matrix given by [44,46],

It should be noted that the characteristic resonances are not necessarily separable using the PCA approach, whereas the MODWT MRA retains the spectral information. Thus, MODWT MRA enables specificity by isolating the petal associated with each resonance into wavelet coefficients at that resonant frequency. The flowchart of Fig. 3 compares the algorithms of MODWT MRA and PCA.

 

Fig. 3. the computational steps of the proposed MODWT MRA technique and the PCA algorithm. The first row shows that a three-dimensional $\varepsilon _{(x,y)}(f)$ is calculated from the measured time-domain THz pulses. The second row shows that the implementation of PCA on $\varepsilon$ fails to preserve the frequency-dependence of $\varepsilon$. The third and fourth rows demonstrate that the frequency information is fully preserved at all stages of the bimodality maximization-based MODWT MRA.

Download Full Size | PDF

This flowchart emphasizes the unique property of retention of the frequency information by MODWT MRA. The first row shows that a spatio-spectral data cube composed of $\varepsilon _{(x,y)}$ is obtained by Fourier transformation and deconvolution of the measured THz-TDS pulses reflected from the sample. Deconvolution removes the response of the measurement setup to extract the sample response function. The second row illustrates the implementation of the PCA on $\varepsilon$. PCA results in a set of two-dimensional PC$_{(x,y)}$ which do not preserve the frequency-dependent information. In contrast, the third row shows that the frequency dimension of $\varepsilon$ is fully preserved at all stages of the bimodality-based MODWT MRA algorithm. In other words, the MODWT vectors, $\tilde {D}_1$, $\tilde {D}_2$,…, $\tilde {D}_J$ and $\tilde {S}_J$, retain the spectral information of the $\varepsilon$ in their 3D structure. This property, in turn, enables the identification of specific petals at their resonant frequencies using the bimodality coefficient optimization process.

3. Results

3.1 Boehler star with multiple spectrally-overlapping samples

In this section, we investigate the robustness of the bimodality-based MODWT MRA of THz reflection images in the presence of scattering. Figure 4 shows the obtained experimental results when each of the five petals of a Boehler star is filled by a different chemical powder. Fig. 4(a) shows the sample target made from $\alpha$-lactose monohydrate, 4-aminobenzoic acid (PABA), 3-hydroxybenzoic acid (3HBA), riboflavin, and high-density polyethylene (HDPE). Figure 4(b) compares the $\varepsilon$ of these five chemicals over the frequency range of $f = 0.2 - 1.1$ THz without a scattering layer in the beam path. The distinct characteristic resonant signatures, including $\alpha$-lactose resonance at 0.53 THz, PABA resonances at 0.6 and 0.8 THz, 3HBA resonance at 0.8 THz, and riboflavin weak resonance at 1.03 THz, are observable in Fig. 4(b). These resonant frequencies are in agreement with the values reported in the literature [47–50]. The extinction coefficients of HDPE (cyan line) and the Boehler star’s resin-based building material (red line) do not contain any resonance absorption lines. In addition to the scattering effects described earlier, the overlap between the resonances of PABA and 3HBA with the broadband absorption of the Boehler star resin itself imposes a further computational challenge in resolving the petal associated with each of these chemicals. Furthermore, as the signal-to-noise ratio (SNR) decreases in the higher frequencies, the $\varepsilon$ of riboflavin, 3HBA, PABA, $\alpha$-lactose, and the resin background overlap, thus further obscuring the image of the riboflavin resonance at 1.03 THz. In this example, maximizing the $\beta (f)$ at 0.53 and 0.6 THz over all unique combinations of $\tilde {D}_j$ reveals that $\sum _{j=3}^6\tilde {D}_j$ yields the highest bimodality value at both resonant frequencies. Figure 4(c) shows the $\beta$ of the images formed using the $\varepsilon$ (red line) in the presence of cloak of 0.3 $\text {g}/\text {cm}^3$ volume density, and their reconstruction from $\sum _{j=3}^6\tilde {D}_j$ (blue line). It can be noticed that in the $\beta$ of $\sum _{j=3}^6\tilde {D}_j$, both resonances are identified as the only local maxima with a bimodality value greater than 5/9. The similarity between the shapes of the resonant signatures at 0.53 and 0.6 THz can explain why the same set of $\tilde {D}_j$ is able to extract them at their respective frequencies. The images shown in Figs. 4(d)-(h) are obtained from a 12$\times$12-mm$^2$ area (48$\times$48-pixels at 0.25-mm pixel size) at the center of the Boehler star. Figure 4(d) shows an image formed by $\sum _{j=3}^6\tilde {D}_j(f)$ at 0.53 THz with a $\beta =0.77$. It can be seen that the petal of $\alpha$-lactose has a higher contrast in comparison to rest of the sample. Figure 4(e) shows an image formed by the same combination at 0.6 THz, yielding a $\beta =0.75$. This image provides a high contrast by identifying only the petal of PABA. Figures 4(f)-(h) show the images formed by MODWT details vectors, including $\tilde {D}_4$(0.8 THz), $\tilde {D}_5$(0.8 THz), and $\tilde {D}_6$(1.03 THz), respectively. Figure 4(f) demonstrates that $\tilde {D}_4(f)$ can extract and isolate the broader resonant feature of PABA at 0.8 THz, resulting in a higher contrast over the petal of PABA, which is comparable to the petal resolved by $\sum _{j=3}^6\tilde {D}_j(f)$ at 0.6 THz. In contrast, Fig. 4(e) exhibits that $\tilde {D}_5(f)$ at 0.8 THz yields a higher contrast on the petal of 3HBA. In other words, $\tilde {D}_j$ of different levels can be employed to distinguish between materials of overlapping resonant signatures at the same resonant frequency. It should be noted that because of the shared similarity in shape between the $\varepsilon$ of 3HBA and the Boehler star resin material, as shown in Fig. 4(b), the difference between the 3HBA petal and the background material is modest in Fig. 4(g). Finally, Fig. 4(h) shows an image formed by $\tilde {D}_6(f)$ at 1.03 THz, resolving the petal of riboflavin.

 

Fig. 4. (a) Boehler star target where each petal is filled by a different chemical powder, including $\alpha$-lactose, HDPE, 3HBA, riboflavin, and PABA. The Boehler star is covered with a diffusive cloak with 0.3 $\text {g}/\text {cm}^3$ volume density. (b) $\varepsilon$ of the chemicals measured without a scattering layer in the beam path. The black arrows mark the resonant frequencies. (c) $\beta$ of the images formed using $\varepsilon$ (orange line) and the images formed using $\sum _{j=3}^6\tilde {D}_j$ (blue line). (d)-(h) Images formed using different combinations of $\tilde {D}_j(f)$ at the resonant frequencies, including $\sum _{j=3}^6\tilde {D}_j$(0.53 THz) in (d), $\sum _{j=3}^6\tilde {D}_j$(0.6 THz) in (e), $\tilde {D}_4$(0.8 THz) in (f), $\tilde {D}_5$(0.8 THz) in (g), and $\tilde {D}_6$(1.03 THz) in (h).

Download Full Size | PDF

The experimental results given by Figs. 4(d)-(h) show that different resonance absorption lines are extracted by MODWT details coefficients of different scales. For example, the image formed using $\tilde {D}_4$(0.8 THz) resolves the petal of PABA, whereas the image created by $\tilde {D}_6$(1.03 THz) captures the area of riboflavin. This observation can be explained by the difference between the spectral shapes of different resonant signatures. Because the resonance of PABA at 0.8 THz is both sharper and narrower in comparison to the resonance of riboflavin at 1.03 THz, it is extracted by lower-level MODWT details vectors, which correspond to a higher-quefrency content in $\varepsilon$. On the contrary, the broader resonant signature of riboflavin has a lower-quefrency content. Therefore, it is extracted by MODWT details vector of a higher level.

3.2 Comparison with alternative hyperspectral technique

Figure 5 shows the results of spectral unmixing using the PCA of $\varepsilon$. Fig. 5(a) shows that each PC is formed by a linear, weighted combination of all spectral images along the measurement bandwidth. The weight associated with each frequency in constructing the PCs is determined by the elements of the eigen-vectors of $\text {cov}(E)$. Figure 5(b) gives the eigen-vectors of the first four PCs in the order of capturing a higher variance in the original set. It can be seen that there is no eminent resonant feature in the eigen-vectors of the first (red line) and second (blue line) PCs. However, the eigen-vector of the third PC (green line) has two distinct peaks at 0.53 and 0.6 THz, revealing the resonant frequencies of $\alpha$-lactose and PABA. Similarly, in the fourth eigen-vector (purple line) the image at 0.53 THz has the greatest value. Therefore, the fourth PC is mainly associated with the resonance of $\alpha$-lactose.

 

Fig. 5. (a) Forming a PC using a linear, weighted combination of the spectral images along the measurement bandwidth. The weight of each image is determined by the eigen-vectors of $\text {cov}(E)$. (b) Eigen-vectors of the first four PC bands in the order of capturing a higher variance. (c)-(f) Images of the first four PCs. The third and fourth PCs yield the petals of $\alpha$-lactose and PABA.

Download Full Size | PDF

In comparison to the results obtained using MODWT MRA, our new approach can identify all constituent chemicals in the presence of scattering, albeit with varying levels of success. In contrast, only the third and fourth PCs are able to identify $\alpha$-lactose and PABA, with loss of specificity between the two. In other words, because each PC is formed by a weighted superposition of all spectral images at different frequencies, it contains contributions from all frequency components. Therefore, the characteristic resonances are not necessarily separable using the PCA approach. In contrast, the wavelet MRA retains the frequency information. Thus, it enables specificity by isolating the petal associated with each resonance into wavelet coefficients at that resonant frequency. This comparison reveals that the scale-frequency decomposition using wavelet transform is advantageous in studying the characteristic resonant frequencies of mixture samples.

3.3 Performance evaluation

In this section, we describe a quantitative approach for comparison between the images of the mixture Boehler star formed using $\varepsilon$, $\tilde {D}_j$, and the PCs. We use the image of the Boehler star obtained in the absence of any cloaks in the beam path, shown in Fig. 6(a), as the reference for delineating the boundary pixels of each petal. The color axis in Fig. 6(a) represents the peak-to-peak amplitude of measured THz signals. Figure 6(b) labels the pixels of chemical using a different color and excludes the diffraction-limited resolution circle at the center. The Boehler star background material is set to black, $\alpha$-lactose to purple, HDPE to blue, 3HBA to green, riboflavin to yellow, and PABA to red. The colored region associated with each material forms the label matrix, as the ground truth, for calculating the mean squared error (MSE) at the resonant frequencies. For example, in the label matrix of $\alpha$-lactose, all the pixels in the purple region of Fig. 6(b) are equal to one, while every pixel elsewhere equals zero. In order to compare the images formed by $\varepsilon$, $\tilde {D}_j$, and the PCs to their corresponding label matrix, we first convert them to a set of binary images. In a binary image, a pixel value is set to one if it is assigned to a given material and is set to zero otherwise. We use an image segmentation approach using a variable threshold over the color axis to form the binary images. The variable threshold is used to create a receiver operating characteristic (ROC) curve for each image [51]. In particular, a new binary image is formed at each cut-off threshold value and is compared to the label matrix for calculating the true positive rate (TPR) and false positive rate (FPR). The ROC curve shows the calculation of TPR or sensitivity on the y-axis versus the FPR or (1-specificity) on the x-axis. A high sensitivity indicates a higher rate of correct assignment of pixels belonging to a material at a specific resonant frequency. A high specificity (1-FPR), on the other hand, implies a higher rate of correct assignment of pixels belonging to the rest of the materials. In each ROC plot in Figs. 6(c)-(g), the diagonal dashed lines show the ROC curve of a random classifier with no predictive power. A higher area under the ROC curve (ROC-AUC) indicates a greater accuracy in differentiation between the pixels associated with a specific resonant mode and the rest of the materials.

 

Fig. 6. (a) Image of the mixture sample obtained without the presence of a cloak in the beam path. The image in (a) is formed using the peak-to-peak amplitude of THz signals measured over a 12$\times$12 mm$^2$ area at the center of the Boehler star. (b) Superposition of the sample label matrices can delineate the pixels associated with each chemical using a different color. (c) ROC curves formed using the images of $\varepsilon$(0.53 THz), $\sum _{j=3}^{6}\tilde {D}_j$(0.53 THz), and PC 4. (d) ROC curves formed using the images of $\varepsilon$(0.6 THz), $\sum _{j=3}^{6}\tilde {D}_j$(0.6 THz), and PC 3. (e)-(g) ROC curves of the images formed using $\varepsilon$ and $\tilde {D}_j$ of the materials detected in Figs. 2(d)-(f), respectively. The diamond markers identify the threshold values maximizing the harmonic mean of sensitivity and specificity. (h) Final aggregate chemical map formed by a weighted combination of the thresholded binary images of MODWT details vectors at the resonant frequencies.

Download Full Size | PDF

Figures 6(c)-(g) show the ROC curves of the images at 0.53 THz ($\alpha$-lactose), 0.6 THz (PABA), 0.8 THz (PABA), 0.8 THz (3HBA), and 1.03 THz (riboflavin), respectively. The solid blue, red, and black lines show the ROC curves obtained using the images formed by $\varepsilon$, $\tilde {D}_j$, and PCA, respectively. It should be noted that in the images of PCs, only the petals of $\alpha$-lactose and PABA are resolved. Therefore, we only include the ROC curves of PC 3 and PC 4 in Figs. 6(c),(d) to compare them with the images obtained by $\varepsilon (f)$ and $\tilde {D}_j(f)$ at 0.53 and 0.6 THz. It can be seen that the ROC-AUC value is always higher for images formed using $\tilde {D}_j$ in comparison to the other approaches. Table 1 gives the MSE values at each resonant frequency for the images formed using $\varepsilon$ and $\tilde {D}_j$. The MSE of the MODWT details vectors is always lower compared to the extinction coefficients, corresponding to an error reduction of 82$\%$ at 0.53 THz, 12$\%$ at 0.6 THz, 71$\%$ at 0.8 THz (PABA), 56$\%$ at 0.8 THz (3HBA), and 71 $\%$ at 1.03 THz. Because of the inherent trade-off between sensitivity and specificity, to generate a final aggregate chemical map, the optimum threshold values are selected to maximize the harmonic mean of sensitivity and specificity. The resultant cutoff threshold for each chemical are shown using the diamond markers in Figs. 6(c)-(g), In Fig. 6(h), we show the aggregate chemical map formed by a weighted combination of the thresholded binary images formed by $\tilde {D}_j$. These images include $\sum _{j=3}^6\tilde {D}_j$(0.53 THz), $\sum _{j=3}^6\tilde {D}_j$(0.6 THz), $\tilde {D}_5$(0.8 THz), and $\tilde {D}_6$(1.03 THz). Comparing Fig. 6(h) to the assigned labels in Fig. 6(b) reveals that the chemical map resultant of the proposed method provides a higher spatial resolution towards the center of the Boehler star. In other words, areas not previously identifiable using the peak-to-peak amplitude or the raw $\varepsilon$ are resolved accurately using our proposed approach. One key limitation of this technique is that materials such as HDPE and the resin background used in the Boehler star structure cannot be differentiated because their dielectric functions do not possess any signature resonant frequencies that can be used to encode the reflected spectra. Future directions of this work include employing time-resolved, polarization-sensitive THz-TDP measurements, which can be used to characterize the scattering events [52]. Furthermore, the use of broadband THz sources, such as air-plasma filaments [53], can enable the implementation of the proposed bimodality coefficient technique at higher frequencies for the identification of even more complex samples.

Table 1. MSE between THz images and the ground truth

View Table

4. Conclusion

We presented experimental results to demonstrate THz spectroscopic reflection imaging through a highly-diffusive cloak. We showed that the scattering effects of a diffusive cloak can obscure the characteristic resonance absorption lines of interest. We presented a multiresolution spectral analysis technique using the wavelet decomposition and reconstruction of THz extinction spectra to mitigate those artifacts. We identified the specific wavelet scales that can extract the obscured resonances of several test materials in a complex sample using the higher-order statistics of the spectral images. In particular, we used the skewness and kurtosis of the spectral images to define a new spectroscopic concept named the bimodality coefficient spectrum. We showed that by using a unique combination of the wavelet details vectors in reconstructing the extinction coefficients, the resonant frequencies can be identified as the local maxima above a defined 5/9 threshold in the bimodality coefficient spectrum. Additionally, we showed that overlapping signature resonances can be separated into wavelet coefficients of different scales. We evaluated the performance of our technique in the identification and labeling the pixels associated with each chemical. We showed that our approach outperforms the raw extinction spectra and principal component analysis. Finally, we demonstrated that our proposed technique resulted in a chemical map with a higher spatial resolution, despite the presence of a scattering layer, as compared to the images formed using raw THz spectra in the absence of a diffusive cloak.