## 1.

## Introduction

Polarization-sensitive optical coherence tomography (PSOCT)^{1} enhances conventional optical coherence tomography (OCT)^{2} by revealing polarization properties, such as retardance and optical axis, in addition to structural information. Such polarization-dependent contrast has potential applications in a variety of areas, such as ophthalmology, dermatology, and dentistry. Due to the “round-trip” signal detection in OCT, conventional PSOCT measurements are integrated results from the sample surface to the measurement depth. Because such cumulative results are not intuitive for interpretation and incorrect for samples with depth-varying optical axis, various algorithms^{3}4.5.6.7.8.^{–}^{9} have been developed to extract depth-resolved local polarization properties. The majority of these reported methods are Jones matrix-based^{5}6.^{–}^{7} OCT systems that require multiple imaging acquisitions with different incident polarization. These algorithms cannot be directly applied to the widely adopted conventional PSOCT systems that use a single circularly polarized incident light. An algorithm^{8} was recently reported to extract local retardance in birefringent samples using a conventional PSOCT system. However, the local optical axis cannot be obtained by using this algorithm. Here, we extended the previous algorithm^{8} and developed a new iterative method that can derive a true depth-resolved local optical axis along with local retardance in conventional PSOCT systems.

## 2.

## Methods

In a conventional PSOCT system, a circularly polarized light is usually used as the incident light, and the horizontally (H) and vertically (V) polarized backscattered light from the sample are measured. The “single-trip” cumulative retardance $\delta $ and “fast” optical axis $\theta $ can then be obtained in addition to the intensity $I=10\text{\hspace{0.17em}}{\mathrm{log}}_{10}({H}^{2}+{V}^{2})$ as:^{8}

## (1)

$$\delta ={\mathrm{tan}}^{-1}(|\mathrm{H}|/|\mathrm{V}|)\phantom{\rule{0ex}{0ex}}\theta =0.5\text{\hspace{0.17em}}{\mathrm{tan}}^{-1}[\mathrm{Im}(\mathrm{H}\times {\mathrm{V}}^{*})/\mathrm{Re}(\mathrm{H}\times {\mathrm{V}}^{*})].$$^{10}and the cumulative retardance $\delta $ is the true integrated result over the imaging path. However, the same interpretation is not true in samples with depth-varying optical axis.

^{8}

To construct the local optical axis and retardance (denoted as $\varphi $ and $\kappa $, respectively) at a particular depth in a birefringent sample, the corresponding $i$’th layer (represented as the $i$’th image pixel in the $A$-scan) is modeled as a linear retarder whose Jones matrix ${\mathbf{J}}_{i}$ is:^{11}

## (2)

$${\mathbf{J}}_{i}={\mathbf{J}}_{\mathrm{L}}({\varphi}_{i},{\kappa}_{i})=\mathbf{R}({\varphi}_{i})\mathbf{\Lambda}({\kappa}_{i})\mathbf{R}(-{\varphi}_{i}),$$^{11}

## (3)

$$\mathbf{R}(\varphi )=\left(\begin{array}{cc}\mathrm{cos}\text{\hspace{0.17em}}\varphi & -\mathrm{sin}\text{\hspace{0.17em}}\varphi \\ \mathrm{sin}\text{\hspace{0.17em}}\varphi & \mathrm{cos}\text{\hspace{0.17em}}\varphi \end{array}\right),\phantom{\rule[-0.0ex]{2em}{0.0ex}}\mathbf{\Lambda}(\kappa )=\left(\begin{array}{cc}{e}^{i\kappa /2}& 0\\ 0& {e}^{-i\kappa /2}\end{array}\right).$$## (4)

$${\mathbf{J}}_{\mathrm{RT}}(n)={[{\mathbf{J}}_{\mathrm{ST}}(n-1)]}^{\mathrm{T}}({\mathbf{J}}_{n}^{\mathrm{T}}{\mathbf{J}}_{n})[{\mathbf{J}}_{\mathrm{ST}}(n-1)],$$## (6)

