Spectral binning for mitigation of polarization mode dispersion artifacts in catheter-based optical frequency domain imaging

Polarization mode dispersion (PMD) has been recognized as a significant barrier to sensitive and reproducible birefringence measurements with fiber-based, polarization-sensitive optical coherence tomography systems. Here, we present a signal processing strategy that reconstructs the local retardation robustly in the presence of system PMD. The algorithm uses a spectral binning approach to limit the detrimental impact of system PMD and benefits from the final averaging of the PMD-corrected retardation vectors of the spectral bins. The algorithm was validated with numerical simulations and experimental measurements of a rubber phantom. When applied to the imaging of human cadaveric coronary arteries, the algorithm was found to yield a substantial improvement in the reconstructed birefringence maps. ©2013 Optical Society of America OCIS codes: (170.4500) Optical coherence tomography; (110.5405) Polarimetric imaging; (170.3010) Image reconstruction techniques; (170.3880) Medical and biological imaging. References and links 1. S. Yun, G. Tearney, J. de Boer, N. Iftimia, and B. Bouma, “High-speed optical frequency-domain imaging,” Opt. Express 11(22), 2953–2963 (2003). 2. M. R. Hee, D. Huang, E. A. Swanson, and J. G. Fujimoto, “Polarization-sensitive low-coherence reflectometer for birefringence characterization and ranging,” J. Opt. Soc. Am. B 9(6), 903–908 (1992). 3. J. F. de Boer and T. E. Milner, “Review of polarization sensitive optical coherence tomography and Stokes vector determination,” J. Biomed. Opt. 7(3), 359–371 (2002). 4. D. K. Kasaragod, Z. Lu, J. Jacobs, and S. J. Matcher, “Experimental validation of an extended Jones matrix calculus model to study the 3D structural orientation of the collagen fibers in articular cartilage using polarization-sensitive optical coherence tomography,” Biomed. Opt. Express 3(3), 378–387 (2012). 5. Y. Yang, A. Rupani, P. Bagnaninchi, I. Wimpenny, and A. Weightman, “Study of optical properties and proteoglycan content of tendons by polarization sensitive optical coherence tomography,” J. Biomed. Opt. 17(8), 081417 (2012). 6. K. H. Kim, M. C. Pierce, G. Maguluri, B. H. Park, S. J. Yoon, M. Lydon, R. Sheridan, and J. F. de Boer, “In vivo imaging of human burn injuries with polarization-sensitive optical coherence tomography,” J. Biomed. Opt. 17(6), 066012 (2012). 7. C. K. Hitzenberger, E. Goetzinger, M. Sticker, M. Pircher, and A. F. Fercher, “Measurement and imaging of birefringence and optic axis orientation by phase resolved polarization sensitive optical coherence tomography,” Opt. Express 9(13), 780–790 (2001). 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). 9. W. Drexler and J. G. Fujimoto, “State-of-the-art retinal optical coherence tomography,” Prog. Retin. Eye Res. 27(1), 45–88 (2008). 10. 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). #185698 $15.00 USD Received 21 Feb 2013; revised 19 Apr 2013; accepted 21 Apr 2013; published 2 Jul 2013 (C) 2013 OSA 15 July 2013 | Vol. 21, No. 14 | DOI:10.1364/OE.21.016353 | OPTICS EXPRESS 16353 11. S. K. Nadkarni, M. C. Pierce, B. H. Park, J. F. de Boer, P. Whittaker, B. E. Bouma, J. E. Bressner, E. Halpern, S. L. Houser, and G. J. Tearney, “Measurement of Collagen and Smooth Muscle Cell Content in Atherosclerotic Plaques Using Polarization-Sensitive Optical Coherence Tomography,” J. Am. Coll. Cardiol. 49(13), 1474–1481 (2007). 12. W.-C. Kuo, N.-K. Chou, C. Chou, C.-M. Lai, H.-J. Huang, S.-S. Wang, and J.-J. Shyu, “Polarization-sensitive optical coherence tomography for imaging human atherosclerosis,” Appl. Opt. 46(13), 2520–2527 (2007). 13. C. E. Saxer, J. F. de Boer, B. H. Park, Y. H. Zhao, Z. P. Chen, and J. S. Nelson, “High-speed fiber based polarization-sensitive optical coherence tomography of in vivo human skin,” Opt. Lett. 25(18), 1355–1357 (2000). 14. W. Y. Oh, S. H. Yun, B. J. Vakoc, M. Shishkov, A. E. Desjardins, B. H. Park, J. F. de Boer, G. J. Tearney, and B. E. Bouma, “High-speed polarization sensitive optical frequency domain imaging with frequency multiplexing,” Opt. Express 16(2), 1096–1103 (2008). 15. M. Yamanari, S. Makita, and Y. Yasuno, “Polarization-sensitive swept-source optical coherence tomography with continuous source polarization modulation,” Opt. Express 16(8), 5892–5906 (2008). 16. E. Z. Zhang and B. J. Vakoc, “Polarimetry noise in fiber-based optical coherence tomography instrumentation,” Opt. Express 19(18), 16830–16842 (2011). 17. M. Villiger, E. Z. Zhang, S. Nadkarni, W.-Y. Oh, B. E. Bouma, and B. J. Vakoc, “Artifacts in polarizationsensitive optical coherence tomography caused by polarization mode dispersion,” Opt. Lett. 38(6), 923–925 (2013). 18. 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). 19. S. H. Yun, G. J. Tearney, B. J. Vakoc, M. Shishkov, W. Y. Oh, A. E. Desjardins, M. J. Suter, R. C. Chan, J. A. Evans, I. K. Jang, N. S. Nishioka, J. F. de Boer, and B. E. Bouma, “Comprehensive volumetric optical microscopy in vivo,” Nat. Med. 12(12), 1429–1433 (2007). 20. S. Yun, G. Tearney, J. de Boer, and B. Bouma, “Removing the depth-degeneracy in optical frequency domain imaging with frequency shifting,” Opt. Express 12(20), 4822–4828 (2004). 21. A. Dal Forno, A. Paradisi, R. Passy, and J. Von der Weid, “Experimental and theoretical modeling of polarization-mode dispersion in single-mode fibers,” IEEE Photon. Technol. Lett. 12(3), 296–298 (2000). 22. E. Ibragimov, G. Shtengel, and S. Suh, “Statistical correlation between first and second-order PMD,” J. Lightwave Technol. 20(4), 586–590 (2002). 23. G. van Soest, M. Villiger, E. Regar, G. J. Tearney, B. E. Bouma, and A. F. W. van der Steen, “Frequency domain multiplexing for speckle reduction in optical coherence tomography,” J. Biomed. Opt. 17(7), 076018 (2012). 24. E. Götzinger, M. Pircher, W. Geitzenauer, C. Ahlers, B. Baumann, S. Michels, U. Schmidt-Erfurth, and C. K. Hitzenberger, “Retinal pigment epithelium segmentation by polarization sensitive optical coherence tomography,” Opt. Express 16(21), 16410–16422 (2008). 25. B. Park, M. Pierce, B. Cense, and J. de Boer, “Real-time multi-functional optical coherence tomography,” Opt. Express 11(7), 782–793 (2003). 26. J. Park, N. J. Kemp, H. N. Zaatari, H. G. Rylander 3rd, and T. E. Milner, “Differential geometry of normalized Stokes vector trajectories in anisotropic media,” J. Opt. Soc. Am. A 23(3), 679–690 (2006). 27. E. Z. Zhang, W.-Y. Oh, M. L. Villiger, L. Chen, B. E. Bouma, and B. J. Vakoc, “Numerical compensation of system polarization mode dispersion in polarization-sensitive optical coherence tomography,” Opt. Express 21(1), 1163–1180 (2013). 28. I. Sharf, A. Wolf, and M. B. Rubin, “Arithmetic and geometric solutions for average rigid-body rotation,” Mechanism Mach. Theory 45(9), 1239–1251 (2010). 29. M. Geissbuehler and T. Lasser, “How to Display Data by Color Schemes Compatible with Red-Green Color Perception Deficiencies,” Opt. Express 21(8), 9862–9874 (2013). 30. E. Götzinger, B. Baumann, M. Pircher, and C. K. Hitzenberger, “Polarization maintaining fiber based ultra-high resolution spectral domain polarization sensitive optical coherence tomography,” Opt. Express 17(25), 22704– 22717 (2009). 31. M. K. Al-Qaisi and T. Akkin, “Swept-source polarization-sensitive optical coherence tomography based on polarization-maintaining fiber,” Opt. Express 18(4), 3392–3403 (2010).


