## Abstract

It is challenging to recover local optic axis orientation from samples probed with fiber-based polarization-sensitive optical coherence tomography (PS-OCT). In addition to the effect of preceding tissue layers, the transmission through fiber and system elements, and imperfect system alignment, need to be compensated. Here, we present a method to retrieve the required correction factors from measurements with depth-multiplexed PS-OCT, which accurately measures the full Jones matrix. The correction considers both retardation and diattenuation and is applied in the wavenumber domain, preserving the axial resolution of the system. The robustness of the method is validated by measuring a birefringence phantom with a misaligned system. Imaging *ex-vivo* lamb trachea and human bronchus demonstrates the utility of reconstructing the local optic axis orientation to assess smooth muscle, which is expected to be useful in the assessment of airway smooth muscle thickness in asthma, amongst other fiber-based applications.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Optical coherence tomography (OCT) measures the amount and path-length of light back-scattered by tissue to reconstruct cross-sectional views of the subsurface microstructure [1]. Many tissues feature similar scattering contrast, which can make their differentiation difficult, and would benefit from additional contrast. Various contrast-enhancing extensions of OCT have been proposed to date, such as optical coherence elastography [2, 3], OCT angiography [4, 5] and polarization-sensitive optical coherence tomography (PS-OCT) [6, 7]. PS-OCT provides additional contrast by measuring the polarization properties of the sample, including birefringence, degree of polarization (DOP), or optic axis orientation. The promise of polarization imaging has been demonstrated in various fields of medicine, including ophthalmology [8, 9], dentistry [10], dermatology [11-13], neurology [14, 15], cardiology [16] and pulmonology [17, 18].

Cumulative phase retardation, or retardance, and cumulative optic axis imaging are the two most straightforward contrasts available from PS-OCT, but since they do not directly provide local sample information, they are difficult to interpret, especially in samples with optic axis orientation that varies with depth. In addition, wrapping of the cumulative phase retardation causes artifacts in cumulative optic axis imaging [19]. It is more desirable, therefore, to assess birefringence and optic axis orientation locally at each depth for imaging samples, especially those composed of differently birefringent layers.

Local optic axis imaging has been demonstrated with free-space PS-OCT [20-24]. Recently, it has also been shown that local optic axis imaging is possible with fiber- and catheter-based PS-OCT employing non-polarization-maintaining optical fiber [25]. This has been achieved by compensating for the transmission through the system and fiber elements, estimated directly from the measurement data, owing to an intrinsic symmetry constraint when measuring polarization properties along identical illumination and detection paths. The specific PS-OCT implementation in that work used sequential illumination with two input polarization states orthogonal to each other on the Poincaré sphere, and relied on the assumption that the imaging system and the sample are free of diattenuation. Whereas dichroism is negligible compared with birefringence in most biological tissues, imperfect system components or inadequate alignment can lead to readily observed diattenuation that impacts the accurate recovery of the optic axis orientation. Further, the processing used therein employs spectral binning and, thereby, sacrifices axial resolution to mitigate polarization mode dispersion (PMD) present in the system components.

One of the recently investigated applications of optic axis orientation imaging is the assessment of airway smooth muscle (ASM) structure. Thickening (remodeling) and contraction of the ASM is considered a primary cause of excessive narrowing and symptoms in diseases such as asthma [26], and the capacity to clearly image the ASM could guide and improve patient therapy [27]. The idea behind imaging the optic axis orientation of the ASM relies on the organization of the airway wall, where the ASM is oriented approximately orthogonal to the surrounding tissue. To date, only one study on PS-OCT measurements of ASM has been presented by Adams *et al*. [18], which used the angle between the apparent optic axis orientation of the similarity transformed local Jones matrix and a reference axis without, however, correcting for the preceding layers. The extent to which such an assessment without depth-correction might compromise ASM assessment in more general airway geometries is as yet unknown.

In this manuscript, we present an improved method for local optic axis orientation imaging that corrects for general system imperfections and maintains the full axial resolution. It is based on a depth-multiplexed PS-OCT system that simultaneously measures the full Jones matrix and, therefore, does not rely on assumptions regarding the system and sample properties. The method is insensitive to system polarization distortions, including wavelength-dependent reference polarization states, and PMD in the system components. The full nominal axial resolution of the employed light source remains available because the corrections are directly applied in the wavenumber (*k*) domain. We validate our method on tissue-like birefringence phantoms, benchmark it against the previous reconstruction approach, and demonstrate imaging of a lamb trachea and human bronchus *ex vivo*.

## 2. Methods

We adapted a PS-OCT system previously used for deep tissue volumetric imaging with needle probes [28], where the needle probes were replaced with a standard scanning microscope configuration. The light source used is a wavelength-swept laser (Axsun Technologies, Bellerica, MA, USA), centered at 1310 nm, with a full sweep range of 100 nm at a 50 kHz sweep rate. 1% of the power of the source, split by a fiber coupler, was directed to a fiber Bragg grating and reflection from the grating was employed to trigger the acquisition and minimize timing jitter. The remaining power of the source, further split by a fiber coupler, was delivered to the sample (80%) and reference (20%) arms via conventional single-mode fibers, respectively. The sampling clock generated by the laser was electronically frequency-doubled to extend the total depth range from 5.0 mm to 10.0 mm. The increased imaging range enabled depth-encoded polarization multiplexing by passively delaying one of two orthogonal polarization states illuminating the sample with a free-space-based polarization delay unit (PDU) [29, 30]. An optical circulator delivered the light from the PDU to a standard scanning microscope configuration, which contains a fiber collimator (F220APC-1310, Thorlabs Inc., New Jersey, USA), a galvanometric x-y scanner (GVS002, Thorlabs Inc.) and a scan lens (LSM03, Thorlabs Inc.), providing approximately 25 μm lateral resolution. For lateral scanning, the distance between adjacent A-lines was adjusted to 5 μm for slight oversampling. In the reference arm, a fiber coupler was used to direct the light to the reference mirror and on to the receiver, imparting a 6 dB loss but avoiding potential polarization effects that could be caused by a circulator. The sample signal and reference signal were combined with a commercial polarization-diverse optical mixer (PDOM-1310, Finisar, Sunnyvale, USA) and detected with a pair of balanced receivers (PDB460C-AC, Thorlabs Inc.), connected to a high-speed dual-channel digitizer (ATS9350, Alazar Technologies Inc., Pointe-Claire, Québec, Canada). The fiberized polarization-diverse optical mixer lacks a linear polarizer on the reference arm channel before interfering the two signals, relying instead on the polarization state of the reference light itself to equally split between the two polarization-diverse channels.

