Multiwavelength phase unwrapping and aberration ... - OSA Publishing

5 downloads 0 Views 402KB Size Report
aberration correction using depth filtered digital holography. Volker Jaedicke,1,* Sebastian Goebel,1 Nektarios Koukourakis,2 Nils C. Gerhardt,1. Hubert Welp,3 ...
4160

OPTICS LETTERS / Vol. 39, No. 14 / July 15, 2014

Multiwavelength phase unwrapping and aberration correction using depth filtered digital holography Volker Jaedicke,1,* Sebastian Goebel,1 Nektarios Koukourakis,2 Nils C. Gerhardt,1 Hubert Welp,3 and Martin R. Hofmann1 1 2

Photonics and Terahertz Technology, Ruhr-Universität Bochum, Universitätsstr 150, 44780 Bochum, Germany

Laboratory for Measurement and Testing Techniques, Technische Universität Dresden, Helmholtzstr. 18, 01062 Dresden, Germany 3

Department of Electrical Engineering and Information Technology, University of Applied Science Georg Agricola, Herner Str. 45, 44789 Bochum, Germany *Corresponding author: [email protected] Received April 2, 2014; revised May 13, 2014; accepted June 9, 2014; posted June 10, 2014 (Doc. ID 209430); published July 9, 2014 In this Letter, we present a new approach to processing data from a standard spectral domain optical coherence tomography (OCT) system using depth filtered digital holography (DFDH). Intensity-based OCT processing has an axial resolution of the order of a few micrometers. When the phase information is used to obtain optical path length differences, subwavelength accuracy can be achieved, but this limits the resolvable step heights to half of the wavelength of the system. Thus there is a metrology gap between phase- and intensity-based methods. Our concept addresses this metrology gap by combining DFHD with multiwavelength phase unwrapping. Additionally, the measurements are corrected for aberrations. Here, we present proof of concept measurements of a structured semiconductor sample. © 2014 Optical Society of America OCIS codes: (090.1760) Computer holography; (090.0090) Holography; (090.1995) Digital holography; (100.3175) Interferometric imaging. http://dx.doi.org/10.1364/OL.39.004160

In this Letter, we present a new approach to processing spectral domain optical coherence tomography (SDOCT) data. To obtain tomographic images in SD-OCT, spectra are point-wise recorded in an interferometric setup. These spectra are Fourier transformed to obtain complex valued depth profiles. In intensity-based imaging, cross-sectional or enfaced images are extracted from these data by taking the absolute value. The axial resolution is determined by the coherence length of the system, which is typically a few microns [1]. In addition, the phase can be used to get height maps of the sample’s surface or interface layers. The accuracies of these height maps are limited by the phase noise of the system (mostly governed by the stability of the interferometer) and can reach subnanometer precision in common path setups [2]. The maximum step height of the sample that can be measured is set by the phase unwrapping problem, e.g., the phase is only known as modulo 2π. This usually limits the maximum resolvable step height of a sample in reflection geometry to half of the wavelength of the system. Therefore, there is a gap for resolvable sample structures, which is determined by the phase unwrapping problem and by the axial resolution of the system in intensity-based imaging. In this Letter, we demonstrate phase imaging of a sample that can be resolved neither with common intensity-based processing nor by regular phase resolved processing using a standard SD-OCT system. Instead, we combined multiple phase maps obtained by depth filtered digital holography (DFDH) [3] with a multiwavelength phase unwrapping algorithm [4]. We corrected the measurement for aberrations in a semi-automated procedure [5]. Previously we used DFDH in a full-field OCT setup in off-axis geometry and filtered out spurious reflections of 0146-9592/14/144160-04$15.00/0

a multilayer sample to obtain a height profile using one phase map [3]. In the original work, aberrations of the measurements were not taken into account. In this Letter, we used a standard scanning SD-OCT system to acquire synthesized holograms in a point-wise manner. In general, the processing scheme of DFDH is similar to spectroscopic OCT. In spectroscopic OCT, the absolute values of depth resolved complex spectra are used to obtain spectroscopic information of biological tissue [6]. In contrast, in DFDH the depth resolved phase data is used to calculate lateral phase maps. For a single backscattering interface, the recorded spectra in SD-OCT can be described by Ikx;y  I S kx;y  I R kx;y q     2 I S kx;y I R kx;y cosk · Δzx;y ;

(1)

