## Abstract

We demonstrate the enhancement of resolution and image quality in terahertz (THz) lens-free in-line digital holography by sub-pixel sampling with double-distance reconstruction. Multiple sub-pixel shifted low-resolution (LR) holograms recorded by a pyroelectric array detector (100 μm × 100 μm pixel pitch, 124 × 124 pixels) are aligned precisely to synthesize a high-resolution (HR) hologram. By this method, the lateral resolution is no more limited by the pixel pitch, and lateral resolution of 150 μm is obtained, which corresponds to 1.26λ with respect to the illuminating wavelength of 118.8 μm (2.52 THz). Compared with other published works, to date, this is the highest resolution in THz digital holography when considering the illuminating wavelength. In addition, to suppress the twin-image and zero-order artifacts, the complex amplitude distributions of both object and illuminaing background wave fields are reconstructed simultaneously. This is achieved by iterative phase retrieval between the double HR holograms and background images at two recording planes, which does not require any constraints on object plane or a priori knowledge of the sample.

© 2016 Optical Society of America

## 1. Introduction

Due to the properties of terahertz (THz) radiations, i.e., non-ionizing and ability to penetrate through non-conducting materials, THz imaging serves as a cutting edge and non-destructive technology in many fields such as biomedical diagnosis, security, quality control, etc [1–3]. Among various THz imaging techniques, including THz near-field imaging [4], THz focal plane imaging with raster scanning [5] or array detectors [6], THz spectroscopic imaging [7] and THz tomography [8], digital holography based on THz coherent radiation is a combination of THz technology and digital holography. During the measurement, the amplitude and phase can be retrieved simultaneously from obtained holograms. It has the advantages of lens-free, high resolution, label-free and potential of real-time imaging; thus has attracted more and more attention during the past few years. Both off-axis [9–14] and in-line [15–21] configurations with different devices have been demonstrated.

Resolution enhancement and twin-image elimination are two main concerns in THz digital holography. According to the Abbe criterion, resolution of a holography imaging system is limited by its numerical aperture (NA) and the wavelength:

*R*, λ, NA,

_{NA}*α*,

_{max}*z*,

*N*,

*∆s*are the best resolution, wavelength, numerical aperture, maximum diffraction angle of the object to the detection area, recording distance between sample and detector, number of pixels, pixel pitch, respectively. It is obvious from Eq. (1) that to improve the resolution of the imaging system, one should increase the value of NA. A common technique is based on synthetic aperture. By multiple overlapping sub-holograms, the lateral resolution reached 200 μm (1.68λ) in reflection for the illuminating wavelength of 118.8 μm (2.52 THz) [13] and 125 μm (1.29λ) in transmission for the illuminating wavelength of 97 μm (3.09 THz) [20]. Another method called self-extrapolation was also proposed to improve the resolution of holography in an in-line scheme by retrieving the wavefront of isolated object larger than the detector [22], and was utilized in THz in-line holography [17,19]. Two methods mentioned above are essentially similar in terms of increasing NA by enlarging the size of detection area. However, when trying to increase the value of NA there is something else we should keep in mind. According to Nyquist sampling theorem, the maximum spatial frequency should not exceed the cut-off frequency of detector, which for in-line scheme with plane waves is written as:

Another concern of in-line holography is the twin-image problem, which is caused by the lack of phase information in the hologram [25]. Therefore, the elimination or reduction of the twin-image artifact can be formulated as a phase retrieval problem, for which various strategies with different constraints on object plane and recording plane have been investigated to iteratively recover the complex amplitude distribution of the sample [25–27]. For different samples, different constraints on object plane have to be applied according to the properties of the samples. In our previous work [17,19], wavefront was reconstructed using the method proposed by Latychevskaia et al [26]. In that scheme, for various samples, only a single hologram is required, and a universal constraint based on positive absorption can be applied on object plane to remove the twin-image, which greatly simplifies the procedure of reconstruction. However, this method also has some limitations as listed below:

