Fiber-based polarization-sensitive OCT of the human retina with correction of system polarization distortions

In polarization-sensitive optical coherence tomography (PS-OCT) the use of single-mode fibers causes unpredictable polarization distortions which can result in increased noise levels and erroneous changes in calculated polarization parameters. In the current paper this problem is addressed by a new Jones matrix analysis method that measures and corrects system polarization distortions as a function of wavenumber by spectral analysis of the sample surface polarization state and deeper located birefringent tissue structures. This method was implemented on a passive-component depth-multiplexed swept-source PS-OCT system at 1040 nm which was theoretically modeled using Jones matrix calculus. High-resolution B-scan images are presented of the double-pass phase retardation, diattenuation, and relative optic axis orientation to show the benefits of the new analysis method for in vivo imaging of the human retina. The correction of system polarization distortions yielded reduced phase retardation noise, and better estimates of the diattenuation and the relative optic axis orientation in weakly birefringent tissues. The clinical potential of the system is shown by en face visualization of the phase retardation and optic axis orientation of the retinal nerve fiber layer in a healthy volunteer and a glaucoma patient with nerve fiber loss.


Introduction
In the past two decades optical coherence tomography (OCT) has become a standard clinical tool in ophthalmology for non-invasive cross-sectional imaging of retinal and choroidal structures [1,2]. OCT is based on optical low-coherence interferometry and provides high axial (2-15 µm) and lateral (5-30 µm) resolutions which makes it ideal for in vivo visualization of single tissue layer morphology and the diagnosis of pathology on a microscopic level [2]. The polarization state of the light changes while passing through several ocular structures such as the cornea [3], the retinal nerve fiber layer (RNFL) [4,5], Henle's fiber layer [6,7], the retinal pigment epithelium (RPE) [7][8][9], the lamina cribrosa [10], and the sclera [11]. Polarization-sensitive OCT (PS-OCT) uses these polarization state changes to measure the depth-resolved polarization properties of individual tissue structures [12][13][14]. In the eye, PS-OCT has proven valuable for the identification of aforementioned tissue structures, the quantitative analysis of RNFL and Henle's fiber layer phase retardation [5,6] In order to facilitate optical alignment and allow the use of fiber-based sample interfaces (e.g. endoscopes) PS-OCT systems are often built with single-mode optical fibers (SMFs) [21]. In SMFs stress-induced birefringence and imperfections in the fiber core symmetry pose a problem by changing the polarization state of the OCT light to an unknown state, which can become aligned in parallel or orthogonal to the sample optical axis [21,22]. Polarization modulation and multiplexing methods have therefore been developed for fiber-based PS-OCT that probe the sample with several distinct (orthogonal) polarization states to solve this problem [21][22][23][24][25][26][27][28]. It was recognized early on that these methods only could achieve optimal results when the changes in the polarization states are wavelength independent, i.e., if polarization mode dispersion (PMD) in the PS-OCT system is negligible [21]. In SMFs PMD causes a differential group delay between the two principal polarization states, which results in axial displacements between the detected polarization states in PS-OCT [21,29,30]. It was shown recently by Villiger et al.
[30] that PMD can add significant amounts of noise and even cause erroneous local retardation changes in PS-OCT measurements. It is therefore evident that PMD and other wavelength-dependent effects should be minimized or corrected. The amount by which the SMF length can be reduced in fiber-based PS-OCT systems is limited and therefore the development of numerical solutions for PMD correction is required. In PS-OCT microscopy on muscle/tendon tissue Zhang et al. [31] showed that a set of three sample arm calibration signals can be used to measure and numerically compensate PMD effects. In a separate study, Villiger et al.
[32] presented a robust Stokes-analysis method that worked directly on endoscopic PS-OCT sample data, and mitigated PMD effects at the expense of a lower axial resolution using a spectral binning method.
In this paper we present a method to correct wavelength-dependent polarization distortions introduced by the setup in ophthalmic PS-OCT measurements of the retina. A passive-component depth-multiplexed swept-source PS-OCT setup at 1-µm was developed and is presented together with its theoretical Jones matrix modeling. The correction method extends the conventional PS-OCT Jones matrix analysis method of Park et al.
[33] with a spectral binning method to measure the polarization states of the retinal surface and deeper located birefringent tissue structures as a function of wavelength. Together with the theoretical model of the setup the analysis of these two tissue structures allowed the determination and correction of the system polarization distortions. This included the simultaneous correction of PMD effects and wavelength-dependent power differences between the polarization states in the interferometer reference and sample arms. Highresolution images are presented of the polarization properties in the macula and optic nerve head areas of a healthy subject and a glaucoma patient.