$${\mathbf{J}}_{n}^{\mathrm{T}}{\mathbf{J}}_{n}={[{\mathbf{J}}_{\mathrm{ST}}^{\mathrm{T}}(n-1)]}^{-1}{\mathbf{J}}_{\mathrm{RT}}(n){[{\mathbf{J}}_{\mathrm{ST}}(n-1)]}^{-1}.$$^{7}

From Eq. (6), both round-trip Jones matrix of the top n layers ${\mathbf{J}}_{\mathrm{RT}}(n)$ and the single-trip Jones matrix of the top ($n-1$) layers ${\mathbf{J}}_{\mathrm{ST}}(n-1)$ are needed to calculate the local Jones matrix ${\mathbf{J}}_{n}$. The ${\mathbf{J}}_{\mathrm{RT}}(n)$ in a birefringent sample is also a linear retarder:^{8}

## (7)

$${\mathbf{J}}_{\mathrm{RT}}(n)={\mathbf{J}}_{\mathrm{L}}({\theta}_{n},2{\delta}_{n})=\mathbf{R}({\theta}_{n})\mathbf{\Lambda}(2{\delta}_{n})\mathbf{R}(-{\theta}_{n}),$$To calculate the single-trip Jones matrix ${\mathbf{J}}_{\mathrm{ST}}(n-1)$, a recursive algorithm can be applied starting from the first layer where ${\mathbf{J}}_{\mathrm{ST}}(1)={\mathbf{J}}_{1}={\mathbf{J}}_{L}({\theta}_{1},{\delta}_{1})$. Once the local retardance and optical axis are determined for the $n$’th layer, the next single-trip cumulative Jones matrix for the top $n$ layers is constructed as:

## (8)

$${\mathbf{J}}_{\mathrm{ST}}(n)={\mathbf{J}}_{n}{\mathbf{J}}_{\mathrm{ST}}(n-1)={\mathbf{J}}_{L}({\varphi}_{n},{\kappa}_{n}){\mathbf{J}}_{\mathrm{ST}}(n-1),$$Instead of calculating the local retardance ${\kappa}_{n}$ from ${\mathbf{J}}_{n}^{\mathrm{T}}{\mathbf{J}}_{n}$ along with the optical axis ${\varphi}_{n}$, we obtained ${\kappa}_{n}$ from the PSOCT measurements of the two neighboring layers as reported previously.^{8} Briefly, from the PSOCT measurements obtained at depth layer $n-1$ (${\delta}_{n-1}$ and ${\theta}_{n-1}$) and depth layer $n$ (${\delta}_{n}$ and ${\theta}_{n}$), the following matrix is constructed:

## (9)

$$\mathbf{N}(n)={\mathbf{J}}_{L}^{-1}({\theta}_{n-1},{\delta}_{n-1}){\mathbf{J}}_{\mathrm{L}}({\theta}_{n},2{\delta}_{n}){\mathbf{J}}_{L}^{-1}({\theta}_{n-1},{\delta}_{n-1}).$$^{7}

## 3.

## Results

The above algorithm was tested on PSOCT images acquired using a bulk-optic frequency domain PSOCT system described elsewhere.^{12} The system has a central wavelength of 844 nm. The depth and lateral resolution of the system are 6.7 and 37 *μ*m in air, respectively. Imaging speed is 50 k $A-\text{lines}/s$. The structure image $I$, cumulative phase retardation $\delta $, relative optical axis $\theta $ were calculated by Eq. (1). A $3\times 3$ (pixel) median filter was applied to the measured H and V components to remove the speckle noise prior to applying the proposed algorithm.

An initial test was conducted to verify the algorithm in a tendon sample with uniform fiber orientation. A piece of chicken tendon was slightly stretched and mounted so that the sample surface was roughly perpendicular to the incident light. A series of PSOCT images were obtained and processed using the aforementioned algorithm while the tendon sample was rotated from around $-10$ to 80 deg with a 10 deg step. The 0 deg was defined as the position when tendon fiber orientation was parallel with the B-scan direction [Fig. 1(a)]. The calculated fast optical axis angle was shifted by 90 deg to display the fibrous orientation, which is along the “slow” optical axis.

