- Email: [email protected]
Effective improvement of depth resolution and reduction of ripple error in depth-resolved wavenumber-scanning interferometry
Optics and Lasers in Engineering 66 (2015) 58–63
Contents lists available at ScienceDirect
Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng
Effective improvement of depth resolution and reduction of ripple error in depth-resolved wavenumber-scanning interferometry Yun Zhang, Yulei Bai, Jinxiong Xu, Weichao Xu, Yanzhou Zhou n Faculty of Automation, Guangdong University of Technology, Guangzhou 510006, PR China
art ic l e i nf o
a b s t r a c t
Article history: Received 28 February 2014 Received in revised form 2 August 2014 Accepted 11 August 2014
The bottleneck problem hindering the development of Depth-Resolved Wavenumber-Scanning Interferometry (DRWSI) is a limited depth resolution due to a finite range of the laser wavenumber scanning. A robust Complex Number Least Squares Algorithm (CNLSA) is proposed in the present study to take full advantage of phase and amplitude information over the entire range of frequencies, and the algorithm performs much better than any window function. Experimental results as well as simulations show that the phase ripple error is reduced and the depth resolution is significantly refined. A CNLSA is likely to become practical and suitable for DRWSI applications. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Interferometric imaging Phase measurement Lasers tunable Resolution
1. Introduction Depth-Resolved Wavenumber-Scanning Interferometry (DRWSI) and Phase-Contrast Spectral OCT are perspective imaging techniques that utilize the amplitude and phase information for all pixels [1–5]. One of the bottleneck problems hindering the development of these methods is a limited depth resolution due to a relatively small range of the wavenumber scanning in the laser's output. To overcome this problem, two possible solutions were recently suggested: (1) using a laser with a wider range of the wavenumber scanning [2–4]; and (2) synthesizing several uncorrelated wavenumber bands into one. The latter method can be extended to synthesize the outputs from several different lasers into a wide range of wavenumber scanning bands [2,3,5]. DRWSI, an extension of the traditional Phase-Shifting Interferometry (PSI), is a linear superposition of all possible interference signals from the scattered lights within the structure. If there are M correlated reflective lights along the depth dimension, there will be M(M 1)/2 amplitude peaks due to the interference [6]. Hibino and co-workers positioned the reference plane properly away from the measured object to make those interference peaks appear in the high frequency band; meanwhile, the signals in the low frequency band were filtered out as noise [6], although these signals carried useful information. As early as 1984, Greivenkamp proposed a method of least squares for the PSI in order to obtain the unwrapped phase map by fitting the
phase shifting data to a sine curve of known period but variable amplitude and phase [7]. Grattan and co-workers proposed a curve fitting signal processing for a interferometric in order to supress the noise [8]. In general, least squares algorithm, suitable for the analysis of functions or series of real numbers, is optimization methods that are widely used in signal and imaging processing. However, after a real number signal in DRWSI is transformed into a series of complex numbers in the frequency domain, its real and imaginary parts contain the depth-resolved information, respectively. Consequently, it is complicated to use least squares method for DRWSI. In the present work, consistent with [9], a new algorithm is proposed for evaluation of wrapped phase maps of interference signals from the depth-scattered multiple reflections. The proposed method takes full advantage of phase and amplitude information in the entire range of frequencies, and has higher depth resolution.
2. Theory Let us assume that the measured object is composed of multiple transparent surfaces S1, S2, ……, SM along its depth dimension; the wavenumbers of the scanned light source are k(n), n¼ 0, 1, …, N 1; for simplicity, the spatial coordinates (x, y) are omitted in the analysis that follows. A k-space Fourier transform of the interference signals from those multiple surfaces M can be written as [2–5]: pffiffiffiffiffiffiffiffiffiffiffi ðN 1Þ U I p U I q ~ f 7 Λpq ; exp½8 i U ð2 U kð0Þ U Λpq þ ϕpq0 Þ U W π 2UΔ k p¼1q¼1 M
~ Þ¼ ∑ Iðf n
Corresponding author. Tel.: þ 86 15918764670. E-mail address: [email protected] (Y. Zhou).
http://dx.doi.org/10.1016/j.optlaseng.2014.08.009 0143-8166/& 2014 Elsevier Ltd. All rights reserved.
M
∑
ð1Þ
Y. Zhang et al. / Optics and Lasers in Engineering 66 (2015) 58–63
59
where I is the light's intensity; Λ is the Optical Path Difference (OPD); the subscripts p and q denote the surfaces p and q, respectively;
ϕpq0 is the initial phase difference between the
surfaces p and q; Δk is the range of the wavenumber scanning; ~ ðf Þ is the k-space Fourier transform of the window function. and W
1/2 ∙(N 1)/ Let the unknown vector X ¼[∑M p ¼ 1 I p ∙(N 1)/Δk, (I1∙I2)
Δk, …, (I(M 1)∙IM)1/2∙(N 1)/Δk, Λ12/π, …, Λ(M 1)M/π, 2∙k(0)∙Λ12 þ ϕ120, …, 2∙k(0)∙Λ(M 1)M þ ϕ(M 1)M0], which represents the amplitude of the DC, the amplitude (Ip∙Iq)1/2∙(N 1)/Δk, the frequency Λpq/π and the phase [2∙k(0)∙Λ12 þ ϕ120] wrapped into the region [ π, π], at the peak of the interference signal Λpq, respectively. The
length of the unkonwn vector X is [3∙M∙(M 1)/2] þ1. ~ Þ are, respectively, the measured spectrum and If I~m ðf Þ and Iðf the spectrum that is modeled by Eq. (1) in the DRWSI, the error E between those spectra can be written by utilizing the definition of the least squares for complex numbers [9]: L
~ Þ I~m ðf Þj2 ; E ¼ ∑ jIðf j j j¼1
ð2Þ
where L is the number of data points in the frequency domain. In order to solve for the unknown vector X, first of all its initial values should be given in terms of their values at the relative peaks Λpq in the measured spectrum. Secondly, the nonlinear least squares fitting method (Levenberg–Marquardt algorithm) is used to minimize the error E and to iteratively compute the vector X. This method is called the Complex Number Least Squares Algorithm (CNLSA) for the DRWSI. Because CNLSA takes full advantage of the information in the entire range of frequencies, it can recover the real OPD Λpq and the phase precisely, even though the OPD between the two surfaces p and q is smaller than the depth resolution of the DRWSI by using Fourier transform with rectangular window. It should be noted that when the range of wavenumber scanning is too narrow, it can cause severe measured spectrum aliasing. In this case, the initial values of X evaluating from the peaks are very far from their real values. Sometimes, the algorithm fails. Thus, choosing the initial values of X properly is the convergence factor of CNLSA. ~ ðf Þ has two undesirable In Eq. (1), the window function W effects on the OPD and phase evaluation in the DRWSI: (1) spoiling its depth resolution by broadening its spectrum; and (2) producing amplitude and phase fluctuations. The depth resolution of a Fourier transform with Hanning window is δΛ ¼4π/Δk and the peak level of its sidelobes is 31 dB, while the depth resolution of a Fourier transform with rectangular window is δΛ ¼ 2π/Δk and the peak level of its sidelobes is 13 dB. Hanning window is characterized by smaller spectral leakage; however, its mainlobe is twice wider than that of the rectangular window. There is no window function that has the narrowest possible mainlobe and the smallest possible sidelobes.
3. Optical setup As shown in Fig. 1, the DRWSI experimental system was designed to be a Michelson interferometer. The light source was a single longitude mode DFB laser diode Toptica LD-0860-0150-DFB-1. Its wavenumber k could be modulated by the temperature without the mode hop. Three batches of wavenumber scanning experiments were performed to verify the phase distribution of the DRWSI. As shown in Fig. 1, the series of the wavenumber scanning k1, k2 and k3 start at the same value of k(0)¼ 8.001 106 m 1, with the corresponding ranges being Δk1 ¼8.182 103 m 1, Δk2 ¼6.168 103 m 1 and Δk3 ¼ 4.307 103 m 1, respectively. The object in the experiments was composed of three flat smooth surfaces S1, S2 and
Fig. 1. The configuration of the optical system: TC—temperature controller; LD— laser diode; LC—laser controller; PC—personal computer; CCD—CCD camera; CBS— cube beam splitter; L—convex lens; OW—optical wedge; AB—the series of the wavenumber scanning k3; AC—the series of the wavenumber scanning k2; and AD—the series of the wavenumber scanning k1.
S3, where S1 and S2 were the two surfaces of an optical wedge with the OPD of Λ120 ¼ 4.68 mm (x¼0 mm, y¼0 mm) and the title angle of 120 . The surface S3 was a flat surface with the OPD of Λ230 ¼ 3.82 mm (x¼0 mm, y¼0 mm). The difference between the OPD Λ120 and Λ230 was dΛ ¼ |Λ230 Λ120|¼ 0.86 mm. Wavenumber monitoring and a Random Sampling Fourier Transform (RSFT) were used to evaluate the wrapped phase maps. A detailed explanation is given in [5].
4. Experimental result Due to the similarity between the OPD Λ12 and Λ23 interference signals, the latter one was chosen to verify the correctness and show the excellent performance of CNLSA. The unwrapped phase maps are performed according to the reference [10]. Ideally, the wrapped and unwrapped phase map at the OPD Λ23 should be smooth, because the interference signals come from the flat smooth surfaces; however, a spectral crosstalk from other interference peaks, in particular from the OPD Λ12, inevitably leads to the appearance of phase ripples at the OPD Λ23. With respect to k1, k2 and k3, because all depth resolutions of the Fourier transform with a Hanning window are much larger than dΛ, all of the wrapped and unwrapped phase maps at the OPD Λ23 are smashed, as shown in Figs. 2(a–c) and 3(a–c), and the phase textures cannot be resolved. With respect to k1, the depth resolution of the Fourier transform with rectangular window of 0.77 mm is a little lower than dΛ. Although there are phase ripples in the wrapped and unwrapped phase map of the OPD Λ23 that is shown in Figs. 2(d) and 3(d), these phase maps are nearly acceptable. With respect to k2 and k3, their depth resolutions of the Fourier transform with rectangular window are 1.02 mm and 1.46 mm, respectively, and are much larger than dΛ; therefore, their phase ripples are larger, as shown in Figs. 2(e, f) and 3(e, f). Because the range of the wavenumber series k3 is the narrowest, its phase maps evaluated by the Fourier transform with rectangular window are totally smashed, as shown in Figs. 2(f) and 3(f), and its phase textures cannot be resolved. In the context below, CNLSA's window function in Eq. (1) is a rectangular window. The initial values of X can be obtained by two steps: (a) recognize the peaks corresponding to DC, surface S23, surface S12 and surface S13 in the measured amplitude spectrum; (b) evaluate the amplitudes,
60
Y. Zhang et al. / Optics and Lasers in Engineering 66 (2015) 58–63
Fig. 2. Wrapped phase maps at the OPD Λ23, for the wavenumber series k1, k2 and k3, respectively; (a–c) Fourier transform with Hanning window; (d–f) Fourier transform with rectangular window; (g–i) CNLSA; (a, d, g) the wavenumber series k1; (b, e, h) the wavenumber series k2; and (c, f, i) the wavenumber series k3.
Fig. 3. Unwrapped phase maps of the OPD Λ23, for the wavenumber series k1, k2 and k3, respectively; (a–c) Fourier transform with Hanning window; (d–f) Fourier transform with rectangular window; (g–i) CNLSA; (a, d, g) the wavenumber series k1; (b, e, h) the wavenumber series k2; and (c, f, i) the wavenumber series k3.
Y. Zhang et al. / Optics and Lasers in Engineering 66 (2015) 58–63
frequencies and wrapped phases at the peaks in the measured spectrum respectively. We use Matlab 2010b on a computer (CPU: Intel Xeon E5645, 2.40 Ghz 12, memory: 8 GB) to do the computation. The time consumption to evaluate CNLSA's wrapped phase maps corresponding to k1, k2 is about 60 h and it becomes about 85 h as the wavenumber scanning is k3. With respect to k1 and k2, no ripples are observed in the wrapped and unwrapped phase maps when CNLSA is utilized, as shown in Figs. 2(g, h) and 3 (g, h). With respect to k3, in spite of poor performance of the other two algorithms, the phase maps evaluated by CNLSA recover, as
Table 1 The maximal unwrapped phase errors at y¼0.14 mm in the wrapped phase maps of the OPD Λ23 evaluated by FT with HW, FT with RW and CNLSA, respectively; WN —wavenumber series; FT—Fourier transform; HW—hanning window; RW—rectangular window; A—absolute maximal phase error; R—relative maximal phase error. WN
k1 k2 k3
FT with HW
FT with RW
CNLSA
A(rad)
R (%)
A (rad)
R (%)
A (rad)
R (%)
88.19 100.31 108.93
321.7 365.9 397.4
1.60 1.99 52.69
5.8 7.3 192.2
– 0.28 2.11
– 1.0 7.7
61
shown in Figs. 2(i) and 3(i). Although there are ripples in the phase maps of CNLSA corresponding to k3, they are still acceptable. In order to compare with the differences among the three algorithms, the section lines at y¼0.14 mm in the unwrapped phase maps of Fig. 3 are chosen to be analyzed. If the unwrapped phase with respect to k1 that is evaluated by CNLSA is taken as the reference, the unwrapped phase differences between the CNLSA, the Fourier transform with Hanning window or rectangular window and the reference are given in Table 1. With respect to k3, the Fourier transform with rectangular window yields a huge phase error, while the CNLSA is still convergent and satisfactory, which unwrapped phase error is similar to that of the Fourier transform with rectangular window with respect to k1 and k2. Another experiment is performed as the series of the wavenumber scanning is k4, where k0 ¼7.308 106 m 1 and Δk¼ 4.503 103 m 1. The object measured is still compose of 3 smooth surfaces S1, S2 and S3, where S1 and S2 are the two surfaces of the optical wedge with the OPD Λ120 ¼ 9.06 mm at x¼0 mm, y¼0 mm and the title angle 60 . Surface S3 is of curvature and set to parallel to the surface S1, with the OPD Λ130 ¼7.69 mm at x¼ 0 mm, y¼0 mm. The interference forms a closed fringe pattern as shown in Fig. 4. Fig. 5(a–f) shows the wrapped and unwrapped phase maps of Λ13 evaluated by Fourier transform with Hanning window, Fourier transform with rectangular window and CNLSA. The depth resolutions of Fourier transform with Hanning window and rectangular window are 2.79 mm and 1.40 mm, respectively, which is larger than dΛ ¼|Λ120 Λ130|¼1.37 mm. Therefore, there are phase ripple errors in the phase maps evaluated by those two methods, as shown in Fig. 5(a, b, d, e). Among the three algorithms, the phase maps evaluated by CNLSA are smoothest, as shown in Fig. 5(c and f). According to the two group experimental results above, it is concluded that the CNLSA is the best among the algorithms for the DRWSI.
5. Simulation
Fig. 4. The closed interference fringe pattern.
The simulations were performed to take into account the optical setup that was described above: the interference signals were from the three surfaces S1, S2 and S3, with the corresponding reflective light intensities of I1 ¼0.58, I2 ¼0.50 and I3 ¼0.21; the initial phases were ϕ120 ¼ ϕ230 ¼ ϕ130 ¼ 0 rad; the wavenumber
Fig. 5. The wrapped and unwrapped phase maps of the OPD Λ13, for the wavenumber series k4; (a, b, c) are the wrapped phase maps evaluated by Fourier transform with Hanning window, Fourier transform with rectangular window and CNLSA, respectively; (d, e, f) are the unwrapped phase maps of (a, b, c), respectively.
62
Y. Zhang et al. / Optics and Lasers in Engineering 66 (2015) 58–63
series was k3. As Λ230 ¼ 3.82 mm, the parameter Λ120 was varied from 7.50 mm to 2.50 mm in steps of 0.005 mm, while the OPD and phase at the OPD Λ230 discrepancies between the preset values and those that were evaluated by a Fourier transform with the Hanning window or rectangular window and CNLSA were simulated. The spectral crosstalk between Λ120 and Λ230 yields the OPD and phase errors by the following mechanisms. (1) Peak drifting occurs as dΛ ¼ |Λ230 Λ120| approaches the depth resolution δΛ. The depth resolutions of Fourier transforms with a rectangular or Hanning window are 1.46 mm and 2.92 mm, respectively. As dΛ approaches its depth resolution from the right side of Fig. 6(a), all three algorithms produce fluctuating errors both in the OPD and phase values. In order to evaluate the fluctuating errors, the ranges of the OPD Λ120 scanning are chosen to be just outside the regions in which double peaks switch to a single peak for all the algorithms,
as shown in Table 2. It is noted that there are relatively larger OPD and phase's offset errors when a Fourier transform with rectangular window is used, because of its poor performance in restraining the sidelobes and the constant OPD between the OPD Λ230 and the DC signal. Among the three algorithms, the offset and standard deviation errors are minimal for the CNLSA, with the finest depth resolution, as shown in Table 2. (2) Peak overlapping occurs as dΛ become smaller than the depth resolution δΛ. Because the phase term 2∙k(0)∙Λpq in Eq. (1) is very sensitive to the Λpq, the double peaks of the OPD Λ230 and Λ120 rapidly switch between the double peaks and a single peak, as dΛ o0.94 mm or dΛ o1.58 mm, for a Fourier transform with rectangular window or Hanning window, respectively. As shown in Fig. 6(a), these cause rapid switching between small and large OPD and phase errors as a result of the automatic peak-searching algorithm. Due to the error in the initial values of X when automatic peak searching is used, the switching error or the depth resolution of the CNLSA occurs at dΛ ¼ δΛ ¼0.80 mm, as shown in Fig. 6. In the experiment described above, with respect to k3, dΛ ¼ 0.86 mm is smaller than the switching point of a Fourier transform with rectangular window (0.94 mm); thus, the wrapped phase map is smashed, as shown in Fig. 2(f). Because dΛ is larger than the switching point of CNLSA (0.80 mm), its wrapped phase map recovers, as shown in Fig. 2(i). Therefore, both the experimental and the simulated results show that CNLSA yields a good OPD and phase performance, which cannot be achieved by using window functions. It is noted that the more correct the initial values X are, the finer the depth resolution of the CNLSA becomes. If the initial values of the OPD Λ230 and Λ120 are given by the preset values in the simulation in spite of the double peaks to a single peak switching, the OPD Λ230 error is 0.07 mm, and the phase error is 0.25 rad (λ/50) at dΛ ¼ δΛ ¼ 0.34 mm, as shown in Fig. 6(b). In this case, the phase error 0.25 rad (λ/50) is similar to the CNLSA error obtained in the experiment with respect to k2, as shown in Table 1; however, the depth resolution is much finer than the experimental one.
6. Conclusion
Fig. 6. The OPD and the wrapped phase errors at the OPD Λ230, simulated by the three algorithms, respectively; green is a Fourier transform with Hanning window; blue is a Fourier transform with rectangular window; red is CNLSA, with respect to the wavenumber scanning series k3; Pha. Err.—the wrapped phase error; OPD Err.— the OPD error; (a) a Fourier transform with Hanning window or rectangular window and CNLSA, automatic peak searching; (b) CNLSA, as the initial values of the vector X are the preset values. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
The spectral leakage from the neighboring peaks is one of the principal error sources in DRWSI. We proposed a robust algorithm, called the Complex Number Least Squares Algorithm (CNLSA), for DRWSI. The algorithm yields a significant improvement of the depth resolution, and suppresses the OPD and phase ripple errors without requiring any increase in hardware complexity. Although the Levenberg–Marquardt algorithm that was used in the present study is a classical method, both the experiment and the simulation suggest that CNLSA performs much better than any window function. Theoretically, it is possible to obtain a better depth resolution by developing a global fitting algorithm that does not depend on the initial values. In conclusion, a CNLSA helps remove the limitation on the range of light source's wavenumber scanning, and is likely to become the principal data evaluating algorithm for DRWSI.
Table 2 The OPD Λ230 and phase errors simulated by FT with HW, FT with RW and CNLSA, respectively, with respect to the wavenumber series k3; std—standard deviation error. Error
OPD (mm) Phase (rad)
FT with HW Λ120 A(6.46 mm, 7.02 mm)
FT with RW Λ120 A (5.00 mm, 5.56 mm)
CNLSA Λ120 A (4.82 mm, 5.38 mm)
Offset
Offset
Offset
7 8.05 10 7 0.03
Std 3
7 4.66 10 7 0.02
3
7 0.14 7 0.70
Std 70.08 70.35
7 3.06 10 7 0.01
Std 3
7 2.50 10 3 7 0.01
Y. Zhang et al. / Optics and Lasers in Engineering 66 (2015) 58–63
Acknowledgment The authors wish to thank the National Natural Science Foundation of China (NSFC) and Provincial Natural Science Foundation of Guangdong (NSFG) for their financial support with the Grants 11072063 and S2012010010327, respectively. References [1] De la Torre-Ibarra MH, Ruiz PD, Huntley JM. Simultaneous measurement of inplane and out-of-plane displacement fields in scattering media using phasecontrast spectral optical coherence tomography. Opt Lett 2009;34(6):806–8. [2] Davila A, Huntley JM, Pallikarakis C, Ruiz PD, Coupland JM. Simultaneous wavenumber measurement and coherence detection using temporal phase unwrapping. Appl Opt 2012;51(5):558–67. [3] Davila A, Huntley JM, Pallikarakis C, Ruiz PD, Coupland JM. Wavelength scanning interferometry using a Ti: sapphire laser with wide tuning range. Opt Laser Eng 2012;50(8):1089–96.
63
[4] Chakraborty S, Ruiz PD. Measurement of all orthogonal components of displacement in the volume of scattering materials using wavelength scanning interferometry. J Opt Soc Am A 2012;29(9):1776–85. [5] Xu J, Liu Y, Dong B, Bai Y, Hu L, Shi C, et al. Improvement of the depth resolution in depth-resolved wavenumber-scanning interferometry using multiple uncorrelated wavenumber bands. Appl Opt 2013;52(20):4890–7. [6] Hibino K, Hanayama R, Burke J, Oreb BF. Tunable phase-extraction formulae for simultaneous shape measurement of multiple surfaces with wavelengthshifting interferometry. Opt Express 2004;12(23):5579–94. [7] Greivenkamp JE. Generalized data reduction for heterodyne interferometry. Opt Eng 1984;23:350–2. [8] Wang Q, Ning YN, Grattan KTV, Palmer AW. A curve fitting signal processing scheme for a white-light interferometric system with a synthetic source. Opt. Laser. Technol 1997;29(7):371–6. [9] Sheppard RJ, Jordan BP, Grant EH. Least squares analysis of complex data with application to permittivity measurements. J Phys D: Appl Phys 1970;3 (11):1759–64. [10] Buckland JR, Huntely JM, Turner SRE. Unwrapping noisy phase maps by use of a minimum-cost-matching algorithm. Appl Opt 1995;34(23):5100–8.
Contents lists available at ScienceDirect
Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng
Effective improvement of depth resolution and reduction of ripple error in depth-resolved wavenumber-scanning interferometry Yun Zhang, Yulei Bai, Jinxiong Xu, Weichao Xu, Yanzhou Zhou n Faculty of Automation, Guangdong University of Technology, Guangzhou 510006, PR China
art ic l e i nf o
a b s t r a c t
Article history: Received 28 February 2014 Received in revised form 2 August 2014 Accepted 11 August 2014
The bottleneck problem hindering the development of Depth-Resolved Wavenumber-Scanning Interferometry (DRWSI) is a limited depth resolution due to a finite range of the laser wavenumber scanning. A robust Complex Number Least Squares Algorithm (CNLSA) is proposed in the present study to take full advantage of phase and amplitude information over the entire range of frequencies, and the algorithm performs much better than any window function. Experimental results as well as simulations show that the phase ripple error is reduced and the depth resolution is significantly refined. A CNLSA is likely to become practical and suitable for DRWSI applications. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Interferometric imaging Phase measurement Lasers tunable Resolution
1. Introduction Depth-Resolved Wavenumber-Scanning Interferometry (DRWSI) and Phase-Contrast Spectral OCT are perspective imaging techniques that utilize the amplitude and phase information for all pixels [1–5]. One of the bottleneck problems hindering the development of these methods is a limited depth resolution due to a relatively small range of the wavenumber scanning in the laser's output. To overcome this problem, two possible solutions were recently suggested: (1) using a laser with a wider range of the wavenumber scanning [2–4]; and (2) synthesizing several uncorrelated wavenumber bands into one. The latter method can be extended to synthesize the outputs from several different lasers into a wide range of wavenumber scanning bands [2,3,5]. DRWSI, an extension of the traditional Phase-Shifting Interferometry (PSI), is a linear superposition of all possible interference signals from the scattered lights within the structure. If there are M correlated reflective lights along the depth dimension, there will be M(M 1)/2 amplitude peaks due to the interference [6]. Hibino and co-workers positioned the reference plane properly away from the measured object to make those interference peaks appear in the high frequency band; meanwhile, the signals in the low frequency band were filtered out as noise [6], although these signals carried useful information. As early as 1984, Greivenkamp proposed a method of least squares for the PSI in order to obtain the unwrapped phase map by fitting the
phase shifting data to a sine curve of known period but variable amplitude and phase [7]. Grattan and co-workers proposed a curve fitting signal processing for a interferometric in order to supress the noise [8]. In general, least squares algorithm, suitable for the analysis of functions or series of real numbers, is optimization methods that are widely used in signal and imaging processing. However, after a real number signal in DRWSI is transformed into a series of complex numbers in the frequency domain, its real and imaginary parts contain the depth-resolved information, respectively. Consequently, it is complicated to use least squares method for DRWSI. In the present work, consistent with [9], a new algorithm is proposed for evaluation of wrapped phase maps of interference signals from the depth-scattered multiple reflections. The proposed method takes full advantage of phase and amplitude information in the entire range of frequencies, and has higher depth resolution.
2. Theory Let us assume that the measured object is composed of multiple transparent surfaces S1, S2, ……, SM along its depth dimension; the wavenumbers of the scanned light source are k(n), n¼ 0, 1, …, N 1; for simplicity, the spatial coordinates (x, y) are omitted in the analysis that follows. A k-space Fourier transform of the interference signals from those multiple surfaces M can be written as [2–5]: pffiffiffiffiffiffiffiffiffiffiffi ðN 1Þ U I p U I q ~ f 7 Λpq ; exp½8 i U ð2 U kð0Þ U Λpq þ ϕpq0 Þ U W π 2UΔ k p¼1q¼1 M
~ Þ¼ ∑ Iðf n
Corresponding author. Tel.: þ 86 15918764670. E-mail address: [email protected] (Y. Zhou).
http://dx.doi.org/10.1016/j.optlaseng.2014.08.009 0143-8166/& 2014 Elsevier Ltd. All rights reserved.
M
∑
ð1Þ
Y. Zhang et al. / Optics and Lasers in Engineering 66 (2015) 58–63
59
where I is the light's intensity; Λ is the Optical Path Difference (OPD); the subscripts p and q denote the surfaces p and q, respectively;
ϕpq0 is the initial phase difference between the
surfaces p and q; Δk is the range of the wavenumber scanning; ~ ðf Þ is the k-space Fourier transform of the window function. and W
1/2 ∙(N 1)/ Let the unknown vector X ¼[∑M p ¼ 1 I p ∙(N 1)/Δk, (I1∙I2)
Δk, …, (I(M 1)∙IM)1/2∙(N 1)/Δk, Λ12/π, …, Λ(M 1)M/π, 2∙k(0)∙Λ12 þ ϕ120, …, 2∙k(0)∙Λ(M 1)M þ ϕ(M 1)M0], which represents the amplitude of the DC, the amplitude (Ip∙Iq)1/2∙(N 1)/Δk, the frequency Λpq/π and the phase [2∙k(0)∙Λ12 þ ϕ120] wrapped into the region [ π, π], at the peak of the interference signal Λpq, respectively. The
length of the unkonwn vector X is [3∙M∙(M 1)/2] þ1. ~ Þ are, respectively, the measured spectrum and If I~m ðf Þ and Iðf the spectrum that is modeled by Eq. (1) in the DRWSI, the error E between those spectra can be written by utilizing the definition of the least squares for complex numbers [9]: L
~ Þ I~m ðf Þj2 ; E ¼ ∑ jIðf j j j¼1
ð2Þ
where L is the number of data points in the frequency domain. In order to solve for the unknown vector X, first of all its initial values should be given in terms of their values at the relative peaks Λpq in the measured spectrum. Secondly, the nonlinear least squares fitting method (Levenberg–Marquardt algorithm) is used to minimize the error E and to iteratively compute the vector X. This method is called the Complex Number Least Squares Algorithm (CNLSA) for the DRWSI. Because CNLSA takes full advantage of the information in the entire range of frequencies, it can recover the real OPD Λpq and the phase precisely, even though the OPD between the two surfaces p and q is smaller than the depth resolution of the DRWSI by using Fourier transform with rectangular window. It should be noted that when the range of wavenumber scanning is too narrow, it can cause severe measured spectrum aliasing. In this case, the initial values of X evaluating from the peaks are very far from their real values. Sometimes, the algorithm fails. Thus, choosing the initial values of X properly is the convergence factor of CNLSA. ~ ðf Þ has two undesirable In Eq. (1), the window function W effects on the OPD and phase evaluation in the DRWSI: (1) spoiling its depth resolution by broadening its spectrum; and (2) producing amplitude and phase fluctuations. The depth resolution of a Fourier transform with Hanning window is δΛ ¼4π/Δk and the peak level of its sidelobes is 31 dB, while the depth resolution of a Fourier transform with rectangular window is δΛ ¼ 2π/Δk and the peak level of its sidelobes is 13 dB. Hanning window is characterized by smaller spectral leakage; however, its mainlobe is twice wider than that of the rectangular window. There is no window function that has the narrowest possible mainlobe and the smallest possible sidelobes.
3. Optical setup As shown in Fig. 1, the DRWSI experimental system was designed to be a Michelson interferometer. The light source was a single longitude mode DFB laser diode Toptica LD-0860-0150-DFB-1. Its wavenumber k could be modulated by the temperature without the mode hop. Three batches of wavenumber scanning experiments were performed to verify the phase distribution of the DRWSI. As shown in Fig. 1, the series of the wavenumber scanning k1, k2 and k3 start at the same value of k(0)¼ 8.001 106 m 1, with the corresponding ranges being Δk1 ¼8.182 103 m 1, Δk2 ¼6.168 103 m 1 and Δk3 ¼ 4.307 103 m 1, respectively. The object in the experiments was composed of three flat smooth surfaces S1, S2 and
Fig. 1. The configuration of the optical system: TC—temperature controller; LD— laser diode; LC—laser controller; PC—personal computer; CCD—CCD camera; CBS— cube beam splitter; L—convex lens; OW—optical wedge; AB—the series of the wavenumber scanning k3; AC—the series of the wavenumber scanning k2; and AD—the series of the wavenumber scanning k1.
S3, where S1 and S2 were the two surfaces of an optical wedge with the OPD of Λ120 ¼ 4.68 mm (x¼0 mm, y¼0 mm) and the title angle of 120 . The surface S3 was a flat surface with the OPD of Λ230 ¼ 3.82 mm (x¼0 mm, y¼0 mm). The difference between the OPD Λ120 and Λ230 was dΛ ¼ |Λ230 Λ120|¼ 0.86 mm. Wavenumber monitoring and a Random Sampling Fourier Transform (RSFT) were used to evaluate the wrapped phase maps. A detailed explanation is given in [5].
4. Experimental result Due to the similarity between the OPD Λ12 and Λ23 interference signals, the latter one was chosen to verify the correctness and show the excellent performance of CNLSA. The unwrapped phase maps are performed according to the reference [10]. Ideally, the wrapped and unwrapped phase map at the OPD Λ23 should be smooth, because the interference signals come from the flat smooth surfaces; however, a spectral crosstalk from other interference peaks, in particular from the OPD Λ12, inevitably leads to the appearance of phase ripples at the OPD Λ23. With respect to k1, k2 and k3, because all depth resolutions of the Fourier transform with a Hanning window are much larger than dΛ, all of the wrapped and unwrapped phase maps at the OPD Λ23 are smashed, as shown in Figs. 2(a–c) and 3(a–c), and the phase textures cannot be resolved. With respect to k1, the depth resolution of the Fourier transform with rectangular window of 0.77 mm is a little lower than dΛ. Although there are phase ripples in the wrapped and unwrapped phase map of the OPD Λ23 that is shown in Figs. 2(d) and 3(d), these phase maps are nearly acceptable. With respect to k2 and k3, their depth resolutions of the Fourier transform with rectangular window are 1.02 mm and 1.46 mm, respectively, and are much larger than dΛ; therefore, their phase ripples are larger, as shown in Figs. 2(e, f) and 3(e, f). Because the range of the wavenumber series k3 is the narrowest, its phase maps evaluated by the Fourier transform with rectangular window are totally smashed, as shown in Figs. 2(f) and 3(f), and its phase textures cannot be resolved. In the context below, CNLSA's window function in Eq. (1) is a rectangular window. The initial values of X can be obtained by two steps: (a) recognize the peaks corresponding to DC, surface S23, surface S12 and surface S13 in the measured amplitude spectrum; (b) evaluate the amplitudes,
60
Y. Zhang et al. / Optics and Lasers in Engineering 66 (2015) 58–63
Fig. 2. Wrapped phase maps at the OPD Λ23, for the wavenumber series k1, k2 and k3, respectively; (a–c) Fourier transform with Hanning window; (d–f) Fourier transform with rectangular window; (g–i) CNLSA; (a, d, g) the wavenumber series k1; (b, e, h) the wavenumber series k2; and (c, f, i) the wavenumber series k3.
Fig. 3. Unwrapped phase maps of the OPD Λ23, for the wavenumber series k1, k2 and k3, respectively; (a–c) Fourier transform with Hanning window; (d–f) Fourier transform with rectangular window; (g–i) CNLSA; (a, d, g) the wavenumber series k1; (b, e, h) the wavenumber series k2; and (c, f, i) the wavenumber series k3.
Y. Zhang et al. / Optics and Lasers in Engineering 66 (2015) 58–63
frequencies and wrapped phases at the peaks in the measured spectrum respectively. We use Matlab 2010b on a computer (CPU: Intel Xeon E5645, 2.40 Ghz 12, memory: 8 GB) to do the computation. The time consumption to evaluate CNLSA's wrapped phase maps corresponding to k1, k2 is about 60 h and it becomes about 85 h as the wavenumber scanning is k3. With respect to k1 and k2, no ripples are observed in the wrapped and unwrapped phase maps when CNLSA is utilized, as shown in Figs. 2(g, h) and 3 (g, h). With respect to k3, in spite of poor performance of the other two algorithms, the phase maps evaluated by CNLSA recover, as
Table 1 The maximal unwrapped phase errors at y¼0.14 mm in the wrapped phase maps of the OPD Λ23 evaluated by FT with HW, FT with RW and CNLSA, respectively; WN —wavenumber series; FT—Fourier transform; HW—hanning window; RW—rectangular window; A—absolute maximal phase error; R—relative maximal phase error. WN
k1 k2 k3
FT with HW
FT with RW
CNLSA
A(rad)
R (%)
A (rad)
R (%)
A (rad)
R (%)
88.19 100.31 108.93
321.7 365.9 397.4
1.60 1.99 52.69
5.8 7.3 192.2
– 0.28 2.11
– 1.0 7.7
61
shown in Figs. 2(i) and 3(i). Although there are ripples in the phase maps of CNLSA corresponding to k3, they are still acceptable. In order to compare with the differences among the three algorithms, the section lines at y¼0.14 mm in the unwrapped phase maps of Fig. 3 are chosen to be analyzed. If the unwrapped phase with respect to k1 that is evaluated by CNLSA is taken as the reference, the unwrapped phase differences between the CNLSA, the Fourier transform with Hanning window or rectangular window and the reference are given in Table 1. With respect to k3, the Fourier transform with rectangular window yields a huge phase error, while the CNLSA is still convergent and satisfactory, which unwrapped phase error is similar to that of the Fourier transform with rectangular window with respect to k1 and k2. Another experiment is performed as the series of the wavenumber scanning is k4, where k0 ¼7.308 106 m 1 and Δk¼ 4.503 103 m 1. The object measured is still compose of 3 smooth surfaces S1, S2 and S3, where S1 and S2 are the two surfaces of the optical wedge with the OPD Λ120 ¼ 9.06 mm at x¼0 mm, y¼0 mm and the title angle 60 . Surface S3 is of curvature and set to parallel to the surface S1, with the OPD Λ130 ¼7.69 mm at x¼ 0 mm, y¼0 mm. The interference forms a closed fringe pattern as shown in Fig. 4. Fig. 5(a–f) shows the wrapped and unwrapped phase maps of Λ13 evaluated by Fourier transform with Hanning window, Fourier transform with rectangular window and CNLSA. The depth resolutions of Fourier transform with Hanning window and rectangular window are 2.79 mm and 1.40 mm, respectively, which is larger than dΛ ¼|Λ120 Λ130|¼1.37 mm. Therefore, there are phase ripple errors in the phase maps evaluated by those two methods, as shown in Fig. 5(a, b, d, e). Among the three algorithms, the phase maps evaluated by CNLSA are smoothest, as shown in Fig. 5(c and f). According to the two group experimental results above, it is concluded that the CNLSA is the best among the algorithms for the DRWSI.
5. Simulation
Fig. 4. The closed interference fringe pattern.
The simulations were performed to take into account the optical setup that was described above: the interference signals were from the three surfaces S1, S2 and S3, with the corresponding reflective light intensities of I1 ¼0.58, I2 ¼0.50 and I3 ¼0.21; the initial phases were ϕ120 ¼ ϕ230 ¼ ϕ130 ¼ 0 rad; the wavenumber
Fig. 5. The wrapped and unwrapped phase maps of the OPD Λ13, for the wavenumber series k4; (a, b, c) are the wrapped phase maps evaluated by Fourier transform with Hanning window, Fourier transform with rectangular window and CNLSA, respectively; (d, e, f) are the unwrapped phase maps of (a, b, c), respectively.
62
Y. Zhang et al. / Optics and Lasers in Engineering 66 (2015) 58–63
series was k3. As Λ230 ¼ 3.82 mm, the parameter Λ120 was varied from 7.50 mm to 2.50 mm in steps of 0.005 mm, while the OPD and phase at the OPD Λ230 discrepancies between the preset values and those that were evaluated by a Fourier transform with the Hanning window or rectangular window and CNLSA were simulated. The spectral crosstalk between Λ120 and Λ230 yields the OPD and phase errors by the following mechanisms. (1) Peak drifting occurs as dΛ ¼ |Λ230 Λ120| approaches the depth resolution δΛ. The depth resolutions of Fourier transforms with a rectangular or Hanning window are 1.46 mm and 2.92 mm, respectively. As dΛ approaches its depth resolution from the right side of Fig. 6(a), all three algorithms produce fluctuating errors both in the OPD and phase values. In order to evaluate the fluctuating errors, the ranges of the OPD Λ120 scanning are chosen to be just outside the regions in which double peaks switch to a single peak for all the algorithms,
as shown in Table 2. It is noted that there are relatively larger OPD and phase's offset errors when a Fourier transform with rectangular window is used, because of its poor performance in restraining the sidelobes and the constant OPD between the OPD Λ230 and the DC signal. Among the three algorithms, the offset and standard deviation errors are minimal for the CNLSA, with the finest depth resolution, as shown in Table 2. (2) Peak overlapping occurs as dΛ become smaller than the depth resolution δΛ. Because the phase term 2∙k(0)∙Λpq in Eq. (1) is very sensitive to the Λpq, the double peaks of the OPD Λ230 and Λ120 rapidly switch between the double peaks and a single peak, as dΛ o0.94 mm or dΛ o1.58 mm, for a Fourier transform with rectangular window or Hanning window, respectively. As shown in Fig. 6(a), these cause rapid switching between small and large OPD and phase errors as a result of the automatic peak-searching algorithm. Due to the error in the initial values of X when automatic peak searching is used, the switching error or the depth resolution of the CNLSA occurs at dΛ ¼ δΛ ¼0.80 mm, as shown in Fig. 6. In the experiment described above, with respect to k3, dΛ ¼ 0.86 mm is smaller than the switching point of a Fourier transform with rectangular window (0.94 mm); thus, the wrapped phase map is smashed, as shown in Fig. 2(f). Because dΛ is larger than the switching point of CNLSA (0.80 mm), its wrapped phase map recovers, as shown in Fig. 2(i). Therefore, both the experimental and the simulated results show that CNLSA yields a good OPD and phase performance, which cannot be achieved by using window functions. It is noted that the more correct the initial values X are, the finer the depth resolution of the CNLSA becomes. If the initial values of the OPD Λ230 and Λ120 are given by the preset values in the simulation in spite of the double peaks to a single peak switching, the OPD Λ230 error is 0.07 mm, and the phase error is 0.25 rad (λ/50) at dΛ ¼ δΛ ¼ 0.34 mm, as shown in Fig. 6(b). In this case, the phase error 0.25 rad (λ/50) is similar to the CNLSA error obtained in the experiment with respect to k2, as shown in Table 1; however, the depth resolution is much finer than the experimental one.
6. Conclusion
Fig. 6. The OPD and the wrapped phase errors at the OPD Λ230, simulated by the three algorithms, respectively; green is a Fourier transform with Hanning window; blue is a Fourier transform with rectangular window; red is CNLSA, with respect to the wavenumber scanning series k3; Pha. Err.—the wrapped phase error; OPD Err.— the OPD error; (a) a Fourier transform with Hanning window or rectangular window and CNLSA, automatic peak searching; (b) CNLSA, as the initial values of the vector X are the preset values. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
The spectral leakage from the neighboring peaks is one of the principal error sources in DRWSI. We proposed a robust algorithm, called the Complex Number Least Squares Algorithm (CNLSA), for DRWSI. The algorithm yields a significant improvement of the depth resolution, and suppresses the OPD and phase ripple errors without requiring any increase in hardware complexity. Although the Levenberg–Marquardt algorithm that was used in the present study is a classical method, both the experiment and the simulation suggest that CNLSA performs much better than any window function. Theoretically, it is possible to obtain a better depth resolution by developing a global fitting algorithm that does not depend on the initial values. In conclusion, a CNLSA helps remove the limitation on the range of light source's wavenumber scanning, and is likely to become the principal data evaluating algorithm for DRWSI.
Table 2 The OPD Λ230 and phase errors simulated by FT with HW, FT with RW and CNLSA, respectively, with respect to the wavenumber series k3; std—standard deviation error. Error
OPD (mm) Phase (rad)
FT with HW Λ120 A(6.46 mm, 7.02 mm)
FT with RW Λ120 A (5.00 mm, 5.56 mm)
CNLSA Λ120 A (4.82 mm, 5.38 mm)
Offset
Offset
Offset
7 8.05 10 7 0.03
Std 3
7 4.66 10 7 0.02
3
7 0.14 7 0.70
Std 70.08 70.35
7 3.06 10 7 0.01
Std 3
7 2.50 10 3 7 0.01
Y. Zhang et al. / Optics and Lasers in Engineering 66 (2015) 58–63
Acknowledgment The authors wish to thank the National Natural Science Foundation of China (NSFC) and Provincial Natural Science Foundation of Guangdong (NSFG) for their financial support with the Grants 11072063 and S2012010010327, respectively. References [1] De la Torre-Ibarra MH, Ruiz PD, Huntley JM. Simultaneous measurement of inplane and out-of-plane displacement fields in scattering media using phasecontrast spectral optical coherence tomography. Opt Lett 2009;34(6):806–8. [2] Davila A, Huntley JM, Pallikarakis C, Ruiz PD, Coupland JM. Simultaneous wavenumber measurement and coherence detection using temporal phase unwrapping. Appl Opt 2012;51(5):558–67. [3] Davila A, Huntley JM, Pallikarakis C, Ruiz PD, Coupland JM. Wavelength scanning interferometry using a Ti: sapphire laser with wide tuning range. Opt Laser Eng 2012;50(8):1089–96.
63
[4] Chakraborty S, Ruiz PD. Measurement of all orthogonal components of displacement in the volume of scattering materials using wavelength scanning interferometry. J Opt Soc Am A 2012;29(9):1776–85. [5] Xu J, Liu Y, Dong B, Bai Y, Hu L, Shi C, et al. Improvement of the depth resolution in depth-resolved wavenumber-scanning interferometry using multiple uncorrelated wavenumber bands. Appl Opt 2013;52(20):4890–7. [6] Hibino K, Hanayama R, Burke J, Oreb BF. Tunable phase-extraction formulae for simultaneous shape measurement of multiple surfaces with wavelengthshifting interferometry. Opt Express 2004;12(23):5579–94. [7] Greivenkamp JE. Generalized data reduction for heterodyne interferometry. Opt Eng 1984;23:350–2. [8] Wang Q, Ning YN, Grattan KTV, Palmer AW. A curve fitting signal processing scheme for a white-light interferometric system with a synthetic source. Opt. Laser. Technol 1997;29(7):371–6. [9] Sheppard RJ, Jordan BP, Grant EH. Least squares analysis of complex data with application to permittivity measurements. J Phys D: Appl Phys 1970;3 (11):1759–64. [10] Buckland JR, Huntely JM, Turner SRE. Unwrapping noisy phase maps by use of a minimum-cost-matching algorithm. Appl Opt 1995;34(23):5100–8.