Terahertz (THz) waves can penetrate diverse non-conducting and optically opaque materials without ionization but are nevertheless strongly absorbed by water molecules and are thus quite suitable for biomedical imaging, materials characterization, and non-destructive testing [1–4]. For samples with low amplitude contrast at the THz band, various phase imaging approaches with different types of sources and layouts have been proposed to reveal the thickness, refractive index, morphology, and internal structure from phase images. Through the utilization of some of these approaches, which use advanced image acquisition and reconstruction strategies, it is possible to achieve real-time imaging [5]. THz pulsed imaging (TPI), which is commonly based on commercial THz time-domain spectroscopy, provides broadband spectral images with both amplitude and phase information [6]. However, it has limited data acquisition speed in the raster scanning mode when using mechanical optical delay lines and THz photoconductive antennas, e.g.,  $\sim 1\mathrm{h}$ for an image of $640\times 480$ pixels [7]. This time-consuming acquisition can be addressed to some extent either by adopting asynchronous sampling with a pair of synchronized femtosecond lasers [8] or by using fast rotary prisms [9]. Full-field TPI is achieved by electro-optic sampling, while the signal-to-noise ratio (SNR) is degraded with the absence of lock-in amplifiers, and the field-of-view (FOV) is restrained by the size of non-linear crystals [10,11]. THz single-pixel imaging can effectively squeeze the long acquisition time by using mask coded illumination and compressive sensing sampling [12,13]. The amplitude and phase changes in a leaf of a terrestrial plant have been tracked in a $32\times 32$ pixels window during its dehydration process with a rate of 6 frames-per-second (fps) [14]. The acquisition time grows linearly with the expansion of the image size. With heterodyne profilometry, it is possible to scan an area of $80\mathrm{mm}\times 80\mathrm{mm}$ within 30 min and retrieve the complex amplitude distribution using two continuous-wave (CW) THz sources with adjacent frequencies [15]. Laser feedback interferometry (also known as self-mixing interferometry) uses quantum cascade lasers as transceivers for THz reflective phase imaging [16,17]. For these point scanning imaging approaches, mechanical inertia is the bottleneck to further boost the imaging speed.

By virtue of the rapid development in THz array detectors, e.g., microbolometers and pyroelectric detectors, CW THz full-field phase imaging has become the new trend in THz imaging. The phase can be retrieved from the recorded intensity of diffraction patterns and holograms, e.g., THz ptychography can simultaneously reconstruct both the incident probe distribution and the complex-valued transmission function of a sample’s partially overlapped illuminated regions using different recording geometries [18–22]. Successive recording of diffraction patterns guarantees the unique solution and convergence speed of this approach, nevertheless restraining its ability to investigate simultaneous phenomena. CW THz digital holography (TDH) is another lensless full-field phase imaging approach [23–32]. Due to its compact in-line scheme, Gabor TDH provides sub-wavelength lateral resolution [23], while the sample size is restrained for the balance of the scattered object beam and the unscattered reference beam [24]. Furthermore, the suppression of the twin image requires time-consuming iterative phase retrieval algorithms [25]. In phase-shifting TDH, the restriction of the sample size is alleviated due to the additional on-axis reference beam in the Mach–Zehnder architecture [27]. It has limited real-time imaging ability since three or more fixed-step phase-shifting holograms need to be recorded to eliminate the twin image. Common-path self-referencing TDH overcomes this issue by recording the off-axis hologram and performing a subsequent spatial Fourier filtering [28]. The FOV of this robust geometry is approximately half of size of the illumination beam. Off-axis Fresnel TDH achieves dynamic imaging in both reflection ($640\times 480$ pixels, 5 fps) [29] and transmission modes ($320\times 240$ pixels, 3.5 fps) [30]. Either Fresnel–Kirchoff or angular spectrum integral is adopted for the diffraction propagation calculation depending on the reconstruction distance. The essential numerical focusing is often time-consuming. Currently, the potential of real-time phase imaging in TDH has not yet been fully released.

