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

Qingyun Li; Karol Karnowski; Peter B. Noble; Alvenia Cairncross; Alan James; Martin Villiger; David D. Sampson

doi:10.1364/BOE.9.005437

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_{tot} (k) = J_{out} (k) \cdot [\int_{0}^{\infty} J_{sam}^{T} (z) \cdot J_{sam} (z) \cdot \exp (i 2 k z) d z] \cdot J_{in} (k) .

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_{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(i2kz). The complex-valued Jones matrix tomogram reconstructed from J_tot(k) can be described as:

J_{tom} (z, p) = F T [J_{tot} (k) \cdot w (k, p)] \approx J_{out} (p) \cdot J_{sam}^{T} (z) \cdot J_{sam} (z) \cdot J_{in} (p),

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₀+ pΔk/(2N), where k₀ 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):

J_{cor} (k) = J_{in}^{T} (k) \cdot J_{out}^{- 1} (k),

leading to a transpose-symmetric interference signal

J_{tot}^{'} (k) = J_{cor} (k) \cdot J_{tot} (k)

. With this correction, the reconstructed tomograms become also transpose-symmetric:

J_{tom}^{'} (z, p) = J_{cor} (p) \cdot J_{tom} (z, p) \approx J_{in}^{T} (p) \cdot J_{sam}^{T} (z) \cdot J_{sam} (z) \cdot J_{in} (p) .

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:

\min_{J_{cor} (p)} \sum_{B-Scan} {‖ J_{tom}^{T} (z, p) \cdot J_{cor}^{T} (p) - J_{cor} (p) \cdot J_{tom} (z, p) ‖}_{2}^{2},

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}^{†} \cdot H \cdot 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) \cdot 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 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 $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].

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

Download Full Size | PDF

3.2 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]:

\begin{matrix} J_{t o t}^{' ​'} (k) = J_{comp}^{T} (k) \cdot J_{tot}^{'} (k) \cdot J_{comp} (k) \\ = J_{in}^{T} (k_{c}) \cdot [\int_{0}^{\infty} J_{sam}^{T} (z) \cdot J_{sam} (z) \cdot \exp (i 2 k z) d z] \cdot J_{in} (k_{c}), \end{matrix}

where

J_{tot}^{' ​'} (k)

is the corrected matrix without system distortions, k_c the central wavenumber and J_comp(k) = J_in⁻¹(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

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 B-scan with sufficient signal (SNR > 20 dB):

\min_{J_{comp} (p)} \sum_{B-Scan} {‖ J_{comp}^{T} (p) \cdot J_{tom,scaled}^{'} (z, p) \cdot J_{comp} (p) - J_{tom,scaled}^{'} (z, k_{c}) ‖}_{2}^{2},

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

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.

Download Full Size | PDF

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

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 A-line in the green box in (a) before (red) and in (b) after spatial averaging (blue).

Download Full Size | PDF

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 $R_{rt} (z) = D \cdot R_{st}^{T} (z) \cdot D \cdot 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:

R_{st} (n \cdot d z) = Q_{n} \cdot \dots \cdot Q_{2} \cdot Q_{1} \cdot Q_{0} .

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:

\begin{array}{l} Q_{0} = \sqrt{R_{rt} (0)} \\ Q_{1} = \sqrt{D \cdot Q_{0} \cdot D \cdot R_{rt} (d z) \cdot Q_{0}^{T}} \\ Q_{2} = \sqrt{D \cdot Q_{1} \cdot Q_{0} \cdot D \cdot R_{rt} (2 \cdot d z) \cdot Q_{0}^{T} \cdot Q_{1}^{T}} \\ ⋮ \\ Q_{n} = \sqrt{D \cdot R_{rt} [(n - 1) \cdot d z] \cdot D \cdot R_{rt} (n \cdot d z) \cdot R_{st}^{T} [(n - 1) \cdot d z]}, \end{array}

where we used the fact that M^T = M⁻¹ 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(kc) is minimally diattenuating, its retardation component still impacts R_rt(z) and defines R_rt(z = 0). Taking its matrix root to compute Q₀ 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.

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

Download Full Size | PDF

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.

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

Download Full Size | PDF

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

Fig. 6 Cross-sectional B-scan views of lamb tracheal samples measured ex vivo. (a), (d), (g) OCT intensity images. Overlay of (b), (e), (h) local phase retardation on intensity images. Overlay of (c), (f), (i) local optic axis orientation on intensity images. (a-c) Region 1A and (df) Region 1B from Sample #1 and (g-i) Sample #2. (Video clips available online: Visualization 2, Visualization 3, and Visualization 4). M: Mucosa; ASM: Airway smooth muscle; C: Cartilage. Scale bar: 500 μm (physical length).

Download Full Size | PDF

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

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

Download Full Size | PDF

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

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

Download Full Size | PDF

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

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:

J = \exp (- i \frac{1}{2} \sum_{n = 1}^{3} r_{n} \cdot σ_{n} + \frac{1}{2} \sum_{n = 1}^{3} d_{n} \cdot σ_{n}),

where n = 1, 2, 3. Real-valued r₁, r₂ and r₃ correspond, respectively, to the Q, U, V components of a retardation vector r= [r₁, r₂, r₃]^T. Similarly, real-valued d₁, d₂ and d₃ correspond, respectively, to a diattenuation vector d = [d₁, d₂, d₃]^T. σ_n is the Pauli basis (n=0, 1, 2, 3):

σ_{0} = [\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}], σ_{1} = [\begin{matrix} 1 & 0 \\ 0 & - 1 \end{matrix}], σ_{2} = [\begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix}], σ_{3} = [\begin{matrix} 0 & - i \\ i & 0 \end{matrix}] .

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:

\begin{matrix} J = I + (\sum_{n = 1}^{3} f_{n} \cdot σ_{n}) + \frac{1}{2!} {(\sum_{n = 1}^{3} f_{n} \cdot σ_{n})}^{2} + \frac{1}{3!} {(\sum_{n = 1}^{3} f_{n} \cdot σ_{n})}^{3} + \dots \\ = \sum_{k = 0}^{\infty} \frac{1}{k!} {(\sum_{n = 1}^{3} f_{n} \cdot σ_{n})}^{k}, \end{matrix}

where I is the identity matrix. Owing to the anticommutator properties of the Pauli basis {σ_m, σ_n} = 2δ_mnI, where δ_mn is the Kronecker delta, the even powers of the above series simplify to:

{(\sum_{n = 1}^{3} f_{n} \cdot σ_{n})}^{2 k} = \underset{c^{2}}{\underset{︸}{{(\sum_{n = 1}^{3} f_{n}^{2})}^{k}}} \cdot I = c^{2 k} \cdot I, where k = 0, 1, 2 \dots

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

\begin{matrix} J = \sum_{k = 0}^{\infty} \frac{{(\sum_{n = 1}^{3} f_{n} \cdot σ_{n})}^{2 k}}{(2 k)!} + \sum_{k = 0}^{\infty} \frac{{(\sum_{n = 1}^{3} f_{n} \cdot σ_{n})}^{2 k + 1}}{(2 k + 1)!} \\ = \sum_{k = 0}^{\infty} \frac{c^{2 k}}{(2 k)!} I + \sum_{k = 0}^{\infty} \frac{c^{2 k + 1}}{(2 k + 1)!} (\sum_{n = 1}^{\infty} \frac{f_{n}}{c} \cdot σ_{n}) \\ = \cosh (c) I + \sinh (c) (\sum_{n = 1}^{3} \frac{f_{n}}{c} \cdot σ_{n}) . \end{matrix}

To retrieve the 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:

q_{n} = T r (σ_{n} \cdot J),

which leads to q₀ = 2cosh(c) and finally leads to:

d_{n} - i r_{n} = \frac{q_{n} c}{\sinh (c)} = \frac{q_{n} \cosh^{- 1} (q_{0} / 2)}{\sinh (\cosh^{- 1} (q_{0} / 2))} .

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₁₁ of J′ is zero, similar to the method of Li et al. [17]:

q_{n} = T r (σ_{n} \cdot \frac{J^{'} \cdot \exp (- i \arg (J_{11}))}{+ \sqrt{\det (J^{'} \cdot \exp (- i \arg (J_{11})))}}) .

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]

Name	Description
Visualization 1	3D rendering of the local optic axis orientation of a birefringence phantom with two layers
Visualization 2	B-scan views of intensity, local phase retardation and local optic axis orientation of the lamb tracheal sample (Region 1A of Sample #1) measured ex vivo.
Visualization 3	B-scan views of intensity, local phase retardation and local optic axis orientation of the lamb tracheal sample (Region 1B of Sample #1) measured ex vivo.
Visualization 4	B-scan views of intensity, local phase retardation and local optic axis orientation of the lamb tracheal sample (Sample #2) measured ex vivo.
Visualization 5	En face views of intensity, local phase retardation and local optic axis orientation of the lamb tracheal sample (Region 1A of Sample #1) measured ex vivo.
Visualization 6	En face views of intensity, local phase retardation and local optic axis orientation of the lamb tracheal sample (Region 1B of Sample #1) measured ex vivo.
Visualization 7	En face views of intensity, local phase retardation and local optic axis orientation of the lamb tracheal sample (Sample #2) measured ex vivo.
Visualization 8	En face views of intensity, local phase retardation and local optic axis orientation of a human airway sample measured ex vivo, with surface flattening in post processing.

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

Abstract

1. Introduction

2. Methods

3. Algorithm

3.1 Symmetrization of the full Jones matrix

3.2 Compensation of system polarization distortions

3.3 Extraction of relative local optic axis orientation

4. Results

4.1 Validation with birefringence phantoms

4.2 Tissue imaging

5. Discussion

6. Conclusions

Appendix

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

Funding

Acknowledgments

Disclosures

References

Supplementary Material (8)

Cited By

Figures (8)

Equations (17)

Biomedical Optics Express