Application of maximum likelihood estimator in nano-scale optical path length measurement using spectral-domain optical coherence phase microscopy

Spectral-domain optical coherence phase microscopy (SDOCPM) measures minute phase changes in transparent biological specimens using a common path interferometer and a spectrometer based optical coherence tomography system. The Fourier transform of the acquired interference spectrum in spectral-domain optical coherence tomography (SD-OCT) is complex and the phase is affected by contributions from inherent random noise. To reduce this phase noise, knowledge of the probability density function (PDF) of data becomes essential. In the present work, the intensity and phase PDFs of the complex interference signal are theoretically derived and the optical path length (OPL) PDF is experimentally validated. The full knowledge of the PDFs is exploited for optimal estimation (Maximum Likelihood estimation) of the intensity, phase, and signal-to-noise ratio (SNR) in SD-OCPM. Maximum likelihood (ML) estimates of the intensity, SNR, and OPL images are presented for two different scan modes using Bovine Pulmonary Artery Endothelial (BPAE) cells. To investigate the phase accuracy of SD-OCPM, we experimentally calculate and compare the cumulative distribution functions (CDFs) of the OPL standard deviation and the square root of the CramérRao lower bound ( SNR 2 / 1 ) over 100 BPAE images for two different scan modes. The correction to the OPL measurement by applying ML estimation to SD-OCPM for BPAE cells is demonstrated. ©2008 Optical Society of America OCIS codes: (110.4500) Optical coherence tomography; (110.0180) Microscopy; (180.3170) Interference microscopy, (170.0110) Medical optics and biotechnology References and links 1. K. Svoboda, C. F. Schmidt, B. J. Schnapp, and S. M. Block, “Direct observation of kinesin stepping by optical trapping interferometry,”Nature 365, 721-727 (1993). 2. J. Farinas and A. S. Verkman, “Cell volume and plasma membrane osmotic water permeability in epithelial cell layers measured by interferometry,” Biophys. J. 71, 3511-3522 (1996). 3. C. Yang, A. Wax, M. S. Hahn, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Phase-referenced interferometer with subwavelength and subhertz sensitivity applied to the study of cell membrane dynamics,” Opt. Lett. 26, 1271-1273 (2001). 4. A. Barty, K. A. Nugent, D. Paganin, and A. Roberts, “Quantitative optical phase microscopy,” Opt. Lett. 23, 817-819 (1998). 5. P. Marquet, B. Rappaz, P. J. Magistretti, E. Cuche, Y. Emery, T. Colomb, and C. Depeursinge, “Digital holographic microscopy: a noninvasive contrast imaging technique allowing quantitative visualization of living cells with subwavelength axial accuracy,” Opt. Lett. 30, 468-470 (2005). 6. G. Popescu, L. P. Deflores, J. C. Vaughan, K. Badizadegan, H. Iwai, R. R. Dasari, and M. S. Feld, “Fourier phase microscopy for investigation of biological structures and dynamics,” Opt. Lett. 29, 2503-2505 (2004). (C) 2008 OSA 27 October 2008 / Vol. 16, No. 22 / OPTICS EXPRESS 17186 #95322 $15.00 USD Received 22 Apr 2008; revised 20 Aug 2008; accepted 6 Oct 2008; published 13 Oct 2008 7. S. Kostianovski, S. G. Lipson, and E. N. Ribak, “Interference microscopy and Fourier fringe analysis applied to measuring the spatial refractive-index distribution,” Appl. Opt. 32, 4744(1993). 8. T. Ikeda, G. Popescu, R. R. Dasari, and M. S. Feld, “Hilbert phase microscopy for investigating fast dynamics in transparent systems,” Opt. Lett. 30, 1165-1167 (2005). 9. M. A. Choma, A. K. Ellerbee, C. Yang, T. L. Creazzo, and J. A. Izatt, “Spectral-domain phase microscopy,”Opt. Lett. 30, 1162-1164 (2005). 10. C. Joo, T. Akkin, B. Cense, B. H. Park, and J. F. de Boer, “Spectral-domain optical coherence phase microscopy for quantitative phase-contrast imaging,” Opt. Lett. 30, 2131-2133 (2005). 11. D. C. Adler, R. Huber, and J. G. Fujimoto, “Phase-sensitive optical coherence tomography at up to 370,000 lines per second using buffered Fourier domain mode-locked lasers,” Opt. Lett. 32, 626-628 (2007). 12. B. White, M. Pierce, N. Nassif, B. Cense, B. Park, G. Tearney, B. Bouma, T. Chen, and J. de Boer, “In vivo dynamic human retinal blood flow imaging using ultra-high-speed spectral domain optical coherence tomography,” Opt. Express 11, 3490-3497 (2003). 13. S. Yazdanfar, C. Yang, M. Sarunic, and J. Izatt, “Frequency estimation precision in Doppler optical coherence tomography using the Cramer-Rao lower bound,” Opt. Express 13, 410-416 (2005). 14. B. Park, M. C. Pierce, B. Cense, S. H. Yun, M. Mujat, G. Tearney, B. Bouma, and J. de Boer, “Real-time fiber-based multi-functional spectral-domain optical coherence tomography at 1.3 μm,” Opt. Express 13, 3931-3944 (2005). 15. A. Papoulis and S. U. Pillai, Probability, Random Variables and Stochastic Processes, (McGraw-Hill, 2002) 4th edition . 16. A. J. Miller and P. M. Joseph, “The use of power images to perform quantitative analysis on low SNR MR images,” Magn. Reson. Imaging 11, 1051-1056 (1993). 17. G. McGibney and M. R. Smith, “An unbiased signal-to-noise ratio measure for magnetic resonance images,” Med Phys. 20, 1077-1078 (1993). 18. B. A. van den, Handbook of Measurement Science, (Wiley, Chichester, England, 1982) Vol. 1, pp. 331-377.

