Exploiting k-space/frequency duality toward real-time terahertz imaging

Imaging at terahertz frequencies has recently received considerable attention because many materials are semitransparent to THz waves. The principal challenge that impedes a widespread use of THz imaging is the slow acquisition time of a conventional point-by-point raster scan. In this work, we present a theoretical formulation and an experimental demonstration of a novel technique for fast compressionless terahertz imaging based on broadband Fourier optics. The technique exploits k-vector/frequency duality in Fourier optics that allows the use of a single-pixel detector to perform angular scans along a circular path, while the broadband spectrum is used to scan along the radial dimension in Fourier domain. The proposed compressionless image reconstruction technique (hybrid inverse transform) requires only a small number of measurements that scales linearly with an image’s linear size, thus promising real-time acquisition of high-resolution THz images. Additionally, our imaging technique handles equally well and on an equal theoretical footing amplitude contrast and phase contrast images, which makes this technique useful for many practical applications. A detailed analysis of the technique’s advantages and limitations is presented, and its place among other existing THz imaging techniques is clearly identified. © 2018 Optical Society of America under


INTRODUCTION
The terahertz frequency range (0.1-10 THz, wavelengths of 3 mm-30 μm) has received considerable attention over the years, thanks to the prospect of numerous advanced imaging applications benefiting from the fact that many materials are semitransparent to THz waves (for example, polymers, plastic packaging, paper, etc.) [1,2].Moreover, unlike X-rays, THz radiation is nonionizing, thus posing no risks to living beings.Many applications have been demonstrated in various applied fields such as security [3], biomedical [4], pharmaceutical [5], food industry [6], and art conservation [7].Despite all the interest and potential, many challenges remain that impede the widespread use of THz imaging.
One of the main limiting factors is the acquisition time of a THz image.Currently, spectral imaging is done using THz timedomain spectroscopy (THz-TDS) systems.The emitted THz pulse contains multiple frequencies and the detection is based on a time-domain sampling of the electric field, which provides direct access, through frequency Fourier transform, to the amplitude and the phase of the picosecond pulse.
There are two main types of broadband THz detectors: 1) THz photoconductive antennas (THz-PCA) [8] and 2) THz detection based on electro-optic sampling (EOS) in nonlinear crystals [9].Both techniques are highly sensitive, but the THz-PCA generally performs better at frequencies below 3 THz with a higher signal-tonoise (SNR) ratio thanks to the use of lock-in amplifiers and single pixel detectors.On the other hand, EOS has a better sensitivity at higher frequencies [9][10][11], but lacks the ability to use lock-in amplification when used together with CCD arrays.For spectral imaging, both methods share many similar challenges.
The first challenge is the slow time-domain pulse sampling typically based on a mechanical optical delay line.A potential solution to this problem is an asynchronous optical sampling that allows to forgo mechanical delay line and features repetition rates of several kHz [12].However, this method is expensive because it uses two synchronized femtosecond lasers.A cheaper solution can be the use of fast rotary optical delay lines.These generally come either in a prism [13][14][15][16] or a mirror [17][18][19][20][21] configuration.
The second challenge is the single-pixel nature of the THZ-PCA detectors.Some attempts were made recently to integrate several antennas on the same semiconductor chip [22][23][24][25].For example, in [23], the authors demonstrated a one-dimensional (1D) linear array of 15 pixels that was able to reduce the acquisition time from 9 h to 36 min.However, low THz signals were reported and were mainly attributed to focusing optics that are complicated to manufacture.Additionally, dense integration of multiple antennas on a chip is problematic due to crosstalk and interference between them.In the EOS technology, on the other hand, it is possible to replace the single-pixel photodiode by a charge-coupled device (CCD) [26].However, due to the impossibility of integrating lock-in amplifiers into the imaging setup, the loss in signal is severe, thus requiring additional data post-processing such as the dynamic subtraction technique [27] or the averaging over multiple frames.
Together, the two above-mentioned challenges are currently the major limiting factors that prevent the proliferation of THz spectral imaging to real-time imaging applications.With a single pixel detector, spectral imaging experiments are generally performed by physically moving the sample in the focal plane of focusing optics.When performing time-domain sampling using an optical delay line (based on a linear micropositioning stage), a single pixel with a spectral resolution of ∼1-3 GHz is typically acquired in 1-5 s.Therefore, even a low spatial resolution 32 × 32 spectral image of 1024 pixels takes ∼1 h to acquire.
Due to the complexity in building high sensitivity THz multipixel arrays, the general trend in the THz community is use advanced signal processing techniques such as compressive sensing (CS) to reduce the number of pixels required to reconstruct an image [28].When applied to imaging, the CS theory is based on the assumption that most objects have a sparse representation in a given basis.This concept allows the efficient reconstruction of a N × N image with less than N 2 measurements.Primarily in CS systems, a single THz frequency is employed because spectral information is typically not required.For example, using Fourier optics and a moving single-pixel detector, a THz image of 4096 pixels was reconstructed with only 500 measurements (12%) at a fixed frequency of 0.2 THz [29].However, this approach required mechanical movement of the detector.In [30], the same authors fixed the single-pixel detector and used a set of binary metal masks in the object plane.These masks formed a basis from which a 1024-pixel image was reconstructed using 300 (29%) to 600 (59%) measurements at a fixed frequency of 0.1 THz.
The possibility of reducing the number of measurements with the use of a single-pixel detector resulted in a spurt of activity in CS applied to THz imaging [31][32][33][34][35]. Recently, researchers have used an optically controlled silicon mask in the THz path [33] to create patterns similar to binary metal masks.A THz image of 1024 pixels was obtained with only 63 measurements in 2 s.In [34], the same group used an electronically controlled spatial light modulator in the THz beam to demonstrate image acquisition at 1 Hz by using 45 masks/s.In [35], THz near-field imaging using CS is demonstrated.Despite these advances, most of the prior research in CS was done on reconstructing amplitude modulated images predominantly in the form of binary metal masks with cutouts.However, in addition to amplitude information, practical applications frequently require phase sensing.For example, in quality control, a scratch on the surface of a plastic sheet will result in a weak amplitude contrast, but a strong phase contrast.
An alternative to CS that samples spatial information at a single THz frequency is to use broadband signals to encode spatial information in their spectra.In fact, using time or spectral data to encode the information of an image is the basis of many imaging techniques at other frequencies.For example, in optical coherence tomography, the depth information is encoded in the spectrum [36].Additionally, originating from radar technology, synthetic aperture imaging (SAI) allows the reconstruction of an image using the echoes of a scattered field measured in time domain.SAI was implemented in the THz range, providing submillimeter resolution.With THz-TDS, the scattered wave can be directly mapped in time domain and direct phase measurement of the electric waveform avoids the use of interference with a reference wave [37].The measurement can also be done with a CW source, where the axial dimension is obtained by sweeping the frequency [38].
Furthermore, dimensionality reduction of an image can be done using spectral information.The serial time-encoded amplified (STEAM) imaging system described in [39,40] is a powerful implementation of the space-to-time image transformation.In the STEAM system, 2D spatial information of the image is encoded into a broadband spectrum using spatial dispersers.Then, using dispersion-compensating fiber and Raman amplification, the spectrum of an optical pulse is mapped into the time domain.By using a single-pixel fast photodetector, one can then reconstruct a 2D image at the laser repetition rate.Similar to the spatial dispersers used in the STEAM system [41], in the THz and millimeter wave range, the echelon diffractive gratings have been used to encode the image in 1D [42] or 2D [43] spatial coordinates.The position in space is then obtained by scanning the frequency in the spectrum.
In this work, we present a compressionless imaging technique based on the k-space/frequency duality in Fourier optics that is capable of reconstructing amplitude and phase contrast THz images using only a small number of measurements proportional to the linear size of the object rather than its area.As noted in [44], in Fourier space a spectral frequency can be equated with a spatial frequency, which opens the possibility to substitute the sampling over a 2D k-space by sampling over a 1D k-space and frequency.Therefore, we first use a single pixel THz-PCA detector and a lock-in amplifier to record the broadband time pulses at some strategically chosen points in the Fourier space.Then, we use a hybrid inverse transform developed in our group to reconstruct both amplitude and phase images from the time traces collected in Fourier space.A detailed analysis of the amplitude and phase image resolution is then presented using experimental measurements and numerical data.Unlike the techniques based on compressive sensing theory, our method is lossless (in term of information), as both spatial and spectral data is used for image reconstruction.Based on solid mathematical formulation, our hybrid inverse transform can be equally well applied both to amplitude and phase imaging.

