Fast focus estimation using frequency analysis in digital holography

Seungtaik Oh; Chi-Young Hwang; Il Kwon Jeong; Sung-Keun Lee; Jae-Hyeung Park

doi:10.1364/OE.22.028926

1. Introduction

Autofocusing is one of the most intriguing features of digital holography(DH). Through an autofocusing process, the focus plane for a digital hologram can be determined without any prior knowledge of recording situation. The most popular autofocusing method in DH is the range search approach based on focus measure. Various focus measures have been used for autofocusing in DH, for instance, self-entropy [1], Fresnelet-sparsity [2], total amplitude [3], the variance of local intensity [4], the Laplacian of Gaussian for multiple objects [5] and recently approximated Tamura coefficient [6]. In general, the range search method consists of two generic procedures: sectional numerical reconstructions and focus analysis. Sectional numerical reconstructions are performed for field propagation along the optical axis within a search range chosen by a user. Then, a focus measure is applied for each sectional reconstruction result to find the correct focus distance. The range search methods usually give highly accurate results, but a main drawback here is that the overall process might be time-consuming if the search range is completely out of the correct focus distance and a very small step size in propagation is chosen for a high accuracy.

In this paper we are going to propose a novel fast frequency-based method to estimate the focus distance from a digital hologram for a single object. Our intuitive idea is that optical rays starting from the hologram plane should converge and form an image at the focus plane. The ray direction is determined by the gradient of phase distribution in an isotropic medium. Unfortunately, the phase derivative is hard to compute exactly since the true phase restoration is not easy. In our method, smoothed-rays are considered instead relying on the energy flow direction which is more general definition of the ray direction. The smoothed-rays are computed by using the local maxima of the local spatial frequency spectrum obtained by windowed Fourier transform(WFT). Then the focus distance is estimated by analyzing the distribution of intersection points of the smoothed-rays. Our method is an intrinsic approach based on local spatial frequency analysis on the hologram plane and furthermore is more efficient than the existing approaches since our method does not require any range searches based on numerical reconstructions and a focus measure.

It is remarked that another intrinsic autofocusing method for optical scanning holography was presented in [7]. The authors modified the optical field through some filtering processes and then extracted the focus distance by using the rotational property of the Wigner distribution function in propagation. In addition, WFT has proved to be an effective tool to analyze the fringe pattern and recover the phase information [8, 9].

2. Theory

2.1. Local spatial frequency and smoothed-ray direction

Let f(x, y) be the complex amplitude of an optical wave field with wavelength λ on the xy-plane and let k = 2π/λ be the wavenumber. It is assumed that the optical wave is propagating toward the positive z-direction. WFT [10] of f for a point (x, y) and a spatial frequency pair (α, β) is defined by

W F_{f} (x, y, α, β) = \iint f (x^{'}, y^{'}) W (x^{'} - x, y^{'} - y) exp [- j 2 π (α x^{'} + β y^{'})] d x^{'} d y^{'},

where W(x, y) is a window function, that is, lowpass signal with ∬W(x, y)dxdy = 1. Using WFT, it is possible to obtain a local frequency spectrum at a given point in space. An important property of WFT is that the spectrogram S_f(x, y, α, β) := ‖WF_f(x, y, α, β)‖² can be interpreted as an energy of f at a given spatial point and frequency [11]. Given a spatial point (x, y), the direction of energy flow can be computed by finding a spatial frequency pair (α̅, β̅) to maximize the spectrogram S_f(x, y, α, β). Specifically, the direction is given by the spatial frequency vector associated with (α̅, β̅), i.e., (λα̅, λβ̅,(1 − λ²(α̅² + β̅²))^1/2). In fact, the direction of energy flow is the same as the ray direction which is normal to the wavefront, that is, iso-contour of the phase distribution [12].