Banded patterns appeared in the images of cumulative retardance $\delta $ and optical axis $\theta $ as expected [Fig. 1(b)]. The images of local retardance $\kappa $ were relative homogeneous, suggesting a uniform distribution of optical retardance in the sample. The calculated images of local optical axis $\varphi $ were also uniform at each rotation angle but had different values at different angles. For quantitative comparison, the mean and standard deviation values of the local retardance and optical axis were calculated for a region of interest (ROI) enclosed by the box in Fig. 1(b). As shown in Fig. 1(c), the mean local retardance remained stable at different rotation angles. However, the mean local axis correctly retrieved the sample orientation. The axis measurements were well-fitted by a linear regression of $y=0.995\text{\hspace{0.17em}}x-3.6\text{\hspace{0.17em}}\text{\hspace{0.17em}}\phantom{\rule{0ex}{0ex}}\mathrm{deg}$ (${R}^{2}=0.997$).

We next tested the algorithm in a sample with depth-varying optical axis. A specimen was constructed by stacking two pieces of chicken tendons at an angle. The upper tendon sample was positioned at $\sim 72\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$, and the lower piece was positioned at about $-13\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$. In the intensity image [Fig. 2(a)], the boundary of the two pieces of tendon can be vaguely identified. The boundary disappeared in the cumulative retardance [Fig. 2(b)] and optical axis [Fig. 2(c)] images that show continuous banded patterns across the boundary. The image of local retardance [Fig. 2(d)] has a uniform appearance suggesting similar retardance in both tendon pieces. However, the two pieces show a distinct contrast in the image of local optical axis [Fig. 2(e)], while the axis appears to be relatively homogeneous within each piece. Figure 2(f) shows the histograms of local retardance and axis obtained from the two small ROI shown in Fig. 2(d) and 2(e). The calculated optical axis within the two ROIs in Fig. 2(f) is $69.5\pm 3.7\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ in the upper piece and $-16.0\pm 6.9\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ in the lower piece; both are consistent with the specimen configuration.

Figure 3 shows the *en face* images extracted along the two dotted lines in Fig. 2(a). Images were filtered by a $5\times 5$ median filter in the *en face* plane for display. Although the two tendon pieces can be discriminated in the structure image at depths where there was a clear separation [Fig. 3(a)], it is difficult to discriminate the thin layer of tissue with a different fiber orientation in the structure image at another depth [Fig. 3(f)]. Compared to other images, the local axis images [Fig. 3(e) and 3(j)] are the most effective in revealing fiber orientation changes at different imaging depths.

The algorithm was also applied to analyze PSOCT images of a plastic specimen. The intensity image [Fig. 4(a)] appears to be homogeneous with little structural variation. However, the cumulative retardance [Fig. 4(b)] shows the typical birefringent band structures. The inhomogeneous banded patterns in the cumulative optical axis [Fig. 4(c)] suggest a varying optical axis, which is clearly confirmed in the local optical axis [Fig. 4(e)]. Figure 4(f)–4(j) are the corresponding *en face* images extracted along the two dashed lines in Fig. 4(a). The images of local axis reveal some polarization features that are not discernible in the cumulative results or the structural images. For example, the “bow tie” pattern in Fig. 4(j2) suggests a twisted optical axis distribution in the *en face* plane that is gradually formed from a shallow depth [Fig. 4(j1)] to a deep depth [Fig. 4(j2)].

## 4.

## Conclusion

We demonstrated an algorithm that can successfully extract depth-resolved local optical axis from conventional PSOCT measurements. Different from the iterative method used previously in Jones matrix-based OCT systems,^{4}5.6.^{–}^{7} we used a reconstructed Jones matrix [Eq. (8)] in the calculation instead of the raw calculated Jones matrix. In addition, the local retardance was calculated using PSOCT results obtained in two consecutive layers.^{8} We believe these modifications helped to control error propagation in the calculation process. This proposed method can be applied to existing PSOCT systems that use a single circularly polarized incident light to map local optical axis and retardance in birefringent samples with negligible diattenuation. This approach can be also extended to Jones matrix-based systems as well to completely characterize samples with strong optical diattenuation.

## Acknowledgments

This project is supported in part by a National Science Foundation Grant CBET-0643190.