A review of algorithms for segmentation of optical coherence tomography from retina.

Optical coherence tomography (OCT) is a recently established imaging technique to describe different information about the internal structures of an object and to image various aspects of biological tissues. OCT image segmentation is mostly introduced on retinal OCT to localize the intra-retinal boundaries. Here, we review some of the important image segmentation methods for processing retinal OCT images. We may classify the OCT segmentation approaches into five distinct groups according to the image domain subjected to the segmentation algorithm. Current researches in OCT segmentation are mostly based on improving the accuracy and precision, and on reducing the required processing time. There is no doubt that current 3-D imaging modalities are now moving the research projects toward volume segmentation along with 3-D rendering and visualization. It is also important to develop robust methods capable of dealing with pathologic cases in OCT imaging.

INTRODUCTION

Optical coherence tomography (OCT) is a recently established imaging technique to describe different information about the internal structures of an object and to image various aspects of biological tissues, such as structural information, blood flow, elastic parameters, change of polarization states, and molecular content.[1] In contrast to OCT technology development which has been a field of active research since 1991, OCT image segmentation has only been more fully explored during the last decade. Segmentation, however, remains one of the most difficult and at the same time most commonly required steps in OCT image analysis. No typical segmentation method exists that can be expected to work equally well for all tasks.[2]

One of the most challenging problems in OCT image segmentation is designing a system to work properly in clinical applications. There is no doubt that algorithms and research projects work on a limited number of images with some determinate abnormalities (or even on normal subjects) and such limitations make them more appropriate for bench and not for the bedside. Moreover, OCT images are inherently noisy, thus often requiring the utilization of 3D contextual information. Furthermore, the structure of the retina can drastically change during disease. Nevertheless, OCT image segmentation is a rapidly growing and important area and a great deal of effortswent into designing algorithms for automatic segmentation of retinal OCTs.

The important steps of OCT image segmentation algorithms may be categorized in to 4 steps: 1) Determining a particular range of OCT datasets (e.g., 2D, 3D, Time Domain, Spectral Domain, macular, ONH, etc.) for which the proposed algorithm will work properly; 2) Allocating proper values for possible parameters of the algorithm; 3) Running the algorithm on determined datasets and acquiring the outcomes (e.g., boundary information, thickness values, classification of normal and abnormal images and etc.); 4) Validating the results by comparing with gold standards or getting the confirmation by an expert.

Here, we review some of the important image segmentation methods for processing retinal OCT images. We may classify the OCT segmentation approaches into five distinct groups according to the image domain subjected to the segmentation algorithm. We define five separate families of segmentation approaches: Methods applicable to A-scan, B-scan, active contour approaches (frequently in 2D), analysis methods utilizing artificial intelligence, and segmentation methods using 3D graphs constructed from the 3D OCT volumetric images. Some details of algorithms representing each class will be discussed in more detail and concluding remarks will be provided to compare the efficiency of different methods in different datasets. It should be noted that intrinsic noisy structure of retinal OCTs (particularly in old OCTs and after diseases) makes simple edge detection algorithms unsuitable for this purpose and researchers have tried to develop new algorithms to overcome such problems.

OTHER RETINAL IMAGING MODALITIES

It should be noted that OCT is not the only possible device for assessing retinal pathologies; the field of ophthalmology was revolutionized in 1851 with the invention of the ophthalmoscope by Hermann von Helmholtz[3] as for the first time detailed examinations of the interior of the eye could be made in living patients. Fundus photography (a low powered microscope attached with a camera),[4] fluorescein Angiography[5] (photographing the retina by injecting fluorescent dyes) and Retinal thickness analyzer (RTA)[6] are other modalities proposed for diagnosis of retinal malfunctions. The latter is capable of rapidly obtaining retinal thickness map covering an area of 3 × 3 mm. The oblique projection of a narrow laser slit beam on retina and recording the backscattered light are the principles of this method.[6] Furthermore, confocal scanning laser ophthalmoscopy (CSLO)[7] provides a three-dimensional topographic representationof the optic disk and peripapillary retina, which is constructed from a series of two-dimensional slices. This three-dimensional representation consists of 256 × 256 (65,536) pixel elements, each of which is a measurement of retinal height at its corresponding location. Three topography images are usually acquired in a single session and thereafter are automatically aligned and averaged to obtain a single mean topography image. Although the CSLO is similar, in many respects, to a CT scan, the light rays used for CSLO cannot penetrate tissue; this limits this modality to depicting the surface topography of the optic disk and para-papillary retina.[8]

PRINCIPLES OF OPTICAL COHERENCE TOMOGRAPHY

Optical coherence tomography utilizes near infrared super luminescent diode light in a fashion similar to the way B-mode ultrasound uses sound to generate two-dimensional images. Although sound has the advantage of penetrating tissue (for example, scanning a fetus in utero), light, with its much shorter wavelength, has the advantage of obtaining significantly higher resolution. A simplified optical setup for a fiber based OCT system utilizing a low coherence source and a Michelson-type interferometer is illustrated in Figure 1. Cross sectional data along an axial line through the sample, called an A-scan, is acquired by axially scanning the position of the reference arm. Interference fringes are acquired at the photodiode detector when the round trip distance from the sample reflection matches that of the reference reflection. The magnitude of the observed fringes is proportional to the reflectivity of the scatterer. A two dimensional profile, or B-scan, is generated by scanning the interrogating beam laterally across the sample and acquiring an axial scan at each lateral location. Subtle differences in adjacent layers are visualized as differences in scattering intensities.[9101112]