Introduction
Polarization sensitive (PS) optical frequency domain imaging (OFDI) is an extension of conventional OFDI [1] or optical coherence tomography (OCT) that determines the change of the polarization state of the probe light, caused by propagation and reflection within the sample [2,3].This capability is particularly useful in biological tissue with regularly arranged fibrous architecture, such as collagen, that may induce form-birefringence and result in a measurable retardation signal.PS-OCT has been used to image collagen-rich tissues such as tendon and cartilage [4,5], the loss of birefringence signal in skin burns [6], and in ophthalmology, where it enables a more detailed interpretation of the retinal structures such as the birefringent retinal nerve fiber layer [7][8][9][10].The birefringence signal measured with PS-OFDI has also been shown to correlate with collagen content of aortic plaques, which could serve as a predictor of plaque vulnerability [11,12].
Unlike free-space propagation, propagation through single mode fibers in general affects the polarization state, and, if uncontrolled, a given input state can be transformed into any possible output state.Nevertheless, in order to benefit from the advantages of fiber-based setups, Saxer et al. introduced a polarization modulation scheme [13], that was adapted in several variations [14,15].It deduces the effect of the sample on the probing light by observing the change of the polarization state of the detected light with sample depth independently for distinct consecutive (or otherwise multiplexed) input polarization states.This removes the requirement to know the precise polarization state incident on the sample and specifically allows the polarization state to change by an unknown amount as the light propagates along the fiber.However, this configuration requires that the polarization change is independent of wavelength and assumes the presence of only a negligible amount of polarization mode dispersion (PMD) [13].Yet, many current clinical applications for OCT and OFDI require the use of catheter-based systems that employ a rotary junction, are composed of a fiber interferometer with optical circulators, and require several meters of fiber for manipulation of the catheter.The combination of these elements creates a significant amount of PMD and is the limiting factor for polarization sensitivity that has hindered the clinical utility of catheter-based PS-OFDI to date [16].In a recent study, we investigated the impact of PMD on the reconstructed image of local retardation [17].The presence of PMD has been found to introduce a systemic bias to the retrieved retardation values and even induce artifacts such as apparent regional birefringent structure.At this point, it is important to recognize the statistical nature of PMD, an excellent review of which is given by Gordon and Kogelnik [18].A given optical system has one nominal value of PMD, but the specific realization of PMD depends on the precise fiber conformation and changes dramatically as a function of the catheter rotation and manipulation.The nominal value corresponds to the average PMD over all possible realizations.A single calibration prior to the measurement to subtract a PMD-induced bias is thus impossible.
In this work, we present a novel reconstruction approach that is more tolerant to instrumentation PMD and thereby mitigates its detrimental effect.The algorithm was validated and tested both with simulations and experimental measurements of a rubber phantom and cadaveric human coronary arteries.The work is organized as follows: first, some general formalism is introduced and the experimental and numerical procedures are described.Next, we present and discuss the different steps of the reconstruction algorithm.Then, the performance of the algorithm is demonstrated with numerical simulations and experimental data of a rubber phantom.Lastly, we demonstrate that this processing scheme improves local retardation images of cadaveric human coronary arteries.

