Open Access
Add to BibSonomyAbstract
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 Js(z) is obtained as
where Jr(z) and Jin 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 where z and z 0 are the depth positions of the point of interest and the surface, respectively. In practice, Jin is a Jones matrix of the input fiber of the interferometer; hence, the effect of fiber birefringence is canceled by Eq. (1).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.
where tr, det, and † denote the trace, determinant, and complex conjugate transpose, respectively. In this paper, the raw phase retardation δ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 Jr(z) and Jin(z), are modeled as the summation of the true Jones matrix and additive complex noise as
where Ss denote the true values of the OCT signals, and Ns 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
SNRs,i and SNRr,i are defined as the ratio of signal energy (signal intensity) of the 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.
Fig. 1 Comparison of the measured phase retardation distribution with different ESNR in the simulation and experiment. The sample is an eighth-waveplate in (a)–(d) and a quarter-waveplate in (e)–(h). (a),(c),(e), and (g) show the results of numerical simulation and (b),(d),(f), and (h) are the experimental results. The marked ESNR values denote the mean ESNR in each experiment or simulation.
Fig. 2 Skewness of the distribution of phase retardation (δM) obtained by Monte-Carlo simulation for several true phase retardations.
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. 213 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, 213 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 221 trials. For δT = 0 and π, 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 δT = 0 and π are not symmetric for any ESNR. For other δT, 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.
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 δM corresponding to several δT values are plotted against ESNR, as shown by the red curves in Figure 3. The curves are obtained by averaging 216 trials of δM 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.
Fig. 3 Estimations of several set true phase retardations using 216 trials of the measured phase retardation. The red and green solidlines denote the estimation results using mean and 4th order MCB estimators, respectively. The black dashed lines denote the true values.
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 δM into transformed phase retardation, δE, for further estimation. The ensemble average of δE approaches the true phase retardation δT as the number of samples increases;
where 〈〉 denotes ensemble average. Note that f(δM) is not only a function of δM, but also a function of ESNR γ and δT because the distribution of δM varies with them.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′(δM), which is a polynomial function of δM and ESNR γ, but not a function of δT, as
where bi(γ) is i-th order polynomial coefficient of the transform function at an ESNR of γ.In the design process of the transform function, bi(γ) is defined to be
To determine bi(γ) for a particular γ, δMs were obtained by Monte-Carlo simulations for several δT s from δT = 0 to π in steps of π/60. If bi(γ) is properly defined, the simulation results would follow the following 61 equations.
where the left-hand-sides are the set true phase retardations δT, δ M,δT is a raw phase retardation obtained by Monte-Carlo simulation with a set true phase retardation of δT, and n is the maximum polynomial order of f′(δM;γ). Here, each ensemble averaged value is obtained via Monte-Carlo simulation with 221 trials; hence, they are known values. Equation (10) can also be written in a vector form as where DT = [0,π/60,···,δT, ···,π]T, B = [b 0(γ),···,bn(γ)]T and And then, the polynomial coefficients bi(γ) are defined as where DM + is the pseudo-inverse of DM, obtained by singular value decomposition.Equation (13) is equivalent to obtaining optimal bi(γ) by a least squares method. Thus, bi(γ) is defined to minimize a squared-sum error defined as
where δT,k = kπ/60 is the true phase retardation of the k-th configuration of the Monte-Carlo simulation and δM,k is the corresponding simulated value of the measured phase retardation.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 216 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, 216 trials of δM or δE 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, π]. 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 δT = 0 and π.
In order to evaluate the performance of the MCB estimator in a more qualitative manner, we defined a mean error energy, ɛ 2, as
where δT,i = iπ/6 is true phase retardation and 〈δM/E,i〉 is the ensemble mean of δM or δE when the true value is δT,i. Figure 5 shows the square-root of ɛ 2 as a function of ESNR, where ɛ was obtained via Monte-Carlo simulation with 216 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.Fig. 5 Simulated precision of mean and MCB estimator for 216 trials. The red, green, and blue curves represent the mean error energy corresponding to a mean estimator, 4th order MCB estimator, and 6th order MCB estimator, respectively.
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 δT and ESNR, where σ 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.
Fig. 7 Estimation result in simulation and experiment. The true values are denoted by black dashed lines, and the mean and MCB estimation results are denoted by red and green lines, respectively. The solid lines denote the mean of a 65536-trial estimation, and the squares, circles, and crosses denote 64-measurement estimations of a quarter wave-plate, one eighth wave-plate, and glass in the experiment.
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.
Fig. 8 B-scan images of in vitro chicken breast muscle. (a) is an intensity image, and (b) is a log-scaled ESNR image. (c) is a single raw phase retardation image. (d) and (e) are phase retardation images obtained from mean and MCB estimators based on 128 B-scans. The white dashed lines denote the positions of the depth signal shown in Fig. 9.
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.
Fig. 9 Plots of A-line signals versus penetration depth. The red dots represent phase retardation values of raw phase retardation (a), obtained by mean estimator (b), and obtained by 4th order MCB estimator (c). Dashed curves represents corresponding ESNR values. The ESNR of (a) is a raw and non-averaged ESNR signal, while the ESNR of (b) and (c) are averaged ESNR of A-lines which have been taken at the same location on the sample and been utilized for the estimation.
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.
Fig. 10 Intensity (a) and ESNR (b) images of an in vivo posterior eye. The second row shows corresponding raw (c) and distribution-transformed (d) phase retardation images, respectively. (e) and (f) are magnified images of the area indicated by the rectangles in (c) and (d).
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 δE that is also associated with the single pixel. On the other hand, the mean operation requires multiple sampled values.
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, DM of Eq. (13) is determined by a Monte-Carlo simulation, and DT is defined by the simulation parameters.
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]
