Fiber-based polarization-sensitive OCT for birefringence imaging of the anterior eye segment

: We demonstrate a prototype system of polarization-sensitive optical coherence tomography (PS-OCT) designed for clinical studies of the anterior eye segment imaging. The system can measure Jones matrices of the sample with depth-multiplexing of two orthogonal incident polarizations and polarization-sensitive detection. An optical clock is generated using a quadrature modulator and a logical circuit to double the clock frequency. Systematic artifacts in measured Jones matrices are theoretically analyzed and numerically compensated using signals at the surface of the sample. Local retardation images of filtering blebs after trabeculectomy show improved visualization of subconjunctival tissue, sclera, and scar tissue of the bleb wall in the anterior eye segment.


Introduction
Optical coherence tomography (OCT) is an imaging modality to measure the two or threedimensional structure of biological samples [1]. In addition to the retinal imaging that is the most popular application of OCT, an anterior segment OCT (AS-OCT) has also been developed to observe the morphometric properties of the anterior eye segment. The technique of the interferometer for AS-OCT was translated from time-domain OCT [2] into spectraldomain OCT and swept-source OCT (SS-OCT) or optical frequency domain imaging (OFDI) [3][4][5][6][7] because of advantages in the sensitivity [8][9][10] and ability of high-speed imaging [11,12].
Conventional OCT utilizes only scattering intensity to form images, and speckle contrast in a broad sense is the sole agent to distinguish different tissues [13]. Polarization-sensitive OCT (PS-OCT) has been developed to extend the function of OCT by exploiting the full nature of polarized light [14,15]. Previous studies of the anterior eye segment using PS-OCT showed promising potential in birefringence imaging for enhancing image contrast in normal tissue [16][17][18] and detecting abnormal birefringence in keratoconus corneas [19,20] and scar tissue [21,22]. In addition, it has been shown that scleral birefringence had statistically significant correlations with scleral elasticity and intraocular pressure in some in vivo and ex vivo studies [23][24][25], suggesting that scleral ultrastructure is partially related to the rheology of the eye. These reports have shown a potential of PS-OCT regarding the application of the anterior eye segment that cannot be obtained by conventional OCT.
Fiber-based PS-OCT is useful to facilitate optical alignment and to enable a flexible probe design for the measurement of Jones matrices or Stokes parameters [26][27][28]. In order to measure these parameters without reducing the effective A-scan rate in SS-OCT, frequencyor depth-multiplexing methods have been developed [17,[29][30][31][32]. Wang et al. [33] applied electric frequency doubling of the sampling clock to extend the axial measurement range of the depth-encoded PS-OCT. Although the electric frequency doubling is effective, its applicability is limited by the bandwidth of the frequency multiplier. If the frequency-swept laser has high nonlinearity of the optical frequency sweep, the multiplied clock can be distorted. In the case of frequency multiplexing with an electric 90-degree phase shifter [34], the applicability is also limited by frequency dependence of the electric phase shifter.
Another challenge in PS-OCT is wavelength-dependent polarization properties of optical components. Recent studies have suggested that polarization mode dispersion (PMD) increases polarimetric noise in PS-OCT [35,36]. Some methods have been developed to overcome this issue using various approaches [37][38][39].
In this paper, we report on our depth-multiplexed PS-OCT using a sampling clock generated by a combination of optical and electric components. We also show our methods to correct wavelength-dependence of polarization-sensitive devices in the system using signals at the surface of the sample and to reduce speckle noise of the Jones matrix using a coherent spatial Gaussian filter. Local retardation images of filtering blebs after trabeculectomy are demonstrated as a clinical application of our system.

Interferometer
The interferometer of our system employs the depth-multiplexing method with passive components [31], but the details were refined and adapted for imaging of the anterior eye segment. Figure 1(a) shows a schematic of the interferometer. Light source is a frequencyswept laser (HSL-20, Santec, Aichi, Japan), which has a center wavelength of 1297 nm, a wavelength range of 112 nm, a sweep rate of 50 kHz, a duty cycle of 51%, and an average output power of 34.0 mW. All of the optical fiber used in this system is single mode fiber (SMF). 10% of the output light is directed to a trigger/clock generator, which will be described below. 90% of the output light is directed at a Mach-Zehnder interferometer. In the sample arm, a polarization delay unit (PDU) was built to apply different optical delays for two orthogonal polarizations. The beam is divided by a polarizing beam splitter (PBS), reflected by right-angle prisms, combined at the same PBS, and coupled to SMF. Since the single PBS ensured the same tilt angle of the PBS for both the incident path and the return path, the beam direction after combining the two beams from two prisms can be easily aligned. This design reduced the effort of the alignment and the implemented size of the PDU, which was 150 × 200 mm 2 using ready-made components. After the PDU, the light passes through a circulator (Oplink Communications Inc., CA, USA), and is directed to a scanner illuminating a sample with an averaged optical power of 7.71 mW. The optical power loss between the light source and the probe was mainly due to the free space implementation of the PDU and the connections of the SMF patch cables. The scanned position on the sample is monitored by a CMOS camera through a dichroic mirror using light-emitting diodes at a wavelength of 940 nm. Backscattered light is coupled to SMF, and directed to a polarizationsensitive detection arm. In the reference arm, the light passes through a variable delay line, and is directed to the polarization-sensitive (PS-) detection unit. In the PS-detection unit, two orthogonal polarizations are separated and detected individually. Fiber Bragg grating (FBG, FBG-SMF-1264-80-0.2-A-{2}60F/E,L = 1M, Tatsuta Electric Wire & Cable Co., Ltd, Osaka, Japan) was used in front of one of two ports of each balanced receiver to embed a peak signal at a constant wavelength in the interference signal. These signals are detected by two balanced receivers that have a bandwidth of 200 MHz (C12792-2555(X), Hamamatsu Photonics, Shizuoka, Japan), filtered by a low-pass filter (SLP-300 + , Mini Circuits, NY, USA), and digitized by a high-speed digitizer (ATS9350, Alazartech, QC, Canada). The lowpass filter was selected so that the nonlinear group delay of the filter did not influence the signal phase of OCT, and was used to eliminate aliasing of high-frequency signals that passed the balanced receivers over their bandwidth. A schematic of the trigger/clock generator is shown in Fig. 1(b) and Fig. 1(c). Interference signals were generated by a Mach-Zehnder interferometer with quadrature modulation as shown in Fig. 1(b). A phase shift of 90 deg was applied between horizontal and vertical polarizations in one arm of the interferometer using a quarter waveplate (QWP). Interference signals of horizontal and vertical polarizations were divided by a PBS and detected individually by two receivers (PDB130C, Thorlabs, NJ, USA). These two signals of 125-235 MHz frequency were bandpass-filtered to suppress unused frequency bands, filtered by comparators to binarize signals at zero-crossing points, and combined with a logical XOR gate to double the signal frequency to 250-470 MHz as shown in Fig. 1(c). The XORed signal was bandpass-filtered, amplified, and supplied to the digitizer as an external clock. This electric processing using XOR gate is similar to Xi et al. except the phase shifting method [34]. One output port of the beamsplitter (BS) that was unused for clock generation in Fig.  1(b) was connected to an FBG. The transmitted light from the FBG was used for monitoring the optical power, and the reflected light from it was used as an optical trigger of the digitizer.
Our optical clock enabled an axial measurement range of 15.28 mm for standard OCT imaging. For depth-multiplexed PS-OCT imaging, the optical delay of the PDU was set to be 6.68 mm, and it was the axial measurement range for PS-OCT imaging.
The sensitivity, axial resolution and signal roll-off of the system were 97.6 dB at 3.41 mm depth, 10.1 µm in tissue and −0.54 dB/mm, respectively.

Prototype system for clinical studies
All of the interferometer except the probe of the sample arm was installed in the interferometer assembly box shown in Fig. 2. The scanning probe was installed on a scanning head, and it was connected to the sample arm of the interferometer with an SMF. All of the hardware, including a computer and a monitor, were mounted on a motorized optical table. This PS-OCT system was set up in the clinic of Tohoku University Hospital. Software for both the measurement and analysis was developed using LabVIEW (National Instruments, TX, USA) so that the system can be operated by orthoptists.

Signal pre-processing
Volumetric raw spectra were acquired using the system described in Section 2. Each spectrum has 3200 pixels that were sampled based on the optical clock. DC offset of the spectrum was calculated by 10th-order polynomial fit to each average of B-scan spectra and subtracted from each spectrum.
In principle, the spectrum has no jitter in the sampled timing owing to the opticallygenerated trigger and sampling clock. In practice, however, we observed the jitter that is in integral multiples of the sampling clock between the horizontal and vertical detection channels and among A-scans. Although the cause of this jitter is unclear, we speculate that it may be because of a stability issue in the electric analog circuit of the analog-to-digital conversion with the external clock in the digitizer. To resolve this issue, a peak signal was superposed in each spectrum at a wavelength of 1264 nm using the FBGs in the PS-detection unit. This approach is an extension of the method by Choi et al.
[40] that compensated for the timing jitter among A-scans of a single detection channel. The following processing was applied in our compensation method; 1) 128 pixels were extracted from each spectrum around the peak signal; 2) a number of the pixels were increased by 16 times using zero-padding technique; 3) the center of rising and falling edges of the peak that were detected at half maximum of the peak was defined as the peak position; 4) all spectra in the volumetric data were shifted to align the peak positions in integral multiples of the original 128-pixel resolution.
Subsequently, the spectrum was apodized by Hanning window and transformed by discrete Fourier transform (DFT). The positive-frequency range was extracted to obtain the complex OCT signal. The depth-multiplexed signals were demodulated by shifting the axial depth in a precision of the discrete pixel resolution, which was 9.55 µm/pixel. Fixed pattern noise in these complex OCT signals was removed by subtracting the complex median at each depth in each B-scan [41]. These OCT signals were transformed into spectral domain by inverse DFT for the successive processing shown in the following section.