Experimental system description and Jones matrix modeling
The developed PS-OCT setup was based on the passive-component depth-multiplexed PS-OCT methods as published by Lim et al. [27] and Baumann et al. [28]. The optical layout of the setup and an example of the acquired OCT images are shown in Fig. 1, and are respectively discussed in sections 2.1 and 2.2. This is followed by the Jones matrix modeling of the PS-OCT setup in section 2.3.

Interferometer and ophthalmic interface hardware
The setup used a 1-µm wavelength swept-source (Axsun Technologies, Inc., MA, USA) with a fixed repetition rate of 100 kHz and an average output power of 23 mW. Light from the swept-source was coupled into a single-mode fiber-based interferometer (left panel Fig. 1) where it was split into a sample and a reference arm in a 70/30 ratio by a fiber coupler.
In the sample arm the light was coupled into a polarization-delay unit (PDU, marked by a dashed box in Fig. 1) where two orthogonal polarization states were created and sent simultaneously to the eye. This was achieved by splitting the incoming light in its two polarization components with a polarizing beam splitter (PBS), after which the two components each travelled a different path length through air before they were recombined by another PBS. This created an optical path delay for the horizontal linear polarized state with respect to the vertical linear polarized state and consequently multiplexed them in depth in the OCT image. Each state was additionally purified by a polarizer (P) to ensure orthogonality of the polarization states, and to suppress inter-state interference as described by Ju et al. [34]. The optical power was set equal for both states using a polarization controller (PC) on the input fiber of the PDU. The light launched from the PDU was divided by a 20/80 coupler in two arms going respectively to the ophthalmic interface for imaging and a calibration mirror (C) for the generation of a reference signal. Returning light from both arms was analyzed in a polarization-sensitive detection unit (described below).
In the reference arm 10% of the light was split off by a fiber coupler to feed a fiber-based Mach-Zehnder Interferometer (MZI). The MZI used two 50/50 couplers with unbalanced airpaths and a balanced detector (1837 GHz-Nirvana, New Focus, CA, USA) to generate a stable interference signal that was used for k-clocking of the data acquisition. The remaining 90% of the reference arm light was transmitted through open air to match its path length and chromatic dispersion to the sample arm.
The polarization-sensitive detection unit consisted of a 50/50 fiber coupler where the light from the sample and reference arms was recombined, and PBS cubes to split the light from both output leads into its orthogonal polarization components. The horizontal and vertical polarization components were recorded by separate balanced detectors (PDB430C, Thorlabs Inc., NJ, USA) denoted as H-det. and V-det. in Fig. 1 respectively. Polarization controllers were installed to align the polarization states of the coupler's output leads, and to distribute the reference arm power equally over both detectors. This alignment was optimized to yield equal optical power and spectral shape for all detector photo-diodes as measured by their fast monitor outputs.
The ophthalmic interface (upper right panel Fig. 1) directed the OCT light (red paths) via a set of galvanometer scanners (GS) and a telescope to the human eye to image the posterior segment. The bundle incident on the cornea had an optical power of 1.5 mW and a diameter (1/e 2 ) of 1.1 mm to provide a theoretical diffraction-limited spot-size (1/e 2 ) of 20 µm on the retina. In order to improve the subject's fixation stability and facilitate instrument alignment a fixation target LED grid (FT) and a pupil camera (CAM) (blue paths) were coupled together via a glass plate (GP) and combined with the OCT light using a dichroic mirror (DM). The DM was placed under 45° in between the scan lens and the ophthalmic lens to ensure a constant incidence angle on the DM during scanning of the OCT beam with negligible (angle dependent) polarization effects. All single-mode fiber components were used at standard lengths (1 -2.5 m lead length) and together had a physical length of ~14 m per interferometer arm. In the sample arm the input light was converted in two depth-multiplexed orthogonal polarization states by a polarization delay unit (PDU) which extended the optical path of the horizontal linearly polarized state. Both states were launched simultaneously into the eye via an ophthalmic interface to record two depth-multiplexed OCT images with a dual-channel polarization-sensitive detection unit (H-det. & V-det.). In the reference arm 10% of the light was fed to a fiber-based Mach-Zehnder interferometer (K-clock) to create a stable interference signal that was used to clock the data-acquisition hardware linearly in wavenumbers. Component abbreviations, PBS: polarizing beam splitter, P: polarizer, PC: polarization controller, C: calibration mirror, GS: galvanometer scanners, DM: dichroic mirror, GP: glass plate, FT: fixation target, CAM: pupil camera.

