Iterative method for zero-order suppression in off-axis digital holography

We propose a method to suppress the so-called zero-order term in a hologram, based on an iterative principle. During the hologram acquisition process, the encoded information includes the intensities of the two beams creating the interference pattern, which do not contain information about the recorded complex wavefront, and that can disrupt the reconstructed signal. The proposed method selectively suppresses the zero-order term by employing the information obtained during wavefront reconstruction in an iterative procedure, thus enabling its suppression without any a priori knowledge about the object. The method is analyzed analytically and its convergence is studied. Then, its performance is shown experimentally. Its robustness is assessed by applying the procedure on various types of holograms, such as topographic images of microscopic specimens or speckle holograms. © 2010 Optical Society of America OCIS codes: (090.1995) Digital holography; (100.2000) Digital image processing; (100.3010) Image reconstruction technique; (100.3175) Interferometric imaging; (120.5050) Phase measurement; (180.3170) Interference microscopy.


Introduction
Digital holography (DH) is a measurement technique based on an interferometric principle, which enables the recovery of both the amplitude and phase of a wavefront. This measurement method has been ported to the investigation of microscopic specimens by combining interferometry and microscopy in digital holographic microscopy (DHM) [1]. This capability led to many different applications, such as high-precision surface topography [2], or biological investigations [3][4][5]. In macroscopic applications, DH enables quality control applications and vibration assessment [6][7][8].
The principle of holography is based on the measurement of the interference between two waves, where one interacted with the specimen, commonly called the object wave, and a reference wave, whose phase profile is precisely controlled. The interference pattern is then measured in intensity on a photographic plate or an electronic system such as a charge-coupled device (CCD) camera. Due to this acquisition process, different terms are recorded: the intensities of both waves, usually referred to as the zero-order term, the object wavefront information, and, due to the real nature of the recorded signal, the complex conjugate of the object wave, called twin-image. This redundant information makes it more difficult to recover an accurate reconstruction of the object signal, and can even perturb it heavily.
Reconstructing the signal containing the information about the object complex wavefront requires separating it from the other terms contained in the interferogram. Historically, two main methods have been employed to accomplish it, either in the spatial or temporal domain. Spatial separation has been performed by using a reference wave which propagation direction makes a slight angle with that of the object wave. In this off-axis configuration, the terms contained in the hologram propagate in different direction, making it possible to isolate the object wavefront [9]. More recently, this spatial separation has been performed in the Fourier domain, either optically [10] or digitally [1].
On the other hand, temporal separation can be achieved through the acquisition of several holograms, giving rise to phase-shifting interferometry (PSI) methods [11]. Classically, PSI requires four holograms with very precise phase-shifts for reconstruction, even though more refined methods were developed, loosening the requirement of high precision phase-shifts [12], or reducing the needed amount of images [13,14]. Recently, hybrid techniques were also developed, working on both temporal and spatial domains [15][16][17].
Different methods were proposed to suppress the zero-order term in off-axis holography. A first method based on multiple measurements uses shutters to directly measure the intensity of the two beams interfering for direct subtraction [18,19]. Intensity terms were considered as constant, making it possible to subtract them with a constant term [20]; then, linear filtering was also used with high-pass filters [21,22]. Non-linear methods were also proposed by employing approximations on the signal along with square operators [23], or wavelet-based transforms by using a Taylor development [24]. The methods cited above always employ an approximation of the interference equation, or restrict the field of application by using strong assumptions on the type of object field. Recently, we also proposed a method providing exact reconstruction by filtering the hologram in the cepstral domain [25], with low requirements on the experimental procedure.
Iterative methods were employed in particular in phase retrieval algorithms [26], and in particular in combination with Fourier interferometry for suppression of the conjugated imaging order [19]. Those methods however can be used for phase-only objects, as this type of algorithms is usually employed for phase recovery from intensity measurements, thus not providing the full complex field.
The method proposed in this article relies on an iterative procedure, which employs the information about the complex object field retrieved during reconstruction to suppress the zeroorder, in order to obtain a new version of the field in which zero-order artifacts are suppressed. Its main advantages are that no approximation is made on the acquisition process, and that it is based on simple physical considerations, improving its robustness.
In Section 2, we present a state-of-the-art reconstruction for off-axis DH, on which part of our proposed method is relying. Then, in Section 3, we present the principle of the iterative technique, and study its convergence. In Section 4, we present the performance of the algorithm with simulated holograms, and analyze the zero-order suppression. Finally, in Section 5, we demonstrate the capabilities of the method with experimental measurements for different situations such as topographic measures of microscopic objects and speckle holograms for low aperture systems.