Modeling of the electric fields with Jones formalism
In depth-multiplexed PS-OCT, four signals are detected from each A-scan using the PDU and two detectors of H and V channels. We start from an equation that was previously described by Braaf et al. in Eq. (8) of [37] to model these electric fields; is a common amplitude of the interfered signals. Since γ is already demodulated as described in Section 3.1, it is replaced with a residual phase γ ′ after the demodulation.
Equation (1) is rewritten with a new symbol Note that the terms of the depth-dependent signal decay, ( ) z β Δ and ( ) z γ β Δ + , are not compensated at this step.
is a Jones matrix of an optical path from the sample's surface up to the input fiber tip of the PS-detection unit, and ( ) i k e ψ is the phase difference between the horizontal and vertical components due to the PS-detection unit, which includes a mismatch of the optical path lengths between the horizontal and vertical components of the electric fields from the PBS to the detectors.

Compensation of the signal decay and k-dependence
In this section, we show our compensation method using the sample's surface. The depth position of the speckle peak at the sample's surface in a B-scan of the volumetric data is detected using an OCT intensity image. A narrow Gaussian window that has a standard deviation of 3 pixels (28.7 µm in air) and the center at the surface was multiplied to the complex OCT A-scan profiles that were calculated in the pre-processing step of Section 3.1, and transformed into spectral domain by inverse DFT. Assuming ( ) is an identity matrix at the sample's surface, the electric fields of the surface are derived from Eq. (3) as 1 2 Equation (4) is demodulated by multiplying a complex conjugate of the normalized and is averaged in the B-scan to increase reliability. The averaged electric fields are described as in Eq. (4) is cancelled, i e χ is a constant common phase and it does not influence the compensation process in the following. Matrix elements of Eq. (5) are individually described for later use as 1 2 Assuming that , relative values of k-dependent parameters in Eq. (5) are derived;  12 12 ( ) ( ) 1 0 : . 0 Equation (3) is rearranged using Eqs. (11) and (12) The k-dependence of ( ) sample k E is corrected using Eqs. (11)-(13) as All other matrix elements of showed k-dependence of the phase that is relative to 1 ( ) arg{ } H surf k E′ as shown in Fig. 3 Fig. 3(d). Similarly, k-and elementdependent amplitudes of ( ) surf k ′ E shown in Fig. 3(a) that were originated from k-and polarization-dependent optical losses in the PDU and the detection unit were outlined in Fig.  3(c). All of these systematic artifacts were corrected in Eq. (15).

