2.1.

Data Collection

In order to accurately measure displacements of the ocular fundus during the cardiac cycle in real time, we made use of a SD-OCT device (Spectralis OCT Plus, Heidelberg Engineering, Germany), the software of which has been custom-modified to record image series in real time, allowing the export of raw images in linear scale. OCT line scans of the same region of the ONH were acquired at 20 Hz for 20 s, during which subjects were asked to look at a fixation light. Since the device is equipped with an active eye tracker that halts acquisition to avoid movement artifacts, subjects are allowed to blink if necessary. Each time series has 401 images ($496\times 768\text{pixels}$) of high resolution ($3.9\mu \mathrm{m}$ axial) acquired over a 15 deg field of view (5 mm approx.) that covers cup, rim, and peripapillary retina, scanned at 45 deg with respect to the fovea to disc axis (Fig. 1). Volunteers covering the glaucoma spectrum, i.e., from glaucoma suspects to advance glaucoma, were recruited at the Maisonneuve-Rosemont Hospital’s Ophthalmology Clinic. Institutional Review Board and Ethics Committee approval was obtained for the study, and all participants gave their informed consent. The research protocol follows the tenets of the Declaration of Helsinki.

Fig. 1

Time series typical images. (a) Fundus image showing the scanned line through the optic nerve head (ONH), in this case at 45 deg. (b) Corresponding optical coherence tomography (OCT) image (two averaged frames).

In order to optimize the parameters of our segmentation pipeline, we presented a set of 25 images that include healthy and several stages of glaucomatous eyes randomly selected from the database, with Spectralis quality scores ranging from 10 to 40, to five observers. Images were segmented varying the values of the main parameters of the algorithm, described below, and raters, masked to parameters and values, were asked to choose the best segmentations for each “unknown” parameter, and the most voted values were selected and kept fixed for all images.

2.2.

Preprocessing

As a first step, a 2-D median filter of size $7\times 7\mathrm{pixels}$ is applied for rough speckle noise removal, then images are submitted to a two-step contrast enhancement technique to suppress features that may lead to erroneous boundary detection. First, intensity exponentiation ($I^2$) is applied to eliminate floaters near the cup and attachment points in case of vitreopapillary detachments. Afterwards, exponentiation $(I^3)+\text{attenuation}$ compensation is required to compensate inhomogeneities due to retinal blood vessels¹⁸ and to sharpen the edges (Fig. 2). In our particular case, since we have previously applied a median filter and intensity exponentiation ($I^2$) to the images, when the exponent value is set to 3, the increment of speckle noise that could be expected according to Ref. 14 is not observed, proving to give the best results for edge detection.

Fig. 2

Preprocessing example. (a) Original ONH image showing floaters in the vitreous. (b) Same image after median filtering and contrast enhancement steps were applied.

2.3.

Inner Limiting Membrane Segmentation

In order to separate the inner layer of the retina from the vitreous, another median filter ($7\times 7\mathrm{pixels}$) is applied. As the image has been previously submitted to contrast enhancement, a threshold of 10% of the maximum intensity is enough to discriminate between the retina and the remaining noise near the edges and is then used to generate the corresponding binary image. An increase of the threshold value results in the loss of the vitreal–retinal boundary information in most cases. The next step is to remove rough edges and bridges by morphological closing with a squared structuring element of 3 pixels width [Fig. 3(a)]. Pixel labeling of eight connected objects is then used to remove the remaining spurious pixels in the vitreous, usually corresponding to highly reflective floaters, by only keeping objects with more than 1000 elements.

Fig. 3

Example of some segmentation steps. (a) Binary image after morphological closing, as the vitreal floaters are not attached to the retina, they are removed after labeling. (b) Final ILM profile, the edge is detected as the first nonzero pixel in each A-scan and interpolated with smoothing splines.

The ILM is set to the most anterior foreground pixel of each A-scan. To exclude misallocated points at the edges of the cup, distances in the axial direction between adjacent contour points $[c(x)]$ are calculated $[\mathrm{\Delta }zi=c(x_i+1)-c(x_i)]$, and we search for distances of $|\mathrm{\Delta }z|$ that are equal to or greater than $27\mu \mathrm{m}$. When a negative targeted value is followed by a positive one and separated by maximum 10 points (transversal direction $x$) then such interval is removed and replaced with a spline-interpolating curve. Figure 3(b) shows the final profile.

2.4.

Deformation Analysis