- 1. Normalization before the reconstruction plays an important role in this scheme. For samples on substrates which are very common in biological applications, if the size of the substrate is large, so that every photon penetrates the substrate before arriving at the detector, an absorption factor must be introduced during normalization. In this case, the non-uniformity of substrates may bring some uncertainty into calculation.
- 2. Reconstruction is based on the assumption that the sample is more absorptive than the background. So the constraint of positive absorption is not valid for samples etched or hollowed on substrates.

Hagemann et al proposed a reconstruction algorithm to retrieve the complex amplitude of the sample by taking holograms at serveral recording planes based on the idea of Allen et al. [28], and applied this to X-ray holography [29]. By this method, constraint on the object plane is not required, and the complex amplitude of the exit wave field at the object plane can be reconstructed without knowing the properties of the sample in advance. In this paper, we demonstrate the first application of this technique on THz in-line holography to improve its reconstruction quality. Furthermore, we also record the background without sample at the same recording planes to reconstruct the complex amplitude of the background wave field at the object plane. By using the complex amplitudes of these two fields, the reconstruction of the object wave can be obtained without the twin-image and zero-order artifacts.

## 2. Experimental setup

Both off-axis and in-line schemes are widely used in holography. The reconstruction in off-axis scheme is free of twin-image. However, in off-axis scheme, spectral filtering is required, which cause a decrease of resolution. Compared with off-axis scheme, in the in-line scheme, the additional reference beam and spectral filtering are not required, and there are no optical elements between sample and detector. This provides the highest resolution which is limited by NA only. However, the reconstructed result suffers from the twin-image problem. Fotunately, this problem can be solved by iterative phase retrieval algorithm. For this reason, the compact Gabor in-line scheme is adopted in our work. The experimental setup is shown in Fig. 1. A CO_{2} pumped continuous-wave THz laser (FIRL 100, Edinburgh Instrument Ltd.) is used as the light source. The average power and frequency of this laser in the experiment are ~80 mW and 2.52 THz, respectively. To know the power instability of this laser, we measure the output power 60 times within 10 minutes. The standard deviation of measured power is 1.7 mW, which is approximately 2.1% of the average value. Since the whole experiment takes less than 5 minutes, we believe the influence on the experimental results caused by laser power fluctuation is rather small, if there is any. The beam is expanded and collimated via two confocal off-axis parabolic mirrors (PM1 and PM2 in Fig. 1) with focal length of 76.2 mm and 152.4 mm, respectively. The sample position is fixed during the whole experiment. A pyroelectric array detector (Pyrocam III, Ophir-Spiricon, Inc.) is used to record the holograms and backgrounds. This detector has 124 × 124 pixels with the size of 85 μm × 85 μm, and the pixel pitch is 100 μm × 100 μm. The detector is mounted on a 3D motorized translational stage with a minimum step size of 3 μm, so that we can move the detector along X and Y direction with sub-pixel amount. In order to reach the highest resolution, the distance between object plane and recording plane 1 is set to 9.1 mm, which is the minimum distance limited by the structure of the detector. For the second distance, in [30], where the similar technique is used in the visible holography, the author mentioned that “the distance between two recorded planes can be selected freely, which does not depend on the sample structure, noise, and other experimental parameters.” However, if the recording planes are too close to each other, the difference between the two holograms is so small that the accuracy of the reconstruction will be infulenced very easily by experimental errors, including the accuracy of translational stage and the noise of the detector, etc. On the other hand, we should also not set the distance to be very large. There are two main reasons: the first one is that power of THz wave decreases fast during propagation in air, and the other is that the spatial frequency conponents in two holograms will be very different. Both facts influence the accuracy of the reconstruction. According to our previous experience, the distance between two recording planes is chosen to be 1 mm. The experiment is performed following the steps listed below:

- (I) Record a LR background image without sample at recording plane 1.
- (II) Place the sample, and keep the detector at recording plane 1. Then, record multiple LR holograms when the positions of detector are (
*X*_{0}+*i∆X*,*Y*_{0}+*j∆Y*), (*i*,*j*= 0, 1, 2), where (*X*_{0},*Y*_{0}) is the original position of the detector, and*∆X*=*∆Y*= 30 μm are the sub-pixel shifts in X and Y directions. Combining the nine LR holograms taken at different positions, we obtain a HR hologram at recording plane 1. - (III) Move the detector to recording plane 2, and repeat step (II) to obtain a HR hologram at recording plane 2.
- (IV) Remove the sample and record another LR background image at recording plane 2.
- (V) The double-distance reconstruction algorithm is applied to iteratively retrieve the complex amplitude of the sample.

## 3. Synthesis of HR hologram by sub-pixel sampling

The basic principle behind the sub-pixel sampling based pixel super resolution synthesis is to combine multiple images containing slightly different views of the same scene with sub-pixel displacement in relation to each other, and use them to construct a high-resolved image with higher sampling rate. It is equivalent to capturing the same scene using a detector with smaller pixels. This technique, referring to a signal processing towards resolution enhancement, is different from the imaging technique overcoming the diffraction limit. There are two challenges to this restoration problem. The first is that, the shifts between the LR images must be estimated precisely which is called image registration, and the second one is the synthesis of HR image from LR images. Over the past two decades, various solutions, based on frequency domain or spatial domain, have been proposed for different applications, such as satellite and aerial imaging, medical image processing, ultrasound imaging and infrared imaging. For a comprehensive overview of the image registration and super resolution synthesis, refer to [23]. The major advantage of this approach applied to THz holography is that sampling, namely pixel pitch here, is no longer a limiting factor, so that the recording distance can be dramatically reduced and imaging with large NA can be achieved, even imaging on chip. Both resolution and the signal-to-noise ratio of the reconstruction will be greatly improved, particularly when the features of the sample are comparable to or smaller than the pixel pitch. Moreover, it's very significant that relatively inexpensive LR array detector can be used to obtain competitive results compared with that by using HR camera.

In our implementation, the sub-pixel registration is based on the cross-correlation peak between LR holograms by use of matrix-multiply discrete Fourier transforms, which achieves the same accuracy as traditional FFT upsampling but with greatly reduced computational time and memory requirements [31]. Since the sample and the detector are not completely parallel in the experiment, the unwanted horizontal interference fringes are very obvious in the hologram shown in Fig. 2(a). The synthesized result is greatly degraded due to the large registration error based on the raw LR holograms. Therefore, a frequency band-pass filtering, which is critical to suppress the influence of the noisy interference fringes, is interactively applied to raw holograms before registration. In addition, inspired by the work of Vandewalle et al. [32], in which a set of aliased images are precisely aligned based on their low-frequency, aliasing-free part, we select the noisy-fringe free area as the registration part in a similar way but in the spatial domain. The pre-processing result on the 150 μm resolution target is illustrated in Fig. 2. The unwanted fringes are greatly suppressed and the red rectangle part in Fig. 2(c) is served as the registration area.

When the noisy-fringe suppressed LR holograms are accurately registered, different LR holograms are transformed into the same coordinate frame of the reference hologram which is usually the first one, via the following equation:

*X*,

_{s_HR}*Y*) are coordinates of points in HR hologram, (

_{s_HR}*X*,

_{s_LR}*Y*) are coordinates in LR hologram,

_{s_LR}*∆X*and

*∆Y*are sub-pixel shifts in column and row respectively and

*F*is the upsampling factor. Then the regular HR grid hologram is synthesized from these irregular sampled data, following the adaptive normalized convolution [33, 34], in which the local neighborhood is approximated by projecting onto a set of basis functions. The polynomial basis functions in our implementation are {

_{us}**1**,

**x**,

**y**,

**x**

^{2},

**xy**,

**y**

^{2}}. Therefore, for a point within a local neighborhood centered at

*s*

_{0}= (

*x*

_{0},

*y*

_{0}), the intensity value at

*s*= (

*x*+

*x*

_{0},

*y*+

*y*

_{0}) can be approximated by a polynomial expansion:

*x*,

*y*) are the local coordinates with respect to

*s*

_{0}and

**P**= [

*p*

_{1}

*p*

_{2}…

*p*

_{6}]

^{T}are the projection coefficients at

*s*

_{0}. The basis functions are accompanied by a so called applicability function ($a$ in Eq. (5)) which is also centered at

*s*

_{0}and represents the appropriate weight to put on the points. Similarly, each input signal has a certainty value (

**c**in Eq. (5)), which has value between 0 and 1 according to its reliability. For example, the certainty value of dead pixels should be zero. The estimation error between the approximation

*f*’(

*s*) and the measurement

*f*(

*s*) of all the

*N*points within a neighourhood

*t*is:By minimizing$\epsilon ({s}_{0})$, the projection coefficients

**P**at an output position

*s*

_{0}can be solved to be:where

**f**is an

*N*× 1 matrix of input intensity

*f*(

*s*),

**B**= [

*b*

_{1}

*b*

_{2}…

*b*

_{6}] is an

*N*× 6 matrix of six basis functions sampled at local coordinates of

*N*points, and

**W**= diag(

**c**) .* diag($a$) is an

*N*×

*N*diagonal matrix from an element-by-element product of the

**c**and $a$. Therefore, the value at output position

*s*

_{0}can be estimated by

*p*

_{1}with local coordinates (

*x*,

*y*) = (0, 0). In case of regularly sampled data, applicability function and basis functions are fixed, the implemention is very efficiently by convolution. For the irregularly sampled data in our situation, the computation of

**P**is required for each output position. The applicability function $a$ is according to [33]:where

*r*is the distance from the center of neighborhood. In our work,

*α*= 2,

*β*= 0, and the radius of local neighborhood

*r*is set to one LR pixel to cover sufficient points. The certainty

_{max}**c**for all the sub-pixel measurement are set to 1 in this paper. If the certainty is estimated before a local polynomial expansion, the influence of outliers caused by misregistration may be alleviated [34].

## 4. Double-distance iterative reconstruction

To remove the twin-image and zero-order artifacts, we reconstruct both the exit wave field *I·O* and background field *I* at the object plane. The wavefront of object *O* at the object plane is then obtained by dividing *I·O* by *I*. In this section, we describe how to extract these information from the HR holograms and LR background images taken at recording planes 1 and 2.

We define the complex amplitude of the exit wave and background fields at recording planes 1 and 2 as *U*_{1}(*I·O*), *U*_{2}(*I·O*), *U*_{1}(*I*) and *U*_{2}(*I*), respectively. Two HR holograms *H*_{1} = |*U*_{1}(*I·O*)|^{2} and *H*_{2} = |*U*_{2}(*I·O*)|^{2} at two recording planes are synthesized from multiple sub-pixel shifted LR holograms. To retrieve the complex amplitude of exit wave field, iteration is performed. The angular spectrum propagation (ASP) method [35] is used to propagate the wave fields back and forth between two recording planes. Assuming the complex amplitude of exit wave field at recording plane *i* at the *k*th iteration is *U _{i, k}* (

*I·O*) (

*i*= 1, 2;

*k*= 0, 1 … n). The iterative details are as follows:

- (II) At the
*k*th iteration, the wave field*U’*(_{2, k}*I·O*) is obtained by forward propagating the*U*(_{1, k}*I·O*) to the recording plane 2. Then measurement constraint is enforced by replacing the amplitude with the square root of the HR hologram*H*_{2}while the phase values remain unchanged. Then the*U*(_{2, k}*I·O*) iswhere φ(

*U*) represents the phase values of the complex field*U*. - (III) Backward propagate
*U*(_{2, k}*I·O*) to the recording plane 1 to get*U’*(_{1, k}*I·O*). Again the measurement constraint is applied and the*U*(_{1, k + 1}*I·O*) isThe updated complex amplitude