Data-acquisition and signal pre-processing
After detection, the signals from H-det. and V-det. were low-pass filtered (SBLP-200+, 200 MHz −3 dB cutoff frequency, Mini-circuits, NY, USA) to reduce high-frequency noise and recorded by a 14-bit resolution dual-channel data-acquisition (DAQ) board (X5-400M, Innovative Integration, CA, USA). The k-clock signal from the MZI was high-pass filtered (SHP-20+ , 15.5 MHz −3 dB cutoff frequency, Mini-circuits) to remove low-frequency signal offsets, and provided the DAQ sampling clock. The MZI air-gap mismatch was set to 8.4 mm to create a k-clock that supported an OCT imaging depth range of 4.2 mm. The sweep-trigger of the swept-source initiated the digitization of A-scans linearly in wavenumbers with 1536 samples for an optical bandwidth of 92 nm centered at 1040 nm. A system sensitivity of 92.2 dB was measured from a mirror near the zero-delay depth as the incoherent intensity sum of the four signals. The signal roll-off with depth was −6 dB over 4.5 mm, and the axial resolution was 8.2 µm in air.
In the lower right panel of Fig. 1 an example is given of the acquired B-scan images for both detectors. By extending the optical path of the horizontal linearly polarized state in the PDU with half the OCT image depth range (2.1 mm), two depth-multiplexed images were recorded per detection channel. In what follows, the corresponding OCT signals are denoted as H 1 and H 2 for H-det., and V 1 and V 2 for V-det. for the horizontal and vertical states respectively. The calibration mirror created two sets of reference signals (C) that were used to phase-stabilize successive A-scans by cross-correlation. This procedure was implemented as described earlier by our group [35] without the need for numerical k-space resampling in this case. Afterwards fixed-pattern noise was subtracted from the data [36] and optical chromatic dispersion was corrected digitally [37] before polarization-sensitive analysis was performed as described in sections 3 and 4.

Jones matrix modeling of the passive-component depth-multiplexed PS-OCT setup
In OCT the signal intensity I of a single scatterer that is measured at the detector is described as a function of wavenumber k by the interference between the complex electric field vectors of the reference | 2 describe the signal intensity from respectively the reference and sample arm, 2Re{ } describes the interference between both arms, and Re{} denotes the real part. In practice the non-interfering terms are removed by balanced detection in swept-source OCT and are therefore omitted from further analysis. The complex electric field can be obtained from Eq. (1) through digital Hilbert transformation [39]: are defined at the detector, and are described by the electric field emitted by the swept-source , and subsequent transformation of its polarization state by the components of each arm as represented by their Jones matrices J ref (k) and J sam (k). While propagating to the detector both fields travel over their own optical distances z ref and z sam , and are attenuated by scalar factors c ref (k) and c sam (k) that are related to optical loss and sample reflectivity: With polarization-sensitive detection the horizontal E H (k) and vertical E V (k) polarization components of ( ) where Δz = z sam -z ref gives the path-length difference between the interferometer arms, c t (k) = 2·c ref (k)·c sam (k), and R H (k) and R V (k) are the horizontal and vertical components of . The sample arm electric field can be further decomposed into a chain of Jones matrices describing individual optical path sections [33,41]: horizontal and linearly vertical polarized. In the ideal case both states have identical optical powers and are perfectly linearly polarized such that can be modeled by an identity matrix and a multiplicative constant [28]. In practice this is not the case and is more accurately modeled by a diagonal matrix that describes the differences between both states, and a common complex electric field component e in (k): In Eq. (7), α(k) describes wavelength-dependent amplitude differences between the two states, β(Δz) is the depth-dependent signal decay (roll-off) from the swept-source instantaneous linewidth, γ describes the depth-displacement introduced by the PDU to the horizontal linear state, and ρ(k) is the phase-difference between the horizontal and vertical field components of the (possibly elliptical) polarization state of ( ) . Together Eqs. (5)-(7) model the electric fields of the four detected signals (H 1 , H 2 , V 1 , and V 2 ): where the common scalar c t (k)·e ikΔz of Eq. (5) is included in e in (k) for simplicity.

Jones matrix analysis to correct system polarization distortions
The principles for tissue polarization analysis using Jones matrix modeling in fiber-optic PS-OCT systems were described by Park et al. [33]. The core element of this analysis is the Jones matrix modeling of the polarization state at the sample surface, which acts as a reference state to assess polarization changes at deeper locations. In this section the analysis of Park et al.
[33] was used as the basis for a more elaborate analysis in which also wavelength-and depth-dependent setup distortions are addressed.