The full Jones matrix of the complex-valued interference signals, **J**_{tot}(*k*), where *k* is the wavenumber, was reconstructed using the steps as follows: subtraction of the recorded background signal; numerical compensation of the chromatic dispersion in the system; Fourier transformation to *z*-domain; compensation of sensitivity roll-off along depth; cropping and coarse alignment of depth-multiplexed signals with integer pixel number; inverse Fourier transformation to the *k*-domain; and precise sub-pixel alignment of the depth-delayed signals by applying a pre-calibrated linear phase ramp in the *k*-domain.

## 3. Algorithm

#### 3.1 Symmetrization of the full Jones matrix

**J**_{tot}(*k*), retrieved from measurements with our PS-OCT system, can be modeled as:

**J**_{in}(*k*) and **J**_{out}(*k*) are *k*-dependent transmission matrices describing the illumination optics (including the PDU), and the detection optics (comprising the effect of non-ideal reference signal states), respectively [31]. **J**_{sam}(*z*) and ${\mathbf{J}}_{\text{sam}}^{T}(z)$ correspond to the propagation from the sample surface to a given depth *z* and back to the surface, respectively, excluding the global phase, which is explicitly specified by exp(*i*2*kz*). The complex-valued Jones matrix tomogram reconstructed from **J**_{tot}(*k*) can be described as:

*w*(

*k*,

*p*), spanning 1/

*N*

^{th}of the full spectral width, Δ

*k*, and moved across the spectrum in equal steps, indexed by

*p*∈ [1, 2

*N*−1],

*N*=8, resulting in 15 spectral bins [32]. The spectral variation of

**J**

_{in}(

*k*) and

**J**

_{out}(

*k*) is sufficiently slow to be assumed constant within each spectral bin, approximated by the central wavenumber of each bin

*k*=

*k*

_{0}+

*p*Δ

*k*/(2

*N*), where

*k*

_{0}indicates the smallest recorded wavenumber. Alternatively, a single Hamming window spanning the entire range is used to reconstruct the full-resolution Jones matrix tomogram

**J**

_{tom}(

*z*) =

*FT*[

**J**

_{tot}(

*k*)].

In a system with symmetric round-trip transmission (i.e., identical illumination path and detection path), both **J**_{tot}(*k*) and **J**_{tom}(*z,p*) would exhibit transpose symmetry [33]. However, because of the generally distinct fiber components in the illumination and detection paths, the retrieved matrices are unlikely to be transpose-symmetric. Yet, the transpose symmetry can be recovered numerically by multiplication on the left-hand side of Eq. (1) with a *k*-dependent correction matrix **J**_{cor}(*k*):

In order to estimate **J**_{cor}(*k*), we first compute the correction matrix **J**_{cor}(*p*) for each spectral bin, following the method described by Villiger *et al.* [25]. **J**_{cor}(*p*) minimizes the sum of the squared Euclidean norm of the difference between the off-diagonal entries of the recovered Jones matrices across an entire B-scan:

**J**

_{cor}(

*p*)) = +1. All points with sufficient signal (SNR > 20 dB) within an entire B-scan are used for this estimation, which can be cast as a matrix eigenvalue problem [25]. In brief, Eq. (5) can be rewritten as a quadratic problem ${\mathbf{j}}_{\text{cor}}^{\u2020}\cdot \mathbf{H}\cdot {\mathbf{j}}_{\text{cor}}$, where

**j**

_{cor}is the vectorized Jones matrix

**J**

_{cor}, the dagger indicates complex conjugation, and

**H**consists of permutations of the sum over the complex conjugate of the inner products of the vectorized measured Jones matrices ${\mathbf{j}}_{\text{tom}}(z)\cdot {\mathbf{j}}_{\text{tom}}^{\u2020}(z)$. The eigenvector of

**H**associated with the smallest eigenvalue, scaled to a determinant of +1, provides the solution to

**J**

_{cor}. Hence, Eq. (5) can be solved directly with an eigenvalue decomposition of the 4 × 4 matrix

**H**for each spectral bin

*p*. Crucially, the recovered

**J**

_{cor}(

*p*) is a general Jones matrix, featuring both retardation and diattenuation, since

**J**

_{tom}(

*z,p*) is not restricted to a pure retardation matrix as in previous work [25].

A transpose-symmetric retardation matrix corresponds to a linear retarder and a transpose-symmetric diattenuation matrix to a linear diattenuator. The effects of retardation and diattenuation are best visualized as vectors using the Stokes formalism and the *Q*, *U*, and *V* coordinates of the Poincaré sphere. The direction of the retardation vector corresponds to the slow optic axis and its length to the amount of retardance. Likewise, the direction of the diattenuation vector corresponds to the polarization state experiencing weakest diattenuation, and its length to the amount of diattenuation between the state opposite on the Poincaré sphere and this state. A linear retarder/diattenuator corresponds to a retardation/diattenuation vector confined to the horizontal *QU* plane of the Poincaré sphere. Alignment of the vector along *V* is referred to as circular retardation/diattenuation. Using the common polar decomposition [34] on a general (both retarding and diattenuating) transpose-symmetric Jones matrix would obscure the underlying symmetry and produce a diattenuation vector no longer confined to the *QU* plane. Instead, we here use simultaneous, or concurrent, decomposition (and synthesis) of a Jones matrix, explained in Appendix A, which preserves this symmetry. Applying this concurrent decomposition to **J**_{cor}(*p*) provides the retardation vector **r**_{cor}(*p*) and diattenuation vector **d**_{cor}(*p*). Using 3^{rd} order polynomial fitting of the discrete **r**_{cor}(*p*) and **d**_{cor}(*p*) vectors, we obtain the continuous variables **r**_{cor}(*k*) and **d**_{cor}(*k*) across the full spectrum, allowing us to reconstruct **J**_{cor}(*k*) (according to Eq. (10) in the Appendix), and avoiding any trade-off with axial resolution as in most previous forms of spectral binning, similar to the method of Braaf *et al.* [31]. Of note, fitting the retardation **r**_{cor}(*p*) and diattenuation **d**_{cor}(*p*) vectors, rather than the entries of the Jones matrix **J**_{cor}(*p*), preserves the determinant of the correction matrix **J**_{cor}(*k*), which is always +1. Because **J**_{cor}(*k*) is stable over multiple image acquisitions, once it is estimated from a single B-scan, it can be directly applied to all remaining B-scans in a three-dimensional dataset.