Lensless Fourier-transform holography (LFH) is another prosperous real-time phase imaging approach that has been validated at other bands, e.g., the soft X-ray, extreme ultraviolet (EUV), and visible-light regime [33–38]. The object beam interferes with the diverging reference beam with the same spherical curvature. The hologram is a sinusoidal fringe of a vector spatial frequency, which makes the most efficient use of the space bandwidth product of the detector [33]. Numerical reconstruction requires a single inverse Fourier transform without additional diffraction propagation, thus providing phase information in almost real-time. To offset the quadratic phase factor, the reference beam is brought to focus at the object plane and is equidistant from the recording plane. The Mach–Zehnder configuration with the intervention of a beam splitter between the sample and the detector is widely adopted in the visible-light domain, but it is too bulky for TDH. Extended recording distance would lead to the decreased lateral resolution in the latter case. In the case of no beam splitter, e.g., the soft X-ray and EUV, an object mask with a pinhole generates a common-path divergent reference beam. The transmission energy would have considerable loss if a mask or a transmissive lens is adopted for THz imaging.

In this study, we propose an off-axis triangular LFH geometry based on oblique illumination and apply it to THz real-time full-field phase imaging. A parabolic mirror (PM) is adopted to generate a spherical reference beam. Successive reflections are suppressed by the triangular layout. The experimental lateral resolution is evaluated using a gold-plated Siemens star. The size of the reconstructed images can be expanded through sub-pixel image registration and image stitching. Finally, the dehydration process of an aquatic leaf is quantitatively monitored.

The schematic diagram of the lensless Fourier-transform TDH (LF-TDH) is shown in Fig. 1. Both the object and the point source of the reference beam are placed at the $x_0–y_0$ plane. The incident angle of the planar illumination beam is assumed to be $\theta$. The array detector is located at the recording plane, which is parallel to the object plane and denoted as the $x–y$ plane. The center of the detector aligns with the center of the illumination beam, and thus is $d\tan \theta$ away from the coordinate center at the recording plane, where $d$ denotes the vertical distance between the object plane and the recording plane. In the proposed triangular configuration, both the object wavefront along the inclined illumination direction and the reference beam can be expressed under the paraxial approximation within the detector, even if the angles $\alpha$ and $\theta$ are not very small practically.

We assume $u(x_0,y_0)$ is the complex amplitude at the object plane. Under the Fresnel approximation, the object wavefront at the recording plane can be simplified as [39] $$\begin{gathered} O(x,y)\propto \exp [\frac{jk}{2d}(x^2+y^2)]\times {\iint }_{-\infty }u(x_0,y_0)\exp \left(j2\pi \frac{\sin \theta }{\lambda }{x}_0\right) \\ \times \exp [\frac{jk}{2d}(x_0^2+y_0^2)]\exp [-\frac{jk}{d}(x{x}_0+y{y}_0)]\mathrm{d}{x}_0\mathrm{d}{y}_0 \\ \propto \exp [\frac{jk}{2d}(x^2+y^2)]U(\frac{x-d\sin \theta }{\lambda d},\frac{y}{\lambda d}), \end{gathered}$$(1)with $${U(p,q)=\Im \{u(x_0,y_0)\exp [\frac{jk}{2d}(x_0^2+y_0^2)\left]\right\}|}_{p=\frac{x}{\lambda d},q=\frac{y}{\lambda d}},$$(2)where $\lambda$ is the wavelength, $k$ represents the wave number, and $k=2\pi /\lambda$. $\Im \{\}$ denotes the Fourier transform. The divergent reference beam can be formulated as follows: $$R(x,y)\propto \exp [\frac{jk}{2d}(x^2+y^2)]\exp \left(\frac{jk}{d}x{x}_r\right),$$(3)in which ($x_r$, 0) is the coordinate of the point source of the reference beam. The reference beam interferes with the object beam with an off-axis angle of $\alpha$. The hologram can be expressed as $$I(x,y)={|O+R|}^2={|O|}^2+{|R|}^2+O^{*}R+O{R}^{*},$$(4)where * denotes the complex conjugate. In the proposed geometry, the incident angle $\theta$ is around half of the off-axis angle $\alpha$ and does not need to be accurately measured. Meanwhile, the effect of the spherical phase factor associated with Fresnel diffraction is eliminated. The fourth term of the right-hand side of Eq. (4) can be simplified as $$O(x,y)R^{*}(x,y)\propto U(\frac{x-d\sin \theta }{\lambda d},\frac{y}{\lambda d})\exp (-\frac{jk}{d}x{x}_r).$$(5)