where I S k is the signal backscattered by the sample and I R k is the signal of the reference arm (DC terms). The last term is the interference of both signals that leads to the characteristic modulation of the spectra by the path length difference Δzx;y , which includes subwavelength variations of the sample. In detail, the DFDH algorithm works as shown in Fig. 1. Step 1: the interference term is isolated by subtracting the DC component, and the spectra are resampled from wavelength (λ) to wavenumber (k) domain. Step 2: an inverse Fourier transform is applied to each spectrum to get complex valued depth profiles. Step 3: the depth of interest (DOI) is filtered by multiplying a window function (i.e., Gaussian window). Step 4: the filtered depth profile is Fourier transformed to again obtain a complex spectrum, which is resampled from wavenumber to © 2014 Optical Society of America

July 15, 2014 / Vol. 39, No. 14 / OPTICS LETTERS

zn  z0n 

Fig. 1.

DFDH signal processing.

wavelength domain. From these depth filtered spectra, lateral phase maps for each discrete wavelength are extracted. The phase maps can be used to calculate a height map for each wavelength using the following equation: zx; yjλ 

λ · φx; yjλ : 4π

(2)

The individual height maps obtained by DFDH may suffer from the phase unwrapping problem. Several methods for obtaining the true, unwrapped phase exist; these can be categorized into numerical and multiwavelength methods [7]. Numerical phase unwrapping algorithms work in an iterative manner by estimating the phase jumps when the phase is wrapped. Nevertheless, in several cases the algorithm will fail due to the sample structure, undersampling, or noise. In contrast, in multiwavelength unwrapping, two or more phase maps at different wavelengths are used to obtain an unambiguous phase map. This concept has been previously used in interferometry [8] and has been adapted for digital holography [4] and OCT [9]. Different algorithms for using more than two phase maps at different wavelengths exist [4,7]. We used hierarchical unwrapping [4], which is wellsuited for the DFDH data due to the ability to use a large number of wavelengths. This method works in an iterative manner with four steps within each of the n iterations. For n iterations, n  1 phase maps are used. In step 1, a new phase map, Φn , is obtained by calculating the difference of two phase maps at different wavelengths and adding 2π when the result is negative, as described in the following equation:  Φn 

φ0 − φn φ0 − φn  2π

if φ0 ≥ φn : if φ0 < φn

(3)

In step 2, the new phase map corresponds to a synthetic wavelength, which can be calculated by the following equation: Λn 

λ0 λn : jλ0 − λn j

(4)

In step 3, the corresponding height map, zn , is calculated by Eq. (2) using the synthetic wavelength, Λn . In step 4 (added from the second iteration on), the height map of the previous iteration step is used as a coarse map to unwrap the height map of the actual iteration step, as described in the following equation:

  Λn z − z0n floor 2 n−1  0.5 : 2 Λn

4161

(5)

Here, z0n is the original height map of the actual iteration, which might be ambiguous. The synthetic wavelength of the first iteration has to be large enough to cover the maximum height difference of the measurement. Since the synthetic wavelength is directly correlated to the noise [Eq. (2)], it is decreased in each iteration step. The last step is added to make sure that the height profile is properly unwrapped with decreasing synthetic wavelengths. This concept can be expanded to an arbitrary number of iteration steps when the noise doesn’t exceed half of the previous synthetic wavelength, i.e., Λn ∕2 > εn−1 . Here Λn is the synthetic wavelength of the nth iteration step, while εn−1 is the noise of the previous height profile. In OCT and DFDH, phase measurements are obscured by aberrations. This is especially problematic when point-wise scanning systems are used. We adapted an approach of previous authors [5], which uses flat areas of the sample and a fit of Zernike polynomials to correct the measurements. Zernike polynomials offer a straightforward and widely used approach to describe and quantify optical aberrations. In detail, the correction algorithm works as follows: first, flat areas of the unwrapped height map (zu ) are isolated (e.g., the zones outside the flat areas are set to zero). Then N coefficients (ai ) of the Zernike polynomials (Z i ) are fitted using Eq. (6) (when zu , is given in polar coordinates as a function of ρ, ϑ): zu ρ; ϑ 

iN X i0

ai Z i ρ; ϑ:

(6)

The coefficients in turn are used to calculate a numerical lens [Eq. (7)] with unit amplitude. The numerical lens is converted to Cartesian coordinates and eventually multiplied with the original complex data for correction, as described in the following equation:  ΓZ ρ; ϑ  exp

−i

N X i0

 ai Z i

 4π : λ

(7)

Our complete algorithm has three stages and combines the three concepts described in the previous paragraphs: first, multiple phase maps are obtained by DFDH; second, the height maps are corrected for aberrations; and third, the final, unwrapped height profile is calculated. The algorithm has a comparable processing time to conventional OCT phase resolved processing with numerical unwrapping. In the first stage, phase maps are obtained by the DFDH algorithm. In DFDH, the window size of the DOI sets the spectral sampling and the separation of the wavelengths of the phase maps. The smallest separation between two consecutive wavelengths in turn sets the maximum achievable synthetic wavelength [Eq. (4)]. The spectral sampling and separation can be reduced by interpolation due to zero padding, but this may introduce additional noise due to highly correlated phase maps. In the second stage, the sample