To illustrate the symmetrization process, we use PS-OCT measurements of an *ex-vivo* trachea sample obtained from a fetal lamb (Fig. 1). In the B-scan intensity image, one can distinguish the mucosa, smooth muscle, and cartilage layer (Fig. 1(a)). The retardation vectors **r**_{cor}(*p*) and diattenuation vectors **d**_{cor}(*p*) of **J**_{cor}(*p*) (*p* ∈[1, 15]) and their polynomial fits across the full spectrum are shown in Fig. 1(b). We also decomposed **J**_{tom}(*z*) and ${\mathbf{J}}_{\text{tom}}^{\prime}(z)$ into retardation and diattenuation vectors using the concurrent decomposition to assess their symmetry. The normalized histograms of the circular (*V*) component of the retardation (Fig. 1(c)) and diattenuation (Fig. 1(d)) vectors before (red) and after (blue) symmetrization, evidence the almost complete cancellation of the circular components of both retardation and diattenuation. This response is consistent with that of the concurrent action of a linear retarder and a linear diattenuator, which are confined to the *QU* plane [35].

#### 3.2 Compensation of system polarization distortions

Although symmetric, ${\mathbf{J}}_{\text{tot}}^{\prime}(k)$ may still exhibit *k*-dependence due to system imperfections, such as PMD, that have to be corrected with a compensation matrix **J**_{comp}(*k*) [25]:

*k*

_{c}the central wavenumber and

**J**

_{comp}(

*k*) =

**J**

_{in}

^{−1}(

*k*)·

**J**

_{in}(

*k*

_{c}). In analogy to the symmetrization of the Jones matrices, we find the compensation matrix

**J**

_{comp}(

*p*) for each spectral bin ${\mathbf{J}}_{\text{tom}}^{\prime \text{}\prime}(z,p)$, and then use 3

^{rd}order polynomial fitting on the concurrently decomposed retardation and diattenuation vectors to get continuous

**r**

_{comp}(

*k*) and

**d**

_{comp}(

*k*) and reconstruct

**J**

_{comp}(

*k*) across the full spectrum.

**J**

_{comp}(

*p*) is found with an iterative method, by computationally minimizing the sum of the squared Euclidean norm of the difference between the

*p*th bin and the central bin across all points within an entire B-scan with sufficient signal (SNR > 20 dB):

**J**

_{comp}(

*p*), it features by construction a determinant of +1, and is defined by its retardation and diattenuation vectors

**r**

_{comp}(

*p*) and

**d**

_{comp}(

*p*). The minimization problem is, hence, simplified to iteratively searching for the six real-valued variables defining

**r**

_{comp}(

*p*) and

**d**

_{comp}(

*p*), rather than the four complex-valued entries of

**J**

_{comp}(

*p*), for each spectral bin

*p*. The iterative search begins with initial values for

**r**

_{comp}(

*p*) and

**d**

_{comp}(

*p*) of [0, 0, 0]

^{T.}The *Q*, *U* and *V* components of the retardation **r**_{comp}(*p*) and diattenuation **d**_{comp}(*p*) vectors of the compensation matrix **J**_{comp}(*p*) and their 3^{rd} order polynomial fits **r**_{comp}(*k*) and **d**_{comp}(*k*) from the same B-scan in Fig. 1 are presented in Figs. 2(a) and 2(b), respectively. Fig. 2(c) shows the average of the squared Euclidean norm of the difference between each bin and the central bin before and after compensation for system polarization distortions..

#### 3.3 Extraction of relative local optic axis orientation

After symmetrization and compensation for *k*-dependence of the system components, ${\mathbf{J}}_{\text{tot}}^{\prime \text{}\prime}(k)$ is expected to be equivalent to Jones matrix measurements with a bulk-optic PS-OCT system, except for the presence of the unknown **J**_{in}(*k*_{c}). Using the full spectral width without any binning, we reconstruct the full-resolution Jones matrix tomogram ${\mathbf{J}}_{\text{tot}}^{\prime \text{}\prime}(z)$ and decompose it into retardation and diattenuation vectors using the concurrent decomposition technique. Any residual circular components along the *V*-direction are set to zero to ensure exact transpose symmetry. Figs. 3(a) and 3(b) display the cumulative round-trip retardation and diattenuation corresponding to the length of the recovered retardation and diattenuation vectors. The diattenuation components are relatively small compared to the retardation and appear dominated by noise. This suggests that both **J**_{in}(*k*_{c}) and the sample feature little diattenuation. To simplify the further processing, we followed the strategy of Fan *et al.* of synthesizing pure retardation matrices and ignoring the diattenuation components [20]. In fact, this corresponds to the least squares solution of approximating a general Jones matrix with a pure retardation matrix. But rather than using SU(2) Jones matrices, we construct SO(3) rotation matrices from the retardation vectors [36] to enable effective incoherent spatial averaging, similar to the averaging of Mueller matrices [28]. Moreover, conversion to SO(3) removes the sign ambiguity of the retardation vector in ${\mathbf{J}}_{\text{tot}}^{\prime \text{}\prime}(z)$. To suppress speckle, the SO(3) matrices are weighted with the tomogram intensity signal (the magnitude of the determinant of the Jones matrix ${\mathbf{J}}_{\text{tot}}^{\prime \text{}\prime}(z)$ before scaling [37]) and then spatially averaged by a three-dimensional Gaussian kernel (FWHM approximately twice the speckle width in both lateral and axial directions), using the same formulation and computationally efficient implementation of the Euclidean mean described previously [25]. In short, the averaged rotation matrices contain depolarization components that have to be removed in order to remain in SO(3). We used three-dimensional averaging throughout this work, but two-dimensional filtering may be necessary for *in vivo* imaging and the presence of motion artifacts. The effect on the retardance of spatially filtered SO(3) rotation matrices is shown in Fig. 3(c). In addition, retardance from the central A-lines of the yellow and green boxes before (Fig. 3(a)) and after (Fig. 3(c)) spatial filtering are plotted in Figs. 3(d) and 3(e), respectively.

