## Abstract

A Monte-Carlo-based phase retardation estimator is developed to correct the systematic error in phase retardation measurement by polarization sensitive optical coherence tomography (PS-OCT). Recent research has revealed that the phase retardation measured by PS-OCT has a distribution that is neither symmetric nor centered at the true value. Hence, a standard mean estimator gives us erroneous estimations of phase retardation, and it degrades the performance of PS-OCT for quantitative assessment. In this paper, the noise property in phase retardation is investigated in detail by Monte-Carlo simulation and experiments. A distribution transform function is designed to eliminate the systematic error by using the result of the Monte-Carlo simulation. This distribution transformation is followed by a mean estimator. This process provides a significantly better estimation of phase retardation than a standard mean estimator. This method is validated both by numerical simulations and experiments. The application of this method to *in vitro* and *in vivo* biological samples is also demonstrated.

©2011 Optical Society of America

## 1. Introduction

Optical coherence tomography (OCT) provides high-resolution depth-resolved images of biological tissues noninvasively [1, 2], hence, OCT is suitable for applications in ophthalmology [3, 4], dermatology [5], dentistry [6], and cardiology [7]. Polarization sensitive OCT (PS-OCT), which possesses all the advantages stated above, is a functional extension of OCT. It enables both conventional backscattering tomography and birefringence tomography [8–11]. Tissue birefringence is strongly associated with the structural properties of biological tissues; hence, PS-OCT has been adopted for imaging skin [12–15], cartilage [16], teeth [17], and the anterior and posterior segments of the eye [18–25].

Phase retardation is an important birefringent property of tissue, and it is widely employed to visualize PS-OCT images. It is a cumulative quantity that rotates in phase along the depth if the sample has birefringence while it remains constant if the sample has no birefringence. In biological tissues, micro-structural changes such as fibrosis, inflammation, and canceration can result in the alteration of phase retardation [21, 23, 25]; hence, the phase retardation image is of significant diagnostic importance. However, as reported in several studies, phase retardation measurements can be drastically affected by detection noise in PS-OCT [26, 27]. A standard mean estimator, i.e., average, cannot provide an appropriate estimation of phase retardation, and it significantly reduces the utility of PS-OCT for quantitative measurement [27].

Almost all PS-OCT algorithms derive phase retardation from multiple complex OCT signals. In the real and imaginary parts of the raw OCT signal, the noise distribution is reasonably assumed to be Gaussian centered at a true value. However, in general, the derivation of phase retardation from the raw OCT signals is complicated and, sometimes, nonlinear. This elaborates derivation process drastically complicates the analytic investigation of noise distribution in phase retardation.

As is widely known, Monte-Carlo simulation is a powerful tool for investigating stochastic processes of elaborate systems, and it has been widely employed for the analysis of the optical scattering property of biological tissues [28], theirs polarization dependency [29], and imaging modalities based on optical scattering, including diffusion tomography [30], photoacoustic tomography [31], and OCT [32]. The Monte-Carlo method has also been employed to investigate the noise property of phase retardation measurement in PS-OCT [27]. In this study, it was shown that the noise distribution in phase retardation is neither Gaussian nor symmetric. This is partially because the measurement range of phase retardation of PS-OCT is typically limited from 0 to *π*; hence, phase retardation signals less than zero or greater than *π* will be aliased in the 0 to *π* range. Because of this asymmetric distribution, both the mean and the mode include a systematic error, and they cannot provide a correct estimation of true phase retardation. Currently, this systematic error is minimized by enhancing the signal-to-noise ratio (SNR) of PS-OCT. However, to have a reasonably small error, a very high SNR, typically more than 20 dB, is required.

The objective of this paper is to propose a method that reduces the systematic error in measurement and to relax the harsh SNR requirement for accurate phase retardation detection. This method involves a Monte-Carlo-based (MCB) phase retardation estimator, and it consists of the following 2 steps. In step 1, the measured raw phase retardation values are transformed by a custom distribution transform function, which is designed on the basis of the noise property of PS-OCT, determined by a Monte-Carlo analysis. In step 2, a mean operation is applied to a set of the transformed phase retardation values measured at the same location in the sample or located within a predefined averaging kernel. By simulation and wave plate experiments, we show that this MCB estimator outperforms the conventional mean estimator in terms of quantitative estimation of phase retardation. Finally, we apply the MCB estimator to phase retardation images of an *in vitro* chicken breast muscle and *in vivo* human retina. The results show improved image contrast and reasonable phase retardation, which are in good agreement with prior knowledge of retinal birefringence.

## 2. Jones Matrix Polarization Sensitive Optical Coherence Tomography

The MCB phase retardation estimator proposed in this paper is specifically designed for Jones matrix PS-OCT. The detailed principle of Jones matrix PS-OCT is described elsewhere [33–35]. Here, we briefly summarize the principle of swept-source Jones matrix PS-OCT [36].