*U*(_{1, k + 1}*I·O*) will serve as the input for the (*k*+ 1)th iteration. - (IV) After
*n*th iteration by repeating step (II) and (III), the results converge, and we take*U*(_{i, n}*I·O*) (*i*= 1, 2) as the value of*U*(_{i}*I·O*).

Similarly, the iterative algorithm within the same framework is also used to retrieve the complex field of the illuminating background. Before the iteration, the size of backround image is changed to be the same as the HR hologram. This is carried out by a simple linear interpolation, since the background is smooth, and does not contain many high frequency components. Typically, an acceptable result can be achieved by 10-50 iterations, and after which, the complex amplitude obtained at the recording plane 1 will be back propagated to the object plane to get the wavefront of exit wave field *I·O* and illuminating background field *I*. The final complex amplitude of object *O* is reconstructed straight forwardly by dividing *I·O* by *I*. The complete process of the double-distance reconstruction is illustrated in Fig. 3. Although the multiple-plane phase retrieval algorithm has been reported before in X-ray imaging [33, 34], it has not been applied in THz region. This work demonstrates its first application in THz holography and it is also the first time that both object *O* and illuminating background field *I* are reconstructed simultaneously in THz holography. In our case, the result reconstructed by only two holograms is already free of twin-image. For this reason, we do not further increase the number of recording planes.

## 5. Results and discussion

The first sample used in the experiment is shown in Fig. 4. It is a fabricated resolution chart of 150 μm, which consists of three 150 μm wide stripes etched on a quartz slide by mechanical processing, with the interval of 150 μm between adjacent lines.

Figure 5 shows the comparison between the LR and HR holograms of the resolution chart. Figure 5(a) is one of the nine LR holograms taken at recording plane 1, which contains 124 × 124 pixels with the pixel pitch of 100 μm, while Fig. 5(c) is the HR hologram synthesized from the LR ones. Since the upsampling factor is 3, Fig. 5(c) contains 372 × 372 pixels with the equivalent pixel pitch of 33 μm. To clearly see the difference between the LR and HR holograms, the blow-ups of the areas in red dash boxes in Figs. 5(a) and 5(c) are shown in Figs. 5(b) and 5(d), respectively. From Fig. 5, it is obvious that the HR hologram shows the details much better than the LR one does, especially in the high frequency region.

With the LR and HR holograms, we reconstruct the complex amplitude of the sample using conventional ASP method by back propagating the holograms to the object plane. The amplitude and phase-contrast results obtained from a single LR hologram are shown in Figs. 6(a) and 6(c). The blow-ups of the red dash boxes in Figs. 6(a) and 6(c) are shown in Figs. 6(b) and 6(d), where the blurred strips cannot be clearly distinguished, while the features reconstructed from the HR hologram in Figs. 6(f) and 6(h) can be resolved much more easily with a sharper edge and a higher contrast. According to the Eq. (1), with the shortest recording distance of 9.1 mm we can make in the experiment, the theoretical resolution reaches ~105 μm. But in practice, it is not easy to clearly distinguish the 150 μm resolution target. The main reason is that, the strip width of 150 μm is comparable to the pixel pitch, only occuping two or three pixels with the pixel pitch of 100 μm × 100 μm so that the reconstructed results are easily disturbed by undersampling and noise. The improved results from the HR hologram demonstrates that the proposed pixel super resolution algorithm based on sub-pixel sampling has positive effects on spatial resolution.

To better understand the sub-pixel sampling, we do a further discussion on the advantages and disadvantages between synthetic aperture and sub-pixel sampling. Both of these two techniques are used to improve the resolution, but they focus on different aspects, specifically detection area and sampling, respectively. Synthetic aperture is an effective and straightforward way to enlarge detection area and consequently increase the NA of the imaging system. But by this method we need to expand the size of illumination beam to be comparable to or even bigger than detection area. As a result, the signal-to-noise ratio of holograms will be decreased due to the reduced energy density. Furthermore, with the increasing of NA by synthetic aperture, the resolution of imaging system will eventually be limited by sampling, which is determined by the pixel pitch. The most obvious advantage of sub-pixel sampling is that, it can break through the limit of the pixel pitch. Therefore, the recording distance can be as short as possible so that imaging with larger NA than synthetic aperture can be achieved. However, sub-pixel sampling cannot be used to increase the field of view while synthetic aperture can.

