## Abstract

Inspection and imaging through scattering layers is a decades-old problem in optics. Unlike x-ray and ultrasonic imaging techniques, time-domain spectroscopy methods can provide detailed chemical and structural information of subsurfaces along with their depth. Although high-resolution time-of-flight measurement in time-domain spectroscopy provides 3D information, unfortunately it also induces an unwanted sensitivity to misalignments of the system and distortion of the layers themselves. Such high sensitivity to alignment and sample surface is a well known problem in time-domain and interferometric imaging, and is a major concern when the alignment error is comparable to the pulse wavelength. Here, we propose and implement an algorithmic framework based on low-rank matrix recovery and alternating minimization to remove such unwanted distortions from time-domain images. The method allows for recovery of the original sample texture in spite of the presence of temporal-spatial distortions. We address a blind-demodulation problem where, based on several observations of the sample texture modulated by undesired sweep distortions, the two classes of signals are separated with minimal damage to the main features. The performance of the method is examined in both synthetic and real data in the case of a terahertz time-domain system, and the successful reconstructions are demonstrated. The proposed general scheme can be implemented to advance inspection and imaging applications in THz and other time-resolved spectral imaging modalities.

© 2016 Optical Society of America

## 1. INTRODUCTION

Due to its fine time resolution and broad spectral coverage, terahertz time-domain spectroscopy (THz-TDS) has become a leading method in THz spectroscopy [1,2], imaging [3], and nondestructive testing [4] of dielectric structures. There is extensive literature on improving THz imaging capability. In general, a large body of literature is focused on improving the hardware, as THz-TDS power levels are usually at sub-microwatt levels [5–7]. This perspective can be extended to improving THz acquisition methodologies such as compressive [8] and wide-field acquisitions [9–12]. On the signal-processing side, the mainstream research has been focused on improvement of signal-to-noise ratio (SNR) [13,14]. Many of these SNR improvement methods have been developed for transmission mode spectroscopy and have appreciable alignment with infrared spectroscopy [15–17]. However, when it comes to imaging, inspection, and content extraction of surfaces and layered structures, THz-TDS can suffer significantly from phase distortions along the sample surface. These sweeping distortions in THz time-domain imaging appear to be the dominant challenge for in-depth imaging and content extraction in densely layered structures, irregular surfaces, or any 2D geometry. Unfortunately, these heavy distortions are directly induced by the nature of the sample and THz-TDS measurement scheme, and cannot be addressed by minor hardware improvement, conventional denoising, or SNR enhancement techniques.

In this paper, we propose and demonstrate a mathematical framework that enables the demodulation of the recorded signal from sweep distortions induced by slight depth variations or the layered structure inter-reflections. We view the problem as a demodulation of the distortion profiles from the sample texture, based on the reflected wave measured at different instances of time. The problem of interest is modeled as a bilinear inverse problem (BIP), which can be addressed using a low-rank matrix-recovery framework. We propose an alternating minimization scheme, where a prior structure is considered for each factor in the BIP. More specifically, we assume the distortion profiles belong to a low-dimensional subspace (which can be extracted from the raw data) and the layer structure is of a binary nature. The model considered for the layer texture is in fact a shape-based approximation (see [18–20] for examples), a phase corresponding to the main texture, and another phase representing the anomalies and inclusions. Using the given priors for each factor in the BIP, we alternatively solve the demodulation problem to exclude the undesired sweep distortions from the THz images.

Our mathematical presentation mainly relies on multidimensional calculus. We use bold characters to denote vectors and matrices. Considering a matrix $\mathit{A}$ and the index sets ${\mathrm{\Gamma}}_{1}$ and ${\mathrm{\Gamma}}_{2}$, we use ${\mathit{A}}_{{\mathrm{\Gamma}}_{1},\text{:}}$ to denote the matrix obtained by restricting the rows of $\mathit{A}$ to ${\mathrm{\Gamma}}_{1}$. Similarly, ${\mathit{A}}_{\text{:},{\mathrm{\Gamma}}_{2}}$ denotes the restriction of $\mathit{A}$ to the columns specified by ${\mathrm{\Gamma}}_{2}$, and ${\mathit{A}}_{{\mathrm{\Gamma}}_{1},{\mathrm{\Gamma}}_{2}}$ is the submatrix with the rows and columns restricted to ${\mathrm{\Gamma}}_{1}$ and ${\mathrm{\Gamma}}_{2}$, respectively.