Figure 1

Block-diagram of optical coherence tomography

The scan speed of the time-domain OCT scanner such as a Stratus OCT (Carl Zeiss Meditec, Inc., Dublin, CA, USA) machine is 400 A-scans/sec. The scan speed is limited by the speed at which the reference mirror can be moved. Therefore, the time-domain OCT scanner can acquire only several cross-sectional images (typically 6 × 128 × 1024 voxels). The slow scan speed causes a bad effect on image quality such as eye movement artifacts [Figure 2a]. The voxel size is 30 × 30 × 2 μm, the voxel depth is 8 bits in gray scale, and the acquisition time for each scan consisting of 6 × 128 × 1024 voxels is 1.92 seconds.[9]

Figure 2

Image quality of time-domain and spectral-domain OCTs (a) Time domain OCT image from Stratus OCT (Carl Zeiss Meditec, Inc., Dublin, CA, USA) (b) Spectral-domain OCT image from Cirrus HD-OCT (Carl Zeiss Meditec, Inc., Dublin, CA, USA)[9]

In time domain OCT the location of scatters in the sample is observed by generation of interferometric fringes at the detector as the reference reflector position is axially translated. In contrast, Fourier domain OCT required the reference arm to be held fixed, and the optical path length difference between sample and reference reflections is encoded by the frequency of the interferometric fringes as a function of the source spectrum. Two configurations have prevailed in Fourier domain systems: spectral domain (SD) OCT uses a grating (network of parallel wires) to spatially disperse the spectrum across an array-type detector, and in swept source (SS) OCT a narrow band laser is swept across a broad spectrum, encoding the spectrum as a function of time [Figure 3].[10111213]

Figure 3

Fourier domain configurations. Spectral domain OCT utilizes a spectrometer in the detection arm, whereas swept source OCT requires a wavelength swept laser and a standard photodiode detector

Regardless of whether the spectrum is sampled in time or in position across an array detector, the frequency of the interferometric fringes as a function of spectrum encodes the location of the scatter, with increasing frequency corresponding to larger optical path length mismatches. A common mathematical tool for extracting the frequency content of a signal is the Fourier transform, however, it must be remembered that the true Fourier transform pair of distance is not wavelength (also units of distance), but in fact wavenumber (with units of inverse distance).[14]

The scan speed of the spectral-domain OCT scanner such as a Cirrus HD-OCT (Carl Zeiss Meditec, Inc., Dublin, CA, USA) machine is 27,000 A-scans/sec. The Cirrus HD-OCT scanner can acquire a retinal volume whose dimension is typically 200 × 200 × 1024 voxels covering 6 × 6 × 2 mm3 [Figure 4]. The voxel size is 30 × 30 × 2 μm, the voxel depth is 8 bits in grayscale, and the acquisition time for each volumetric scan consisting of 200 × 200 × 1024 voxels is 1.48 seconds. The image quality of the spectral-domain OCT is much better than that of the time-domain OCT [Figure 2b].[9]

Figure 4

Example dimensions of a 3-D spectral OCT image. A typical spectral OCT images covers a volume of 6×6×2 mm3 (200×200×1024 cubic voxels) using one of the protocols on the Cirrus machine by Carl Zeiss Meditec[15]

The spectral-domain OCT scanner can provide true 3-D views of the retinal structure. Two kinds of OCT scans, the scan centered at the fovea (macular scan) related to central vision and the scan centered at the optic nerve head (ONH-centered scan) related to peripheral vision, are mainly acquired to examine patients’ eyes [Figure 5].[9]

Figure 5

Spectral-domain volumetric OCT scans of the retina (a) Color fundus photograph of the retina. The green square shows the position of the scan centered at the fovea (macular scan), and the blue square represents the position of the scan centered at the ONH (ONH-centered scan), (b) X-Y image of the macular scan, (c) X-Z image of the macular scan corresponding to the green line in image (b,) (d) Y-Z image of the macular scan corresponding to the red line in image (b), (e) X-Y image of the ONH-centered scan (f) X-Z image of the ONH-centered scan corresponding to the green line in image (e), (g) Y-Z image of the ONH-centered scan corresponding to the red line in image (e)[9]

THE LAYERS OF THE RETINA

Covering the inside of most of the eye [Figure 6], the retina is a multilayered structure responsible for transforming light energy into neural signals for further use by the brain. In very general terms, the processing of light starts with the light sensitive photoreceptor cells (rods and cones), which are actually located in the outer portion of the retina (away from the incoming light).[15] These cells convert the light signal into action potentials that are transmitted by the bipolar neurons in the central layers of the retina to the ganglion cells of the inner retina. It is the axons of the ganglion cells that eventually exit the eye to form the optic nerve.[15] Other cells in the retina, such as horizontal cells, amacrine cells and interplexiform neurons, also help in the processing of the neural signal at a local level. Neuroglial cells (such as Muller cells) provide structure and support.[16] Many of the cells of the retina are illustrated in Figure 6a. Based on its appearance from light microscopy [Figure 6b], the retina is traditionallyconsidered to be composed of the following ten major “layers” (starting with the outermost layer):[16]

