1. Introduction

Polarization-sensitive optical coherence tomography (PS-OCT) [1–3] is an extension of OCT [4] that can be used to detect the polarization properties of the imaging sample, including diattenuation [5,6], depolarization [7–9], and birefringence [3,10,11]. In biological tissue, sensitivity to birefringence is particularly relevant due to the form birefringence of numerous types of tissue, including collagen, elastin, and muscle. Due to the prevalence of these types of tissue a number of potential biomedical applications have been identified in a wide variety of organ systems, including the eye [12], skin [13], and luminal organ systems such as the cardiovascular system [14–16] and the lung [17,18]. Applications involving the eye or the skin can be tailored to using a free-space (“bulk”) imaging scheme [3,19], but many applications necessitate the use of a fiber optic catheter. The simplest and most common type of such a catheter involves a length of optical fiber with side viewing optics that is proximally connected to a rotary joint, allowing the catheter to spin independently of the rest of the system [20]. This configuration permits measurements of sample retardance, but the optic axis (OA) of the sample is distorted by the polarization effects of the system [2].

Access to the OA of the birefringent tissue could have significant applications in both differentiating and characterizing tissue, such as the in vivo assessment of airway smooth muscle in asthma [18]. In this work we describe a method for obtaining the OA of the imaging sample in endoscopic PS-OCT. This method requires no hardware modifications and imposes minimal computational load. Additionally, no pre-calibration procedure is required and as a result the technique can be applied retroactively to data already obtained.

In the following sections we outline the theory of both the problem and our proposed solution, and couple this with the results of simulations that demonstrate the impact of this theory on obtained data. We then present our experimental results, including our validation experiments in which we assess the accuracy and consistency of our calibration technique, as well as the outcome of applying our calibration to human airway data that was obtained in vivo. These experimental results illustrate both the validity and robustness of the technique itself, as well as the potential of the technique to impact future studies.

2. Measuring the OA of a birefringent sample in endoscopic PS-OCT

In the Stokes formalism [21] in which we operate for analyzing our PS-OCT data, the polarization state measured at the detectors of any implementation of a PS-OCT system is determined by the polarization effects of three distinct components on the input state: the optical path from source to sample (“system input”), the sample, and the path from sample to detector (“system output”). Mathematically, these components can be represented by a product of Mueller matrices:${\text{S}}_{\text{o}}(z)={\text{M}}_{\text{B}}{\text{M}}_{\text{SD}}(z){\text{M}}_{\text{S}}(z){\text{M}}_{\text{A}}{\text{S}}_{\text{i}}$, where ${\text{S}}_{\text{o}}$ is the measured state, ${\text{S}}_{\text{i}}$ the input state,${\text{M}}_{\text{B}}$ the system output matrix, and${\text{M}}_{\text{A}}$ the system input matrix. The sample matrix is a function of the image depth z and has been further divided into the forward component${\text{M}}_{\text{S}}(z)$ and the reflected component${\text{M}}_{\text{SD}}(z)$, where the D subscript here and throughout the manuscript indicates the relationship${\text{M}}_{\text{D}}=\text{D}{\text{M}}^T\text{D}$, with $\text{D}=\text{diag}(1,1,1,-1)$ [9]. From this point onward we make the following two assumptions: 1) that the system elements are pure retarders, and 2) that the sample can generally be treated as a series (in depth) of linear retarders. The former can be made a more realistic approximation by compensating for polarization mode dispersion (PMD) with the use of spectral binning [22], while the latter assumption is consistent with observation for many types of tissue [23].