The remainder of this paper is organized as follows. In Section 2, we overview the physics of the problem and discuss the main demodulation problem to be addressed. In Section 3, we present the main algorithmic framework to address the demodulation problem using an alternating minimization scheme. Finally, in Section 4, the sensitivity and performance of the algorithm are assessed with synthetic data. We further experimentally demonstrate the blind demodulation of the data recorded from both single and multilayered structures and report the algorithm outcomes for the experimental data. We conclude this section with some remarks and future directions of research.

## 2. PHYSICS OF THE PROBLEM

Consider an electric field that is linearly polarized in the $x$ direction, propagating along the $z$ direction in the free space. The waveform is considered to be a finite-duration THz pulse $\chi (t)$, for which the traveling field along every point $(x,y)$ is

and $c$ is the wave speed. Due to the confocal nature of the measurement at each $(x,y)$ position, the problem of interest can be considered as extracting the contents of a dielectric slab placed perpendicular to the wave-propagation axis, based on analyzing the reflected field ${\mathrm{E\u20d7}}_{0}^{-}(z,t)$ [Figs. 1(a) and 1(b)].For a homogeneous layer of width $d$, with reflection coefficient $\rho $ and refraction index ${n}_{\rho}$, the returned signal can be analytically expressed by convolving the pulse $\chi (\xb7)$ with a train of impulse functions with decaying coefficients. More specifically [21],

whereThe response in Eq. (1) can still be reasonably accurate when the dielectric slab is homogeneous along the $z$ axis, i.e., $\rho (x,y,z)={\rho}_{0}(x,y)$ for $0\le z\le d$. It only requires plugging the point-wise values of $\rho $ and ${\tau}_{\rho}$ in the formulation.

In the case of perfect or close to perfect reflection, where $\rho \approx 1$, the first impulse term in Eq. (1) dominates all the other terms and the returned wave is a simple modulation of $\rho $ with the pulse waveform. When $0<\rho \ll 1$, the most dominant terms in $u(\tau )$ correspond to the first and second impulses and the remaining terms decay exponentially fast. Specifically, in many practical scenarios where the pulse width is sufficiently small and the slab width $d$ is large, the first reflected waveforms do not overlap with the subsequent reflections. In this case, still a modulation of $\rho $ with the pulse waveform is observed in different time intervals.

For multilayer dielectric slabs, the impulse response can still be cast as an impulse train; however, deriving closed-form expressions for the coefficients is a sophisticated task and beyond the scope of this paper. When the reflection coefficients of the layers are small (and hence transmission is the leading phenomenon), the dominant terms associated with different layers remain reasonably distinct from one another. If the pulse width is sufficiently small, for each layer the corresponding dominant reflection is still proportional to the reflection coefficient of that layer.

In all the situations just stated, for a measurement location $(x,y,{z}_{0})$ with a fixed $z$ component and sampling time $t={t}_{0}$, the reflected signal is $\rho (x,y)u({t}_{0}+{z}_{0}/c)$, which is a constant multiple of $\rho (x,y)$. In other words, the reflected images should contain the same pattern as $\rho $. However, due to nonideal configurations such as a tilted or uneven sample, the effective measurement points are $(x,y,{z}_{0}+\u03f5(x,y))$, where $\u03f5(x,y)$ is a function of small magnitude—for instance, in the case of a tilted sample $\u03f5(x,y)={\alpha}_{1}x+{\alpha}_{2}y$, where ${\alpha}_{1}$ and ${\alpha}_{2}$ are small constants.

Technically, the demodulation problem of interest in this paper corresponds to extracting the sample structure $\rho (x,y)$, based on observing the reflected signal ${\mathrm{E\u20d7}}_{0}^{-}(x,y,z,t)$ at a fixed $z$ coordinate and time samples $t={t}_{1},\dots ,{t}_{M}$. Due to the nonideal configurations just stated, the returned signal, sampled at a time ${t}_{j}$, is in the form of $\rho (x,y){u}_{j}(x,y)+{n}_{j}(x,y)$, where ${u}_{j}(x,y)$ depends on the pulse waveform, sampling time, and $\u03f5(x,y)$. The additive term ${n}_{j}(x,y)$ models the noise uncertainty and undesired electromagnetic interactions. We will refer to the ${u}_{j}(x,y)$ factor as a sweep distortion profile, which is often a slowly varying function when the misconfiguration is small.