PS-OFDI setup
Throughout this work, a polarization modulated, fiber-based optical frequency domain imaging (OFDI) setup was considered (Fig. 1) [1].An electro-optic modulator (EOM) transformed the source into linear and circular polarization for alternating A-lines.The light was directed towards the sample by a circulator (AC Photonics), and either coupled through a rotary junction to a side-looking optical probe for intravascular imaging [19], or scanned over the sample with a bench-top galvanometric mirror scanner.The reference light was split before the EOM so that the reference light polarization state remained constant and a polarization controller and linear analyzer were inserted prior to the receiver.An acoustooptic modulator (AOM) was applied for resolving depth degeneracy [20].The receiver stage mixed the reference and the sample signal and detected independently both orthogonal polarization states in a dual balanced configuration.The wavelength-swept laser source swept over a range of 115 nm, centered at 1320 nm, at a repetition rate of 54 kHz.The signal was digitized at 85 MHz, recording 1536 samples per sweep.Assuming a Hanning-window for the spectral shape, this resulted in a full width at half maximum (fwhm) of the intensity point spread function (PSF), i.e. the squared norm of the reconstructed tomogram of a single reflector, of f Z = 9.4 μm, in tissue, assuming a refractive index of n = 1.33.In the lateral direction, both the fiber probe and the bench-top system had a fwhm of the intensity PSF of f X = 30 μm.

Signal formalism
From the signal recorded with OFDI or another OCT modality, it is generally possible to isolate one of the two complex conjugate interference terms.As a function of the wavenumber k, and taking into account the evolution of the polarization states according to the Jones formalism, this intensity signal is proportional to the sample electric field and can be expressed as Bold capitals are matrices and bold lower case characters designate vectors.The two components of the measurement vector f for each location and wavenumber are the x-and ypolarization of the signal at the detector.J A and J B are the system Jones matrices from the source to the sample, and from the sample to the receiver, respectively.J S is the sample Jones matrix, α(k) is the source power spectrum, and e q in = [e x , e y ] T = e x e x + e y e y is the q-th input polarization state.According to the modulation scheme with processing in the Stokes domain, e 1  in = [1, 0] T and e 2 in = 2 -1/2 [1, i] T , alternating between adjacent A-lines.The wavenumber-dependence of J A and J B defines the magnitude of their PMD.System PMD refers to the combined PMD of the two elements J sys = J B J A .This is the amount of PMD to which light reflected off a mirror at the sample location is subject.The mean (across random states of J A and J B , varied by gently moving the sample fiber) PMD of the experimental system described in Fig. 1 was measured to be 85 fs.
Other than the wavenumber-dependence of the system matrices, the system is assumed free of any polarization dependent or general attenuation, i.e. all system components are modeled as unitary matrices.
The tomogram t(z,x) is reconstructed from Eq. ( 1) by Fourier transformation to obtain the spatially resolved Jones vector of the sample reflectivity.This Jones vector can be cast in the Stokes formalism, where the complex valued two-component vector is expressed as a real valued four-component vector s = [I,Q,U,V] T .The four Stokes components are, however, not independent, and for an entirely coherent signal I is equal to the Euclidean norm of [Q,U,V] T , and it is sufficient to express these three components.

Simulation of tomograms and PMD
The wavenumber dependence of J A and J B was described in terms of PMD and the sample Js was modeled as a large number of point-like scatterers, distributed on a sub-resolution scale and embedded in a birefringent medium: Here, Δn is the sample birefringence, which is alternatively cast as local retardation α, expressing the amount of retardation per path length at the central wavenumber k c .The realvalued rotation matrix P aligns the sample optic axis with the laboratory frame.For each birefringence value, 64 independent A-lines were simulated.
The system matrices J A and J B were constructed in a manner similar to [21] by multiplying 50 linear birefringent elements with random orientation.The amount of linear birefringence of each element was drawn from a normal distribution.The standard deviation of this distribution was adjusted to yield specified nominal amounts of system PMD.The simulated PMD was consistent with theoretical models for both first and second order PMD [22].J A and J B were generated independently.
For each simulated sample J S , 40 realizations of PMD matrices J A and J B with the same nominal amount of system PMD were generated.For each combination, the tomogram was reconstructed, followed by the processing to extract local retardation as discussed in the remainder of this paper, to identify the influence of independent realizations of PMD.