It is noted that a linear phase shift of $\sin \theta /\lambda$ in the frequency domain is equivalent to a translation of $d\sin \theta$ in the space domain [39]. Considering the orientation of the detector in this geometry, the center of the hologram coincides with the center of the Fourier transform, which inherently offsets the phase shift from the oblique illumination. It is noted that as long as the angle $\theta$ is smaller than or equal to 15°, the value of $\sin \theta$ is close to the value of $\tan \theta$ with the relative error below 4%, which can be reasonably negligible. Let us use the coordinates ($x_1$, $y_1$), denoting the location where the detector is placed, where $x_1=x-d\sin \theta$ and $y_1=y$, so that Eq. (5) becomes $$O(x_1,y_1)R^{*}(x_1,y_1)\propto U(\frac{x_1}{\lambda d},\frac{y_1}{\lambda d})\exp (jk\sin \theta x_r)\exp (-\frac{jk}{d}{x}_1{x}_r),$$(6)where the first exponential term in Eq. (6) is a constant term that can be omitted in the phase image, while the second exponential term is a linear phase factor that is going to be accounted for when the order of $u_{+1}$ is cropped in the reconstructed image. The real image can be reconstructed using a single inverse Fourier transform ${\Im }^{-1}\{\}$ of the hologram: $$\begin{gathered} u_{+1}(x_0,y_0)={\Im }^{-1}\{O(x_1,y_1)R^{*}(x_1,y_1)\} \\ \propto u(x_0-x_r,y_0)\exp [\frac{jk}{2d}{(x_0-x_r)}^2+y_0^2]. \end{gathered}$$(7)

Meanwhile, the DC term $u_0(x_0,y_0)={\Im }^{-1}\{{|O|}^2+{|R|}^2\}$ locates at the center of the reconstructed object image, and the twin image $u_{-1}(x_0,y_0)={\Im }^{-1}\{O^{*}R\}$ is placed in ($-x_r$, 0), which is centrosymmetric to the real image. With a proper distance $x_r$ and certain off-axis angle $\alpha$, the real image is not superposed with the DC term and the twin image.

The double exposure method [40] is used to compensate for the additional quadratic phase factor in Eq. (7). The corrected phase is obtained by subtracting the phase distribution, which is reconstructed from a background hologram. This hologram is recorded using identical experimental parameters without the sample. Phase unwrapping is subsequently conducted if necessary by using the least squared method [41]. Theoretically, this method is also available for correcting phase aberration [42], which is not adopted because the reconstruction quality is not as good as the double exposure method in our reconstruction.

The diagram of the experimental setup is depicted in Fig. 2. A far-infrared ${\mathrm{CO}}_2$ laser pumped ${\mathrm{CH}}_3\mathrm{OH}$ vapor to emit a CW THz beam (295-FIRL-HP, Edinburgh Instruments) at a wavelength of 118.83 μm (2.52 THz) with a maximum output power of up to 500 mW. The diameter of the emitted beam was $\sim 7\mathrm{mm}$. It was expanded to 22 mm in diameter using two gold-coated off-axis PMs, ${\mathrm{PM}}_1$ and ${\mathrm{PM}}_2$, with effective focal lengths of 50.8 mm and 152.4 mm, respectively. The collimated beam was partitioned into two beams by an HRFZ-Si beam splitter ($n=3.4287$ at 2.5197 THz), the diameter and the thickness of which are 50.8 mm and 3 mm, respectively. The splitting ratio between the transmissive beam and the reflective beam was 54%:46%. The reflective portion went through the thin sample and generated the object beam. Another beam was reflected by ${\mathrm{PM}}_3$ with an off-axis angle of 60° and became a spherical reference beam with a focal length of 101.6 mm. The off-axis angle $\alpha$ between the object beam and reference beam was $\sim 26.5°$.