The recent application of SD-OCT to phase measurement has resulted in significant improvements in phase stability, sensitivity, and speed compared with those of time-domain OCT based systems [12].In many applications, however, especially those with low signal to noise ratio (SNR), the phase sensitivity decreases, making it difficult to measure nanometerscale path length and refractive index differences that are required to characterize organelle structure and function.
In this paper, we formulate a theory for the probability distribution function (PDF) for the phase and intensity in spectral-domain optical coherence phase microscopy (SD-OCPM) [10] and demonstrate a good agreement between the theoretical and experimental PDFs.Moreover, we theoretically and experimentally depict the phase sensitivity of SD-OCPM as a function of SNR.Previously, the fundamental uncertainty limits on frequency/phase estimation precision in Doppler-OCT/OCT in the case of additive noise have been reported by using either the Cramér-Rao lower bound (CRLB) or phasor noise analysis [13][14].We show that the phase sensitivity approaches the square root of CRLB at high SNR; however, the square root of the CRLB is not valid for predicting the phase sensitivity either at low SNR or for an optical path length (OPL) equal to an integer number of half the center wavelength.In addition, we have developed a maximum likelihood (ML) estimator for optimal estimation of phase, intensity, and SNR in SD-OCT.We show via simulation that the ML estimator outperforms the conventional mean estimator in terms of phase precision.We present ML estimated Bovine Pulmonary Artery Endothelial (BPAE) cell intensity, SNR, and OPL images for two different scan modes.To investigate phase precision of our SD-OCPM using two different scan modes, the cumulative distribution functions (CDFs) of OPL standard deviation and the square root of the CRLB over 100 images are calculated and compared.Finally, we validate our proposed ML estimator by acquiring 100 quantitative phase contrast images of a BPAE cell using SD-OCPM for two scan modes and show the improved measured phase by the ML estimator.

Probability density functions of intensity and phase
Fourier transform of the interference spectrum received from a SD-OCPM system is intrinsically complex-valued, and represents the magnitude (ρ) and phase (θ) at a certain location and time.The actual complex value (S) can be corrupted by the addition of a white (circular symmetric) complex Gaussian variable (N), giving the measured complex value X expressed as, with ρ and θ the magnitude and phase of the actual complex value S. The joint probability distribution function of the real parts, X R (S R ,σ), and the imaginary part, X I (S I ,σ), of the corrupted data is given by, with σ characterizing the magnitude of the noise.To find the PDFs of the amplitude and phase of the corrupted data, we define the real and imaginary parts in the polar coordinate system ) cos(φ By using a bivariate transformation, Eq. ( 4) can be written as The marginal PDF of the intensity (I=R 2 ) is given by [15], The marginal PDF of the measured signal phase (0≤φ<2π) is extracted by, .) 2 By defining the random variable u=(R-ρ cos(φ-θ))/σ, Eq. ( 9) can be written as (11) The marginal PDF of the phase is given by, where the signal to noise ratio is defined as SNR=ρ 2 /(2σ 2 ), and I 0 (.), and Q(.) are the zeroth order modified Bessel function and Q function, respectively.Using upper and lower bounds on the Q function [15], we can show the marginal PDF of the phase is bounded by two analytical expressions as follows: When (φ-θ) is equal to π/2 or3π/2, the marginal PDF of the phase is equal to (exp (-SNR))/2π.Using Eq. ( 12), the marginal PDF of the phase approaches a uniform distribution over [0, 2π] for low SNR values.For high SNR values, the derived inequalities (13)(14) show that the marginal PDF of the phase can be approximated as follows: .
Figure 1(a)-(b) shows the derived PDF of the corrupted phase (φ) in Eq. ( 12) for different true phase values (θ) at two different SNR values.As shown in Figs.1(a)-1(b), when either SNR decreases or the true phase value approaches zero or 2π, the PDF broadens due to the branch cut at 0 modulo 2π and causes an increase in the measured phase variance (or standard deviation).At a SNR equal to 20 dB (CRLB=0.005), the calculated phase variances using Eq. ( 12) were 0.005 and 8.35 rad 2 for the true phase values π/4 and π/100, respectively.When the SNR was decreased to 6 dB (CRLB=0.126), the phase variances were calculated to be 0.74 and 4.97 rad 2 for the true phase values π/4 and π/10, respectively.In this scenario, there is a significant difference between the exact phase sensitivity (standard deviation) and the reported fundamental uncertainty limits on the estimated phase ( ) using the CRLB [13].