Specifically, inspired by the applications presented in this paper, we are interested in single or multilayer composite slabs of a binary nature where, for each layer,

Figures 1(c) and 1(d) show examples of the observed reflection from a binary slab (letter “M” printed on a card) at two different time instances. We are willing to make a binary characterization of the sample content by demodulating and excluding the multiplicative sweep distortion profiles affecting the desired image. In the remainder of this section, we will present an inversion scheme to characterize the slab contents based on such multiplicative observations.

## 3. METHODS

#### A. General Setup

For a more concise formulation, we use the more compact spatial notation $\mathit{x}=(x,y)$. Based on the discussions in the previous section, our observation is in the form of

The ultimate goal is recovering the image $\rho (\mathit{x})$ and the distortion profile ${u}_{j}(\mathit{x})$ based on the observation of ${y}_{j}(\mathit{x})$. Our strategy to address the problem is considering known subspace models for the signals ${u}_{j}(\mathit{x})$. More specifically, ${u}_{j}(\mathit{x})\in {\mathcal{S}}^{j}$, where ${\mathcal{S}}^{j}=\mathrm{span}({s}_{1}^{j}(\mathit{x}),\dots ,{s}_{{N}_{j}}^{j}(\mathit{x}))$ and ${N}_{j}$ is the $j$th subspace dimension. In Section A.3, we introduce an algorithm to construct a set of low-dimensional subspaces ${\mathcal{S}}^{1},\dots ,{\mathcal{S}}^{M}$ from the set of observations ${y}_{1}(\mathit{x}),\dots {y}_{M}(\mathit{x})$. Also, based on the structure of the problem, we may benefit from prior assumptions about $\rho (\mathit{x})$, as will be detailed in the remainder of this section.

### 1. Problem Reformulation as a Low-Rank Recovery Scheme

Consider discretizing the domain $D$ into a collection of pixels ${\{{\mathit{x}}_{i}\}}_{i=1}^{P}$. Based on the model in Eq. (4), a natural way of inverting the data for the image and sweep-distortion profiles is through the following minimization:

It is an interesting fact that under certain conditions on $\mathcal{A}$ (e.g., see [22]), the approximate solution of Eq. (10) can coincide with the solution of Eq. (9). Once a rank-one solution ${\mathit{X}}^{*}$ is available, the factors $\mathit{\alpha}$ and $\mathit{\beta}$ can be determined up to a constant multiple ($\mathit{\beta}={\mathit{X}}_{\text{:},1}^{*}$ and ${\mathit{\alpha}}^{T}={\mathit{X}}_{1,\text{:}}^{*}/{\mathit{X}}_{\mathrm{1,1}}^{*}$).

Other than the relaxation technique, there are other methods reported in the literature to address instances of Eq. (9). We are specifically interested in the alternating minimization approach [23] since it allows the incorporation of prior information into the reconstruction of each bilinear factor.

### 2. Alternating Minimization: A Maximum a Posteriori Framework

In minimizing the nonconvex program [Eq. (5)], a reliable technique would be to set up an alternating minimization scheme. Basically, by fixing any of the operands in the bilinear term in Eq. (5) (either $\rho $ or ${u}_{j}$), the problem turns into a standard least squares problem. The process would be to initialize one of the operands in the bilinear model, say $\rho $; with this quantity fixed, the resulting least squares in terms of ${u}_{j}$ is solved. By plugging in the acquired solutions for ${u}_{j}$’s, we can solve another least squares problem in terms of $\rho $ and proceed with such alternation until convergence.

This alternation approach is capable of approximating the solution to Eq. (9). Again, under certain incoherence conditions on $\mathcal{A}$ (slightly stricter than the ones stated for nuclear norm minimization [22,23]), the approximate solution coincides with the solution of Eq. (9). The major advantage of using an alternation scheme over the nuclear norm surrogate is the possibility of incorporating prior information beyond the subspace constraint on each factors of the problem in Eq. (4) (e.g., see [24]).

Basically, the proposed scheme corresponds to an alternation
between the maximum likelihood (ML) estimates of
$\rho $ and ${u}_{j}$. When some level of prior
knowledge about the $\rho $ and/or ${u}_{j}$ exists, the framework could be
generalized to an alternation between the maximum *a
posteriori* (MAP) estimates of
$\rho $ and ${u}_{j}$, as will be detailed in the
remainder of this section.

