1.

Introduction

In recent years, the use of terahertz (10¹² Hz) radiation has been the subject of much research in areas as diverse as semiconductor and inorganic glass characterization,^{1, 2, 3, 4} pharmaceuticals and solid state chemistry,^{5, 6, 7, 8} catalysis and carbonaceous materials,^{9, 10, 11, 12} and process quality control and medicine.^{13, 14, 15} Of particular interest has been a technique known as terahertz time-domain-spectroscopy (THz-TDS) and its associated quasi three-dimensional analogue: terahertz pulsed imaging (TPI).¹⁶ THz-TDS and TPI both use broadband (typically 0.1to 3 THz) sources of pulsed radiation, and the time-varying electric field is measured, from which both frequency dependent amplitude and phase information can be derived. This allows the frequency dependent complex dielectric permittivity to be directly calculated at terahertz frequencies. In addition to this useful ability, pulsed sources afford good axial resolution (30 μm has been reported,¹⁷ however this is dependent on the refractive index of the medium). Although in-plane resolution is classically diffraction limited to approximately 200 μm (at 1.5 THz, dependent on frequency), there have been numerous studies exploring near-field techniques to overcome this limitation and increase resolution.^{18, 19, 20, 21}

Within the area of medical imaging and spectroscopy, research into the use of pulsed sources of terahertz radiation has been driven by a number of key considerations that make this technology of particular interest to the medical field. First, the energy of terahertz frequency photons is low (especially compared to high frequency radiations such as ultraviolet and x-rays) making it nondestructive and nonionizing.²² However, heating effects could still be an issue, and so safe levels of exposure have been set in the milliwatt regime.^{23, 24, 25} This level is orders of magnitude above the typical powers of pulsed terahertz systems, where common average output powers are reported to be on the order of 100 nW.²⁶ Intermolecular bonds, such as hydrogen bonds, have been found to have resonances extending into the terahertz region. Furthermore, the dielectric relaxation of polar liquids that dominates the microwave response also extends into this region. Therefore, water has a significant dielectric response in the terahertz region and absorbs the radiation efficiently.^{27, 28, 29} Additionally, recent work by Havenith and co-workers has shown that free, bulk, and interfacial water appear to have different optical properties at terahertz frequencies, which they probed by investigating the protein-water interaction,^{30, 31} which is a significant result in terms of biological applications of terahertz radiation.

Most biological tissues have high water content and so the terahertz absorption of the polar liquid through relaxational and resonant responses, coupled with the low power of the pulsed terahertz sources, limits the penetration depth of terahertz radiation into biological tissue to the order of a few millimeters at best. The precise value is dependent upon both the signal-to-noise ratio (SNR) of the system and the specific tissue type under investigation; for example, the penetration depth for adipose tissue is over 3 mm for a system with an SNR of 1000 at a frequency of 1 THz, whereas for skin tissues this typically drops to 500 μm at most.³² Such characteristics limit any in vivo studies to utilize reflection geometries only. While this may appear initially to be a significant drawback, it is this high sensitivity to water (and other intermolecular interactions) that makes terahertz radiation of such interest in the medical field as it can be used to detect subtle changes in the tissue composition, in particular, in epithelial tissues which are not detectable with techniques such as MRI. It is these changes in the structure and composition of tissues that can be relevant for identifying abnormalities such as cancer. For example, pulsed terahertz images of freshly excised breast cancer have been able to show contrast between healthy and cancerous tissue, even when the cancer has been of the in situ noncalcified form.³³ This is of particular interest as noncalcified tumors are often missed during breast conserving surgery (as often they are not picked up on x-rays and they are not palpable).