LF-TDH holograms were recorded by a pyroelectric detector (Pyrocam IV, Ophir-Spiricon, Inc.), with a pixel pitch of $80\mathrm{\mu m}\times 80\mathrm{\mu m}$ featuring $320\times 320$ pixels. Because CW THz beams must be chopped to create a changing signal, the detector contains an integrated chopper before the sensor. The chopping frequency was set to 50 Hz. To enhance the contrast of the holographic fringes, multiple frames were accumulated for each hologram. The accumulation numbers are different in the following experiments. In our imaging system, the Fresnel number was 4.48, so it could be considered that the detector was located in the Fresnel diffraction zone [39]. It is noted that the point source of the reference beam has to be placed at the object plane. An auxiliary thermal sensitive plate was placed next to the sample at the object plane blocking the path of the reference beam. Once the point source was on the plate, a tiny hole was burnt by the converging reference beam. Meanwhile, the reconstructed image would be out of focus if the point source and the sample were not on the same plane. The distance $d$ was $\sim 72\mathrm{mm}$, which was the same for the following experiments.

To determine the spatial resolution of the LF-TDH setup, a Siemens star with an outer diameter of 20 mm was adopted. A photograph of the sample is given in Fig. 3(a). The sample had a spoke angle of 7.5° and a total of 48 spokes. The opaque spokes were fabricated by sputtering an Au coating on a high-resistivity float zone Si (HRFZ-Si) substrate via photolithography. The thicknesses of the substrate and Au coating are 500 μm and 50 nm, respectively. The resolutions at the outer diameter and the center circle of this target are 0.38 and 7.66 lp/mm, respectively. The illuminated region of the sample is inside the white dashed-line box in Fig. 3(a). For the reconstruction, 500 frames were recorded and accumulated via Gaussian fitting to be a synthesized hologram [25]. Figure 3(b) shows the accumulated hologram. The fringe contrast is the metric to evaluate the quality of the hologram [43], the value of which is 0.682, as shown in the red rectangle in Fig. 3(b).

Figure 3(c1) shows the reconstructed amplitude image of the object after an inverse Fourier transform of the hologram. For better display, it was resized from the original $320\times 320$ pixels into $1200\times 1200$ pixels by the nearest interpolation method. This approach was adopted for the other two samples. It can be seen that the zero-order image occupies too much energy, which leads to blurred details in the real image. The Laplacian filter was applied to the hologram to filter out the DC term [44], and the result is shown in Fig. 3(c2). It is noted that the Tukey window function was multiplied by the filtered hologram to suppress the negative influence of the aperture diffraction [45]. The cropped region illustrated as the red-line box in Fig. 3(c2) is $220\times 220$ pixels.

In this experiment, the gold-coated sample has high reflectivity and fine features. To further suppress the minor residual successive reflections, it was not perfectly parallel to the detector. Therefore, the distortion known as anamorphism should be considered. It is shown in Fig. 3(d1) that the aspect ratio of the length to the width in the reconstructed Siemens star was 1.43. The anamorphism was numerically corrected by adopting a tilted plane correction method [21]. The results can be seen in Fig. 3(d2). Theoretically, the lateral resolution is anisotropic due to the rectangular detector. The highest resolution was $\delta =\lambda d/\sqrt{2}M\mathrm{\Delta }x=236\mathrm{\mu m}$ in the diagonal direction [39], where $M$ and $\mathrm{\Delta }x$ are the pixel number and pitch, respectively. It is illustrated in the enlarged sub-image of Fig. 3(d2) that the experimental resolution was $\sim 346\mathrm{\mu m}$ in the diagonal direction.

It is illustrated in Fig. 4(a) that a polypropylene (PP) plate patterned with letters was adopted as the sample. The length and width of the patterned region are $\sim 45\mathrm{mm}$ and $\sim 5\mathrm{mm}$, respectively. The thickness of the patterned and unpatterned areas of the plate was measured to be $\sim 411\mathrm{\mu m}$ and $\sim 337\mathrm{\mu m}$ by a screw-thread micrometer, respectively. The absorption coefficient and the refractive index are $\sim 3.9{\mathrm{cm}}^{-1}$ [46] and 1.5151 at 2.5197 THz (TeraPulse-4000, TeraView), respectively. The sample was mounted parallel to the detector. The influence of the non-uniform illumination of the THz beam on the fidelity of reconstruction cannot be neglected. To suppress this effect, two holograms in the presence and absence of the object were recorded and reconstructed, separately. Each hologram was accumulated by 500 frames via Gaussian fitting [25]. Theoretically, the amplitude difference with or without the embossed letters is $\sim 1.5\%$, and thus this PP plate can be categorized as a pure phase object. However, the letters can be clearly recognized in both the amplitude and phase images. There are two possible reasons to explain the mechanism of the silhouette in the reconstructed amplitude images. First, LF-TDH is sensitive to the thickness difference of the sample, i.e., the defocused regions. Second, THz waves have stronger diffraction effect than the coherent visible light due to the much longer wavelength. The cut-profile in Fig. 4(b4) represents clear edges of the embossed pattern and the plate. The mean phase difference between the patterned and unpatterned areas was $\sim 2.47\mathrm{rad}$, and the corresponding thickness was $\sim 78\mathrm{\mu m}$ with relative error of 5.4%.

