Three-dimensional polarization sensitive OCT imaging and interactive display of the human retina.

Polarization sensitive OCT has recently been shown to provide tissue specific contrast, enabling direct identification of retinal layers based on the intrinsic properties of their interaction with light. However, the capabilities of displaying and analyzing 3D datasets in scientific publications were rather limited. Within the framework of the Interactive Science Publishing project, we present new ways of displaying and analyzing 3D sets of various polarization parameters recorded in healthy and diseased human retinas. These datasets can be interactively explored by the reader. Furthermore, we provide data of the 3D distribution of backscattered Stokes vectors to allow the reader to develop and test their own data processing algorithms.


Introduction
Optical coherence tomography (OCT) [1][2][3] is now a well-established tool for highresolution cross sectional imaging of the human retina. While the original OCT technique was based on mechanically scanning a reference mirror to perform A-scans in time domain, spectral domain (SD) OCT [4][5][6] has caused a paradigm change in the OCT community after the discovery of its huge sensitivity advantage [7][8][9]. Various applications of SD-OCT to image the human retina with high speed and high resolution in 2 and 3 dimensions have been successfully demonstrated [10][11][12].
Conventional, intensity based OCT achieves a high depth resolution of a few μm and can resolve several retinal layers, however, it cannot directly differentiate different tissues. Polarization sensitive (PS) OCT takes advantage of the additional polarization information carried by the reflected light [13,14] and can thereby reveal important information on the tissue that is unavailable in conventional OCT. Thereby, PS-OCT can generate tissuespecific contrast and also provide quantitative information.
Different methods of PS-OCT have been reported in literature. While early work measured only reflectivity and retardation [13,14], newer work extended the measurements to various other parameters like Stokes vectors [15,16], Müller [17] and Jones matrix [18], birefringent axis orientation [19] and diattenuation [20][21][22]. We developed a method that combines the PS low coherence interferometry setup first devised by Hee et al. [13] with a phase sensitive recording of the interferometric signals in the two orthogonal polarization channels [19], thus allowing to measure three parameters, reflectivity, retardation, and birefringent axis orientation simultaneously, requiring only a single measurement per sample location. Originally developed for time domain OCT, the technique was later adapted to SD-OCT, enabling the acquisition of 20000 A-lines/s with a sensitivity of 98 dB [23].
Previous work on SD PS-OCT has already demonstrated its suitability for acquiring highresolution 2D cross sections and 3D datasets. However, the possibilities to display 3D datasets were limited to fly-through movies and movies of animated volume renderings. The full information present in the 3D datasets was not always accessible or difficult to comprehend. Within the framework of the Interactive Science Publishing pilot project, it is the purpose of this paper to present the first full 3D PS-OCT datasets of the human retina that can be interactively investigated and visualized by the reader in various ways, thus fully exploiting the information present in these 3D datasets. Examples of healthy retinas and retinas with different disorders are presented.