A. Hybrid Image Reconstruction Algorithm Using a Generalized Fourier Optics Approach
In what follows, we first introduce a hybrid image reconstruction algorithm based on a generalization of a Fourier optics approach to broadband pulses.The Fourier optics theory states that, for a fixed frequency ν, a field profile Sx; y; ν generated by an object that is placed in the front focal plane (object plane) of a convex lens is Fourier-transformed according to [45]  (2) The spatial frequencies, also known as components of the k-space, are related to the ξ; η coordinates as When using broadband pulses, a raster scanning in the Fourier plane [Fig.1(b)] results in the acquisition of a hyperspectral cube where, for each frequency, there is an image in the k-space [Fig.1(c)].However, as the k-components (3) are proportional to the frequency ν, one quickly arrives at the intriguing idea of sampling the k-space using the broadband nature of the THz-TDS pulse, rather than by a mechanical scanning of the k-space by displacing the detector.Indeed, by fixing a detector at a fixed position with coordinates ξ 0 ; η 0 and using a broadband light ν ∈ ν min ; ν max one can sample a linear segment of the k-space described by By changing the ratio η 0 ∕ξ 0 , the whole k-space can be sampled.The simplest way to change this ratio is to measure several points along a circle of fixed radius ρ 0 in the ξ; η plane [Figs.1(d) and 1(e)].Mathematically, for further consideration, it is more convenient to write the Fourier transform (1) in polar coordinates as where ρ is the position vector in the observation point in the Fourier plane, while r is the position vector in the object plane.
The inverse Fourier transform ( 2) is then written as As can be seen from the integral in Eq. ( 6), the frequency ν of the probing wave and the distance of the detector from the origin ρ in the Fourier plane are entering the formulation in a symmetrical fashion.Therefore, in principle, the integral over the spatial coordinate in the Fourier plane ρ at a fixed frequency ν 0 can be replaced by an integral over the frequency ν at a fixed radius ρ 0 .In other words, instead of fixing the frequency of the probing wave and recording the Fourier image by 2D scanning of the point detector, we can instead only scan along a single 1D circular path of radius ρ 0 , while using the broadband spectrum of the probing pulse to sample the k-space along the radial direction.In this way, the full 2D raster scan at a fixed frequency over the Fourier plane is avoided and is substituted by a time-domain scan along a single circle.
Mathematically, we, therefore, define a new hybrid inverse transform by substituting the integration over the 2D Fourier space in (6) by a hybrid integration over a 1D spatial coordinate and frequency, It is important to note that the image Sr reconstructed using Eq. ( 7) is different from the original image Sr; ν, as given by the standard inverse Fourier transform Eq. ( 6).Indeed, Sr; ν is a hyperspectral image that can be different for different frequencies ν.In contrast, Sr is a compounded image that incorporates information from all the frequencies sampled by the pulse.Therefore, one must recognize that although the definition of the hybrid inverse transform (7) is derived using physical arguments, its final form should be considered as a mathematical abstraction, and therefore can be generalized even further as where U ref ν is a certain frequency-dependent reference function responsible for the normalization of the integral ( 8).As we will see later, this normalization step is crucial for correct phase and amplitude retrieval, as well as for correct interpretation of the image given by (8).In what follows, we demonstrate that, in several important cases of amplitude and phase masks, the properly normalized hybrid transform (8) recovers the original amplitude and phase information.