Following these assumptions, all of the Mueller matrices involved are reduced to rotation matrices, and a reference measurement (typically the reflection from the sample’s surface) can then be used to obtain a similarity transformation of the sample matrix from which the sample retardance can be readily obtained [2]. The OA of the sample, however, is rotated by the unknown M_B. This is further complicated when a rotating fiber optic imaging catheter is employed, causing a component of the system to be dependent on the angle of the catheter: ${\text{S}}_{\text{o}}(z,\theta )={\text{M}}_{\text{B}}{\text{M}}_{\text{PD}}(\theta ){\text{M}}_{\text{S}}{(z)}^2{\text{M}}_{\text{P}}(\theta ){\text{M}}_{\text{A}}{\text{S}}_{\text{i}}$, where$S_0(z,\theta )S_0{(z_{s0},\theta )}^{-1}={\text{M}}_{\text{B}}{\text{M}}_{\text{PD}}(\theta ){\text{M}}_{\text{Q}}{\text{M}}_{\text{S}}{(z)}^2{\text{M}}_{\text{Q}}^T{\text{M}}_{\text{PD}}^{}{(\theta )}^T{\text{M}}_{\text{B}}^T$ is the Mueller matrix of the catheter as a function of its angular position$\theta$. We note that here we are treating the sample as being uniform both azimuthally and in depth. The former we do only for simplicity and clarity in describing the theory. Sample inhomogeneity in depth, however, can introduce additional transformation factors that also need to be compensated for [24].

We now make a third assumption that$S_0(z,\theta )={\text{M}}_{\text{B}}{\text{M}}_{\text{PD}}(\theta ){\text{M}}_{\text{Q}}{\text{M}}_{\text{S}}{(z)}^2{\text{M}}_{\text{Q}}{\text{M}}_{\text{P}}(\theta ){\text{M}}_{\text{A}}{S}_i$ can be decomposed such that${\text{M}}_{\text{P}}(\theta )={\text{M}}_{\text{P}}\text{R}(\theta )$, where$\text{R}(\theta )$ represents the fiber rotation; i.e., that the fiber birefringence is constant along the rotation. This formulation imposes an active transformation of the sample along with the fiber optic catheter and in doing so maintains a consistent coordinate system along the circumference of the sample. This coordinate system is depicted in Fig. 1. The final modification that we make to our equation for the measured Stokes vector is to separate out the fixed component of the catheter sheath from the Mueller matrix of the sample:

 

Fig. 1 Depiction of the coordinate system at the sample surface, in Stokes space. This system results in rotational invariance for rotations in the plane perpendicular to the axis of the catheter.

Download Full Size | PDF

From Eq. (2) it is apparent that in order to solve for the OA of the sample at each angular position of the catheter it is first necessary to determine the product of the three rotation matrices$M_B$,${\text{M}}_{\text{P}}(\theta )$, and${\text{M}}_{\text{Q}}$. Specifically, we wish to determine this product and apply it to Eq. (2) in order to resolve the physical OA in the local coordinate system of the sample, as depicted in Fig. 1. We note that while this coordinate system is not absolute in that it depends on the position of our catheter (in contrast with that of the “laboratory” coordinate system), it nevertheless would allow us to obtain the information we desire; that is, the intrinsic characteristic of the sample itself. We have performed simulations that demonstrate the compounding effects of these individual components on sample OA measurements, which we summarize in Fig. 2. This data was generated using randomized values for the retardance and OA of the individual system components included in Eq. (1), with the “measured” OA values obtained from the similarity transformation of the respective matrices. For example, the “measured” sheath OA [Fig. 2(d)] was calculated according to the product ${\text{M}}_{\text{B}}{\text{M}}_{\text{PD}}(\theta ){\text{M}}_{\text{Q}}^2{\text{M}}_{\text{PD}}^T(\theta ){\text{M}}_{\text{B}}^T$ . These results highlight several key features of the effect of the preceding components on the ones that follow them: the axis of${\text{M}}_{\text{PD}}(\theta )$ is rotated from the V = 0 plane by$M_B$, the axis of${\text{M}}_{\text{Q}}$ rotates in a circle parallel to${\text{M}}_{\text{PD}}(\theta )$ and offset in part by the magnitude of the fibers birefringence, and the planar (V = 0) sample value OAs become non-planar and have the appearance of being randomly scattered on the Poincare sphere.

 

Fig. 2 Simulation data demonstrating the effect of the system on sample OA measurements. Each data point represents an angular position of the catheter fiber, with the range encompassing 0 to 2π. (a) Actual and (b) measured catheter OA. (c) Actual and (d) measured sheath OA. (e) actual and (f) measured sample OA. The randomly generated sample OA values are shifted from the Q-U plane to seemingly arbitrary locations on the full Poincare sphere.