The recovered SO(3) retardation matrices after spatial filtering, denoted as **R**_{rt}(*z*), express the retardation accumulated during the round trip to depth *z* and back. If **R**_{st}(*z*) is the single trip retardation matrix to depth *z*, we have ${\mathbf{R}}_{\text{rt}}(z)=\mathbf{D}\cdot {\mathbf{R}}_{\text{st}}^{T}(z)\cdot \mathbf{D}\cdot {\mathbf{R}}_{\text{st}}(z)$, where **D** is a diagonal matrix diag(1, 1, −1). In SO(3), **D**·**M*** ^{T}*·

**D**corresponds to the reverse transmission through element

**M**, in analogy to the simple transpose of Jones matrices [33].

**R**

_{st}(

*z*) can be described as the sequential transmission through individual tissue layers

**Q**

_{n}of thickness

*dz*, the axial sampling distance:

Each layer is assumed to act as a linear retarder. Its retardation vector defines the local retardation and local optic axis orientation that we are attempting to recover. The sequence of linear retarders defining **R**_{st}(*z*) results in a general retarder. However, we only measure the round-trip matrix **R**_{rt}(*z*), which cancels the circular retardation component. Hence, we follow the recursive reconstruction of Fan *et al.* [20] to isolate the round-trip transmission through each individual layer and compute its square root to recover the single-pass transmission and construct **R**_{st}(*z*) for the subsequent layer, starting at *z* = 0:

**M**

*=*

^{T}**M**

^{−1}in SO(3). At pixel locations with insufficient signal (SNR < 5dB),

**Q**is replaced with the identity matrix. Finally, the retardation vector of

**Q**is retrieved and axially filtered with a rectangular filter equal in width to the axial FWHM of the spatial filter previously used to average the SO(3) retardation matrices. The length and orientation of this filtered vector provide the local retardation and the local optic axis orientation, respectively.

Although we confirmed that **J**_{in}(*k*c) is minimally diattenuating, its retardation component still impacts **R**_{rt}(*z*) and defines **R**_{rt}(*z* = 0). Taking its matrix root to compute **Q**_{0} only recovers the linear retardation of **J**_{in}(*k*_{c}). The non-detectable circular retarding component in **J**_{in}(*k*_{c}) offsets the recovered optic axis orientations of the sample, making the optic axis measurement relative.

## 4. Results

#### 4.1 Validation with birefringence phantoms

