Topography of hidden objects using THz digital holography with multi-beam interferences

We present a method for the separation of the signal scattered from an object hidden behind a THz-transparent sample in the framework of THz digital holography in reflection. It combines three images of different interference patterns to retrieve the amplitude and phase distribution of the object beam. Comparison of simulated with experimental images obtained from a metallic resolution target behind a Teflon plate demonstrates that the interference patterns can be described in the simple form of three-beam interference. Holographic reconstructions after the application of the method show a considerable improvement compared to standard reconstructions exclusively based on Fourier transform phase retrieval. © 2017 Optical Society of America OCIS codes: (040.2235) Far infrared or terahertz; (090.1995) Digital holography; (110.6795) Terahertz imaging; (260.3160) Interference. References and links 1. D. Mittleman, Sensing with Terahertz Radiation (Springer, 2013). 2. W. L. Chan, J. Deibel, and D. M. Mittleman, “Imaging with terahertz radiation,” Rep. Prog. Phys. 70(8), 1325– 1379 (2007). 3. Y.-S. Lee, Principles of Terahertz Science and Technology (Springer Science & Business Media, 2009). 4. M. Tonouchi, “Cutting-edge terahertz technology,” Nat. Photonics 1(2), 97–105 (2007). 5. A. Redo-Sanchez, N. Laman, B. Schulkin, and T. Tongue, “Review of terahertz technology readiness assessment and applications,” J. Infrared Millim. Terahertz Waves 34(9), 500–518 (2013). 6. O. Morikawa, K. Yamamoto, K. Kurihara, M. Tani, F. Kuwashima, and M. Hangyo, “Sub-terahertz imaging using time-domain signals obtained with photoconductive spiral antennas,” J. Opt. Soc. Am. B 33(9), 1940–1948 (2016). 7. S. Wang, B. Ferguson, D. Abbott, and X.-C. Zhang, “T-ray imaging and tomography,” J. Biol. Phys. 29(2-3), 247–256 (2003). 8. Y. Zhang, W. Zhou, X. Wang, Y. Cui, and W. Sun, “Terahertz digital holography,” Strain 44(5), 380–385 (2008). 9. U. Schnars and W. Jueptner, Digital Holography (Springer, 2005). 10. E. Hack and P. Zolliker, “Terahertz holography for imaging amplitude and phase objects,” Opt. Express 22(13), 16079–16086 (2014). 11. P. Zolliker and E. Hack, “THz holography in reflection using a high resolution microbolometer array,” Opt. Express 23(9), 10957–10967 (2015). 12. M. Locatelli, M. Ravaro, S. Bartalini, L. Consolino, M. S. Vitiello, R. Cicchi, F. Pavone, and P. De Natale, “Real-time terahertz digital holography with a quantum cascade laser,” Sci. Rep. 5(1), 13566 (2015). 13. Z. D. Taylor, R. S. Singh, D. B. Bennett, P. Tewari, C. P. Kealey, N. Bajwa, M. O. Culjat, A. Stojadinovic, H. Lee, J.-P. Hubschman, E. R. Brown, and W. S. Grundfest, “THz medical imaging: in vivo hydration sensing,” IEEE Trans. Terahertz Sci. Technol. 1(1), 201–219 (2011). 14. P. K. Rastogi and E. Hack, Optical Methods for Solid Mechanics: A Full-Field Approach (John Wiley & Sons, 2012). 15. G. Rajshekhar, S. S. Gorthi, and P. Rastogi, “Estimation of multiple phases from a single fringe pattern in digital holographic interferometry,” Opt. Express 20(2), 1281–1291 (2012). 16. R. Kulkarni and P. Rastogi, “Simultaneous measurement of in-plane and out-of-plane displacements using pseudo-Wigner-Hough transform,” Opt. Express 22(7), 8703–8711 (2014). 17. R. Kulkarni and P. Rastogi, “Digital holographic moiré for the direct and simultaneous estimation of strain and slope fields,” Opt. Express 22(19), 23192–23201 (2014). #290207 https://doi.org/10.1364/OE.25.011038 Journal © 2017 Received 8 Mar 2017; revised 28 Apr 2017; accepted 28 Apr 2017; published 3 May 2017 Vol. 25, No. 10 | 15 May 2017 | OPTICS EXPRESS 11038