Consider $f(\rho ,{\{{u}_{j}\}}_{j=1}^{M}|{\{{y}_{j}\}}_{j=1}^{M})$ to be the joint probability density function (pdf) of $\rho $ and ${u}_{j}$ given the measurements ${y}_{j}$. The MAP estimates of $\rho $ and ${u}_{j}$ correspond to the maximizing parameters of the posterior distribution:

In the remainder of this section, we elaborate on how to implement the proposed scheme using certain facts about the problem. We would like to note that we make some reasonable assumptions and provide some slight modifications in implementing the steps outlined in Algorithm 1.

### 3. Iterative Reconstruction of the Sweep Distortion Profiles

In this section, we mainly focus on addressing the first stage of Algorithm 1 (line 3), which is determining the distortion profiles ${u}_{j}$ based on the observations ${\{{y}_{j}\}}_{i=1}^{M}$ and the knowledge of $\rho $.

To address this problem we only assume that, based on the physics of the problem, the signals ${u}_{j}$ reside in a low-dimensional subspace. In other words, ${\mathit{u}}_{j}={\mathit{S}}^{j}{\mathit{\alpha}}_{j}$ for the discrete representation. However, we make no prior assumptions about the coefficient vector ${\mathit{\alpha}}_{j}$. Since $\mathit{\rho}\odot {\mathit{u}}_{j}=\mathrm{diag}(\mathit{\rho}){\mathit{u}}_{j}$, we simply consider solving the least squares problems of the form

The choice of subspace plays a key role for this problem. In the remainder of this subsection, we present a one-time process which generates embedding subspaces ${\mathit{S}}^{j}$ based on the observations ${\mathit{y}}_{j}$.

From a technical standpoint, ${y}_{j}(\mathit{x})\approx \rho (\mathit{x}){u}_{j}(\mathit{x})$ where, as stated earlier, ${u}_{j}$ is a rather smooth function and, for our application, $\rho (\mathit{x})$ is a function of almost binary nature. Therefore, if we have a multi-scale representation of ${y}_{j}$, there should be a good overlap between ${u}_{j}$ and the course-level approximation of ${y}_{j}$. The high-frequency components of ${y}_{j}$ should be mainly in hold of $\rho $ and the noise.

Based on this argument, consider a wavelet representation of the signal ${y}_{j}$ as

### 4. Iterative Reconstruction of the Image Profile

Inspired by the second stage proposed in Algorithm 1, the main objective in this section is finding the MAP estimate for $\rho (\mathit{x})$. In order to maintain a more compact formulation and closed-form expressions, we make some independence assumptions about the pixel values as will be detailed in the remainder of this section.

Based on the binary model considered, we assume that each pixel in the image $\rho $ belongs to one of the two classes 0 or 1, and denote the pixel class by $\mathcal{C}(\mathit{x})$. Following the categorization in Eq. (3), the pixel values in class 0 concentrate around ${\rho}^{0}$ and the pixel values associated with class 1 concentrate around ${\rho}^{1}$. Consider $\mathit{\rho}={[{\rho}_{1},\dots ,{\rho}_{P}]}^{T}$ and $\mathcal{C}={[{\mathcal{C}}_{1},\dots ,{\mathcal{C}}_{P}]}^{T}$ to be vectors containing the pixel values and the class labels for the discrete image.

For our estimation problem, we are willing to address the maximization problem

whereGenerally speaking, for arbitrary random variables $\rho $, $y$, $u$, and $\mathcal{C}$ (discrete), using the Bayes’ rule we have

To further proceed with simplifying the MAP objective, we assume that the pixel values in $\mathit{\rho}$ are independent of each other. Based on this assumption and the i.i.d nature of the noise, we can state that

Clearly, based on Eq. (4), $f({y}_{j,i}|{\rho}_{i},{\mathcal{C}}_{i})\sim \mathcal{N}({\rho}_{i}{u}_{j,i}^{(k)},{\sigma}^{2})$, which is an immediate result of the normal noise model. We also assume that

*a priori*, a reasonable assumption is ${p}_{0}={p}_{1}=0.5$.