In this system, the incident polarization is continuously modulated by an electro-optic modulator (EOM). This modulation multiplexes two OCT signals that correspond to two orthogonal polarization axes of the EOM in its carrier frequency. Thus, the OCT signal corresponding to one of the axes of the EOM has a high-modulation frequency along its spectral scan, while the spectral interference signal of the other axis is not modulated. These modulated and non-modulated spectra were detected by a polarization diversity (PD) detection unit, in which horizontal and vertical polarization components of the spectral interference signals are independently detected by two detectors. The polarization axes of the EOM and those of PD detection units are configured to be inclined at 45 degrees; hence, the modulated and non-modulated polarization components are detected by both the horizontal and the vertical detectors. The modulated and non-modulated signals were numerically demultiplexed after detection. The demultiplexed spectra were processed using a standard SS-OCT algorithm, and four OCT signals were obtained. Here, the OCT signals are denoted by *I*
_{0,}
* _{H}* (

*z*),

*I*

_{1,}

*(*

_{H}*z*),

*I*

_{0,}

*(*

_{V}*z*), and

*I*

_{1,}

*(*

_{V}*z*), where the subscripts 0 and 1 denote non-modulated and modulated signals and the subscripts

*H*and

*V*denote horizontal and vertical polarization states of the PD detector, respectively. Then, the cumulative Jones matrix at a particular point in the sample

*J*(

_{s}*z*) is obtained as

*J*(

_{r}*z*) and

*J*are the raw Jones matrix measured at the point of interest and the Jones matrix measured at the surface of the sample, respectively [36]. By using the OCT signals, these matrixes are defines as

_{in}*z*and

*z*

_{0}are the depth positions of the point of interest and the surface, respectively. In practice,

*J*is a Jones matrix of the input fiber of the interferometer; hence, the effect of fiber birefringence is canceled by Eq. (1).

_{in}The round-trip phase retardation, *δ _{M}*(

*z*), is obtained by matrix diagonalization [27] or by a trace method [37]. In this paper, we use the trace method represented by Eq. (4) rather than matrix diagonalization because it enables faster calculation.

*δ*(

_{M}*z*) is expressed as the measured phase retardation.

## 3. Noise in Phase Retardation Measurement

#### 3.1. Noise Model and Monte-Carlo Simulation

Numerical simulations are performed to investigate the error and noise properties of phase retardation measurement. In the simulations, the noise is described using the same model as that used in a previous study [27]. Thus, the raw OCT signals, including *J _{r}*(

*z*) and

*J*(

_{in}*z*), are modeled as the summation of the true Jones matrix and additive complex noise as

*S*s denote the true values of the OCT signals, and

*N*s denote complex noise in each channel. The real and imaginary parts of the noise follow zero-mean Gaussian distributions. The noises are totally decorrelated with each other, even though all of them are assumed to have identical standard deviations.

As shown in Eqs. (34)–(39) and Fig. 2 of Ref. 27, the noise property is described as a function of the effective SNR (ESNR; *γ*) rather than the SNR of each detection channel. The ESNR is defined as

*SNR*and

_{s,i}*SNR*are defined as the ratio of signal energy (signal intensity) of the

_{r,i}*i*-th incident light and noise energy of the

*i*-th detection channel at the point of interest (denoted by the subscript

*s*) and surface of the sample (denoted by the subscript

*r*), respectively. Because of this property, we utilize ESNR rather than the SNR of each detection channel in the following analysis.

It should be noted that here we follow the definition of SNR of Jones matrix OCT described in Section 3.1.1 of Ref. 27. In this definition, the signal is the summation of the signal energies of all detections for a single incident polarization state, while the noise is defined as the noise level in each single detection. In our Jones matrix PS-OCT, two detectors, i.e. for horizontal and vertical polarization states, are utilized for polarization diversity detection. When we measure a void region, the signal energy detected by a single detector is equivalent to the noise energy. According to the above mentioned definition, signal is the summation of the signal energies of the two channels, each of which is equivalent to the noise energy in this case. Hence the ESNR measured at a void region (noise region) becomes around 3 dB. This pseudo increasing of SNR is appeared only at a very low signal region.

#### 3.2. Noise Property

The phase retardation measurements and numerical simulation were performed for a one-eighth wave-plate (EWP) and a quarter wave-plate (QWP), whose round-trip phase retardations are *π*/2 and *π*, respectively. For the wave-plates measurement, a tunable neutral density filter was located in a sample arm to control the incident power and hence SNR. 2^{13} A-lines were acquired for each SNR configuration. The intensity of OCT signal in air was averaged to obtain the noise floor of each detection channels, and the SNR and ESNR were calculated for the successive estimation of phase retardation. For each simulation, 2^{13} Monte-Carlo trials were performed. Figure 1 shows the distribution of the raw measured/simulated phase retardation, *δ _{M}*, obtained from numerical simulations (left column) and experiments (right column) using an EWP (first and second rows) and a QWP (third and fourth rows). The histograms show good agreement between the simulations and the experiments. These results validate the proposed Monte-Carlo model.

Figure 1 also shows that the distributions are asymmetric. This asymmetry is shown more quantitatively in Fig. 2, where the skewness of the distributions of the true phase retardations (*δ _{T}* = 0,

*π*/6,

*π*/3,

*π*/2, 2

*π*/3, 5

*π*/6,

*π*radians) is plotted against ESNR. The skewness is defined as the third moment about the mean divided by the third power of standard deviation, and is one of the statistical measures of asymmetric level of prbability density function. Each distribution was obtained using a Monte-Carlo simulation with 2