The effective FOV of this layout, restrained by the size of the illumination beam, was approximately $22\mathrm{mm}\times 22\mathrm{mm}$. It was not sufficient for the whole sample width. Therefore, the sample was mounted on a motorized two-axis translation stage (LTS300, Thorlabs). It was scanned along the $x$ axis in steps of 1.65 mm to ensure that the adjacent positions had an area overlapping ratio of $\sim 85\%$. To obtain complete distribution of the object, sub-pixel image registration [47] and image stitching [31] algorithms were adopted after the reconstruction of the individual holograms. First, the relative transverse displacement of adjacent reconstructed images was calculated between the overlap region, which is necessary for image registration and intensity correction [47]. Second, the section images were stitched to form the expanded image. After that, the weighted smoothness was utilized to reduce the stitching error. The enlarged amplitude and phase images synthesized by 20 individual images are illustrated in Figs. 4(c) and 4(d), respectively. The image width was expanded from 18.3 mm to 44.2 mm, which was two times larger than the effective FOV.

Leaves are the main organ of photosynthesis and transpiration in plants. The changes in their water contents indicate the physiological status of plants. Different from commonly used destructive methods, e.g., gravimetry, THz spectroscopy and imaging have been validated for non-invasively quantifying the water contents of leaves [14,48–52]. Compared with terrestrial herbaceous plants, aquatic leaves have higher water contents, which dehydrate much more rapidly after leaving the water environment. A leaf of Hottonia inflata was selected as the sample and mounted on a 3-mm-thick quartz slide, which was parallel to the detector. It is shown in Fig. 5(a1) that this aquatic leaf has 11 pinnate deep-parted alternate leaflets. The image shows the leaf veins with a submillimeter resolution, which are responsible for the delivery of water and dissolved nutrients. Frames were recorded successively at a chopping frequency of 50 Hz for 300 s. A photograph of the dehydrated leaf after data acquisition is shown in Fig. 5(a2).

Considering the rapid dehydration process of the aquatic leaf, each hologram was obtained by accumulating 30 frames without Gaussian fitting, since 30 frames are statistically insufficient to extract the expectancy of the intensity from the Gaussian statistics. Therefore, the time interval of the adjacent holograms was 0.6 s. The computational time for the complete reconstruction of each hologram, on a desktop personal computer (PC) with 3.4 GHz i5-7500, 8 GB random access memory (RAM), and R2016B MATLAB, was 0.187 s, which is $\sim 1/3$ of the time interval of the successive holograms. The calculation includes the inverse Fourier transform and other processing steps, e.g., interpolation, Laplacian filtering, apodization, phase aberration correction, and phase unwrapping. It is noted that anamorphism correction is unnecessary, since the sample was mounted parallel to the detector.