Off-axis digital holography reconstruction
We present in this section a state-of-the-art method for off-axis hologram reconstruction. This method will be employed during the rest of the article in order to compare performances with the proposed iterative method. Moreover, the principle of the technique described in Section 3 is strongly based on the equations described below.
When acquiring the interference pattern between two waves, the different terms recorded can be expressed as where the two square terms correspond to the intensities of respectively the object and reference waves, or * is the image term, and o * r is its complex conjugate, i.e. the twin image. If the reference wave is considered as a plane wave with a propagation vector subtending an angle with the optical axis, the reference can be expressed as r = R exp [−iϕ(x, y)], where ϕ(x, y) is a linear phase function defining the tilt, and R a real constant value defining the overall amplitude of the plane wavefront. In this case, the different terms of Eq. (1) will be modulated in the spectral domain as where ⊗ represents the convolution operator, ˆ the Fourier transform operator, and (ω 0,x , ω 0,y ) corresponds to the frequency shift induced by the tilt function ϕ(x, y). The reference having a constant intensity, its spectrum results in a Dirac function δ . The result of the interference is illustrated in Fig. 1, where a simulation of a hologram has been performed, resulting in the image shown in Fig. 1(a), to which corresponds the spatial spectrum shown in Fig. 1(b), where the different terms of Eq. (1) are identified. The addition of terms in Eq. (2) shows the potential overlap which can require zero-order suppression. However, the overlap in the spectral domain and in the spatial domain can be different, so that one should note that historically, the off-axis configuration has been proposed in order to provide separation of the different terms in the spatial domain. Therefore, one could always ensure the separation of the terms by increasing the propagation distance of the wavefield between the hologram and image plane. However, the digital propagation induces numerical artifacts which increase with longer propagation distances, and having a long distance between the image plane and the detector can induce vignetting. This explains why most present implementations of DHM employ short reconstruction distances, that make it necessary to use new methods for zero-order suppression.
Under those conditions, the zero-order term is composed of a Dirac function situated at the origin and of the object autocorrelation, which has a spectral bandwidth of twice the one of the imaging term. In the case of DHM, the bandwidth of the imaging term will be limited by the resolution of the microscope objective (MO), that we define here as ω c . In this case, the object autocorrelation will have a spectral bandwidth of [−2ω c , 2ω c ], which can lead to spectral overlap with the imaging order. However, as it can be seen in Fig. 1(b), the spectral bandwidth can appear as smaller in the case of smooth objects. Finally, the twin image has a bandwidth identical to the imaging order, but is modulated in the opposite side of the spectrum. Therefore, in the case the modulation is large enough compared to the bandwidth of the imaging orders, i.e. ω o > ω c , this term will not overlap with the imaging term in the off-axis configuration, and thus can be filtered out. This condition can be in fact usually readily fulfilled in most implementations, as using one quadrant of the Fourier spectrum for the imaging order makes it possible to reach diffraction-limited resolution when any magnifying optics is employed [25].
We consider first in this Section the case where no spectral overlap is occurring between the terms in Eq. (1), leading therefore to perfect reconstruction regards to zero-order filtration. Under this condition, the imaging order can be filtered in the Fourier plane, and can therefore be described as whereŴ (ω x , ω y ) is a window function to select the imaging term. After having selected the imaging term, one can employ a phase mask, commonly called numerical parametric lens (NPL) [27,28], to demodulate the signal, and compensate for aberrations. In the case the hologram was not acquired in the image plane, as it is the case in lensless holography or Fresnel holography, the object wavefront can then be propagated in the image plane by employing the Fresnel integral [29,30]. One should note that nowadays, a transform based on a convolution kernel is commonly preferred to the integral formalism, since it preserves the sampling distance along propagation [30,31]. However, thorough this article, we will use the integral formalism for didactic purposes, as its model is closer to a physical propagation.
We have considered in the basic theoretical development of this Section that no spectral overlap is occurring between the different terms of Eq. (1). However, this implies very strict conditions on the spatial bandwidth of the imaging order, and can lead to a loss in spatial resolution [25]. Suppressing the zero-order makes it possible to avoid being limited by the camera resolution, by loosening the constraints on bandwidth and required tilt between the beams, thus providing optically-limited resolution of reconstruction, by letting one full quadrant to be employed for the imaging order in the case of off-axis DHM.