For a general optical wave, there exist multiple local peaks in its spectrogram. It can be assumed that each local peak represents the direction of the local energy flow associated with a source of object wave. Given a spatial point (x, y), let n be the number of the local maxima of the spectrogram S_f at (x, y). (Let ${(α_{i} β_{i})}_{i = 1}^{n}$ be the frequencies which are the local maximizers of S_f at (x, y). It is sensible to find an averaged direction of energy flow for the multiple local peaks and so smoothed-ray direction is defined by

d (x, y) = \frac{w_{1} {\tilde{d}}_{1} + \dots + w_{n} {\tilde{d}}_{n}}{‖ w_{1} {\tilde{d}}_{1} + \dots + w_{n} {\tilde{d}}_{n} ‖},

where w_i = S_f (x, y, α_i, β_i) and

{\tilde{d}}_{i} = (λ α_{i}, λ β_{i}, {(1 - λ^{2} (α_{i}^{2} + β_{i}^{2}))}^{1 / 2})

. By definition, the smoothed-ray direction is the normalized weighted sum of the spatial frequency vectors associated with local peaks. For convenience, the energy flow direction chosen by the global maximum is called non-smoothed ray direction and clearly the smoothed-ray direction is exactly the same as the non-smoothed ray direction when the spectrogram has a single peak, for example, in a spherical wave case.

2.2. Focus distance estimation based on analysis of intersections of rays

Now, we explain how to compute the focus distance from given smoothed-rays. Let m be the number of the smoothed-rays. Let ${p_{k} = (x_{k}, y_{k}, 0)}_{k = 1}^{m}$ and ${d_{k} = (a_{k}, b_{k}, c_{k})}_{k = 1}^{m}$ be the ray sampling points and the associated smoothed-ray directions, respectively. Let ${L_{k}}_{k = 1}^{m}$ be smoothed-rays given by {p_k = (x_k, y_k, 0)} and {d_k = (a_k, b_k, c_k)}. First, we need to obtain all possible intersections of {L_k}. An intersection is calculated for three mutually non-parallel smoothed-rays by applying a least square approach, i.e., finding the minimizer of the sum of squared distances to the three rays, which leads to solve Ax = b such that

A = [\begin{array}{l} 3 - \sum a_{k_{i}}^{2} & - \sum a_{k_{i}} b_{k_{i}} & - \sum a_{k_{i}} c_{k_{i}} \\ - \sum a_{k_{i}} b_{k_{i}} & 3 - \sum b_{k_{i}}^{2} & - \sum b_{k_{i}} c_{k_{i}} \\ - \sum a_{k_{i}} c_{k_{i}} & - \sum b_{k_{i}} c_{k_{i}} & 3 - \sum c_{k_{i}}^{2} \end{array}] and b = [\begin{array}{l} \begin{array}{l} \sum (x_{k_{i}} - x_{k_{i}} a_{k_{i}}^{2} - y_{k_{i}} a_{k_{i}} b_{k_{i}}) \\ \sum (y_{k_{i}} - x_{k_{i}} a_{k_{i}} b_{k_{i}} - y_{k_{i}} b_{k_{i}}^{2}) \\ \sum (- x_{k_{i}} a_{k_{i}} c_{k_{i}} - y_{k_{i}} b_{k_{i}} c_{k_{i}}) \end{array} \end{array}],

where each summation is taken for i = 1, 2, 3.

Once all intersection calculations are done, we collect the positive z-coordinates of intersections {z₁, ..., z_N} to be analyzed for focus distance extraction. Let $z_{\min} = {min}_{i = 1}^{N} {z_{i}}$ , $z_{\max} = {max}_{i = 1}^{N} {z_{i}}$ and $\bar{z} = (\sum_{i = 1}^{N} z_{i}) / N$ . The three quantities z_min, z_max and z̅ provide fundamental information for the focus distance, but for more accurate and sophisticated estimation, we further consider a density distribution 𝒟 for {z_i} over an interval I = [z_min, z_max]. In order to compute 𝒟, a uniform grid G needs to be constructed on I. Each grid cell of G has a density which is the number of z_i’s contained in the grid cell. Intuitively, the focus distance is closely related to a high density area. Based on this, we focus on the local maxima {d₁, ..., d_M} of 𝒟. Let {c₁, ..., c_M} be the centers of grid cells giving local maxima of 𝒟. Finally our focus distance estimator z̃ is defined by the density-weighted average of {c_i}:

\tilde{z} = \frac{\sum_{i = 1}^{M} d_{i} c_{i}}{\sum_{i = 1}^{M} d_{i}} .

3. Numerical results and discussions

In all our numerical examples, we have used M₄ kernel function [13] $W_{h} (x, y) = \frac{1}{h^{2}} \tilde{W} ({(x^{2} + y^{2})}^{1 / 2} / h)$ , where

\tilde{W} (r) = \frac{10}{7 π} {\begin{array}{l} 1 - \frac{3}{2} r^{2} + \frac{3}{4} r^{3} & : r \leq 1 \\ \frac{1}{4} {(2 - r)}^{4} & : 1 < r \leq 2 . \\ 0 & : r > 2 \end{array}

Our window function is computationally efficient since W̃(r) consists of cubic polynomials and has a finite support. Note that the window size is 4h in this case. In order to analyze the local features of an optical signal, the window size needs to small enough. However, if the window size is too small, then many important details in the spectrogram are smoothed out, and the accuracy of our method may be degraded. In all our simulations, the ratio of the window size to the hologram size is set to be 4%.

3.1. Spherical wave

Let us consider a converging spherical wave field with wavelength λ = 633nm defined on the xy-plane: f (x, y) = exp(−jkr)/r, where $r = \sqrt{x^{2} + y^{2} + z_{0}^{2}}$ and z₀ > 0. The wave field f (x, y) is discretized on a computational domain H with N × N pixels and a pixel size Δp, where N = 2048 and Δp = 9.765625μm. The horizontal or vertical size of the square domain H is L = NΔp = 20mm. We choose 16 sampling points uniformly distributed on H for ray computation. For each ray sampling point, the smoothed-ray direction is computed by Eqs. (1) and (2). With the smoothed-rays, we calculate the focus distance by Eqs. (3) and (4). For this spherical wave case, the density distribution is similar to a delta function defined at the source distance.

The errors between the original source distance and the estimated distance predicted by our method are shown in Fig. 1. It is seen that for 0.3 ≤ z₀ ≤ 5, our estimation is very sharp and the error given by |z̃ − z₀|/z₀ is less than 0.25%(see Fig. 1(a)). The Nyquist sampling condition implies that the source distance z₀ needs to be larger than the Nyquist distance z_nyquist = L/(λΔp) = 0.3085m. If z₀ < z_nyquist, we have a noisy wave field on H and the accuracy of our estimation is degraded, which is shown in Fig. 1(b). For z₀ > 10, our method has a tendency that the accuracy becomes worse as the source distance increases. In order to make an intersection at far distance, we need rays with very steep slopes, but the error in numerical computation of intersection for steep rays is considerable and furthermore the slopes of rays are bounded in our discrete situation.

Fig. 1 Focus estimation error for spherical wave with the source at (a) 0.3 ≤ z₀ ≤ 20 and (b) short distances z₀ ≤ 0.4.

Download Full Size | PDF

3.2. Holograms for virtual object

We are going to test our estimation method for a planar cross-shaped object with the normal parallel to the z-axis. The same computation domain H as the previous case is used. The object size is set to be 10mm which is half of the size of H. The object wave is computed by setting a random phase configuration on the object surface. The amplitude and the phase profiles of the object wave on H are shown in Figs. 2(a) and 2(b). We again use the same ray sampling points as the previous case. For the object z-depth z₀ = 1, our method estimates the focus distance z̃ = 1.00736 with an error |z̃− z₀|/z₀ = 0.73%. The numerical reconstruction result at the estimated focus distance z̃ is shown in Fig. 2(c). The distribution of the intersection density has a noisy bell shape whose center is located at the original focus distance as shown in Fig. 3.

Fig. 2 The wave profiles, (a) the amplitude and (b) the phase, for a hologram of for a planar cross-shaped object with the focus distance z₀ = 1 and (c) the numerical reconstruction at the estimated focus distance z̃ = 1.00736.

Download Full Size | PDF

Fig. 3 The intersection density for a planar cross-shaped object.

Download Full Size | PDF

In addition, comparison results between the smooth-ray and the non-smoothed ray are presented in Fig. 4(a). An additional result for the twice scaled object is also shown in Fig. 4(b) to see how the object size affects the accuracy of our method. The smoothed ray provides more stable and accurate estimation results than the non-smoothed ray for most z-depths. It is emphasized that for the larger object, the non-smoothed ray produces extremely bad estimations as shown in Fig. 4(b), but our smoothed-ray produces fairly good estimations the errors of which are almost the same as the smaller object case.

Fig. 4 Estimation error comparison between the smoothed-rays and the non-smoothed rays for planar cross-shaped objects of size (a) 10mm and (b) 20mm.

Download Full Size | PDF

For z₀ = 1.3 in Fig. 4(b), the estimation error is reduced from 65.28% to 2.62% by using the smoothed-ray. The error reducing is because the smoothness in ray direction change is improved significantly by our smoothed-ray. This can be verified in Fig. 5 which is the projections of ray direction vectors to the hologram plane. In addition, the same projections for the smaller object case is presented in Fig. 6.

Fig. 5 Projected rays to the hologram plane for the planar cross-shaped object of size 20mm with the Z-depth at z = 1.3: (a) smoothed ray and (b) non-smoothed ray.

Download Full Size | PDF

Fig. 6 Projected rays to the hologram plane for the planar cross-shaped object of size 10mm with the Z-depth at z = 1.3: (a) smoothed ray and (b) non-smoothed ray.

Download Full Size | PDF

Big errors for the smoothed-ray(Fig. 4(b)) are observed around z₀ = 0.4, which results from the fact that the Nyquist distance which is inversely related to the object size is determined as z_nyquist = 0.617 and so the object is located too close to the hologram plane.

3.3. Phase-shifting hologram for a real coin

A phase-shifting hologram of a real coin is tested as well. The amplitude and phase profiles of the coin hologram are shown in Figs. 7(a) and 7(b). The wavelength is λ= 532nm and the hologram has Δp = 2.2μm and 1940 × 1940 pixels. The parameters are determined by the recording situation. 16 smoothed rays are used to estimate the focus distance. The correct focus distance of the hologram is 36cm. Our method estimates the focus distance at 35.6507cm and the numerical reconstruction result at this distance is shown in Fig. 7(c). The intersection density of the smoothed-rays is shown in Fig. 8.

Fig. 7 The wave profiles, (a) the amplitude and (b) the phase, for a phase-shifting hologram of a real coin and (c) the numerical reconstruction result at the estimated focus distance z̃ = 0.356507(brightness enhanced).

Download Full Size | PDF

Fig. 8 Intersection density for a phase-shifting hologram of a real coin.

Download Full Size | PDF

As an additional remark, we need to explain what distance is computed by our method for a hologram of two objects with different z-depths z₁ and z₂(z₁ < z₂). Let d = z₂ − z₁. If d is small enough, our method predicts a single focus distance around the averaged distance z̅ = (z₁ + z₂)/2. But if d is large, the estimated focus distance is shifted toward the near object since the energy of the optical flow of the far object becomes much weaker than that of the near object.

3.4. Computation time

The most expensive part in our method is the WFT calculation to find ray directions. The computation time for the other processes such as ray intersection and density analysis is negligible. So, the computational complexity of our method is almost equal to that of m operations of FFT, where m is the number of rays. For all simulations presented in this paper, our focus estimation computation is done in 6 seconds.

In order to see how fast our method is, a comparison result with a range search method is also presented. The Tamura coefficient studied in [6] is used as the focus measure. For the same examples, a fixed search range [0.3, 3] and a step size Δz = 0.01 are considered and it takes about 190 seconds to compute the focus distance. It is noted that the computation time of a range search method mainly depends on the number of field propagations and therefore the computation time of the test range search method may increase if the search range is enlarged or Δz is reduced for better accuracy.

At the current stage, all computations in our method are carried out on CPU and it is highly expected that a significant speed up can be obtained if our method is implemented on GPU.

4. Conclusion

A novel fast intrinsic method to compute the focus distance from digital hologram has been proposed. Our method computes the object depth by analyzing the intersection points of optical rays starting from the hologram. The smoothed-ray direction is introduced instead of the pure ray direction. The smoothed-ray direction is based on the local maxima of the local spatial frequency spectrum obtained by the windowed Fourier transform and, as a result, it can avoid the computation of derivatives of phase distribution of the optical wave field. Since the local spatial frequency information on the hologram plane is used, our method does not require a time-consuming range search based on sequential numerical reconstructions and focus measure both of which are essential in the previous auto-focusing methods for digital hologram. At present, our current method cannot estimate multiple focus distances for a hologram with multiple objects. To show effectiveness of our method, numerical computations have been carried out for holograms for virtual and real objects. The results show that our method give highly accurate focus estimations for practically sensible distances in a high speed.

Acknowledgments

The authors would like to thank all anonymous reviewers for their valuable comments and suggestions. This work was supported by the ‘Cross-Ministry Giga KOREA Project’ of the Ministry of Science, ICT and Future Planning, Republic of Korea [ GK14C0100, Development of Interactive and Realistic Massive Giga-Content Technology].

References and links

1. J. Gillespie and R. A. King, “The use of self-entropy as a focus measure in digital holography,” Pattern Recogn. Lett. 9, 19–25 (1989). [CrossRef]

2. M. Liebling and M. Unser, “Autofocus for digital fresnel holograms by use of a fresnelet-sparsity criterion,” J. Opt. Soc. Am. A 21, 2424–2430 (2004). [CrossRef]

3. F. Dubois, C. Schockaert, N. Callens, and C. Yourassowsky, “Focus plane detection criteria in digital holography microscopy by amplitude analysis,” Opt. Express 14, 5895–5908 (2006). [CrossRef] [PubMed]

4. L. Ma, H. Wang, Y. Li, and H. Jin, “Numerical reconstruction of digital holograms for three-dimensional shape measurement,” J. Opt. A: Pure Appl. Opt. 6, 396–400 (2004). [CrossRef]

5. M. Tachiki, M. Itoh, and T. Yatagai, “Simultaneous depth determination of multiple objects by focus analysis in digital holography,” Appl. Opt. 47, D144–D153 (2008). [CrossRef] [PubMed]

6. P. Memmolo, C. Distante, M. Paturzo, A. Finizio, P. Ferraro, and B. Javidi, “Automatic focusing in digital holog-raphy and its application to stretched holograms,” Opt. Lett. 36, 1945–1947 (2011). [CrossRef] [PubMed]

7. T. Kim and T.-C. Poon, “Autofocusing in optical scanning holography,” Appl. Opt. 48, H153–H159 (2009). [CrossRef] [PubMed]

8. Q. Kemao, “Windowed fourier transform for fringe pattern analysis,” Appl. Opt. 43, 2695–2702 (2004). [CrossRef] [PubMed]

9. G. Rajshekhar and P. Rastogi, “Estimation of the phase derivative using an adaptive window spectrogram,” J. Opt. Soc. Am. A 27, 69–75 (2010). [CrossRef]

10. F. Hlawatsch and G. Boudreauz-Bartels, “Linear and quadratic time-frequency signal representations,” IEEE Signal Process. Mag. 9, 21–67 (1992). [CrossRef]

11. H. M. Ozaktas, Z. Zalevsky, and M. A. Kutay, The Fractional Fourier Transform with Applications in Optics and Signal Processing (John Wiley & Sons, 2000).

12. J. W. Goodman, Introduction to Fourier Optics (Roberts & Company Publishers, 2005), 3rd ed.

13. K. Toraichi, M. Kamada, S. Itahashi, and R. Mori, “Window functions represented by b-spline functions,” IEEE T. Acoust. Speech 37, 145–147 (1989). [CrossRef]

Fast focus estimation using frequency analysis in digital holography

Abstract

1. Introduction

2. Theory

2.1. Local spatial frequency and smoothed-ray direction

2.2. Focus distance estimation based on analysis of intersections of rays

3. Numerical results and discussions

3.1. Spherical wave

3.2. Holograms for virtual object

3.3. Phase-shifting hologram for a real coin

3.4. Computation time

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (8)

Equations (5)

Optics Express