## Abstract

We demonstrate a novel computational method for high resolution image recovery from a single digital hologram frame. The complex object field is obtained from the recorded hologram by solving a constrained optimization problem. This approach which is unlike the physical hologram replay process is shown to provide high quality image recovery even when the dc and the cross terms in the hologram overlap in the Fourier domain. Experimental results are shown for a Fresnel zone hologram of a resolution chart, intentionally recorded with a small off-axis reference beam angle. Excellent image recovery is observed without the presence of dc or twin image terms and with minimal speckle noise.

©2013 Optical Society of America

## Corrections

Kedar Khare, P. T. Samsheer Ali, and Joby Joseph, "Single shot high resolution digital holography: Erratum," Opt. Express**21**, 5634-5634 (2013)

https://www.osapublishing.org/oe/abstract.cfm?uri=oe-21-5-5634

## 1. Introduction

Digital Holography (DH) has been established as an important imaging modality with wide range of applications like bio-imaging, microscopy, phase contrast and quantitative phase imaging, 3D displays, nondestructive imaging techniques, etc. DH involves recording of a hologram or interference pattern between a reference beam and an object beam – both derived from the same source - on a sensor array with subsequent image reconstruction by numerical algorithms [1,2]. With the advent of sensor technology and fast computers DH has opened new possibilities for 3D imaging. While digital array detectors (CCD or CMOS) have replaced the photographic film as the preferential recording medium for DH, most of the current image reconstruction methods used in DH are still equivalent to the physical hologram reconstruction process. As a result the recovered images show artifacts such as the dc and the twin image formation when a single hologram frame is used for image recovery. Further the image resolution in off-axis DH is considered limited by the minimum reference beam angle condition [23, Section 9.4.3] which states that the minimum allowable angle between the reference beam and the nominal object beam direction is given by${\theta}_{\mathrm{min}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\mathrm{arcsin}(3B\lambda )$, with $B$being the highest spatial frequency in the object. The pixel sizes of the current CCD and CMOS arrays used in DH are typically in the range of 2-5 µm. When using visible light wavelengths, this pixel size range along with the minimum reference beam angle condition gives rise to an artificial resolution limit on DH imaging, although the digital detector array may itself be capable of higher resolution imaging. We point out that overcoming the above limitations on DH imaging requires a new a computational approach for image recovery as described here.

The hologram pattern recorded as a result of the interference of an object beam $O(x,y)$ and a reference beam $R(x,y)$ in the detector plane may be described by:

In an off-axis digital holography setup [3] the reference beam may be described as $R\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{R}_{0}\text{\hspace{0.17em}}\mathrm{exp}(i\text{\hspace{0.17em}}k\text{\hspace{0.17em}}x\text{\hspace{0.05em}}\text{\hspace{0.17em}}\mathrm{sin}\theta )$ with ${R}_{0}$being the constant amplitude of the reference plane wave and $k\text{\hspace{0.17em}}=\text{\hspace{0.17em}}2\pi \text{\hspace{0.17em}}/\text{\hspace{0.17em}}\lambda $with$\lambda $denoting the wavelength of illumination. The first two terms on the right hand side are centered at the zero frequency in 2D Fourier transform space and are referred to as the dc terms. The location of the cross terms is decided by the carrier frequency corresponding to the off-axis angle $\theta $of the reference beam. One of the important tasks in DH numerical processing is the suppression of the dc terms and the twin image term$O*R$. The remaining term $OR*$when processed appropriately gives the image of the object that is of interest.

The removal of the dc and the twin image terms has been achieved in the DH literature using several different techniques. In early digital holography work, Kreis and Juptner [4] suggested the subtraction of the average intensity from the recorded hologram for suppression of the dc term. Takaki *et al.* [5] used two or more recorded images either using blocking shutters or using phase shifters to eliminate the zero-order and twin image terms. Phase shifting methods in digital holography have also been studied by Yamaguchi and Zhang [6]. Cuche *et al.* [7] used spatial filtering approach in the frequency domain for eliminating the dc and twin image terms. A sparse sampling scheme for recovering the object function from recorded hologram was studied by Khare and George [8] which required the elimination of the dc terms by acquiring them separately and subtracting from the recorded hologram. Chen *et al.* [9] have suggested a procedure for approximately generating the dc terms from the recorded hologram and then subtracting them. Ma *et al.* [10] designed a complex FIR filter which is then used with single hologram frame for eliminating the dc and the twin image contributions. Liebling *et al.* [11,12] have proposed two algorithms- one based on non-linear Fresnelet approximation and another based on local weighted least squares - for estimation of amplitude and phase from a single hologram. Zhong *et al.* [13,14] made use of 1D &2D Gabor wavelets for the recovery of the complex object fields from an off axis digital hologram. We would like to remark that most of the dc and twin image removal techniques in above references depend either on the assumption that the dc and cross terms in the hologram are well separated in the Fourier domain or require multiple hologram frames as in phase shifting methods. As already stated above, this requirement of separation of the dc and the cross terms is often considered as a key factor governing the resolution achievable in single frame DH imaging. In this paper, we describe a new computational approach for image recovery that gives high quality image recovery from a single hologram frame even when the dc and the cross terms have substantial overlap with each other. We treat the problem of recovery of the complex object field $O(x,y)$in the hologram plane as a constrained optimization problem with a smoothness constraint on the object field. The method bears some similarity to the recent algorithm developments in the area of compressive sensing [15,16] although the problem here is somewhat different in nature.

The paper is organized as follows. In Section 2, we describe our computational approach for recovering the amplitude and phase of the complex object field from a recorded hologram. The experimental results are described in Section 3. In order to illustrate the novelty of this approach we have intentionally recorded a Fresnel zone hologram with a small reference beam angle so that the dc and the cross terms overlap in Fourier domain. The image recovery results are shown in Section 3 with comparison with a few standard Fourier filtering techniques. One simulation experiment with a synthetically computed hologram is also shown in order to illustrate the accuracy of our approach. Finally in Section 4 we provide concluding remarks.

## 2. Image recovery in DH as a constrained optimization problem

We model the problem of reconstruction of the complex object field $O$ in the hologram plane as a constrained optimization problem [17]. For a digitally recorded hologram$H$, a known reference beam $R$ and an unknown complex object function $O$ in the hologram plane, we define a cost function to be minimized as

*q*denotes some neighborhood ${N}_{p}$of the pixel denoted by

*p*. The weight function ${w}_{pq}$will in general be a decreasing function of the distance between pixels

*p*and

*q*. These types of penalty functions are often used in medical imaging literature for iterative image reconstruction [18]. A solution to the minimization problem in Eq. (2) will lead to a complex valued function that is locally smooth over the scale of the neighborhood window size ${N}_{p}$that is used to define$\psi (O,O*)$. Small values of $\alpha $correspond to a solution which is close to the L2-norm error data fit; whereas large values of $\alpha $correspond to overly smooth solutions that may have large L2-error with respect to the given data. The introduction of this smoothness constraint may be justified on the following grounds: (a) In a Fresnel zone hologram configuration, the complex object function $O(x,y)$ is a result of diffraction of light waves from the object to the hologram plane. Any sharp edges in the object are thus expected to be blurred due to diffraction. (b) The solution $O(x,y)$ is not expected to be modulated by the carrier fringes or any of their harmonics. We solve the optimization problem above iteratively by a gradient descent scheme [17]. Since the cost function in Eq. (2) is a function of both$O$and$O*$, the steepest gradient direction is computed with respect to $O*$ [19]. The gradient of the cost function in Eq. (2) with respect to $O*$ may be calculated as

^{th}iteration by ${O}^{(n)}$ the corresponding update step is given by

*t*is a positive constant representing the step size which may be selected by standard line-search methods [20]. Operationally, we have implemented the iterative method in Eq. (4) by alternating between the least square term and the smoothness term. The L2-error is reduced by using Eq. (4) with ‘$\alpha $’ set to zero. We have implemented the smoothness constraint by convolving the resultant updated solution with an averaging filter of appropriate size. The size of the averaging filter will correspond to the neighborhood window size in the definition of $\psi (O,O*)$. Sophisticated penalty functions for minimizing local differences in the complex object function

*O*may certainly be used [18,21], however, a simpler approach as used here is sufficient for demonstrating the main idea. While we do not provide any convergence proofs, such alternating iterative procedures have been used successfully in the X-ray image reconstruction literature recently [22]. The iterative solution implemented in this paper with the averaging filter denoted as $G$ is thus as follows:

*G*may thus be selected to suppress the carrier-frequency component and its second harmonic that may get introduced into the solution. We observe that using an averaging filter with a spatial extent of approximately equal to half the carrier-fringe period (corresponding to the second harmonic of the carrier fringe) is sufficient for this purpose as we illustrate in the following section.

The conventional off-axis hologram replay process involves illuminating the hologram with the conjugate reference beam followed by back-propagation to the image plane. This process is similar to the demodulation of AM/FM signals in communication problems. The various spatial filtering approaches for recovering the object field from a recorded hologram may be considered equivalent to the physical replay process. The demodulation of the hologram pattern achieved by the iterative optimization procedure described above is, however, quite different in nature. As long as the carrier-fringes are sampled sufficiently by the digital array, the object function $O(x,y)$ in the hologram plane is seen to be recoverable even in case of overlap in the dc and cross terms.

The result of the iterative procedure in Eq. (6) gives a smooth complex object function $O(x,y)$ without any modulation which may then be Fresnel back-propagated [23] to the desired plane to obtain the required image. This is done by convolving the recovered function $O(x,y)$with the back-Fresnel impulse response given by

## 3. Experiments on single shot image recovery in Digital Holography

We have tested the image recovery method described in Section 2 in a Fresnel zone hologram setup as shown in Fig. 1 . The Mach-Zehnder configuration used here is suitable for recording the hologram with a small angle between the object and the reference beams. We selected a small nominal angle (~0.6 degrees) between the reference and the object beam. For our setup this makes sure that the dc and the cross terms in the hologram overlap with each other in the Fourier domain. This case is of importance as it is then difficult to separate these terms by spatial filtering methods, and thereby deciding the achievable resolution in a DH imaging experiment. Figure 2(a) shows the Fresnel zone hologram $H(\text{\hspace{0.17em}}x,\text{\hspace{0.17em}}y\text{\hspace{0.17em}}\text{\hspace{0.05em}})$recorded on an array detector (Imaging Source, DMK 72BUC02, CMOS, 1944 x 2592 pixel array, 2.2 µm square pixels). A part of the USAF resolution chart was used as the object. The amplitude of the Fourier transform of the hologram pattern is shown in Fig. 2(b) on a logarithmic scale in order to enhance the dynamic range of the display. The three peaks in this image correspond to the dc and the cross terms as in Eq. (1). The spread of the energy around each of the peaks shows overlap of the cross terms with the dc term. In order to determine the reference beam $R\text{\hspace{0.17em}}(\text{\hspace{0.17em}}x,\text{\hspace{0.17em}}y\text{\hspace{0.17em}}\text{\hspace{0.05em}})$ we removed the object in Fig. 1 and recorded an interference of two plane waves at the nominal angle. This interference pattern can be demodulated easily using Fourier transform techniques to obtain the complex function$R\text{\hspace{0.17em}}(\text{\hspace{0.17em}}x,\text{\hspace{0.17em}}y\text{\hspace{0.17em}}\text{\hspace{0.05em}})$. Determination of the reference beam may be considered as a calibration step which needs to be performed only once. The same $R\text{\hspace{0.17em}}(\text{\hspace{0.17em}}x,\text{\hspace{0.17em}}y\text{\hspace{0.17em}}\text{\hspace{0.05em}})$ may be used for any object as long as the setup remains the same.

In what follows we will show image recovery using the single recorded hologram pattern in Fig. 2(a) with a few different approaches to clearly illustrate the advantages of the constrained optimization procedure outlined in Section 2. The knowledge of the complex valued function $R\text{\hspace{0.17em}}(\text{\hspace{0.17em}}x,\text{\hspace{0.17em}}y\text{\hspace{0.17em}}\text{\hspace{0.05em}})$ as obtained from the calibration step will be assumed.

#### 3.1 Image recovery by simulating the physical hologram replay process

In this most straightforward approach we multiply the recorded hologram pattern $H(\text{\hspace{0.17em}}x,\text{\hspace{0.17em}}y\text{\hspace{0.17em}}\text{\hspace{0.05em}})$ as in Fig. 2(a) by the conjugate reference beam and Fresnel back-propagate the resultant function to the image plane. The corresponding recovery is shown in Fig. 3 . We observe that the dc and the two cross terms produce overlapping patterns and resulting in poor image quality.

#### 3.2 Image recovery by spatial filtering approach

Spatial filtering of the recorded hologram pattern is another standard method for processing digital holograms. We have already noted that the overlapping of the terms makes it difficult to select an appropriate spatial filter. We illustrate this point using three different spatial filters with which we attempt to separate out the cross term$OR*$. In the Fourier transform image Fig. 2(b), the separation of the dc and the sideband peak is 122 pixels. First we choose a spatial filter of size 40x40 pixels centered on the carrier frequency peak. The resulting image recovery after Fresnel back propagation is shown in Fig. 4(b) . The recovery does not have much interference from the dc term but is seen to have low resolution as expected. A larger sized 122x122 spatial filter (filter window size is made equal to the carrier frequency peak distance from center) can improve the image resolution as shown in Fig. 4(d) but the recovery now has interference from the dc term causing introduction of speckle noise which obscures the high resolution features. For an even larger filter size of 244x244 pixels as in Fig. 4(e), the corresponding reconstruction in Fig. 4(f) is completely obscured by contribution from the dc term. The useful resolution has thus not improved due to the speckle artifacts and/or obscuration due to the dc term when the filter window size is increased. Although simple square shaped filters are used here, any windowing on the filters is not expected to solve the dc term interference problem in this type of approach.

#### 3.3 Image recovery by solving minimum L2-norm problem

As a first step in DH recovery by an optimization approach, we follow the iterative scheme in Eq. (6) without any smoothing operation. The iteration may be described as:

*t*was selected using line search during each iteration and was found to be in the range 0.01 to 0.001. We observe that the iteration is fully in the hologram plane domain and requires purely algebraic operations and only a few seconds are required for completing 20 iterations for a hologram frame size of 1944x2592. The computation was performed using MATLAB (our code is not necessarily optimized in any way) on a PC with 3.3 GHz Xeon processor. We observe that the recovered field is modulated by the carrier fringes and their harmonic components. When this field is Fresnel back-propagated to the image plane, the recovery as shown in Fig. 5(c), 5(d) has artifacts due to this modulation. Figure 5(d) is a magnified version of a part of Fig. 5(c). In particular multiple overlapping orders are seen in the recovery as a result of the modulation. It is however important to note that the contribution of the dc terms is seen to be removed and some of the high frequency features in the center of the resolution chart appear better resolved as compared to those in Fig. 4(d).

#### 3.4 Image recovery by solving minimum L2-norm problem with smoothness constraint

Although the iterative solution in Eq. (8) provides a minimum L2-norm solution for the complex field $O(x,y)$ from a single hologram frame, the solution is modulated by the carrier fringes and their harmonics. The modulation on the solution needs to be averaged out in the iterative process to obtain the desired solution which we show next. The iterative procedure as in Eq. (6) is now followed. The carrier fringe period in the recorded hologram (Fig. 1(a)) corresponded to 24 pixels on our CMOS array detector. For the averaging filter$G$we use a 12x12 pixel template with all ones which is normalized by overall factor of 144. Once again 20 iterations starting with zero image as the initial guess were performed. The step size *t* was selected using line search during each iteration and was again found to be in the range 0.01 to 0.001. The computation time required was approximately 10 seconds due to the additional step of convolution with the averaging filter. Once again we would like to point out that the use of specialized hardware or optimization of code could result in much better timing performance although this is not the main focus of the present work. Figures 6(a)
and 6(b) show amplitude and phase of the recovered field $O(x,y)$in the hologram plane and Figs. 6(c), 6(d) show the image recovery after back propagating this field to the object plane. Figure 6(d) is a magnified version of a part of Fig. 6(c). We observe a high quality image recovery without much corruption due to speckle noise. Further, there is no interference from the twin image or the dc terms and as a result the high frequency features in the resolution chart are much more clearly visible. The experimental results demonstrate the advantage of our constrained optimization approach for image recovery from a single digital hologram. We believe that this is the first demonstration of high quality image recovery from single hologram in DH, where the limit posed by overlapping of the dc and cross terms has been overcome. All the image recoveries in Figs. 3–6 have been shown on the same window/grayscale level.

#### 3.5 Complex object wave and image recovery for a synthetic hologram

In order to evaluate our iterative approach for accuracy we test our procedure for a synthetically computed hologram. An object consisting of an image of the letter “A” is used as in Fig. 7(a) and Fresnel diffraction computation is performed to obtain the amplitude (Fig. 7(b)) and the phase (Fig. 7(c)) of the diffracted wave. The object extent is 5 mm and the Fresnel propagation distance used is 5 cm and a wavelength of 0.5 µm is assumed. A hologram of this object is computed numerically (Fig. 7(d)) by interfering the diffracted field with a small angle (0.6 degree) off-axis plane wave as in our experiment. The iterative procedure is applied to this numerically computed hologram and the amplitude and phase of the complex object wave recovered from the hologram are shown in Figs. 7(e), 7(f) respectively. For the numerical values considered, the step size selected using line search was found to be in the range of 0.1-0.001 over 20 iterations. We observe that the recovered fields in Figs. 7(e), 7(f) appear slightly blurred compared to those in Figs. 7(b), 7(c) respectively. This is expected due to the averaging filter used in the iterative procedure. The recovered complex object field is seen to match the computed Fresnel diffraction field (Figs. 7(b), 7(c)) numerically with approximately 10% L2-norm error. The back-Fresnel propagation of this recovered object field yields excellent image recovery as shown in Fig. 7(g) thus validating our procedure for this synthetic case. An absolute difference image between the original object (Fig. 7(a)) and the recovery (Fig. 7(g)) is shown in Fig. 7(h). We observe that the only significant difference is at the edges of the object indicating slight blurring of the recovered object with respect to the original object in Fig. 7(a) that we started with.

## 4. Conclusion

Digital Holography (DH) and recording of interference patterns in general is a very important problem in optics and the ability to perform single shot high resolution DH imaging as shown here has multiple potential applications. The computational part of the DH imaging as it is performed in the current literature is equivalent to the physical hologram reconstruction as in film based holography. As a result, although the DH data is available numerically, the performance of single shot DH imaging is believed to be limited by the same factors such as the dc and twin image formation and the minimum reference beam angle condition. We believe that this is the first demonstration which suggests that the numerical processing if done in an unconventional way may overcome the above limits on DH imaging. We have presented a novel computational method where the complex object field is obtained from a single hologram frame by solving a constrained optimization problem. Experimental results show that high quality image recovery is possible even when the dc and the cross terms in a hologram overlap in the Fourier domain. While we have considered the case of off-axis holography, our initial simulation studies suggest that such procedure may be usable for more general configurations as long as the reference beam has some known spatially varying form. The recovery of both amplitude and phase of the object field in the hologram plane as achieved here suggests that this approach may provide a powerful tool for single shot quantitative phase imaging applications. The resolution limit and quantitative accuracy achievable by the constrained optimization approach as shown here is an aspect that needs further evaluation, although we have presented brief analysis using a synthetic hologram. We will report on these aspects of the present problem and explore potential applications of this technique in future.

## References and links

**1. **U. Schnars and W. P. Jüptner, *Digital Holography–Digital Hologram Recording, Numerical Reconstruction and Related Techniques* (Springer-Verlag, Berlin, 2005).

**2. **L. Yaroslavsky, *Digital Holography and Digital Image Processing* (Kluwer Academic, 2010).

**3. **E. Leith and J. Upatnieks, “Reconstructed wavefronts and communication theory,” J. Opt. Soc. Am. **52**(10), 1123–1128 (1962). [CrossRef]

**4. **T. M. Kreis and W. P. Juptner, “Suppression of the dc term in digital holography,” Opt. Eng. **36**(8), 2357–2360 (1997). [CrossRef]

**5. **Y. Takaki, H. Kawai, and H. Ohzu, “Hybrid holographic microscopy free of conjugate and zero-order images,” Appl. Opt. **38**(23), 4990–4996 (1999). [CrossRef] [PubMed]

**6. **I. Yamaguchi and T. Zhang, “Phase-shifting digital holography,” Opt. Lett. **22**(16), 1268–1270 (1997). [CrossRef] [PubMed]

**7. **E. Cuche, P. Marquet, and C. Depeursinge, “Spatial filtering for zero-order and twin-image elimination in digital off-axis holography,” Appl. Opt. **39**(23), 4070–4075 (2000). [CrossRef] [PubMed]

**8. **K. Khare and N. George, “Direct coarse sampling of electronic holograms,” Opt. Lett. **28**(12), 1004–1006 (2003). [CrossRef] [PubMed]

**9. **G. L. Chen, C. Y. Lin, M. K. Kuo, and C. C. Chang, “Numerical suppression of zero-order image in digital holography,” Opt. Express **15**(14), 8851–8856 (2007). [CrossRef] [PubMed]

**10. **L. Ma, H. Wang, Y. Li, and H. Zhang, “Elimination of zero-order diffraction and conjugate image in off-axis digital holography,” J. Mod. Opt. **56**(21), 2377–2383 (2009). [CrossRef]

**11. **M. Liebling, T. Blu, and M. A. Unser, “Complex-wave retrieval from a single off-axis hologram,” J. Opt. Soc. Am. A **21**(3), 367–377 (2004). [CrossRef] [PubMed]

**12. **M. Liebling, T. Blu, and M. A. Unser, “Non-linear Fresnelet approximation for interference term suppression in digital holography,” Proc. SPIE **5207**, 553–559 (2003). [CrossRef]

**13. **J. Zhong, J. Weng, and C. Hu, “Reconstruction of digital hologram by use of the wavelet transform,” in *Digital Holography and Three-Dimensional Imaging*, OSA Technical Digest (CD) (Optical Society of America, 2009), paper DWB16.

**14. **J. Zhong and J. Weng, “Reconstruction of digital hologram by use of the wavelet transform,” in *Holography, Research and Technologies*, J. Rosen, ed.(InTech, 2011)

**15. **E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory **52**(2), 489–509 (2006). [CrossRef]

**16. **K. Khare, C. J. Hardy, K. F. King, P. A. Turski, and L. Marinelli, “Accelerated MR imaging using compressive sensing with no free parameters,” Magn. Reson. Med. **68**(5), 1450–1457 (2012). [CrossRef] [PubMed]

**17. **M. Bertero and P. Boccacchi, *Introduction to Inverse Problems in Imaging* (IOP, 1998).

**18. **J. A. Fessler, “Penalized weighted least square image reconstruction for positron emission tomography,” IEEE Trans. Med. Image. **13**(2), 290–300 (1994). [CrossRef]

**19. **D. H. Brandwood, “A complex gradient operator and its application in adaptive array theory,” IEE Proc. H Microwaves Opt. Antennas **130**, 11–16 (1983).

**20. **S. Boyd and L. Vandenberghe, *Convex Optimization* (Cambridge University Press, 2004).

**21. **S. Geman and D. Geman, “Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images,” IEEE Trans. Pattern Anal. Mach. Intell. **PAMI-6**(6), 721–741 (1984). [CrossRef] [PubMed]

**22. **E. Y. Sidky, Y. Duchin, X. Pan, and C. Ullberg, “A constrained, total-variation minimization algorithm for low-intensity x-ray CT,” Med. Phys. **38**(S1Suppl 1), S117–S125 (2011). [CrossRef] [PubMed]

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