Speckle reduction of Jones matrices with coherent spatial Gaussian filter
Speckle has both attributes as a fundamental carrier of information and noise source for OCT [13]. In PS-OCT, Jones matrix averaging methods have been proposed with global phase cancellation to reduce the speckle noise [21,37,42]. In this section, we describe our speckle reduction technique using a coherent spatial Gaussian filter.
The Jones matrix at each pixel in a B-scan has an unpredictable global phase offset because of a coherent combination of densely packed unknown scatterers. Assuming that the sample has a uniform optic axis of birefringence without diattenuation in a small rectangular spatial kernel, the Jones matrix in the kernel can be modeled instead of Eq. (16) as The axial pixel separation δz was set to be 4 pixels (27.5 µm in tissue).

Results of PS-OCT imaging
One healthy eye and three glaucomatous eyes of four human volunteers were involved in this study. The study protocol adhered to the tenets of the Declaration of Helsinki and was approved by the institutional review board of Tohoku University Graduate School of Medicine. The anterior eye segment was scanned by the PS-OCT prototype system with a raster scan protocol and a radial scan protocol. In the raster scan protocol, a lateral region of 12 × 12 mm 2 was scanned by 800 × 128 A-scans. In the radial scan protocol, a lateral range of 12 mm was scanned by 800 A-scans and the B-scan direction was rotated in a counterclockwise direction to cover 360 degrees with 128 B-scans. The best B-scan was selected from these scan patterns to show representative properties of the each subject.