We validated the method using a custom-made birefringence phantom (Sample #1). We cut thin strips from a polylactic acid (PLA) 3D printer filament. The PLA filament exhibits intrinsic birefringence with its fast optic axis parallel to the filament’s axis [38]. The strips were prepared so the longer edge is parallel to the optic axis of the filament. We arranged multiple thin strips at different angles with respect to the laboratory frame and used adhesive tape to fix the assembly. No contrast between strips is visible for the OCT intensity, local phase retardation, and depolarization index images (Figs. 4(a-c)), where the depolarization index was computed following Lippok *et al*. [39]. The local optic axis visualization, however, strongly correlates with the orientation of each strip (Fig. 4(d)). Note that since the measured optic axis orientation is relative to an unknown laboratory frame, an offset angle is manually applied here to the measured optic axis orientation to align it with the optic axis wheel, purely for ease of comparison.

We next performed experiments to demonstrate that the local optic axis orientation measured at greater depths is not affected by the optic axis orientation of the preceding layers. We made an additional PLA phantom (Sample #2) with two layers of PLA strips (total thickness ≈2 mm) immersed in ultrasound gel and attached to a piece of adhesive tape, and generated *en face* and cross-sectional local optic axis orientation images (Figs. 4(e-i)). In the top (Fig. 4(e)) and bottom (Fig. 4(f)) layers, multiple strips are oriented at different angles. Figure 4(g) and Fig. 4(i) are cross-sectional images of the green and pink dashed lines in Fig. 4(e) and Fig. 4(f), respectively. Figure 4(e) and Fig. 4(f) correspond to *en face* views of the red and blue dashed layers in Fig. 4(g) and Fig. 4(i), respectively. Figure 4(h) is a snapshot of the 3D rendering of the optic axis orientation overlaid on the intensity of the phantom (Visualization 1) by ImageJ [40]. These results demonstrate that our methodology can provide reliable local optics axis orientation without the influence of the preceding layers.

We next show here that, through symmetrization and compensation of system PMD, our method becomes robust to variations in the system. We do this by deliberately introducing distortions into the system. We repeatedly measured the same PLA phantom (Figs. 5(a-d)), whilst varying the distortions in between measurements by adjusting the polarization controller in the reference arm (Fig. 5(a)). Because the polarization mixer lacks a linear polarizer on the reference input, varying the reference polarization state is akin to introducing system diattenuation [31]. To characterize the distortion, we analyzed the intensity in each reference polarization state by summing the square of the magnitude of the demodulated complex-valued fringe signals (dominated by the reference signal) over the two input states. The optic axis orientations remain almost unchanged regardless of the distortions (Fig. 5(b)). In addition, a benchmark against the method published by Villiger *et al.* [25] is also performed, the results of which (Fig. 5(c)) are expected to be suboptimal due to its neglect of system diattenuation. Briefly, in Villiger *et al.* [25], the imaging system uses two sequential probing Stokes vectors orthogonal to each other on the Poincaré sphere, and the processing assumes pure retardation of the sample and the system. To ease the comparison with a common reference, the measured relative optic axis orientations in Figs. 5(b) and (c) are individually offset, such that the averaged orientations of the dashed green boxes are 6°, matching with their physical orientations and, thus, in effect, converting the relative OA into an absolute OA. Since there is an intrinsic π ambiguity with the optic axis orientation, we confine the absolute optic axis between 0° and 180°, and the acute and obtuse mean value and the standard deviation of the absolute OA of three ROIs (solid boxes in (b) and (c)) are plotted in Figs. 5 (d) and (e), respectively. Apparently, by means of the method described in this manuscript, the absolute OA shows more consistency across different distortions of the reference signal, suggesting better robustness of this method.

#### 4.2 Tissue imaging

In addition to phantom measurements, we performed a pilot study with biological samples. Fetal (128 d, term = 150 d) lamb trachea was excised, incised at the anterior surface and laid flat with the luminal surface exposed for imaging. Animal tissue was sourced via approved institutional tissue sharing processes. Two separate trachea samples (Sample #1 and Sample #2) were scanned by PS-OCT, and Sample #1 was scanned twice over different regions (Region 1A and Region 1B). Cross-sectional images were taken perpendicular to the trachea longitudinal dimension. In Fig. 6, cross-sectional intensity images from different scans are shown in the first column (Figs. 6(a), (d), (g)), local phase retardation images, overlaid on intensity, are shown in the second column (Figs. 6(b), (e), (h)), and local optic axis orientation images overlaid on intensity are shown in the third column (Figs. 6(c), (f), (i)). Some contrast between tracheal smooth muscle and the surrounding tissue is visible on the OCT intensity images for Region 1A (Figs. 6(a-c)) and Region 1B (Figs. 6(d-f)) of Sample #1, but barely visible for Sample #2 (Figs. 6(g-i)). However, when optic axis orientation is considered, we observe excellent contrast for all regions, clearly delineating muscle (green) from the inner mucosa (purple).

Some indication of the orientation of the fibrous structures in different layers of the lamb trachea is observable in the OCT intensity images when *en face* projection (without flattening in post-processing) is used (Figs. 7(a), (d), (g)). Local phase retardation mapping offers rather homogeneous images of the same *en face* layers (Figs. 7(b), (e), (h)); whereas, local optic axis mapping provides excellent contrast (Figs. 7(c), (f), (i)).

We also imaged *ex-vivo* an airway specimen (bronchus) obtained from an adult smoker with normal lung function. Tissue was collected after clinically-indicated lung resection surgery and was approved by the Department of Health Human Ethics Research Committee. Following PS-OCT imaging, the sample was fixed in formalin, embedded in paraffin, sectioned (5 μm) and stained with hematoxylin and eosin (H&E). The OCT intensity image (Fig. 8(a)) provides no contrast between different layers at all. The local phase retardation image (Fig. 8(b)) shows higher birefringence in some regions, however, it remains challenging to tell apart the ASM layer from the mucosa. In comparison, the local relative optic axis orientation image (Fig. 8(c)) clearly contrasts the ASM layer from the mucosa and other tissue layers and shows good correlation with the H&E histology (Fig. 8(g)), where the ASM layer is marked with arrows. Figures 8(d-f) and (g-i) are, respectively, the *en face* views of intensity, local phase retardation and optic axis at approximately 180 μm and 420 μm deep into the sample from the tissue surface. Surface flattening is used purely for ease of interpretation, since the sample is curved in geometry. The dashed red lines in Figs. 8(d-f) indicate the location of the B-scan in Figs. 8(a-c).

## 5. Discussion

Polarized light microscopy (PLM) has long been used to reveal precious sample contrast by leveraging intrinsic sample properties [41, 42]. PLM uses carefully designed and aligned optical components to manage the polarization states of light. In comparison, fiber-based PSOCT employs components that alter the polarization state of the light and act differently across the employed spectrum. The method presented in this manuscript estimates wavelength-dependent correction terms, thereby improving on these deficiencies. The estimation is performed directly on sample data, without specific requirements on sample properties and without the need for specific calibration signals.

Comparisons of the optic axis orientation at three distinct regions of a birefringence phantom suggest better robustness of the method in this manuscript (red in Figs. 5(d-e)) than the one in Villiger *et al.* [25] (blue in Figs. 5(d-e)). The advantages of the approach proposed here are: i) the corrections are made in the *k*-domain, so the axial resolution of the system is preserved, which is especially important in resolving thin structures of thickness close to the OCT axial resolution; ii) the diattenuation of the system is taken into account and, therefore, the errors introduced by ignoring diattenuation are avoided and the effort of aligning the reference signal is minimized; iii) unlike those systems that require probing beam(s) with circular polarization states, the method requires no alignment of input polarization states on the sample, and its extraction of the relevant polarization parameters is robust, which might be beneficial for future clinical applications, where careful alignment procedures are impractical.

According to Park *et al.* [35], for fiber-based PS-OCT, there exists an inherent ambiguity in the sign of the relative optic axis orientation. Although there is a unique solution to the symmetrization of the polarization measurements, there remains an ambiguity in the sign of the matrix root when computing **Q**_{0}, which corresponds to a ‘flip’ of the *QU*-plane of the Poincaré sphere. A straightforward calibration step allows removing this sign ambiguity, whereby a birefringence phantom with at least two distinct, non-orthogonal, and a *priori* known optic axis orientations is imaged. Matching the sense of the relative axis orientations readily provides the required sign to correctly interpret the axis orientations, as illustrated in Fig. 4, and remains consistent until the system retardance is altered.

Once the sign ambiguity is removed, it is possible to convert the relative local optic axis orientation, as measured here, to absolute local optic axis orientation if the orientation of any layer is known and, thus, serves as a reference for the entire sample, as in Fig. 5. Yet, in many cases, the relative optic axis orientation provides sufficient information to differentiate layers of the sample and renders unnecessary the need for absolute optic axis orientation. Additionally, we found that in our system the relative orientation of the optic axis remained unchanged with respect to the laboratory frame, even over periods of several days, as long as the system components and fibers remained untouched. Therefore, we anticipate that calibration with a polarization element once per measurement session may well be sufficient to provide a measure of absolute optic axis orientation.