B. Amplitude Masks
First, we consider the case of an object in the form of an amplitude mask.The resultant image Sr; ν is assumed to be spacefrequency separable and, in the object plane, is defined as A typical example of this imaging modality would be a binary mask represented by an opaque screen with a cutout pattern (image) that is illuminated with a pulsed light.In this case, Sr is a function that has only two values (1 or 0) depending whether it corresponds to the cutout or the opaque part of the screen.Here, Eν is the frequency dependent field of the probing pulse.Furthermore, as a reference function U ref ν in the hybrid inverse transform (7), we take the normalized THz-TDS trace as measured by a point detector in the center of the Fourier plane (ρ 0 0) with the amplitude mask still present in the object plane, In practice, we first perform a time-domain measurement in the origin of the Fourier plane, then U 0; ν is found using a frequency Fourier transform of the measured time-domain pulse.Using Eq. ( 5), one can also write (10) as In this case, the hybrid inverse transform becomes (12) Note that the reference U ref ν is chosen to contain an additional multiplicative factor jcf ∕ν that is necessary to ensure that the reconstructed image Sr; ϕ is proportional to the original one Sr; ν.In fact, from the detailed mathematical analysis presented in Supplement 1, Section 2.2, it follows that, in the case of amplitude masks, the hybrid inverse transform results in the original image normalized by a constant factor proportional to the area of the image, Sr; φ Sr; φ RR Sr; φdφrdr : As an example, in Fig. 2, we experimentally demonstrate the imaging of a metallic mask with a cutout in the form of a Canadian maple leaf [Fig.2(a)].Our intent here is to demonstrate the reconstruction of a complex image using the broadband spectrum and very few pixels in the Fourier plane.For comparison, we first perform regular imaging by recording the field distribution in the Fourier plane by raster scanning with a point detector on a 102 × 102 mm grid with 1.5 mm resolution (4624 pixels).At each pixel, a full-frequency spectrum is recorded, thus a complete hyperspectral image is acquired.As an example, in Figs.2(b) and 2(c) we present image amplitude and phase distributions in the Fourier plane at the particular frequency of 0.57 THz.The reconstruction of the original image is then performed using the standard inverse Fourier transform in Cartesian coordinates [Eq.( 2)].The result contains both amplitude and phase information for every frequency.In Fig. 2(d), we present the amplitude distribution of the reconstructed image at 0.57 THz that corresponds to the amplitude and phase distributions in the Fourier plane presented in Figs.2(b) and 2(c).
We now show that the hybrid inverse transform presented in this section [Eq.( 12)] can reconstruct a simple binary image with considerably fewer pixels.Thus, in Figs.2(e) and 2(f ), we present the recorded amplitude and phase distributions of the image in the k-space using only 180 pixels in the Fourier plane.These points are sampled on a circle of radius ρ 0 25 mm with the center at the origin.The traces are acquired at equal angular intervals.The data is then normalized by the complex spectrum recorded at the origin of the k-space [expression (11)].The k-space at Figs. 2(e) and 2(f ) is constructed by interpreting the spectral data as components of the k-space with Eq. ( 3).As expected, this k-space compares remarkably with the k-space obtained through raster-scanning [Figs.2(b) and 2(c)].Finally, the image is reconstructed numerically using the integration in polar coordinates of Eq. ( 12) [Fig.2(g)].Additionally, in Figs.2(h) and 2(i), we study the quality of the reconstructed image when using a different number of pixels on the circle (between 45 and 20).We note that even when using as little as 20 pixels, the maple leaf can still be clearly recognized.