Efficacy of the processing methods in a healthy human eye
The healthy anterior eye segment at the superior angle of 64-year-old female was measured by PS-OCT, and the processed results are shown in Fig. 4. The OCT intensity image in Fig.  4(b) shows the angle of the anterior chamber without artificial signal fading due to birefringent tissue [18]. The image quality of local retardation images depending on the kdependence correction and the coherent spatial Gaussian filter is shown in Fig. 4(c). A color map of Ametrine was used to display the local retardation images [44]. Without applying both processing methods, polarimetric speckle noise was too high to see any difference between subconjunctival tissue and sclera as shown in Fig. 4(c1). Although the polarimetric speckle noise was reduced by applying the coherent spatial Gaussian filter in Fig. 4(c2), the local retardation image still had noise, which decreased the contrast between the subconjunctival tissue and sclera. With the k-dependence correction but without the coherent spatial Gaussian filter, the contrast between the subconjunctival tissue and sclera was improved, yet fine speckle noise obscured fine structures in the sclera as shown in Fig. 4(c3). The speckle noise in Fig. 4(c3) was caused by the calculation method of local birefringence shown in Eq. (19), where the multiplication of two pixels could enhance the speckle noise. By applying both processing methods, the best image quality of the local retardation was obtained as shown in Fig. 4(c4). The major structures near the angle were indicated in Fig. 5(a) and Fig. 5(b) that were enlarged images of Fig. 4(b) and Fig. 4(c4) at the light-blue boxed region in Fig. 4(b). The sclera had unevenly high local retardation. The region near the scleral spur and the Schwalbe's line including the trabecular meshwork had high local retardation, which would be useful for diagnosis of the angle-closure glaucoma as suggested previously [18,45]. Iris pigment epithelium, ciliary body and choroid indicated by white arrows in Fig. 4(c4) appeared to have high local retardation because of polarization scrambling of melanin [16,17,46]. Corneal stroma that was closely perpendicular to the probing beam appeared to have high local retardation as indicated by a light-blue arrow in Fig. 4(c4), which was also reported previously [21]. The coherent spatial Gaussian filter has an additional benefit at the boundary between the tissue and void regions. Figure 6 shows magnified images at the red and green boxed regions of Fig. 4(b). In the local retardation images without the coherent spatial Gaussian filter, the local retardation was artificially high at the surface of the sample indicated by a light-blue arrow of Fig. 6(c) and at the boundary between the tissue and void space indicated by the white arrows of Fig. 6(c) and Fig. 6(d). These artifacts were caused by the calculation method of the local birefringence using separated pixels [43]. They were mitigated by the coherent spatial Gaussian filter as shown in Fig. 6(e) and Fig. 6(f). Since this filter has the effect of extrapolating the tissue signal to the void region, the resultant local retardation becomes low at the boundary.  , local retardation with k-dependence correction without coherent spatial Gaussian filter (c), (d) and local retardation with kdependence correction and coherent spatial Gaussian filter (e), (f) from the red and green boxed regions in Fig. 4(b), respectively.

Measurements of known samples
To validate the processing methods further, we measured 2% intralipid and a quarter waveplate, which have known values of birefringence. Figure 7 shows the images of the 2% intralipid with or without applying the processing methods. As Fig. 4(c3), Fig. 7(b3) had remaining speckle noise caused by the calculation method of local birefringence shown in Eq. (19). The coherent spatial Gaussian filter was effective to reduce it as shown in Fig. 7(b4). To show the statistical distribution of the images in Fig. 7, the histograms of the local retardation images were plotted in Fig. 8. It was clear that the measured local retardation became closer to zero in Fig. 8  A quarter waveplate (WPQ05M-1310, Thorlabs, NJ, USA) was measured as a birefringent sample. For this measurement, the lateral scan was not applied. In the processing, the coherent spatial Gaussian filter was disabled, because the sample does not have the speckle noise. A phase retardation from the front surface to the back surface of the sample and relative orientation of the optic axis were calculated without or with the k-dependence correction, and were plotted in Fig. 9. Each plot was calculated from a single A-scan. The orientation was left to be wrapped at 90 degrees. In both of Fig. 9(a) and Fig. 9(b), the results with the k-dependence correction were closer to the theoretical values, demonstrating the efficacy of the correction.