The optic axis identified by our processing corresponds to a right-handed rotation in the Poincaré space with right-hand circularly polarized light (looking against the source) located at [*Q*, *U*, *V*]* ^{T}* = [0, 0, 1]

*, and identifies the slow, i.e., retarding, axis of the birefringent material [36]. Biological birefringent tissues, such as collagen and muscle, feature positive birefringence and, hence, the slow axis corresponds to their extraordinary refractive index aligned along their fiber direction [43]. Optic axis imaging in such tissues, as in Figs. 5–8, provides excellent contrast between layers with different orientations. Local phase retardation imaging, however, might fail to contrast such layers if they have similar birefringence. Additionally, tensile or compressional stress of the sample, occurring in preparation or*

^{T}*in situ*, might influence the amount of birefringence with unknown effect, but should, in principle, not change the local optic axis orientation of any layer in the sample. For tissues without birefringence, the local SO(3) becomes close to the identity matrix, whose optic axis orientation is undefined. Practically, the residual noise ensures that the resulting optic axis orientation is, indeed, randomly oriented. The higher the birefringence of a sample, the more trustworthy is the measured optic axis orientation. Hence, it is possible to use the local phase retardation as a mask to suppress regions with low/no birefringence. However, our samples were sufficiently retarding throughout that we did not use this strategy. In Fig. 7, the inner mucosa has, indeed, lower birefringence than the smooth muscle, but still features sufficient birefringence to produce reliable optic axis maps.

With the current unoptimized code running on MATLAB R2016b, the processing time is on average 2.0 s per B-scan (consisting of 1000 A-lines and computing 1.2 mm in depth) on a desktop equipped with an Intel CPU (i7-7700 @3.6GHz) and with 32GB RAM. The symmetrization correction matrix **J**_{cor}(*k*) and system polarization distortion compensation matrix **J**_{comp}(*k*) are both estimated from only a single B-scan at the beginning of the processing and directly applied to all the other B-scans in the data set. It takes approximately 19 s to find **J**_{cor}(*k*) and **J**_{comp}(*k*). Total time to process a 3D data set (1000 × 1000 A-lines) currently is approximately 30 min for computing 1.2 mm in depth and 63 min for computing 5.0 mm in depth. To maximize the computing efficiency, we first crop the data volume based on its 3D intensity, and then compute the optic axis only to a certain depth, which is typically less than 2.0 mm, depending on the thickness of the measured sample. We expect that the reconstruction could be significantly accelerated by improving the numerical implementation and using parallel processing with a graphics card.

## 6. Conclusions

We have developed a robust method for local phase retardation and local relative optic axis orientation imaging that corrects for general system imperfections and maintains the full axial resolution, based on a depth-multiplexed PS-OCT system that simultaneously measures the full Jones matrix. We validated this method with birefringence phantoms and demonstrated imaging of lamb trachea and human airway samples *ex vivo*. We demonstrated it as a promising method for segmenting ASM from other wall structures, which, with further validation, may be a useful tool to study ASM remodeling in asthma.

## Appendix

## Concurrent decomposition of a general Jones matrix into simultaneously acting retardation and diattenuation vectors

A general Jones matrix is defined, to within an arbitrary phase and scaling term, as:

*n*= 1, 2, 3. Real-valued

*r*

_{1},

*r*

_{2}and

*r*

_{3}correspond, respectively, to the

*Q*,

*U*,

*V*components of a retardation vector

**r**= [

*r*

_{1},

*r*

_{2},

*r*

_{3}]

*. Similarly, real-valued*

^{T}*d*

_{1},

*d*

_{2}and

*d*

_{3}correspond, respectively, to a diattenuation vector

**d**= [

*d*

_{1},

*d*

_{2},

*d*

_{3}]

*. σ*

^{T}*is the Pauli basis (*

_{n}*n*=0, 1, 2, 3):

This decomposition assumes that retardation and diattenuation act simultaneously, or concurrently, in contrast to the polar decomposition [44], which groups retardation and diattenuation into a sequential order. The matrix exponential in Eq. (10) above can be expanded in terms of *f _{n} =* (

*d*)/2 into the series:

_{n}− ir_{n}**I**is the identity matrix. Owing to the anticommutator properties of the Pauli basis {σ

*, σ*

_{m}*} = 2δ*

_{n}

_{mn}**I**, where δ

*is the Kronecker delta, the even powers of the above series simplify to:*

_{mn}Identifying the series expansions of the cosine hyperbolic and the sine hyperbolic allows rewriting Eq. (12) as:

To retrieve the *r _{n}* and

*d*parameters from a given Jones matrix, we use the trace property

_{n}*Tr*(σ

*·σ*

_{m}*) = 2δ*

_{n}*of the Pauli basis and define:*

_{mn}*q*

_{0}= 2cosh(

*c*) and finally leads to:

In Eq. (15), it is assumed that det(**J**) = +1, which is the case if parameterized as indicated in Eq. (10), because the determinant of a matrix exponential is identical to the exponential of the matrix trace, which is zero in this case. However, to decompose a general Jones matrix **J**′, it needs to be divided by the square-root of its determinant. To limit complications with the ambiguous sign of the square-root, we rotate the arguments of **J**′ such that the phase of the first entry *J*_{11} of **J**′ is zero, similar to the method of Li *et al.* [17]:

Whereas traditional polar decomposition relies on computationally demanding matrix operations, including eigenvalue decomposition, the concurrent decomposition, which assumes the concurrent effect of retardation and diattenuation, offers computationally efficient implementation and is straightforward to vectorize. Crucially, decomposing a transpose-symmetric Jones matrix in the concurrent fashion results in both retardation and diattenuation vectors that are confined to the *QU*-plane. Polar decomposition obscures this symmetry for diattenuation.

## Funding

Australian Research Council (DP160104969); National Institutes of Health (P41EB-015903).

## Acknowledgments

Q. Li would like to acknowledge the financial support of the Australian Government International Research Training Program (RTP) Fee Offset Scholarships.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

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

**2. **K. V. Larin and D. D. Sampson, “Optical coherence elastography - OCT at work in tissue biomechanics [Invited],” Biomed. Opt. Express **8**(2), 1172–1202 (2017). [CrossRef] [PubMed]