Next, in order to investigate the effect of the double-distance phase retrieval algorithm on quality improvement of reconstruction, we do the experiment on the sample of ‘THz’ pattern with line width of 100 μm, which is fabricated in a similar manner to the 150 μm resolution target. The complex amplitudes reconstructed by the proposed double-distance reconstruction algorithm and the conventional ASP method are shown in Fig. 7. The top row shows the amplitude and the bottom row shows the phase-contrast of the reconstruction.The amplitude and phase-contrast distributions reconstructed from a single HR hologram by ASP are shown in Figs. 7(a) and 7(b), where the results obviously suffer from the out of focus twin-image and the zero-order artifacts. As a comparison, after 20 iterations of the presented double-distance reconstruction algorithm applied to the double HR holograms at two recording planes, the complex amplitude of exit wave *I·O* are given in Figs. 7(c) and 7(d). The results are improved without twin-image artifacts, but its amplitude in Fig. 7(c) is still influenced by the non-uniform illuminating background. After dividing by the complex field of illuminating background, the final reconstructions are shown in Figs. 7(g) and 7(h), where the illuminating background is removed in the amplitude reconstruction with a more uniform background and a higher contrast. Figures 7(e) and 7(f) show the reconstruction of illuminating background *I*. Typically, the result can be converged after 10-50 iterations. It is easy to extend the algorithm to multiple recording planes so that a similar result may be obtained with less iterations. In our work, the reconstructed result from two planes is already twin-image free, so only two planes are adopted for reducing the recording time. It is important to stress again that the proposed method is robust and easy to realize due to the fact that there are not any constraints on the object plane or limitations to sample.

## 6. Conclusion

We demonstrate the enhancement of both the spatial lateral resolution and reconstruction quality in THz in-line digital holography by combining sub-pixel sampling and double-distance phase retrieval algorithm. We apply the sub-pixel sampling method to bypass the limit of the pixel pitch, by which the lateral resolution is enhanced to at least 150 μm (1.26λ) corresponding to the illuminating wavelength of 118.8 μm with the PYIII detector whose pixel pitch is 100 μm × 100 μm. The experimental result is better than already reported resolution to date when taking into account the wavelength. Relatively inexpensive array detector of large pixel pitch can be used to obtain competitive results compared with that by HR camera. Before the HR hologram is synthesized, multiple LR holograms are pre-processed by frequency filtering to reduce the unwanted fringes and only the noisy-fringe free area is served as the registration part for precisely aligning the LR holograms. To suppress the twin-image and zero-order artifacts, we propose a double-distance reconstruction algorithm which iteratively recover the phase information between the two recording planes. It does not require any constraints on object plane or a priori knowledge of the sample, making the procedure robust and easy to implement. The reconstructed results demonstrate the algorithm’s effectiveness and it is the first time in THz digital holography that the wavefront of both object and illuminating background are reconstructed simultaneously. The double-distance algorithm can be easily extended to multiple planes so that less iteration is needed to reconstruct a similar results, but with more recording time. To sum up, we report on high resolution THz in-line holography by sub-pixel sampling with double-distance reconstruction. The high lateral resolution and the high quality reconstructions, combining the unique properties of THz wave, make the presented method has a great potential for security, industrial and biomedical applications.

## Funding

National Natural Science Foundation of China (NSFC) (61475011); National Key Scientific Instrument and Equipment Development Project (2011YQ13001802); CAEP THz Science and Technology Foundation (201406, MT20150903, MT20150806).

## References and links

**1. **M. Tonouchi, “Cutting-edge terahertz technology,” Nat. Photonics **1**(2), 97–105 (2007). [CrossRef]

