## Abstract

Spectral-domain optical coherence phase microscopy (SD-OCPM) 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ér- Rao lower bound $\frac{1}{\sqrt{2\mathit{SNR}}}$ 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

## 1. Introduction

The development of phase imaging modalities may permit quantitative measurements on the structure and dynamics of cellular specimens [1-3]. Several phase imaging methods have been investigated including: 1) noninterferometric methods [4], 2) digital holographic microscopy [5], 3) full-field phase microscopy based on a programmable spatial light modulator [6], 4) Fourier fringe analysis [7], 5) and Hilbert transform [8]. Several quantitative phase imaging schemes in reflection by utilizing either time-domain [3] or Fourier-domain (swept source/spectral-domain) optical coherence tomography (SS/SD-OCT) have also been proposed [9-11].

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 nanometer-scale 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 SDOCPM for two scan modes and show the improved measured phase by the ML estimator.

## 2. 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*), and the imaginary part,

_{R},σ*X*(

_{I}*S*), of the corrupted data is given by,

_{I},σ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

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,

By defining the random variable *u=(R-ρ cos(ϕ-θ))/σ*, Eq. (9) can be written as

The marginal PDF of the phase is given by,

where the signal to noise ratio is defined as *SNR*=*ρ ^{2}*/(

*2σ*), and

^{2}*I*(.), and

_{0}*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 $\frac{1}{\sqrt{2\mathit{SNR}}}$ using the CRLB [13].

## 3. 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-21) are given by the vector summation of the N measurements in the complex plane.

## 4. 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:

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.

## 5. 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.

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 $\frac{1}{\sqrt{2\mathit{SNR}}}$ 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.

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).

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. Figure 8 and Figs. 5(c)-5(d) show that the corrections were mostly applied to the pixels located at OPLs equal to an integer number of half the center wavelength (ring pattern in the cell image) with low SNR where the noise dominated the signal strength. The resulting pattern in Fig. 8 also shows that the ML estimator improved the transverse position error which was due to repeatability of the B-scan pattern in the BM scan mode.

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 MB-scan modes may enable us to image dynamic samples with different large time constants (slow processes).

## Acknowledgment

This research was supported in part by the National Institute of Health under grant number NCRR 19768, and CIMIT, the Center for the Integration of Medicine and Innovative Technology with funding from the Department of the Army under its Cooperative Agreement DAMD17-02-2-0006. S.M.R. Motaghiannezam is currently with the California Institute of Technology, Pasadena, CA. C. Joo is currently with General Electric Global Research Lab, Niskayuna, NY. J.F. de Boer is currently with VU University, De Boelelaan 1081, T-0.67, Amsterdam, The Netherlands, jfdeboer@few.vu.nl.

## 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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef]

**4. **A. Barty, K. A. Nugent, D. Paganin, and A. Roberts, “Quantitative optical phase microscopy,” Opt. Lett. **23**, 817–819 (1998). [CrossRef]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**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). [CrossRef] [PubMed]

**17. **G. McGibney and M. R. Smith, “An unbiased signal-to-noise ratio measure for magnetic resonance images,” Med Phys. **20**, 1077–1078 (1993). [CrossRef] [PubMed]

**18. **B. A. van den, *Handbook of Measurement Science*, (Wiley, Chichester, England, 1982) Vol. 1, pp. 331–377.