^{21}trials. For

*δ*= 0 and

_{T}*π*, the skewness is not zero in nearly the entire ESNR range. Since skewness becomes zero for a symmetric distribution, these non-zero values of the skewness indicate that the distributions of

*δ*= 0 and

_{T}*π*are not symmetric for any ESNR. For other

*δ*, the skewness is nearly zero, i.e., the distribution is symmetric for high ESNR. However, the distribution becomes asymmetric for the ESNR less than around 20 to 25 dB. These asymmetries may be attributed to the aliasing at the perimeters of the measurement range and to the nonlinear relationship between the Jones matrix and the phase retardation.

_{T}Because of these asymmetries, the mean of the distribution cannot give a reasonable estimation of the true phase retardation *δ _{T}*. The mean of the simulated

*δ*corresponding to several

_{M}*δ*values are plotted against ESNR, as shown by the red curves in Figure 3. The curves are obtained by averaging 2

_{T}^{16}trials of

*δ*obtained by the Monte-Carlo simulation. The mean of

_{M}*δ*deviates from the true value (dashed lines) and approaches around 2.15 rad as ESNR decreases.

_{M}## 4. Monte-Carlo-Based Phase Retardation Estimator

Conventionally, mean and maximum likelihood estimators have been utilized for quantitative estimation of polarization [24]. However, the asymmetric distribution presented in Section 3.2 implies that these standard estimators do not provide appropriate estimations. This is because these estimators assume symmetric distribution of measured values. To overcome this problem, we propose an MCB estimator that involves a two-step estimation algorithm consisting of a non-linear distribution transform and a conventional mean estimation.

#### 4.1. Non-Linear Distribution Transform and MCB Estimator

The first step of our MCB estimator is the distribution transform. We assume a transform function *f*(*δ _{M}*) that transforms measured raw phase retardation

*δ*into transformed phase retardation,

_{M}*δ*, for further estimation. The ensemble average of

_{E}*δ*approaches the true phase retardation

_{E}*δ*as the number of samples increases;

_{T}*f*(

*δ*) is not only a function of

_{M}*δ*, but also a function of ESNR

_{M}*γ*and

*δ*because the distribution of

_{T}*δ*varies with them.

_{M}In practical phase retardation estimation, *δ _{T}* cannot be utilized as prior information, while ESNR can be obtained prior to the estimation directly from the measured OCT signals. Hence, we define a set of the suboptimal but practical transform function

*f*′(

*δ*), which is a polynomial function of

_{M}*δ*and ESNR

_{M}*γ*, but not a function of

*δ*, as

_{T}*b*(

_{i}*γ*) is

*i*-th order polynomial coefficient of the transform function at an ESNR of

*γ*.

In the design process of the transform function, *b _{i}*(

*γ*) is defined to be

To determine *b _{i}*(

*γ*) for a particular

*γ*,

*δ*s were obtained by Monte-Carlo simulations for several

_{M}*δ*s from

_{T}*δ*= 0 to

_{T}*π*in steps of

*π*/60. If

*b*(

_{i}*γ*) is properly defined, the simulation results would follow the following 61 equations.

*δ*,

_{T}*δ*

_{M,δT}is a raw phase retardation obtained by Monte-Carlo simulation with a set true phase retardation of

*δ*, and

_{T}*n*is the maximum polynomial order of

*f*′(

*δ*;

_{M}*γ*). Here, each ensemble averaged value is obtained via Monte-Carlo simulation with 2

^{21}trials; hence, they are known values. Equation (10) can also be written in a vector form as

**D**= [0,π/60,···,

_{T}*δ*, ···,

_{T}*π*]

*,*

^{T}**B**= [

*b*

_{0}(γ),···,

*b*(

_{n}*γ*)]

*and*

^{T}*b*(

_{i}*γ*) are defined as

**D**

_{M}^{+}is the pseudo-inverse of

**D**, obtained by singular value decomposition.

_{M}Equation (13) is equivalent to obtaining optimal *b _{i}*(

*γ*) by a least squares method. Thus,

*b*(

_{i}*γ*) is defined to minimize a squared-sum error defined as

*δ*=

_{T,k}*kπ*/60 is the true phase retardation of the

*k*-th configuration of the Monte-Carlo simulation and

*δ*is the corresponding simulated value of the measured phase retardation.

_{M,k}A transform function *f*′(*δ _{M}*;

*γ*) is now defined for particular

*γ*. For practical applications of phase retardation measurement,

*f*′(

*δ*;

_{M}*γ*) is defined for

*γ*of 0 dB to 50 dB in steps of 1 dB.

This transform function can be utilized to correct the contrast of phase retardation tomography images. Further, to obtain a qualitative phase retardation value, one of standard estimators can be applied after this distribution transform. According to the assumption made in the design of the transform function, a mean estimator would be reasonable. In this paper, we denote this combination of the distribution transform and mean estimator as an MCB estimator.

#### 4.2. Performance of Monte-Carlo-Based Estimator

### 4.2.1. Estimation Error

The detailed performance of the MCB estimator was evaluated via Monte-Carlo simulations. Figure 3 shows the estimated phase retardation at several ESNR and true phase retardations, where 2^{16} trials of Monte-Carlo simulation were performed for each configuration. The green curves were obtained using the MCB estimator with the 4th order transform function, and the red curves were obtained using a mean estimator, as discussed in Section 3.2. As shown in this figure, the MCB estimator gives reasonable estimation even with a low ESNR of around 5 dB, while a standard mean estimator suffers from a significant error, even at an ESNR of 20 dB.

For further understanding of this error property, the estimation error is plotted as a function of the true phase retardation and ESNR, as shown in Fig. 4. Figures 4(a) and 4(b) show the estimation error of mean estimator and 4th order MCB estimator, respectively. In these estimations, 2^{16} trials of *δ _{M}* or

*δ*were averaged. Fig. 4(a) clearly shows that the mean estimator includes a significant error if the true value is not close to 2.15 rad, and the error becomes larger as the true value approaches the perimeter of the measurement range of phase retardation, i.e. [0,

_{E}*π*]. In contrast, the error is well controlled by the MCB estimator, as shown in Fig. 4(b). When the ESNR is higher than 5 dB, the error is less than

*π*/50 for most of the part of the plot. The maximum error is less than

*π*/20 rad, which appears at

*δ*= 0 and

_{T}*π*.

In order to evaluate the performance of the MCB estimator in a more qualitative manner, we defined a mean error energy, *ɛ*
^{2}, as

*δ*=

_{T,i}*iπ*/6 is true phase retardation and 〈

*δ*〉 is the ensemble mean of

_{M/E,i}*δ*or

_{M}*δ*when the true value is

_{E}*δ*. Figure 5 shows the square-root of

_{T,i}*ɛ*

^{2}as a function of

*ESNR*, where

*ɛ*was obtained via Monte-Carlo simulation with 2

^{16}trials, followed by a mean estimator (red curve), 4th order MCB estimator (green curve), and 6th order MCB estimator (blue curve). For very high ESNR, e.g., higher than 25 dB, all the estimators provide a reasonably small amount of error. However, if the ESNR becomes lower than 25 dB, the error of the mean estimator increases rapidly, while the errors of the MCB estimators remain reasonably small until around 8 dB.

### 4.2.2. Randomness of Estimation

Thus far, we discussed the estimation error obtained with a sufficiently large number of trials. However, in practice, the number of trials/measurements we can perform at a single location is limited. Similarly, the kernel size of a local averaging filter that may be applied to a phase retardation image should be small. If the distribution of the transformed phase retardation is broad, i.e., randomeness is high, a large number of measurements are required to restrain the randomness of the result. Hence, the randomness of estimation is also of particular interest.

For qualitative evaluation, we define randomness as the root mean square error (RMSE) of the estimation as

In Fig. 6, *σ* is plotted as a function of ESNR and *δ _{T}*. We consider

*σ*as a performance criterion of precision in estimation. Figures 6(a)–6(c) show

*σ*against

*δ*and ESNR, where

_{T}*σ*was obtained with mean, 4th order MCB, and 6th order MCB estimation, respectively.

*σ*of all of these three methods increase as ESNR decreases. The 4th and 6th order MCB estimators possess a similar level of randomness if the ESNR is higher than 15 dB. However, the randomness in the 6th order MCB becomes significantly higher than that in the 4th order MCB if the ESNR is less than 10 dB. Therefore, we have used the 4th order transform function rather than a higher order transform function. In summary, the randomness in estimation increases as ESNR decreases, and the mean and higher order MCB generate the lowest and highest random error.

## 5. Experimental Validation

#### 5.1. Standard Samples

The experimental performance of the MCB estimator is quantitatively evaluated by measuring the phase retardations of a glass plate, a one-eighth wave-plate, and a quarter wave-plate, whose round trip phase retardations are 0, *π*/2, and *π*, respectively.

For this measurement, PS-SS-OCT with a 1.3 *μ*m probe was employed. The details of the principle of this system are described in Section 2 and Ref. 36, and the details of the hardware configuration are described in Ref. 38. In short, this system has a depth resolution of 8.3 *μ*m in air and a measurement speed of 30,000 A-lines/s. The measurements are performed with several ESNR configurations, which are controlled by a neutral density filter placed in front of the sample. For each ESNR configuration, 64 A-lines were successively obtained and utilized for the estimation.

The estimations are shown in Fig. 7, where the red dots and green dots represent the estimations obtained by a mean estimator and MCB estimator, respectively. The red and green curves denote the corresponding estimations obtained by the Monte-Carlo simulations, and they are identical to those in Fig. 3. The experiment is in good agreement with the Monte-Carlo simulation, and the MCB estimator enables better estimation than the mean estimator.

The estimation was not performed for ESNR less than 8 dB because of the difficulty in the experimental determination of ESNR, which arises from the non-uniform SNR of each detection channel of the PS-SS-OCT. Because the modulation efficiency of the source polarization of the system is not perfect, the SNR of the modulated channels are 8 to 10 dB lower than those of the non-modulated channels. The ESNR (*γ*), which is a significant parameter for the estimation, becomes unreliable in this low ESNR range; hence, *δ _{E}* also becomes unreliable.

#### 5.2. In Vitro Chicken Breast Muscle

The MCB estimator was applied to an *in vitro* chicken breast muscle, which is commonly utilized as a standard sample of PS-OCT evaluation. A 10-mm region on the sample was scanned to obtain a B-scan that consists of 512 A-lines. The B-scan measurements were successively performed 128 times at the same location on the sample; hence, 128 measurements were performed at a single point of an OCT image. Before applying the estimator, data points whose ESNR is lower than 8 dB were discarded because of their unreliability (see the last paragraph of Section 5.1). The remaining data points were substituted by the mean and 4th order MCB estimators to form estimated-phase retardation images.

Figure 8 shows the intensity image (a), ESNR image (b), a raw (non-avaraged) phase retardation image (c) and estimated phase retardation images formed by the mean estimator (d) and MCB estimator (e). Although the mean estimator significantly reduces the noise as shown in Figs. 8(c) and 8(d), the contrast of the fringe shown in Fig. 8(d) is still erroneously decreases in a deeper region. This is because, with low ESNR, the mean estimator underestimates the phase retardation when the true phase retardation is close to *π*, and overestimates the phase retardation when the true phase retardation is close to zero, as discussed in Section 4.2.1. On the other hand, the MCB estimator provides higher contrast of the fringe, as shown in Fig. 8(e). It should be noted that with the MCB estimator, the fringe contrast decreases in a very deep region. This may be attributed to OCT signal crosstalk due to multiple scattering and the limitations of MCB discussed in Section 4.2.1.

The red dots in Figs. 9(a)–9(c) respectively show examples of depth resolved phase retardation signals extracted from a raw phase retardation image (a) and from phase retardation images obtained by the mean estimator (b) and the MCB estimator (c). The transversal locations of these signals are indicated by dashed lines in Figs. 8(c)–8(e). The dashed curves present corresponding depth-resolved ESNR signals in a logarithmic scale. The ESNR of Fig. 9(a) is a non-averaged ESNR signal, while the ESNR of Figs. 9(b)–9(c) are averaged ESNR of A-lines which have been taken at the same location on the sample and been utilized for the estimation. It should be noted that the ESNR in the air is not 0 dB but around 3 dB because of the reason discussed in Section 3.1. The phase retardation is cumulative along depth; hence, a clear fringe pattern is observed in the estimated phase retardation by both methods. A clear decay was observed in the contrast of phase retardation fringe in the mean estimator image, caused by ESNR decreasing along the penetration depth; however, this did not happen in estimation using the MCB method.

#### 5.3. In Vivo Posterior Eye Measurement

To show the utility of the distribution transform for correcting the contrast of phase retardation images without successive mean averaging, we applied this transform to a phase retardation image of an *in vivo* human eye. Here, only a single phase retardation image is utilized. The measurement was performed by PS-SS-OCT with a 1.0 *μ*m probe beam with a depth resolution of 11.0 *μ*m in tissue and a scanning speed of 30,000 A-lines/sec [39]. A 6-mm region on the retina was scanned once, and a B-scan with 2048 A-lines were obtained. The distribution transform was then applied to the phase retardation image.

Figure 10 shows the intensity (a), ESNR (b), raw phase retardation (c), and distribution-transformed phase retardation tomography (d). In the phase retardation images, the pixel which has lower SNR than 5dB in the intensity image was masked and appeared as a black. It is known that the sensory retina at the macular region has a very small or no birefringence because the nerve fiber layer is very thin in the parafoveal region or absent at the fovea. However, in the raw phase retardation tomography (Fig. 10(c) and its magnified version of Fig. 10(e)), the inner retina mainly appears as blue, which indicates moderate phase retardation in the sensory retina. In contrast, the same part of the retina appears as green in the distribution-transformed phase retardation image (Fig. 10(d) and its magnified version of Fig. 10(f)), which indicates the absence of birefringence in the sensory retina. The distribution-corrected phase retardation image can describe the birefringence properties of the eye more reasonably.

## 6. Discussion

#### 6.1. Full-Step MCB Estimation of In Vivo Tomography

As discussed above, the MCB estimator involves two steps, the distribution transform and ensemble mean operation. The distribution transform is a point-to-point function, i.e., the inputs of the distribution transform function are ESNR and *δ _{M}* of a single pixel, and the output is a single value of

*δ*that is also associated with the single pixel. On the other hand, the mean operation requires multiple sampled values.

_{E}Because of this property, the mean operation requires multiple phase retardation images measured at exactly the same location on the sample, while the distribution transform function requires only one image. This is an undesirable requirement for *in vivo* measurement.

This limitation can be overcome by three strategies. The first strategy is to apply only the transform function, as explained in Section. 5.3. This strategy is useful for qualitative observation; however, is not sufficient for quantitative assessment.

The second strategy is the use of the local mean as an alternative to the ensemble mean, in which the phase retardation values within a predefined special extent, known as a kernel, are averaged.