Finally, with reference to the term $f({\rho}_{i}|{\mathcal{C}}_{i})$, knowing that a pixel belongs to
class 0 (or 1), we know that its value is likely to be close to
${\rho}^{0}$ (or ${\rho}^{1}$). In practice,
$f({\rho}_{i}|{\mathcal{C}}_{i})$ should model the uncertainty on
how the values of each class concentrate around the mean class
value. While we can use different distributions to model such
concentration, for simplicity and in order to obtain closed-form
expressions, we assume that $f({\rho}_{i}|{\mathcal{C}}_{i}=0)$ and $f({\rho}_{i}|{\mathcal{C}}_{i}=1)$ are in the form of
*truncated normal distributions*, taking
positive arguments and mainly concentrating around
${\rho}^{0}$ and ${\rho}^{1}$, respectively. By definition, a
random variable $\rho $ truncated to values greater than
zero takes the pdf

Figure 2 shows the underlying truncated distributions. The positivity assumption in Eq. (17) is imposed by the physics of the problem, where the reflection coefficient of the slab is always considered to be a positive quantity.

Having the constituting terms modeled in Eq. (16), we can proceed with the pixel value and class label MAP estimation. To maximize $G(\mathit{\rho},\mathcal{C})$, we may minimize the negative log function

Based on the proposed distribution models, one can easily verify that

While the parameters ${\sigma}_{0}$ and ${\sigma}_{1}$ have statistical interpretations, from an optimization standpoint they can be considered as free parameters controlling the algorithm’s performance. When ${\sigma}_{0}$ and ${\sigma}_{1}$ are set to be small quantities, the algorithm converges faster at the expense of more sensitivity to the initialization (possibility of recovering a local minimizer). Assigning relatively larger values to these quantities paves the path toward identifying the global minimizer in a greater number of iterations. This is also a more reliable parameter selection in the case of noisy observations.

## 4. SIMULATION AND EXPERIMENTS

In this section, we assess the performance of the proposed technique in demodulating simulated and real data. The simulation results mainly highlight the performance of the method for various numbers of frames and noise levels in the data. The second set of experiments demonstrates the method’s success in removing multiplicative sweep distortion profiles from actual THz measurements.

#### A. Simulated Data

We simulate the reflected signal from a dielectric slab with the $x$-$y$ profile depicted in Fig. 3(a). The dielectric width is $d=100\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}$. For the binary dielectric slab, the reflection coefficients are ${\rho}^{0}=0.3$ and ${\rho}^{1}=0.1$. The sample is probed with a bipolar THz pulse (simply derivative of a Gaussian) as

In recovering the binary profile, for the subspace construction step associated with this example we use ${N}_{1}=\dots ={N}_{10}=100$ most dominant wavelet basis. We use the symlet wavelet family consistently through this example and the remaining experiments.

Figure 3(c.1) reports the reconstructed reflection profile $\rho $ using our proposed scheme. The algorithm converges in less than 20 iterations and, as may be observed, the result is no more in hold of the distortion traces. A simple rounding of the result to the closest binary value (between ${\rho}^{0}$ and ${\rho}^{1}$) yields Fig. 3(c.2), which is reasonably close to the original profile depicted in Fig. 3(a).

Figure 3(d.1) reports the reconstruction using the nuclear norm framework proposed in Section A.1. While some level of distinction between the two binary phases is achieved, the fact that our proposed scheme efficiently uses the binary prior helps us outperform the nuclear norm reconstruction. A similar binary approximation of the nuclear norm reconstruction is demonstrated in Fig. 3(d.2), which clearly fails to characterize some of the main details in the original binary profile.

To more extensively analyze the performance of our algorithm, we
proceed by the sensitivity analysis of the main parameters affecting
the reconstruction. For this purpose, two main scenarios are
considered. First, an ideal case where the exact sweep distortion
subspace is known *a priori*. From a theoretical
standpoint, one way to characterize a subspace for the sweep
distortion profiles is to run exactly the same experiment with a
homogeneous slab and take the observations as the subspace basis.
While experimentally such accurate characterization of the subspace is
hard, here we present it as a hypothetical case for comparison
purposes and to compare the results against an ideal setup.

For the second class of experiments, we use the proposed wavelet thresholding framework to determine the distortion subspaces. For each proposed scenario, we test the algorithm for various values of $M$ (number of available frames) and the SNR in the observations. Both the reconstructed $\rho $ as well as a binary rounded version are compared with the reference image depicted in Fig. 3(a), and the mean squared error (MSE) is reported.