Download Full Size | PDF

3. Description of technique for obtaining the sample OA

As mentioned in Section 2, solving for the OA of the sample requires solving the product of the rotation matrices$M_B$,${\text{M}}_{\text{PD}}(\theta )$, and${\text{M}}_{\text{Q}}$ in order to recover the coordinate system described in Fig. 1. We have previously demonstrated that we can measure${\text{M}}_{\text{Q}}$ using the inner and outer surfaces of the sheath in the image, which we know is oriented along the –Q direction in the coordinate system of the sample [18]. Any attempt to determine${\text{M}}_{\text{PD}}(\theta )$ introduces two significant complications: 1) a priori knowledge of${\text{M}}_{\text{PD}}(\theta )$ is limited and cannot be calibrated or fixed due to catheter dynamics, and 2) there is no reference point within the image range from which to obtain the similarity transformation of${\text{M}}_{\text{PD}}(\theta )$. We can partially overcome these difficulties by making the following observation:

 

Fig. 3 Experimental results for catheter and sheath OA measurements. (a) Catheter and (b) sheath measurements over one rotation. (c) Catheter and (d) sheath measurements from a 3 cm pullback performed in vivo, with the OA values shifted and scaled to the range of [0,1] and depicted in true RGB (where the Q component is mapped to the R channel, and the U and V components the G and B channels, respectively). These images demonstrate the consistency of these measurements even during in vivo imaging.

Download Full Size | PDF

The measurements we are now able to make provide us with equations in the following forms:

In practice, polarimetry noise prevents all the measured OA values of${\text{M}}_{\text{Q,meas}}(\theta )$from aligning perfectly to a single axis. To compensate for this, the final rotation${\text{M}}_{\text{V}}$is determined from which rotation aligns the mean vector to Q = −1. Alternatively, additional filtering and/or interpolation may be performed in order to yield more consistent results in cases where the signal from the ball lens or the sheath is not well-defined. Once applied, this method offers a significant advantage over our previously published work in that it allows us to recover unambiguously the orientation of all linear OA values of our sample, and not just those limited to the Q-axis as we previously reported [18].

4. Experimental results

The experimental data presented here was obtained using a swept-source PS-OCT system with an electro-optic modulator for input polarization diversity. This is a standard system configuration which we have previously described [18]. Imaging was conducted using a number of different catheters with either ball lens or GRIN lens distal optics. No differences were observed in the performance of our calibration technique using the different catheters. To mitigate the effect of polarimetry noise the Stokes data that was used for calibrating the sample OA was filtered along the azimuthal direction using a Gaussian filter with an angular width of 4.2° at 1 standard deviation. The sample data was acquired from Stokes data filtered as above and using a depth offset of 20 µm for the differential birefringence calculations [6]. All data processing was performed using Matlab 2015a.

In order to first determine the accuracy of our calibration algorithm, we performed imaging using a birefringent rubber phantom. The rubber phantom, when stretched, exhibits a high degree of OA uniformity along its stretched axis [18]. It is therefore a useful alternative to imaging a waveplate, and is a closer analog to the imaging of biological tissue. We imaged the phantom using two different approaches in order to investigate the accuracy and consistency of our calibrations: first, the phantom was imaged at a fixed position with respect to the catheter, and rotated about its axis from 0 to 180° in 10° increments [Fig. 4(a)]. In the second approach, the phantom was maintained at a fixed orientation (with its stretched axis at 36° with respect to the long axis of the catheter) and rotated around the catheter during imaging over one complete rotation [Fig. 4(b)]. The results of these two experiments are depicted in Fig. 5 and Fig. 6, respectively. The latter experiment allowed us to ascertain the consistency with which a given orientation is measured at each angular position of the catheter. In both experiments the PS-OCT system and catheter were operated using the same parameters we apply to our in vivo imaging (with the probe rotating approximately 33 times per second) and with no special handling or arrangement that would preclude translating this technique to an in vivo setting.

 