Introduction
Imaging with terahertz (THz) radiation [1,2] has been an active field of research for the last two decades, thanks to the continuous improvement of THz sources and detectors.The unique capabilities of THz radiation when it interacts with matter [3] made it applicable to nondestructive testing, security inspection, and studies on biological samples [4,5].THz radiation is known to penetrate deep into most non-conducting dielectric materials, e.g.plastics, wood, paper and ceramics, which are usually opaque at optical frequencies.
Coherent THz imaging provides additional information about the sample, because the phase as well as the amplitude of the THz electric field is accessible, by employing either phase sensitive detectors or interferometric setups.Scanning focused THz pulses across the sample is still a valuable option, especially when using photoconductive antennas [6], but it comes at the expense of acquisition time.Wavefront time-domain techniques using THz collimated pulses obviate the need to move the beam relative to the sample [7,8].However, time-domain techniques typically require synchronization between the pump pulse interacting with the sample, and the probe pulse triggering the coherent acquisition.
Full-field, continuous wave (CW) imaging within a digital holography (DH) setup [9] is a well-established technique which preserves phase sensitivity while delivering reconstructions from a single interferogram.Holographic reconstructions allow retrieving the complex field of the beam scattered from an object, using the superposition with a reference beam [9].Being developed first at optical wavelengths, they are now commonly employed in the THz spectral region, too [10][11][12].A reflection holographic setup [11] is particularly suited for samples with considerable absorption in the THz region, for example biological specimens [13].THz imaging in reflection of objects hidden behind THz-transparent materials was demonstrated in simple cases [12,13].However, no systematic studies were reported.There are situations where the covering material, be it a medical gauze in contact with injured skin, or part of a multi-layered sample, introduces severe artifacts in the interference pattern which prevent effective holographic reconstructions of the hidden object.
In this paper, we model the interference patterns recorded with our holographic setup by a multiple-beam interference.The problem is mathematically equivalent to the simultaneous measurement of multidimensional deformation components via DH at optical wavelengths used in mechanical engineering [14].Several methods have been proposed to achieve signal separation from a multi-wave DH setup with multiple illuminating object beams and a single reference beam [15][16][17][18][19][20].Extraction of the interfering object beams in the spatial frequency domain is challenging, due to overlapping components of the signals.Rajshekhar et al. [19] solved the problem by introducing a spatial carrier in one of the object beams.In [15], an alternative approach based on a piecewise multicomponent polynomial phase formulation was provided, thus avoiding the need for a spatial carrier.Kulkarni et al. [20] proposed an iterative procedure which consists of dividing the interference field in blocks where the phase has a simple form and applying an augmented matrix based reconstruction to each block.
However, most of these methods exploit the fact that the intensity of the two object beams can be independently controlled in the experimental setup.This is not feasible in our case, since the multiple object beams come from the covering material and the hidden object, and the ratio of their intensities as well as their phases is governed by their structural and optical properties like the absorption coefficient and the refractive index at THz frequencies.Therefore, separation approaches based on phase shifting of the interfering beams fail, as well.
Multiple-beam analysis is also an issue when testing the surface flatness of transparent plates with interferometric methods.A three-beam interference is measured between a reference beam and two object beams coming from the front surface and the slightly tilted back surface of the plate.In order to unwrap the three-beam interferogram, Sunderland et al. [21,22] proposed a method applicable in a Twyman-Green or in a Fizeau setup.They subtracted a two-beam interferogram obtained by blocking the reference beam from the threebeam interferogram.This led to an intensity distribution out of which the phase distribution at the front and back surfaces can be separated through continuous wavelet and pseudo-Wigner-Hough transforms.
We present a method for signal separation which retrieves the amplitude and phase distribution of the beam scattered from a hidden object in a DH setup.Like some of the previous approaches it is based on difference measurements.However, our experimental setup allows us to record a third interferogram.After isolation of the object beam, its phase is retrieved with a standard Fourier demodulation technique.The method is not iterative and it does not require prior knowledge about the amplitudes of the interfering beams nor their interference angle.
The next section illustrates the theory of the method and provides simulations of interference patterns out of which the different signals are separated.Section 3 describes the experiment.Section 4 provides validation of the method by comparing experimental results with simulations, showing that the experimental images can be reproduced with a three-beam interference.Holographic reconstructions performed on data treated with and without the method are compared.The phase reconstruction allows us to retrieve the profile of the hidden object.A brief discussion of the relevant points and a summary conclude the paper.