Figure 4 shows the MSE values for various numbers of available frames and a fixed SNR value of 10 dB. The dark solid curve corresponds to the case of exact subspace characterization and the dashed curve standing close to it is the MSE value for the reconstructed $\rho $ rounded to the closest binary value. We can see that for only 7 frames (or more), a perfect or close to perfect recovery of the binary profile is possible. For fewer numbers of available frames, the MSE still remains at a controlled level.

The lighter solid line (and the dashed line standing close to it) report the MSE values of similar experiments when the sweep distortion subspaces are characterized through the wavelet thresholding. We can see that, aside from a slight performance gap, the MSE values follow a similar reduction pattern as the ideal case. Based on the reported MSE values, the recovery is still close to the reference image. The slight performance gap between this case and the ideal case is the result of error in exact characterization of the sweep distortion subspaces.

Figure 5 shows the MSE
values for the case of a fixed number of frames
($M=10$) and varying observation noise. For
this experiment, in the case of known distortion subspace, a perfect
recovery of the sweep distortion profile is possible for SNR values
approximately exceeding 10 dB (by rounding to the closest
binary value). In the case of characterizing the sweep distortion
subspace through wavelet thresholding, the MSE values follow a similar
trend as the case of knowing the distortion subspace *a
priori*. The only difference is a slight performance gap
which is again linked to the error in exact characterization of the
sweep distortion subspaces.

#### B. Experimental Data

To experimentally apply the demodulation process, we used a standard THz-TDS system in reflection geometry. We exploited a bipolar pulse with 2.0 THz bandwidth and an effective center frequency of 1.0 THz. In the first experiment (Fig. 6), we used a binary sample with metallic surface and a cross-shaped pattern carved into a $15\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 15\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ metal surface with 2 mm depth. As noted in Figs. 6(a.1)–6(a.3), the spatial sweep distortions are dominantly present in the time domain observations. For this experiment, $M=25$ frames were available, three of which are shown here.

Figures 6(b.1)–6(b.3) and 6(c) report the decoupled sweep distortions from the binary cross profile. We can observe a good characterization of the cross profile thanks to the large number of available frames and high SNR observations. A more challenging experiment with fewer available frames and noisier data is presented next.

A more challenging experiment uses a three-layer paper sample with the letters M-I-T printed on each page from front to back, respectively. The page dimensions are $15\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}\times 30\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ and the thicknesses are 3.0 μm, stacked flush next to each other ($\sim 30\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{\mu m}$ gap), similar to the structure of a book. Figures 7(a.1)–7(a.3) show the reflected signal from the first layer at three different instances of time. Figures 7(b.1)–7(b.3) and 7(c.1)–7(c.3) show reflection samples from the second and third layers, respectively. The available number of frames corresponding to the first, second, and third layer are $M=25$, 6, 8, respectively.

Followed by the observations in each column of Fig. 7, the demodulation results are presented for the three letters at different pages of the sample, as the THz pulse travels deeper into the pages. The rows 4–6 in Fig. 7 correspond to the sweep distortion profiles decoupled from the observed data. The last row corresponds to the recovered character profiles.

For this multilayer setup, the SNR drops after the first layer and a larger number of inter-reflections is induced; thus, the recovered letters “I” and “T” in panels Figs. 7(b.7) and 7(c.7) have more noise compared to the recovered letter “M” on the first page [panel Fig. 7(a.7)].

The clean results from the first experiment (Fig. 6) on the metallic surface are anticipated due to high SNR and the perfectly flat surface of the polished steel; however, for the paper, the reflection is rather small ($<0.15$), the distortion contrast is at the same level of the signal contrast itself, and the pages have slight depth variations across them. The results from the paper sample in Fig. 7 indicate the true robustness of the decoupling technique for practical samples.

It is worth noting that the major burden in characterizing the sample inhomogeneity profiles (characters in this example) is the varying signal level in the observed images. For instance, while the contrast between the character “M” and the background can be visually detected in Fig. 7(a.3), the pixel values drastically vary from one portion of the letter to the other. As a result, characterizing the letter based on the pixel values becomes an inconceivable task. However, once a binary profile is reconstructed through the proposed algorithm, other post-processing schemes can be employed to characterize the reconstructions, especially in the case of the more noisy reconstructions such as Fig. 7(c.7). The interested reader is referred to [26–29], where advanced shape-composition techniques are employed to accurately identify binary inclusions in noisy images. This class of techniques can accurately characterize the inclusion in an image like Fig. 7(c.7). Even when the noise is higher, parts of the object are missing due to signal occlusions, and we have clutter or overlapping characters caused by the inter-reflections from neighboring layers.

