Retinal pigment epithelium segmentation by polarization sensitive optical coherence tomography

We present a new method for identifying and segmenting the retinal pigment epithelium (RPE) in polarization sensitive optical coherence tomography (PS-OCT) images of the human retina. Contrary to previous, intensity based segmentation algorithms, our method uses an intrinsic tissue property of the RPE: its depolarizing, or polarization scrambling effect on backscattered light. Two different segmentation algorithms are presented and discussed: a simpler algorithm based on retardation data, and a more sophisticated algorithm based on local variations of the polarization state calculated from averaged Stokes vector elements. By using a state of the art spectral domain PS-OCT instrument, we demonstrate the method in healthy and diseased eyes.


Introduction
Optical coherence tomography (OCT) [1][2][3] is now a well-established tool for highresolution cross sectional imaging of the human retina. Driven by its unique imaging capabilities, initial applications of OCT concentrated on qualitative imaging. In recent years, quantitative imaging became more and more important, such as measuring the thickness of the retina and of its individual layers, and generating 2D thickness maps of these layers, derived from 3D data sets. An important step for these quantitative diagnostic applications of OCT is the segmentation of retinal layers. Several algorithms have been developed to identify the boundaries of retinal layers [4][5][6][7]. These algorithms are essentially based on variations of backscattered intensity between layers or at the boundaries of layers. Intensity based algorithms are, however, sensitive to various factors like illumination conditions, presence of vessels, layer continuity, or phase washout caused by ocular motions.
Among the most important retinal layers are the photoreceptor layer and the retinal pigment epithelium (RPE), lying in close proximity at the posterior side of the retina. The photoreceptor layer converts the light energy into electrical signals, and the RPE plays a decisive role in the metabolism of the photoreceptors. Any disturbance of the RPE can lead to severe visual impairment, e.g. in age related macular degeneration (AMD) [8]. A direct visualization and segmentation of the RPE would therefore be of considerable diagnostic value. However, segmentation of the RPE by intensity based algorithms can be problematic: OCT images usually show three closely spaced boundaries at the posterior side of the retina that have approximately equal backscattering strength: the boundary between inner and outer photoreceptor segments (IS/OS boundary), the end tips of the photoreceptors (ETPR, or Verhoeff's membrane), and the RPE [9,10]. Because of their close proximity and similar backscattering, they are difficult to distinguish in intensity based images, especially if their normal appearance is disturbed by some retinal disorder.
Polarization sensitive (PS) OCT is a functional extension of OCT that provides information beyond backscattered intensity [11,12]. Different methods of PS-OCT have been reported in literature, providing access to various quantities like retardation [11,12], Stokes vectors [13,14], Müller [15] and Jones matrices [16], birefringent axis orientation [17] and diattenuation [18][19][20]. We developed a method that combines the PS low coherence interferometry setup first devised by Hee et al. [11] with a phase sensitive recording of the interferometric signals in the two orthogonal polarization channels [17], thus allowing to measure three parameters, reflectivity, retardation, and birefringent axis orientation simultaneously. Originally developed for time domain OCT, the technique was later adapted to spectral domain (SD) OCT, enabling the acquisition of 20,000 A-lines/s with a sensitivity of 98 dB [21].
Various applications of PS-OCT to imaging of ocular structures like the cornea [22], the anterior chamber [23], and the retina [9,24,25] have been reported. Using PS-OCT, structures of the ocular fundus could be classified into polarization preserving (e.g. photoreceptor layer), birefringent (e.g. retinal nerve fiber layer (RNFL), Henle's fiber layer, sclera, scar tissue) [9,21,[24][25][26][27][28][29][30], and polarization scrambling or depolarizing (e.g. retinal pigment epithelium (RPE)) [9,21,27,28,30]. Polarization preserving tissue can be recognized in PS-OCT images by constant retardation and axis orientation with depth. Birefringence manifests itself by an increase in retardation with depth, while the axis orientation is usually rather constant or varies on a scale considerably larger than the speckle size. Polarization scrambling tissue can be identified in PS-OCT images by random retardation and axis orientation values, varying from speckle to speckle.
In this paper, we demonstrate the use of PS-OCT for identifying and segmenting the RPE. We use its polarization scrambling property as an intrinsic, tissue specific contrast mechanism that differentiates RPE tissue from other retinal structures. We present two different algorithms, a simple, faster algorithm that uses only retardation information, and a second, more sophisticated method that evaluates local polarization state fluctuations from averaged Stokes vector elements, thus exploiting also phase information. We demonstrate our method in the macula of healthy volunteers and patients with AMD and pseudovitelliform dystrophies and discuss advantages and drawbacks of both segmentation algorithms.

Methods
We used an SD PS-OCT setup, whose technical details are published elsewhere [21]. The instrument records 20,000 spectra (or A-scans) per second in each polarization channel. 3D data of the human retina region, covering a scan field of 15° × 15° and consisting of 60 Bscans (1000(x) × 1024(z) pixels) were acquired within 3 seconds. Imaging depth is ~ 3 mm. The measured system sensitivity was 98 dB. The theoretic depth resolution within the retina (assuming a refractive index of 1.38) is 4.5 μm. The instrument needs only one measurement per sample location to retrieve information on intensity, retardation, and axis orientation.
After standard SD-OCT data preprocessing (fixed pattern noise removal, zero padding, and rescaling of the spectral data from wavelength to wavenumber space) the influence of anterior segment birefringence was compensated by a numerical method [31]. The corrected polarization data are then used for RPE segmentation. We used two different algorithms for this purpose. The first, simpler segmentation algorithm is based on the following idea: Light backscattered from most of the retinal layers is in a well-defined polarization state. Only light backscattered from the RPE is depolarized. This is observed in PS-OCT images as a random polarization state that varies from speckle to speckle, giving rise to randomly varying retardation values. To avoid erroneous segmentation caused by noise (which also gives rise to random retardation values in PS-OCT imaging), we first apply a thresholding procedure based on the intensity data to gate out areas with low signal intensity (i.e., only polarization data from areas where the intensity exceeds the noise level by a factor of 3 or more is further processed for RPE segmentation; the variation of intensity noise with imaging depth obtained from calibration measurements is taken into account in this procedure). Within the areas where the intensity threshold criterion is met, the retardation images are analyzed by use of a floating window that moves over the entire image. (Different window sizes were evaluated, cf. Fig. 3(d); a size of 15(x) × 6(z) pixels, or ~70(x) × 18(y) μm turned out to provide best results and was used to process the data shown in this paper.) We calculate a histogram of retardation values within the window for each window position. From the histogram, we obtain the standard deviation (SD) of the distribution of retardation values. If the SD exceeds a certain value (we used an empirically determined threshold of 16°), the area in the corresponding window is regarded as displaying a random polarization state, and the corresponding tissue is classified as RPE.
The second algorithm is based on Stokes vector analysis and has relations to the degree of polarization DOP. Since OCT is based on a coherent detection method, the light detected from a single speckle is always fully polarized, i.e. DOP = 1 within a single speckle [32]. However, in case of a depolarizing layer, the polarization state of adjacent speckles is uncorrelated, i.e. it varies randomly. This random variation is exploited in the algorithm: After intensity thresholding (see above), we calculate the Stokes vector S of each pixel: (1) where I, Q, U, V denote the four Stokes vector elements, A x and A y denote the amplitudes of the two orthogonal polarization channels x and y, and Δ is their phase difference. Then we average Stokes vectors over adjacent pixels by calculating the mean value of each Stokes vector element within a rectangular window (same size as with first algorithm) and derive a quantity that we call the degree of polarization uniformity DOPU: (2) where the indices m indicate mean Stokes vector elements. The quantity DOPU has relations to the degree of polarization DOP well known from polarization optics, as can be seen by its definition (Eq. (2)). 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 the apparent degree of polarization obtained by temporal averaging [32] and to the quantity that describes the local correlation of polarization states and was recently used for detection of multiply scattered light by OCT [33].
In case of a polarization preserving or birefringent tissue, the value of DOPU should be approximately 1, in case of a depolarizing layer, DOPU should be lower than 1. Similar to the first method, we defined a threshold (DOPU = 0.65, empirically determined), where values below this threshold were regarded as depolarized light, and the corresponding tissue was classified as RPE.
All data acquisition and processing routines were written in LabVIEW (National Instruments). The total computation time of the segmentation algorithms (without preprocessing and cornea compensation), for a 3D data set consisting of 60 B-scans (1000 A-scans each) is 8.3 minutes for algorithm 1 and 31 minutes for algorithm 2 (AMD Athlon 64×2 Dual Core Processor 4800+, 2.41 GHz, 2 GB RAM). Figure 1 shows a PS-OCT B-scan across the fovea of a healthy human eye. In the intensity image of Fig. 1(a), it can be seen that the three posterior layers/boundaries of the retina (IS/ OS, ETPR, RPE) have rather similar reflectivity. However, the polarization sensitive images (b-d) clearly show that the most posterior of these layers, the RPE, changes the polarization state of backscattered light, while the other retinal layers render the polarization state essentially unchanged. The retardation image ( Fig. 1(b)) and the image of Stokes vector element U (1(c)) reveal the random nature of the polarization state backscattered by the RPE (images of Stokes vector elements Q and V have a rather similar appearance and are therefore not shown here). The image of the degree of polarization uniformity DOPU ( Fig.  1(d)) quantitatively shows the extent of the depolarization caused by the RPE, while simultaneously showing the high DOPU values retained by the other layers. It should be noted that only light directly backscattered by the RPE is depolarized. Light that traverses the RPE and is backscattered from deeper structures is still in a well defined polarization state.

RPE segmentation in a healthy fovea
The two algorithms used for RPE segmentation are illustrated with the help of the histograms in Fig. 2. Figures 2(a) and 2(b) show histograms of retardation values obtained from windows within the polarization maintaining IS/OS boundary (2(a)) and within the polarization scrambling RPE (2(b)), respectively. The corresponding window positions are indicated by white rectangles in Fig. 1(b) (window size not drawn to scale). The larger width of the histogram taken from the RPE is clearly observed. The retardation SD is 8.5° in Fig.  2(a) and 18.5° in Fig. 2(b). Figures 2(c) and 2(d) show similar histograms of the Stokes vector elements Q, U, V taken from the window within the IS/OS (2(c)) and within the RPE (2(d)). The wider histograms in the latter case are clearly observed. From these histograms, DOPU values of 0.94 and 0.39, respectively, are derived. Figure 3 illustrates the performance of the two RPE segmentation algorithms. The same data set as in Fig. 1 is used. For segmentation, the window from which the histograms are calculated, floats over the entire image. Areas where the retardation SD or the DOPU value meet the threshold criterion (cf. methods section) are marked by red color. Figure 3(a) was segmented by algorithm 1 (retardation SD), Fig. 3(b) by algorithm 2 (Stokes vectors). Figure  3(c) compares the results of the two segmentation algorithms: areas identified as RPE by both algorithms are shown in blue, areas identified only by algorithm 1 are shown in red, while those identified only by algorithm 2 are shown in green. In this case, both algorithms perform almost identical. Therefore, we used the simpler and faster algorithm 1 in the

Europe PMC Funders Author Manuscripts
Europe PMC Funders Author Manuscripts majority of cases, especially for time consuming 3D segmentation (see also discussion in section 4 for a comparison of algorithm performance in more complicated cases). Figure 3(d) compares the effect of the evaluation window size on the performance of algorithm 1 (the effect is very similar for algorithm 2, therefore, the corresponding results are not shown here). From left to right, window sizes of 4(x) × 2(z) pixels (~19(x) × 6(y) μm), 8 × 4 pixels (~37 × 12 μm), 15 × 6 pixels (~70 × 18 μm), 30 × 12 pixels (~140 × 36 μm), and 45 × 20 pixels (~210 × 60 μm) were used. While small window sizes (left) provide too few points for reliable statistics (causing gaps in the segmented layer), large windows (right) generate coarse segmentations that smear out finer details. Therefore, we used the window size of 15 × 6 pixels (~70 × 18 μm) (center) for evaluating the data sets of this work, as it provides the best compromise. Figure 4 shows results obtained in an AMD patient. A single cross section taken from a 3D data set is displayed. Figure 4(a) is a reflectivity image, 4(b) a retardation image, 4(c) shows the degree of polarization uniformity DOPU, and 4(d) a reflectivity image overlaid with the RPE segmented by the retardation SD based algorithm. The retina shows a large RPE atrophy, and a druse. The polarization scrambling RPE is visible in the retardation and in the DOPU image. Its contour is clearly segmented in the overlay image ( Fig. 4(d)). The RPE is intact on the right hand side of the image, lying anterior to the druse. From the image center to the left, the large atrophy is clearly visible. Some segmented patches (red) in the area below the RPE atrophy might indicate pigmented choroidal tissue [8] (Fig. 4(d)). Figure 5 shows an en face image of the thickness of the structure obtained by the segmentation process (retardation SD algorithm) across the 3D data set (optical thickness measured along A-scans, the white line indicates the position of the B-scan corresponding to Fig. 4). The boundary between the intact RPE and the atrophic area is clearly visible. The small patches of segmented tissue within the atrophic area might correspond to pigmented choroidal tissue (cf. Fig. 4). Figure 6(a) shows a movie of an animated volume rendering of the same 3D data set. The movie shows the backscattered intensity in gray, with the overlaid segmentation in red color. The highlighting of the segmented tissue is activated during the animation, giving it the impression of inserting the red color marking the RPE into the reflectivity data over time; at the end of the movie, the intensity data are switched off, and only the segmented layer remains. Figure 6(b) shows another movie of the same data set where only the segmented RPE is shown. Figure 7 shows a cross sectional image taken from a 3D PS-OCT data set of a patient with a pseudo-vitelliform pattern dystrophy (adult-onset foveomacular vitelliform dystrophy). The intensity image (Fig. 7(a)) shows a lesion with thickening of the retina and an accumulation of tissue of unknown type beneath the photoreceptor layer. The retardation image ( Fig. 7(b)) and the image of DOPU (Fig. 7(c)) show that this tissue is depolarizing and therefore probably consists of RPE cells and pigment loaded macrophage (inflammatory cells). Outside the lesion, the RPE layer is still intact, while at the lesion border atrophic areas can be observed. Figure 7(d) shows that the segmentation algorithm clearly detects the region of accumulated depolarizing tissue. Figure 8 shows an en face image of the thickness of the segmented structure obtained from the 3D data set. The thickening of the depolarizing layer within the pseudo-vitelliform lesion, surrounded by atrophic areas, can be clearly seen. Figure 9 shows two movies from the animated volume rendered 3D data set (same data set as in Fig. 8). In Fig. 9(a), only a part of the volume is visible, showing the intensity in gray, with the overlaid segmentation in red. In Fig. 9(b), only the segmented RPE (of the whole 3D data set) is shown. The increased RPE thickness within the lesion, and the atrophic areas are clearly observed.

Discussion
We have demonstrated the use of PS-OCT for an important application in quantitative retinal imaging: the segmentation of the RPE. Previous algorithms for segmenting retinal layers, including all algorithms that are used in commercial OCT instruments, are based on intensity images. While they can work well if retinal layers are just thickened or bended, they are based on certain assumptions like layer continuity. Therefore, these algorithms frequently have problems if layers are interrupted, discontinuous, or partly missing. E.g., they frequently fail in identifying RPE atrophies, taking Bruch's membrane for the RPE in atrophic areas. The two segmentation algorithms presented here are the first that are not based on backscattered intensity but on polarization properties of backscattered light, i.e. they exploit an intrinsic tissue property. A main advantage is therefore that these algorithms do not rely on any assumptions about the order of layers and do not require that a layer is contiguous. Layer interruptions, e.g. by atrophies or shadowing by blood vessels, do not disturb the algorithm, the RPE is found on both sides of such interruptions.
We demonstrated two different algorithms: a simpler algorithm based on the retardation data, and a more sophisticated algorithm based on averaged Stokes vector elements. While the first algorithm uses only the amplitudes measured in the two polarization channels, the latter also takes the phase difference of the signals into account. Since there is more information available in the latter case, it is worthwhile to consider if this algorithm can provide superior results. Although both algorithms perform very similar in simple cases (healthy eye, Fig. 3) there are cases where the Stokes vector based algorithm is superior. Figure 10 shows a B-scan through the red patch within the RPE atrophy visible on the left side of Fig. 5. Both segmentation algorithms are compared. Figure 10(a) shows the results obtained from algorithm 1 (retardation SD), while in Fig. 10(b) algorithm 2 (Stokes vectors) was used. Figure 10(c) shows the corresponding retardation image, Fig. 10(d) a direct comparison of the segmentation results provided by the two algorithms. Although the two algorithms yield rather similar results in most areas, there are also differences. In the region below the RPE atrophy, where light, due to the absence of the RPE, penetrates down to the strongly birefringent sclera, the algorithm based on retardation SD falsely segments some regions in the sclera (indicated with a blue circle in Fig. 10). The strong retardation in this area obviously fools the algorithm by increasing the width of the retardation SD (which causes the red patch mimicking a thickened RPE in Fig. 5). The algorithm based on Stokes vectors is more reliable in this case. Therefore, there is a trade-off between algorithm speed and reliability. Optimized routines written in a more efficient programming language and parallel processing will overcome this problem and the Stokes vector algorithm might be favored for further work.
An interesting feature of our method is the possibility to quantify the thickness of the pigmented tissue. However, the interpretation of the thickness maps (cf. Figs. 5 and 8) deserves some caution. The RPE consists of a cell monolayer of only a few μm thickness. The layer thickness obtained by the segmentation process is larger than the value expected for the thickness of the healthy RPE. The reason is that the segmented structure is the result of a convolution of the real tissue extension with the evaluation window function. Deconvolution might improve thickness measurements, however, would require a deeper knowledge of the statistics of the polarization state distribution caused by the RPE. There is a trade-off between the number of data points available for retardation and Stokes vector element statistics (and thus the reliability of the algorithms) on the one hand, and achievable resolution on the other hand. The window size presently used (~70 μm (x) × 18 μm (z)) is a compromise. Improvements might be achieved by improving the depth resolution (larger source bandwidth) and/or by 3D windowing (which would require a larger sampling density in y-direction), providing more data points for histogram calculations and/or smaller inplane window sizes and thus better resolution of segmented RPE.
Although the thickness of the normal, healthy RPE cannot be accurately measured, the thickness maps are valuable tools for detection of RPE loss or of accumulations of depolarizing tissue to a thickness well beyond the evaluation window extension (cf. Fig. 8).
None of the commercial OCT instruments presently available can detect such accumulations, and no gold standard for any measurements of pigmented tissue thickness presently exists.
The algorithms presented here might be useful for different tasks of quantitative evaluation of the human retina. On the one hand, mapping of total retinal thickness might be improved. Recent studies with intensity based OCT have shown that errors of retinal thickness measurement were common, with at least some error present in 92% of scans [34]. The inability to clearly identify the RPE in intensity based OCT images makes it rather difficult to create a reliable algorithm for automated retinal thickness measurements. The segmentation of the RPE via its polarization properties might therefore improve retinal thickness measurements. On the other hand, thickening of the RPE can be visualized and the presence and size of atrophies can be quantified with our method and displayed in en face maps. This en face visualization of the RPE can then be compared with other en face imaging methods, like fluorescein angiography or auto fluorescence imaging, and perhaps help to avoid invasive imaging methods at least in some cases. Another important application of RPE segmentation could be automated screening of the ocular fundus for small RPE atrophies that are presently very difficult to detect. These lesions could be important for the treatment decision and outcome prediction of modern antiangiogenetic therapy in neovascular AMD, where the presently evaluated morphometric parameters like retinal thickness are known to be unreliable predictors [35]. A reliable quantification of these RPE lesions would be especially important for follow up studies and treatment control.  Histograms of retardation (a, b) and Stokes vector element (c, d) distributions in selected evaluation windows (cf. Fig. 1(b)). (a, c) Data from window A (IS/OS); (b, d) data from window B (RPE).    En face map of thickness of segmented layer obtained from a retina with AMD and RPE atrophy (same data set as Fig. 4). Color bar: 0 -150 μm.