**3. **B. F. Kennedy, P. Wijesinghe, and D. D. Sampson, “The emergence of optical elastography in biomedicine,” Nat. Photonics **11**(4), 215–221 (2017). [CrossRef]

**4. **C.-L. Chen and R. K. Wang, “Optical coherence tomography based angiography [Invited],” Biomed. Opt. Express **8**(2), 1056–1082 (2017). [CrossRef] [PubMed]

**5. **P. Gong, S. Es’haghian, K. A. Harms, A. Murray, S. Rea, B. F. Kennedy, F. M. Wood, D. D. Sampson, and R. A. McLaughlin, “Optical coherence tomography for longitudinal monitoring of vasculature in scars treated with laser fractionation,” J. Biophotonics **9**(6), 626–636 (2016). [CrossRef]

**6. **J. F. de Boer, C. K. Hitzenberger, and Y. Yasuno, “Polarization sensitive optical coherence tomography - a review [Invited],” Biomed. Opt. Express **8**(3), 1838–1873 (2017). [CrossRef] [PubMed]

**7. **B. Baumann, “Polarization sensitive optical coherence tomography: a review of technology and applications,” Appl. Sci. (Basel) **7**(5), 474 (2017). [CrossRef]

**8. **B. Cense, T. C. Chen, B. H. Park, M. C. Pierce, and J. F. de Boer, “In vivo depth-resolved birefringence measurements of the human retinal nerve fiber layer by polarization-sensitive optical coherence tomography,” Opt. Lett. **27**(18), 1610–1612 (2002). [CrossRef]

**9. **M. Pircher, C. K. Hitzenberger, and U. Schmidt-Erfurth, “Polarization sensitive optical coherence tomography in the human eye,” Prog. Retin. Eye Res. **30**(6), 431–451 (2011). [CrossRef] [PubMed]

**10. **T. Louie, C. Lee, D. Hsu, K. Hirasuna, S. Manesh, M. Staninec, C. L. Darling, and D. Fried, “Clinical assessment of early tooth demineralization using polarization sensitive optical coherence tomography,” Lasers Surg. Med. **42**(10), 738–745 (2010). [CrossRef]

**11. **W. C. Lo, M. Villiger, A. Golberg, G. F. Broelsch, S. Khan, C. G. Lian, W. G. Austen Jr., M. Yarmush, and B. E. Bouma, “Longitudinal, 3D imaging of collagen remodeling in murine hypertrophic scars in vivo using polarization-sensitive optical frequency domain imaging,” J. Invest. Dermatol. **136**(1), 84–92 (2016). [CrossRef] [PubMed]

**12. **P. Gong, L. Chin, S. Es’haghian, Y. M. Liew, F. M. Wood, D. D. Sampson, and R. A. McLaughlin, “Imaging of skin birefringence for human scar assessment using polarization-sensitive optical coherence tomography aided by vascular masking,” J. Biomed. Opt. **19**(12), 126014 (2014). [CrossRef] [PubMed]

**13. **E. Li, S. Makita, Y.-J. Hong, D. Kasaragod, and Y. Yasuno, “Three-dimensional multi-contrast imaging of in vivo human skin by Jones matrix optical coherence tomography,” Biomed. Opt. Express **8**(3), 1290–1305 (2017). [CrossRef] [PubMed]

**14. **H. Wang, A. J. Black, J. Zhu, T. W. Stigen, M. K. Al-Qaisi, T. I. Netoff, A. Abosch, and T. Akkin, “Reconstructing micrometer-scale fiber pathways in the brain: multi-contrast optical coherence tomography based tractography,” Neuroimage **58**(4), 984–992 (2011). [CrossRef] [PubMed]

**15. **H. Wang, T. Akkin, C. Magnain, R. Wang, J. Dubb, W. J. Kostis, M. A. Yaseen, A. Cramer, S. Sakadžić, and D. Boas, “Polarization sensitive optical coherence microscopy for brain imaging,” Opt. Lett. **41**(10), 2213–2216 (2016). [CrossRef] [PubMed]

**16. **M. Villiger, K. Otsuka, A. Karanasos, P. Doradla, J. Ren, N. Lippok, M. Shishkov, J. Daemen, R. Diletti, R.-J. van Geuns, F. Zijlstra, G. van Soest, P. Libby, E. Regar, S. K. Nadkarni, and B. E. Bouma, “Coronary plaque microstructure and composition modify optical polarization: a new endogenous contrast mechanism for optical frequency domain imaging,” JACC: Cardiovasc. Imag., Dec., epub ahead of print (2017).

**17. **J. Li, F. Feroldi, J. de Lange, J. M. A. Daniels, K. Grünberg, and J. F. de Boer, “Polarization sensitive optical frequency domain imaging system for endobronchial imaging,” Opt. Express **23**(3), 3390–3402 (2015). [CrossRef] [PubMed]

**18. **D. C. Adams, L. P. Hariri, A. J. Miller, Y. Wang, J. L. Cho, M. Villiger, J. A. Holz, M. V. Szabari, D. L. Hamilos, R. Scott Harris, J. W. Griffith, B. E. Bouma, A. D. Luster, B. D. Medoff, and M. J. Suter, “Birefringence microscopy platform for assessing airway smooth muscle structure and function in vivo,” Sci. Transl. Med. **8**(359), 359ra131 (2016). [CrossRef] [PubMed]

**19. **C. Fan and G. Yao, “Mapping local retardance in birefringent samples using polarization sensitive optical coherence tomography,” Opt. Lett. **37**(9), 1415–1417 (2012). [CrossRef] [PubMed]

**20. **C. Fan and G. Yao, “Imaging myocardial fiber orientation using polarization sensitive optical coherence tomography”,” BiomedOpt. Express **4**(3), 460–465 (2013). [CrossRef]

**21. **M. Todorović, S. Jiao, L. V. Wang, and G. Stoica, “Determination of local polarization properties of biological samples in the presence of diattenuation by use of Mueller optical coherence tomography,” Opt. Lett. **29**(20), 2402–2404 (2004). [CrossRef]