Figure 6

Different views of the retina (a) Schematic illustration of cellular layers of retina, (b) Light micrograph of a vertical scan through central human retina, (c) Cross section of the eye with illustration of the retina, (d) OCT view of macular retina[15]

Reginal pigment epithelium (RPE): single layer of pigmented hexagonal cells

Photoreceptor layer: the outer (containing the light-sensitive discs) and inner segments of rods and cones

External (or outer) limiting membrane (ELM or OLM): intercellular junctions between photoreceptor cells and between photoreceptor and Muller cells (not an actual membrane)

Outer nuclear layer (ONL): rod and cone cell bodies

Outer plexiform layer (OPL): synapses between photoreceptor cells and cells from the inner nuclear layer

Inner nuclear layer (INL): cell bodies of bipolar cells, horizontal cells, amacrine cells, interplexiform neurons, Muller cells, and some displaced ganglion cells

Inner plexiform layer (IPL): synaptic connections between bipolar cell axons and ganglion cell dendrites

Ganglion cell layer (GCL): mostly ganglion cell bodies

Nerve fiber layer (NFL): ganglion cell axons

Internal limiting membrane (ILM): innermost membrane of retina separating the retina from the vitreous.

Intraretinal layers are also visible from optical coherence tomography images as shown in the example image of [Figure 6c and d]. Note that the precise anatomical correspondence of the layers visible in OCT images is not known. The anatomical labels are our current presumption based on comparisons with histology and images from higher resolution OCT scanners published in the literature.[17]

CLINICAL APPLICATIONS OF OCT SEGMENTATION

The most motivating point in clinical applications of OCT segmentation is providing a quantitative tool to help the ophthalmologists to manage the high complexity of the OCT data and to create a tool for better observation of different boundaries and individual layers. There is no doubt that in new modalities of OCT, the 3D data is overly abundant which makes the interpretation really difficult for the observer; but, a correct segmentation can provide useful information about the layers and makes the ophthalmologist able to perform better diagnosis and treatment. A variety of applications may be defined with this approach which can be categorized as:

Earlier detection of ganglion cell loss in cases of concurrent optic nerve swelling can be possible. In cases of optic disc edema, it is possible that information in the macular scans, such as the thickness of the ganglion cell layer from an intraretinal segmentation, be able to detect such loss at an earlier stage[15]

The thinning of the RNFL is an important marker used by ophthalmologists in the diagnosis of the disease. It has also been used to differentiate between glaucomatous and normal eyes[18192021]

Volumetric measurements in the whole scan, within selected sectors/quadrants and radial distances can be useful in detection of retinal pathologies

The comparison of retinal thickness during the treatment can become more simple and quantitative using intraretinal segmentation

Calculation of similarities between right and left eye can be assisted using intraretinal segmentation and such similarities seem to be a promising sign of healthiness in individuals

Different pathologies may make particular layers thinner (or even eliminates them), which can be found through layer segmentation approaches.

REVIEW OF ALGORITHMS FOR OCT IMAGE PREPROCESSING

Preprocessing is the first and one of the most important parts of an image processing system. OCT images suffer from the intrinsic speckle noise, which decreases the image quality and complicates the image analysis. This particular noise is the foundation of existing problems in the precise identification of the boundaries of the various cellular layers of the retina and other specific retinal features present in the OCT tomograms.[2] Usually all OCT image analysis methods proposed in the literature consist of a preprocessing step before performing any main processing steps. Table 1 shows a relatively complete classification of denoising algorithms employed in OCT segmentation. As it can be seen in this table, median filter and non-linear anisotropic filter are the most popular methods in OCT image denoising. The basic problem associated with most of the denoising algorithms is their intrinsic consequencein decreasing the image resolution. The popularity of methods like non-linear anisotropic filters and wavelet diffusion can be justified through their ability in preserving the edge information. It is also important to know that a great number of newly developed algorithms which utilize graph based algorithms are literally independent from noise and do not use any particular denoising algorithm.[39404142]

Table 1

Classification of preprocessing algorithms employed in researches on OCT segmentation

REVIEW OF ALGORITHMS FOR OCT IMAGE SEGMENTATION

OCT image segmentation is expected to help the ophthalmologists in diagnosis and treatment of retinal pathologies; however, it should be mentioned that structural changes due to diseases usually changes the OCT appearance, considerably. Therefore, a low portion of available researches have focused on structurally elaboratedillnesses and most of papers deal with pathologies that despite making changes in thickness or placement of layers, would retain the logical organization of retinal layers. The other common problem in image segmentation is unavoidable speckle noise consideration which was thoroughly reviewed in previous section. The only point to be discussed in noise management is that new versions of OCT imaging like Spectralis Heidelberg are equipped with high speed of imaging which makes the capable of obtaining many pictures in each instant and averaging them to reduce the unwanted speckle noise. The intensity fluctuation is the third tip to be considered in OCT layer separation. Different researchers have diverse ideas about this problem; some investigations show that this fluctuation is only due to noise of the imaging system, but a great deal of papers assume it as a result of absorption or scattering in retinal layers, where the intensity diminishes with increasing the depth of imaging in retina. The problem of blood vessels should also be considered in OCT images, which makes discontinuities in boundaries of different layers; The robustness of the algorithm in presence of blood vessel artifacts is discussed in many papers,[41424344] and some researchers tried to propose an algorithm in preprocessing step to compensate for the effect of these vessels.[44] Furthermore, motion artifacts can decrease the quality of images and mutilate the operation of segmentation method; such problems are recently solved in new versions of OCT imaging systems like Spectralis Heidelberg which has an eye tracking systems to automatically compensate the eye movements.