## 5. CONCLUSION

The main outcome of this research is the development of a general demodulation scheme to process reflected signals conveniently applicable to THz imaging. While a careful processing of the THz data can produce high-resolution images thanks to its high frequency, dealing with phenomena such as inter-reflections, noise, and modulated distortions is inevitable. We showed that reformulation of the problem as a bilinear inverse problem and making practical assumptions, such as the binary prior, can assist us with a successful inversion of THz measurements in THz-TDS systems.

Using basic statistical assumptions about the image and the observation, we were able to develop an iterative inversion scheme which is computationally efficient, distributable, and scalable to big data sets. In fact, the main computational load at every step of the algorithm is a least squares solve. The remaining computations are in the form of thresholding-type operations.

Since the main objective in many applications (such as water profilometry, extracting coded structures, and extracting subsurface cracks) is a binary characterization of the inhomogeneity, we used a binary prior in the modeling in our samples. The algorithm can be naturally generalized to the case of polyadic priors, where more than two levels are considered for the image profile. Basically, the proposed algorithm can be further modified or combined with post- or pre-processing tools to handle a larger class of imaging applications.

## Funding

Massachusetts Institute of Technology Media Laboratory Consortia (2746038); National Science Foundation (NSF) (RI 1527181).

## REFERENCES

**1. **Y.-C. Shen, “Terahertz pulsed
spectroscopy and imaging for pharmaceutical applications: a
review,” Int. J. Pharm. **417**, 48–60
(2011). [CrossRef]

**2. **M. Tonouchi, “Cutting-edge
terahertz technology,” Nat.
Photonics **1**, 97–105
(2007). [CrossRef]

**3. **C. Seco-Martorell, V. López-Domnguez, G. Arauz-Garofalo, A. Redo-Sanchez, J. Palacios, and J. Tejada, “Goya’s
artwork imaging with terahertz waves,”
Opt. Express **21**, 17800–17805
(2013). [CrossRef]

**4. **A. Redo-Sanchez, N. Laman, B. Schulkin, and T. Tongue, “Review of terahertz
technology readiness assessment and
applications,” J. Infrared Millimeter
Terahertz Waves **34**, 500–518
(2013). [CrossRef]

**5. **B. Heshmat, H. Pahlevaninezhad, Y. Pang, M. Masnadi-Shirazi, R. Burton Lewis, T. Tiedje, R. Gordon, and T. E. Darcie, “Nanoplasmonic
terahertz photoconductive switch on GaAs,”
Nano Lett. **12**, 6255–6259
(2012). [CrossRef]

**6. **B. Heshmat, H. Pahlevaninezhad, and T. Darcie, “Carbon
nanotube-based photoconductive switches for THz detection: an
assessment of capabilities and
limitations,” IEEE Photon. J. **4**, 970–985
(2012). [CrossRef]

**7. **B. Heshmat, M. Masnadi-Shirazi, R. B. Lewis, J. Zhang, T. Tiedje, R. Gordon, and T. E. Darcie, “Enhanced terahertz
bandwidth and power from GaAsBi-based
sources,” Adv. Opt. Mater. **1**, 714–719
(2013). [CrossRef]

**8. **W. L. Chan, K. Charan, D. Takhar, K. F. Kelly, R. G. Baraniuk, and D. M. Mittleman, “A single-pixel
terahertz imaging system based on compressed
sensing,” Appl. Phys.
Lett. **93**, 121105 (2008). [CrossRef]

**9. **W. Withayachumnankul and D. Abbott, “Terahertz imaging:
compressing onto a single pixel,”
Nat. Photonics **8**, 593–594
(2014). [CrossRef]

**10. **C. Li, J. Grant, J. Wang, and D. R. Cumming, “A nipkow disk
integrated with Fresnel lenses for terahertz single pixel
imaging,” Opt. Express **21**, 24452–24459
(2013). [CrossRef]