Functioning principle of the iterative method
The method we propose in this article is based on the principle of approximating the zero-order through the reconstructed signal of Eq. (3) which can be obtained with the method described in Section 2, and to subtract it from the measured hologram in order to get a zero-order-free signal. The main purpose of this approach is to enable object zero-order suppression without any knowledge about the original object wave, so that no further recording is required for suppression. This approach makes it possible to keep the one-shot feature of off-axis DHM.

Iterative method principle
The zero-order can be divided in two terms, which correspond to the intensities of the two waves employed for the interference. The reference intensity |r| 2 can simply be recorded by the camera by blocking the object wave, for instance, in order to get an estimator |r exp | 2 . As the reference wave is commonly well controlled and constant in time, this operation can be considered as a calibration step before the experiment. The fact of using a measurement enables the suppression of the experimental noise occurring on the reference wave. It is however more complicated to compensate for the object wave intensity term |o| 2 , since its profile is specimen-dependent.
In order to get an estimator of the object wave intensity, one can use the reconstructed signal of Eq. (3). This can be done by calculating the autocorrelation of the calculated field, and normalizing it with the experimental reference intensity where it is considered in first approximation that |r exp | 2 ≈ R 2 . This equation provides an estimation of the object wave intensity, so that it is possible to get a hologram with an attenuated zero-order by subtracting those terms as However, as described in Section 2, the different terms of the hologram are often overlapping in the spectral domain, so that Eq. (4) will be disturbed by some part of the zero-order. The approach shown in Eq. (5) has been presented for example in Refs. [18,19], where both the object and reference intensity were measured for further subtraction, thus requiring three images for reconstruction. The principle of the proposed method is therefore to use an iterative approach by employing I of Eq. (5) at step k in order to get a more accurate estimator |o est,k+1 | 2 for further suppression.
In this fashion, the iterative algorithm can be summarized as follows: 1. Compensate for the reference intensity by subtracting the term |r exp | 2 .
4. Compute the hologram with attenuated zero-order I k (x, y) according to Eq. (5).
The performance of the algorithm will mainly depend on how strong the spectral overlap between the terms is, and the proposed method has clearly a usability range. In the case no spectral overlap occurs, the reconstruction will be exact without using zero-order suppression methods. On the other hand, if the spectral overlap is too important, and if the zero-order intensity is rather strong, one can intuitively understand that the algorithm will not provide an exact reconstruction.