In this section we focus on different segmentation methods, introduced in OCT image segmentation researches. We may classify the OCT segmentation approaches into five distinct groups according to the image domain subjected to the segmentation algorithm. We define five separate families of segmentation approaches: Methods applicable to A-scan, B-scan, active contour approaches (frequently in 2-D), analysis methods utilizing artificial intelligence, and segmentation methods using 3D graphs constructed from the 3D OCT volumetric images.

A-Scan methods were firstly introduced by Hee[10] and were popular until 2005.[2223] The method introduced in Hee[10] was on the basis of variations in intensity and measured retinal and RNFL thickness. Instead of using simple threes holding, this method was based on one-dimensional edge detection in each A-Scan and was looking for the two most effective edges using peak detection.

Huang[1] used a similar method in hereditary retinal degenerations in experimental animals and humans.

In 2000, George[22] used the similar dual threshold to segment and choriocapillaries from OCT images without taking into account any spatial characteristics of OCT images.

In designing a retinal boundary detector, Koozekanani[23] used retinal anatomy and the principles of OCT operation to make various assumptions about the image boundary characteristics. For instance, the normal retina has smooth boundaries without discontinuities or gaps and the inner boundary is always above the outer boundary. Because the OCT acquires each A-scan separately, they applied one-dimensional edge detection to each image column (A-scan) individually, to overcome the multiplicative speckle noise in 2-D images which could make edge detectors problematic. Each A-scan penetrates both retinal boundaries and so they assumed that every image column intersects exactly two boundaries. They marked the positively sloped, leading edges of peaks, as they are more consistent and easier to detect. They modeled the boundary displacement between adjacent scans as a Mth order Markov sequence. Boundary deviations from one A-scan to the next were due to the actual slope of the retinal surface, patient motion relative to the OCT machine and image noise. They represented the boundary of interest (inner or outer) as a weighted sum of the neighboring A-scans, making the boundary model Mth order. Thus, they used an autoregressive mathematical model for this Markov process. Retinal thickness was calculated with an error comparable to the 10μm resolution of the OCT system used, representing a substantial improvement over clinical measurements provided by the Humphrey 2000 OCT built-in algorithm. Since this model relied on simply connecting 1D points, it was sensitive to noise. Thus, it demanded to apply special rules to correct the errors in the extracted layer borders.

Herzog[24] proposed a method for retinal layer segmentation from axial TDOCTs through the optic nerve head (ONH). The basis of the method was looking for curves with maximum number of located points on boundaries and with minimum rate of change. An example of retinal choroid boundary identification in Herzog's[24] method is shown in Figure 7a and the best fit model used for extracting the cup limits is demonstrated in Figure 7b. Given the retinal-vitreal boundary in Herzog,[24] the retina and optic cup can be modeled through the parametric set of equations. The left and right retinal-optic disk regions were modeled using straight line segments while the optic cup was modeled by a parabolic segment. The goal was to find where the breakpoints (the extents of the optic cup) were located. If two points were picked arbitrarily for the breakpoints, the model parameters in Equation could be found via least-squares fit.

Figure 7

(a) An example of retinalchoroid boundary identification in Herzog's[24] (b) shows the best fit model used for extracting the cup limits

Shahidi[25] and Ishikawa[20] proposed similar algorithms based on intensity variations. It is worthy of mentioning that the work of Ishikawa[20] was the first report demonstrating that the thickness of the innermost layers in the macula had diagnostic power comparable with that of circum papillary nerve fiber layer (cpNFL) in glaucoma studies[2] and the number of segmented retinal layers in this study were four.

A new approach in preprocessing was also introduced by Ishikawa:[20]

Aligned z-offset (starting location of the meaningful signal on each sampling line) by cross-correlation (shifting sampling lines so that the sum of the products of adjacent pixels is maximized; Figure 8a).

Figure 8

Preprocessing of an OCT image by Ishikawa.[20] (a) A raw OCT macular scan image was aligned by cross-correlation, (b) A modified mean filter was applied to the aligned image above

Equalized the histogram of pixel intensity on each line by scaling the pixel intensities to the same minimum and maximum values.

Applied a modified mean filter (kernel size 7 × 5) to remove speckles Figure 8b.

Bagci[31] in 2007 used correlation of axial A-scans to find six retinal layers and many papers[3145] used the same method for retinal thickness analysis in OCT datasets.

In 2008 Shrinivasan[26] modified the algorithm proposed by Koozekanani[23] to work on UHR OCT. In this study the thick scattering region of the outer retina previously attributed to the RPE was shown to consist of distinct scattering bands corresponding to the photoreceptor outer segment tips, RPE, and Bruch's membrane. Six retinal layers were totally segmented in this algorithm.

A-Scan methods lacked the contribution from 3D image context and suffered from excessive computation time and lack of layer detection accuracy; however, additional A-scan approaches have been introduced recently.[464748] Fabritius[46] incorporated 3D intensity information to improve the intensity based segmentation and segmented the ILM and RPE directly from the OCT data without massive pre-processing in about 17-21 seconds with error smaller than 5 pixels in 99.7% of scans.