Algorithm
The principle strategy behind our algorithm relies on the fact that the impact of PMD depends on the ratio of its magnitude with the system resolution.An interferometer with a mean PMD of τ PMD = 50 fs results in far smaller artifacts when used with a light source providing an axial resolution of τ c = 2f Z /c = 80 fs, where f Z = 12 μm is the axial resolution and c is the speed of light, than with a system having a twice better resolution of τ c = 40 fs [17].Using a spectral window narrowed by a factor N thus allows reducing the influence of the system PMD on the polarization image.Instead of rejecting the measurements outside the window, the idea is to perform spectral multiplexing, akin to a recent report [23].The spectral interference pattern is split into multiple, overlapping, spectrally narrowed windows.From each window it is possible to retrieve a measure of local retardation.In the Stokes formalism, local retardation is represented as the 3x3 rotational matrix operator: Here, ω defines the rotation axis, and φ the amount of rotation.p and q form an orthonormal set together with ω, but are ambiguous otherwise.The rotation can thus be represented as a rotation vector conveying both the rotation axis and the amount of rotation.The rotation vectors (Ω) of the different spectral bins all represent the same tissue birefringence, and it would now be tempting to compute their arithmetic mean as an average measure of local retardation.Indeed, the spectral narrowing reduced the impact of PMD on each measure of local retardation, and an additional averaging over the bins would further reject noise contributions.However, the orientation of the rotation vectors (ω) is subject to an unknown amount of polarization rotation, induced by the system matrices J A and J B .And in the presence of PMD, these matrices induce a different amount of rotation for each spectral bin.Yet, it is possible to estimate the relative orientation between the spectral bins and correct for it.This corrective rotation is specific to each bin, but identical along an entire A-line.Even if the sample birefringence varies with depth and generates rotation vectors (Ω) that vary along z, after the corrective rotation, they are superimposed and describe the identical depth variation.With the ability to use an entire A-line to estimate the correction, enough data points are, in typical imaging experiments, available for a robust estimation.This correction generally aligns the rotation vectors (Ω) of the various spectral bins to the central bin, but maintains the vectorial nature of the local retardation (and its noise), which is beneficial for the final averaging and retrieval of the absolute amount of local retardation.
Table 1 summarizes the 6 steps that make up the proposed spectral binning algorithm.The input is the complex valued fringe pattern f q (k,x) for both input polarization states q, and the output is the local retardation.Here we describe the purpose and implementation of each of these six steps.
Step 1: Reconstruction of Jones tomogram.The first step performs the tomogram reconstruction of each spectral bin, defined by the integer ratio N of the original spectrum, i.e.Δk/N, where Δk is the spectral support of the wavelength swept laser source, spanning from k 0 to k 1 .We use a Hanning function as the window, shifted by half its width between bins.This results in 2N-1 bins, as in [23].
Step 2: Transformation to Stokes space and spatial filtering.The reconstructed Jonestomograms are cast in the Stokes formalism.The intensity-based Stokes formalism enables efficient filtering, which is necessary to reduce the impact of signal-low regions in speckle and to reduce a PMD-induced bias [17].Because of the spectral binning, the tomograms are significantly oversampled along the axial direction and axial filtering over a few pixels provides no benefit.Along the lateral direction, however, averaging over several A-lines helps to reduce speckle.A Gaussian window of fwhm w X was used for experimental data.For the numerical simulations, where the A-lines feature independent speckle realizations, a simple mean filter was applied.
Step 3: DOPU calculation and Stokes normalization.Before normalizing the filtered Stokes vectors, the ratio of the Euclidian norm of the filtered Q, U, and V components with the filtered intensity I is computed for each spectral bin, and then these ratios are averaged across the bins to obtain a measure of degree of polarization uniformity (DOPU).The DOPU is used in combination with the threshold value th dopu to generate a mask and serve as a pixelby-pixel inclusion/exclusion criteria for step 5.The measure of DOPU as introduced by Gotzinger et al. [24] takes the mean of unfiltered normalized Stokes vectors over a rectangular window with an additional intensity threshold to exclude signal-low pixels.If all the Stokes vectors are perfectly aligned, then this mean vector has a unitary norm.Instead, we average spatially filtered Stokes vectors before normalization, and compare the norm of this filtered [Q,U,V] T to the filtered I, which likewise amounts to one if all the original Stokes vectors were perfectly aligned, providing a similar measure of DOPU.In practice, both measures of Reconstruction of Jones tomogram with 2N-1 spectral bins Transformation to Stokes space and spatial filtering