In order to measure displacements of the neuroretinal tissue due to pulsatility, in addition to segmentation, an analysis algorithm has been developed to track the position of specific retinal regions in time (frame to frame) and to quantify tissue deformation. For alignment, we have used the Matlab intensity-based automatic image registration function with mean-squared error metric and one-plus-evolutionary optimizer. The registration is limited to rigid transformations (translation and rotation), taking as a reference the first image of the series.

Once the segmentation has been performed, the result is a set of 401 ILM profiles per time series. With the purpose of discarding the curves automatically where the ILM delineation has failed, all the profiles are shifted vertically so that their minima coincide [Fig. 4(a)], and the average and standard deviation (S.D.) of the depth ($z$) values at each A-scan position are calculated. Outliers are then tagged in two steps. First, we identify outlier pixels in each one of the profiles, defined as those points located at a distance of 3 S.D. or more from the average value in each A-scan. Next, the S.D. of the number of outlier pixels per curve is determined, and profiles with more outliers than such S.D. are discarded [Fig. 4(b)]. In general, the number of discarded profiles per time series is usually between 0 and 10 but not higher than 15.

Fig. 4

(a) Profiles resulting from segmentation of the whole time series plotted together after vertical shift. Some of the profiles show inaccurate segmentation at the bottom of the cup or the peripapillary retina; hence, two-step outlier identification is carried out. (b) Final profiles after outlier elimination process (10 profiles discarded).

The peripapillary retinal deformation can now be determined from the time series by quantifying changes in the distance measured from the prelaminar tissue to a selected region on the retina. Since the major blood vessels cross the optic nerve through the retina on the nasal region, we have decided to investigate the temporal side, where vessels are less dense, selecting three regions of $200\mu \mathrm{m}$ ($T$, $T2$, $T3$) starting at the Bruch’s membrane opening, as shown in Fig. 5(a). As the bottom of the cup of all profiles is already located at $z=0$, the distance that we want to measure is simply $z(x)$, and deformation is then calculated as the S.D. of such distances averaged within each region [Fig. 5(b)]. It is important to emphasize that although the deformation value obtained per time series could be on the order of the axial resolution of the device, the actual range of measured distances in the time series is between five and ten times larger.

Fig. 5

Deformation measurements: the green line in (a) sets the reference at the deepest point of the ONH, and the red arrow represents the distance measured to determine retinal deformation. Three distance values are calculated per frame, corresponding to the regions T, T2, and T3. (b) Deformation is defined as the standard deviation of the $z(x)$ values of all profiles in each one of the regions (red rectangles).

As a summary, Fig. 6 shows a flow diagram of the segmentation and analysis algorithms proposed.

Fig. 6

Schematic representation of the ILM segmentation and deformation analysis algorithms to measure peripapillary retinal deformation from video-rate OCT image series of the optic nerve.

3.1.

Assessment of the Segmentation Algorithm

To evaluate the performance of the segmentation algorithm, we compared our method with two other established methods: a modified Canny filter¹⁹ (MCF) and a hierarchic approach²⁰ (HA) briefly described below. For each method, the parameters were chosen to give the best results.

MCF: the preprocessing steps include median filtering with a squared mask of 13 pixels width and intensity normalization. Then, the first derivative of a Gaussian function multiplied by a Gaussian is used to create two convolution masks to estimate the gradients in both directions ($\sigma x=0.2$, $\sigma y=4$), and the resulting gradient image is thresholded with $T=0.2*(I_{\max }-I_{\min })+I_{\min }$, with $I_{\max }$ and $I_{\min }$, the maximum and minimum intensity values of the image being processed. Next, linear interpolation is used in a $3\times 3$ neighborhood to find two pixels that straddle the gradient direction; the central pixel is then defined as an edge point if its gradient magnitude is greater than those of the interpolated points. Finally, the ILM is set as the first nonzero pixel per A-scan.

HA: intensity normalization, median filtering with a $5\times 5\mathrm{pixels}$ mask, and image decomposition to a lower resolution with a 16-pixel block using the median constitute the preprocessing. Then, an auxiliary matrix that contains the intensity differences between adjacent depth positions in each A-scan is calculated and employed to define two binary images. The first one contains only the pixels with the highest difference value per A-scan, while a threshold arbitrarily set at 0.03 is applied to the auxiliary matrix to generate the second binary image. These two images are used to identify the borders of the retina, and to correct erroneous recognition, a set of constrictions must be fulfilled. The next steps consist on reducing the decomposition area by half on each iteration to increase accuracy; the positions of the layers determined in the previous iteration are connected with linear interpolation and individual point matching. As final steps, the layers are rescaled, a $5\times 5$ median filter is applied to remove points misallocated in shadowed areas, and interpolation gives the resulting segmented profile.