Maximum likelihood estimator
The most intuitive way of estimating data is mean estimation without a priori knowledge of the data PDF.While several estimation methods have been proposed [15][16], ML estimators are known to be consistent and asymptotically efficient.In addition, if there exists an unbiased estimator of which the variance attains the CRLB, it is given by the ML estimator [17].
Considering a set of N independent, Gaussian distributed, complex data points, the joint PDF of the complex data, p c is simply the product of the PDFs of these data points: The resulting PDF is called the Likelihood function L [18].The ML estimate of each parameter is found by maximizing L with respect to that parameter [18].To simplify maximizing L, we may work with the natural logarithm of p c , which is a monotonic function.Thus, The ML estimates of amplitude (intensity), phase, and variance are found by maximizing Log L, with respect to these variables.At the maximum, the first order derivatives of Log L with respect to these variables are zero.From the resulting equations, ML estimators of the amplitude, phase, and variance are found to be: It can be seen that the ML estimates in Eqs.(18)(19)(20)(21) are given by the vector summation of the N measurements in the complex plane.

Experimental and simulation results
Figure 2 shows a schematic diagram of the SD-OCPM system employed in the experiment.A broadband 800 nm (λ 0 ) Kerr-lens mode-locked laser (~130 nm in FWHM, FemtoLasers, Austria) was employed for high resolution OCPM.The OCPM system was described in detail in Ref [10].Briefly, the OCPM beam passed through the XY beam scanners, and was introduced into the inverted microscope (Axiovert 200, Carl Zeiss) through its back port.The beam was then magnified by a telescope composed of the scan and tube lenses, and delivered to the specimen through the microscope objective (NeoFluar 20×, NA: 0.5, Carl Zeiss).For OCPM imaging, the reflection from the bottom surface of a coverslip served as the reference, whereas the back-scattered waves from the focal volume inside the specimen were the measurement fields.The back-scattered beams were reflected by the dichroic mirror, and coupled back to the fiber for the interference spectrum measurement.
The interference signals were detected by a high-speed spectrometer, where A-line rates were set to 5 and 10 kHz with different duty cycles (integration times).The measured spectrum was then Fourier-transformed to obtain complex data (intensity and quantitative phase information) at the focal location.SD-OCPM generated the structural and phase images of the specimen, as the beam scanned the specimen with the galvanometer XY scanners and the piezo-electric transducer (PZT).The OPL, optical path length difference between the reference and sample beams, is given by: ) , ( 4) , ( 0 The lateral resolution of SD-OCPM was measured as 0.75 μm at the full width at half maximum (FWHM), and the axial resolution, which is the combination of confocal and coherence gatings, was found as 2.58 μm through the measurement of the intensity at the corresponding depth as a mirror surface moves along the optical axis.
To show and validate the PDF of the OPL (phase) and the exact OPL sensitivity for each SNR value, we acquired 655360 data points by illuminating a stationary cover slip at an Aline rate of 10 kHz.The integration time (τ) of CCD camera per A-line was set at different values to control the SNR. Figure 3(a) shows a shifted histogram of the measured OPLs and theoretically derived PDF (Eq.( 12)) at SNR~58 dB (τ=25 μs).While the mean value of the measured OPLs was 222.396 nm, the histogram of the measured OPLs and the derived PDF were offset to zero nm. Figure 3(b) shows the experimental and theoretical results for the OPL standard deviation as a function of the measured SNR for OPL~222 nm. Figure 3(c) shows the theoretically predicted standard deviation of the OPL as a function of SNR for different OPLs.It is clear the sensitivity degrades when SNR decreases.In addition, when the true value of the OPL approaches an integer number of half the center wavelength (branch cuts= {0 nm, 400 nm,…}) the OPL standard deviation increases.At very small OPL values (for example OPL=3 nm), the standard deviation increased with increasing SNR, which can be explained by a PDF of the phase that became narrower around 0 and 2π with a mean phase value of approximately π, as depicted in Figs.1(a)-1(b).The standard deviation gradually decreased again when the PDF of the phase became asymmetric with a mean shifting towards zero when the SNR increased further.Figure 3(d) shows that the OPL standard deviation deviates from the square root of the CRLB at low SNR and OPL values close to an integer number of half the center wavelength.In statistics, the precision of an estimator is determined by the spread of the estimates when the experiment is repeated under identical conditions, and is represented by the standard deviation or the variance of the estimator.We considered the root mean-squared error (RMSE=√MSE) as a performance criterion.To find the precision of the ML and mean estimators, we simulated and calculated the RMSE of the two estimators for several phase mesurements.The sample data points (simulated OPL values) were divided into M sets of N data points, and we applied the ML and mean estimators to each set of data.Finally, we calculated the RMSE of each estimator by using the estimated OPLs for each set and the true OPL value for at least 65533 datasets.
The simulation results in Fig. 4 show the precision of the ML and mean estimators as a function of data points (N).At SNR=10 dB, true OPL value of 25 nm (θ=π/8), and N=100, the RMSE of the ML and mean estimators were 1.4 nm and 18.7 nm, respectively.The simulation results show that the ML method outperforms the mean method.For example, the precision of the ML and mean estimators were 2.2 nm and 58.2 nm at SNR=6 dB, true OPL value of 25 nm, and N=100, respectively.

Images
We assessed the performance of the ML estimator by imaging a prepared BPAE cell using the experimental set up depicted in Fig. 2. The cells were fixed and did not exhibit any motion during the measurements.The system magnification and scanning speed gave an image size 256×256 pixels over a 172 μm×172 μm field of view (FOV).In order to acquire several images, we employed two different scanning modes, which will be referred to as the MB-scan and the BM scan.In OCT, an A-scan refers to a single reflectivity depth scan in the sample.A B-scan is the 2D cross sectional OCT image constructed from multiple A-scans taken over different transverse positions.An M-scan refers to multiple A-scans taken over time at the same transverse location.Using the same terminology, an MB-scan is a term used to describe multiple M-scans taken over different transverse positions to create a 2D OCT image of M scans.A BM-scan on the other hand is multiple B-scans taken over time for the same transverse scan region.In the MB-scan, the beam was positioned at a desired location of the specimen, and 120 spectra were acquired for 24 ms at an A-line rate of 5 kHz.Since the settling time of the scanners was ~2 ms, the latter 100 out of 120 spectra were chosen to avoid the effect of scanner settling dynamics on the phase measurement.In the MB-scan mode, the total acquisition time was 26 min for an image size of 256×256 pixels at 5 kHz A-line rate.In the BM-scan, the total acquisition time was 22 min for the same MB-scan image size and 5 kHz A-line rate.To examine the spatial phase stability for the MB and BM scan modes, we obtained 100 phase images of a coverslip and measured the standard deviation across a field of view of 43μm×43μm with 70×70 pixels.The measured standard deviations across the field of view were measured to be 1.7 nm and 2.0 nm in air at an SNR~30 dB for the MB and BM scan modes, respectively.The measured standard deviation increases in BM scan mode, which is due to reliance on the repeatability of the scanner pattern.) and the OPL standard deviation at each pixel were calculated using the ML estimated SNR and OPL in Eq. ( 21) and Eqs.(19, 22), respectively.
Figures 6(b)-6(c) shows CDFs for two different modes.The probability of an OPL standard deviation less than 2 nm was ~0.9, which was 0.05 less than the probability of √CRLB in the MB scan mode.However, the probability of an OPL standard deviation less than 2 nm was 0 in the BM scan mode.As shown in Fig. 6(b), the OPL standard deviation of the BM scan had an offset ~ 3 nm due to the transverse position error of the scanner.The probability of √CRLB < 2 nm in the BM scan mode was ~0.8, which is 0.15 less than the MB scan.To quantify the OPL (phase) correction using the ML and mean estimators, we calculated the magnitude OPL differences of the ML estimated images and one of the BPAE images in each scan mode.Figure 8(a)-8(b) represents the magnitudes of the OPL correction using ML estimators in the BM and MB scan modes.Using ML estimation over 100 images, the mean OPL corrections were 7 nm and 10 nm in the MB and BM-scan modes over 256×256 pixels, respectively.The mean corrections increased by a factor of 3 over one of the BPAE cells.We believe this method is applicable for quantitative phase imaging methods where high phase sensitivity is required under low SNR conditions.Subsurface phase imaging can be one example where the structure of interest is imbedded in a turbid or transparent media.This method is not appropriate for dynamic samples with small time constants (rapid processes).Different time separations between two phase measures over a given pixel using BM and MBscan modes may enable us to image dynamic samples with different large time constants (slow processes).