B-Scan methods allowed dealing with 2D noise by incorporating better denoising algorithms during the preprocessing step. However, the dependency of these algorithms on noise reduction required very complicated and time-consuming denoising methods which made these algorithms too weak from the speed point of view.[202830] Additionally, the underlying intensity based methods and the relevant threshold selection was a difficult problem that made the methods case-dependent.

In 2006 Boyer[28] introduced a parabolic model of the cup geometry and used a modified version of the Markov model corporate by Koozekanani.[23] In this work, they extracted the parameters needed for clinical evaluations like the cup-to-disk ratio for the first time in TDOCT images.

One year later, in 2007, Baroni[30] used edge likelihood function maximizing for boundary detection. Along the A-scans, median filtering was applied (as 1D filter of 5 pixels); this eliminated isolated noisy pixels without losing small details. In the second step, the reflectivity values were summed row by row, to obtain a longitudinal cumulative grey level histogram [Figure 9a–d]. Detection of these peaks allowed an approximate delimitation of retina. Boundary detection was performed in the third step. The image was considered as overlapping cross-sectional stripes, 10 pixels wide. Each stripe was averaged to obtain a cross sectional grey level profile, which was filtered using the Gaussian gradient. The maxima and minima of the filtered profiles were detected by a simple peak detector and considered significant if their value was greater than the average value. The most important drawback of this method was that it assumed to have completely aligned images whose HRC layer was completely horizontal. In this study, for the first time, they reported the potential of thickness, densitometry, texture and curvature in TDOCT images for identifying retinal diseases.

Figure 9

Illustration of the computer method for the identification of three main retinal layers in normal eye (a) Regions of interest detected by horizontal cumulative histogram, (b) Approximate contours of ILM, inner and outer retina interface, and retina-pigment epithelium interface (c) Final contours; note the ILM contour smoothed for computing tortuosity index (see text), (d) RNFL strip used for reflectivity and texture analysis (it is dashed only for display purpose)

In 2008, Tan[49] utilized dynamic programming and 2D gradient information to extract boundaries and showed the relation between the thickness of retinal intra layers like RNFL, GCL and IPL and probability of suffering from glaucoma.

Recently in 2010, Kajic[50] used a Dual-tree complex wavelet (DTCW) denoising algorithm in preprocessing and utilized a statistical model based on texture and shape that captured the variance of the training data to segment eight layers on unseen data.

Active contours approaches for OCT image segmentation were first proposed by Cabrera Fernαndez[37] and modified by Yazdanpanah.[39] Unfortunately the required processing time and exact error reports are not available for any of the mentioned papers, which make such methods difficult to compare with other published methods. Regardless, active contour algorithms surpass the performance of intensity based B-scan approaches, both in resistance to 2D noise and in accuracy.

Cabrera Fernαndez[37] in 2004 applied active contours (a gradient vector flow (GVF) snake Model) to extract fluid-filled regions in the retinal structure of AMD patients.

One year later, Mujat[51] used deformable splines to calculate the thickness of RNFL in Two volumetric SDOCT data sets of the same eye from a single subject. The algorithm performed well in 350 frames with only few and isolated boundary detection errors. The advantage of this automated snake methodology was that it was able to provide larger area maps of the RNFL thickness facilitating the correct registration of ROIs with visual field defects which could allow better longitudinal evaluation of RNFL thinning in glaucoma studies.

In 2009, Yazdanpanah[39] introduced a modified Chan-Vese's energy-minimizing active contour algorithm in a multi-phase framework to segments SDOCT data from rodent models. Using the anatomical knowledge of expert, a circular shape prior was defined to avoid initial training. The approach had three main features. First, it could segment all intra-retinal layers simultaneously due to the multi-phase property of the algorithm. Second, they incorporated a shape prior term that enabled the algorithm to accurately segment retinal layers, even where the region-based information was missing, such as in inhomogeneous regions. Finally, the algorithm was robust and avoided the re-initialization problem that was associated with the level set approach.

In the same year, Mishra[32] used OCT data from healthy and diseased rodent retina and speckle noise and other typical artifacts in OCT images were handled by using an adaptive vector-valued kernel function in the precise layer boundary optimization step. For segmentation purpose, a modified active contour algorithm was utilized by using sparse dynamic programming method and two-step kernel based optimization scheme.

Recently in 2010, Mayer[52] used the minimization of an energy function consisting of gradient and local smoothing terms on 72 scans from glaucoma patients and 132 scans from normal subjects. A mean absolute error per A-Scan of 2.9 μm was achieved on glaucomatous eyes, and 3.6 μm on healthy eyes which proved that the approach provided a reliable tool for extracting diagnostic relevant parameters from OCT B-Scans for glaucoma diagnosis.

In 2011, Ghorbel[53] proposed a method for the segmentation of eight retinal layers in Heidelberg spectralis SDOCT images. The approach was based on global segmentation algorithms, such as active contours and Markov random fields. Moreover, a Kalman filter was designed to model the approximate parallelism between the photoreceptor segments and to detect them. The performance of the algorithm was tested on a set of 700 retinal images acquired in-vivo from healthy subjects.

