Segmentation and quantification of retinal lesions in age-related macular degeneration using polarization-sensitive optical coherence tomography

During the last two decades, optical coherence tomography (OCT) evolved as a noncontact technique for imaging transparent and translucent samples with micrometer-scale resolution.¹ OCT is still subject to permanent enhancements in order to increase imaging performance, and a wide range of applications has opened in biomedical imaging.^2- 3 Mainly in ophthalmology, OCT has emerged as an indispensable tool for diagnosis of eye disease and for monitoring during therapy.⁴ Recent advances in terms of sensitivity and imaging speed using spectral domain (SD) OCT allow one to record three-dimensional (3-D) data sets comprising several tens of thousands of depth profiles within a few seconds.^5- 7

In clinical routine, apart from 2-D imaging modalities such as fundus photography, scanning laser ophthalmoscopy, fluorescence angiography, and fundus autofluorescence, OCT is used for retinal imaging in three dimensions. Since reviewing and interpretation of large 3-D data sets may be costly in terms of time and effort, several attempts have been made to ease their analysis. Anatomical characteristics and parameters such as retinal nerve fiber layer thickness or the macular volume can be determined from OCT data sets recorded in the posterior eye segment. Retrieval of such parameters alleviates the extraction and evaluation of possible pathological changes in ocular structures. Therefore, in recent years, segmentation of features in OCT data sets has gained particular interest.^8- 15 Most of the approaches for segmentation of retinal layers are based on reflectivity dissimilarity of different tissue types. However, while the layered structure of the retina is well ordered in healthy eyes, it is sometimes hard to differentiate structures solely based on their reflectivity in eyes with retinal disease. Nevertheless, so-called functional extensions of OCT have been shown to provide additional contrast by being able to record not only reflectivity profiles but also depth-resolved information about flow velocity, changes in the probe light's wavelength spectrum, or the polarization properties of tissues as in Doppler OCT, spectroscopic OCT, and polarization-sensitive (PS) OCT, respectively.^16- 20 In the human retina, birefringent, polarization-preserving, and depolarizing structures can be distinguished by PS-OCT.^21- 25 Recently, segmentation of the retinal pigment epithelium (RPE) was demonstrated by PS-OCT based on that layer’s—in contrast to most of the other layers—depolarizing character.²⁶ Using this tissue-specific contrast, our group showed the potential of PS-OCT to image the RPE in patients with age-related macular degeneration (AMD).^27- 28

AMD is a multifactoral disease that results in progressive vision loss in the central region of the retina (i.e., the macula). AMD is the leading cause of blindness and severe vision loss in the industrialized world's elderly population.²⁹ In the early stages of the disease, metabolic endproducts accumulate in Bruch's membrane and build up a hydrophobic barrier between the posterior vascular layers and the RPE and, therewith, the retina. RPE detachments can result from these deposits, which are called drusen.³⁰

At late stages, AMD is classified as wet (neovascular, exudative) or dry (nonexudative) AMD. The latter branch, which accounts for 80–90% of all diagnosed patients, is characterized by sharply delineated areas of severe depigmentation or apparent absence of the RPE with a minimum diameter of 175 μm. Wet AMD—although accounting for only 10—20% of all patients—is responsible for 80% of vision loss by AMD.³¹ It is characterized by RPE detachments, subretinal or sub-RPE neovascular membranes, scar tissue, subretinal hemorrhages, and related hard exudates (lipids).³⁰

The RPE plays the key role in all types of AMD. In this paper, we demonstrate the potential of PS-OCT as a tool for not only 3-D high-resolution retinal imaging with additional contrast but also for quantitative assessment of retinal disorders in AMD. On the basis of the additional contrast provided by PS-OCT, software algorithms have been developed for quantification of clinically relevant parameters such as the area and volume of drusen and the area of atrophic zones in geographic atrophy. After a description of the PS-OCT device and the segmentation algorithms in Section 2, results from measurements in patients, including a repeatability analysis, are presented in Section 3 and discussed in Section 4.