t t t t t t t t t t s s s s z x m z x m h x h x w
DOPU calculation and Stokes normalization ˆ, , Extraction of local retardation matrices R(z,x,m) for each spectral bin Estimation of correction matrices P corr (x,m) DOPU result in similar masks when used with a suitable th dopu .For computational efficiency the definition in Table 1 is preferable, because the norm of [Q,U,V] T can directly be used for normalizing the Stokes vectors, and because it does not require an additional intensity threshold value.
Step 4: Extraction of local retardation matrices R(z,x,m) for each spectral bin.The local retardation matrices are calculated by comparing the Stokes vectors of each bin at depths offset by dz.The original algorithm for this extraction is based on geometric reasoning [25].It was devised to find an accurate measure of the retardation angle, and put little emphasis on the orientation and consistency of the optic axis.While this does not affect the conventional technique for extracting local retardation, this ambiguity would interfere with the vector averaging we propose in step 6.For this vector averaging, it is important that the retrieved rotation angle describes a right-handed rotation around the accompanying rotation axis.This can be implemented with a simple sign check.To further improve on the accuracy of the retrieved rotation axis, we formulated an alternative algorithm for calculating the local retardation matrix.Indeed one is trying to solve the matrix equation where u q and v q are the Stokes vectors at the -dz/2 and dz/2 offsets for the input polarization state q.With only a 3x2 measurement matrix available, the system is underdetermined.If R is assumed to represent a pure rotation, i.e. ignoring diattenuation, with only three free parameters, it becomes feasible to construct a solution for R. With a large number of additional measurements available, it is in principle possible to extract a more general R [26], but here we limit ourselves to the case of pure birefringence.Following the geometrical reasoning it is possible to find the unique rotation axis around which u 1 is mapped onto v 1 at the same time as u 2 is mapped onto v 2 .Yet, the two rotation angles differ in general, and Park et al. [25] used a relative weighting to account for the orientation of the Stokes vectors with this rotation axis and the increased susceptibility to noise for more parallel orientation.When using this original method, the final rotation angle provides a good estimation of the underlying rotation value, if compared in absolute value.Figure 2(a) (blue circles) displays the estimation accuracy for a large number of randomly generated rotation vectors, extracted from two random orthogonal normalized Stokes vectors u 1 and u 2 and their rotated versions v 1 and v 2 .Independent normally distributed noise vectors n 1,2 and m 1,2 with mean amplitude 0.1 were added to these vectors: This amount of Stokes-noise corresponds to a SNR of 10. Figure 2(b) displays the orientation accuracy of the retrieved optic axis, and Fig. 2(c) the sign consistency, when compared to the ground truth.The sign of the optic axis ω was adjusted to define a right handed system, by multiplication with the sign of (u q × v q )•ω, for the state q exhibiting the larger angle with ω, and where × is the vector cross-product and • the scalar product.To improve the rotation vector's orientation accuracy, we solve for R(z) in the least square sense rather than using the original approach.In this case, we look for the single rotation axis and the single rotation angle that minimize the overall least square error only two vectors to find the least square solution, it is possible to reformulate a geometric reasoning, that leads to the identical least square solution.Constructing two orthonormal sets from the original Stokes vector pair, and reapplying the same geometric reconstruction, it is guaranteed to find a single optic axis with identical rotation angles: ( ) ( ) In Figs.2(a)-2(c), we display the performance of this least squares solution (red squares) relative to the original approach (blue circles).In Figs.2(a) and 2(b) we can see that the accuracy of the rotation angles suffers slightly for the new formulation while the orientation of the optic axis experiences an improvement, especially in terms of the standard deviation.
Likewise, the sign-consistency with this new formulation is improved, which is crucial for the vector averaging of step 6.
Step 5: Estimation of correction matrices P corr (x,m).This step corrects for the PMDinduced variation of the tissue's apparent optic axis between the spectral bins.The rotation matrices R(z,m) of the bins are related through a similarity transformation to the ground truth R'(z): , ' As explained in [27], finding the algebraic solution to Eq. ( 7) is possible if at least two independent R(z,m) and R'(z) are available.In the present case, we have an entire A-line available to find P(m).However, the values of R(z,m) are subject to noise, and we can, at best, find an estimate of P(m).By limiting both R(z, m) and R'(z) to pure rotation matrices, the problem can be simplified to finding the rotation matrix P(m) that maps the local retardation vectors Ω(z,m) onto the retardation vectors of the central bin Ω(z,N) in the least square sense: The optimization is performed over several adjacent A-lines, defined by w x , the width of the filtering kernel h(x).Only the points x sel and z sel that have sufficient DOPU are taken into account, to mask points with random Stokes values.A threshold value of th dopu = 0.6 was used.In the special case where R'(z) defines a single optic axis along the entire tissue (or is the identity matrix in other regions), Eq. ( 7) is inherently ambiguous, as any P(m) multiplied with a rotation matrix with the same rotation axis as R'(z) results in the identical R(z,m).This is, however, of little concern, as the ambiguity is reflected only in P(m) which is not further used except to retrieve R'(z), which on the other hand is not affected itself by this ambiguity.
Step 6: Estimation of local retardation α(z,x).With the corrected rotation vectors calculated for each spectral bin, the final step computes the average rotation vector.Algebraically, this can be shown to be a first order approximation of both the Riemannian and Euclidian mean of the rotation matrix operators [28], but it is easier to compute and sufficiently accurate.The norm of the mean rotation vector, divided by dz provides the overall estimation α(z,x) of local retardation.It is the vectorial averaging of this last step that makes the spectral binning efficient.

Results
Throughout the results section, the conventional Stokes-domain processing approach [25] served as reference against which the spectral binning algorithm was benchmarked.To guarantee a fair comparison and use an equal amount of averaging in both algorithms, the conventional algorithm first computed the Stokes tomogram using the entire available spectral support, but then applied both lateral and axial spatial averaging with a kernel h(x,z).In order to be matched with the result of the spectral binning algorithm using the fraction N of the original spectrum, the shape of h(x,z) was defined as where W(k,N) is a Hanning window of width Δk/N, and h X (x) is the same lateral kernel as for the spectral binning.This way, the image of a single point-like scatterer would result in the identical intensity image for both averaging strategies.From the filtered Stokes vectors, the conventional algorithm retrieved the corresponding measure of DOPU, followed by normalization and extraction of the local retardation.Each local retardation value obtained in this conventional fashion benefited from the same amount of averaging as the corresponding value of the spectral binning algorithm.The conventional algorithm is summarized in Table 2.