Artificial intelligence based approaches were presented in[3329] and relied on a multiresolution hierarchical support vector machine (SVM) or on fuzzy C-means clustering techniques. The first one reported to have low ability in detection (six pixels of line difference and 8% of thickness difference) and a high time complexity (two minutes). But, the latter reported to have better results by only two pixels of linear difference and 45 seconds of time complexity. Overall, these methods cannot be categorized as established standard approaches since later-introduced methods like graph-based approaches can surpass them both in accuracy and time complexity.

In 2007 Fuller,[33] developed a semiautomatic segmentation system in which the morphology of retinal structures could be discovered and refined by a clinician. The clinician interactively specified the location of a retinal layer on a few select slices of the volume. This selection was then extrapolated throughout the entire volume using a SVM classifier in order to create segmentation. Once segmented, they provided visualizations and measurements of the resulting segmentation to aid in disease diagnosis. The main visualization interface was an interactive 3D volume rendering of the segmented portions of the volume. Furthermore, SVM considered a voxel's mean value and variance across multiple resolutions [Figure 10] in order to gracefully handle the speckle noise and to give the SVM a global perspective over feature shapes. They used a radial basis function kernel since they assumed it can represent the feature space well and also that speckle noise is normally distributed across the data. Additionally, it allows non-linear separation of the space. They reported that 68% of the thickness differences between the SVM segmentation and the manual segmentation fell below six voxel units.

Figure 10

Fuller[33] constructed a mipmap-like hierarchy in order to compute varying levels of data distributions to sample as input to a SVM algorithm. Topleft image shows a slice through a level having a resolution of 475×150×48. Remaining images show levels having a resolution half the one before it. (The change in the “bump” is due to the influence of neighboring slices in front and behind the shown slice.)

One year later in 2008, Mayer[29] calculated the RNFL thickness in circular SDOCT B-scans from 5 normal and 7 glaucoma eyes using Fuzzy C-means clustering technique without the need of parameter adaptation for pathological data. They reported that 97% of the upper and 74% of the lower RNFL layer boundary points lied within a two pixel range from the manual segmentation of the evaluation data set.

Vermeer[54] in 2011 presented a method for three-dimensional retinal layer segmentation in OCT images by a flexible method that learned from provided examples. Parts of representative OCT scans were manually segmented and used by the algorithms to learn from. Learning and classification of pixels was done by a support vector machine. Smoothness of the detected interfaces was guaranteed by the level set regularization that was applied after the pixel classification. The procedure was the same for all layers, except for the manually segmented data used to train the classifier. RMS errors for the top and bottom of the retina were between 4 and 6 μm, while the errors for intra-retinal interfaces were between 6 and 15 μm.

3D graph-based methods seem so far to be best suited for the task in comparison to the above-mentioned approaches. Their time requirements can be reduced to about 45 seconds per 3D volume (480 × 512 × 128 voxels) and they routinely achieve high accuracy with about 2.8 μm of layer-surface segmentation error. Such methods take advantage of newly developed 3D imaging systems, which provide better visualization and 3D rendering of the segmented boundaries.[27344042] By design benefitting from contextual information represented in the analysis graph, these methods are robust to noise and do not require advanced noise reduction techniques in the preprocessing steps. While there is no theoretical limit on the number of layers that can be simultaneously segmented by these approaches, up to 11 layers are routinely identified in retinal OCT images, performance that is unavailable to the other above-referenced algorithms.

In 2008, Garvin[36] developed an automated 3-D segmentation approach for the division of the retina on macular optical coherence tomography (OCT) scans into five layers. She compared its performance on OCT scans of patients with unilateral anterior ischemic optic neuropathy to that of human experts. Figure 11b and c shows an example of the six surfaces (labeled 1-6) they desired to find on each 3-D composite image Figure 11a. The 3-D composite image associated with each eye in Garvin's[36] work was created in two major steps. In the first step [Figure 12], raw scans for a particular angular location (e.g., all the vertical scans) were individually aligned so that boundary six (the retinal pigment epithelium) appeared approximately straight in the aligned image. Each scan was aligned by first finding boundaries 1, 5, and 6 simultaneously using an optimal graph search approach performed in 2-D. To ensure smoothness, a least-squares spline was fit to boundary six. The columns were then translated so that this spline would be a straight line in the aligned image. In the second step of this stage, each aligned image was registered to the first image in its location set by exhaustively searching for the best whole-pixel translation. After the determination of surfaces 1, 5, and 6, the remaining surfaces were found sequentially (allowing the utilization of other surface locations in the cost functions) in the following order: surface 4, surface 3, and finally, surface 2. The graph search approach guaranteed that the optimal feasible (satisfied smoothness and interaction constraints) surfaces would be found with respect to the designed cost functions. Figure 13 demonstrates three example results reflecting the best, median, and worst performances according to the overall unsigned border positioning error in garvin.[36] They reported overall mean unsigned border positioning error to be 6.1 ± 2.9 μm, a result comparable to the interobserver variability (6.9 ± 3.3 μm).

Figure 11