Several reconstructed amplitude and phase images were extracted at an interval of 40 s, as shown in Figs. 5(b) and 5(c). Obviously, the morphology of the leaf at the initial moment was distinctly different from that in the final moment. At $t=81\mathrm{s}$, there was no sign of a droplet attached to the leaf. The contrast of the amplitude image was primarily generated from the transmittance difference due to the water content in the main and minor veins. At $t=161\mathrm{s}$, it is clearly shown in Fig. 5(b) that the absorption of THz radiation by the tips and margins of leaflets was significantly reduced, indicating that the leaf was in a dehydrated state. At $t=281\mathrm{s}$, the amplitude image shows that the tips and margins of the leaflets had higher THz transmission than the leaf main vein. In contrast, the phase distribution at that moment was clearer than at the other moments. It is known that the unwrapped phase reflects the morphology of the leaf. During dehydration, the thickness of the leaf became smaller. The dynamic change process of 500 successive images with a time interval of 0.6 s can be viewed in the video (Visualization 1). The evaporation of the water surrounding the leaf can be seen in Visualization 1, which changed the focusing position and blurred the images at the initial moments. Figures 5(d1) and 5(d2) reveal that the water content of the leaf decreased distinctly, which is due to the chloroplasts contained inside the cell collapsing and reducing the leaf thickness [53]. Further analysis of the change in water content and its calculation using the Lambert–Beer attenuation law can be seen in Appendix A, where the whole leaf as well as its small parts, i.e., the leaf tip, the rachis, and the area where a leaflet connects to the rachis, is investigated and compared.

We proposed LF-TDH, in which a triangular geometry based on tilted illumination was built. Similar to other TDH approaches based on CW THz sources and array detectors, e.g., in-line TDH and off-axis Fresnel TDH, LF-TDH is a full-field phase imaging approach. However, unlike these TDH approaches, the complex calculation of diffraction propagation is replaced by a single inverse Fourier transform of the hologram so that the speed of reconstruction is obviously accelerated. The proposed layout is different from other lensless holographic geometries at other frequencies. Compared with current lensless Fourier-transform digital holographic layouts, the proposed layout not only avoids drilling holes at the object mask and wasting THz radiation, but also simplifies the interference architecture and squeezes the recording distance for attainable high resolution. In this system, both the object and recording planes are not perpendicular to the optical axis, the emergent beam from the sample was reflected to the opposite direction by the sensor of the detector, and vice versa. Therefore, successive reflections between the object and the detector known as the Fabry–Perot effect can be greatly suppressed. If the object is not parallel to the detector, e.g., the Siemens star, a tilted plane correction method is needed with a frequency spectrum coordinate transformation.

The theoretical and experimental lateral resolutions of the proposed geometry are 236 μm and 346 μm, respectively. The pixel pitch of the THz array detector, i.e., microbolometers and pyroelectric detector, is smaller than the wavelength of the THz beam. Therefore, the aperture of the detector rather than the pixel pitch is the bottleneck of the lateral resolution of LF-TDH. The use of synthetic aperture is worthy of being investigated to enhance the resolution. Sub-pixel registration and the stitching of overlapped reconstructed images is an effective method to expand the image size. The time interval of temporal adjacent reconstructed images was 0.6 s, while the data acquisition speed, i.e., the chopping frequency of the detector, was 50 Hz. Currently, the bottleneck for the time resolution is the low SNR of the single-frame hologram. For quantitative phase imaging, 30 frames were the minimum accumulation number to get a synthesized THz digital hologram with good fringe contrast. The speed can be easily boosted by replacing the pyroelectric detector with a microbolometer. Both the radiation power fluctuation of the THz laser and the uneven distribution of the THz illumination beam had a negative influence on the fidelity of the individual reconstruction images and the synthesized ones. An appropriate THz spatial filter will help to produce a uniform Gaussian beam.

In summary, LF-TDH has a compact and robust lensless recording layout and simplified reconstruction procedure. Full-field amplitude and phase images can be reconstructed simultaneously. Although the study of dehydration of terrestrial plant leaves is not new at THz frequencies, dynamic detection of an aquatic leaf undergoing rapid dehydration is realized at the THz band for the first time, to the best of our knowledge. The proposed approach is capable of measuring the water variations on the order of sub-seconds in different parts of the leaf from successive amplitude images. When the leaf becomes more transparent to the THz radiation during its dehydration, the changes of the morphology are still clearly revealed by a series of reconstructed phase images. This method provides unprecedented insights into the water dynamics of plants and renders wide applications in THz biomedical imaging. It is noted that the natural evaporation process of the leaf moisture cells may be accelerated by the strong CW THz radiation, which was $\sim 0.284{\text{mW/}\mathrm{mm}}^2$. With the development of high-power THz sources, the thermal effect of THz waves should be considered for biomedical imaging. Our future work will involve the calculation of the absolute water content values of the leaf as well by establishing an acropetal mass transport model for water.