Methods
The principles and technical details of the SD PS-OCT setup used in this work are published elsewhere [23]. Here is a brief summary: Light emitted from a super luminescent diode (Superlum, Moscow; center wavelength 839 nm, FWHM bandwidth 53 nm) illuminates, after being vertically polarized, a bulk optics Michelson interferometer, where it is split by a non polarizing beam splitter (NPBS) into a sample and a reference beam. The reference light transmits a quarter wave plate (QWP) oriented at 22.5°, and is reflected by a mirror. After double passage of the QWP, the orientation of the polarization plane is at 45° to the horizontal, providing equal reference power in both channels of the polarization sensitive detection unit. The sample beam passes a QWP oriented at 45°, which provides right-handed circularly polarized light to the sample. An x-y galvanometer scanner scans the beam over the retina (fast scan direction: x = horizontal; 1/e 2 diameter at cornea ~ 2 mm; diffraction limited focal spot size at retina ~ 12 μm).
After recombination of the two beams at the NPBS, light is directed towards a polarization sensitive detection unit, where it is split into orthogonal polarization states by a polarizing beam splitter. The two orthogonally polarized beam components are coupled into two polarization maintaining fibers (PMFs) and directed towards two separate spectrometers which are designed identically, consist of similar components, and have to be aligned with sub-pixel accuracy with respect to each other (for details of alignment requirements, see ref. [23]). Each spectrometer consists of a reflection grating (1200 lines/mm), a camera lens with a focal length of 200 mm, and a 2048 element line scan CCD camera with a pixel size of 14 × 14 μm 2 (Atmel Aviiva M2 CL 2014). To reduce the data transfer rate, only 1024 pixels of each camera were read out (for better image quality, the data were expanded again to 2048 pixels per spectrum by zero padding in postprocessing). The maximum line rate of the camera is 29 kHz, and via camera link and a high speed frame grabber board (PCI 1428 National Instruments) data could be transferred continuously to a personal computer. The resolution of the camera is 12 bit per pixel. The sensitivity of our system was 98 dB with an integration time of 50 μs and a power of 700 μW onto the sample. The sensitivity decay, due to the finite spectrometer resolution, was 14 dB (equal in both channels) along the measurement range of 3 mm.
Our system was operated at an A-scan rate of 20k A-lines/sec. 3D data of the human retina, covering a scan field of 15° × 15° and consisting of 60 B-scans (1000(x) × 1024(z) pixels) are acquired within ~ 3 seconds. Imaging depth (in air) is ~ 3 mm. The theoretic depth resolution within the retina (assuming a refractive index of 1.38) is 4.5 μm. Pixel spacing of original datasets is ~ 4.5 μm (x) × 75 μm (y) × 3 μm (z, measured in air). Because of a trigger delay between x-scanner and B-scan data acquisition, images were truncated by 50 pixels in x-direction at one side of the image, leaving 950 transverse pixels (or ~ 14.25°). To improve 3D data processing and interactive image manipulation speed with OSA ISP software, the data were further reduced (after all the image processing steps described below) before generating the final OSA ISP data files: in x-direction, data were downsampled by a factor of 2; in z-direction, data were downsampled by a factor of 2, then the images were cropped in z-direction to remove areas that contain no information (vitreous, noisy areas below the choroid). The final z-extension of the images was 1.5-1.8 mm (in air, details are given in figure captions). Table 1 summarizes pixel spacings, pixel counts, and approximate physical extensions of the OSA ISP datasets.
After data collection the following data pre-processing steps were performed: At first fixed pattern noise, originating from the camera readout, was eliminated. This procedure consists of subtracting a mean spectrum (averaged over 1000 A-scans) from each spectral dataset, inverse Fourier transforming the dataset, removing two remaining sharp frequencies generated by the camera, and Fourier transforming the data back to obtain a spectrum free of fixed pattern noise. Afterwards, zero padding to 2048 pixels was performed, the data were rescaled into k-space (including residual dispersion compensation), and the inverse FFT of both signals was calculated, finally providing the pre-processed A-scan data of the horizontal and vertical polarization channels. Prior to calculation of the polarization data, the influence of anterior segment birefringence was compensated by a numerical method [33].
In a next step, datasets of polarization independent reflectivity R, retardation δ, and optic axis orientation θ were derived from the corrected horizontal and vertical channel data as previously described [19,23]: with A being the amplitude and Φ the phase of the respective channel, and the indices H and V denoting the horizontal and the vertical polarization channel, respectively. ΔΦ = Φ V − Φ H is the phase difference between the two channels. The unambiguous measurement ranges are 90° for δ and 180° for θ.
Apart from retardation, there is another light-tissue interaction mechanism of considerable interest for retinal imaging: We recently demonstrated that the RPE scrambles the polarization state of backscattered light, i.e., acts as a depolarizing layer [27]. This can be used for identifying and segmenting the RPE [34]. Since OCT is a coherent method, the light detected from a single speckle is always fully polarized, i.e., the degree of polarization (DOP) cannot be directly measured by OCT [35]. However, depolarization can be observed in PS-OCT images because it leads to uncorrelated polarization states of adjacent speckles, i.e., their polarization state varies randomly. In a recent paper, we defined a new quantity, the degree of polarization uniformity DOPU [34] that is closely related to the DOP and measurable by OCT. DOPU exploits the random variation of the polarization state of adjacent speckles and is obtained in the following way: At first, we calculate the Stokes vector S of each pixel: (4) where I, Q, U, V denote the four Stokes vector elements. Then we average Stokes vectors over adjacent pixels by calculating the mean value of each Stokes vector element within a floating rectangular window (size 15(x) × 6(z) pixels, or ~70(x) × 18(y) μm) and derive DOPU within this window by the following equation: (5) where the indices m indicate mean Stokes vector elements. As can be seen by Eq. (5), DOPU is closely related to DOP well known from polarization optics. However, we want to emphasize the statistical origin of DOPU, i.e., it can only be derived by local averaging of Stokes vector elements. It might therefore also be regarded as a spatially averaged DOP and is closely related to the apparent degree of polarization obtained by temporal averaging [35] and to the quantity that describes the local correlation of polarization states and was recently used for detection of multiply scattered light by OCT [36]. (It should be mentioned that DOPU is different from the degree of circular polarization (DOCP) that was recently discussed in the context of PS-OCT [37]. Since DOCP is based on only one Stokes vector element (V), it does not really provide a measure for the degree of polarization; e.g., retardation can transform a circularly polarized beam into a linear polarized beam with DOCP = 0 but DOP and DOPU still equal to 1.) In case of a polarization preserving or birefringent tissue, the value of DOPU is approximately 1, in case of a depolarizing layer, DOPU is lower than 1. To avoid erroneous data points caused by noise (which also gives rise to random Stokes vector elements), we first apply a thresholding procedure based on the intensity data to gate out areas with low signal intensity (threshold value: 3 times local intensity noise). To reduce computation time, software compensation of corneal retardation can be omitted for the calculation of DOPU (since randomness of polarization states can be judged also from uncorrected data).
Prior to generation of OSA ISP ready data, motion artifacts were corrected by correlation of B-scans (using ImageJ module "StackReg" [38]). Reflectivity data are displayed on a logarithmic gray scale. Polarization data δ, θ, and DOPU are displayed on a false color scale. For display of polarization data in 3D volume renderings and in cross sectional images through the 3D datasets by OSA ISP software, intensity and polarization datasets are fused to a combined dataset where the polarization values are encoded by color while the intensity data form the alpha channel (i.e., the opacity of a data point corresponds to the intensity) [39]. This method of display ensures that areas of low reflectivity appear transparent in a volume rendering, allowing an unobstructed view of the polarization data in high-reflectivity areas. Furthermore, erroneous polarization data caused by noise in low-intensity areas are gated out of 2D cross sections by this novel display method.
More than 200 eyes of healthy subjects and patients with various diseases were imaged by SD PS-OCT. All measurements were approved by the ethics committee of the Medical University of Vienna and followed the tenets of the Declaration of Helsinki. Informed consent was obtained from the subjects and patients after explanation of the nature and possible consequences of the study. In this paper, a few selected cases are shown to demonstrate the 3D imaging and display capabilities. Full 3D datasets are provided to allow the reader to interactively view and explore the datasets, vary aspect angles, change transfer functions and coloring schemes. Finally, Stokes vector data with full resolution that were only pre-processed (i.e., prior to motion correction, cornea compensation, downsampling) are provided in the Appendix to give readers the opportunity to develop and test their own algorithms for processing polarization data and compare their results with the images shown here. Figure 1 shows a 3D dataset obtained in the optic nerve head (ONH) region of the left eye of a male healthy volunteer (age 35). The figure is a tableau of snapshots from OSA ISP analysis sessions. Four different data are shown in the different snapshots: intensity ( Fig.  1(a)), retardation ( Fig. 1(b)), axis orientation ( Fig. 1(c)), and degree of polarization uniformity DOPU (Fig. 1(d)). Each snapshot contains four images: volume rendering (top left), en face section (x-y; top right), horizontal cross section (x-z; bottom left), and vertical cross section (y-z; bottom right). (Because of the strongly anisotropic sampling ratio, the resolution is lower in the y-direction).