The third strategy is the use of slope fitting. Depth-oriented slope fitting of phase retardation has been widely utilized for the quantification of birefringence [12, 26, 40]. The slope fitting, i.e., least square linear fitting, is a maximum likelihood estimation with an assumption of symmetric distribution of noise, and the mean is a special case of linear fitting in which the slope is predefined as zero. Both the slope fitting and the mean are maximum likelihood estimations with the assumption of symmetric noise distribution; hence, the slope fitting would be applicable as an alternative to the mean operation. In the proposed scheme, the slope fitting of phase retardation is applied along the depth after the distribution transform. This would provide better estimation of birefringence of a biological specimen than the standard slope fitting without a distribution transform.

#### 6.2. Applicability of MCB Estimator to Other PS-OCT Algorithms

It is known that a similar systematic error exists in other PS-OCT algorithms including the widely utilized Hee’s algorithm [8, 26]. Although the MCB estimator presented in this paper cannot be directly applied to other PS-OCT algorithms, it may be possible to design a tailored MCB estimator for the other algorithms by using Eqs. (11)–(13). In this designing process, the noise model of the PS-OCT should be properly customized for the PS-OCT algorithm of interest. Then, **D _{M}** of Eq. (13) is determined by a Monte-Carlo simulation, and

**D**is defined by the simulation parameters.

_{T}In our Jones matrix PS-OCT, the SNR of each detection channel does not independently affect the error property; however, the ESNR *γ* dominates this property. Thus, **B** was determined for each *γ* value. However, this ESNR dependency is not warranted in other PS-OCT algorithms. The independent factors of the error property should be carefully considered in the design process.

## 7. Conclusion

We proposed a nonlinear method to estimate a correct phase retardation value from raw phase retardation values measured via PS-OCT. This estimator involves of two operations, distribution transform and mean operation. The distribution transform function was designed to eliminate the asymmetric distribution of phase retardation, which causes a systematic error in the mean estimation of phase retardation. The transform function was designed on the basis of Monte-Carlo analysis of the error property of PS-OCT. The superior performance of the MCB estimator, as compared to that of a standard mean estimator, has been demonstrated by numerical simulations and experiments. The MCB estimator was also applied to an *in vitro* chicken breast muscle, and it showed higher contrast of the phase retardation fringe than a standard mean estimator. The distribution transform was applied to an *in vivo* human posterior eye, and a reasonable phase retardation image was obtained. MCB can potentially improve the phase retardation image quality of PS-OCT; moreover, it can further will enhance the ability of quantitative phase retardation measurement of PS-OCT.

## Acknowledgments

This study is partially supported by the Japan Science and Technology Agency through the contract of the development program of advanced measurement systems, and a research grant from Tomey Corporation.

## References and links

**1. **D. Huang, E. Swanson, C. Lin, J. Schuman, W. Stinson, W. Chang, M. Hee, T. Flotte, K. Gregory, C. Puliafito, and J. Fujimoto, “Optical coherence tomography,” Science **254**, 1178–1181 (1991). [CrossRef] [PubMed]

**2. **A. Fercher, W. Drexler, C. Hitzenberger, and T. Lasser, “Optical coherence tomography - principles and applications,” Rep. Prog. Phys. **66**, 239–303 (2003). [CrossRef]

**3. **M. Hee, J. Izatt, E. Swanson, D. Huang, J. Schuman, C. Lin, C. Puliafito, and J. Fujimoto, “Optical coherence tomography of the human retina,” Arch. Ophthalmol. **113**, 325–332 (1995). [CrossRef] [PubMed]

**4. **J. G. Fujimoto, W. Drexler, J. S. Schuman, and C. K. Hitzenberger, “Optical coherence tomography (OCT) in ophthalmology: Introduction,” Opt. Express17, 3978–3979 (2009), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-17-5-3980. [CrossRef]

**5. **J. Welzel, “Optical coherence tomography in dermatology: a review,” Skin Res. Tech. **7**, 1–9 (2001). [CrossRef]

**6. **L. Otis, M. Everett, U. Sathyam, and B. Colston, “Optical coherence tomography: A new imaging technology for dentistry,” J. Am. Dent. Assoc. **131**, 511+ (2000). [PubMed]

**7. **I. Jang, B. Bouma, D. Kang, S. Park, S. Park, K. Seung, K. Choi, M. Shishkov, K. Schlendorf, E. Pomerantsev, S. Houser, H. Aretz, and G. Tearney, “Visualization of coronary atherosclerotic plaques in patients using optical coherence tomography: Comparison with intravascular ultrasound,” J. Am. Coll. Cardiol. **39**, 604–609 (2002). [CrossRef] [PubMed]

**8. **M. Hee, D. Huang, E. Swanson, and J. Fujimoto, “Polarization-sensitive low-coherence reflectometer for birefringence characterization and ranging,” J. Opt. Soc. Am. B9, 903–908 (1992), http://www.opticsinfobase.org/abstract.cfm?URI=josab-9-6-903. [CrossRef]

**9. **J. de Boer, T. Milner, and J. Nelson, “Determination of the depth-resolved Stokes parameters of light backscattered from turbid media by use of polarization-sensitive optical coherence tomography,” Opt. Lett.24, 300–302 (1999), http://www.opticsinfobase.org/ol/abstract.cfm?URI=ol-24-5-300. [CrossRef]