4162

OPTICS LETTERS / Vol. 39, No. 14 / July 15, 2014

measurement is hierarchically unwrapped using the phase maps obtained by DFDH. Here one has to determine the appropriate set of wavelengths of the phase maps for the hierarchical unwrapping algorithm. This includes the starting wavelength, the step size between the wavelengths, and the last wavelength. The maximum synthetic wavelength that is necessary to cover the maximum height variations of the measurement depends not only on the sample’s structures but also on the aberrations introduced by the system (i.e., tilt between the sample and reference arm). Since noise is the most critical factor for the algorithm, the appropriate phase maps have to be carefully selected. For aberration correction, flat areas of the sample are used to fit Zernike polynomials and to calculate the numerical lens. Eventually all phase maps, φx; yjλ , are corrected with this numerical lens. In the third stage, the final height map of the sample is obtained by hierarchical unwrapping of the corrected phase maps. The OCT system we used is a commercially available Thorlabs Callisto SD-OCT system with 100 nm, 3 dB spectral bandwidth and a central wavelength of 930 nm. The lateral and axial resolutions were 8 and 7 μm, respectively. The beam was scanned in the lateral dimension by two galvanometric-driven mirrors, while at each position a spectrum was recorded with a 1024-pixel line scan camera with a 1 kHz frame rate. The spectral sampling set the total available depth range to 1.7 mm in air and the depth sampling to 3.32 μm. First we analyzed the noise in our system by measuring a flat mirror substrate. The substrate had a roughness of less than 1 Å and flatness of λ∕10 at 632.8 nm. We recorded a data set of the mirror substrate with 512 × 512 data points with 1 mm lateral dimensions. First the data were processed with the DFDH algorithm. The window width for the DOI was set to 64 points (212.5 μm). We obtained 64 phase maps at wavelengths ranging from 868 to 992 nm that were numerically unwrapped. For this purpose we tried different approaches and selected Flynn’s method, which was implemented using the code provided by the book by Ghiglia and Pritt [10]. The aberrations of the measurement were corrected (8 Zernike coefficients) as described above using the full area of the substrate. We then calculated the individual noise of each of the 64 phase maps and took the median value. The remaining root mean square (rms) error was 0.088 μm. To compare the noise behavior between hierarchical and numerical unwrapping, we further calculated the rms noise for an increasing number of iterations using hierarchical unwrapping. As expected, the noise dropped continuously for an increasing number of iteration steps (minimum of 0.1 μm for eight phase maps). This is because the noise decreased with decreasing synthetic wavelengths. After falling to this minimum, the noise increased again, which can be explained by two mechanisms: spikes occur due to pixels with high noise, and phase jumps and thus discontinuities in the height maps occur. The latter is the case when the whole profile cannot be unwrapped properly. This analysis shows that the hierarchical unwrapping algorithm is suited for unwrapping without considerably increasing the noise. Furthermore, due to the decreasing noise when multiple

Fig. 2. Top: microscope image of the GaAs semiconductor sample. Bottom: profilometer scans obtained with stylus device.

iterations are used, hierarchical unwrapping can be superior to unwrapping based on only two phase maps. In semiconductor metrology, nanometer accuracy across a micrometer range is mandatory. Therefore, a semiconductor sample was chosen for a proof of concept measurement. The selected GaAs semiconductor sample was produced by molecular beam epitaxy and was surface structured by etching. Thus its dimensions are well-known from the process parameters. In addition, the heights of the structures were verified using profilometer scans from a stylus device (bottom of Fig. 2). The blue line and the red line in the microscope image in top of Fig. 2 indicate the positions where the line scans were obtained. A data set of the semiconductor sample was recorded with the same parameters as the mirror substrate. First we used conventional OCT-based phase imaging and numerically unwrapped the phase map obtained from the complex OCT depth profiles using the depth of the highest intensity. The result is shown in Fig. 3. The x and y profiles were taken across the same positions as in Fig. 2. As expected, the unwrapping failed, since the structures of the sample were far beyond the fundamental limit of half of the wavelength of the system. In contrast, we applied our three-stage algorithm. First we used the DFDH algorithm and calculated 100 phase

Fig. 3. Numerically unwrapped height profiles of the semiconductor sample.

July 15, 2014 / Vol. 39, No. 14 / OPTICS LETTERS