**22. **L. Azinfar, M. Ravanfar, Y. Wang, K. Zhang, D. Duan, and G. Yao, “High resolution imaging of the fibrous microstructure in bovine common carotid artery using optical polarization tractography,” J. Biophotonics **10**(2), 231–241 (2017). [CrossRef]

**23. **Y. Wang, K. Zhang, N. B. Wasala, D. Duan, and G. Yao, “Optical polarization tractography revealed significant fiber disarray in skeletal muscles of a mouse model for Duchenne muscular dystrophy,” Biomed. Opt. Express **6**(2), 347–352 (2015). [CrossRef] [PubMed]

**24. **Y. Wang, M. Ravanfar, K. Zhang, D. Duan, and G. Yao, “Mapping 3D fiber orientation in tissue using dual-angle optical polarization tractography,” Biomed. Opt. Express **7**(10), 3855–3870 (2016). [CrossRef] [PubMed]

**25. **M. Villiger, B. Braaf, N. Lippok, K. Otsuka, S. K. Nadkarni, and B. E. Bouma, “Optic axis mapping with catheter-based polarization-sensitive optical coherence tomography, ” Optica (accepted).

**26. **Y. S. Prakash, “Emerging concepts in smooth muscle contributions to airway structure and function: implications for health and disease,” Am. J. Physiol. Lung Cell. Mol. Physiol. **311**(6), L1113–L1140 (2016). [CrossRef] [PubMed]

**27. **G. M. Donovan, J. G. Elliot, F. H. Y. Green, A. L. James, and P. B. Noble, “Unravelling a clinical paradox - why does bronchial thermoplasty work in asthma?” Am. J. Respir. Cell Mol. Biol. **59**(3), 355–362; e-pub ahead of print (2018). [CrossRef] [PubMed]

**28. **M. Villiger, D. Lorenser, R. A. McLaughlin, B. C. Quirk, R. W. Kirk, B. E. Bouma, and D. D. Sampson, “Deep tissue volume imaging of birefringence through fibre-optic needle probes for the delineation of breast tumour,” Sci. Rep. **6**(1), 28771 (2016). [CrossRef] [PubMed]

**29. **B. Baumann, W. Choi, B. Potsaid, D. Huang, J. S. Duker, and J. G. Fujimoto, “Swept source/Fourier domain polarization sensitive optical coherence tomography with a passive polarization delay unit,” Opt. Express **20**(9), 10229–10241 (2012). [CrossRef] [PubMed]

**30. **Y. Lim, Y.-J. Hong, L. Duan, M. Yamanari, and Y. Yasuno, “Passive component based multifunctional Jones matrix swept source optical coherence tomography for Doppler and polarization imaging,” Opt. Lett. **37**(11), 1958–1960 (2012). [CrossRef] [PubMed]

**31. **B. Braaf, K. A. Vermeer, M. de Groot, K. V. Vienola, and J. F. de Boer, “Fiber-based polarization-sensitive OCT of the human retina with correction of system polarization distortions,” Biomed. Opt. Express **5**(8), 2736–2758 (2014). [CrossRef] [PubMed]

**32. **M. Villiger, E. Z. Zhang, S. K. Nadkarni, W.-Y. Oh, B. J. Vakoc, and B. E. Bouma, “Spectral binning for mitigation of polarization mode dispersion artifacts in catheter-based optical frequency domain imaging,” Opt. Express **21**(14), 16353–16369 (2013). [CrossRef] [PubMed]

**33. **N. Vansteenkiste, P. Vignolo, and A. Aspect,“Optical reversibility theorems for polarization: application to remote control of polarization,” J. Opt. Soc. Am. A **10**(10), 2240–2245 (1993). [CrossRef]

**34. **S.-Y. Lu and R. A. Chipman, “Interpretation of Mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A **13**(5), 1106–1113 (1996). [CrossRef]

**35. **B. H. Park, M. C. Pierce, B. Cense, and J. F. de Boer, “Optic axis determination accuracy for fiber-based polarization-sensitive optical coherence tomography,” Opt. Lett. **30**(19), 2587–2589 (2005). [CrossRef] [PubMed]

**36. **J. P. Gordon and H. Kogelnik, “PMD fundamentals: polarization mode dispersion in optical fibers,” Proc. Natl. Acad. Sci. U.S.A. **97**(9), 4541–4550 (2000). [CrossRef] [PubMed]

**37. **J. Li and J. F. de Boer, “Coherent signal composition and global phase determination in signal multiplexed polarization sensitive optical coherence tomography,” Opt. Express **22**(18), 21382–21392 (2014). [CrossRef] [PubMed]

**38. **Y. Furuhashi, Y. Kimura, and H. Yamane, “Higher order structural analysis of stereocomplex-type poly (lactic acid) melt-spun fibers,” J. Polym. Sci., B, Polym. Phys. **45**(2), 218–228 (2007). [CrossRef]

**39. **N. Lippok, M. Villiger, and B. E. Bouma,“Degree of polarization (uniformity) and depolarization index: unambiguous depolarization contrast for optical coherence tomography,” Opt. Lett. **40**(17), 3954–3957 (2015). [CrossRef] [PubMed]

**40. **J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld, B. Schmid, J.-Y. Tinevez, D. J. White, V. Hartenstein, K. Eliceiri, P. Tomancak, and A. Cardona, “Fiji: an open-source platform for biological-image analysis,” Nat. Methods **9**(7), 676–682 (2012). [CrossRef] [PubMed]

**41. **R. Oldenbourg, “A new view on polarization microscopy,” Nature **381**(6585), 811–812 (1996). [CrossRef] [PubMed]

**42. **R. A. Carlton, “Polarized Light Microscopy,” in *Pharmaceutical Microscopy*, R. A. Carlton, ed. (Springer, 2011), pp. 7–64. [CrossRef]

**43. **M. Wolman, “Polarized light microscopy as a tool of diagnostic pathology,” J. Histochem. Cytochem . **23**(1), 21–50 (1975). [CrossRef] [PubMed]

**44. **T. Tudor, “Vectorial Pauli algebraic approach in polarization optics. I. Device and state operators,” Optik **121**(13), 1226–1235 (2010). [CrossRef]