Functioning and nonfunctioning filtering blebs
Trabeculectomy is an ophthalmic surgery to reduce the intraocular pressure (IOP) of the glaucomatous eye by creating a filtering bleb that drains the aqueous flow. Figure 10(a), 10(c) and 10(e) are images of a functioning filtering bleb in a 78-year-old female whose IOP was 10 mmHg 10 days after trabeculectomy on the left eye for primary open angle glaucoma (POAG). The OCT intensity image in Fig. 10(c) shows a large cleft in the filtering bleb. The local retardation image in Fig. 10(e) shows negligibly low birefringence in the bleb wall. An exception is indicated by a white arrow in Fig. 10(e), which is birefringence of a suture as seen in the color photo of Fig. 10(a). The sclera and scleral flap shows birefringence as well as the sclera of the healthy eye. Figure 10(b), 10(d) and 10(f) show images of a nonfunctioning filtering bleb in a 53-year-old male whose IOP was 36 mmHg 5 years after trabeculectomy on the left eye for POAG. The OCT intensity image in Fig. 10(d) shows a narrow cleft in the filtering bleb and partial adhesion of the scleral flap to the bleb wall. The local retardation image in Fig. 10(f) shows apparent birefringence encapsulating the cleft near the inner boundary of the bleb wall as indicated by white arrows. This indicates the scar tissue that prevented the drainage of the aqueous flow from the cleft.

Bleb reconstruction surgery
Bleb reconstruction is an ophthalmic surgery to reopen the drainage path of the aqueous flow by removing the scar tissue and reconstructing the filtering bleb. A filtering bleb that had decreased functioning in the right eye of an 85-year-old female, which was diagnosed as exfoliation glaucoma and had an IOP of 22 mmHg, was measured by PS-OCT before and after the bleb reconstruction surgery. Figure 11(a), 11(c) and 11(e) show images of the filtering bleb before the surgery. The OCT intensity image in Fig. 11(c) showed narrow cleft of the filtering bleb and partial adhesion of the scleral flap to the bleb wall at the location indicated by a red arrow. In the local retardation image of Fig. 11(e), highly birefringent scar tissue was found at the location of the adhesion indicated by a white arrow. Although the cleft was narrow, it was not fully covered by the scar tissue, which would explain the moderately high IOP of 22 mmHg. Two weeks after the bleb reconstruction surgery for this patient, the IOP decreased to 4 mmHg. A color photo and PS-OCT images after the surgery are shown in Fig. 11(b), 11(d) and 11(f). In the OCT intensity and local retardation images of Fig. 11(d) and 11(f), it was observed that the scar tissue was extracted as indicated by a white arrow and the region of the cleft was enlarged. It is noted that the location of the scar tissue cannot be distinguished from that of the bleb wall without scarring if only the OCT intensity images are available. Local retardation images are invaluable to identify scar tissue objectively. Volume-rendered video clips of the local retardation before and after the bleb reconstruction surgery for the same patient of Fig. 11 are shown in Fig. 12 and Fig. 13, respectively. The movies were created using Voreen 4.4 (voreen.uni-muenster.de) [47]. In these movies, the cleft of the filtering bleb and the scleral flap were manually segmented and visualized with light-blue and red colors, respectively. Before the surgery shown in Fig. 12, the scar tissue covered a large area above the scleral flap, and the cleft region was smaller than the scleral flap. After the surgery shown in Fig. 13, a large part of the scar tissue was removed, and the cleft region was enlarged, indicating successful treatment for this filtering bleb.

Discussion
resultant local retardation images were shown in Fig. 14. Interestingly, there were no major qualitative differences among these methods in both rectangular and Gaussian window functions. Hence, these methods would work similarly for the purpose of spatial filtering except extreme cases, such as highly unbalanced SNRs of the matrix elements. Figure 14 also suggests that a deviation of the experimental data from the assumption that Eq. (17) should be unitary does not cause unstable results. We developed PS-OCT for the anterior eye segment, and demonstrated in vivo measurements of the filtering blebs after trabeculectomy. Our PS-OCT could measure highquality local retardation images of filtering blebs with the techniques of k-clock doubling, kdependence correction and coherent spatial Gaussian filter. The patient with good control of IOP showed low birefringence in the bleb wall, and the patients with poor control of IOP showed high birefringence of scar tissue in the bleb wall. These results were consistent with previous studies [21,22]. Furthermore, we demonstrated PS-OCT imaging of the filtering bleb before and after bleb reconstruction surgery. These results suggested that PS-OCT would be able to help surgeons plan needling revision or reconstruction surgery of the filtering bleb.