Pulsed terahertz imaging/spectroscopy can be performed in two basic geometries; transmission and reflection. Because of the time domain nature of the measurements, both time and frequency analyses can be employed. By means of the Fourier transform, the time-domain data can yield the absorption coefficient and refractive index of the sample, and they have been measured for several tissue types including breast cancer,³⁴ skin cancer,³⁵ colon cancer,³⁶ liver cirrhosis,³⁷ scar tissue from burning,³⁸ and corneal tissues.^{38, 39} It was found that for all these tissue types that the abnormal tissue has a higher absorption coefficient than the corresponding healthy tissue at terahertz frequencies. Additionally, for the liver tissue a study was devised to quantify how much of the difference is due to changes in water content, and how much is due to structural changes. It was found that for frequencies above 0.4 THz, structural changes contributed to over half of the change in absorption coefficient, and thus water was not the sole source of contrast in the terahertz properties of liver cirrhosis.³⁷ This is a very significant finding and opens up further avenues of research as previously it had been thought that water was the dominant source of image contrast.

In time-domain analyses, the data can be thought of as being analogous to an ultrasound A-scan. Reflections of the terahertz transient off different layers result in an observable pulse train that can be used to determine the structure of the sample at various depths. Some of the present authors have investigated the ability of pulsed terahertz radiation to observe the different layers of the skin using a reflection system, where it was possible in some cases to resolve the pulses from the skin surface and from the stratum corneum/epidermis interface which can have important implications in skin hydration.⁴⁰ Time-domain analysis has also been used in studies of skin cancer with terahertz radiation where it was found that cancerous areas exhibited larger reflections of the terahertz pulse, indicative of their increased refractive index values when compared to their healthy analogues.⁴¹ Most biological tissues have a high water content and so both the terahertz resonance of the hydrogen bonds and dielectric relaxation of the polar molecules mean that the penetration depth of terahertz radiation into biological tissue is very shallow (of the order of millimeters at most depending on the tissue type). As a result, in vivo imaging is performed in reflection geometry. Since skin is the most easily accessible organ of the body, we have been using it to conduct in vivo experiments with a handheld probe designed by Teraview Ltd (Cambridge, UK). In order to extract as much useful information as possible from the probe scans there is a need for improved data processing techniques to increase the axial resolution and resolve closer interfaces. In this paper various data processing techniques for such a probe system are presented and explored, using data taken from human skin in vivo as examples.

2.

Materials and Methods

A handheld probe system [Fig. 1a] was used to take the terahertz in vivo measurements of human skin reported in this paper. The system was originally envisaged to facilitate intrasurgical measurements in the operating theatre and so its design has prioritized flexible positioning and fast data acquisition. The handheld probe consists of a molded plastic outer casing with a 1.5-cm long quartz window [Fig. 1b]. Inside the probe are two photoconductive devices (the emitter and detector) coupled via a fiber-optic to the 800 nm ultrafast pulsed laser (Vitesse, Coherent, Santa Clara, California) contained within the main body of the unit. The emitted terahertz radiation passes through a Risley beam steerer before being guided by the silicon waveguide onto the quartz plate, with the reflected signal following a similar route to the detector [Fig. 1c]. The Risley beam steerer scans the terahertz pulse along the long axis of the quartz plate, allowing one-dimensional line scans to be produced of a sample placed on the quartz plate. This geometry allows the probe to be “wanded” over areas of interest; as a result the probe gives additional flexibility to sample geometries and sites, but at the trade-off of decreased SNR and dynamic range (DR) due to the losses induced by the fiber optics and the multiple reflections off the components in the terahertz beam path [Fig. 1c]. As an example, the frequency-domain DR (as defined in Ref. 42) is on average around 13 dB lower across the usable frequency range when compared to a more conventional reflection time-domain system (see Fig. 2).

Fig. 1

(a) Photograph of the terahertz imaging probe system from TeraView Ltd., UK, the laser light used to generate the terahertz is fed along an umbilical from the main unit to the probe head. (b) The imaging window is made of quartz and is 1.5-cm long. (c) Schematic of the terahertz optics within the probe head. The Risley prisms counter rotate at the same speed, causing the terahertz beam to transcribe a line scan along the quartz window.

Fig. 2