Fig. 1 .
Fig. 1.The theoretical PDF of the corrupted phase (φ) for different true phase values at (a) SNR= 6 dB and (b) SNR=20 dB.

Fig. 4 .
Fig.4.The simulated precision of the ML and mean estimators as a function of data points at SNR=10 dB and true OPL=25 nm.

Fig. 5 .
Fig. 5. ML estimated intensities of the BPAE cell in the (a) BM and (b) MB scan modes.ML estimated SNRs of the BPAE cell in the (c) BM and (d) MB scan modes.

Figures 5 (
Figures 5(a)-5(d) represents the ML estimated intensities and SNRs of the BPAE cell by taking 100 images using the BM and MB scan modes.The degradation in the measured intensity and SNR in the BM scan mode could be attributed to the transverse position error of the scanner.Using the SNR map, we are able to identify those cell areas where the phase measurement accuracy was degraded by poor SNR.Figure 6(a) shows the quantitative ML estimated OPL image of the BPAE cell by taking 100 images in the MB scan mode.To determine the accuracy of our OPL images, we calculated the cumulative distribution functions (CDFs) of the OPL sensitivity based on the square root of the CRLB and the standard deviation over all pixels of 100 images for the BM and MB scan modes.The square root of the CRLB ( SNR 2 / 1) and the OPL standard deviation at each pixel were calculated using the ML estimated SNR and OPL in Eq. (21) and Eqs.(19, 22), respectively.Figures6(b)-6(c) shows CDFs for two different modes.The probability of an OPL standard deviation less than 2 nm was ~0.9, which was 0.05 less than the probability of √CRLB in the MB scan mode.However, the probability of an OPL standard deviation less than 2 nm was 0 in the BM scan mode.As shown in Fig.6(b), the OPL standard deviation of the BM scan had an offset ~ 3 nm due to the transverse position error of the scanner.The probability of √CRLB < 2 nm in the BM scan mode was ~0.8, which is 0.15 less than the MB scan.

Figure 6 (
Figures 5(a)-5(d) represents the ML estimated intensities and SNRs of the BPAE cell by taking 100 images using the BM and MB scan modes.The degradation in the measured intensity and SNR in the BM scan mode could be attributed to the transverse position error of the scanner.Using the SNR map, we are able to identify those cell areas where the phase measurement accuracy was degraded by poor SNR.Figure 6(a) shows the quantitative ML estimated OPL image of the BPAE cell by taking 100 images in the MB scan mode.To determine the accuracy of our OPL images, we calculated the cumulative distribution functions (CDFs) of the OPL sensitivity based on the square root of the CRLB and the standard deviation over all pixels of 100 images for the BM and MB scan modes.The square root of the CRLB ( SNR 2 / 1) and the OPL standard deviation at each pixel were calculated using the ML estimated SNR and OPL in Eq. (21) and Eqs.(19, 22), respectively.Figures6(b)-6(c) shows CDFs for two different modes.The probability of an OPL standard deviation less than 2 nm was ~0.9, which was 0.05 less than the probability of √CRLB in the MB scan mode.However, the probability of an OPL standard deviation less than 2 nm was 0 in the BM scan mode.As shown in Fig.6(b), the OPL standard deviation of the BM scan had an offset ~ 3 nm due to the transverse position error of the scanner.The probability of √CRLB < 2 nm in the BM scan mode was ~0.8, which is 0.15 less than the MB scan.

Fig. 6 .)
Fig. 6. (a).The quantitative ML estimated OPL image of the BPAE cell using the BM scan mode.Cumulative distribution functions of OPL sensitivity based on the square root of the CRLB ( SNR 2 / 1 ) in blue and standard deviation over 100 images in red for (b) BM and (c) MB scans.

Figures 7 (
Figures 7(a)-7(b) shows the OPL standard deviation image over of a BPAE cell using the BM and MB scans, respectively.The high standard deviation area matches the low SNR value pattern in Figs.5(c)-5(d).

Fig. 7 .
Fig. 7.The quantitative OPL standard deviation image of the BPAE cell using the ML estimator in (a) BM and (b) MB scan modes.