Signal roll-off correction and depth-delay demodulation
The Jones matrix analysis started with the selection of a single B-scan from an in vivo retinal data set for which the data was pre-processed as described in section 2.2. Subsequently the complex electric field information was obtained via the digital Hilbert transform step of Eq.
(2) with the addition of a division by the amplitude roll-off curve in the depth domain to correct for β(Δz) and β(Δz + γ) in Eq. (8). The amplitude roll-off curve was obtained in a separate experiment by measuring the signal amplitude on a mirror that was displaced over the full OCT imaging depth. The intensity roll-off difference between the two depthmultiplexed states was maximally −3.6 dB, which is equivalent to an amplitude factor of 0.66. In addition, the depth-displacement γ between the two launched polarization states was measured from the calibration mirror signals by determining the first order phase difference (slope) as a function of wavenumber. The depth-delayed state was demodulated (frequencyshifted) on top of the undelayed state by applying an opposing first order phase-ramp. This typically corrected γ up to a quarter of a pixel, since PMD elsewhere in the setup caused slight under-or over-correction.
Rewriting Eq. (8) for these corrections, the sample data is described by E sample (k):

K-dependent sample surface state determination
In the next step the sample surface state was measured as a function of wavenumber using a spectral binning method, for which the principle is shown in Fig. 2(A). In spectral binning the OCT fringe data as a function of wavenumber is sequentially multiplied with a series of shifted narrow Gaussian apodization windows to select different spectral bands. The spectral bands are individually processed to obtain spectral properties from the OCT signal such as variations in amplitude and phase. Spectral binning has been successfully demonstrated in a variety of OCT signal processing methods such as spectroscopic OCT analysis [42], the reduction of speckle noise [43], the measurement of optical chromatic dispersion [44,45], and PMD-robust Stokes-based processing of PS-OCT data [32].
The PS-OCT data was first Fourier transformed without spectral binning and the location of the retinal surface pixels were determined by image segmentation. The PS-OCT data was then separately processed with spectral binning for 15 overlapping spectral bands of 384 samples (~23 nm) in width that were equally spaced over the OCT spectrum with interval steps of 82 samples (~5 nm). For every spectral bin the data was Fourier transformed and the pixels at the retinal surface were extracted. The surface state was then created by complex averaging over all surface pixels for each of the four electric field components (E H1 , E H2 , E V1 , and E V2 ). In the averaging process first every pixel was normalized to (divided by) E H1 to account for the amplitude and phase offset of each pixel (see also section 4.4) [46]. In Fig.  2(B) the surface state components E H1 , E H2 , E V1 , and E V2 are plotted respectively for their real and imaginary parts at the sample number of each spectral bin center. It can be clearly seen that E H2 , E V1 , and E V2 varied gradually but with significant amplitude across the spectrum while E H1 was constant as a consequence of the normalization step before the complex averaging. The smooth variation of the electric field components enabled the use of spline extra-and interpolation to extend the surface state to all wavenumbers of the OCT spectrum (solid lines in Fig. 2(B)). Since the sample surface pixels were not (or minimally) affected by the polarization properties of the retinal tissues, the surface state E surf (k) was modeled according to Eq. (8) with exclusion of J s (k) and only contained terms describing the setup components and the cornea: Equations (9) and (10) were used to calculate the similarity transformation of J s (k): In Eq. (11) J out (k) and the reference arm diagonal matrix are cancelled out, which removed the optical path sections with the longest SMF lengths (together ~21 m) and therefore eliminated the largest contributors of PMD. The preservation of J in (k) was further preferred over more common methods that preserve J out (k) [33], since J in (k) is non-diattenuating and can be described as a unitary matrix of the SU(2) group. This was not the case in the current setup for J out (k) and the reference arm diagonal matrix, since polarization dependent losses within the bulk optics of the polarization-sensitive detection circuit caused deviations from SU (2) [47] and complicated their exact modeling.