Simulations
For the numerical simulations, the signal of scattering samples with uniform sample birefringence and corrupted by system PMD were computed as outlined in the section 2.3.
The sample birefringence was adjusted to yield a local retardation of 0 to 1 deg/μm, corresponding to values reported for aortic tissue [11].Noise corresponding to an average signal-to-noise ratio (SNR) of 15 dB was added.For each combination of sample (J S ) and PMD (J A , J B ) matrices, the local retardation was retrieved according to the spectral binning and the conventional reference algorithm.Lateral averaging was performed over 5 A-lines and dz was set to dz = 5.26 f Z .The number of bins was varied from 1 to 5. The mean and standard deviation of the retrieved local retardation in a region of interest spanning from 500 μm below the sample surface to 900 μm and over the 64 simulated A-lines was evaluated.For each sample birefringence 40 realizations of PMD matrices with the same nominal PMD value were generated to obtain histograms that approximate the probability density functions of the mean and the standard deviation of the local retardation.This point is crucial to account for the statistical nature of PMD.Not only the median value of these distributions is important, but also their width, as this indicates how much the local retardation is likely to change for a different realization of PMD, which is induced in the experimental setting by fiber movement.Because these distributions are frequently skewed, we report the median together with the 10% and 90% quantiles in the results.Figure 3 displays the simulation results for a sample without any intrinsic birefringence, and different nominal values of system PMD, as function of the spectral fraction N. PMD is reported as ratios of τ PMD /τ c .In the absence of sample birefringence, the influence of system PMD is most detrimental.The adverse effect of PMD is clearly visible in Fig. 3.For N = 1, both algorithms show a comparable performance.For an increasing N, however, a dramatic reduction both in the mean value as well as in the standard deviation is obvious in case of the spectral binning algorithm.The benefit of an increased N can be expressed as where f represents either the mean or the standard deviation of the local retardation.For the spectral binning, this improvement follows very accurately a 1/N 2 law, both for the mean and the standard deviation of the local retardation.Stated differently, by selecting N = 5, the impact of PMD can be reduced by a factor of 25.The reduction of the spectral width results in a degradation of the axial resolution.This can, however, be tolerated within certain limits, as will be discussed in section 5.The conventional processing likewise improves with increased N, but to a far lesser degree.
Transformation to Stokes space and spatial filtering , has width and height

t t t t t t t t t t s s s s z x z x h x z h x z w w
Extraction of DOPU and normalization of Stokes vectors Extraction of local retardation matrix R(z,x) and local retardation α(z,x) , For a sample with a finite birefringence, the detrimental impact of PMD becomes in general less severe.Figure 4 shows the distributions of the mean and the standard deviation of the local retardation for different amounts of sample birefringence, for the case N = 5, and again comparing the spectral binning with the conventional algorithm.For the lowest birefringence values, the signal is ultimately limited by SNR, and keeps a finite offset from the true value.Importantly, however, the spectral binning successfully narrows down the spread of the probability distribution for a given nominal amount of system PMD.This way, the result of the algorithm is rendered independent of the specific PMD realization and provides the identical result in a robust fashion, over the entire range of simulated sample birefringence values.The standard deviation still maintains a noticeable spread, but it is significantly reduced when compared to the corresponding conventional retrieval.Next, to find experimental demonstration of the advantages of the spectral binning algorithm, a rubber-phantom, cut from a fiber dust cap and strained by stretching, was measured with the bench-top scanner.The sample provided a homogenous backscattering signal similar to biological tissue, with a varying amount of stress induced birefringence controlled by the amount of stretching.Using the bench-top scanner maintained the identical realization of system PMD throughout the measurement, defined by the specific conformation of the fiber interferometer.

Rubber Phantom
Figure 5 displays tomograms and local retardation maps for conventional and spectral binning algorithms applied across varying induced strain.The processing parameters were N = 5, w X = 5 f X and dz = 5.12 f Z .In case of the conventional processing and higher strains, banding artifacts induced by system PMD appeared (indicated by the white arrows in Fig. 5(b)).The origin of these artifacts is described in detail in [17].The local retardation measurements obtained by spectral binning appear smoother, exhibit no banding, and find lower retardation values for the low strain configuration.For higher strain, an increase in local retardation with sample depth is found.This can be explained by the geometry of the straining: the rubber-phantom has the shape of a short tube into which two small metal rods are inserted and then pulled apart.The inner surface obtains thus the highest strain, whereas the outer surface remains slightly more relaxed, due to compression of the rubber band in the radial direction and shear forces.Part of the increase in local retardation with depth could also be attributed to the decreasing SNR.
The right hand side of Fig. 5 displays the means and standard deviations of the local retardation over three regions of interest distributed along depth and in function of the applied strain.The regions of interest spanned across the entire displayed sections, and 11 B-scans in the out of plane direction.The values of the spectral binning exhibit superior standard deviations and linearity with strain, when compared to the results of the conventional algorithm.The experimental setting suffers from a nominal PMD to coherence time ratio τ PMD /τ c of close to unity.Judging from the numerical simulations, an even more significant improvement could be expected.Experimental limitations such as polarization-dependent scattering that were not considered in the numerical simulations can account for this difference.
Despite the improvement in linearity between the strain and the local retardation, it is offset to a finite retardation in the absence of strain.It is possible that the rubber cap exhibits birefringence even in the absence of external strain due to residual internal stress.In order to confirm this hypothesis, the birefringence would have to be measured with an alternative polarimetry instrument, or the phantom could be cut from a larger block of rubber.The current analysis nevertheless clearly demonstrated an improvement in linearity and a reduction in noise obtained by spectral binning.To stress the potential of improved polarization measurements for clinical application, coronary arteries of a human cadaveric heart were measured.The protocol was approved by the institutional review board of Massachusetts General Hospital, and the time between death and imaging was less than 72h.The heart was catheterized, and pullbacks in the three coronary arteries were recorded.Figure 6 shows a typical cross-sectional view of a fibrous plaque in the left circumflex indicating, with a white arc of a circle, a region of intimal hyperplasia, and with a red arc of a circle, the region of the fibrous plaque.These structures could be clearly identified using the data from the entire pullback, but are more difficult to perceive in an individual cross-sectional intensity view.The processing parameters for computing the local retardation were N = 5, w X = 5 f X and dz = 5.12 f Z .With the experimental setting suffering from a nominal PMD to coherence time ratio τ PMD /τ c close to unity, the local retardation map obtained with conventional processing, shown in Fig. 6(b), fails to reveal any meaningful structure and demonstrates again the deleterious effect of system PMD on the polarization signatures.In contrast, in the local retardation image obtained with spectral binning (Figs.6(d) and 6(e)), a layered architecture is observed in the region of intimal hyperplasia (white arc of a circle in Fig. 6(a)), showing a first layer of lower local retardation, followed by a region of higher retardation, and then dropping back to a lower value.This matches well with the architecture of the vessel wall, which is composed of the intima as the innermost layer, followed by the media, and finally the adventitia, as indicated by the labels in Fig. 6(e).In the region of the plaque (red arc of a circle), a higher birefringence is detected in the intimal layer, consistent with the presence of fibrosis.The exact origin of the observed birefringence signal and correlation with vessel morphology merits further investigation.However, the observed polarization features agree well with the known facts that the media is rich in smooth cells (SMC) and elastin, while regions of significant fibrosis feature more abundant and thicker collagen fibers, and all these components have been shown to cause an increased birefringence signal in PS-OCT [11].Interestingly, the intensity view only shows very poor contrast for these structures, which demonstrates the complementary nature of the local retardation signal.

