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
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 , and nondestructive testing  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  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 and the index sets and , we use to denote the matrix obtained by restricting the rows of to . Similarly, denotes the restriction of to the columns specified by , and is the submatrix with the rows and columns restricted to and , 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 direction, propagating along the direction in the free space. The waveform is considered to be a finite-duration THz pulse , for which the traveling field along every point is1(a) and 1(b)].
For a homogeneous layer of width , with reflection coefficient and refraction index , the returned signal can be analytically expressed by convolving the pulse with a train of impulse functions with decaying coefficients. More specifically ,2) corresponds to a reflection. Particularly, the first term corresponds to the reflection from the front surface of the slab, the second term () corresponds to the back surface reflection, and the remaining terms correspond to the returned waves after a number of inter-reflections within the slab.
The response in Eq. (1) can still be reasonably accurate when the dielectric slab is homogeneous along the axis, i.e., for . It only requires plugging the point-wise values of and in the formulation.
In the case of perfect or close to perfect reflection, where , the first impulse term in Eq. (1) dominates all the other terms and the returned wave is a simple modulation of with the pulse waveform. When , the most dominant terms in 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 is large, the first reflected waveforms do not overlap with the subsequent reflections. In this case, still a modulation of 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 with a fixed component and sampling time , the reflected signal is , which is a constant multiple of . In other words, the reflected images should contain the same pattern as . However, due to nonideal configurations such as a tilted or uneven sample, the effective measurement points are , where is a function of small magnitude—for instance, in the case of a tilted sample , where and are small constants.
Technically, the demodulation problem of interest in this paper corresponds to extracting the sample structure , based on observing the reflected signal at a fixed coordinate and time samples . Due to the nonideal configurations just stated, the returned signal, sampled at a time , is in the form of , where depends on the pulse waveform, sampling time, and . The additive term models the noise uncertainty and undesired electromagnetic interactions. We will refer to the 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,1(b)].
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.
A. General Setup
For a more concise formulation, we use the more compact spatial notation . Based on the discussions in the previous section, our observation is in the form of
The ultimate goal is recovering the image and the distortion profile based on the observation of . Our strategy to address the problem is considering known subspace models for the signals . More specifically, , where and is the th subspace dimension. In Section A.3, we introduce an algorithm to construct a set of low-dimensional subspaces from the set of observations . Also, based on the structure of the problem, we may benefit from prior assumptions about , as will be detailed in the remainder of this section.
1. Problem Reformulation as a Low-Rank Recovery Scheme
Consider discretizing the domain into a collection of pixels . 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:4). Consider to be a matrix representation of the subspace and to be a subspace in which the signal lives (in the case of no such assumption, can be simply the identity matrix spanning the canonical basis). Under this assumption, the minimization in Eq. (5) can be cast as 6) is nonlinear in terms of and . However, since the outer product of two nonzero vectors is a rank-one matrix, an equivalent reformulation of Eq. (6) is the minimization 9) makes it a combinatorial problem. As a remedy, a well-known approximation strategy is to use the nuclear norm of , denoted by , as a convex proxy to :
It is an interesting fact that under certain conditions on (e.g., see ), the approximate solution of Eq. (10) can coincide with the solution of Eq. (9). Once a rank-one solution is available, the factors and can be determined up to a constant multiple ( and ).
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  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 or ), the problem turns into a standard least squares problem. The process would be to initialize one of the operands in the bilinear model, say ; with this quantity fixed, the resulting least squares in terms of is solved. By plugging in the acquired solutions for ’s, we can solve another least squares problem in terms of 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 (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 ).
Basically, the proposed scheme corresponds to an alternation between the maximum likelihood (ML) estimates of and . When some level of prior knowledge about the and/or exists, the framework could be generalized to an alternation between the maximum a posteriori (MAP) estimates of and , as will be detailed in the remainder of this section.
Consider to be the joint probability density function (pdf) of and given the measurements . The MAP estimates of and 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 based on the observations and the knowledge of .
To address this problem we only assume that, based on the physics of the problem, the signals reside in a low-dimensional subspace. In other words, for the discrete representation. However, we make no prior assumptions about the coefficient vector . Since , we simply consider solving the least squares problems of the form12) is 12)].
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 based on the observations .
From a technical standpoint, where, as stated earlier, is a rather smooth function and, for our application, is a function of almost binary nature. Therefore, if we have a multi-scale representation of , there should be a good overlap between and the course-level approximation of . The high-frequency components of should be mainly in hold of and the noise.
Based on this argument, consider a wavelet representation of the signal as25]. For a fixed subspace order , our proposed subspace selection for corresponds to the top scaling/wavelet basis functions with the largest coefficients (in magnitude). In other words, we sort the wavelet coefficients in descending magnitude order and select the basis functions associated with the top coefficients. The selected basis (in discrete form) corresponds to the columns of . This process is performed only once during the entire reconstruction and once the subspaces are determined, they remain fixed throughout the alternation scheme.
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 . 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 belongs to one of the two classes 0 or 1, and denote the pixel class by . Following the categorization in Eq. (3), the pixel values in class 0 concentrate around and the pixel values associated with class 1 concentrate around . Consider and 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
Generally speaking, for arbitrary random variables , , , and (discrete), using the Bayes’ rule we have13), the vectors are known and deterministic and our intention to derive Eq. (14) is to make the MAP estimation in terms of probability density functions which all make such prior assumption. Also, since the denominator term in Eq. (14) is independent of and , in Eq. (15) we neglected the corresponding term as it will be a constant factor in the maximization.
To further proceed with simplifying the MAP objective, we assume that the pixel values in are independent of each other. Based on this assumption and the i.i.d nature of the noise, we can state that16) can now be modeled based on the problem setup.
Clearly, based on Eq. (4), , which is an immediate result of the normal noise model. We also assume that
Finally, with reference to the term , knowing that a pixel belongs to class 0 (or 1), we know that its value is likely to be close to (or ). In practice, 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 and are in the form of truncated normal distributions, taking positive arguments and mainly concentrating around and , respectively. By definition, a random variable 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 , we may minimize the negative log function
Based on the proposed distribution models, one can easily verify that15). The minimization objective is separable in terms of the cost for each pixel and we can minimize each term individually. Also, minimizing can be performed in a serial manner by first minimizing with respect to and then with respect to the class label . More concretely, 18), a simple calculation yields 20) if feasible, or 0 otherwise. Basically, 20) in some sense generates the optimal solution , but leaves us with an ambiguity about the class label. The label can be obtained based on the comparison just shown to find the ultimate MAP estimator for the image profile. Algorithm 2 summarizes the steps sketched in this section and presents the overall process to perform the demodulation task.
While the parameters and have statistical interpretations, from an optimization standpoint they can be considered as free parameters controlling the algorithm’s performance. When and 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 - profile depicted in Fig. 3(a). The dielectric width is . For the binary dielectric slab, the reflection coefficients are and . The sample is probed with a bipolar THz pulse (simply derivative of a Gaussian) as3(b) are three sample images where an almost vertically moving sweep distortion across the frames is observable. The samples are synthetically corrupted with white noise ().
In recovering the binary profile, for the subspace construction step associated with this example we use 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 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 and ) 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 (number of available frames) and the SNR in the observations. Both the reconstructed 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 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 () 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 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, 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 and the thicknesses are 3.0 μm, stacked flush next to each other ( 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 , 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 (), 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.
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.
Massachusetts Institute of Technology Media Laboratory Consortia (2746038); National Science Foundation (NSF) (RI 1527181).
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).