A spectral domain PS-OCT setup was used for imaging, a scheme of which is shown in Fig. 1. A superluminescent diode (Superlum Diodes, Inc.) centered at 839 nm with a FWHM bandwidth of 58 nm was used as a light source. The light was inserted into a Michelson interferometer in a vertically polarized state using a bulk optics polarizer. After passing a polarizing beamsplitter (PBS), the light was split into a sample and a reference beam by means of a nonpolarizing beamsplitter (BS). In the sample arm, a quarter wave plate (QWP) was used in order to provide circularly polarized light onto the eye. Further sample arm optics included an X-Y galvanometer scanner and a relay of two identical achromat lenses. A probe power of 700 μW on the cornea was used, which is well below the limits drawn by American National Standards Institute (ANSI) and International Electrotechnical Commission (IEC).^32- 33 Half of the light beam returning from the sample arm was again directed into the source arm. There, its horizontally polarized portion was detected by an avalanche photodiode (APD), in order to be able to record scanning laser ophthalmoscope (SLO) images using the same sample arm optics.

In the reference arm, the beam passed a variable neutral density filter, glass prisms for dispersion compensation, and was reflected by a mirror mounted on a motorized translation stage. Light returning from sample and reference arm, respectively, interfered at the nonpolarizing BS. Note that the QWP present in the reference arm of the original PS interferometer scheme devised by Hee et al. ¹⁹ was omitted in our design in order to avoid saturation of the APD by light returning from the reference mirror. This altered layout required an additional half-wave plate (HWP) oriented at 22.5 deg in the detection arm, which rotated the beam's polarization state by 45 deg in front of the PBS in the detection arm, such that equal reference arm power was provided in both orthogonal polarization channels. In a software compensation step, this rotation was compensated for by multiplying the OCT signals from the two channels with the Jones matrix of a HWP oriented at –45 deg. Now, the previously published method for computing PS-OCT signals could be applied.³⁴

From the interferometer exit, the signals of the two orthogonal polarization channels were guided to two identical spectrometer units employing diffraction gratings with 1200 lines/mm and achromat lenses (f = 200 mm). The spectral interference signals were read out simultaneously at two line-scan cameras (e2v AViiVA M2 CL 2014) operating at 20 kHz. Each readout of the two cameras yielded the spectral data necessary to compute not only reflectivity but also phase retardation, axis orientation, and Stokes vectors for one depth profile.

Two operational modes were available for imaging, which were timed and controlled using a field-programmable gate array (FPGA). For alignment of the eye under investigation, SLO images and OCT B-scan images were recorded and displayed in an alternate fashion to provide an overview in three dimentions. SLO and OCT images were recorded in 165 and 20 ms, respectively. In the real-time display, each SLO image was followed by four OCT images. Once the eye was aligned properly, 3-D data sets covering a retinal cube of 6.2×6.7×3.3 mm³ (x×y×z) were recorded. One out of three sampling patterns (64×1024, 128×512, 256×256) could be selected. For this study, only data sets sampled 128×512 (that is, 128 B-scans each consisting of 512 A-lines) were used.

The system was used for imaging in the eyes of patients suffering from various ocular diseases. Compliant patients diagnosed with drusen and advanced dry AMD (geographic atrophy), respectively, were selected and imaged for this study. The study protocol was approved by the ethics committee of the Medical University of Vienna and followed the tenets of the declaration of Helsinki.