Intracoronary imaging
Figure 6 also displays the DOPU for both processing schemes.Using the entire spectral support, the DOPU is destroyed in many regions due to PMD.Computing the DOPU independently over the spectral bins overcomes this limitation and provides a high DOPU to about 1.5 mm inside the tissue, well beyond the media.
Because the retardation signal is complementary to the intensity tomogram, it is interesting to combine both in a single representation.To this end, we used a near isoluminescent colormap ("Ametrine" from Geissbuehler et al. [29]) to display the birefringence in color, and the intensity tomogram in the luminescence.Note the subtle difference with conventional hue-saturation-value mapping.The isoluminescent map only uses the CMYK color space and takes into account the luminance rather than the value.The combined image is displayed in Fig. 6(d).The weighting with the intensity signal in the luminescence channel efficiently masks the random retardation values in signal-low regions and provides a direct co-registration of polarization features with the tissue structure.

Discussion
Spectral binning offers a simple tool to significantly improve local retardation measurements and mitigate system PMD.Although the catheter rotation and manipulation of the catheter changes the specific realization of PMD, this no longer impacts the retrieved local retardation.The narrowing of the spectral window by factor N obviously increases the axial width of the PSF by the same factor.Yet, the local retardation is extracted by comparing the Stokes vectors at two different depths, separated by dz.This axial distance ultimately limits the axial resolution of the local retardation map.As long as N*f Z < dz, no additional sacrifice in axial resolution is made.Adjusting dz is thus critical as it also has an important impact on the accuracy of the retrieved birefringence values [17].Besides the spectral binning, averaging in the lateral direction is needed, to reduce the impact of low SNR regions in speckle pattern.We found that N = 5, in combination with a matched dz and a lateral averaging over 5 independent speckle patterns (i.e.> 5 f X ) provided satisfying results.
With spectral binning, it is possible to recover local retardation maps with a system suffering from PMD that closely match the measurements that would be possible on a system free of any PMD.As seen in Fig. 3, in the absence of PMD, lateral averaging alone would suffice and one could set N = 1.The depth spacing for differential measurements, dz, and hence the axial resolution of birefringence measurements, cannot be changed, though, without introducing an additional error on the local retardation, and spectral binning indeed allows to approach the situation of a system without any PMD.
Instead of local retardation, many prior research reports on PS-OCT have displayed the accumulated retardance, which is the retardation matrix between the Stokes vector of a reference signal -usually the tissue surface -and the depth dependent Stokes vectors.This removes the dependency on dz and works well on homogenous tissues, leaving it to the viewer to take the local derivative by eye and conclude on the birefringence of the tissue.For more involved tissue architectures with different layers of varying birefringence and changing orientation of the optic axis, such as coronary arteries, reading cumulative retardation maps becomes difficult.Extracting the local retardation provides a far more intuitive view in this case.The current work focused on local retardation and restricted the sample matrix to a pure rotational operator.In principle, a more general Muller-matrix could be assumed to extract additional polarimetry parameters, such as diattenuation [26].
An alternative solution to resolving the challenges imposed by PMD is to re-engineer the optical system and reduce the physically present amount of PMD.The specifications on differential group delay (DGD) of current fiber-circulators is about between 50 and 100 fs for each port combination, depending on the manufacturer, accounting for sufficient system PMD to pose significant difficulty in obtaining accurate birefringence measurements.Compared to the amount of PMD originating from the circulator, the contribution from the sample arm is small.With its location between the two passages through the circulator, even a small variation of its transmission matrix due to fiber motion or catheter rotation generates, however, a large variation of the overall system PMD realization.Improved DGD specifications of fiber-circulators would obviously help.Replacing the fiber circulators with a 50:50 fiber coupler would drastically reduce the amount of system PMD, but come in hand with a 5-6 dB signal penalty, which is critical for clinical imaging.
The use of polarization maintaining (PM) fibers is yet another way to circumvent the limitations of PMD in fiber-based PS systems, and has been adapted by several groups [30,31].Most current catheter-based clinical systems, however, rely on a rotary junction for helical scanning of the fiber probe, which precludes the use of PM fibers.
In a recent work [27], we introduced a calibration method that accurately determined the system matrices J A and J B individually for each A-line.With these matrices known, the detrimental effect of PMD can be entirely compensated numerically.While this method directly and in principle fully removes the effect of PMD from the system, it requires hardware modifications to generate the calibrating signals.Because spectral binning does not require changes to the imaging system hardware, it offers an interesting paradigm for reducing the impact of PMD on local birefringence.
The spectral binning algorithm works on data from conventional polarization modulated OCT data.It can be applied to already existing data sets and new data can be acquired with current, state of the art polarization modulation OFDI systems, without very restrictive requirements on the employed optical circulators and overall system PMD.This simplifies the route to investigating the potential of polarimetry measurements on clinical data.This is of particular interest for intracoronary measurements.The ability of polarization imaging to assess collagen and smooth muscle cells was previously demonstrated on ex vivo aortic plaques with free-space PS-OCT [11].With spectral binning, it now becomes possible to perform intravascular studies in animal models of atherosclerosis and in human patients.