Discussion on convergence
This Section does not intend to give a formal demonstration of convergence, since the full mathematical development is out of the scope of this article. However, the equations derived below make it possible to deduce the necessary conditions for convergence of the iterative algorithm, and provide an intuition on the causes of the eventual error sources which can be identified when employing this method for zero-order suppression.
In the case the different terms contained in Eq. (1) are not totally separated in the spectral domain, the filtering operation of Eq. (3) will contain parts of zero-order, as it is usually the case. The estimator of Eq. (4) will not be totally accurate in these conditions, so that the estimated imaging order is The part of zero-order contained in the second term of the right-hand of Eq. (6) will be reflected in the calculation of the estimator, giving rise to several terms. If we express the filtered object intensity as O W (x, y) = |o| 2 ⊗W (x, y), the zero-order estimator at the first iteration becomes One can see four terms in Eq. (7), resulting from the autocorrelation operation. The first is the estimator employed for suppressing the object zero-order as in Eq. (4), while the three others are error terms induced by the presence of some object intensity in the filtering window. The second term corresponds to the square of the object intensity, while the two last ones are modulated cross-terms. The different components of Eq. (7) are represented in Fig. 2 for the hologram of Fig. 1 at first iteration, showing their relative intensity through the multiplying factors used to make the terms visible.
At each iteration, the terms of Eq. (7) will be modified by the filtering and modulus operations. By recurrence, one can deduce a general formula for the error between the estimator and the object zero-order for the iteration step k as where O( f ) is a polynomial function. One can note that the parasitic terms of Eq. (8) are powers of K = o R . This demonstrates that a first condition for the algorithm to converge is that the reference amplitude has to be stronger than that of the object, i.e. R > |o|, so that the summation in Eq. (8) does not diverge. Indeed, the error in the estimation at first iteration will be once again multiplied by the factor K at the next iteration, making them thus decrease.
On the other hand, one can note that the efficiency of the object zero-order suppression will strongly depend on the accuracy of the estimation of the reference intensity |r exp | 2 , in order to minimize the error induced by the first term on the right hand of Eq. (8).
Finally, the error estimation shows that the iterative technique enables efficient zero-order suppression, although it will not fully suppress it. The summation in the error expression will not converge to zero, implying that some artifacts will remain even at convergence. However, as it will be shown in the following Sections, those artifacts can fairly be neglected in practical cases. One should note also that the reconstruction quality will also depend on the reference wave shape. In our model, we have considered a plane reference wave, which leads to the results presented in Section 4. However, a reference wave having a strongly irregular phase profile will induce errors in the reconstruction, since the estimator of the object autocorrelation is based on the imaging order, which is multiplied by the reference wave. The procedure proposed in this article prevents some errors to propagate by normalizing with an experimental reference, but this will not compensate for strongly aberrated reference waves.

Results on synthesized holograms
In this Section, we present the applicability of the proposed method on simulated holograms, which enable analyzing the accuracy and the speed with which the algorithm converges to the solution. The synthetic holograms were created in accordance with the experimental conditions presented in Section 5, so that comparison of efficiency between ideal and real conditions could be made.

Topographic holography
The hologram employed in this case is a representation similar to a topographic image of a USAF 1951 test target, which contains specific spatial frequencies. The hologram and its spectrum are shown in Fig. 1. The optical arrangement generates a spectrum in which some overlap occurs between the imaging order and the object zero-order [cf. Fig. 1(b)]. However, since the object is rather smooth in this case, most of the spectral information is given by frequencies close to the origin, implying that the spectral bandwidth of the object autocorrelation is mostly contained in the [−ω c , ω c ] range.
The synthetic hologram is generated by firstly applying a low-pass filter to the phantom, to account for the aperture of the MO, taken in this case as a 20×, NA = 0.4 MO, with λ = 661 nm and a camera sampling of Δx = 6.45 μm. The object is then propagated to the camera plane (d = 5 cm), where its interference with a tilted plane wave has been calculated to generate the hologram. The ratio K has been taken in this case as 1/5. The reconstruction for the standard method described in Section 2 is shown in Fig. 3(a), where the Fourier filtering step has been performed with a square selecting the proper quadrant of the spectrum [cf. dashed square in Fig. 3(c)]. The zero-order can be identified as a ghost in the bottom-right corner of the figure, highlighted by a dashed square in Fig. 3(a), 3(b). The iterative reconstruction was applied to this hologram, so that the object zero-order is suppressed, as it can be seen from the spectrum of the wavefield in the hologram plane shown in Fig. 3(c). Accordingly, the reconstruction of the hologram shows an imaging order free of ghost in amplitude [cf. Fig. 3(b)]. The results were obtained in this case with 10 iterations, which is plentifully sufficient in those conditions to suppress the zero-order. However, depending on the recording conditions, such as the value of K or the spectral overlap between terms, more iterations could possibly be needed in order to get sufficient extinction. The speed of convergence for this hologram is shown in Fig. 4 by representing the sum of ε-defined in Eq. (8)-on the field of view for each iteration which has been normalized by the object intensity.
One can see that the error converges to a value of 3%, which corresponds to the residue of the algorithm, as described in Section 3.2.