4163

problem in phase-based imaging and by the resolution in intensity-based imaging. By means of DFDH, we obtain depth filtered phase maps at different wavelengths. This in turn enables the use of multiwavelength phase unwrapping techniques. In contrast to the approach of previous authors [9], our concept uses more than two phase maps for unwrapping. Using more than two phase maps can be beneficial in terms of noise, as we have shown in our analysis. Furthermore, the use of numerical lenses and Zernike polynomials allows us to correct the aberrations introduced by the system. We demonstrated the benefit of our concept by analyzing the topological structure of an etched semiconductor sample. The structures of the sample are far beyond half of the maximum wavelength and well below the axial resolution of the system used. Therefore, they cannot be resolved by conventional phase or intensity imaging, but can be resolved with our approach. Due to the depth filtering, the proposed technique can also be used for measurements in multilayer semitransparent samples [3,11]. As a next step, we want to develop an automated algorithm that takes the noise of the measurement and the information provided by intensity into account. So far the different parameters for the algorithm have to be found in a heuristic search. In addition, more advanced filter methods [12] and improved system design, such as a common path, can improve the results. Fig. 4. Height map and profiles of the hierarchically unwrapped semiconductor sample measurement.

maps of a DOI with a width of 199 μm (including zero padding of 40 points). Second, a height map for aberration correction was calculated by using hierarchical unwrapping, employing four phase maps (wavelengths of 942, 943, 944, and 945 nm). Therefore, the largest and smallest synthetic wavelengths were 887.36 and 296.73 μm, respectively. All 100 phase maps were subsequently corrected with a numerical lens (8 Zernike coefficients), which was calculated from flat areas of the sample between the structures. Finally an unambiguous height profile was obtained from the corrected phase maps by hierarchical unwrapping with a new set of wavelengths (967, 973, 980, and 986 nm). This set the smallest synthetic wavelength to 50.695 μm, which was considerably smaller than the previously used smallest synthetic wavelength and thus reduced the noise. The height map of the sample and profile plots are shown in Fig. 4. These were low pass filtered to suppress spike noise. In contrast to the conventional processed phase data, the structures of the sample can be correctly measured. As a drawback, the results suffer from worse lateral resolution than the profilometer scans. In addition, some unwrapping errors can be seen in noisy parts of the measurement, and the y profile is slightly tilted. In conclusion, we demonstrated a three-stage algorithm based on DFDH to obtain topographic height maps of surfaces or interfaces using a standard SD-OCT system. Our approach addresses the metrology gap of resolvable sample structures set by the unwrapping

We acknowledge Arne Ludwig for providing the semiconductor sample, Felix Mitschker for the help with the profilometer scans, and the Research School Plus, the RWTÜV Stiftung, and the European Space Agency ESA for funding. References 1. J. Walther, M. Gaertner, P. Cimalla, A. Burkhardt, L. Kirsten, S. Meissner, and E. Koch, Anal. Bioanal. Chem. 400, 2721 (2011). 2. M. Choma, A. K. Ellerbee, C. Yang, T. L. Creazzo, and J. A. Izatt, Opt. Lett. 30, 1162 (2005). 3. N. Koukourakis, V. Jaedicke, A. Adinda-Ougba, S. Goebel, H. Wiethoff, H. Höpfner, N. C. Gerhardt, and M. R. Hofmann, Opt. Express 20, 22636 (2012). 4. C. Wagner, W. Osten, and S. Seebacher, Opt. Eng. 39, 79 (2000). 5. T. Colomb, E. Cuche, F. Charrière, J. Kühn, N. Aspert, F. Montfort, P. Marquet, and C. Depeursinge, Appl. Opt. 45, 851 (2006). 6. V. Jaedicke, S. Agcaer, F. E. Robles, M. Steinert, D. B. Jones, S. Goebel, N. C. Gerhardt, H. Welp, and M. R. Hofmann, Biomed. Opt. Express 406, 35 (2013). 7. M. K. Kim, SPIE Rev. 1, 018005 (2010). 8. C. R. Tilford, Appl. Opt. 16, 1857 (1977). 9. H. C. Hendargo, M. Zhao, N. Shepherd, and J. A. Izatt, Opt. Express 17, 5039 (2009). 10. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software (Wiley, 1998). 11. S. Goebel, V. Jaedicke, N. Koukourakis, H. Wiethoff, A. Adinda-Ougba, N. C. Gerhardt, H. Welp, and M. R. Hofmann, Proc. SPIE 8589, 85891J (2013). 12. R. Seara, A. A. Go Alves, and P. B. Uliana, Appl. Opt. 37, 2046 (1998).