Example composite image (from Garvin[36]) with labeled intralayer segmentation and 3-D visualization of three surfaces (top and bottom of images have been cropped to aid in visualization) (a) Composite image (b) Six surfaces (labeled 1-6) and five corresponding intralayers (labeled A-E). The anatomical correspondence is our current presumption based on histology and example images from higher-resolution research OCT scanners: (a) nerve fiber layer (NFL), (b) ganglion cell layer and inner plexiform layer (GCL + IPL), (c) inner nuclear layer and outer plexiform layer (INL+OPL), (d) outer nuclear layer and photoreceptor inner segments (ONL + IS), (e) photoreceptor outer segments (OS). (c) Example 3-D visualization of surfaces 1, 3, and 4[36]

Figure 12

Individual scan alignment (top and bottom of images have been cropped to aid in visualization)[36]

Figure 13

Three example results reflecting the best, median, and worst performances according to the overall unsigned border positioning error in garvin.[36] (a) Best case composite image, (b) Best case composite image with segmented borders, (c) Best case composite image with average manual tracing, (d) Median case composite image, (e) Median case composite image with segmented borders, (f) Median case composite image with average manual tracing, (g) Worst case composite image, (h) Worst case composite image with segmented borders, (i) Worst case composite image with average manual tracing[36]

In 2009, Abràmoff[40] combined a multiscale 3D graph search algorithm and a voxel column classification algorithm using a k-NN classifier to segment the ONH cup and rim in SDOCT data from 34 glaucoma patients. This preliminary study showed for the first time a high correlation between segmentation results of the ONH cup and rim from SDOCT images and planimetry results obtained by glaucoma experts on the same eye. They reported the correlation of algorithm to three independent experts was 0.90, 0.87, and 0.93, respectively.

In 2010, Lee[27] presented an improved and fully automatic method based on a similar methodology using graph search combined with a k-NN classifier that employed contextual information combined with a convex hullbased fitting procedure to segment the ONH cup and rim on 27 SDOCT scans from 14 glaucoma patients. They reported unsigned error for the optic disc cup was 2.52 ± 0.87 pixels (0.076 ± 0.026 mm) and for the neuroretinal rim was 2.04 ± 0.86 pixels (0.061 ± 0.026 mm).

In the same year, Yang[41] proposed a method in which detection of 9 boundaries followed a two-step segmentation algorithm Figure 14. First, a customized canny edge detector was used to create a map showing local main edges. A complementary gradient map in the axial direction was also acquired with a larger kernel. Second, a graph was built based on a combination of the axial intensity gradient and the canny edge maps. The layer boundary was then extracted bythe shortest path search applied to the graph using a dynamic programming algorithm. Overall the ICC of each boundary was above 0.94, the mean coefficient of variation was less than 7.4%, and the mean standard deviation was less than 2.8 μm. Algorithm was able to segment low intensity and low contrast OCT images in a very short time without degrading the accuracy. In addition, pre extraction of vessel locations, which is not a trivial operation, was unnecessary in this method, as shown in Figure 15.

Figure 14

Illustration of nine boundary segmentation flow in Yang.[41] Panel a is the original OCT image acquired using Topcon 3D OCT-1000 equipment. The image was first aligned as shown in b and the ILM and IS/OS were detected as in c; d, e and f illustrate the BM/Choroid, OS/RPE, IPL/INL, NFL/GCL, GCL/ IPL, INL/OPL and ELM were detected in order; in the end, all nine boundaries were converted back to the original OCT image coordinates as shown in g[41]

Figure 15

Illustration of segmentations on images with blood vessel artifacts[41]

Again in 2010, Quellec[34] introduced a novel ten-layer automated segmentation of OCT data from 13 normal eyes and from 23 eyes with CNV, intra-, and subretinal fluid and pigment epithelial detachment and twenty-one textural features were implemented in 3-D and measured locally for each layer. The variations of texture and thickness across the macula in these ten layers, averaged over thirteen normal eyes, defined the normal appearance of maculae in SD-OCT scans. A machine learning approach that classifies the retinal pathology based on feature- and layer-specific properties in comparison with the normal appearance of maculae was reported. In this study, they reported an extended method for automated segmentation of 10 intraretinal layers identified in [Figure 16] from 3-D macular OCT scans using multiscale 3-D graph search technique. The basic concept of this approach was to detect the retinal surfaces in a subvolume constrained by the retinal surface segmented in a low-resolution image volume. The cost functions for the graph searches, capable of detecting the retinal surfaces having the minimum costs, were inverted gradient magnitudes of the dark-to-bright transition from top to bottom of the OCT volume for surfaces 1, 3, 5, 7, 9, 10 and those of the bright-to-dark transition for surfaces 2, 4, 6, 8, 11 [Figure 16]. They reported that mean unsigned surface positioning errors were less than 6 μm and confirmed that useful 3-D textural information can be also extracted from SD-OCT scans to aid local retinal abnormality detection.

Figure 16

Segmentation results of 11 retinal surfaces (10 layers) (a) X-Z image of the OCT volume (b) Segmentation results, nerve fiber layer, ganglion cell layer, inner plexiform layer, inner nuclear layer, outer plexiform layer, outer nuclear layer, outer limiting membrane, inner segment layer, connecting cilia, outer segment layer, Verhoeff's membrane, and retinal pigment epithelium. The stated anatomical labeling is based on observed relationships with histology although no general agreement exists among experts about precise correspondence of some layers, especially the outermost layers. (c) 3-D rendering of the segmented surfaces (N: nasal, T: temporal)

Furthermore in 2010, chiu[55] reported a skillful approach based on graph-based theory and dynamic programming thatsignificantly reduced the processing time required for image segmentation and featureextraction. This methodology was able to address sources of instability such asthe merging of layers at the fovea, uneven tissue reflectivity, vessel hypo-reflectivity and thepresence of pathology. Interestingly, the approach incorporates an automatic initialization thatbypasses the need for manual endpoint selection.

In 2012, Kafieh[42] proposed a segmentation method capable of detecting 12 retinal boundaries using diffusion map based segmentation algorithm. In contrast to recent methods of graph based OCT image segmentation, this approach did not require edge-based image information and rather relied on regional image texture. Consequently, the method demonstrated robustness in situations of low image contrast or poor layer-to-layer image gradients. Diffusion map was applied on 2D and 3D OCT datasets and in each application; the procedure was composed of two steps (one for partitioning the data to important and less important sections, and another one for localization of internal layers). In the first step, the data pixels/voxels were grouped in rectangular/cubic sets to form a graph node. The weights of graph were also calculated based on geometric distances of pixels/voxels and differences of their mean intensity. The first diffusion map clustered the data into three parts, second of which was the area of interest and the two other sections were eliminated from the next calculations. In the second step, the remaining area went through another diffusion map algorithm and the internal layers were localized based on their similarity of texture Figure 17. Figure 18shows an example of extracted OCT layers using the mentioned diffusion-map based method.[42] The method was tested on 23 datasets from two patient groups (10 datasets, 3D-OCT data from patients diagnosed with glaucoma and 13 datasets, 3D-OCT obtained from normal eyes without pathologies). The mean unsigned border positioning errors (mean ± SD) was 8.52 ± 3.13 and 7.56 ± 2.95 micrometers for the 2D and3D approaches, respectively. In addition, pre extraction of vessel locations was unnecessary in this method, as shown in Figure 19.

Figure 17

Construction of graph nodes from a 3D OCT[42]

Figure 18

(Up-left) Surface 1, (Up-right) First and 7th surfaces, (Down-left) Surfaces 9 to12, (Down-right) Surfaces 3 to7[42]

Figure 19

Robustness of the proposed algorithm to blood vessel artifacts[42]

Table 2 shows a brief look at the mentioned approaches and compares the OCT systems, preprocessing method, error range, and computation time.

Table 2

A brief look at the mentioned approaches and compares the OCT systems, preprocessing method, error range, and computation time

There is no doubt that algorithms and research projects work on a limited number of images with some determinate abnormalities (or even on normal subjects) and such limitations make them more appropriate for bench and not for the bedside. The automatic segmentation will never replace the physicians but the improved segmentation methods and the ability of these methods in providing more information in less complicated data will offer valuable help to ophthalmologists in better diagnosis and treatment of retinal pathologies.

CONCLUSION

In contrast to OCT technology development which has been a field of active research since 1991, OCT image segmentation has only been more fully explored during the last decade. Segmentation, however, remains one of the most difficult and at the same time most commonly required steps in OCT image analysis. No typical segmentation method exists that can be expected to work equally well for all tasks.[2] In this paper, we tried to cite most related works from 1997 to 2012, however, this is in no way complete. It should also be noticed that the number as reported in Tables cannot be used for direct comparison of the relative performances, since different settings are utilized in each method.

Current researches in OCT segmentation are mostly based on improving the accuracy and precision, and on reducing the required processing time. There is no doubt that current 3-D imaging modalities are now moving the research projects toward volume segmentation along with 3-D rendering and visualization. It is also important to develop robust methods capable of dealing with pathologic cases in OCT imaging.

Eventually, it should be emphasized that the automated segmentation methods can be categorized as an undeniable assistant to physicians and there is no doubt that such methods can never replace a trained doctor.

BIOGRAPHIES

Rahele Kafieh received the BS degree from Sahand University of Technology, Iran, in 2005 and the MS degree from the Isfahan University of Medical Sciences, Iran, in 2008, both in biomedical engineering. She is currently working toward the PhD degree in the Department of Biomedical Engineering at Isfahan University of Medical Sciences. Her research interests are in biomedical image processing, computer vision, graph algorithms, and sparse transforms. She is a student member of the IEEE and the IEEE Signal Processing Society

E-mail: r_kafieh@yahoo.com

Hossein Rabbani is an Associate Professor at Isfahan University of Medical Sciences, in Biomedical Engineering Department also Medical Image & Signal Processing Research Center. Involved research topics include medical image/volume processing, noise reduction and estimation problem, image enhancement, blind deconvolution, video restoration probability models of sparse domain&s coefficients especially complex wavelet coefficients. He is a member of IEEE, Signal Processing Society, Engineering in Medicine and Biology Society, and Circuits and Systems Society

E-mail: h_rabbani@med.mui.ac.ir

Saeed Kermani obtained his BS from the Department of Electrical Engineering of Isfahan University of Technology in Isfahan, Iran, 1987, and he received the MS in Bioelectric Engineering from Sharif University of Technology, in 1992 and his PhD in Bioelectric Engineering at AmirKabir University of Technology, Tehran, Iran, in 2008. He is Assistant Professor of Medical Engineering at the Department of Medical Physics and Medical Engineering in the School of Medicine of Isfahan University of Medical Sciences, Iran. His research interests are in biomedical signal and image processing techniques

E-mail: kermani@med.mui.ac.ir