Speckle hologram simulations
In this Section, we present simulation results based on the reconstruction of speckle holograms. In the case speckle is occurring, the object wave contains very high spatial frequencies, created by the random phase pattern generated by the specimen roughness. This can be identified in the spectrum of the generated hologram as an approximately constant spectral energy in the aperture of the MO, as presented in Fig. 5(a).
The synthetic hologram is based on the same object as in Section 4.1, where the phase has been set to a random value uniformly distributed between [−π/2, π/2]. This simulation considers a low aperture optical system, composed of a 2 f lens system. The field diffracted by the specimen is propagated to the first lens of the optical system, where a low-pass filter is applied to account for the optical component aperture. The field is next propagated to the second lens where its curvature is compensated. Then, it is propagated to the hologram plane where the interference with a tilted reference wave is computed. The optical system was taken for a numerical aperture of 0.045, a wavelength of λ = 760 nm and Δx = 6.45 μm. The results of the iterative reconstruction are shown in Fig. 5, with both the spectrum of the field in the hologram plane [cf. Fig. 5(b)] and the reconstructed amplitude in the image plane [cf. Fig. 5(c)], for 10 steps of iteration. As it can be seen in the spectrum of the image, the iterative reconstruction efficiently suppresses the zero-order. However, some spectral energy can still be identified at the center of the spectrum, although not seen in Fig. 5(a). This is due to the residues after convergence, which are larger in this case compared to the smooth hologram of Section 4.1, since the spectral overlap is more important in this case. Nevertheless, one can see from the amplitude reconstruction that the residues are still negligible, giving an imaging order free of artifacts from the zero-order.

Experiments
In this Section, we present experimental results based on the iterative reconstruction. In accordance with Section 4, the measurements were carried on a USAF 1951 test target and speckle measurements based on the reflection on a coin.

Topographic measurements
A USAF 1951 test target has been measured with a standard reflection DHM setup [1], based on a Michelson interferometer. Parameters are similar to the simulation of Section 4.1, so that the specimen is illuminated with a plane wave of wavelength λ = 661 nm, through a 20× MO (NA = 0.4). The reflected light then interferes with a reference wave, which profile is well controlled and curvature matched with the object wave. The interference pattern is recorded by a CCD camera, with a pixel size of Δx = 6.45 μm, for a field of view of 512×512. The target was placed in a configuration where the elements of the target generate strong spatial frequencies in the direction of the interference fringes, giving rise to a spectrum similar to the one used in Section 4.1. The results of the reconstruction are shown in Fig. 6, with the same type of filtering as in Section 4.1, consisting in a quadrant filtering as shown in Fig. 3(c), and a propagation along d = 4.8 cm. The amplitude of the reconstruction for the standard method [cf. Fig. 6(a)] contains a strong part of zero-order, which is suppressed with the iterative technique as shown in Fig. 6(b).
The corresponding phase image is presented in Fig. 6(c) for the iterative case, where it can be identified that the noise on the signal is approximately constant, which is an indicator that the zero-order has been suppressed, so that the phase signal is not disrupted by the phase artifacts of the zero-order. We show in Fig. 7 a phase profile of the steps generated by the chromium layer [cf. dashed line in Fig. 6(c)] for both standard and iterative reconstructions. One can see that the iterative technique keeps the quantitative information of the phase, by providing the same mean phase change on steps. However, the iterative technique is free from a strong phase noise which can be identified on the phase profile of the standard method, due to artifacts generated by the zero-order term.