Healthy subject
The intensity images ( Fig. 1(a)) show the usual features like optic disk cupping with the lamina cribrosa at the bottom of the cup, vessel shadows (x-y image), and various retinal layers and vessel cross sections (x-z and y-z images).
The retardation images ( Fig. 1(b)) show several interesting features that are not discernible in the intensity images. For clearer demonstration of some of the details, the volume rendering was cropped, the boundaries of the volume in Fig. 1(b) (top left) are indicated in the three cross sectional views. The following features can be observed (marked by arrows): the increase in retardation caused by the birefringent inferior RNFL bundle (color change from blue to green), the strongly birefringent scleral ring (rim of the scleral canal), and the thin polarization scrambling RPE layer (appearing in green color). Furthermore, the strongly birefringent sclera can be seen in some areas. (In this subject, the penetration of the probing light into deeper layers is stronger than usual, so the sclera can partly be observed. This gives rise to so-called "atypical scanning laser polarimetry (SLP) patterns" in some subjects imaged with SLP, a well known artifact of SLP imaging technology [40,41]).
Figure 1(c) shows axis orientation data. The anatomically varying axis orientation of the nerve fiber bundles around the ONH can clearly be observed (for a better view of this feature, the x-y cross section (top right) is replaced by a slightly oblique section that is approximately parallel to the retinal surface). Figure 1(d) shows the degree of polarization uniformity DOPU. In most retinal layers, a high value of DOPU is observed (orange to red color). DOPU remains rather constant within the birefringent RNFL, and is only slightly decreased in layers below the RPE. The pronounced decrease of DOPU within the RPE (green color) is clearly observed. It should be noted that reduced DOPU values are also observed in parts of the lamina cribrosa and in the scleral ring. This is probably an artifact caused by the strong and locally varying birefringence that changes Stokes vectors significantly within the evaluation window in these areas. Figure 2 shows a movie of the animated, volume rendered retardation dataset of Fig. 1(b), showing all the above mentioned features from various aspect angles. In addition, all four full 3D datasets can be interactively explored by using the OSA ISP viewer. Figure 3 shows a 3D dataset of the macula of the right eye of a 62 years old male patient with age related macular degeneration (AMD). The clinical diagnosis was choroidal neovascularization with extensive RPE atrophy, visual acuity was 0.1. Figure 3(a) provides intensity data and Fig. 3(b) shows DOPU data. For comparison, Fig. 3(c) shows a fluorescein angiogram. The OSA ISP snapshots contain the same sub-images as those of Fig. 1. The intensity dataset shows large cysts in the retina (cf. cross sectional images) and atrophic areas that can be recognized by increased penetration of the probing light into subretinal tissue (providing a locally increased imaging depth). However, the location of the atrophic lesions is only coarsely discernible, small atrophies are not clearly visible. In the DOPU images, atrophic areas are easily recognized by the missing depolarizing RPE layer (cross sectional images). Moreover, the inclined upwards view towards the bottom of the (cropped) volume rendered dataset (top left of tableau of Fig. 3(b)) clearly shows atrophic areas (yellow-orange) between the green-blue areas where the RPE is still present. Interactive manipulation of the 3D datasets allows these structures and features to be explored in a more comprehensive way than the 2D snapshots. Figure 4 shows a 3D dataset of the macula of the right eye of a female AMD patient (age 77 years) with cystoidic macular edema and scarring. Visual acuity was 0.05. Large cysts and multiple smaller cysts are observed in the neurosensory retina (intensity images, Fig. 4(a)). The regularly stacked structure of the photoreceptor layer at the posterior side of the retina is missing, instead, a strong hyper reflective band is observed at the posterior boundary of the retina. From the intensity images, the type of tissue corresponding to this hyper reflective band is unclear, it cannot be told whether the RPE is still present or not. The DOPU images ( Fig. 4(b)) clearly show that only small islands of depolarizing tissue (RPE; green and blue patches) are left, major parts of hyper reflective tissue show high DOPU values (orangered), indicating that they likely consist of fibrotic scar tissue, conforming to an end stage of AMD. This can be seen in more detail in the movie of the animated 3D DOPU dataset (Fig.  5) and by interactive exploration of the 3D dataset by the OSA ISP viewer. Figure 6 shows a 3D dataset of the retina of a 72 years old female patient with a flat choroidal nevus superior to the nerve head (12 o'clock position; a benign tumor; no change in size for the last years). The intensity dataset ( Fig. 6(a)) shows a normally layered retina in most parts of the field of view. However, a small focal pigment epithelium detachment overlying a solid extrusion of the choroid can be seen (arrow) surrounded by an area of slightly increased backscattering below the retina. The retardation images (Fig. 6(b)) show an increased level of retardation in this area, however, it is unclear from this figure if the observed green color is caused by true retardation or by polarization scrambling. The DOPU images (Fig. 6(c)) clearly show strong depolarization in this area in addition to the thin depolarizing band of the RPE. Figures 6(d) and (e) show additional views of the bottom (inclined upwards and straight upwards) of the 3D rendered DOPU dataset that show the large extension of the depolarizing nevus. Figure 6(f) shows a conventional color fundus photo for comparison. More details can be found by interactively exploring the 3D datasets by the OSA ISP viewer.

Discussion
We have shown that the polarization state of backscattered light can be used as an intrinsic, tissue specific contrast means for retinal imaging. While retardation images exploit the birefringence of the RNFL and the sclera for tissue identification and quantification (with axis orientation providing additional information on the orientation of RNFL bundles), the polarization scrambling properties of pigmented tissue like RPE and nevus can be used to clearly identify these tissues and locate areas of RPE lesions in AMD.
Different versions of PS-OCT have been reported, the main differences are the use of bulk optic versus fiber optic interferometers, and the number of input polarization states required to obtain the polarization information. We use a bulk optics interferometer and illuminate the sample with circularly polarized light which has the advantage that only a single input state is needed for acquiring the polarization data shown in this paper. The drawback of our method is that it is not straightforward to be implemented in fiber optics, and that it relies on the assumption of negligible diattenuation. This assumption, however, seems to be largely warranted because diattenuation measured in various tissues by other multi-input state PS-OCT systems is very small [20][21][22]. The main advantage of using a single input state is, apart from reduced measurement time and system complexity, that the system is not sensitive to phase changes that can occur between adjacent A-scans. Furthermore, it is less influenced by noise originating from uncorrelated speckle fields: PS-OCT methods based on fiber optics usually need at least four data points: two from the tissue surface and two from within the tissue, to eliminate the unknown fiber birefringence, and to derive the true (depth integrated) tissue birefringence corresponding to one point in the tissue. Images computed from these four datasets are more sensitive to noise degradation than images just computed from a single dataset.
The comprehensive 3D acquisition of intensity and polarization data provides new diagnostic opportunities and simultaneously poses new challenges in terms of appropriate data display. On the one hand, display of conventional structural data in the form of reflectivity images and volume renderings is required; on the other hand, an unobstructed view of polarization data is desirable. Our solution to this challenge was to combine intensity and polarization data to a combined dataset that encodes the polarization quantity into color and the reflectivity into opacity. In this way, low-intensity areas that cause large polarization noise (e.g., data from areas in the vitreous, where no light is reflected and therefore the corresponding polarization data consist entirely of random noise) appear totally transparent and provide an unobstructed view on the polarization data originating from true tissue structures.
The new OSA ISP 3D data visualization tool is already very powerful and capable of displaying these combined datasets in various ways, as cross sections of arbitrary orientation and position, as well as in form of volume rendered datasets and movies derived from animations of these volume renderings. While the cross sectional views allow to switch between polarization data and reflectivity data (by right-clicking on a cross sectional image and selecting "Color Mapped -Intensity Modulated" (= polarization) or "Opacity" (=intensity) from the "Color Mapping" sub-menu; for optimum view of the respective modality, the contrast (image transfer function) has to be modified accordingly), the volume rendering of such a combined dataset can presently display only the polarization state. To show a volume rendered reflectivity dataset requires an additional dataset containing just 3D intensity data. This has the drawback of not being able to switch directly between the two views which would allow for detailed comparison of tissue features in polarization and reflectivity under exactly the same viewing conditions (3D angle, magnification, cropping). Future versions of the OSA ISP viewer will hopefully provide this direct switching feature.
Apart from visualization and display issues, also the data processing and the selection of polarization parameters that are to be displayed and analyzed deserve some discussion. In previous papers, we (and others) focused on retardation (and the associated axis orientation) to characterize birefringent tissue like the RNFL. Recently, we introduced the parameter DOPU [34], which is closely related to the DOP and is of great interest for characterizing depolarizing tissue like the RPE. However, other parameters, or other algorithms to derive these parameters, might turn out to be even more relevant or appropriate. For example, the evaluation window used for obtaining DOPU might not be optimal in certain cases (see, e.g., the image artifacts near strongly birefringent areas like the scleral ring, Fig. 1(d)). Other window shapes, 3D windows, or more sophisticated evaluation kernels [36] might provide more reliable results in such difficult situations.
To give the readers the opportunity to develop and test their own algorithms for processing polarization data, we provide the most general form of these data, the 3D distribution of Stokes vectors. Figure 7 shows an OSA ISP snapshot of the Stokes vector dataset corresponding to Fig. 1. The data were only preprocessed in the intensity regime (fixed pattern noise removal, zero padding, re-scaling into k-space, dispersion compensation, Fourier transform to real space, alignment of B-scans to reduce motion artifacts; to avoid mathematically inconsistent Stokes vector elements that can be caused by interpolation between pixels, a simpler motion correction algorithm that only corrects for axial motions by integer pixel numbers was used in this case). No corneal birefringence compensation and no averaging of Stokes vector elements were performed to allow the reader to develop all the polarization state manipulations themselves. The data are encoded in the following way: intensity is encoded in the alpha channel (opacity), the Stokes vector elements Q, U, V are encoded into red, green, and blue colors of RGB color space, respectively (8 bit scales).
We would like to point out that, although variations of the polarization state can be observed by changes of the colors in Fig. 7, a direct interpretation of these colors is not meaningful; instead, extraction of physically meaningful parameters by the readers developing their own data processing algorithms is encouraged. These can be compared to the parameters we extracted (cf. Fig. 1). Since especially the calculation of DOPU data benefits from the maximum possible resolution, we provide Stokes vector datasets of full resolution, i.e., before downsampling, and also before motion correction (the only size reduction is the limitation of the z-extension to 750 data points by cropping areas that contain no information). These datasets can be downloaded via these links, Data001; Data003; Data004; Data006, where Data001 corresponds to the data of Fig. 1, Data003 to those of Fig. 3, etc. Each link opens a web page that contains the Stokes vector elements I, Q, U, and V as a stack of 60 16bit TIFF images that correspond to the 60 B-scans of each 3D dataset. Each image, individually viewed, shows the 2D distribution of the corresponding Stokes vector element, so, e.g., image no. I00 of Data001 shows the first Intensity B-scan of the dataset corresponding to Fig. 1, Q59 of Data003 shows the image of the Q-distribution of the last (60th) B-scan of the dataset corresponding to Fig. 3, etc. Intensity data are on a linear scale; Q, U, and V are also on a linear scale, with the lowest value corresponding to −1, the highest to +1.  PS-OCT retardation image of healthy nerve head. Movie of animated volume rendering (same dataset as Fig. 1(b)) (Media 1).        3D PS-OCT dataset of healthy human retina in the nerve head area (same original data as in Fig. 1). The 3D distribution of the backscattered Stokes vector is shown (View 12). Opacity is encoded by intensity, Stokes vector elements Q, U, V are encoded in red, green, and blue colors of RGB color space (0 corresponds to −1, 255 to +1 for each of the variables that are encoded by color values). Image size: ~ 14.25°(x) × 15°(y) × 1.8 mm(z, in air).