C. Phase Masks
Second, we study the case of an object in the form of a dispersionless phase mask with no absorption.As an example, we can consider imaging scratches on the surface of a transparent material plate/film surrounded by air [see Fig. 3(a)].In this case, the optical path of the probing light in the object plane is given by Δr Δ 0 − μr L a n a L m n m − n m − n a hr; (14) where n a and n m are the frequency-independent refractive indices of the air and the material, L a and L m are the distances travelled in the air and in the material, while hr is the spatially dependent depth of the scratch.In this case, the field distribution in the object plane can be written as Sr; ν SrEν exp j2πνΔ 0 − μr∕c; (15) where Sr defines the slow amplitude variation of the aperturelimited beam in the object plane, Eν is the frequency dependent field of the probing pulse, and μr n m − n a hr is the optical path differential (a scratch, for example) across the object plane that we want to image.
In the case of phase masks, we define the reference U ref ν in the hybrid inverse transform (8) using a THz pulse recorded at the origin of the k-space ρ 0 0 and measured using a reference sample without the scratch hr 0, Note that unlike in the case of [11] used for amplitude masks, there is no division by the frequency in (16).This is an important difference with the case of amplitude masks that is detailed in the Supplement 1.The hybrid inverse transform (8) then becomes In Supplement 1, Section 2.3, we demonstrate that the imaginary part of the reconstructed image ( 17) is directly proportional to the optical path differential in the object plane, where the term Sr∕ RR drSr represents the slow-varying amplitude of the incident aperture-limited beam that can be measured independently (if desired) using the binary mask method described earlier.Moreover, we find experimentally, that our imaging approach given by the hybrid inverse transform (17) is so sensitive that it can readily detect nonuniformities in the phase distribution across the imaging plate.Such nonuniformities can be due to defects in the substrate geometry and composition, substrate misalignment, or even due to phase variation in the wavefront of the probing beam.Mathematically, this means that instead of a constant phase Δ 0 used in (15), we have to assume a somewhat nonuniform phase distribution Δ 0 r across the surface of a reference substrate.Therefore, to improve the quality of phase imaging, we find it beneficial to substract the phase image of a pure substrate from the phase image of a phase mask written on a similar substrate.The phase image of a substrate  18) and a supplementary beam amplitude measurement.
Imf S0 rg is reconstructed in exactly the same fashion as the phase image Imf Srg of a phase mask written on a similar substrate, by performing exactly the same steps and using the same reference (in both cases) as described above.
As an example, we demonstrate the imaging of the Greek letter π inscribed as a 100 μm deep engraving in a L m 1 mm thick slab of photosensitive resin (PlasCLEAR) printed using a stereolithography 3D printer (Asiga Freeform PRO2) [Figs.3(a) and 3(b)].This polymer is transparent in the THz region and has an almost constant refractive index of n m 1.654 [46].The amplitude and the phase distributions of the image in the k-space are presented in Figs.3(c) and 3(d).As in the case of the amplitude masks, these were inferred using the k-space/frequency duality in Fourier optics.In what follows, our goal is to map the height hr of the engraving.First, we compute the imaginary part of the hybrid inverse transform Imf Srg (18), which is shown in Fig. 3(e).There, the π engraving is clearly visible, but it is surrounded by a halo due to the somewhat nonuniform substrate used in the experiment.Two additional steps are required to extract the absolute value of the scratch depth hr.First, as described above, we retrieve the phase image of the somewhat nonuniform substrate without the scratch Imf S0 rg, which is shown in Fig. 3(f).Then, the phase variation across the reference substrate is removed from the image with the scratch by subtracting the two phase images Imf Srg − Imf S0 rg, and the result is shown in Fig. 3