Segmentation of the RPE was performed based on its depolarizing (i.e., polarization scrambling), character as reported in Ref. ²⁶. For each B-scan, the signals recorded at the two exits of the PS-OCT instrument were processed using regular SD-OCT processing, including subtraction of an averaged spectrum, rescaling, numerical dispersion compensation, and Fourier transform.³⁵ From the resulting signals, I_1,2(z) = A_1,2(z) exp[iΦ_1,2(z)], the Stokes vector S = (I Q U V)^T was computed for each image pixel. In order to gate out noisy pixels, for further processing, only pixels with intensities above a depth-dependent threshold I_thr(z) = I_thr(z = 0 μm) −z·1.7×10⁻³ dB/μm were considered. I_thr(z = 0 μm) was set to a value of ∼11 dB above a level equal to two standard deviations of the average noise amplitude in each A-scan.³⁶

Each of the Stokes vector elements Q, U, and V was normalized by I and then averaged in a sliding window sized 8(x)×5(z) pixels (∼100×16 μm²). Thereof, the degree of polarization uniformity (DOPU) was calculated as Display Formula

Now, from each DOPU B-scan image, depolarizing structures can be derived by segmenting pixels with low DOPU values, that is, pixels with DOPU values below a threshold DOPU_thr. Segmentation results may vary, depending on the choice of DOPU_thr, and the optimal DOPU_thr value can be different, depending on the respective eye and on the image quality. For the results and figures in this paper, a fixed value of 0.8 was chosen for DOPU_thr.

Figure 2 shows PS-OCT B-scan images of a human retina in the fovea region of a healthy volunteer's left eye. In Fig. 2, the reflectivity B-scan image is displayed where healthy layers can be distinguished based on the intensity of backscattered light. The respective DOPU image is shown in Fig. 2. While most layers—those appearing in reddish color—do not considerably alter the sample beam's polarization state, polarization-scrambling manifests in lower DOPU values displayed in green and blue. Note that, due to the computation of DOPU by averaging in a sliding evaluation window, the spatial resolution is somewhat reduced in this image. By considering only pixels with low DOPU values, depolarizing structures can be segmented as shown in Fig. 2, where low-DOPU-valued pixels are overlaid in red. In Fig. 2, the position of the internal limiting membrane (ILM), which was found by intensity thresholding, is overlaid on the reflectivity image in blue. Furthermore, the position of the pixel with the lowest DOPU value in each A-line (i.e., the pixel exhibiting the highest depolarization), which in the healthy retina is located in the RPE, is displayed in red. A fundus image computed from the OCT data set by axial summation of the intensity for each A-line is shown in Fig. 2, and in Fig. 2, a retinal thickness map derived by calculating the axial distance between ILM and RPE position can be seen. An average refractive index of 1.38 was assumed in the retina.

Drusen were segmented by first locating the actual position of the RPE by its DOPU value. This line, found by tissue-specific contrast, acts as a “backbone” for the following calculation of the position where the RPE should be in a healthy eye. Finally, the drusen are segmented and quantified by calculating the difference between the actual position and the normal “should-be” position.

For each B-scan, the starting point of the calculation is the axial position Z_RPE(i) of the minimum DOPU value in each A-line, that is, the red line in Fig. 2. The index i indicates the index of the A-line in the respective B-scan and Z_RPE(i) is measured as the distance from the zero delay, which is located at the upper border of all shown B-scan images, and Z_RPE(i) increases in the posterior direction. First, missing values of Z_RPE(i) are interpolated in A-scans, where no pixel is below the threshold DOPU_thr. Then, outliers of Z_RPE(i) are cut and interpolated. Outliers are defined as values Z_RPE(i), which differ by >33 μm (i.e., 10 pixels) either from the mean of Z_RPE(i − 1) and Z_RPE(i − 1) or from the mean of Z_RPE(i + 1) and Z_RPE(i + 2), or from both. If Z_RPE(i) is an outlier, then Z_RPE(i) is replaced by the average of its four neighboring values Z_RPE(i + m), where m = −2, −1, +1, +2. The subsequent function is smoothed using a Savitzky–Golay filter (polynomial order: 3, filter length: 6).

The result of these first steps is a function

The new function A_j is used in the following second step to compute the output B_j: Display Formula