**2. **W. L. Chan, J. Deibel, and D. M. Mittleman, “Imaging with terahertz radiation,” Rep. Prog. Phys. **70**(8), 1325–1379 (2007). [CrossRef]

**3. **X. Yin, B. W.-H. Ng, and D. Abbott, *Terahertz Imaging for Biomedical Applications: Pattern Recognition and Tomographic Reconstruction* (Springer, 2012).

**4. **S. Hunsche, M. Koch, I. Brener, and M. C. Nuss, “THz near-field imaging,” Opt. Commun. **150**(1–6), 22–26 (1998). [CrossRef]

**5. **B. B. Hu and M. C. Nuss, “Imaging with terahertz waves,” Opt. Lett. **20**(16), 1716–1718 (1995). [CrossRef] [PubMed]

**6. **A. W. M. Lee and Q. Hu, “Real-time, continuous-wave terahertz imaging by use of a microbolometer focal-plane array,” Opt. Lett. **30**(19), 2563–2565 (2005). [CrossRef] [PubMed]

**7. **Y. C. Shen, T. Lo, P. F. Taday, B. E. Cole, W. R. Tribe, and M. C. Kemp, “Detection and identification of explosives using terahertz pulsed spectroscopic imaging,” Appl. Phys. Lett. **86**(24), 241116 (2005). [CrossRef]

**8. **B. Ferguson, S. Wang, D. Gray, D. Abbot, and X.-C. Zhang, “T-ray computed tomography,” Opt. Lett. **27**(15), 1312–1314 (2002). [CrossRef] [PubMed]

**9. **R. J. Mahon, J. A. Murphy, and W. Lanigan, “Digital holography at millimetre wavelengths,” Opt. Commun. **260**(2), 469–473 (2006). [CrossRef]

**10. **M. S. Heimbeck, M. K. Kim, D. A. Gregory, and H. O. Everitt, “Terahertz digital holography using angular spectrum and dual wavelength reconstruction methods,” Opt. Express **19**(10), 9192–9200 (2011). [CrossRef] [PubMed]

**11. **S. H. Ding, Q. Li, Y. D. Li, and Q. Wang, “Continuous-wave terahertz digital holography by use of a pyroelectric array camera,” Opt. Lett. **36**(11), 1993–1995 (2011). [CrossRef] [PubMed]

**12. **E. Hack and P. Zolliker, “Terahertz holography for imaging amplitude and phase objects,” Opt. Express **22**(13), 16079–16086 (2014). [CrossRef] [PubMed]

**13. **P. Zolliker and E. Hack, “THz holography in reflection using a high resolution microbolometer array,” Opt. Express **23**(9), 10957–10967 (2015). [CrossRef] [PubMed]

**14. **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**, 13566 (2015). [CrossRef] [PubMed]

**15. **K. Xue, Q. Li, Y. D. Li, and Q. Wang, “Continuous-wave terahertz in-line digital holography,” Opt. Lett. **37**(15), 3228–3230 (2012). [CrossRef] [PubMed]

**16. **Q. Li, K. Xue, Y. D. Li, and Q. Wang, “Experimental research on terahertz Gabor inline digital holography of concealed objects,” Appl. Opt. **51**(29), 7052–7058 (2012). [CrossRef] [PubMed]

**17. **L. Rong, T. Latychevskaia, D. Wang, X. Zhou, H. Huang, Z. Li, and Y. Wang, “Terahertz in-line digital holography of dragonfly hindwing: amplitude and phase reconstruction at enhanced resolution by extrapolation,” Opt. Express **22**(14), 17236–17245 (2014). [CrossRef] [PubMed]

**18. **J. Hu, Q. Li, and S. Cui, “Research on object-plane constraints and hologram expansion in phase retrieval algorithms for continuous-wave terahertz inline digital holography reconstruction,” Appl. Opt. **53**(30), 7112–7119 (2014). [CrossRef] [PubMed]

