Optical coherence tomography (OCT) is a noninvasive imaging technique for recording high-resolution cross-sectional images of transparent and translucent tissues.1, 2, 3 Conventional OCT measures spatially resolved backscattered intensity with a resolution on the order of a few . Polarization-sensitive OCT (PS-OCT) is a functional extension that takes advantage of the additional polarization information carried by the reflected light.4, 5 Thereby, PS-OCT can reveal important information about biological tissue that is unavailable in conventional OCT. Several ocular structures can change the light’s polarization state and are therefore of interest in terms of PS-OCT imaging: in the retina, the nerve fiber layer6 and Henle’s fiber layer7 are birefringent, which is of great importance for glaucoma diagnostics.8 The optic nerve head is surrounded by the birefringent scleral rim,9 which might be used as a landmark in studies of optic disk anatomy. Furthermore, a polarization scrambling layer is located near the retinal pigment epithelium (RPE) which might be useful for diagnostics of age-related macular degeneration (AMD).9, 10 In the anterior segment, the cornea is birefringent,11 which might be of interest for analysis and diagnosis of diseases like keratoconus and corneal scars.12 The iris has polarization preserving anterior layers and a polarization scrambling pigment epithelium.13
Different methods of PS-OCT have been reported in the literature. While early work on PS-OCT4, 5 measured only reflectivity and retardation of a sample, in recent years, many proposals have been made to extract more information, e.g., on Stokes vectors,14, 15 Müller16 and Jones matrix17 distribution, birefringent axis orientation,14, 18 and diattenuation.19, 20, 21 We developed a method that combines the PS low coherence interferometry setup first devised by Hee 4 with a phase-sensitive recording of the interferometric signals in the two orthogonal polarization channels,18 thus allowing us to measure the three most important parameters, reflectivity, retardation, and birefringent axis orientation, simultaneously. We implemented our PS-OCT technique in several types of OCT systems: classical A-scan-based time domain OCT (TD-OCT),18 transversal scanning TD-OCT,13 and spectral domain OCT (SD-OCT).22 The functionality of the systems was demonstrated in several different ocular tissues: cornea,12, 23 iris,13 and retina.9, 10, 22
Compared to other methods that probe the sample with light beams of different polarization states, our method has the advantage of requiring only a single input polarization state, which makes the instrument simpler and faster. However, our method also has a drawback: it requires that the beam incident on the sample is in a well-defined (circular) polarization state. While this is no problem for probing superficial tissues like cornea or skin, it can be a problem if the tissue of interest is located behind a birefringent layer. This is exactly the case in retinal imaging: the retina is imaged through the birefringent cornea. In this case, the beam probing the retina is no longer circularly polarized but in an elliptical polarization state defined by the corneal birefringence. The result is that the measured polarization patterns will be distorted.
A similar problem arises in scanning laser polarimetry (SLP),6 a technique that measures the thickness of the birefringent retinal nerve fiber layer (RNFL) in the region around the optic nerve head for purposes of glaucoma diagnostics.8 It has been shown that a careful compensation of the individual corneal birefringence by a variable retarder is required to provide reliable RNFL thickness measurements by SLP.24
It is the purpose of this paper to report on a method for corneal birefringence compensation that can be used with our single-input-state PS-OCT technique. Our method measures the polarization state of light backscattered from the retinal surface to calculate the corneal retardation and axis orientation. This information is used to compensate for the corneal birefringence by a software algorithm. The method is equally suited for TD- and SD-OCT. We demonstrate the method in combination with a state of the art SD-PS-OCT instrument in a test target measured through a well-defined retarder and in the macula and nerve head region of the living human eye.
Details of our SD-PS-OCT system have been published elsewhere.22 This section contains a brief summary and a description of the polarizing elements in the interferometer to explain the principle of corneal birefringence compensation.
Figure 1 shows a sketch of the optical setup. It is based on the original polarization-sensitive low-coherence interferometer first described by Hee 4 Light emitted from a superluminescent diode (SLD) (fiber pigtailed, collimator focal length , center wavelength , bandwidth ) illuminates, after being vertically polarized, a Michelson interferometer, where it is split by a nonpolarizing beamsplitter (NPBS) into an object and a reference beam. The reference light passes through a variable density filter, a dispersion compensation unit, and a quarter wave plate (QWPr), oriented at , and is reflected by a mirror. After double passage of QWPr, the orientation of the polarization plane is at to the horizontal, providing equal reference power in both channels of the polarization sensitive detection unit. The sample beam passes quarter wave plate QWPs (oriented at ), which provides circularly polarized light to the sample via an galvanometer scanner.
After recombination of the two beams at the NPBS, light is directed toward a polarization-sensitive detection unit, where it is split into orthogonal polarization states by a polarizing beamsplitter. The two orthogonally polarized beam components are guided by polarization maintaining fibers toward two identical spectrometers that record the spectral interferograms corresponding to the horizontal and vertical polarization channels.
Each spectrometer consists of a reflection grating (1200 lines/mm), a camera lens with a focal length of , and a 2048-element line scan CCD camera (Atmel Aviiva M2 CL 2014) with a resolution of 12 bits per pixel. The cameras were operated with an integration time of per A-scan, which provided a sensitivity of at a power of onto the sample. The sensitivity decay was along the measurement range of . To reduce the amount of data, only 1024 pixels of each camera where read out. Data post-processing was performed as previously described.22
Corneal Birefringence Compensation
The polarization state of a fully polarized light beam traversing polarizing optical elements can be calculated by the Jones formalism.25 Each light beam is characterized by a two-element Jones vector whose top and bottom elements describe the beams horizontal and vertical electric field components, respectively. Each optical element is described by a Jones matrix. The elements of these vectors and matrices are generally complex. If a light beam traverses optical elements, the Jones vector of the resulting output beam is found by multiplying the input Jones vector by the Jones matrices of the optical elements in the order they are traversed by the beam.
The Jones vectors of reference and sample beams corresponding to Fig. 1 were previously analyzed in detail.18, 26 It was shown that the reference beam, after double-passing the elements in the reference arm, is in a plane polarization state with the polarization plane oriented at to the horizontal, thus providing equal reference power to both channels, with no phase delay between these components. The sample is illuminated by a circularly polarized beam exiting from the QWPs. The Jones vector of the sample beam, after double-passing the entire sample arm, is:is the Jones matrix of the QWPs, is the Jones matrix of a sample with retardation and (cumulative) optic axis orientation , and is the reflectivity of the sample. On deriving Eq. 1, a vertically polarized input state with and input field amplitude equal to 1 was assumed, and diattenuation was neglected.
The right side of Eq. 1 shows that the exiting sample arm Jones vector is complex valued, the magnitude depending on , the phase terms depending on and . The phase difference, however, depends only on . Therefore, the information on sample retardation and axis orientation are decoupled: while can be obtained from the amplitude ratio of the interferometric signals, is entirely encoded in their phase difference .
In the case of SD-OCT, amplitude and phase information are easily obtained. Since real valued intensity data are (inverse) Fourier transformed, a complex signal is directly obtained22:is the spectral intensity as a function of wavenumber , the structure term, the depth coordinate, the (optical) path length difference between reference and sample beams, the amplitude, and the phase of the interferometric signal in real space. The indices and denote the horizontal and vertical polarization channels, respectively. From the analytic signals on the right side of Eq. 2, the quantities , , and can be directly obtained: for and for .
Let us now assume that the sample is not directly illuminated by the circularly polarized beam exiting the QWPs, but, as in the case of retinal imaging, at first traverses an additional birefringent element, the cornea, that has a Jones matrix . In this case, Eq. 1 is changed to:is the sample arm output Jones vector measured through the cornea. A simple decoupling of and is, in general, no longer possible, preventing a direct solution for and as in Eqs. 4, 5.
If the birefringence parameters and of the cornea are known, the problem can be solved. One possible solution is the method used in SLP, where a variable retarder (VR) with a Jones matrix is placed in front of the eye. Applied to our PS-OCT system, this would mean that Eq. 6 is expanded and we obtain the sample arm output Jones vector measured through the VR and the cornea:and have to be measured. In SLP, this is done by setting and analyzing the polarization pattern of light backscattered from the macula region of the human eye.27 A very characteristic so-called bow-tie pattern is observed that is caused by the combined birefringence effects of the Henle’s fiber layer and cornea. From the orientation of the bow tie, and under the assumption that the macula is normal, and can be derived. Once and are known, the birefringence parameters of the VR are set to values and . In the this case, the Jones matrices and cancel each other. In case of our PS-OCT system, Eq. 7 is now reduced to Eq. 1, i.e., the influence of the corneal birefringence is eliminated, and a corrected output Jones vector is obtained from which the retinal birefringence parameters can directly be derived by Eqs. 4, 5.
To avoid the additional complexity of a variable retarder, and to be able to compensate corneal birefringence parameters that vary with the corneal position through which the retina is imaged, we use a software-based solution that is related to the VR method used in SLP. In a first step, two-dimensional (2D) B-scan images or three-dimensional (3D) data sets of the retina are recorded through the birefringent cornea. From the measured Jones vectors , , and , values are obtained by Eqs. 4, 5. According to Eq. 6, these values correspond to a combined retarder consisting of cornea and retina and therefore cannot be used directly. However, and obtained at the surface of the retina correspond to those of the cornea ( and ). These surface values are extracted and used to model the situation corresponding to Eq. 7 in software, and new, artificial sample arm output Jones vector is calculated:to to indicate that the measured sample is now the retina. is therefore the Jones matrix of the retina whose birefringence parameters and are to be determined. If is chosen so that and , cancels , and the result of Eq. 8 reduces to Eq. 1. The inputs to the right side of Eq. 8 are well known, and can be calculated as a function of position in the retina. From the thus-calculated Jones vector, amplitudes and phases are extracted and and are obtained by Eqs. 4, 5.
It should be mentioned that other computational algorithms for compensating the corneal birefringence are conceivable. While a direct subtraction of the corneal retardation is possible only if the axes of cornea and retina are parallel, algorithms based on Stokes vector analysis or Poincaré sphere rotation might be used. Such algorithms could be adapted from those published for PS-OCT using multiple input states.15, 21
In a first step, the compensation method was demonstrated in a test target consisting of a lens and a stack of birefringent tape resembling the retina. In front of the lens, a liquid crystal variable retarder (LCVR) was mounted to simulate measurements through a birefringent cornea.
We recorded B-scan images of the artificial retina with the LCVR set to zero retardation and to a value mimicking the corneal birefringence (measured value: or for our wavelength of ). Figure 2 shows the result. Figure 2a is a reflectivity image. The individual layers of the stack of tape are clearly seen. Figures 2b, 2c, 2d show retardation images, recorded with (2b), (2c), and , but with cornea compensation [Fig. 2d], as described in Sec. 2.2. In Fig. 2b, near-zero retardation is observed at the surface of the sample, increasing gradually with depth. In the uncompensated measurement through the retarder [Fig. 2c], the surface already shows higher retardation that further increases with depth. To obtain and of the artificial cornea, we extracted the surface of the sample by a peak detection algorithm working on the reflectivity image and calculated and from the interferometric signals corresponding to the peaks found in the A-scans. For the compensation algorithm, we used mean values of and derived from all A-scans. After compensation, the surface shows again blue color [Fig. 2d], indicating near-zero retardation.
For a better quantitative evaluation of the results, we calculated histograms of the values obtained at the sample surface and histograms of corresponding to the second band of the tape. (Because in the case of , the interferometric signal in the vertically polarized detection channel is very low, the phases are rather random and no meaningful axis can be derived; therefore, we used the signal from the second band where a better signal strength is obtained.) Fig. 3 shows the results.
Figures 3a, 3b, 3c show the surface retardation values derived from Fig. 2b, 2c, 2d, respectively. In Fig. 3a, the peak of the histogram is at , very low as expected from the fact that zero retardation was adjusted. (The noise and imperfect polarization components cause the residual value). When switching the retarder to the value mimicking the cornea, a larger surface value of is observed [Fig. 3b]. After compensation [Fig. 3c], the value is again close to zero ; however, due to noise, the distribution is somewhat broadened. The axis histogram obtained without retarder [Fig. 3d] has a peak at . On measuring through the retarder, the value is changed to [Fig. 3e]. After compensation, the value is [Fig. 3f], close to its original value. These results in the test target demonstrate that our software compensation method works well to restore the original birefringence patterns.
In the next step, we performed PS-OCT imaging in a healthy eye of a human volunteer, after full informed consent was obtained. We recorded 3-D data sets of the macula area and cross-sectional circumpapillary scans (nerve head area). A 3-D data set consists of 60 B-scans, each B-scan containing 1000 A-scans. The scanned area covers , or . A circumpapillary scan consists of 4000 A-scans. In each B-scan and circumpapillary scan, the retinal surface and the inner/outer photoreceptor segment boundary (IS/OS boundary) were extracted by searching for the corresponding intensity maximum, and from the interferometric signals corresponding to these maximum values, retardation and axis orientation values were determined.
Figure 4 shows the uncompensated data obtained in the macula. Figures 4a and 4b show the retardation and slow axis orientation (referenced to the vertical instrument axis), respectively, measured across the retinal surface at the inner limiting membrane (ILM). These maps resemble the cornea’s birefringence, i.e., Figs. 4a and 4b show the distribution of and . Figures 4c and 4d show the distribution of and at the IS/OS boundary, respectively. These patterns resemble the combined birefringence of the cornea and retina. For comparison, Fig. 4e shows a depth-integrated reflectivity image of the same area, and Fig. 4f a conventional (reflectivity) B-scan across the fovea centralis, where individual retinal layers can be identified. Since the RNFL is rather thin in the imaged area, the only retinal structure that has a noticeable birefringence is the Henle’s fiber layer. In this layer, the fibers are arranged in a pattern running radially away from the fovea centralis. Therefore, the Henle’s layer acts as a retarder with an axis orientation that changes linearly with azimuth angle around the fovea (the slow axis is parallel to the fiber orientation). If such a retarder is combined with a retarder of approximately constant axis orientation and retardance (as the cornea), the combined retardation pattern will show areas of increased retardation (where the axes of cornea and Henle’s layer are parallel) and areas of reduced retardation (where the axes are orthogonal). The result is a bow-tie pattern well known from SLP.27 This bow-tie pattern is clearly visible in Fig. 4c. The dark arms of the bow-tie pattern are parallel to the fast axis of the cornea, the bright arms are parallel to its slow axis. The pattern orientation indicates a corneal slow axis oriented at nasally downward near the foveal center. The interpretation of the uncompensated axis orientation pattern in Fig. 4d is not straightforward.
The polarization states measured at the retinal surface were used to compensate for the corneal birefringence, as described in Sec. 2.2. As can be seen in Figs. 4a and 4b, and are not constant across the entire image. There is a slow variation across the images, overlaid by noise. To allow for the slow variation of the surface states and to reduce the noise, we used a floating average based filter algorithm to obtain the surface states that entered into our software compensation algorithm [Eq. 8]. The window size of the filter was data points, corresponding to .
Figure 5 shows the compensated images corresponding to Fig. 4. Figures 5a and 5b show retardation and slow axis at the surface, Figs. 5c and 5d the values at the IS/OS boundary, respectively. The values obtained at the surface of the retina [Fig. 5a] are now rather constant across the image and lower than in the uncompensated case. The slightly higher values observed near the center of the fovea are probably caused by the rather poor signal intensity obtained in that area at the retinal surface (due to the curvature of the retinal surface near the foveal pit, the backscattered light in this area is usually rather weak). A low signal power is known to distort PS-OCT signals, usually biasing very low toward higher values in this type of PS-OCT configuration.26 The axis orientation at the surface is rather random, as expected if the corneal influence is removed (no axis is defined in that case).
To better quantify the effectiveness of corneal compensation, we calculated histograms of and as measured at the retinal surface for the uncompensated case [Figs. 6a and 6b ] and for the compensated case [Figs. 6c and 6d]. Before compensation, the width of the distribution is broad, stretching from to [Fig. 6a]; after compensation, a rather narrow distribution with a peak at is observed [Fig. 6c]. Ideally, the compensated should be zero; however, this value is not reached because of noise. (The algorithm used to calculate [Eq. 4] always provides positive values, so the mean will be larger than zero for a noisy signal.) The distribution shows a very pronounced peak in the uncompensated case [Fig. 6b], i.e., a rather well defined corneal slow axis at (referenced to the vertical instrument axis). This corresponds to an orientation of nasally downward, in good agreement with the obtained from the bow-tie pattern. After compensation [Fig. 6d], is much more evenly distributed, however, the compensation is not perfect.
At the IS/OS boundary, the bow-tie pattern is removed after compensation [Fig. 5c], as expected. A slightly increased retardation around the central fovea is observed, probably indicating the birefringence of the Henle’s layer. The compensated axis orientation image now shows a pronounced pattern around the fovea: two full-color oscillations from blue over green, yellow to red, corresponding to an axis change of , are observed along a circle around the fovea, indicating the radially oriented Henle’s fibers.
Figure 7 shows a circumpapillary scan with a diameter of . Figure 7a is a reflectivity image, Figs. 7b and 7c show the uncorrected retardation and axis orientation, respectively. In the thick superior bundle of the RNFL (left) an increase of retardation is observed with depth [color change from dark blue to green; Fig. 7b]. A similar effect is expected in the inferior RNFL bundle (right); however, due to corneal birefringence varying with scan position, an increased retardation is already observed at the anterior surface of the retinal tissue, distorting the expected retardation pattern. The axis orientation is expected to vary by a total of as we go along the circumpapillary scan (because the nerve fibers emerge approximately radially from the nerve head). This would mean two full-color oscillations from blue over green, yellow, orange, to red along the scan. Although a color change can be observed, some colors are missing.
After software compensation of the corneal birefringence (the surface state was averaged over a window size consisting of 400 adjacent A-scans), the artifacts are largely removed: the retardation measured at the surface of the retina is low (blue color) along the entire scan [Fig. 7d], and the axis orientation now shows the expected two full-color oscillations from left to right [Fig. 7e].
The PS-OCT method we use is based on the assumption that the sample is illuminated by a circularly polarized beam. While this condition can usually be met with a free-space PS-OCT setup and samples that are directly accessible, it is usually not fulfilled in retinal imaging where the sample beam first traverses the birefringent cornea. This is no problem if only qualitative polarization information is required, e.g., to identify a depolarizing (or polarization scrambling) layer as the RPE. However, other applications require quantitative information on tissue birefringence, e.g., glaucoma diagnostics via measuring RNFL retardation. In this case, the birefringence of the cornea has to be compensated.
In SLP, a similar problem arises, and it was solved by placing a variable retarder in front of the cornea that cancels the corneal birefringence. The advantage of this method is that it is straightforward and fast: after one measurement of corneal birefringence in the macular region, the variable retarder is set appropriately and the nerve head area is scanned. In principle, this method can also be used for PS-OCT. However, the method has some disadvantages: it requires additional hardware and an additional calibration measurement (in the macula), and it assumes that the birefringence of the cornea is constant, regardless of the scanning position.
In a previous paper,23 we have shown that this is not necessarily the case. If the sampling beam deviates from perpendicular incidence, the birefringence encountered by the beam in the cornea changes. The assumption of constant corneal birefringence therefore requires a very precise alignment of the instrument with the measured eye. Furthermore, certain diseases such as corneal scars and keratoconus can cause a very irregular birefringence pattern in the cornea.12, 28 To allow for these variations of corneal birefringence with scan position, we decided to use a software-based compensation method. This method models the hardware-variable retarder by an additional Jones matrix whose components are derived from the polarization state measured at the surface of the retina. Since the surface state can be measured locally, the corneal birefringence variations with scanning angle can be accommodated. However, because of signal noise, a certain averaging is required to determine reliable surface states. The window size we used was for 3-D data sets and 1/10 of the circumpapillary arclength for the 2-D circumpapillary scans. The latter is coarser than the averaging used in a previous study,29 where the circumpapillary arc was divided into 48 sectors in which tissue birefringnece was obtained from an average of 32 adjacent A-scans. However, in our case, only the surface state is averaged, not entire A-scans. (The retardation is calculated A-scan wise.) Since the surface state varies only slowly [at least with normal corneas; see Figs. 4a and 7b], we think that this coarser averaging is justified because surface speckle is more efficiently removed. However, our averaging window might be too coarse for measuring through a cornea with keratoconus. (In this case, the aberration-induced image degradation might anyway prevent meaningful OCT imaging.)
By using the retinal surface polarization state as input for corneal birefringence compensation, our method has similarities to the Stokes vector–based method that uses the polarization state observed at the retinal surface as a reference state.29, 30 Compared to that technique, our method has the advantage of requiring only one input polarization state. However, the disadvantage is that our method assumes that the retarder anterior to the retina has a value between 0 and . While this is true if there is only the cornea as an anterior retarder, this assumption is in general not true if additional retardation is introduced by other optical elements, e.g., a single mode fiber. Therefore, our method cannot easily be extended to a fiber-based PS-OCT setup.
A disadvantage of software-based methods over hardware solutions is certainly the increased computational effort. The compensation of a 3-D data set such as that in Fig. 5 presently takes several minutes (Pentium 4, , LabView code). However, since our software code consists essentially of matrix multiplications performed for each data point, optimized parallel code could greatly reduce computation time.
Financial support from the Austrian Science Fund (FWF Grant No. P16776-N02) is acknowledged.