(g).
There, we can see that the phase image has been considerably improved and only the letter π is visible.To extract the absolute value of the scratch depth, we need to further divide the phase image Imf Srg − Imf S0 rg by the slow varying amplitude of the aperture-limited beam Sr∕ RR drSr.As mentioned earlier, this value can be found by performing a hybrid inverse transform [Eq.( 12)] while treating the substrate without a scratch as an amplitude mask.Finally, the height distribution can be extracted using Eq. ( 18) and, as seen from the Fig. 3(h), the π symbol is an engraving of depth around 100 μm.

A. Image Resolution
First, we quote the resolution of a standard Fourier optics technique when using raster scanning in the Fourier plane.In this case, the minimal achievable resolution in the object plane follows the Fourier transform properties.Particularly, using the Nyquist theorem, the minimal resolution achievable is δr 0.5∕k max , where k max is the maximal spatial frequency sampled in the k-space.When using raster scanning, the resolution is therefore limited by the grid size Δξ in the Fourier plane.For example, in the x direction, with the definitions in Figs.4(a) and 4(b), the resolution is dx 0.25∕Δk ξ 0.25λF ∕Δξ, where Δξ is the grid size in the Fourier plane.Therefore, the resolution improves when using smaller wavelengths and a larger grid size.In the meantime, the grid spacing in the k-space d k x limits the field of view (maximal image size) such as Δx 1∕d k x λF ∕d ξ, where d ξ is the grid spacing in the Fourier plane.
In the hybrid inverse transform, k max is a function of the maximal radial position of the detector ρ max and the maximal THz frequency ν max of the pulse, so that k max ν max ρ max ∕cF .Therefore, the minimal achievable resolution is set by the Nyquist theorem to δr 0.5cF ∕ν max ρ max 0.5λ min F ∕ρ max , where λ min c∕ν max .
In the case of amplitude masks, as shown in Supplement 1, Section 3.1, the resolution of the hybrid inverse transform algorithm [Eq.(12)] closely follows the Nyquist theorem limit δr 0.5λ min F ∕ρ 0 , where ρ 0 is the radius of the circle along which the time-domain data is sampled.Particularly, we find that the hybrid inverse transform acts as a linear smoothing filter that simply averages the image inside a circle of radius ∼δr defined above.
In the case of phase masks, as shown in Supplement 1, Section 3.2, the resolution of the hybrid inverse transform algorithm [Eq.(17)] is somewhat more complex.First, the resolution in the object plane is limited by the Nyquist theorem.Second, a spatially dependent correction term proportional to the local optical path variation due to surface inhomogeneity (presence of the engraving) must be added to the resolution, δrr 0.5λ min hrn m − n a • F ∕ρ 0 , where hr is the local height of the scratch, while n m and n a are the refractive indices of the substrate material and surrounding air.To experimentally demonstrate the effect of both the radial position of the detector and the THz bandwidth on the image resolution, we present in Fig. 4 imaging of the paper cutout in the form of a snowflake fabricated using a commercial paper puncher [Fig.4(c)].The paper cutout can be considered as a phase mask due to the subwavelength thickness of the paper and its low absorption, therefore the hybrid inverse transform in the form of Eq. ( 17) is used.In Figs.4(d)-4(i), we summarize our measurements by presenting the reconstructed phase images as a function of increasing THz frequencies (ν max ) and radial positions (ρ 0 ).From these images, we conclude that the resolution of the system can be enhanced either by using spectrally broader THz pulses or by positioning the detector as far from the origin of the Fourier plane as possible.Indeed, by increasing the THz bandwidth, the resolution of the image clearly improves, as it can be seen qualitatively by comparing Figs.4(d)-4(f) and 4(g)-4(i).Moreover, by positioning the point detector at a larger radial position ρ 0 , the image resolution again is enhanced, as seen by comparing Figs.4(d)-4(g), 4(e)-4(h), and 4(f)-4(i).There, the small dendrites in the snowflake become resolvable when the radial position ρ 0 is increased.We note that the improvement in the image resolution stops when the radial position of the detector (ρ 0 ) is increased beyond the size of the limiting aperture of the optical system (lens size).