DR for a conventional “flatbed” terahertz reflection imaging system (black solid line) compared to the probe system (dashed red line) measured from the reflection off a mirror reference. The probe system displays both a lower bandwidth and DR when compared to the conventional imaging system.

Skin is composed of two primary layers, the epidermis and the dermis. The epidermis itself is subdivided into five further layers, the uppermost layer being the stratum corneum, a layer of dead cells with a low water content varying in thickness from 10 to 200 μm depending on body site, being the thickest on the soles and palms. It is the thickness and composition of these layers that determines the form of the measured terahertz waveform due to the superposition of reflected pulses from the interfaces [see Fig. 3a]. As a result of the unique geometry of the probe, our group has needed to develop a number of processing steps (detailed below) which have allowed us to extend the earlier in vivo studies of skin^{40, 43} to include a variety of skin types from different regions of the body. In particular, this was made possible due to the ease of positioning the terahertz probe in a variety of different positions [see Fig. 3b].

Fig. 3

(a) Schematic of the layers of skin, along with the expected reflections from each interface. (b) Photograph highlighting the increased flexibility of the probe system.

2.2.

Improved Baseline Subtraction

Despite its usefulness regarding signal alignment, the reflection off the inner silicon-quartz interface causes interference with the main peak as the transient pulse resulting from this interface is still nonzero when the second reflection reaches the detector. Figure 4a shows a typical unprocessed time domain waveform recorded by the terahertz probe. Pulse A is the reflection from the lower silicon-quartz interface and it is this pulse that is used to align the measurements. Pulse B is the reflection from the quartz-sample interface (in this case the sample is air) and it is this pulse that is used to extract the complex dielectric permittivity of the sample. To extract accurate material properties, the reflection from the sample needs to be isolated and thus it is important to remove any contribution from the first pulse. This unwanted contribution is called the “baseline.”

Fig. 4

Time domain spectroscopy data for isopropanol. (a) Representative sample time-domain waveform from the probe system (unprocessed data). (b) The quartz-sample reflected waveform before (solid light line) and after (dashed dark line) baseline extraction. The subtracted baseline is also shown (solid dark line). (c) Absorption coefficient and (d) refractive index data for isopropanol calculated using reflection measurements both before (solid light line) and after (solid dark line) baseline extraction. The transmission spectrum (dashed dark line) and literature data from Yomogida (Ref. 44) (crosses) are shown for comparison.

We have devised a method to calculate a baseline to remove unwanted reflections and hence, transient pulse features from the sample reflection. A measurement of water is made as well as a reference measurement of air and is described in more detail in Ref. 45. Since water makes perfect contact with the quartz window and its terahertz properties are well known, this measurement can be used to deduce the system baseline. The effect of the baseline extraction in the time domain pulse is shown in Fig. 4b along with a typical baseline response. The baseline subtraction flattens out the terahertz response post peak, with the processed waveform from 19 ps appearing noticeably flatter than the unprocessed waveform, which is as expected for a terahertz pulse of this type. The effect on the frequency-dependent optical constants is pronounced, as shown in Figs. 4c, 4d, where the absorption coefficient of isopropanol recorded using the reflection system agree far more closely to the transmission system data across the frequency range after utilization of the baseline subtraction method.

3.

Results and Discussion

3.1.

Comparison of Skin Properties from Different Areas

To exploit the probe flexibility, the terahertz time-domain response from five different body sites was measured. The five sites were the forehead, cheek, chin, dorsal forearm, and the palm of the hand. Ten subjects, 7 male and 3 female between the ages of 24 and 29 were measured over 5 days at each of the 5 different skin sites. Figure 5 shows the average terahertz impulse responses recorded from the five positions using the hand-held probe, after performing a double Gaussian deconvolution (DGIF) on the recorded time-domain waveforms.

Fig. 5

Average terahertz impulse response recorded by the hand-held terahertz probe for the five different skin sites. The error bars represent the statistical deviation of the sample dataset of 50 measurements taken from 10 subjects during the 5 day period.