Speckle hologram
In this Section, we present experimental results based on the iterative method for speckle holograms. We measured a part of coin surface on a reflection setup similar to the one employed in Section 5.1, but with low aperture. In this case, a magnification of 3× is employed, provided by an achromatic lens, with a source wavelength of 760 nm. Apart from those elements, all other experimental parameters are identical.
The results of the reconstruction are shown in Fig. 8, with the hologram spectrum [cf. Fig. 8(a)], in which the uniform spectral energy E generated in the bandwidth of the aperture can be clearly identified, giving a similar signal as in the simulations of Section 4.2. This behavior is due to the roughness of the coin, which generates approximately random phasors in some parts of the field of view.
The reconstructed amplitude of the letter "A" engraved on a Swiss coin for the standard and iterative reconstructions are shown respectively in Fig. 8(b) and 8(d), demonstrating also in this case the suppression of the zero-order. The spectrum of the signal after using the iterative procedure is shown in Fig. 8(c), where the zero-order is very strongly attenuated.
In order to quantify the efficiency of the method, we measured the mean spectral energy in parts of the spectrum for respectively the standard (std) and iterative (iter) techniques, on regions of the zero-order (zero), and in parts containing measurement noise (noise) outside the imaging orders. Those regions are identified in Fig. 8(c). One can estimate the efficiency of the zero-order suppression by calculating The meaning of Eq. (9) is the calculation of the ratio of the spectral energy of the zero-order before and after employing the iterative method, giving an indicator of the relative attenuation of the zero-order on an energy point of view. The subtraction of the spectral energy of a zone outside of the interference terms (where no information is contained) accounts for the electronic noise generated by the camera. This computation gives an efficiency of S = 97.8% for the measurements shown in Fig. 8, indicating that most of the spectral energy carried by the zeroorder has been suppressed.

Conclusions
We presented in this paper an iterative technique for suppressing the zero-order term in off-axis digital holography. The method relies on employing the object field obtained through standard spatial filtering methods to subtract the object autocorrelation term. It is based on the simple physical model of intensity recording, thus making it robust for many different types of signal.
We demonstrated the applicability of the method with different types of holograms, such as topographic phase maps in the context of digital holographic microscopy, or speckle holography for low aperture systems. The functioning principle of the method has been shown through simulations, with which convergence and robustness of the zero-order suppression could be assessed. The efficiency of the method was then demonstrated experimentally, showing a good  agreement with simulated results. Some discrepancies were however identified, as the zeroorder is not perfectly suppressed due to experimental noise and optical aberrations of the system, but only strongly attenuated.
Suppressing the zero-order enables employing a larger spatial bandwidth for the imaging term and therefore makes it possible to increase the spatial resolution of the reconstruction by suppressing the need of minimizing the aperture of the system because of artifacts which can occur due to the object intensity. We also proposed a suppression of the reference term through a calibration measurement, in order to suppress the spatial frequencies coming from the coherent noise of the reference wave. The method employing an iterative reconstruction, it is not suitable for real-time reconstruction. However, it keeps the one-shot capability of offaxis digital holography, as it needs only one hologram for reconstruction. Furthermore, the reconstruction method is independent of the type of measured object, as no a priori information is employed to retrieve the complex field.
This method can be applied by satisfying simple conditions, such as employing a reference wave whose intensity is globally stronger than the one of the object wave. The speed of convergence of the algorithm will essentially depend on how the different terms encoded in the hologram are overlapping, since a larger spectral overlap will make it more difficult for the algorithm to converge. Suppression of the zero-order in our experiments resulted in more accurate phase profile for topographic measurements, or ghost-free images in speckle holography.