B. Advantages and Limitations of the Hybrid Image Algorithm
Now, we would like to comment on the concept of applying the presented imaging method in the context of industrial-strength, real-time THz imaging systems.The main advantage of our method is its ability to reduce significantly the image acquisition time, while retaining the high SNR offered by the TDS-THz setups based on photoconductive antennas, albeit with the loss of spectral data.In fact, the acquisition time of our method scales proportionally to the object linear size rather than its area, which is the case in standard 2D raster scanning setups.Moreover, as our system is based on the THz-TDS setup, the same imaging system can perform both hybrid imaging of the whole object (with the loss of the spectral data), followed by a more precise hyperspectral imaging using a slow raster scan of a smaller target area of an object.Both amplitude imaging modality and phase imaging modality are supported by the hybrid imaging setup, thus allowing efficient imaging of the objects with high intensity contrast (amplitude modality) encountered for example in objects with deep carvings or cutouts, as well as the low intensity contrast (phase modality) encountered for example in objects with small scratches or imperfections on their surface.
Additionally, we note that hybrid imaging based on the k-space/frequency duality opens a realistic approach to completely forgo mechanical scanning in high-speed 2D THz imaging systems.Indeed, we note that the main reason why the standard 2D raster scanning is slow is the mechanical scanning of the sample.Thus, with the current state of the art in micropositioning stages, a reliable line sampling (back and forth along a single line of a ∼5 cm long image) can be done with ∼1 Hz rates.Therefore, even a modestly large image containing 100 lines would take several minutes to acquire.At this point, we are not even talking about the hyperspectral imaging, but rather about imaging with a monochromatic THz source.Within our approach, we require the acquisition of full THz time traces along a single circle in the Fourier space.Currently, this acquisition is accomplished using a 1D scan with a point detector.Alternatively, mechanical scanning of a point detector along a circle in the k-space can be forgone by adopting a 4f system layout (instead of the 2f system used in our work), while using absorbing photomasks in the Fourier plane for sampling [44].
In this arrangement, a single pixel detector would be placed at a focal distance of the second lens, while dynamic absorbing masks in the Fourier plane [32,33] can be used to sample along the circle path.The fastest method, however, would be to use a circular array of THz photoconductive antennas, which have been already demonstrated by several research groups [22][23][24][25], thus resulting in a simultaneous interrogation of all the spatial points along a circular path in the k-space.
Finally, according to our methodology, the sampling in the second dimension is performed via interpretation of the frequency spectrum of a registered broadband THz pulse.Currently, a slow linear delay line is used for the THz pulse acquisition, thus limiting per pixel pulse acquisition rates to 1 Hz.Slow linear delay lines can be substituted by the much faster rotary delay lines that were recently demonstrated by several research groups [13][14][15][16][17][18][19][20][21] and that can enhance the acquisition rate to ∼10-100 Hz.In fact, mechanical optical delay lines can be replaced altogether by opto-electronic delay lines integrated into the fs lasers used in the THz-TDS systems, thus resulting in full THz spectral scanning rates of ∼100 Hz per pixel [12].In any case, even if the optical delay line can be made infinitely fast, the spectrum acquisition time in the THz-TDS systems is limited by the lock-in averaging time constant of ∼10 ms, which is a minimal time required per frequency to get a reliable intensity reading with a SNR of 20-40 dB in amplitude.Using a modest spectral resolution of ∼30 GHz, a full THz spectrum of up to 3 THz can then be acquired in ∼1 s.Alternatively, tunable THz systems based on frequency difference generation in photomixers can also be used to detect spectral information [47].While using technology similar to the photoconductive antennas, and featuring spectral acquisition rates of ∼1 Hz, such signals offer considerably higher SNRs and signal intensities.Therefore, the use of linear THz antenna arrays together with fast delay lines (specifically, opto-electronic delay lines) negates the need for any moving parts in the THz imaging systems based on photoconductive antennas.This would result in ∼1 Hz acquisition rates of medium resolution (∼100 × 100), THz images with very high SNRs, either in amplitude contrast or phase contrast modalities.These acquisition rates can be further improved either by sacrificing SNR while using faster optical delay lines, or not sacrificing SNR by using more powerful (albeit considerably more expensive) broadband THz sources such as those based on the two-color fs laser mixing and plasma filamentation (see [48], for example).