Eigenvector analysis of deeper located birefringent tissue structures
In this section the remaining system polarization distortions are determined by simultaneously characterizing the unknown setup parameters J in (k) and α(k)·e i(kγ'+ρ(k)) . This is only possible for Eq. (11) when J s (k) has distinct eigenvalues and eigenvectors (i.e. J s (k) is not defective) which prevents Eq. (11) from cancelling itself to an identity matrix. In practice this condition is satisfied when the measurement is subject to birefringence or diattenuation. The analysis described in this section was therefore performed only on selected tissue pixels with significant phase retardation (>0.1 rad) and signal-to-noise ratio (SNR >10 dB). A large quantity of pixels in a single retinal B-scan is normally suited, since the birefringent RNFL and Henle's fiber layer impose phase retardation on all deeper tissue layers. These pixels were obtained and analyzed for different spectral bands by applying the spectral binning procedure of section 3.2 to Eq. (11). The eigendecomposition of Eq. (11) gives the left eigenvectors: where R s represents a unitary change of basis for J s (k), and s 1 and s 2 are two complex scalars resulting from the eigendecomposition algorithm (LAPACK, eig() function, Matlab 2010b, The Mathworks, Inc.) to normalize the length of the eigenvectors to 1. R s is modeled as a standard rotation matrix and describes the tissue optic axis as shown later in section 4.1. In order to use ν(k) for the correction of the remaining system polarization distortions, s 1 and s 2 have to be removed since they are artificial and do not represent any physical parts of the system. This is achieved by decomposing ν(k) in its individual terms and reconstructing the most common realization of ν(k) without s 1 and s 2 . J in (k) describes the sample arm optical path from the exit of the PDU through the setup and the cornea, and can be separated in two Jones matrices: J in (k) = J c (k)·J sys (k). In this decomposition J c (k) describes the local corneal polarization properties and J sys (k) describes the setup components between the PDU and the eye. The matrices J c (k), J sys (k) and R s are members of the SU(2) group and for simplicity J in −1 (k)·R s −1 = J sys −1 (k)·J c −1 (k)·R s −1 can therefore be replaced in Eq. (12) by a wavelength-dependent SU(2) matrix in the form [x(k), -y*(k); y(k), x*(k)] with |x(k)| 2 + |y(k)| 2 = 1. Rewriting Eq. (12) according to these steps gives: Subsequently Eq. (13) was used to determine α(k) and the amplitudes of s 1 , s 2 , x(k) and y(k): It was found that the parameters |ν 11 |, |ν 12 |, |ν 21 |, |ν 22 |, and φ ν (k) were described by Gaussian-like distributions that varied around their mean due to (speckle) noise and local variations in the corneal polarization properties J c (k) and tissue optic axis orientation R s . In Fig. 3(A) examples of histograms of all five parameters are shown for the center spectral bin. Clear distinct distributions were observed from which the position of the histogram maximum (the mode) represented the most common value for each parameter over all pixels. The modes of |ν 11 |, |ν 12 |, |ν 21 |, |ν 22 |, and φ ν (k) for all spectral bins were used for the calculations of Eqs. across wavenumber for these parameters is shown in Fig. 3(B). The upper graph shows the amplitudes of |x(k)| and |y(k)| with a maximum variation of ~0.1 which was related to wavelength-dependent polarization rotations of J sys (k). The power difference factor α(k) was found to deviate from 1 due to imbalance in the optical power distribution by the PDU between the two polarization states. In the lower graph of Fig. 3(B) the blue curve shows the variation of φ ν (k) across wavenumber which consists of a small linear ramp caused by the residual demodulation factor γ' and a sinusoidal-like component that was related to the fiber orientation of J sys (k). The parameters were extended to all wavenumbers (solid lines).
The obtained values for α(k), |x(k)|, |y(k)|, and φ ν (k) were used to construct a matrix Q(k) that represents an estimate of the most common realization of ν(k), and could be rewritten according to the matrix terms of Eq. (13) without s 1 and s 2 : where φ xy (k) = arg(x(k)·y(k)) and the overbars denote the most common realization of a variable or matrix. In general the polarization properties of the setup components were constant for all A-scans in a single data set and therefore the product between Q −1 (k) and ν(k) cancelled them out: Equation (19) is close to an identity matrix, but the variation of J c −1 and R s −1 around their most common realizations caused considerable amplitude in the off-diagonal elements for the majority of the realizations of κ(k), and still allowed the determination of φ xy (k) from each column vector of κ(k) via: In Fig. 3(A) the lower right histogram shows the distribution of φ xy (k) for the center spectral bin, which had a clear distinct peak although wider and lower in peak-height than was observed for φ ν (k). The lower graph in Fig. 3(B) shows the variation across wavenumber of φ xy (k) in red, which was small since the larger part of the system distortions was already corrected.

Defining system correction factors
The analysis derived in sections 3.2-3.3 allowed the determination of system correction factors from a single B-scan of a volumetric data set that subsequently could be applied to all other B-scans. Individual correction factors were obtained for setup distortions of both the optical path into the eye, as well as for the optical path out of the eye. The former is obtained from Eqs. (18) and (20) and is described as: and describes the most common realization of ν(k) multiplied with a phase factor. The latter is obtained out of the surface state description of Eq. (10) and Q in (k) as: Using Q in (k) and Q out (k) the original sample data E sample (k) as given by Eq. (9) of every Bscan was corrected for all system distortions: where the depth-resolved Jones matrices are written as bold italic variables, J rc and s R are constant with depth, and the convolution with the swept-source spectral shape is omitted for simplicity. In Eq. (24) e in (z) is a complex scalar that describes a common amplitude and phase which are respectively related to the reflectivity and the random axial position of the scatterers within each image pixel.

Intensity, phase retardation and diattenuation
The intensity of every pixel was calculated from Eq. (24) as the sum of the intensities for the individual field components:

Spatial electric field averaging for SNR improvement and speckle reduction
Speckle noise and low SNR can significantly degrade the quality of measured tissue polarization parameters [46,53]. A spatial averaging method for complex electric fields was therefore developed to simultaneously improve the SNR and reduce speckle noise, at the expense of a small reduction in resolution. This method uses similar basic concepts as described for the Jones matrix averaging methods of Lim et al. [46] and Ju et al. [34].
Each pixel is subject to an individual offset in amplitude and phase as given by e in (z) in Eq. (24). The phase offset is seemingly random for every pixel and spatial complex averaging between pixels is therefore only effective when these phase offsets are corrected. The simplest way to do so is to divide the electric field components of each pixel through their individual Ê H1 (or any other) field component, which removes the offsets in amplitude and phase and aligns Ê H1 for all pixels. This simple method is however not optimal since the alignment on a single field component is subject to (speckle) noise. The method is therefore extended by including the other field components and neighboring pixels in the alignment process. A spatial kernel was defined including n pixels. For every pixel the phase angle between the reference electric field component Ê H1 and the other components Ê V1 , Ê H2 , and Ê V2 were calculated. Since these angles depended on the local polarization properties of the tissue, they had a strong spatial consistency, and were therefore averaged over all n pixels in the kernel with a weighting by their amplitudes to suppress (speckle) noise influences: The phase angles of Eq. (31) were used to correct Ê V1 , Ê H2 , and Ê V2 to the phase angle of Ê H1 , and subsequently the phase angle of the vector sum over all field components was calculated: which gave the estimated noise-free phase angle of Ê H1 . The angle ψ cor was used to correct the four original electric field components of the center pixel of the kernel. This procedure was applied to all pixels in a B-scan to obtain a uniform alignment of the phase offset. Afterwards, spatial complex averaging was applied per electric field component using the same spatial kernel size to improve the SNR and reduce the speckle noise for every pixel.
In practice, the current method used one of the diagonal electric field components Ê H1 or Ê V2 as the reference field component, since in the nearly similar matrix form of Eq. (24) these have the largest amplitudes and give therefore the best performance. The method is however easily adapted to use other reference components in case the amplitudes of Ê H1 or Ê V2 are insufficient.

System performance and imaging results
The in vivo measurements in humans adhered to the tenets of the Declaration of Helsinki and were approved by the Erasmus Medical Center's Institutional Review Board. Informed consent was obtained from each subject.

Technical system evaluation
In order to evaluate the system polarization distortions a healthy volunteer was measured and the polarization state at the retinal surface was determined as a function of wavenumber as described by Eq. (10). The sample surface state effectively describes the polarization transfer function of the system including the cornea, and can be used to estimate the depth shift between the input polarization states due to PMD in the system [30]. Ten data sets were acquired with different orientations of the sample arm fiber to determine the average PMD depth shift of the system. An average PMD depth shift of 1.6 ± 0.44 µm (0.29 ± 0.08 pixel) was measured which is equivalent to a time delay of 11 ± 2.9 fs. This value is in agreement with earlier measurements of Villiger et al. [30] in a swept-source OCT system at 1310 nm without fiber-optic circulators.
The compensation of the system polarization distortions by the correction steps described in sections 3.2-3.4 is illustrated in Fig. 4. In row (A) the retinal surface polarization state was plotted as a function of wavenumber in normalized Stokes vectors on the Poincaré sphere. The first (blue), the middle (red), and the last (green) B-scans from a typical volumetric data set of the retina were selected, for which the vertical linearly polarized input state was plotted. The B-scans were acquired 3.3 s apart from each other over a total volume acquisition time of 6.6 s. The left Poincaré sphere of Fig. 4(A) shows the Stokes vector traces for the case without correction, which traveled over an arc of approximately a quarter of the globe from the -U to the +V pole. The nearly identical traces indicate that the system polarization distortions were constant over the acquisition time of the data set. Small shifts in the global location of the traces are attributed to changes in corneal polarization properties. The right Poincaré sphere of Fig. 4(A) shows the same Stokes vector traces after applying the correction for system polarization distortions. The system correction factors were calculated from the middle B-scan and subsequently applied to the other B-scans. It can be clearly seen that after correction the Stokes vector traces reduced to single locations on the Poincaré sphere, which indicates the complete correction of the system polarization distortions. In addition, the locations near the -Q pole were consistent with vertical linearly polarized light, while the offsets from this pole are attributed to different corneal birefringence for these Bscans. The horizontal linear input state showed the same behavior at the exact other side of the Poincaré sphere (not shown). These results confirmed that it is possible to apply the system correction factors obtained from a single B-scan to the whole data set, while maintaining optimal correction. The same analysis was performed for the polarization state of birefringent RNFL tissue that was obtained slightly inferior from the optic nerve head at a depth of 159 µm. The analysis was performed on a cross-sectional tissue area of 276 µm × 36 µm (width × depth) with an average DPPR of 0.53 rad. In row (B) of Fig. 4  function of wavenumber without correction for system polarization distortions. In this case the polarization state of the RNFL tissue follows a trace similar as was observed for the uncorrected polarization state at the retinal surface. The right Poincaré sphere shows the same Stokes vector trace after correction for the system polarization distortions. It can be clearly seen that after correction the Stokes vector trace reduced to a small location on the Poincaré sphere, which indicates a successful correction of the system polarization distortions. These results confirm the applicability of the developed correction method to in vivo birefringent tissues. In addition a quarter-wave plate (10RP34-1064, Newport, CA, USA) was measured at different orientation angles to evaluate the accuracy in determining DPPR and optic axis orientation [25,33]. The measured DPPR with and without correction for system polarization distortions was respectively 3.05 ± 0.07 rad and 3.05 ± 0.06 rad, for a manufacturer specification of 3.08 rad at 1040 nm. The difference between the set and measured axis orientation with and without correction was respectively 1.8 ± 1.6 deg and 1.5 ± 1.2 deg. The equivalent results with and without correction indicated that the system polarization distortions have minimal impact on the measurement of strong birefringent samples, and accurate values were obtained for the phase retardation and optic axis orientation in both cases.

Cross-sectional PS-OCT imaging in a healthy volunteer
In Fig. 5 the effect of correcting system polarization distortions is shown in cross-sectional PS-OCT images of the macula in a healthy volunteer. The images were obtained slightly inferior to the foveal center, consist of 2000 A-scans, and are displayed over 1.3 mm in depth (in tissue) and 6.6 mm (22°) in width. Spatial electric field averaging was performed over 3 × 7 pixels (depth × width) to reduce (speckle) noise and improve the SNR. Pixels outside the tissue were detected by an intensity threshold (SNR < 8 dB) and are displayed in gray in the images of the polarization parameters. In Fig. 5(A) the intensity information visualizes the structures of the fovea (right side) and a part of the inferior nerve fiber bundle (left side). Images of the DPPR without and with correction for the system polarization distortions are shown respectively in Figs. 5(B) and 5(C). In both images a clear increase in DPPR was observed for the inferior nerve fiber bundle and the sclera. Especially the DPPR of the sclera increased rapidly with depth, and saturated the image dynamic range. The most dominant difference was the more uniform (less noisy) appearance of the DPPR after correction. Also the increased DPPR caused by Henle's fiber layer was shown more clearly after correction, which was observed in the inner/outer segment junctional layer (IS/OS) flanking the center of the fovea. These differences were even better appreciated when the dynamic range of the images was scaled back from π/2 to π/4 in Figs. 5(D) and 5(E) to visualize small changes in DPPR without and with correction. In these images it is evident that the noise on the DPPR was lower after correction, which was confirmed by the average measured DPPR of 0.20 ± 0.19 rad without and 0.14 ± 0.11 rad with correction for a 10 pixel thick tissue layer starting from the retinal surface. In addition the average DPPR of Henle's fiber layer was measured as 0.22 ± 0.08 rad without and 0.28 ± 0.06 rad with correction. The latter is in close agreement with results reported by Cense et al. [6] who used a PS-OCT system with minimal SMF length and consequently should have been minimally affected by PMD.
The diattenuation is shown in Figs. 5(F) and 5(G) without and with correction for the system polarization distortions, respectively. In these images it can be seen that the correction minimized the diattenuation to nearly zero and significantly reduced its noise. This is shown as a darker and more uniform representation of the diattenuation in Fig. 5(G) and is especially appreciated in the strong reflecting layers of the RNFL and IS/OS. The median and interquartile ranges of both these layers combined were 0.09 (0.04 -0.17) without and 0.06 (0.03 -0.12) with correction; and were chosen here over the mean and the standard deviation due to the severe skewness of the diattenuation distribution. These results are consistent with previous studies which reported a negligible diattenuation for biological tissues [33,54,55].
Figures 5(H) and 5(I) show the relative optic axis orientation without and with correction for the system polarization distortions, respectively. The relative optic axis orientation was noisy for pixels with a low SNR or a low phase retardation since the eigenvectors are poorly defined in these cases. This occurs especially at the retinal surface where the surface state correction removed all polarizing effects, but was also evident around the fovea where polarizing tissue structures were absent. In the image without correction there is however a false indication of a uniform optic axis orientation around the fovea, which was due to the data deviating from the SU(2) form in the analysis of section 4.3. The optic axis orientation was however well defined in structures with strong DPPR such as the RNFL inferior nerve bundle and the sclera, and shows equivalent results with and without correction. This is in contrast with the axis orientation cast by Henle's fiber layer in the IS/OS layer, which only revealed the expected radial distribution after correction and is seen here in cross-section.
In Fig. 6 images are shown for the superior optic nerve bundle with correction for system polarization distortions. The images consist of 2000 A-scans and are displayed over 0.9 mm  Fig. 6(A) clearly demonstrates the negligible diattenuation for this tissue in Fig. 6(B). In addition, the nerve fiber bundles created strong DPPR and are shown in clusters separated by blood vessels in Fig. 6(C). In Fig. 6(D) the relative optic axis orientation visualizes the direction of the fiber bundles in leaving the optic nerve head area. The nearly twofold color cycle in the lateral direction indicates a full rotation of the optic axes in the Poincaré sphere. This is consistent with the nerve fibers leaving the optic nerve head in directions distributed over a 180° angle in real space and is in agreement with earlier PS-OCT findings of Pircher et al. [7].

En face retinal PS-OCT imaging in a healthy volunteer and a glaucoma patient
The performance of the developed PS-OCT system for clinical measurements is demonstrated by imaging the macula and optic nerve head areas from a healthy volunteer and a glaucoma patient. For each area a volumetric PS-OCT data set was obtained with 300 B-scans and 2000 A-scans / B-scan in an acquisition time of 6.6 s. Compiled en face PS-OCT maps were created for the intensity, the DPPR, and the relative optic axis orientation. The intensity en face map was calculated by summing the logarithmic scaled intensity over depth including both the retina and choroid. The en face maps of the polarization parameters were obtained from the IS/OS since this layer displays the polarization effects for the complete retina with high SNR [7,17]. In addition, the RNFL thickness of both volunteers was measured with Scanning Laser Polarimetry (SLP) using a GDxPRO (Carl Zeiss Meditec Inc., CA, USA). This instrument scans the retina with polarized light at 785 nm and converts the detected DPPR to a measure of RNFL thickness using a fixed factor of 0.67 nm/µm (single pass phase retardation to thickness) [56]. The SLP thickness maps are plotted in their native clinically used color map, and cropped and positioned according to the PS-OCT field of view.
The results for the right eye of the healthy volunteer are shown in Fig. 7 for a field of view of 5.3 × 9.6 mm 2 (17.8° × 32.1°). In Fig. 7(A) the en face intensity image displays retinal structures such as several blood vessels leaving the optic nerve head. In Fig. 7(B) the en face map of the retinal DPPR shows the strong retardation of the nerve fiber bundles leaving the optic nerve head in the superior and inferior directions. Around the fovea a (donut) ring of increased retardation can be seen caused by Henle's fiber layer. The optic axis orientation is shown in Fig. 7(C), which clearly visualizes the orientation of the nerve fibers leaving the optic nerve head in all radial directions, while the radial extension of the Henle fibers is clearly observed in the fovea. Note that in PS-OCT the measured optic axis orientation angles in the Stokes convention are subject to a pi-ambiguity in the Poincaré sphere [51]. This leads to a wrapping of the optic axis orientation every 90° in real space, and consequently gives a fourfold wrapping for a full 360° rotation as is shown for the fovea and optic nerve head [6]. The RNFL thickness map as measured by SLP is given in Fig. 7(D) and shows an equivalent retardation pattern as observed with PS-OCT including the retardation ring from Henle's fiber layer and the nerve bundles superior and inferior of the optic nerve head.
The results for the left eye of the glaucoma patient are shown in Fig. 8 for a field of view of 5.8 × 8.9 mm 2 (19.3° × 29.7°). Similar to the healthy subject, Fig. 8(A) shows the en face intensity image with the retinal structures. Figure 8(B) shows the en face DPPR map with significantly lower retardation in the superior nerve bundle and an almost absent retardation in the inferior bundle in comparison with the healthy volunteer. In contrast the retardation ring of Henle's fiber layer was well preserved. The optic axis orientation in Fig. 8(C) shows the severe loss in the inferior nerve bundle as an absence in orientation which is uniformly colored in cyan in the lower left image corner. The SLP thickness map shows a similar pattern with a significantly decreased thickness for the inferior RNFL bundle as can be seen in Fig. 8(D). This is in agreement with the PS-OCT findings, although seemingly more retardation was observed by SLP than with PS-OCT. This difference might be attributed to an imperfect corneal compensation for the SLP [19]. The retardation ring of Henle's fiber layer was however observed similarly to the PS-OCT.