Theory and simulations
We use an off-axis reflection DH setup with a 45° geometry like the one shown in Fig. 1(a).Let the coordinate system be such that the detector plane defines the xy-plane and let the wave vectors of all beams mostly lie in the yz-plane.Assume that the object we aim to image is hidden behind a covering material, which reflects a non-negligible amount of the incoming beam.Including this additional field, depicted in Fig. 1(b), the overall irradiance measured by the detector can be written as: Here, each complex field Ẽ i = E i exp(iφ i ), i ∈ {R, O, A} is described by an amplitude E i and a phase φ i , where the subscripts R, O, A denote the reference beam, the beam scattered from the object, and the contribution from the covering material, respectively.The time dependence exp(-iωt) has been discarded.We do not set any constraints on the intensity of Ẽ A compared to the one of Ẽ O .Equation ( 1) can be cast in the following form: In order to illustrate the information contained in Eq. ( 2), simulations of three-beam interference patterns were performed.A 9-spoked reflective Siemens Star (Fig. 2 The parameters used in the simulations reproduced those from a typical experiment (see section 3).The angle α between the hidden object and the covering material (Fig. 1(b)) was set to 0° without loss of generality in the model.The three complex fields were propagated through the Fresnel-Kirchhoff integral [24] and rotated to the detector plane as explained in [11].A vector sum in the phasor domain was computed to simulate the interference pattern of Eq. ( 2) at the detector position.
Figure 2(d) shows the resulting three-beam interference pattern I ROA .For sake of simplicity but without loss of generality, the wave vectors of all beams were chosen in the yzplane.
To reconstruct the object, the contribution due to the object field Ẽ O needs to be separated from the contribution Ẽ A .In common off-axis digital holography the interference pattern between the object beam and the reference beam can be demodulated in the Fourier domain with the help of the carrier frequency, which has to be large enough to separate the signal from the low-frequency disturbances [25].Interference pattern I RA and (i) its spectrum.(j) Intensity distribution of the subtraction performed to retrieve the object's complex beam (Eq.( 5)) and (k) its spectrum.f y is the spatial frequency along the y axis.All the spectra are displayed with a common scale.
In our case, Ẽ O and Ẽ A are modulated by the carrier frequency to the same spectral region, since their scattering surfaces are parallel to each other (α = 0°).Thus, they cannot be separated by simply windowing the spectrum.This is made clear by plotting the Fourier spectrum of I ROA (Fig. 2(e)), where the terms proportional to cos(φ R -φ A ) and cos(φ R -φ O ) overlap around ± f c .The characteristic rays of the diffraction pattern of the Siemens Star cannot be clearly seen due to the overlying diffraction pattern from the covering material.
As the signals cannot be separated in the spectrum, separation is achieved with difference measurements.With our experimental setup, we can record two more interference patterns, namely: ( ) ( ) I OA is the interference pattern obtained by blocking the reference beam.Similarly, I RA is the interference pattern obtained by moving the object out of the detector field of view.The intensities corresponding to Eqs. ( 3) and ( 4) are shown in Figs.2(f) and 2(h) along with their corresponding spectra in Figs.2(g) and 2(i).We notice that ( ) The subtraction of the two-beam interference patterns from the three-beam interference pattern reveals the interference term between the object and reference signals.Figures 2(j) and 2(k) show the distribution of the subtracted intensities and its spectrum, respectively.We immediately see that the diffraction pattern of the Siemens Star is much more evident in Fig. 2(k) than in Fig. 2(e).Now, the object signal can be safely extracted from the spectrum in Fig. 2(k) as long as the spectrum of E A 2 is centered around the zero peak.It should be noted, that the subtraction performed in Eq. ( 5) eliminates the speckle signal contained in Ẽ A even in cases where it is comparable to the object signal, because all the contributions are assumed to be stationary in time.In fact, the effectiveness of Eq. ( 5) to separate the object signal mainly depends on the spatial frequency content of E A 2. The larger the speckle size, the more the spectrum of E A 2 is confined to low spatial frequencies, thus making the separation more effective.E R E O and (φ R -φ O ) are obtained from Eq. ( 5) with a Fourier fringe demodulation algorithm.In an analogous way, E R E A and (φ R -φ A ) are obtained using Eq. ( 4).Subtraction of the fringe pattern 2E R E A cos(φ R -φ A ) from I RA leaves E R 2 + E A 2 .Knowing the sum and the product of E R 2 and E A 2 , the separation of E R and E A is straightforward.Finally, the object amplitude E O = (E R E O )/E R is obtained.Since the phase of the reference beam at the detector plane φ R is constant over the beam size, the phase differences (φ R -φ O ) and (φ R -φ A ) correspond to φ O and φ A , respectively.In this way, all the distributions E R , E O , E A , φ A , and φ O at the detector plane were retrieved.Alternatively, E R can be directly measured by blocking the beam before it interacts with the covering sample and the hidden object.Similarly, E A can be measured after blocking the reference beam and removing the hidden object.Although not strictly necessary, such measurements can serve as a comparison with the retrieved amplitudes.
The amplitude and the phase of the complex field Ẽ O at the detector plane are shown in Figs.3(a) and 3(b), respectively.Back-propagation to the corresponding position of the hidden object led to the amplitude and phase reconstructions shown in Figs.3(c) and 3(d), respectively.The reconstructed phase distribution is flat where the amplitude is non-zero because the simulated object is a flat object.On the background, where the amplitude is very low, the phase is uniformly distributed in the range [-π, π] rad.Therefore, for better visibility the phase map was thresholded with an amplitude-based criterion.

Experiment
The general principle discussed in the previous section was experimentally demonstrated using the off-axis DH setup of Fig. 1(a).The THz source was a CW far-infrared gas laser optically pumped by a CO 2 laser (FIRL100, Edinburgh Instruments, Livingston, Scotland).Its wavelength can be chosen from many emission lines between 60 and 500 μm, which can be considered monochromatic for our purpose.We used the emission line at λ = 96.5 μm caused by a rotational transition in CH 3 OH molecules and delivering several tens of mW.A quartz beam splitter transmits about 70% of the incident power, while a metal mirror directs the beam towards the two objects.The hidden object was a 100 μm-thick metallic 9-spoked Siemens Star fabricated by laser ablation with an inner diameter of 4 mm.It was placed about 17 mm from the detector plane and made an angle of about 45° to the detector plane, which ensured a reliable determination of the phase of the object beam through Fourier fringe demodulation (Eqs.( 4) and ( 5)).The covering sample was a 2 mm-thick Teflon (PTFE) plate, which is known to be highly transparent and weakly scattering in the low THz range [3].The distance between the two objects was estimated to 2 mm.By reflecting part of the incident beam from its front and back surface, the Teflon plate provided the contribution Ẽ A .Interferometric signals were recorded with an uncooled microbolometer array based on amorphous silicon featuring 480 × 640 pixels on a 17 μm pixel pitch (Gobi-640-GigE, Xenics, Leuven, Belgium).Although it was designed for infrared radiation between 8 and 14 μm, it was shown to work at THz wavelengths and to have better lateral resolution than currently available dedicated THz bolometers and pyroelectric detectors [26].
Because of the incoherent infrared emission from the environment at room temperature, the incident beam was modulated with a shutter.Each interference pattern is the difference of averaging 20 images with the shutter open and 20 images with the shutter closed.This procedure was automated using a LabVIEW-based software and took a couple of seconds per interference pattern corresponding to Eqs. ( 2)-(4).

Results
The experimental interference patterns described by Eqs. ( 2)-(4), i.e.I ROA , I OA and I RA , are shown in the top row of Fig. 4. The low frequency oblique modulation observed in Figs.4(a) and 4(b) is due to a non-zero angle α between the object and the covering material.The experimental interference pattern I RA (Fig. 4(c)) contains the carrier frequency due to the interference between Ẽ R and Ẽ A and interference effects along the rim caused by scattering of the reference beam from the rectangular detector housing.Our 45° geometry implies that all the experimental images contain more noise on their left-hand side, where the object was farther from the detector.
Simulations of the experiment are shown in the bottom row of Fig. 4. We found that the best parameter values reproducing the experimental images were a propagation distance of 17.7 mm; interference angles θ of θ yz = 48° and θ xz = 4° on the yz-and xz-plane, respectively; an angle of α yz ≈α xz = 7°.These values, particularly α yz and α xz , should be understood as effective or average parameters, because the Teflon plate had an undetermined curvature, and its orientation could not be simply described by two angles with respect to the surface of the Siemens Star.Taking the high frequency signal as background noise, the experimental intensity of the spatial frequency peak corresponding to the modulation frequency in the I OA pattern was obtained with a simulated amplitude ratio E A / E O ≈0.5.It must be noted that this result is consistent with our 45° setup geometry and the optical properties of our samples, namely: the Teflon plate had a thickness of 2 mm, a refractive index of 1.4 and an absorption coefficient of 2 cm −1 at 3 THz [3]; 80% of the total area of the Siemens Star was 100% reflecting; the incident radiation was linearly polarized along the x-axis (s-polarization).A considerable improvement is obtained after application of our subtraction method.The interference effects due to the Teflon plate are reduced and the structure of the Siemens Star can be clearly seen in Figs.5(d) and 5(e).
The backpropagation distance was optimized so to have sharp reconstructions.For reconstructions in Figs.5(a) and 5(b), a sample-to-detector distance of 17.7 mm was used.Although it was not changed when placing the Teflon plate, optimization led to a distance of 16.7 mm in the case of Figs.5(d) and 5(e).Furthermore, the image was displaced by roughly 0.5 mm along the horizontal axis upon insertion of the plate.A minor shift along the vertical axis can also be observed.We ascribe these effects to refraction through the plate.These differences in the reconstruction parameters are consistent with our acquisition geometry and the optical properties of the plate and the Siemens Star.
The reconstructed phase distributions are not flat because the Siemens Star was slightly bent.By unwrapping the phase along the black lines shown in Figs.5(b) and 5(e), we could provide an estimate of the curved profile, which is plotted in Figs.5(c) and 5(f) for the reconstructions without and with Teflon plate, respectively (black symbols: experimental data points; red line: parabolic fit).The two fits led to the same curvature to within a relative error smaller than 10%.A maximum variation of about 20 μm along the profile of the Siemens Star was estimated.When converting the wrapped phase to profile height, if no previous knowledge about the sample topography is available, dual wavelengths methods can be applied [27].

Conclusions
In this paper, a multiple-beam analysis of interference patterns of an object hidden behind a THz-transparent material was presented in the framework of THz DH in reflection.Simulations and experimental results confirm the applicability of a three-beam interference formulation in order to explain the features observed in the experimental images.
An algorithm for the separation of the hidden object beam was presented.The method basically differs from other signal separation algorithms used to determine multidimensional deformation components with two object beams and a reference beam.The idea of difference measurements is also used in interferometric investigations of the surface flatness of transparent plates.By combining three interference patterns, we have shown that the amplitude and phase of the object beam can be retrieved, without further knowledge about the intensity of the beams or their interference angles.Holographic reconstructions of the object topography substantially improved after the application of the method, thus confirming that contributions from a scattering material like a Teflon plate were properly taken into account.A profile of the hidden object was retrieved and was shown to have the same curvature as measured without covering material.
The study confirms the potential of THz radiation as a unique tool for the characterization of structures embedded in THz-transparent materials.In its current state, our method is applicable to contact experiments, where a sample is pushed against a rigid and THztransparent covering material.In this case, the covering material can be measured without the hidden sample.Potential applications in industrial inspection would still require that the object is not encapsulated in a case, or that one sample can be extracted in order to measure the diffracting properties of the housing alone.However, the presented indications on how to interpret multi-beam interference patterns and how to extract information about the sample remain valid and interesting for several applications.

Fig. 1 .
Fig. 1.(a) Schematic of a typical off-axis digital holography setup.(b) Insight into the threebeam interference arising from scattering from the object and the covering material.See the main text for an explanation of the symbols.
(a)) was used as the hidden object.Without loss of generality, the reference beam (Fig.2(b)) was approximated by a wave with a Gaussian intensity distribution and flat phase front.The contribution Ẽ A (Fig.2(c)) is supposed to hinder the reconstruction of the hidden object.For sake of illustration, a set of Gaussian beamlets with random phases on top of a larger Gaussian distribution was chosen in order to mimic a speckle pattern.Intensity and phase of the beamlets were simulated according to speckle statistics, i.e. an exponential distribution of intensities and uniform distribution of phases within [-π, π] rad[23].The spectrum of Ẽ A extends to spatial frequencies (inversely proportional to the speckle size) lower than the spatial carrier introduced by the tilted reference beam at f c ≈2 sin(θ/2)/λ.Here, θ is the interference angle between the reference beam and the object beam (Fig.1(b)), while λ is the wavelength.The amplitudes in Figs.2(a) and 2(c) were modulated by a Gaussian envelope to simulate the incident beam profile.

