Robust reconstruction of local optic axis orientation with fiber-based polarization-sensitive optical coherence tomography

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


Introduction
Optical coherence tomography (OCT) measures the amount and path-length of light backscattered 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][12][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][21][22][23][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 wavelengthdependent 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.

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

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 J T sam (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(i2kz).The complex-valued Jones matrix tomogram reconstructed from J tot (k) can be described as: where we use spectral binning with a Hanning window 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, 2N−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/(2N), 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): leading to a transpose-symmetric interference signal J′ tot (k) = J cor (k)•J tot (k).With this correction, the reconstructed tomograms become also transpose-symmetric: 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: (5) subject to det(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 j † cor •H•j 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 j tom (z)•j † tom (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 transposesymmetric 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 J′ tom (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].

Compensation of system polarization distortions
Although symmetric, J′ tot (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]: where J″ tot (k) is the corrected matrix without system distortions, k c the central wavenumber and In analogy to the symmetrization of the Jones matrices, we find the compensation matrix J comp (p) for each spectral bin J′ tom (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 pth bin and the central bin across all points within an entire Bscan with sufficient signal (SNR > 20 dB): where J′ tom,scaled (z,p) is J′ tom (z,p) scaled to a determinant of +1, as explained in Appendix A, which is necessary to remove any difference in the global phase and intensity between the spectral bins.Employing the concurrent decomposition to synthesize 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 complexvalued 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.

Extraction of relative local optic axis orientation
After symmetrization and compensation for k-dependence of the system components, J″ tot (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 J″ tot (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 J″ tot (z).To suppress speckle, the SO(3) matrices are weighted with the tomogram intensity signal (the magnitude of the determinant of the Jones matrix J″ tot (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.( ) 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: where we used the fact that M T = M −1 in SO(3).At pixel locations with insufficient signal (SNR < 5 dB), 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.

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

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

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 PS-OCT 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.
In our system, the polarization-diverse optical mixer is a compact commercial product, which prevents us from inserting a linear polarizer to filter the reference signal to have identical power and phase for the orthogonal detection channels.As a result, it is the dominant source of apparent system diattenuation.However, in a system where the reference signal is polarized and correctly aligned, we suggest making the Jones matrix symmetric with a correction matrix on the right-hand side of J tot (k), canceling J in (k) instead of J out (k).This would correct more efficiently for an imbalance in the PDU, where the ratio of the power in the two orthogonal polarized illumination beams might drift over time.In our case, we carefully aligned the PDU before taking measurements, such that the two illumination beams had equal power onto the sample.As a result, J in (k) featured minimal diattenuation that could be simply ignored.To further improve the robustness to system drift and polarization dependent loss, it would be possible to estimate the residual system diattenuation present in J in (k) by minimizing the cumulative diattenuation of J″ tot (z) across an entire B-scan, under the assumption that the sample itself is free of diattenuation.In this work, we also implicitly assumed that the sample retardation is independent of the wavelength, i.e., the same phase retardation is measured across all spectral bins.The same assumption is made when reconstructing conventional PS-OCT images without spectral binning.Yet, the correction for the k-dependence of the system components relies on their effects being constant in time.kdependence of the sample retardation, if sufficiently varying in orientation across a B-scan, should not significantly impact the ability to correct for the system transmission.
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] T , 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 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.

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 A. 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: 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 n − ir n )/2 into the series: where I is the identity matrix.Owing to the anticommutator properties of the Pauli basis {σ m , σ n } = 2δ mn I, where δ mn is the Kronecker delta, the even powers of the above series simplify to: (13) Identifying the series expansions of the cosine hyperbolic and the sine hyperbolic allows rewriting Eq. ( 12) as: which leads to q 0 = 2cosh(c) and finally leads to: ( ) ( ) ( ) ( ) sinh sinh cosh 2 n n n n q q q c d ir c q 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]: ' exp arg det ' exp arg 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.

Fig. 1 .
Fig. 1.Example of symmetrization of full Jones matrix with an ex-vivo lamb trachea sample.(a) Cross-sectional intensity image of the sample.M: Mucosa; ASM: Airway smooth muscle; C: Cartilage.(b) Plots of retardation vectors r cor (p) (circles) and diattenuation vectors d cor (p) (stars) of J cor (p) and their 3 rd order polynomial fits r cor (k) (solid lines) and d cor (k) (dashed lines) across the full spectrum.k c is the central wavenumber.Normalized histogram of the circular (V) component of the retardation vectors (c) and diattenuation vectors (d) of the full Jones matrix before (red, J tom (z)) and after (blue, J′ tom (z)) symmetrization.Images of the retardation vectors of the full Jones matrix before (e) and after symmetrization (f).Images of the diattenuation vectors of the full Jones matrix before (g) and after symmetrization (h).Scale bar: 500 μm (physical length assuming refractive index of 1.40 of biological samples here and throughout this manuscript).

Fig. 2 .
Fig. 2. Compensation of system polarization distortions.(a) Retardation r comp (p) (circles) and (b) diattenuation d comp (p) (stars) components of the compensation matrix J comp (p) and their 3 rd order polynomial fits across the full spectrum, r comp (k) and d comp (k) (solid lines).(c) The average of the squared Euclidean norm of the difference between the Jones matrices of each spectral bin J′ tom,scaled (z,p) and those of the central bin J′ tom,scaled (z, k c ) before (red) and after (blue) compensation of system polarization distortions.

Fig. 3 .
Fig. 3. Generation and spatial filtering of cumulative SO(3).Image of cumulative retardance (a) and (b) diattenuation of the sample decomposed from the general Jones matrix.(c) Image of retardance of SO(3) after spatial averaging.(d) Retardance of the central A-line in the yellow box in (a) before (red) and in (b) after spatial averaging (blue).(e) Retardance of the central Aline in the green box in (a) before (red) and in (b) after spatial averaging (blue).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 haveR rt (z) = D•R T st (z)•D•R 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:

Fig. 4 .
Fig. 4. Validation with custom-made birefringence phantoms (Phantom #1 for (a-d), Phantom #2 for (e-i)).En face view of (a) OCT intensity image.Overlays on intensity image, respectively, of: (b) local phase retardation; (c) depolarization index; and (d) local optic axis orientation.(e) & (f) En face views of the local optic axis orientation of the top layer (e) and the bottom layer (f), corresponding to sections indicated by the red and blue dashed lines in (g) and (i), respectively.(g) Cross-sectional view through the green dashed line in (e) and (f).The white region in the upper left corner and the pink region in the upper right corner of (f) show the adhesive tape that fixes the phantom.(h) A frame from a 3D rendering (available online: Visualization 1) of the local optic axis orientation of Phantom #2.(i) Cross-sectional view through the pink dashed line in (e) and (f).Scale bars: 2 mm (physical length).

Figure 4 (
g) and Fig. 4(i) are cross-sectional images of the green and pink dashed lines in Figs.4(e) and Fig. (f), respectively.Figure 4(e) and Fig. 4(f) correspond to en face views of the red and blue dashed layers in Figs.4(g) and Fig. (i), respectively.

Fig. 5 .
Fig. 5. Local optic axis orientations with different system distortions.(a) The intensity of the reference signal in vertical (solid) and horizontal (dash) directions in different measurements.The optic axis orientation overlaid on the intensity images, processed, respectively, with the method in this manuscript (b) and the method in Villiger et al. [25] (c), where the OA orientations of the dashed green boxes are offset to match with their physical orientation.Plots of mean value and standard deviation of the acute (d) and obtuse (e) absolute optic axis orientation in ROIs of (b, red) and (c, blue).Scale bars: 1 mm (physical length).

Fig. 7 .
Fig. 7. En face views of lamb trachea samples measured ex vivo.OCT intensity, local phase retardation and local optic axis for: Regions 1A (a-c) and 1B (d-f) of Sample #1 and Sample #2 (g-i), corresponding to the layers indicated by the dashed yellow lines in Fig. 6, respectively.The dashed cyan lines indicate the locations of the B-scan in Fig. 6.The localized blue region in (b) and small white regions in (e) and (h) are artifacts due to strong reflections from the sample surface.M: Mucosa; ASM: Airway smooth muscle.(Video clips available online: Visualization 5, Visualization 6, and Visualization 7).Scale bar: 500 μm (physical length).

Fig. 8 .
Fig. 8. Human airway sample measured ex vivo.B-scan view of (a) OCT intensity image, (b) local phase retardation and (c) local optic axis orientation.En face views of (d, g) OCT intensity image, (e, h) local phase retardation and (f, i) local optic axis orientation of the layer approximately 180 μm and 420 μm deep into the sample from the tissue surface, respectively.(Full-depth en face view video clip available online: Visualization 8).(a-c) correspond to the dashed red lines in (d-i).(j) H&E histology image with ASM indicated by arrows.M: Mucosa; ASM: Airway smooth muscle.Small blue regions in (e) are artifacts due to strong reflections from the sample surface.Scale bar: 500 μm (physical length).

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 ] T .Similarly, real-valued d 1 , d 2 and d 3 correspond, respectively, to a diattenuation vector d = [d 1 , d 2 , d 3 ] T .σ n is the Pauli basis ( r n and d n parameters from a given Jones matrix, we use the trace property Tr(σ m •σ n ) = 2δ mn of the Pauli basis and define: