1. Introduction

SDOCT [1, 2] was recently established as a powerful, real-time technique for investigating the depth structure of biomedical tissue with the purpose of non-invasive optical diagnostics. SDOCT is a much faster technique [3] with superior sensitivity [4, 5] as compared to the commercially available time-domain OCT systems and allows for video-rate OCT scans. In ophthalmic applications, it has been suggested that OCT may be helpful for diagnosing glaucoma by measuring the thickness of the RNFL [6–9]. In these reports the RNFL thickness was evaluated with time-domain OCT commercial instruments for only a small number of circular scans, in general three, and not as a full map of the retina. In this paper we present an automated and efficient method for determining the thickness of the RNFL. The method has been applied with the same set of parameters to a large number of retinal scans including the optic nerve head (ONH) and the fovea, and for different measurement sessions. The display modality, also introduced here, combines the RNFL thickness map and a reflectivity map (a fundus-type image), together with the cross-sectional images of the retina (OCT movie), facilitating clinical interpretation of the OCT data. Our analysis has the advantage of providing large area maps of the RNFL thickness that would allow correct registration of the area of interest, which is very important in the longitudinal evaluation of RNFL thinning in glaucoma. In contrast, the currently available individual circular scans need to be measured at precisely the same retinal location, which is very difficult.

Boundary detection has been studied since the early days of computer vision and image processing, and different approaches have been proposed. Segmentation algorithms have also been applied to retinal imaging either for estimating the thickness of various retinal layers [10] or for evaluating the thickness of the retina [11]. The procedure described here is based on a deformable spline (snake) algorithm [12, 13]. As the snake seeks to minimize its overall energy, its shape will converge on the image gradient contour. However, in general, the snake is not allowed to travel too much, and proper initialization is required. The snake parameters (elasticity, rigidity, viscosity, and external force weight) are set to allow the snake to follow the boundary for a large number of retinal topographies. Deformable spline algorithms are widely used in medical imaging. For clarity we briefly describe the principles of the algorithm here. A deformable curve can be thought of as a set of nodes connected together by springs. As the nodes of the contour move, the springs are stretched and flexed, leading to internal restoring forces within the snake which increase the “energy” of the curve. In edge detection with snakes, each contour node is displaced based on a “force” vector derived from the local gradient of the intensity edge map (corresponding to the normalized gradient magnitude of the image). The deformable curve is evolved iteratively until equilibrium is reached between the internal forces arising from the snake stiffness and rigidity and the external forces arising from the image intensity gradient. At this point, the total “energy” due to internal and external forces is minimized. Due to the presence of noise and other artifacts, the minimization process is sensitive to local minima and therefore, contour initialization must be close to the true boundary locations. We use Gaussian blurring of the image when computing the edge map and its gradient in order to reduce sensitivity to local minima and also to increase the capture range of true boundaries within the image.

2. Data processing

The experimental setup used here for video-rate OCT scans has been previously described [3] except that a Ti:Sapphire laser (Integral OCT, Femtolasers, Austria) was used instead of an SLD as light source for the presented measurements. The spectral bandwidth of the laser was about 147 nm. The depth range and axial resolution of the system were 1.6 mm and 2.6 μm in tissue (n = 1.38), respectively. The cross-sectional images consist of 1000 A-lines acquired at a rate of 29 kHz. Two data sets are analyzed in this paper: a) a scan volume of 9.67×6.4×1.6 mm³ that includes both the ONH and the fovea, and b) a scan volume of 6.4×6.4×1.6 mm³ centered on the foveal pit. The pixel size is 9.6×32×1.56 μm³, and 6.4×32×1.56 μm³, respectively. The signal-to-noise ratio (SNR) evaluated as the maximum intensity per frame relative to the noise floor varied for the two data sets from 31 dB to 36 dB, and from 36 dB to 44 dB, respectively. The difference between the two data sets can mainly be attributed to the system’s sensitivity decay with depth since the depth position of the analyzed part of the retina with respect to the zero-delay point changes between measurement sessions.

The depth profile in SDOCT is obtained as the Fourier transform (FFT) of the spectral interference in a Michelson interferometer [1, 2]. The data processing steps to create a good quality structural SDOCT image were described previously [14–16]. To summarize, a DC component was determined by averaging a large number of spectra corresponding to different depth profiles (this averages out the spectral modulation). This DC component is removed from each spectrum recorded by the linear CCD to create a typical, zero-mean interferogram. DC removal eliminates the fixed pattern noise as well as the strong reference arm autocorrelation peak at zero delay that would otherwise create a non-uniform background in the structural image. Interpolation to equally spaced samples in k-space was performed after zero padding by a factor of 4 [17]. Various dispersion compensation schemes for OCT have been investigated [18–21]. We implemented a previously described procedure [14], as follows: the dispersion mismatch between the reference and sample arm can be obtained as the phase of the depth profile for either a strong retinal specular reflection in the center of the fovea in-vivo, or for a mirror in a model eye prior to the in-vivo measurement. The dispersion is compensated by multiplying the spectrum in k-space with the complex conjugate of this dispersion mismatch. In short, the preliminary steps in data processing are: DC subtraction, zero padding, interpolation to k-space, and dispersion compensation. Then, the dispersion compensated spectrum in k-space is Fourier transformed to create the depth profile.

Another important step in data processing is the removal of motion artifacts. During the measurement time, the retina moves due to various reasons, such as involuntary head and eye motion, and cardiac pulsation. These motion artifacts within each frame are significantly reduced in SDOCT as compared to the slower time-domain OCT, however, they are still present from frame to frame in video-rate movies and have to be corrected for proper reconstruction of the retinal topography. In hardware, longitudinal motion (along the OCT beam) can be compensated using adaptive ranging [22], while lateral motion can be balanced using an eye-tracker [23]. Currently, we correct the motion artifacts in software. Two-dimensional cross-correlation of pairs of consecutive frames provides the relative shift of one image with respect to the second image, and we use this information to shift back the second image, both lateral and axial. It is an automated correction and does not require user input.

2.1 Anterior boundary

Data processing for determining the RNFL thickness involves two steps of edge detection. The analysis is performed at this point in 2D, frame by frame, by identifying the edges corresponding to the anterior and posterior boundaries of the RNFL.

The first step is to identify the anterior boundary (AB) of the RNFL. A few preliminary steps are performed in determining an initial guess for AB. The image is blurred using a Gaussian kernel with a standard deviation radius of 4 pixels for smoothing out the edges present in the image while suppressing spurious artifacts due to imaging noise. The edge image is calculated as the magnitude of the image gradient, and rescaled between 0 and 1 by subtracting its minimum value and normalizing to its maximum value. The edge image is then converted to a binary image by keeping all edges above a threshold value of 0.07 (i.e., 7% of the strongest edge). This threshold is determined such that AB remains a continuous line in the edge image. Certain areas in the OCT scan, mainly around the ONH and the fovea, contain weak edges that still need to be identified and this dictates how low the threshold should be set. Some false edges are, however, preserved this way due to the noise in the image. They are eliminated by removing any object in the binary image that has an area smaller than 500 pixels (about 0.07% of the total image size). The purpose is to preserve only continuous lines across the image, and therefore this value can be set based on the size of the image. The result is shown in Fig. 1(a). The initial guess for AB is then determined as the first unity pixel from the top in the binary edge image along each column. To avoid holes in the identified AB, the columns that contain no edge are removed and then AB is linearly interpolated over the removed A-lines. A median filter with a moving window of 50 pixels is applied to smooth out the AB.

 

Fig. 1. Intermediate steps in finding the RNFL thickness. (a) binary edge image obtained from the gradient magnitude of the retinal cross-section; (b) AB - red line, obtained from the multiresolution deformable spline; (c) smoothed field f; (d) edge field s calculated as the rescaled (between 0 and 1) magnitude of the image gradient; (e) binary intensity mask obtained using mean(f) - std(f) as threshold; (f) mean RPE (red dots) estimated between the blue dots and filtered mean RPE (black line) for identifying the blood vessels’ position; (g) dilated and eroded edge field in the interest area; (h) thinned edge field without vertical edges and with only positive edges; (i) initial guess of PB (red dots) and its filtered version; (j) valid A-lines after the deformable spline algorithm; (k) interpolated and median filtered PB; l) retinal cross-section with RNFL boundaries, AB - red line, PB - blue line.

Download Full Size | PDF

This initial guess of AB is used as initialization for a multiresolution deformable spline algorithm. The external force field for the snake algorithm is obtained as the gradient of the edge image that was generated as described above for a Gaussian kernel with a standard deviation radius of 2 pixels. The purpose of the Gaussian blur, as described earlier, is to make sure the initial guess for AB is within the capture range of the true boundary at AB. This also determines the resolution in finding AB. The algorithm is then repeated with a smaller value of the standard deviation of the kernel (Fig. 1(b)), 0.5 pixels, therefore a better resolution, using as initialization the value of AB obtained from the previous run of the algorithm with a coarse resolution.

The anterior boundary of the RNFL can be used to create and display the 3D topography of the retinal surface.

2.2 Posterior boundary

The second step in determining the RNFL thickness is to identify the posterior boundary (PB) of the nerve fiber layer. Due to the fact that the transition between the RNFL and the ganglion cell layer is not as sharp and uniform as between the vitreous humor and the retina, a more elaborate algorithm is required for identifying the PB. Additional challenges are found in this case, such as the shadow cast by the blood vessels from the RNFL that generate “holes” in the PB, and the fact that there is no RNFL in the ONH area (the nerve fibers are along the optic nerve, perpendicular to the retina, and therefore parallel to the OCT beam). If the nerve fibers are perpendicular to the incident laser beam they can be detected as a layer, however, if they are parallel to the OCT beam (as it happens in the optic nerve) they cannot be detected, there is no “horizontal” boundary (perpendicular to the laser beam) to reflect the light.

A few preliminary steps are required before actual estimation of the PB. First, everything above the AB is removed and the image is realigned to AB based on the assumption that the anterior and posterior boundaries are relatively parallel and the algorithm in determining PB is more stable for predominantly horizontal edges. A depth of 350 points, corresponding to a layer thickness of 552 μm, is kept along each A-line. This depth is selected to include the relevant part of the depth profiles. All images in Fig. 1 from (c) to (k) have AB as the top edge, and represent the steps in determining PB.

Second, of major importance is the processing of the image to generate the smoothed field f (Fig. 1(c)) and the edge field s (Fig. 1(d)) by using an algorithm for joint anisotropic smoothing and edge-preservation. The edge field s is the rescaled (between 0 and 1) magnitude of the gradient of f, as described above. The method is based on a conjugate gradient based iterative minimization [24] of a cost functional which balances a data fidelity term against prior models for image smoothness and edge density [25, 26]. The cost functional is given by

where the notation ∥.∥_k represents an l_k norm (k = 1,2), and the integration is done over the entire image. The first term represents the data fidelity and controls the degree to which the smoothed field f resembles the original image g. The second term, the smoothness constraint, penalizes large gradients in f except where edges exist (s~1), generating the greatest anisotropic smoothing far away from the edges. The last two terms, the edge penalty, control the edge width and prevent the minimization process from placing edges everywhere in the edge field s. The real positive scalars α (75), β (0.03), and ρ (0.01) adjust the relative weighting between these competing terms. The solution of Eq. (1) is obtained by iterative minimization of the cost functional until convergence to a predefined tolerance level is achieved or the maximum number of iterations (50) is exceeded. All the subsequent steps in the image processing are based on f and s, and not the original image.

Identification of blood vessels position in the RNFL is also required before the estimation of the PB. This is done based on the analysis of the retinal pigment epithelium (RPE). Scattering and absorption of the OCT beam due to the blood significantly reduce the ability of the OCT beam to probe behind the blood vessels and the structural image appears to have a “shadow” underneath the blood vessels (Fig. 1). This creates “holes” in the PB as well as in the RPE. We need first to identify the RPE in the image and estimate the index of the A-lines corresponding to these “holes”. These indices, declared as invalid, will later be used in the analysis of the PB.

The identification of the RPE is based on both the magnitude of the smoothed field f and its gradient, the edge field s. An intensity mask (Fig. 1(e)) is generated from the smoothed field f based on the statistical properties of f, mean and standard deviation, using mean(f) - std(f) as threshold. A binary image is also created by thresholding the edge field s. The threshold value is selected as described in the previous section (0.05 in this case), to preserve even weaker edges. Small patches with area smaller than 50 points are removed from the binary edge image to clear out spurious points. The lowest (deepest) pixel with value one along each A-line is selected from the two binary images as the posterior boundary of RPE. This boundary is then smoothed out with a 55 points moving window median filter. A band with a predefined thickness (25 points) above the posterior RPE boundary is then selected from the smoothed field f (delimited by the blue dots in Fig. 1(f)) and is averaged along the A-lines to generate a mean value of the RPE reflectance along each A-line (red dots in Fig. 1(f)). The mean RPE reflectance is then filtered out using a Savitzky-Golay FIR smoothing filter (black line in Fig. 1(f)) that fits the data with a third order polynomial in a 201 pixels moving window. The index of an A-line is declared invalid if the mean RPE value drops on that A-line by more than a certain predefined value (5 in reflectance units) as compared to the filtered mean RPE. As mentioned, these invalid indices correspond to the “holes” in the RPE generated by the shadow of the blood vessels. They also correspond to the ONH area where there is no RPE. The invalid indices are extended by a number of 10 A-lines on both sides of the shadows, as identified above, just to make sure that the vertical edges corresponding to the boundary of the shadows are removed from the image. Generally, the blood vessels have a larger diameter than the holes identified in the RPE.

The RPE area is removed from the smoothed field f, the intensity mask, and the binary edge image, and the rest of the processing is focused on the posterior boundary of the RNFL.

The intensity mask is applied to the edge field in order to remove the edges that are outside the interest area. To avoid broken lines, the binary edge image is dilated and then eroded (Fig. 1(g)). A morphological thinning operation is then performed on the binary edge image. Vertical and negative edges are removed since we are looking for predominantly horizontal edges with the correct slope from the RNFL to the ganglion cell layer (Fig. 1(h)). An initial guess of the PB is estimated as the first pixel of value one from the top along each A-line (red dots in Fig. 1(i)). For the A-lines that have no pixel of value 1 in the binary edge image, the PB is set to coincide with AB. This situation corresponds to the areas around the ONH and fovea. The columns with invalid indices are removed and the PB is linearly interpolated over the invalid regions. To make sure that the PB is relatively smooth, the smoothing procedure mentioned above for the RPE analysis employing a Savitzky-Golay FIR smoothing filter is used here as well with a moving window of 51 points (black line in Fig. 1(i)). This step is intended to remove points that are too far from the smoothed version of PB.

At this point we are ready to apply the deformable spline algorithm. The intensity mask is applied to the original edge field s, and the edge field is then blurred with a Gaussian kernel of 2 pixels standard deviation radius. The external forces are calculated as gradient of the rescaled edge field, and they are set to zero for the A-lines with invalid indices. The A-lines with invalid indices are removed from the result of the snake algorithm (Fig. 1(j)) and PB is linearly interpolated over the invalid regions. A final median filter with a moving window of 11 points is applied to the calculated PB (Fig. 1(k)).

Since the image was realigned to the anterior boundary, the result obtained here is actually the RNFL thickness. Coming back to the original image, one has to add this result to the AB values to obtain the true PB (Fig. 1(l)).

The results of the algorithm described above are shown in Fig. 2 by displaying the identified AB and PB on the corresponding structural image. The video-rate OCT movie has been corrected in software for motion artifacts.

 

Fig. 2. Movie of the OCT scan showing the anterior (red) and posterior (blue) boundary of the RNFL a) including the ONH and the fovea (1.99 MB), and b) centered on the foveal pit (2.39 MB). The movies consist of 170 (a) and 180 (b) frames, respectively, displayed in a reversed-contrast logarithmic gray-scale at 30 fps. Each frame has a size of 8.85×1.2 mm² (a) and 6.07×0.786 mm² (b), respectively. The vertical size was increased by a factor of 4.65. (14.4 MB version (a) and 14.5 MB version (b))

Download Full Size | PDF

All the parameters in the program, such as threshold values, snake parameters for anterior and posterior boundaries, and RPE thickness, are set based on a large number of OCT scans in different areas of the retina to account for a statistically significant variability in the boundaries’ characteristics.

The snake parameters need to be set different for AB and PB given the different properties of the two boundaries. All data sets presented here were processed with fixed settings for AB and PB. The processing of a single image took 62 seconds on a 3.2GHz Pentium 4 processor.

3. Results

Generally, a video-rate OCT movie is hard to follow and interpret. It is not always clear what happens with the blood vessels moving laterally in the image and with the continuous and sometimes sudden change in their apparent diameter. If the orientation of a blood vessel changes from frame to frame, it appears in the OCT-movie that the vessel’s diameter changes. There are situations in the OCT-movie when the vessel is parallel to the frame and it “appears” very thick, and then in the next frame is very thin since it changed its orientation with respect to the frame. Additional information is required for correct interpretation of OCT-movies. Recently, a visualization was demonstrated where a fundus-type image was shown as a still picture simultaneously with the OCT movie [27]. A line across this picture indicated the position of the OCT scan for orientation along the retinal surface. This picture was obtained by integrating the depth profiles (the reflectance) and was displayed on a logarithmic scale. This operation creates enough contrast to differentiate between the A-lines that correspond to blood vessels’ location and those that do not intersect blood vessels. Clearly, the scattering and absorption on blood reduce the total reflectance in the area corresponding to blood vessels and the vessels appear darker than the surroundings. We obtained a smoother image by integrating the logarithmic depth profile and displaying it in linear scale as shown in the right image of Fig. 3.

For the dataset corresponding to the results shown in Fig. 2(a), Fig. 3, Fig. 5, and Fig. 6, the size of the acquired scan area is 9.6×6.4 mm² and includes both the ONH and the foveal pit. The images were cropped to 922 A-lines to eliminate sides’ artifacts and cover a depth of 1.2 mm (761depth points), while the retinal maps are made of 170 scans representing a scan area of 8.85×5.73 mm² acquired in 5.9 seconds. For the results shown in Fig. 2(b), Fig. 4 and Fig. 7, the size of the acquired scan area is 6.4×6.4 mm² and is centered on the foveal pit. The images were cropped to 949 A-lines by 500 depth points for 180 scans representing a scan volume of the size 6.07×5.76×0.786 mm³ acquired in 6.2 seconds.

The integrated reflectance image is compared in Fig. 3 and Fig. 4 to images of the same area in the same eye obtained with different methodologies, fundus imaging and fluorescein angiography, respectively.

 

Fig. 3. Fundus image (left), and integrated reflectance image (right) for the same eye. The size of the image is 8.85×5.73 mm².

Download Full Size | PDF

As compared to a fundus image of the same eye (Fig. 3), the integrated reflectance image shows the same structure of the blood vessels with a very good quality, approaching that of a fundus photo.

 

Fig. 4. Angiogram (left) and integrated reflectance image (right) for the same eye. The size of the image is 6.07×5.76 mm².

Download Full Size | PDF

Figure 4 shows an angiogram (left) and the integrated reflectance image (right) for the same area around the foveal pit. Both Fig. 3 and Fig. 4 demonstrate that the integrated reflectance map can be used as a reliable representation of the retinal vasculature. The data sets used for generating these integrated reflectance maps have been corrected in software for motion artifacts (as described in Section 2), however, slight discontinuities in the blood vessels can still be noticed at the top part of both images. These discontinuities could be due to a temporary loss of fixation.

3.1 RNFL thickness map

The experimental measurements presented here were performed on the right eye of a healthy volunteer. Figure 5 shows a map of the RNFL thickness as obtained following the procedure described in the previous Section. The RNFL thickness map was smoothed out with a 10×10 pixels median filter corresponding to an area of 96×320μm². On the color map, the darkest blue areas represent the absence of the RNFL, corresponding to the ONH and fovea. Lighter blue corresponds to thinner areas of the RNFL, while darker red represents thicker RNFL areas. The maximum RNFL thickness here is 177 μm. This RNFL map is consistent with known normal retinal anatomy, as the RNFL thickness is greater superior and inferior to the ONH. The classic bow tie pattern (in red) is seen (Fig. 5). This map is a quantitative assessment of the RNFL thickness and provides a comprehensive evaluation of large retinal areas as opposed to a limited number of circular scans measured with the current commercial instruments.

 

Fig. 5. RNFL thickness map. The dark blue areas correspond to the ONH (right) and fovea (left), and indicate no RNFL thickness. The dark red areas indicate a maximum thickness of 177 μm. The color bar is scaled in microns.

Download Full Size | PDF

Figure 6 displays the combined retinal information analyzed here: the top left picture shows the integrated reflectance as a map of the retina, while the bottom left image shows the RNFL thickness map. The right side movie is a video-rate OCT scan corrected for motion artifacts. Using this display modality, one can follow much easier in the movie the structure of the blood vessels pattern. The position of each depth scan, indicated by the horizontal red line across the maps, can be related more explicit to retinal morphology, i.e. the ONH and the fovea. The position of the blood vessels across the cross-sectional images, indicated by their “shadow”, can easily be correlated with the intersection of the horizontal line with the vasculature evident in the integrated reflectance map. The integrated reflectance map illustrates also the orientation of the blood vessels with respect to the cross-sectional scans, allowing for a clear interpretation of the continuous and sometimes sudden change in the apparent vessels’ diameter. The association of the integrated reflectance map and of the RNFL thickness map with the OCT movie gives clinicians a more intuitive means of interpreting the OCT data for diagnosing retinal diseases such as glaucoma.

 

Fig. 6. Movie (1.91 MB) with combined integrated reflectance map (top left), RNFL thickness map (bottom left), and retinal cross-sectional images corrected for motion artifacts (11 MB version). The color scheme for the RNFL thickness map is scaled in microns, dark blue meaning no thickness, and dark red being a maximum of 177 μm, as shown in Fig. 5. The movie consists of 170 frames displayed at 30 fps. Each map has a size of 8.85×5.73 mm². The vertical size of the cross-sectional images (1.2 mm) was increased by a factor of 4.65.

Download Full Size | PDF

A similar OCT scan centered on the foveal pit is shown in Fig. 7 together with the corresponding integrated reflectance map and RNFL thickness map. It shows the vasculature around the fovea and the retinal depth structure with a greater detail as compared to Fig. 6.

 

Fig. 7. Movie (1.33 MB) with combined integrated reflectance map (top left), RNFL thickness map (bottom left), and retinal cross-sectional images corrected for motion artifacts (7.33 MB version). The color map for the RNFL thickness map is scaled in microns, dark blue meaning no thickness, and dark red being 105 μm. The movie consists of 180 frames displayed at 30 fps. Each map has a size of 6.07×5.76 mm². The vertical size of the cross-sectional images (0.786 mm) was increased by a factor of 4.65.

Download Full Size | PDF

The dark blue band on the center left side of the RNFL thickness map shown in Fig. 7 corresponds to the temporal raphe [28], a structure located temporal to the macula. Since our system cannot distinguish individual fibers, the structure and the direction of the retinal nerve fibers cannot be seen in the thickness map. However, the RNFL thickness is small in the raphe area since there are not too many fibers. Going away from the raphe, more and more fibers comprise the RNFL and the thickness increases. Our RNFL thickness map is consistent with the RNFL distribution pattern in Vrabec’s report [28] that was confirmed later on by ophthalmoscopy [29].

4. Conclusions

To summarize, we have presented an automated and efficient method for determining a map of the RNFL thickness, based on anisotropic noise suppression and edge detection using deformable splines. The method has been applied with the same set of parameters to a large number of retinal scans including the ONH and the fovea, and for different measurement sessions on a single subject. The algorithm performed well in 350 frames within an SNR range of 31–44 dB, where the boundary was incorrectly identified in only a few and isolated places. The performance of the automated algorithm in patient measurements could be adversely affect by ocular opacities or diseases that reduce the SNR, or by a reduced contrast between RNFL and ganglion/inner plexiform layer. The RNFL thickness map is a quantitative assessment and provides evaluation of large retinal areas as compared to a limited number of circular scans measured with the current commercial instruments. The data was presented in a correlated manner that includes the RNFL thickness map and an integrated reflectance image together with an ultra-high resolution OCT movie. This representation provides a comprehensive picture to clinicians and we anticipate that it will be a valuable tool for interpreting the OCT data for diagnostic purposes. The RNFL thickness maps can potentially be used for a thorough evaluation of the RNFL thickness in longitudinal studies of glaucoma progression. These studies require large area RNFL thickness maps, which may allow for more accurate correlations of RNFL thinning with visual field defects as opposed to individual circular scans that need to be measured at precisely the same retinal location and that give less information.