**10. **G. Yao and L. V. Wang, “Two-dimensional depth-resolved Mueller matrix characterization of biological tissue by optical coherence tomography,” Opt. Lett.24, 537–539 (1999), http://www.opticsinfobase.org/ol/abstract.cfm?URI=ol-24-8-537. [CrossRef]

**11. **S. Jiao and L. V. Wang, “Jones-matrix imaging of biological tissues with quadruple-channel optical coherence tomography,” J. Biomed. Opt. **7**, 350–358 (2002). [CrossRef] [PubMed]

**12. **B. Park, C. Saxer, S. Srinivas, J. Nelson, and J. de Boer, “In vivo burn depth determination by high-speed fiber-based polarization sensitive optical coherence tomography,” J. Biomed. Opt. **6**, 474–479 (2001). [CrossRef] [PubMed]

**13. **Y. Hori, Y. Yasuno, S. Sakai, M. Matsumoto, T. Sugawara, V. Madjarova, M. Yamanari, S. Makita, T. Yasui, T. Araki, M. Itoh, and T. Yatagai, “Automatic characterization and segmentation of human skin using three-dimensional optical coherence tomography,” Opt. Express14, 1862–1877 (2006), http://www.opticsinfobase.org/abstract.cfm?URI=oe-14-5-1862. [CrossRef] [PubMed]

**14. **S. Sakai, M. Yamanari, A. Miyazawa, M. Matsumoto, N. Nakagawa, T. Sugawara, K. Kawabata, T. Yatagai, and Y. Yasuno, “In vivo three-dimensional birefringence analysis shows collagen differences between young and old photo-aged human skin,” J. Invest. Dermatol. **128**, 1641–1647 (2008). [CrossRef] [PubMed]

**15. **S. Sakai, N. Nakagawa, M. Yamanari, A. Miyazawa, Y. Yasuno, and M. Matsumoto, “Relationship between dermal birefringence and the skin surface roughness of photoaged human skin,” J. Biomed. Opt. **14**, 044032 (2009). [CrossRef] [PubMed]

**16. **W. Drexler, D. Stamper, C. Jesser, X. Li, C. Pitris, K. Saunders, S. Martin, M. Lodge, J. Fujimoto, and M. Brezinski, “Correlation of collagen organization with polarization sensitive imaging of in vitro cartilage: implications for osteoarthritis,” J. Rheumatol. **28**, 1311–1318 (2001). [PubMed]

**17. **A. Baumgartner, S. Dichtl, C. Hitzenberger, H. Sattmann, B. Robl, A. Moritz, Z. Fercher, and W. Sperr, “Polarization-sensitive optical coherence tomography of dental structures,” Caries Res. **34**, 59–69 (2000). [CrossRef]

**18. **Y. Yasuno, M. Yamanari, K. Kawana, T. Oshika, and M. Miura, “Investigation of post-glaucoma-surgery structures by three-dimensional and polarization sensitive anterior eye segment optical coherence tomography,” Opt. Express17, 3980–3996 (2009), http://www.opticsexpress.org/abstract.cfm?URI=oe-17-5-3980. [CrossRef] [PubMed]

**19. **F. Fanjul-Velez, M. Pircher, B. Baumann, E. Goetzinger, C. K. Hitzenberger, and J. L. Arce-Diego, “Polarimetric analysis of the human cornea measured by polarization-sensitive optical coherence tomography,” J. Biomed. Opt. **15** (2010). [CrossRef] [PubMed]

**20. **Y. Yasuno, M. Yamanari, K. Kawana, M. Miura, S. Fukuda, S. Makita, S. Sakai, and T. Oshika, “Visibility of trabecular meshwork by standard and polarization-sensitive optical coherence tomography,” J. Biomed. Opt. **15**, 061705 (2010). [CrossRef]

**21. **B. Cense, T. Chen, B. Park, M. Pierce, and J. de Boer, “Thickness and birefringence of healthy retinal nerve fiber layer tissue measured with polarization-sensitive optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. **45**, 2606–2612 (2004). [CrossRef] [PubMed]

**22. **E. Götzinger, M. Pircher, and C.K. Hitzenberger, “High speed spectral domain polarization sensitive optical coherence tomography of the human retina,” Opt. Express13, 10217–10229 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?uri=oe-13-25-10217. [CrossRef]

**23. **M. Pircher, E. Goetzinger, O. Findl, S. Michels, W. Geitzenauer, C. Leydolt, U. Schmidt-Erfurth, and C. K. Hitzenberger, “Human macula investigated in vivo with polarization-sensitive optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. **47**, 5487–5494 (2006). [CrossRef] [PubMed]

**24. **M. Yamanari, M. Miura, S. Makita, T. Yatagai, and Y. Yasuno, “Phase retardation measurement of retinal nerve fiber layer by polarization-sensitive spectral-domain optical coherence tomography and scanning laser polarimetry,” J. Biomed. Opt. **13**, 014013 (2008). [CrossRef] [PubMed]