The impulse responses for the three positions on the face are quite similar, with a small peak in the signal followed by a large trough. Similarly, the impulse response recorded for the dorsal forearm is similar in shape, but with a more pronounced peak. In contrast, the impulse response measured from the palm shows a distinct second feature; a further trough appearing after the main trough. For the case of the palm it appears that the stratum corneum is thick enough for the reflections from the two interfaces [see Fig. 3a] to be resolved using the waveforms processed using a DGIF. For the other four sites, variation in the peak heights and troughs was observed. In a previous in vivo study of the dorsal forearm⁴⁰ it was postulated that the height and appearance of a small peak before the main trough was indicative either of a dry surface layer of skin, existence of an air gap due to hairs, or a combination of both. This is because dry skin (and air) has a lower refractive index than that of quartz and therefore no phase change is observed in the pulse upon reflection. If the size of the small peaks observed in the measurements is related to surface dryness it would suggest that the dorsal forearm was drier than the other sites measured, and the chin was the wettest surface on the face. Previous hydration studies of the face and volar forearm found that the hydration level of the forehead, cheek, chin, and volar forearm areas were 77.9±6.6, 70.4±7.1, 76.1±6.5, and 61.9±4.1, respectively.⁴⁶ These appear to agree with the relative peak heights for the terahertz measurements, as the dorsal forearm (not the volar forearm in the present case) has a much larger peak height than the facial sites. In addition, a qualitative analysis of the terahertz impulse responses for the three facial sites would suggest that the cheek is drier than the forehead and chin areas, in agreement with the previously reported hydration study. However, quantitative analysis is not possible with these datasets due to the large statistical variation within both the terahertz and hydration measurements. Recent work by Brown and co-workers have used a Bruggeman mixing model to calculate the hydration level of ex vivo porcine corneal tissues.^{39, 47} They demonstrated a noise equivalent water concentration sensitivity of 0.19% for their terahertz spectroscopy system (TPI3000 CP, TeraView, UK), double the sensitivity estimated for optical and ultrasound pachymetry.

In an attempt to improve the resolution of the impulse response function and allow us to resolve the lower reflection from thinner types of the stratum corneum, our group has investigated a number of alternative processing techniques. These are wavelet denoising and sparse representation.

3.2.

Wavelet Deconvolution to Remove Noise

Deconvolution of the system response in the time domain is usually achieved by division of a sample and reference pulse in the frequency domain. Given a sample pulse g, reference pulse h, and baseline b, the impulse response of the sample f can be calculated using the following equation:

Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\rm FFT}(f) = \frac{{{\rm FFT}({g - b})}}{{{\rm FFT}({h - b})}} \Rightarrow f = {\rm FFT}^{ - 1} \left[ {\frac{{{\rm FFT}({g - b})}}{{{\rm FFT}({h - b})}}} \right], \end{equation}\end{document} $$\mathrm{FFT}\left(f\right)=\frac{\mathrm{FFT}\left(g-b\right)}{\mathrm{FFT}\left(h-b\right)}\Rightarrow f={\mathrm{FFT}}^{-1}\left[\frac{\mathrm{FFT}\left(g-b\right)}{\mathrm{FFT}\left(h-b\right)}\right],$$

where FFT and FFT⁻¹ denote the fast Fourier transform and its inverse, respectively. To suppress the amplified noise effects that Eq. 1 can cause, a bandpass filter can be used. Usually, including in the study reported in Sec. 1.1, a double Gaussian filter has been used:

Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\rm filter} = \frac{1}{{{\rm HF}}}\exp \left({ - \frac{{t^2 }}{{{\rm HF}^2 }}} \right) - \frac{1}{{{\rm LF}}}\exp \left({ - \frac{{t^2 }}{{{\rm LF}^2 }}} \right), \end{equation}\end{document} $$\mathrm{filter}=\frac{1}{\mathrm{HF}}\exp \left(-\frac{t^2}{{\mathrm{HF}}^2}\right)-\frac{1}{\mathrm{LF}}\exp \left(-\frac{t^2}{{\mathrm{LF}}^2}\right),$$