These three approaches were applied to 25 ONH OCT images randomly chosen from the data set described in Sec. 2.1, with Spectralis quality scores ranging from 11 to 35. In order to validate our algorithm quantitatively, as there is no gold standard for ONH OCT images to compare with, five specialists were asked to delineate the ILM manually for the same set of images. Manual segmentation was performed using a tablet equipped with a stylus. For each image, the average manually segmented ILM was calculated and used to compare the accuracy of every specialist and automated method. Absolute deviation from the mean trace for every A-scan of all images is shown in Fig. 7 as violin plots.

Fig. 7

Quantitative evaluation of the proposed segmentation algorithm. A set of 25 ONH OCT images was manually segmented by five specialists, and the average ILM was used as gold standard for validation. The accuracy of segmentation was evaluated by the absolute deviation from the gold standard of each A-scan for all images. Results for each one of the automated methods are shown as violin plots. The white dot represents the mean value of each distribution.

Our algorithm shows variability comparable to the specialists’, highlighting its robust performance on low signal images. On the other hand, MCF failed to segment the vitreal-retinal boundary due to high speckle content on the images. Notwithstanding HA showed a more robust performance to noise than MCF, its border detection was considerably affected by shadow effects that reduce visibility of some retinal regions. Although MCF and HA show good results for high-quality images ($Q>25$), where a significant contrast at the vitreal-retinal boundary is observed, its accuracy on low signal images is not always satisfactory. Figure 8 shows some examples of the segmented results along with the original images.

Fig. 8

ILM automated segmentation results using different algorithms: a modified Canny filter (MCF), hierarchic approach (HA), and the proposed method. The second column shows the average manual ILM used as gold standard. Images shown have spectralis quality scores of 11, 15, 18, and 30 (from top to bottom).

The time of analysis, together with satisfactory results, is of crucial importance if we take into account that each time series comprises hundreds of images. The average ILM segmentation time for MCF, HA (dividing individual images on blocks of $16\times 16$, $8\times 8$, $4\times 4$), and our algorithm is 9.4 s, 2.4 s, and 0.40 s, respectively, calculated over a set of 100 images ($496\times 768$) analyzed in a PC (3.30 GHz Intel Core i3 processor, 4 GB of RAM).

3.2.

Peripapillary Retinal Deformation

Results from 29 healthy and 12 glaucoma eyes with functional trabeculectomy blebs were obtained using the proposed segmentation and analysis algorithms. In this case, two time series of the same eye were recorded per subject, at 45 deg (inferotemporal region) and 135 deg (superotemporal) with respect to the fovea to disc axis. These particular angles were chosen as they correspond to the areas of major insult in glaucoma. All eyes were transformed to right-eye format.

Recorded time series were excluded due to fixation problems, poor quality images (Spectralis quality $\text{score}<10$ for most of the images), and large blood vessels governing the movement of the neuroretinal tissue. The latter exclusion criterion is vital as we could be measuring directly the pulsation of the vessel instead of the tissue displacement due to choroidal pulsatility, and includes time series where the ONH has a blood vessel located at the base of the cup.

For the sake of simplicity, peripapillary retinal deformation values per eye at each one of the angles correspond to the average value calculated over the three regions defined along the temporal side (T, T2, and T3) on the time series (see Fig. 5). The healthy group shows higher average deformation for both regions, $9.0\pm 3.6\mu \mathrm{m}$ inferotemporal and $8.8\pm 6.6\mu \mathrm{m}$ superotemporal, compared with $4.8\pm 1.8\mu \mathrm{m}$ and $6.9\pm 3.5\mu \mathrm{m}$ in the glaucoma bleb group. However, statistical difference was only found at the inferotemporal region ($P=0.005$) as can be seen in Fig. 9. Group differences were analyzed with the Mann–Whitney test. Mean group deformation values are listed in Table 1.

Fig. 9

Comparison of retinal deformation values between healthy and bleb eyes. Statistical difference was found only at the inferotemporal region ($P=0.005$).

Table 1

Mean peripapillary retinal deformation per group.a

This preliminary study shows that retinal displacement is different for healthy and treated diseased eyes. A more comprehensive characterization could help us to understand how neuroretinal tissue deformation is related to glaucoma progression.