**25. **M. Miura, M. Yamanari, T. Iwasaki, A. E. Elsner, S. Makita, T. Yatagai, and Y. Yasuno, “Imaging polarimetry in age-related macular degeneration,” Invest. Ophthalmol. Vis. Sci. **49**, 2661–2667 (2008). [CrossRef] [PubMed]

**26. **M. Everett, K. Schoenenberger, B. Colston, and L. Da Silva, “Birefringence characterization of biological tissue by use of optical coherence tomography,” Opt. Lett.23, 228–230 (1998), http://www.opticsinfobase.org/ol/abstract.cfm?URI=ol-23-3-228. [CrossRef]

**27. **S. Makita, M. Yamanari, and Y. Yasuno, “Generalized Jones matrix optical coherence tomography: performance and local birefringence imaging,” Opt. Express18, 854–876 (2010), http://www.opticsexpress.org/abstract.cfm?URI=oe-18-2-854. [CrossRef] [PubMed]

**28. **L. V. Wang, “Mechanisms of ultrasonic modulation of multiply scattered coherent light: a monte carlo model,” Opt. Lett.26, 1191–1193 (2001), http://www.opticsinfobase.org/ol/abstract.cfm?URI=ol-26-15-1191. [CrossRef]

**29. **G. Yao and L. V. Wang, “Propagation of polarized light in turbid media: simulated animation sequences,” Opt. Express7, 198–203 (2000), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-7-5-198. [CrossRef]

**30. **S.-P. Lin, L. Wang, S. L. Jacques, and F. K. Tittel, “Measurement of tissue optical properties by the use of oblique-incidence optical fiber reflectometry,” Appl. Opt.36, 136–143 (1997), http://ao.osa.org/abstract.cfm?URI=ao-36-1-136. [CrossRef] [PubMed]

**31. **Z. Xie, L. V. Wang, and H. F. Zhang, “Optical fluence distribution study in tissue in dark-field confocal photoacoustic microscopy using a modified monte carlo convolution method,” Appl. Opt.48, 3204–3211 (2009), http://ao.osa.org/abstract.cfm?URI=ao-48-17-3204. [CrossRef] [PubMed]

**32. **D. Smithies, T. Lindmo, Z. Chen, J. Nelson, and T. Milner, “Signal attenuation and localization in optical coherence tomography studied by Monte Carlo simulation,” Phys. Med. Biol. **43**, 3025–3044 (1998). [CrossRef] [PubMed]

**33. **S. Jiao, W. Yu, G. Stoica, and L. Wang, “Optical-fiber-based mueller optical coherence tomography,” Opt. Lett.28, 1206–1208 (2003), http://ol.osa.org/abstract.cfm?URI=ol-28-14-1206. [CrossRef] [PubMed]

**34. **B. H. Park, M. C. Pierce, B. Cense, and J. F. de Boer, “Jones matrix analysis for a polarization-sensitive optical coherencetomography system using fiber-optic components,” Opt. Lett.29, 2512–2514 (2004), http://ol.osa.org/abstract.cfm?URI=ol-29-21-2512. [CrossRef] [PubMed]

**35. **M. Yamanari, S. Makita, V. D. Madjarova, T. Yatagai, and Y. Yasuno, “Fiber-based polarization-sensitive fourier domain optical coherence tomography using b-scan-oriented polarization modulation method,” Opt. Express14, 6502–6515 (2006), http://www.opticsinfobase.org/abstract.cfm?URI=oe-14-14-6502. [CrossRef] [PubMed]

**36. **M. Yamanari, S. Makita, and Y. Yasuno, “Polarization-sensitive swept-source optical coherence tomography with continuous source polarization modulation,” Opt. Express16, 5892–5906 (2008), http://www.opticsexpress.org/abstract.cfm?URI=oe-16-8-5892. [CrossRef] [PubMed]

**37. **S.-Y. Lu and R. A. Chipman, “Homogeneous and inhomogeneous Jones matrices,” J. Opt. Soc. Am. A11, 766–773 (1994), http://josaa.osa.org/abstract.cfm?URI=josaa-11-2-766. [CrossRef]

**38. **Y. Lim, M. Yamanari, and Y. Yasuno, “Polarization sensitive corneal and anterior segment swept-source optical coherence tomography,” Proc. SPIE **7550**, 75500O (2010). [CrossRef]

**39. **M. Yamanari, Y. Lim, S. Makita, and Y. Yasuno, “Visualization of phase retardation of deep posterior eye by polarization-sensitive swept-sourceoptical coherence tomography with 1-μm probe,” Opt. Express17, 12385–12396 (2009), http://www.opticsexpress.org/abstract.cfm?URI=oe-17-15-12385. [CrossRef] [PubMed]

**40. **W.-C. Kuo, N.-K. Chou, C. Chou, C.-M. Lai, H.-J. Huang, S.-S. Wang, and J.-J. Shyu, “Polarization-sensitive optical coherence tomography for imaging human atherosclerosis,” Appl. Opt. **46**, 2520–2527 (2007), http: //ao.osa.org/abstract.cfm?URI=ao-46-13-2520. [CrossRef] [PubMed]