Conclusion
In this work, we presented a processing strategy to mitigate detrimental effects of system PMD and reconstruct improved local retardation maps in a robust fashion, insensitive to fiber motion and probe rotation.By spectral binning, the impact of PMD onto each bin is reduced to a level where the usual spatial averaging is sufficient to provide a good estimate of tissue retardance.The main benefit of the proposed algorithm emerges from maintaining the matrix formalism and representing the result of each spectral bin as a depth-dependent rotation matrix.The average of the rotation matrices of the spectral bins is then approximated as the arithmetic mean of their rotation vectors.Before taking this average, however, the rotation vectors of the individual spectral bins are corrected for their PMD specific rotation error, which can be estimated from the data itself.
The algorithm was first validated with extensive numerical simulations.Experimental measurements of a rubber phantom provided further confirmation of its performance.And its application to the imaging of coronary arteries provided a very strong qualitative demonstration of its benefits.We believe that this work, in combination with ongoing and future efforts will enable catheter-based tissue polarimetry and will significantly help the advancement of OFDI in a range of clinical applications.

Fig. 2 .
Fig. 2. Estimation of (a) rotation angle, and (b) orientation of the rotation axis for the original (following [25], displayed as blue circles) and the least squares (red squares) estimation.Errorbars indicate standard deviations.(c) Frequency of correct sign estimation, i.e. the sign of the scalar product of the estimated and ground truth rotation axis.

Fig. 3 .
Fig. 3. Performance assessment of spectral binning with numerical simulation in the absence of sample birefringence.Median and 10% and 90% quantiles -indicated by error bars -of the estimated probability density functions of the mean and the standard deviation of the local retardation as function of binning N or the equivalent filtering for the conventional algorithm for different nominal amounts of system PMD.(a) Mean of the local retardation for spectral binning.(b) Mean of the local retardation for conventional processing.(c) Benefit of increased N for spectral binning (full lines) and conventional processing (dashed lines).The black thick line indicates 1/N 2 .(d)-(f) Same information for the standard deviation of the local retardation.

Fig. 4 .
Fig. 4. Performance assessment of spectral binning with numerical simulation for different amounts of sample birefringence.Median and 10% and 90% quantiles -indicated by error bars -of the estimated probability density functions of the mean ((a) and (b)) and the standard deviation ((c) and (d)) of local retardation, for spectral binning ((a) and (c)) and the conventional processing ((b) and (d)) and N = 5.The points at zero sample birefringence correspond to the rightmost values in panels (a),(b),(d),(e) of Fig. 3.

Fig. 5 .
Fig. 5. Measurement of rubber phantom with stress induced birefringence.(a) Log-scaled intensity images, averaged out of plane over 11 B-scans, for 0%, 22.2%, 44.4%, 66.6%, and 88.8% of strain.(b) Corresponding local retardation retrieved with conventional processing, with the white arrows indicating the banding artifact, and (c) spectral binning, likewise averaged over 11 B-scans out of plane.The scale bars show 500μm along depth and 1mm in the lateral direction.(d)-(f) Analysis of mean and standard deviation, indicated by the errorbars, computed over three regions of interest (over the entire displayed sections and the 11 Bscans out of plane) at increasing depths visualized in (a)-(c) by the black rectangles: (d) depth 1, (e) depth 2, and (f) depth 3. The conventional processing results in noisy measurements with oscillations that mask the linear relation between stress and local retardation found with spectral binning.

Fig. 6 .
Fig. 6.PS imaging of cadaveric human coronary artery.(a) Intensity image, indicating a region of intimal hyperplasia (white arc of a circle) and a fibrous plaque (red arc of a circle).(b) Conventional reconstruction of the local retardation and (c) the corresponding DOPU.(d) Overlay of the intensity signal as luminance and the local retardation obtained with spectral binning in isoluminescent colormap.(e) Local retardation reconstructed with spectral binning, identifying the layered architecture as intima (i), media (m) and adventitia (a).The white arrow points to the region of increased birefringence within the fibrous plaque.(f) The DOPU obtained with spectral binning.Scale bars: 1mm.

Table 1 . Spectral multiplexing algorithm for mitigation of PMD Parameters
N, h X (x), dz, th dopu Input