CONCLUSION
In conclusion, we have demonstrated a hybrid Fourier imagingbased reconstruction method that uses both spatial and spectral information to map the spatial frequencies in the Fourier plane.Then, using solid mathematical foundations, we develop theoretically and demonstrate numerically and experimentally the application of several algorithms for compressionless reconstruction of high-contrast amplitude images and low-contrast phase images.This work is motivated by the need for a real-time, high SNR THz imaging system capable of medium-to-high spatial resolution.We believe that the proposed hybrid spatial/spectral image acquisition modality in the Fourier plane can enable such systems with already existing technologies.
Funding.Canada Research Chairs; Canada Foundation for Innovation (CFI) (Project 34633) in Ubiquitous THz Photonics.
)where F is the lens focal distance and c is the speed of light.The Fourier transform U ξ; η; ν of the original field can be measured directly in the back focal plane of the lens (also known as the Fourier plane), where ξ; η are the Cartesian coordinates of the observation point [see the schematics in Fig.1(a)].Moreover, the original field distribution Sx; y; ν generated by the object can be reconstructed using the inverse Fourier transform, Sx; y; ν jν cF ZZ dξdη • U ξ; η; ν exp j2πν cF xξ yη :

Fig. 1 .
Fig. 1.Hybrid image reconstruction algorithm: (a) schematic of the object plane and the Fourier plane; (b) raster scanning on a 2D Cartesian grid in the Fourier plane; and (c) corresponding hyperspectral k-space cube.(d) Hybrid reconstruction algorithm with a 1D circular scan in the Fourier plane and (e) corresponding k-space inferred using spectral information.

2 .
Reconstruction of a binary amplitude image in the form of a maple leaf cutout in the metallic plate.(a) Schematic of the maple leaf cutout.Standard raster scanning: (b) amplitude and (c) phase of the k-space at a single frequency of 0.57 THz (4624 pixels).(d) Image reconstruction using the standard inverse Fourier transform (2).Image reconstruction using the hybrid inverse transform (12): (e) inferred k-space amplitude and (f) phase distribution using spectra of the THz time traces acquired at 180 pixels positioned along a circle of radius ρ 0 25 mm around the origin of the Fourier plane.(g) Image reconstructed using hybrid inverse transform (12) with 180 pixels, (h) 45 pixels, and (i) 20 pixels.

Fig. 3 .
Fig. 3. Reconstruction of a phase contrast image in the form of the shallow engraving of the Greek letter π onto a slab of transparent plastic.(a)-(b) Schematic of the sample n a 1, n m 1.654, and L m 1 mm.Image reconstruction using the hybrid inverse transform (17): (c) inferred k-space amplitude and (d) phase distribution using spectra of the THz time traces acquired at 180 pixels positioned along a circle of radius ρ 0 25 mm around the origin of the Fourier plane.(e) Reconstructed phase image of a substrate with engraving Imf Srg and (f) without the engraving Imf S0 rg.(g) Improved phase image of the engraving Imf Srg − Imf S0 rg.(h) Reconstructed depth of the engraving from Eq. (18) and a supplementary beam amplitude measurement.

Fig. 4 .
Fig. 4. Impact of the THz bandwidth and the radial position of the detector on phase image resolution.Schematics of (a) the k-space and (b) the object plane showing relations between resolutions and image sizes.(c) Photograph of the phase mask in the form of a snowflake cutout in very thin (∼100 μm) paper.Reconstructed phase images using Eq.(17) with ρ 0 25 mm and (d) ν max 0.46 THz, (e) ν max 0.66 THz, and (f) ν max 0.86 THz.Reconstruction with ρ 0 30 mm and (g) ν max 0.46 THz, (h) ν max 0.66 THz, and (i) ν max 0.86 THz.