where t represents the time axis with zero in the middle, and HF and LF relate to the high and low frequency cut-offs, respectively. Equation 2 is then incorporated into Eq. 1 as follows:

Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} f = {\rm FFT}^{ - 1} \left[ {{\rm FFT}({{\rm filter}})\frac{{{\rm FFT}({g - b})}}{{{\rm FFT}({h - b})}}} \right]. \end{equation}\end{document} $$f={\mathrm{FFT}}^{-1}\left[\mathrm{FFT}\left(\mathrm{filter}\right)\frac{\mathrm{FFT}\left(g-b\right)}{\mathrm{FFT}\left(h-b\right)}\right].$$

Equation 3 represents the DGIF method used in Sec. 1.1. However, when considering the impulse responses recorded for samples using the terahertz probe, the decreased DR and SNR when compared to the conventional flatbed system, especially at higher frequencies, limits its ability to resolve distinct reflections, such as the quartz/skin interface and stratum corneum/epidermis interface.

An alternative approach to denoising terahertz signals has been proposed by our group that uses Wiener filtering and wavelet shrinkage, called frequency-wavelet domain deconvolution (FWDD), to increase noise rejection at higher frequencies and enable the probe to resolve additional interfaces.^{48, 49} Figure 6 shows the impulse response from a terahertz pulse recorded from a part of the palm with a thin stratum corneum thickness denoised either with the FWDD algorithm (a) or with the DGIF algorithm (b). The advantage of the wavelet transformation technique is in its superior pulse-preserving capabilities due to the similarity in the shape of the terahertz impulse response to typical wavelet basis sets. The regularization parameter β is used to modulate the effect of the initial Wiener filtering step which is represented by the following equation [a modified form of Eq. 3 above]:

Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} f_{{\rm Wiener}} = {\rm FFT}^{ - 1} \left\{ {\frac{{{\rm FFT}({g - b})}}{{{\rm FFT}({h - b})}}\left[ {\frac{{| {{\rm FFT}({h - b})} |^2 }}{{| {{\rm FFT}({h - b})}|^2 + \beta \frac{{N\sigma ^2 }}{S}}}} \right]} \right\}, \end{equation}\end{document} $$f_{\mathrm{Wiener}}={\mathrm{FFT}}^{-1}\left\{\frac{\mathrm{FFT}\left(g-b\right)}{\mathrm{FFT}\left(h-b\right)}\left[\frac{{\left|\mathrm{FFT}\left(h-b\right)\right|}^2}{{\left|\mathrm{FFT}\left(h-b\right)\right|}^2+\beta \frac{N{\sigma }^2}{S}}\right]\right\},$$

where g, h, and b are as defined previously, N is the number of points in the time domain, σ² is the estimated noise variance, and S is the mean power spectral density of the baseline subtracted measurements, defined as follows:

Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} S = \frac{{\left\| {({g - b}) - {\rm mean}({g - b})} \right\|_2^2 - N\sigma ^2 }}{{\left\| {({h - b})} \right\|_2^2 }}, \end{equation}\end{document} $$S=\frac{{‖\left(g-b\right)-\mathrm{mean}\left(g-b\right)‖}_2^2-N{\sigma }^2}{{‖\left(h-b\right)‖}_2^2},$$

where |X|₂ signifies the l ₂ norm of X. Typically a small value of β is used, between 0.001 and 0.05 in our previous work, in order to accurately preserve the pulse shape, at the expense of leaving in some noise. This noise can be dealt with by the wavelet shrinkage step performed upon the output of Eq. 4. It is this pulse preserving property that makes it possible to resolve closer interfaces, as is demonstrated in Ref. 48, to which the reader is directed for further information. Key to the denoising capabilities of the FWDD algorithm (and most wavelet denoising techniques) is the wavelet shrinkage step, more specifically the setting of the threshold value. If this is correctly estimated, and the appropriate wavelet basis is used, the low coefficient bases should contain mostly noise information, and so setting them to zero will “denoise” the data. However, a poor choice of wavelet basis or a too high (or low) threshold value will result in a poorly filtered dataset and could result in additional peaks being attributed to the data. An alternative approach is to consider that the number of reflections will be sparsely represented in the time domain, which is summarized below.