Fig. 4 Experimental configurations for rubber phantom imaging. (a) Configuration in which the phantom is rotated about its axis and the angle $\theta$ is measured. (b) Configuration in which the orientation of the phantom in the imaging plane is fixed as indicated by ${\theta }_f$, and maintained at this angle as the entire phantom is rotated around the catheter ($\phi$).

Download Full Size | PDF

 

Fig. 5 Experimental results obtained by rotating a phantom from 0 to 180° in 10° increments. (a) Measured versus actual orientation and (b) measured V component. The mean absolute error for the orientation was 1.4° and the mean V component was −0.02.

Download Full Size | PDF

 

Fig. 6 Rotation of the phantom (at 36° fixed orientation) 360° around the catheter. (a) The individual, uncorrected Stokes parameters. (b) The uncorrected angle data combined with the structural image. (c) The Stokes parameters and (d) angle data post-correction. The mean measured angle was 36.7° (S.D 2.64°) and the mean V component was −0.1 (S.D 0.118). Scale bar, 1 mm.

Download Full Size | PDF

The combination of these two experiments allows us to quantify both the accuracy and consistency of our calibrations. The results of these data demonstrate that the V component of the sample OA is minimized following the calibrations, as expected. The angular measurements for the OA orientation were then obtained in the Q-U plane. A small region within the phantom encompassing approximately 3 speckles was averaged to produce each data point, with no out of frame averaging performed. From these experiments we determined that with our calibration technique we can measure the OA orientation both accurately (mean absolute error: 1.4°) and consistently over the full rotation of the catheter (S.D: 2.64°).

To verify that our method can be applied in data acquired in a clinical setting, we retroactively applied the technique to data that was acquired of human airways in vivo during bronchoscopy. These data sets were obtained over a period of years and with a variety of catheters. The airway wall is both useful for demonstrating the OA calibration and an important biological application of the technique. The elasticity and the need for caliber regulation of the airway wall means it is comprised of a number of birefringent tissues, with varying orientations. This includes the longitudinally oriented fiber bundle network in the submucosa [25] and the circumferentially-oriented airway smooth muscle, which wraps around the airway wall at angles up to 20° [26].

Figure 7 demonstrates the results of applying the technique to an airway. The OA image [Fig. 7(a,b)] is depicted in HSV format with the OA data encoded in the hue, the structural image the value, and constant saturation. By isolating a +/− 20° range within the corrected OA image we can easily identify the circumferentially-oriented airway smooth muscle [Fig. 7(c)]. In contrast, the retardance image alone [Fig. 7(d)] highlights the increase in information that is obtained by analyzing the OA data: the retardance of the airway wall layers is both weak and homogeneous, and as a result it is difficult to obtain any new information by analyzing the retardance alone. However, by masking the retardance image with the ASM image in Fig. 7(c), we can further isolate the retardance associated with the ASM band [Fig. 7(e)]. This data can then be further analyzed to assess muscle tone, which is relevant to certain diseases and disorders of the airway [18].

 

Fig. 7 OA correction technique applied to a cross-sectional image of an airway acquired in vivo from a healthy human volunteer. (a) Uncorrected and (b) corrected OA images combined with the structural image. The range of corrected values can be restricted to identify the airway smooth muscle (white arrow) (b). Using this result as a mask we can extract from the retardance image (c) the retardance associated with the muscle alone (d). Scale bars, 500 µm.

Download Full Size | PDF

5. Conclusion

In conclusion, we have developed a calibration technique for obtaining the sample OA in endoscopic PS-OCT data sets. This technique is computationally efficient, robust, and does not require any pre-calibration procedures or hardware modifications. Using a uniform retarder we have demonstrated that this technique provides both accuracy and consistency in measuring the OA direction along the full rotation of the endoscopic catheter. We have retroactively applied this calibration technique to in vivo endoscopic imaging data sets and have verified that the calibration approach is viable in the clinical setting. We anticipate that this technique may be of significant use in differentiating and characterizing tissue across multiple organ systems and applications.