Fig. 2 .
Fig. 2. Simulation results.Amplitude distribution of (a) the Siemens Star, (b) the reference beam, and (c) the contribution Ẽ A .Their interference pattern is shown in (d), along with its corresponding Fourier spectrum (e).(f) Interference pattern I OA and (g) its spectrum.(h)Interference pattern I RA and (i) its spectrum.(j) Intensity distribution of the subtraction performed to retrieve the object's complex beam (Eq.(5)) and (k) its spectrum.f y is the spatial frequency along the y axis.All the spectra are displayed with a common scale.

Fig. 3 .
Fig. 3. Simulation results.(a) Retrieved amplitude and (b) phase distributions of the hidden object at the detector plane.(c) Reconstructed amplitude and (d) phase distributions after backpropagation to the position of the object.

Fig. 4 .
Fig. 4. Top row, experimental interference patterns: (a) I ROA , (b) I OA , (c) I RA .Bottom row, corresponding simulated patterns: (d) I ROA , (e) I OA , (f) I RA .Amplitude and phase reconstructions are shown in the first and second row of Fig. 5, respectively.The reconstructions shown in Figs.5(a) and 5(b) were obtained without the Teflon plate in front of the Siemens Star, with a standard Fourier transform phase retrieval algorithm.The wavy reconstruction of some of the spokes is probably due to spurious waves scattered to the detector from the surroundings of the object.The same procedure, with the sample hidden behind the Teflon plate, leads to the reconstructions shown in Figs.5(g) and 5(h).The low frequency oblique modulation is still present and it severely compromises the reconstruction of the hidden Siemens Star.A considerable improvement is obtained after application of our subtraction method.The interference effects due to the Teflon plate are reduced and the structure of the Siemens Star can be clearly seen in Figs.5(d) and 5(e).The backpropagation distance was optimized so to have sharp reconstructions.For reconstructions in Figs.5(a) and 5(b), a sample-to-detector distance of 17.7 mm was used.Although it was not changed when placing the Teflon plate, optimization led to a distance of 16.7 mm in the case of Figs.5(d) and 5(e).Furthermore, the image was displaced by roughly 0.5 mm along the horizontal axis upon insertion of the plate.A minor shift along the vertical axis can also be observed.We ascribe these effects to refraction through the plate.These differences in the reconstruction parameters are consistent with our acquisition geometry and the optical properties of the plate and the Siemens Star.The reconstructed phase distributions are not flat because the Siemens Star was slightly bent.By unwrapping the phase along the black lines shown in Figs.5(b) and 5(e), we could provide an estimate of the curved profile, which is plotted in Figs.5(c) and 5(f) for the reconstructions without and with Teflon plate, respectively (black symbols: experimental data points; red line: parabolic fit).The two fits led to the same curvature to within a relative error smaller than 10%.A maximum variation of about 20 μm along the profile of the Siemens Star was estimated.When converting the wrapped phase to profile height, if no previous

Fig. 5 .
Fig. 5. Reconstructions of the Siemens Star.Amplitude (a), (d), (g), wrapped phase (b), (e), (h), and profile (c), (f).(a), (b), (c): the Siemens Star was not hidden; (d), (e), (f): the Siemens Star was hidden behind a Teflon plate and the signal separation was applied to eliminate the contribution from the plate; (g), (h): the Siemens Star was hidden behind the Teflon plate and a standard phase retrieval algorithm was applied.As in Fig. 3, the reconstructed wrapped phase distributions were arbitrarily set to zero where the amplitude is within the noise level.In the plots (c) and (f), black symbols represent experimental data along the cross-sections highlighted in (b) and (e), while red lines show a parabolic fit.(i): colorbar for wrapped phase reconstructions.