Fig. 6

(a) A typical probe impulse response from the palm of the hand when treated with the FWDD algorithm. (b) The same impulse response filtered using the DGIF with varying filter parameters.

3.3.

Sparse Deconvolution

Deconvolution is usually an ill-posed problem, especially when the impulse response of the system is frequency-selective: additional assumptions on the underlying signal, i.e., before convolution, have to be made in order to retrieve it with reasonable confidence. In the case of the reflection terahertz measurements, it is assumed that what is generating the measured pulse train is actually “sparse” (i.e., containing only a few non-negligible values): indeed, the measured terahertz waveform essentially results from the strong reflections caused by large refractive index changes both at the interfaces and within the sample. It is thus reasonable to look for a sparse signal f(t) that fits the measurements g(t).

This kind of problem has attracted a lot of attention in the last decade. A popular approach is called “Basis Pursuit” and consists in minimizing a “sparsifying” norm (typically, the l ¹ norm), under a quadratic constraint that specifies the quality of the data fitting.⁵⁰ Formally

Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} &&\mathop {\min }\limits_f \int \!\!\!{| {f(t)}|} dt,\,\,{\rm under}\,\,{\rm the}\,{\rm constraint}\,\,{\rm that} \nonumber\\ &&\quad \int {| {({h * f})(t) - g(t)} |^2 dt \le E}, \end{eqnarray}\end{document} $$\begin{array}{ccc}& & \underset{f}{min}\int \left|f\left(t\right)\right|dt,\mathrm{under}\mathrm{the}\mathrm{constraint}\mathrm{that}\hfill \\ & & \int {\left|\left(h*f\right)\left(t\right)-g\left(t\right)\right|}^{2}dt\le E,\hfill \end{array}$$

where E is some fitting error that we are prepared to accept (E is either zero for the genuine Basis Pursuit approach or the variance of the noise in the case of noisy measurements). Solving a discretized version of this problem (i.e., integrals become sums) can be done using an iterated reweighted least-squares algorithm which is quite fast in the current problem.⁵¹

A typical result obtained by applying this technique is illustrated in Fig. 7. Using the sparse deconvolution response, the time difference between the two interfaces is estimated to be 0.55 ps, which is comparable to the time difference observed for the two interfaces observed in similar palm data after processing with the FWDD technique.

Fig. 7

(a) The measured terahertz signal [g(t)] (red dotted line) and the reconstructed terahertz signal (blue solid line) obtained by convolving the sparse deconvolution result f(t) with an estimate of the system impulse response h(t) (not shown here). (b) The sparse deconvolution result f(t).

4.

Conclusion

In order to continue to drive forward new applications of pulsed terahertz imaging and spectroscopy in vivo, novel measurement geometries are required. In particular, the handheld probe was conceived as a way to conveniently acquire data of epithelial tissues in vivo. However, new data analysis techniques have been required in order to extract meaningful data from the probe measurements. This paper has detailed a number of processing steps that have been developed by our group to facilitate this extraction process. The baseline subtraction algorithm substituting a water measurement for the previously used quartz plate has allowed us to eliminate any further reflections reproducibly by ensuring good contact with the quartz plate through the use of a liquid medium. As a result, spectroscopic results that are comparable to transmission geometry results are now possible. The DGIF, FWDD, and sparse deconvolution techniques have all shown promise in removing the effects of noise and bandwidth limitations on resolving interfaces. In the future we hope to continue to improve these techniques to increase the potential application areas of terahertz spectroscopy and imaging.