Optics and Lasers in Engineering 50 (2012) 322–327
Contents lists available at SciVerse ScienceDirect
Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng
Complex-amplitude-based phase unwrapping method for digital holographic microscopy Shuo Liu n, Wen Xiao, Feng Pan, Fanjing Wang, Lin Cong School of Instrument Science and Optoelectronics Engineering, Beihang University, Beijing 100191, China
a r t i c l e i n f o
a b s t r a c t
Article history: Received 6 September 2011 Received in revised form 26 October 2011 Accepted 8 November 2011 Available online 21 November 2011
A complex-amplitude-based phase unwrapping approach for digital holographic microscopy is proposed in this paper. A quality map is derived directly from the reconstructed complex amplitude distribution of object wave to evaluate noise influence and phase reliability in the wrapped phase image. Quality-guided phase unwrapping algorithm is then implemented with the quality map to retrieve continuous phase profile. Unwrapping errors caused by unreliable phase data are successfully suppressed. The effectiveness of this method is demonstrated with simulation and experimental results. & 2011 Elsevier Ltd. All rights reserved.
Keywords: Phase unwrapping Digital holographic microscopy Quality map
1. Introduction Digital holographic microscopy (DHM) is a useful method permitting nondestructive, marker-free and high-resolution phase contrast imaging, which has been applied in many fields, including living cells investigation [1–3] and microstructure topography [4,5]. In DHM, the object wave transmitted or reflected by the specimen is magnified through a microscope objective. A better resolution is achieved as the magnified image of the specimen is captured instead of the specimen itself. A digital camera is adopted to record the hologram formed by the interference of object wave and reference wave. The object wavefront is numerically reconstructed from the digital hologram. The phase image is then extracted from the complex amplitude image of object wave, which could be used to analyze 3D profile of the specimen. Unlike its natural performance, the calculated phase image is typically wrapped in ( p, p) radians, making 2D phase unwrapping a necessary procedure. Phase unwrapping in DHM is usually challenging as the phase data in wrapped phase image are possibly degraded due to low signal-to-noise ratios (SNRs), which might be caused by severe noise corruption or low signal intensity. Numerous methods have been proposed to minimize the unwrapping errors induced by unreliable phase data. One strategy is based on comparing multiple phase images wrapped in different manners, which could be obtained by using multiple wavelengths [6,7] or varying numerical reconstruction distance [8]. These methods demand additional
n
Corresponding author. Tel./fax: þ 86 10 82339736. E-mail address:
[email protected] (S. Liu).
0143-8166/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.optlaseng.2011.11.006
experimental or computational efforts. Some phase unwrapping algorithms developed for other coherent techniques could be adopted in DHM as well. Most of these methods rely on single phase map and exhibit more efficiency, e.g., branch-cut algorithms [9,10], region-growing algorithm [11], least-squares methods [12,13], etc. Among them, an approach known as quality-guided phase unwrapping algorithm [14–21] is widely used, as it makes a good trade-off between noise robustness and execution time [22,23]. This algorithm utilizes the quality map depicting the goodness of phase data to optimize the unwrapping path, whose performance exceedingly depends on the accuracy of quality map. Conventional quality maps are derived from phase data statistics, e.g., phase derivative variance [22], which sometimes fail to evaluate phase reliability accurately. Quality maps could be defined with the original data obtained in measurement, such as fringe modulation [19], correlation coefficient [22], etc. Unfortunately, these quality measures are developed for specific applications thus not available in DHM. It is important to generate more effective quality maps for DHM phase evaluation and phase unwrapping. DHM enables to record and reconstruct the wavefront of object wave. The signal retrieved in DHM is the complex amplitude image, which directly reflects the noise influence. In this paper, a complex-amplitude-based method for DHM phase unwrapping is proposed. A novel quality map is defined taking into account the relative variation of reconstructed complex amplitude values. The data qualities in different regions of the reconstructed image are successfully evaluated. Quality-guided phase unwrapping is performed with this quality map to recover the continuous phase profile. The availability of this approach is validated with simulation and experimental results.
S. Liu et al. / Optics and Lasers in Engineering 50 (2012) 322–327
2. Principle of the method 2.1. Phase retrieval in DHM Based on the principle of digital holography, the implementation of DHM consists of two main procedures. At first, a digital camera is adopted to record the hologram Hx,y formed by the interference of reference wave Rx,y and object wave Ox,y 2
Hx,y ¼ 9Ox,y þRx,y 9 ,
ð1Þ
where x and y are the spatial coordinates in the camera target plane. Subsequently, the hologram Hx,y is processed with numerical reconstruction algorithms to obtain the complex amplitude distribution of object wave in the image plane. By using convolution algorithm [24], this procedure could be expressed as
Cm,n ¼ DFT
1
fDFTff ilterfHx,y UR0x,y ggUGg,
ð2Þ
where m and n denote the spatial coordinates in the image plane; Cm,n is the reconstructed complex amplitude image; R0x,y denotes the digital reference wave for reconstruction; G denotes discrete transfer function; filter{} denotes digital filtering process; DFT{} and DFT 1{} denote Discrete Fourier Transform and its inverse transform, respectively. Both the amplitude image Am,n and the phase image jm,n could be calculated from Cm,n , given as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Am,n ¼ 9Cm,n 9 ¼ RefCm,n g2 þ ImfCm,n g2 , ð3Þ
jm,n ¼ arctan
ImfCm,n g , RefCm,n g
ð4Þ
where arctan is the four-quadrant arctangent operator, Re{} and Im{} denote the real and imaginary parts, respectively. The phase image jm,n indicates optical path difference, which is involved with the profile information and refractive index distribution of the specimen. As jm,n is wrapped in the range of ( p, p) radians after arctangent calculation, 2D phase unwrapping must be implemented for further analysis. 2.2. Complex-amplitude-based quality map
where jkm,n is the phase of Ckm,n . Eq. (7) indicates that the proposed quality measure is determined by the phase difference between Cm,n and Ckm,n . In a severely noisy area, the differences would be considerable due to the random variations of complex amplitude values, resulting in low quality values. In contrast, in a less noisy area, the complex amplitude values mostly reflect the specimen information, giving rise to smaller differences and higher quality values. Eq. (6) is therefore an evaluation of noise influence and phase data reliability. Unlike conventional quality measures defined with phase data or intensity values, this quality map is derived from the complex amplitude image reconstructed in DHM. It is therefore referred to as complex-amplitude-based (CA) quality map in the following sections. 2.3. Quality-guided phase unwrapping After the extraction of CA quality map, quality-guided phase unwrapping is realized to eliminate 2p discontinuities in the wrapped phase image jm,n . In this procedure, phase pixel with highest quality value is unwrapped first, while pixels with lower quality values are unwrapped subsequently until all pixels are unwrapped. This path-following algorithm is similar regardless of which quality measure is used. A detailed description of this procedure could be found in [22]. With the guidance of CA quality map, the influence of unreliable phase data could be considerably suppressed.
3. Simulation analysis Computer simulation is performed to analyze the capability of proposed method. As for DHM, it has been demonstrated that the real noise sources involved in recording give rise to complex Gaussian noise after numerical reconstruction [25], which is similar to the property of speckle noise [26]. The complex amplitude image Cm,n retrieved in DHM is hereby simulated by disturbing object wavefront Sm,n with additive complex Gaussian noise Nm,n , expressed as
Cm,n ¼ Sm,n þ Nm,n ¼ am,n expðjjm,n Þ þ LðRm,n þ jIm,n Þ,
In DHM, the phase data in jm,n are easily degraded due to noise corruption, which would cause phase errors in unwrapping procedure. Both the specimen information and noise corruption are involved with the complex amplitude distribution Cm,n , only exhibiting different properties. Generally, the spatial variation of Cm,n induced by noise is disordered, while that caused by specimen information is more regular. This difference is hereby utilized to estimate the noise influence. The phase quality at the pixel (m, n) is derived from the variation of complex amplitude values in the neighborhood of the pixel. At first, the complex amplitude values Ci,j within the k k window centered at (m, n) are summed up, given as follows:
Ckm,n ¼
X
k1 m k1 2 r i rm þ 2
i,j
k1 n k1 2 rj r n þ 2
Ci,j ,
,
ð5Þ
where k is an odd integer. Subsequently, Ckm,n is compared with the complex amplitude value at the center pixel (i.e.,Cm,n ), by calculating the normalized inner product of them. The quality value is defined as Z m,n ¼
ReðCkm,n Cnm,n Þ 9Ckm,n 99Cm,n 9
,
ð6Þ
where the asterisk * denotes the complex conjugate. Eq. (6) is equivalent to Z m,n ¼ cosðjkm,n jm,n Þ,
ð7Þ
323
ð8Þ pffiffiffiffiffiffiffi where j ¼ 1; am,n and jm,n are the amplitude and phase distributions associated with simulated object wave, respectively; L denotes the noise level; Rm,n and Im,n are the real and imaginary parts of simulated noise, respectively, which obey Gaussian random statistics with mean values of 0 and standard deviations of 1. The object phase distribution jm,n is constructed taking advantage of computer processing, resulting in continuous phase profile with values in (0, 30) radians, as illustrated in Fig. 1(a). Since the reconstructed amplitude image is usually nonuniform in practice, am,n is simulated by peaks function in MATLAB with values in (1, 10), as shown in Fig. 1(b). Fig. 1(c) and (d) present the wrapped phase image and amplitude image calculated from the resultant complex amplitude image Cm,n with noise level L¼1.5, respectively, which simulate the results reconstructed in DHM. The proposed method is employed in phase unwrapping of Fig. 1(c) to confirm its availability. The quality-guided methods based on two other quality measures are also utilized for the purpose of comparison. One of these two quality maps is the phase derivative variance (PDV) map, which is considered as one of the most effective phase-based quality measures [22]. As the amplitude image in Fig. 1(d) is the square root of reconstructed wave field intensity, it reflects data reliability to some extent and could be used as quality measure as well, which is referred to as AMP quality map. Fig. 2(a), (b) and (c) illustrate the unwrapped phase images obtained with the methods based on CA quality map (k¼3), PDV quality map (k¼3) and AMP quality map,
324
S. Liu et al. / Optics and Lasers in Engineering 50 (2012) 322–327
respectively. The regions indicated by the dashed rectangles in Fig. 2(a), (b) and (c) are enlarged and shown in Fig. 2(d), (e) and (f), respectively. Compared with Fig. 1(a), it is clear that the continuous phase profile of simulated object wave is successfully
Fig. 1. Simulated DHM result: (a) and (b) are the phase image and amplitude image of simulated object wave, respectively; (c) and (d) are the wrapped phase image and amplitude image with noise corruption (L¼ 1.5), respectively. All of them are 500 500 pixels.
restored with the proposed CA method. The results obtained with PDV and AMP quality maps are involved with remarkable phase discontinuities and unwrapping errors in the regions with low SNRs. All the phase unwrapping calculations were carried out on computer with a Core 2 DUO E7200 processor. The execution times corresponding to CA (k¼3), PDV (k ¼3) and AMP methods are 1.2709 s, 1.4698 s and 1.0788 s, respectively. Furthermore, the performances of these three approaches are compared in the case of different noise levels. The phase errors are quantitatively estimated by calculating the averaged absolute phase differences between the unwrapped phase images and the simulated phase map in Fig. 1(a). Fig. 3 presents the averaged phase errors associated with these three methods under noise level L¼1.0,
Fig. 3. Averaged phase errors associated with the methods based on CA (k¼3), PDV (k¼3) and AMP quality maps when noise level L¼ 1.0, 1.1, 1.2 y 2.0.
Fig. 2. Unwrapped phase maps of Fig. 1(c) based on the methods: (a) CA (k¼ 3); (b) PDV (k¼ 3); (c) AMP; (d), (e) and (f) are the enlarged regions marked by the dashed rectangles in (a), (b) and (c), respectively. (a), (b) and (c) are 500 500 pixels, (d), (e) and (f) are 100 100 pixels.
S. Liu et al. / Optics and Lasers in Engineering 50 (2012) 322–327
1.1, 1.2 y 2.0. It is noticed that the proposed approach is involved with the smallest phase errors in these three methods, which proves its effectiveness in phase unwrapping.
Fig. 4. DHM experimental result: (a) digital hologram; (b) wrapped phase image retrieved from (a); (c) enlarged region of interest in (b); (d) amplitude image corresponding to (c). (a) and (b) are 1024 1024 pixels, (c) and (d) are 400 400 pixels.
325
4. Experimental verification The availability of proposed method is further verified with DHM experimental results. The living samples of mouse osteoblastic MC3T3-E1 cells were investigated in transmission mode with an improved DHM imaging system [3]. A frequency-doubled Nd:YAG laser (532 nm, 50 mW) was used as the light source. An infinity-corrected microscope objective (20 , NA¼0.4) was adopted to improve lateral resolution. A CMOS camera (1024 1024, pixel size 6.7 mm 6.7 mm) was utilized to record the hologram. The phase aberrations induced by microscope objective and liquid culture medium were compensated in numerical reconstruction. Fig. 4(a) shows a captured digital hologram of 1024 1024 pixels. Fig. 4(b) is the wrapped phase image retrieved from Fig. 4(a) with reconstruction distance of 17.8 cm. The region of interest marked by the dashed rectangle in Fig. 4(b) is enlarged and illustrated in Fig. 4(c), which is 400 400 pixels. Although the living cells are usually considered as pure phase objects, it is found that the reconstructed amplitude image is not homogeneous in this experiment. Fig. 4(d) is the amplitude distribution corresponding to Fig. 4(c). The proposed CA method (k ¼3) is then used to unwrap the phase map in Fig. 4(c). Similar to that in simulation section, the PDV quality map (k ¼3) and AMP quality map, which takes advantage of Fig. 4(d) in this case, are also adopted to guide phase unwrapping procedure. The unwrapped phase images obtained with these three methods are shown in Fig. 5. It can be seen that phase discontinuities still exist in all unwrapped images, especially in the areas with low amplitude modulation, as shown in Fig. 5(d)–(f). The unweighted L0 measure values are calculated to evaluate the averaged numbers of phase discontinuities in these results [22]. The unweighted L0 measure values corresponding to Fig. 5(a), (b) and (c) are 0.0090, 0.0107 and 0.0128, respectively. It means that the proposed CA approach is involved with the fewest phase discontinuities. Fig. 6 presents another experimental result of
Fig. 5. Unwrapped phase images of Fig. 4(c) based on the methods: (a) CA (k¼ 3); (b) PDV (k¼ 3); (c) AMP; (d), (e) and (f) are the enlarged regions marked by the dashed rectangles in (a), (b) and (c), respectively. (a), (b) and (c) are 400 400 pixels; (d), (e) and (f) are 150 150 pixels.
326
S. Liu et al. / Optics and Lasers in Engineering 50 (2012) 322–327
Fig. 6. DHM experimental result: (a) wrapped phase image; (b) amplitude image corresponding to (a); unwrapped phase images obtained with the methods: (c) CA (k¼ 3); (d) PDV (k¼3); (e) AMP; (f), (g) and (h) are the enlarged regions marked by the dashed rectangles in (c), (d) and (e), respectively. (a)–(e) are 300 300 pixels; (f), (g) and (h) are 150 150 pixels.
MC3T3-E1 cells obtained with the same DHM system. The reconstruction distance is 18.4 cm. Fig. 6(a) and (b) illustrate the wrapped phase image and its corresponding amplitude image, respectively. Fig. 6(c), (d) and (e) are the unwrapped phase images acquired with CA (k¼3), PDV (k¼3) and AMP methods, respectively. The unweighted L0 measure values calculated with Fig. 6(c), (d) and (e) are 0.0081, 0.0102 and 0.0126, respectively. It also validates that the proposed CA method is useful in suppressing phase unwrapping errors. All the results prove that the proposed CA method offers an effective tool for phase quality evaluation and phase unwrapping in DHM.
5. Conclusion In this paper, a DHM phase unwrapping method is realized by using a quality map extracted from the reconstructed complex amplitude image. The difference between specimen information and noise distribution is utilized to estimate phase reliability at each pixel, with which quality-guided unwrapping approach is performed. Computer simulation demonstrates that this method
is effective in suppressing unwrapping errors induced by unreliable phase data. Experimental results further prove its feasibility in DHM phase unwrapping in the case of noise corruption. This method provides an alternative tool to evaluate phase quality and recover phase profile of the specimen in DHM. It might be improved and employed in more application fields in the future.
Acknowledgments This work was supported by the National Science Foundation of China (No. 31000387 and No. 61177006). References [1] Rappaz B, Marquet P, Cuche E, Emery Y, Depeursinge C, Magistretti P. Measurement of the integral refractive index and dynamic cell morphometry of living cells with digital holographic microscopy. Opt Express 2005;13(23): 9361–73. ¨ [2] Kemper B, Carl D, Schnekenburger J, Bredebusch I, Schafer M, Domschke W, et al. Investigation of living pancreas tumor cells by digital holographic microscopy. J Biomed Opt 2006;11(3). 034005-1-8.
S. Liu et al. / Optics and Lasers in Engineering 50 (2012) 322–327
[3] Liu S, Pan F, Wang Z, Wang F, Rong L, Shang P, et al. Long-term quantitative phase-contrast imaging of living cells by digital holographic microscopy. Laser Phys 2011;21(4):740–5. ¨ [4] Charrie re F, Kuhn J, Colomb T, Montfort F, Cuche E, Emery Y, et al. Characterization of microlenses by digital holographic microscopy. Appl Opt 2006;45(5):829–35. [5] Miccio L, Finizio A, Grilli S, Vespini V, Paturzo M, De Nicola S, et al. Tunable liquid microlens arrays in electrode-less configuration and their accurate characterization by interference microscopy. Opt Express 2009;17(4):2487–99. [6] Gass J, Dakoff A, Kim MK. Phase imaging without 2p ambiguity by multiwavelength digital holography. Opt Lett 2003;28(13):1141–3. [7] Parshall D, Kim MK. Digital holographic microscopy with dual-wavelength phase unwrapping. Appl Opt 2006;45(3):451–9. [8] Khmaladze A, Epstein T, Chen Z. Phase unwrapping by varying the reconstruction distance in digital holographic microscopy. Opt Lett 2010;35(7): 1040–2. [9] Goldstein RM, Zebker HA, Werner CL. Satellite radar interferometry: twodimensional phase unwrapping. Radio Sci 1988;23(4):713–20. [10] Zheng D, Da F. A novel algorithm for branch cut phase unwrapping. Opt Laser Eng 2011;49(5):609–17. [11] Baldi A. Phase unwrapping by region growing. Appl Opt 2003;42(14):2498–505. [12] Pritt MD, Shipman JS. Least-squares two-dimensional phase unwrapping using FFTs. IEEE Trans Geosci Remote Sensing 1994;32(3):706–8. [13] Chyou JJ, Chen SJ, Chen YK. Two-dimensional phase unwrapping with a multichannel least-mean-square algorithm. Appl Opt 2004;43(30):5655–61.
327
[14] Bone DJ. Fourier fringe analysis: the two-dimensional phase unwrapping problem. Appl Opt 1991;30(25):3627–32. [15] Quiroga JA, Gonzalez-Cano A, Bernabeu E. Phase-unwrapping algorithm based on an adaptive criterion. Appl Opt 1995;34(14):2560–3. [16] Xu W, Cumming I. A region growing algorithm for InSAR phase unwrapping. IEEE Proc 1996;4:2044–6. [17] Pritt MD. Phase unwrapping by means of multigrid techniques for interferometric SAR. IEEE Trans Geosci Remote Sensing 1996;34(3):728–38. [18] Li W, Su XY. Phase unwrapping algorithm based on phase fitting reliability in structured light projection. Opt Eng 2002;41(6):1365–72. [19] Quan C, Tay CJ, Chen L, Fu Y. Spatial-fringe-modulation-based quality map for phase unwrapping. Appl Opt 2003;42(35):7060–5. [20] Yu C, Peng Q. A correlation-based phase unwrapping method for Fourier transform profilometry. Opt Laser Eng 2007;45(6):730–6. [21] Kemao Q, Gao W, Wang H. Windowed Fourier-filtered and quality-guided phase-unwrapping algorithm. Appl Opt 2008;47(29):5420–8. [22] Ghiglia DC, Pritt MD. Two-dimensional phase unwrapping: theory, algorithms, and software. Wiley; 1998. [23] Zappa E, Busca G. Comparison of eight unwrapping algorithms applied to Fourier-transform profilometry. Opt Laser Eng 2008;46(2):106–16. [24] Kreis T. Handbook of holographic interferometry. Wiley-VCH; 2005. [25] Pandey N, Hennelly B. Effect of additive noise on phase measurement in digital holographic microscopy. 3D Research 2011;2(1). 01006-1-6. [26] Goodman JW. Speckle phenomena in optics: theory and applications. Roberts & Co; 2007.