**11. **H. Shen, L. Gan, N. Newman, Y. Dong, C. Li, Y. Huang, and Y. Shen, “Spinning disk for
compressive imaging,” Opt.
Lett. **37**, 46–48
(2012). [CrossRef]

**12. **A. W. Lee and Q. Hu, “Real-time,
continuous-wave terahertz imaging by use of a microbolometer
focal-plane array,” Opt.
Lett. **30**, 2563–2565
(2005). [CrossRef]

**13. **Y. Chen, S. Huang, and E. Pickwell-MacPherson, “Frequency-wavelet
domain deconvolution for terahertz reflection imaging and
spectroscopy,” Opt.
Express **18**, 1177–1190
(2010). [CrossRef]

**14. **R. Galvão, S. Hadjiloucas, J. Bowen, and C. Coelho, “Optimal
discrimination and classification of THz spectra in the wavelet
domain,” Opt. Express **11**, 1462–1473
(2003). [CrossRef]

**15. **Y. Deng, Q. Sun, F. Liu, C. Wang, and Q. Xing, “Terahertz
time-resolved spectroscopy with
wavelet-transform,” in *3rd
International Congress on Image and Signal Processing
(CISP)* (IEEE,
2010), Vol. 7,
pp. 3462–3464.

**16. **A. D. Burnett, W. Fan, P. C. Upadhya, J. E. Cunningham, M. D. Hargreaves, T. Munshi, H. G. Edwards, E. H. Linfield, and A. G. Davies, “Broadband terahertz
time-domain spectroscopy of drugs-of-abuse and the use of
principal component analysis,”
Analyst **134**, 1658–1668
(2009). [CrossRef]

**17. **X. Yin, B. W.-H. Ng, and D. Abbott, *Terahertz Imaging for Biomedical
Applications: Pattern Recognition and Tomographic
Reconstruction*
(Springer,
2012).

**18. **A. Aghasi, I. Mendoza-Sanchez, E. L. Miller, C. A. Ramsburg, and L. M. Abriola, “A geometric
approach to joint inversion with applications to contaminant
source zone characterization,”
Inverse Prob. **29**, 115014 (2013). [CrossRef]

**19. **A. Aghasi, M. Kilmer, and E. L. Miller, “Parametric level
set methods for inverse problems,”
SIAM J. Imaging Sci. **4**, 618–650
(2011). [CrossRef]

**20. **E. L. Miller, L. M. Abriola, and A. Aghasi, “Environmental
remediation and restoration: hydrological and geophysical
processing methods,” IEEE Signal
Process. Mag. **29**(4),
16–26
(2012). [CrossRef]

**21. **R. W. Scharstein, “Transient
electromagnetic plane wave reflection from a dielectric
slab,” IEEE Trans. Educ. **35**, 170–175
(1992). [CrossRef]

**22. **B. Recht, M. Fazel, and P. A. Parrilo, “Guaranteed
minimum-rank solutions of linear matrix equations via nuclear norm
minimization,” SIAM Rev. **52**, 471–501
(2010). [CrossRef]

**23. **P. Jain, P. Netrapalli, and S. Sanghavi, “Low-rank matrix
completion using alternating minimization,”
in *Proceedings of the Forty-Fifth Annual ACM Symposium on
Theory of Computing*
(ACM, 2013),
pp. 665–674.

**24. **K. Lee, Y. Wu, and Y. Bresler, “Near optimal
compressed sensing of sparse rank-one matrices via sparse power
factorization,” arXiv:1312.0525
(2013).

**25. **G. Strang and T. Nguyen, *Wavelets and Filter Banks*
(SIAM,
1996).

**26. **A. Redo-Sanchez, B. Heshmat, A. Aghasi, S. Naqvi, M. Zhang, J. Romberg, and R. Raskar, “Terahertz
time-gated spectral imaging for content extraction through layered
structures,” Nat. Commun.
(to be published).

**27. **A. Aghasi and J. Romberg, “Convex cardinal
shape composition,” SIAM J. Imaging
Sci. **8**, 2887–2950
(2015). [CrossRef]

**28. **A. Aghasi and J. Romberg, “Sparse shape
reconstruction,” SIAM J. Imaging
Sci. **6**, 2075–2108
(2013). [CrossRef]

**29. **A. Aghasi and J. Romberg, “Learning shapes by convex
composition,” arXiv:1602.07613
(2016).