**19. **L. Rong, T. Latychevskaia, C. Chen, D. Wang, Z. Yu, X. Zhou, Z. Li, H. Huang, Y. Wang, and Z. Zhou, “Terahertz in-line digital holography of human hepatocellular carcinoma tissue,” Sci. Rep. **5**, 8445 (2015). [CrossRef] [PubMed]

**20. **H. Huang, L. Rong, D. Wang, W. Li, Q. Deng, B. Li, Y. Wang, Z. Zhan, X. Wang, and W. Wu, “Synthetic aperture in terahertz in-line digital holography for resolution enhancement,” Appl. Opt. **55**(3), A43–A48 (2016). [CrossRef] [PubMed]

**21. **J. Q. Hu, Q. Li, and G. H. Chen, “Reconstruction of double-exposed terahertz hologram of non-isolated object,” J. Infrared Milli. Terahz. Waves **37**(4), 328–339 (2016). [CrossRef]

**22. **T. Latychevskaia and H.-W. Fink, “Resolution enhancement in digital holography by self-extrapolation of holograms,” Opt. Express **21**(6), 7726–7733 (2013). [CrossRef] [PubMed]

**23. **S. C. Park, M. K. Park, and M. G. Kang, “Super-resolution image reconstruction: a technical overview,” IEEE Signal Process. Mag. **20**(3), 21–36 (2003). [CrossRef]

**24. **W. Bishara, T. W. Su, A. F. Coskun, and A. Ozcan, “Lensfree on-chip microscopy over a wide field-of-view using pixel super-resolution,” Opt. Express **18**(11), 11181–11191 (2010). [CrossRef] [PubMed]

**25. **L. Denis, C. Fournier, T. Fournel, and C. Ducottet, “Twin-image noise reduction by phase retrieval in in-line digital holography,” Proc. SPIE **5914**, 148–161 (2005). [CrossRef]

**26. **T. Latychevskaia and H.-W. Fink, “Solution to the twin image problem in holography,” Phys. Rev. Lett. **98**(23), 233901 (2007). [CrossRef] [PubMed]

**27. **L. Rong, F. Pan, W. Xiao, Y. Li, and F. J. Wang, “Twin image elimination from two in-line holograms via phase retrieval,” Chin. Opt. Lett. **10**(6), 060902 (2012). [CrossRef]

**28. **L. Allen and M. Oxley, “Phase retrieval from series of images obtained by defocus variation,” Opt. Commun. **199**(1–4), 65–75 (2001). [CrossRef]

**29. **J. Hagemann, A. L. Robisch, D. R. Luke, C. Homann, T. Hohage, P. Cloetens, H. Suhonen, and T. Salditt, “Reconstruction of wave front and object for inline holography from a set of detection planes,” Opt. Express **22**(10), 11552–11569 (2014). [CrossRef] [PubMed]

**30. **Y. Zhang, G. Pedrini, W. Osten, and H. Tiziani, “Whole optical wave field reconstruction from double or multi in-line holograms by phase retrieval algorithm,” Opt. Express **11**(24), 3234–3241 (2003). [CrossRef] [PubMed]

**31. **M. Guizar-Sicairos, S. T. Thurman, and J. R. Fienup, “Efficient subpixel image registration algorithms,” Opt. Lett. **33**(2), 156–158 (2008). [CrossRef] [PubMed]

**32. **P. Vandewalle, S. Süsstrunk, and M. Vetterli, “A frequency domain approach to registration of aliased images with application to super-resolution,” EURASIP J. Appl. Signal Process. **2006**, 71459 (2006). [CrossRef]

**33. **H. Knutsson and C.-F. Westin, “Normalized and differential convolution,” in Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition (IEEE, 1993), pp. 515–523. [CrossRef]

**34. **T. Q. Pham, L. J. Van Vliet, and K. Schutte, “Robust fusion of irregularly sampled data using adaptive normalized convolution,” EURASIP J. Adv. Sig. Pr. **2006**, 1–12 (2006). [CrossRef]

**35. **J. W. Goodman, *Introduction to Fourier Optics* (McGraw-Hill, 1996).