Drusen can now be found by computing the difference between actual RPE position Z_RPE and normal RPE position Z_NRPE for each A-line. In the case of very large drusen or drusenoid pigment epithelium detachments in retinas appearing heavily bent in the OCT B-scan image, Z_NRPE may be displaced in anterior direction (i.e., between the actual position of Bruch's membrane and the RPE). This will not affect the calculated drusen area, but the volume may be underestimated. We observed this effect in confluent drusen with a diameter of ∼3 mm. For the results presented in this study, which included patients with drusen classified as AREDS score 2 and 3, this effect was not apparent.

Using PS-OCT, atrophic zones in dry AMD can be detected by detecting holes in the depolarizing RPE layer.

In a straightforward manner, one could simply sum the number of depolarizing pixels [i.e., the red pixels in Fig. 2] for each A-line of a 3-D PS-OCT data set. The resulting thickness map of the depolarizing layer now reveals the existence and absence of polarization scrambling tissue in an en face view.

However, depending on the subject's pigmentation, depolarization may also be apparent in the choroid. Signals from this layer usually are rather weak when the probe beam anteriorly has to penetrate the highly scattering and absorbing RPE. Contrarily, in eyes with geographic atrophy, the missing RPE usually results in higher penetration depths (sometimes down to the sclera) and in stronger OCT signals from the choroid. Dependent on the subject, choroidal tissue exhibits a more or less polarization scrambling character. In some cases, depolarizing spots with DOPU values as low as those in RPE tissue may be noticeable. If, in such eyes, the number of all depolarizing pixels is integrated over depth in order to compute a thickness map of the RPE, then polarization scrambling pixels in the choroid may add to those in the RPE. This can result in mimicking of healthy structures and masking of RPE defects, if the RPE is absent or thinned as in the case of geographic atrophy.

To prevent from integrating choroidal pixels into computing of the thickness map of the depolarizing layer, only depolarizing pixels located in a shallow band close to the photoreceptor layer and Bruch's membrane should be contributing.

For this purpose, the ILM position is used as a starting point. The initial ILM position is smoothed similar to the RPE position in the previous section by using a third-order polynomial SG filter with a window length of 200, yielding the function

By summing the number of depolarizing pixels within the evaluation band along each A-line, an RPE thickness map is generated. In order to assess the size of atrophic zones, the thickness map is binarized and smoothed by removing isolated pixels using an erosion filter. In a last step, the algorithm automatically detects patches of connected pixels, which can now be converted into atrophic areas by scaling them with the known pixel area.

Figure 3 shows results from drusen segmentation in a PS-OCT data set recorded in the right eye of a female patient at age 70. Reflectivity B-scan images are shown in Figs. 3. In Figs. 3, an overlay of the segmented ILM and RPE as well as of the computed normal RPE position is displayed in blue, red, and green color, respectively. The fundus projection image is shown in Fig. 3; the positions of the B-scans shown in Figs. 3 are indicated by red arrows. A retinal thickness map computed as the axial distance between ILM and segmented actual RPE position is shown in Fig. 3, while for the thickness map in Fig. 3, the axial distance between ILM and the approximated normal RPE position were used. In these maps, single pixels exhibiting thickness values of >150 pixels (corresponding to a retinal thickness of 350 μm) were linearly interpolated between neighboring pixels with values of <350 μm. Finally, a drusen thickness map displaying the axial deviation of the actual RPE position from the computed normal RPE position is shown in Fig. 3. Deviation values of <8 pixels (25.6 μm) were tolerated as small variations of the axial DOPU_min-position within the RPE layer and are displayed gray.

From the drusen thickness map [Fig. 3], drusen dimensions can be assessed quantitatively. Namely, their area and volume can be computed from the known area of a pixel in the en face image [∼12.5(x)×50.2(y) ≈ 630 μm²] and the axial extension of one voxel (∼3.2 μm, optical distance). For the data set shown in Fig. 3, the total area of drusen in the scanned region of the retina was calculated to be 5.527 mm² with a volume of 0.188 mm³. An average refractive index of 1.38 was assumed for drusen volume computation.

Results from applying the segmentation algorithms to the right eye of an 80-year-old female patient diagnosed with dry AMD are shown in Fig. 4. In Fig. 4, an en face fundus image computed from the 3-D data set is shown. The red arrows indicate the position of the reflectivity B-scans of Figs. 4. In Figs. 4, an overlay of depolarizing structures inside and outside of the evaluation band is shown in red and green, respectively. The retinal thickness computed as axial distance between ILM and DOPU_min position is mapped in Fig. 4. In the map shown Fig. 4, the overall thickness of depolarizing tissue [i.e., red and green pixels in Figs. 4] is displayed, while for the thickness map of Fig. 4, only depolarizing pixels within the evaluation band (i.e., the red pixels in the B-scan images) were considered as the depolarizing layer. After processing this image as described in Section 2d, two atrophic zones, sized 0.570 and 3.841 mm², were found automatically. Zones that were connected to the image border [cf. lower left and right corner of Fig. 4] were excluded by the algorithm from the calculation of the total area. The total area of geographic atrophy was 4.412 mm² for this data set.

In order to check the repeatability of the areas and volumes determined for drusen and atrophic zones, respectively, repeated measurements during a single visit were performed in several patients. Five data sets each, from seven eyes diagnosed with drusen and from five eyes diagnosed with geographic atrophy were evaluated. For both the drusen and geographic atrophy data sets, a fixed value of DOPU_thr = 0.8 was chosen for RPE segmentation. For every eye diagnosed with drusen, the drusen thickness maps were registered and cropped to the same retinal section using the OCT fundus image and Image J. Then, the area and volume of drusen were computed and analyzed. In the case of geographic atrophy, only such atrophic zones were integrated into area calculation that did not exceed the OCT device's scan area in all of the five maps. For that reason, one atrophic zone was excluded in one of the investigated eyes.

Results from repeatability measurements can be found in Tables 1- 2 as well as in Figs. 5- 6 for drusen, and in Table 3 and Fig. 7 for the evaluation of geographic atrophy segmentation. Drusen areas in the range of 3.225–9.802 mm² were measured. The average area was 5.735 mm². Measured areas varied by 0.169–0.721 mm² [standard deviation (sd)] with a mean variation of 0.419 mm². The coefficient of variation (sd/mean) was in average 0.077, ranging from 0.032 to 0.126. For the same eyes, drusen volumes ranging from 0.084 to 0.456 mm³ (mean value = 0.211 mm³; mean sd = 0.014 mm³) were observed. Variations between 0.007 and 0.023 mm³ (sd) were computed. The average coefficient of variation was determined to be 0.075 (range: 0.040–0.107).

Segmentation of geographic atrophy yielded areas between 0.620 and 7.788 mm² (mean value = 3.132 mm², mean sd = 0.326 mm²). The variation of determined areas resulted in standard deviations in the range of 0.056–0.979 mm² and in coefficients of variation between 0.036 and 0.140, with an average of 0.092.

The box plots shown in Figs. 5- 7 illustrate the distribution of measured values. The upper and lower edges of the box indicate the 25 and 75% quartile, whereas the whiskers show the minimum and maximum of the measured areas and volumes, respectively. The bar inside the box represents the median.

In the results shown in Section 3, the performance of the introduced segmentation algorithms may be appreciated. Although the method seems to perform well in the eyes shown in Figs. 3- 4, the repeatability measurements show variability values of 7.7 and 7.5% for determined drusen areas and volumes, and a coefficient of variation of 9.2% for areas of geographic atrophy. These variations correspond to 1% (661 pixels) and 0.8% (514 pixels) of the total scanned area for drusen area and geographic atrophy segmentation, respectively. In the case of drusen volume assessment, the COV of 7.5% corresponds to 0.01% (6702 voxels) of the total scanned volume. We believe that these variations are mainly due to eye motion during the measurement. For the PS-OCT instrument used for this study, it takes roughly 3.3 s to record one 3-D data set. Any fixation instability and (transverse) eye motion during that period will inevitably result in distortion of the reconstructed en face images, such as OCT fundus images, drusen thickness maps, or thickness maps of depolarizing tissue. Hence, any quantitative data, such as areas or volumes retrieved from these maps, will be affected. In order to increase repeatability even more, either transverse eye motion has to be compensated for or acquisition speed has to be increased. Transverse motion compensation (especially along the slow scanning Y direction) may be achieved for instance by registering the recorded 3-D OCT data to SLO images from the same fundus region during postprocessing, because SLO images can be recorded much faster than 3-D OCT data sets. A hardware solution for transverse motion compensation would be using a retinal tracker for stabilizing the probe beam with respect to the retina under investigation.³⁷ Recent developments in light source and data-acquisition technology have enabled line rates of several 100,000 A-lines per second, which is more than 10 times the line rate of the PS-OCT instrument we used, which can further reduce motion artifacts.^38- 39 However, even the variation of geographic atrophy area determined from our repeatability measurements [i.e., ∼0.33 mm² (cf. Table 3)] is still significantly smaller than the area growth rate of 1.52 mm²/yr and 4.57 mm²/2yr determined by using FAF and fundus photography reported in literature.^40- 41

One advantage of PS-OCT–based segmentation as compared to—clinically eshablished—imaging modalities, such as fundus photography, FAF, or scanning laser polarimetry, is its ability for depth-resolved imaging. Although all of the aforementioned techniques have proven to be able to image and discriminate pigmented structures, the images will only be the results of a depth integration. Therefore, signals from structures lying on each other will add in the resulting image, which can lead to masking or mimicking of pathology as well as to misinterpretation of retinal features. In contrast, PS-OCT enables discrimination and localization of pigmented structures like the RPE in three dimensions.

It also must be pointed out that, although the area and volume data presented in this paper showed the overall pathology areas and volumes in the respective 3-D data set, the drusen and geographic atrophy maps could also be used to assess the dimensions of single features such as individual druses or atrophic patches.

The results presented in this study were based on two-dimensionally operating segmentation algorithms; that is, segmentation of the RPE was performed B-scan by B-scan. Evaluation in three dimensions could even improve segmentation results because they would enable one to apply smaller evaluation kernels for DOPU computation as well as probably even smoother area and volume assessment in three dimensions. Such further improvements of the segmentation algorithms might also increase data-evaluation speed. Currently, using a state-of-the-art personal computer (Intel Core 2 Quad processor, 4×2.50 GHz, 4 GB RAM, Windows XP) data evaluation from raw data to saving PS-OCT B-scan images and maps takes ∼1.4 and ∼2.3 s per B-scan for segmenting and evaluating drusen and geographic atrophy, respectively (including loading data from a USB hard drive).

In this paper, we presented new applications of PS-OCT for segmentation and quantitative assessment of retinal pathologies in AMD. Area and volume of macular drusen were determined as well as areas of atrophic zones in late stage dry AMD. The segmentation algorithms were based on the polarization scrambling characteristics of the RPE, which allows distinguishing this layer from other retinal layers. The assessed pathology dimensions showed good repeatability and, once again, presented PS-OCT as a versatile method for ophthalmic imaging. PS-OCT might be a promising technique for imaging of retinal pathology in everyday ophthalmology.

The authors thank Martin Wurm and David Stifter at RECENDT GmbH as well as Christoph Wölfl at the Medical University of Vienna for their excellent technical support. Financial support by the Austrian Science Fund (FWF Grant No. UNSPECIFIED P19624-B02) and from the European Union project FUN OCT (FP7 HEALTH, Contract No. UNSPECIFIED 201880) is gratefully acknowledged.