|Home | About | Journals | Submit | Contact Us | Français|
This paper presents a generalized framework for segmenting closed-contour anatomical and pathological features using graph theory and dynamic programming (GTDP). More specifically, the GTDP method previously developed for quantifying retinal and corneal layer thicknesses is extended to segment objects such as cells and cysts. The presented technique relies on a transform that maps closed-contour features in the Cartesian domain into lines in the quasi-polar domain. The features of interest are then segmented as layers via GTDP. Application of this method to segment closed-contour features in several ophthalmic image types is shown. Quantitative validation experiments for retinal pigmented epithelium cell segmentation in confocal fluorescence microscopy images attests to the accuracy of the presented technique.
Fast, accurate, and objective detection and quantification of imaging biomarkers are crucial for the study and diagnosis of ophthalmic diseases. Due to the time consuming and subjective nature of manual segmentation, considerable work has been done in recent years to automate the segmentation of ocular structures. For example, many layer segmentation algorithms have been developed to delineate the retinal [1–6], choroid-scleral [7,8], and corneal layers [9,10] on optical coherence tomography (OCT) images .
However, in addition to layered structures, there is also a need to automatically segment and quantify closed-contour features in ophthalmic images. Examples of closed-contour anatomical and pathological structures include retinal pigmented epithelial (RPE) cells in confocal fluorescence microscopy images of flat-mounted retina , intra-retinal cysts in OCT images of patients with diabetic macular edema (DME) [13,14], and photoreceptors in adaptive optics scanning ophthalmoscope (AOSO) images of retina [15–18]. Previous works have been presented to segment closed-contour features, including fluid or detachments on OCT images of the retina [19–21], drusen in fundus photographs [22,23] or in OCT images , geographic atrophy in fundus auto-fluorescence images , vessel lumens in intravascular OCT images , corneal cells [27–31], and photoreceptors on AOSO images [32–37]. However, these segmentation algorithms were often developed for specific ophthalmic applications, rather than as a general framework applicable to both layer and closed-contour features. Furthermore, in some cases such as with photoreceptor or corneal cell identification, cell counts or densities were achieved but the actual structure size and shape were not necessarily determined [27,29,31–35,37].
To segment circular biological structures, transformation of images into the polar domain has proven to be an effective tool within the medical community [38–41]. This is due to the observation that circles in the Cartesian domain are represented as lines in the polar domain. Using this transform, the mathematical operations and manipulations required to segment the image are greatly simplified since these operations no longer need to be performed radially. However, for non-circular, closed-contour structures, the polar transform is no longer adequate since oblong or convoluted features deviate from a flat path.
In this paper, we present a general framework to segment arbitrarily shaped, closed-contour features. We previously developed a framework for layer segmentation based on graph theory and dynamic programming (GTDP)  to segment retinal and corneal layers on spectral domain optical coherence tomography (SDOCT) images [5,6,10]. Section 2 of this paper extends the application of our layer segmentation technique to include closed-contour segmentation. The algorithm is based on the projection of closed-contour features in the Cartesian domain into lines in a virtual domain, which we term the quasi-polar domain. In the quasi-polar domain, features of interest are then segmented as layers using our GTDP technique. In Section 3, we apply the framework to segment RPE cells in confocal microscopy images of an age-related macular degeneration (AMD) mouse model, and in Section 4.1 we evaluate the algorithm’s effectiveness in detail. Lastly, to highlight the flexibility of the presented framework, Section 4.2 illustrates the segmentation of closed-contour features for other ophthalmic features, including cysts on SDOCT and cone photoreceptor cells on AOSO images.
This section introduces a generalized framework for closed-contour feature segmentation. The method is an extension of our previously presented technique for layer segmentation via GTDP . The core steps are outlined in Fig. 1 and discussed in detail in the following subsections. Section 3 then provides the specific algorithmic steps required to implement this framework for RPE cell segmentation in confocal fluorescence microscopy images.
The first step is to attain pilot estimates of all structures to segment. These estimates should contain information about both the location and shape of the structures. Structure localization prevents the algorithm from confusing the target structure with similarly-shaped extraneous features, while shape estimation results in a more accurate quasi-polar transform, as will be discussed in Section 2.2.
For a grayscale image in the Cartesian domain, the pilot estimates can be represented as a binary image where
and denotes the pixel at the ith row and jth column of the image. Many basic techniques such as thresholding, contours, local minima, or the multiscale generalized likelihood ratio test  can be used to attain this estimate, and Section 3.2 describes one such method for pilot structure estimation. To isolate a single structure, and can be cropped to smaller image blocks both where and
This section explains the quasi-polar transform in depth. In its simplest form, the quasi-polar transform can be described by a polar transform followed by a flattening step. Figure 2 shows an illustration, where a circle, an oval, and an arbitrary shape in the Cartesian domain (Figs. 2(a-c)) are first transformed into the polar domain (Figs. 2(d-f)). Information about the shape of the pilot estimate is then used to flatten the structure into a line in the quasi-polar domain (Fig. 2(g)). Once transformed into the quasi-polar domain, the structure can be segmented as if it were a layer using the GTDP technique.
It can also be observed from Fig. 2 that a simple polar transform is sufficient to map closed-contour shapes (Figs. 2(a-c)) into layered structures (Figs. 2(d-f)). Motivation for a further transformation into the quasi-polar domain (Fig. 2(g)) stems from the tendency for shorter geometric paths to be found when using the GTDP technique, especially in low signal-to-noise ratio (SNR) imaging scenarios. This is due to the observation that a cut (segmentation) with fewer traversed nodes often yields a lower total weight. As a result, oblong features are disadvantaged since their paths in the polar domain are longer than the shortest geometric path (Figs. 2(e, f)). By performing an extra flattening step, the path of the oblong structure now reflects the shortest geometric distance.
The quasi-polar transform is implemented using the following three steps:
Since Step 3 flattens the pilot estimate instead of the structure itself, in general the structure is not perfectly flat in This also implies that a pilot estimate shape closer to the actual structure results in a flatter structure in the quasi-transform domain.
Once the image is transformed into the quasi-polar domain, it can be segmented using our GTDP-based automatic layer segmentation algorithm . In summary, this technique represents the image as a graph containing nodes (i.e. image pixels), and edges, which connect adjacent pixels. We assign weights to each of the edges based on a priori information about the structure to prefer paths through the structure boundaries. Our GTDP technique includes automatic endpoint initialization, where an additional column is added to either side of the image with minimal vertical edge weights. As we detailed in , the start and end nodes can then be arbitrarily defined anywhere within these two columns. Finally, any optimization method that finds the minimum weighted path from the start node to the end node (e.g. Dijkstra’s method ) can be used to segment the structure.
Once the structure is segmented in the quasi-polar domain, the segmentation information is transformed back into the Cartesian domain by applying the inverse of Steps 1 and 3 from Section 2.2.
After the first structure is segmented, Steps 1-3 in Section 2.2 are repeated for each subsequent structure (e.g. for each individual cell or cyst). We can achieve this either by sequential segmentation of individual features, or by parallelizing the segmentation to make this process more efficient.
AMD is the leading cause of irreversible blindness in Americans older than 60 years . In a recent study, Ding et al. showed that anti-amyloid therapy protects against RPE damage and vision loss in an APOE4 mouse model of AMD . In , the average size and density of the cells were used as biomarkers for disease progression. To attain these quantitative metrics, 22,495 cells were analyzed using customized semi-automatic segmentation software. To achieve a drastically faster segmentation rate, this section describes an implementation of the generalized closed-contour segmentation algorithm from Section 2 to fully automatically segment RPE cells in confocal microscopy images of flat-mounted retina of normal and AMD mouse models.
In the following, we describe our fully automatic segmentation method as tailored to the data set from Ding et al.’s report published in PNAS . However, the adaption of our method to any other data set is straightforward. The data set in  included 23 confocal microscopy fluorescence images (0.64 × 0.64 mm2) of central RPE flat mounts from 17 mice with approximately 1000 cells per image. RPE cells were stained with a rabbit antibody against ZO-1 (40–2200, 1:100; Invitrogen) and imaged on a Nikon Eclipse C1 microscope.
We first normalized the intensities of the image from  We then smoothed using an 11 × 11 pixel Gaussian filter with standard deviation pixels. Next, we computed the extended-minima transform  that found all minima and suppressed those with a depth less than 0.05. To generate the pilot estimate for all cells, we set for all points z lying within the convex hull  of each minima. To localize the desired cell and its corresponding pilot estimate, we first cropped and to 201 × 201 pixel image blocks corresponding to the largest conceivable cell size in our model. Figure 3(a) shows the cropped cell image, and Fig. 3(b) shows the cropped pilot estimate image. Both images were centered at the computed for the desired pilot estimate. Since contained other cell estimates, we then set all pixels outside the desired estimate in to zero.
We next transformed (Fig. 3(a, d, f)) by following Steps 1-3 in Section 2.2 with the following implementation details:
While Fig. 3(f) shows the resulting quasi-polar transformed cell image Fig. 3(e) shows the quasi-polar transformed pilot estimate. As can be seen, the transform yielded a fully flattened pilot estimate and an approximately flattened cell structure.
Once the quasi-polar transformed image was generated, we used our GTDP layer segmentation framework to segment the cell boundary. To do so, we defined the following three cost functions:
We implemented these three cost functions using the following MATLAB (The MathWorks, Natick, MA) notation, where and are corresponding columns in and respectively, such that
We then combined these cost functions such that
Finally, we calculated the edge weights using
where is the weight assigned to the edge connecting pixels and and The notation indicates a normalization of the values x to range from y to z
Prior to finding the shortest path, we discarded all edges from the graph containing nodes in previously segmented cells to avoid segmenting them again. We also excluded the edges that allowed those cells to be included within a newly segmented cell. We then found the minimum-weighted path p using Dijkstra’s method, where and represent the row and column axes of respectively. Since in the quasi-polar domain, must equal in practice however, this is not guaranteed when using automatic endpoint initialization (Section 2.3). Therefore, if we set both endpoints equal to and found the corresponding minimum-weighted path We then set both endpoints equal to and found the path The ultimate selected path p (Fig. 3(f), magenta) was the one that passed through the brighter overall path, where
We transformed the segmentation back into the Cartesian domain (Fig. 3(g)) by following the method in Section 2.4. However, since the cell was segmented using only a cropped section of the image, we shifted the coordinates of the segmented layer back to its appropriate location on the original-sized image The result was a binary edge image where all pixels coincided with the segmentation such that For all pixels and that were unconnected, we used linear interpolation to connect the dots and create a closed-contour.
To segment all other cells, we iterated the steps outlined in Sections 3.3 through 3.6 for all pilot estimates of the cells in the image. To prevent gaps from occurring between adjacent cells, we created a preference for already-segmented cell borders by brightening their intensity values, where and for all z where
The gold standard segmentation from  used in our validation study (Section 4.1) separated cells by drawing a line through the middle of the cell border with a single pixel thickness. Since the presented algorithm did not necessarily yield a single-pixel cell border, this section describes the post-processing steps required to match our automatic segmentation with the gold standard protocol, and to remove any erroneous cell segmentations.
First, we removed all cells smaller than 50 pixels by setting for all z lying within the cell. We then generated the skeleton  of and removed all spurs . Next, we looked at individual cell edges by setting the branch points in equal to zero. We defined the minimum edge threshold to be where was the set of all pixels where If all pixels for a given cell edge fell below t, or if 80% of the pixels fell below t and there were at least 15 such pixels, then the cell edge was removed from Once we removed all invalid cell edges, we restored all branch points on and removed any newly created spurs or pixels no longer connected to the cell edge.
This section shows the closed-contour segmentation results obtained using the procedures described in Section 3. Section 4.1 evaluates in detail the RPE cell segmentation algorithm described in Section 3, while Section 4.2 briefly illustrates the segmentation of closed-contour features for other ophthalmic applications.
To validate our fully automatic segmentation results against the gold standard, we compared the cell count and mean cell size for each of the 23 images in . To prevent from biasing our results, these images were not used during algorithm development stage. Instead, a single confocal microscopy fluorescence image of an AMD mouse model not included in the validation data set was used to train the algorithm.
We did not alter the data in Ding’s study in any shape or form, with one slight exception. In , closed-counter features smaller than 100 pixels were absorbed into the neighboring cells, since such regions were likely a result of incorrect manual segmentation. In our replication of the results from , we instead discarded these regions from our quantitative comparison since such regions were negligible. In , regions of the image considered invalid due to imaging errors, flat-mount preparation errors, or cells being cut off by image borders, were marked separately and discarded from the analysis. In this study, we used the markings from  to ignore the exact same invalid regions in our analysis. Note that, the above two types of outliers occupied a very small percentage of the total area segmented.
Quantitative results for the RPE cell segmentation algorithm described in Section 3 are shown in Table 1 . The average error between automatic segmentation and the gold standard was 1.49 ± 1.44% for the cell count and 1.32 ± 1.53% for the mean cell area. The median error was found to be 0.97% and 0.71% for the detected number of cells and mean cell area, respectively. Qualitative results for this validation study are shown in Fig. 4 , where Fig. 4(a) (corresponding to Image 16 in Table 1) is the automatically segmented confocal fluorescence image of an AMD flat-mounted mouse retina, and Figs. 4(b-g) are zoomed-in cell segmentation results. Figure 5 shows two confocal images of RPE cells (Figs. 5(a, d)) along with their gold standard (Figs. 5(b, e)) and automatic (Figs. 5(c, f)) segmentation results. The image shown in Figs. 5(a-c) corresponds to Image 9 in Table 1 with approximately median error, and Figs. 5(d-f) show Image 10 from Table 1 with maximum error. Lastly, Fig. 6 shows zoomed-in images of Images 6 and 9 from Table 1 comparing the gold standard (Figs. 6(a-e)) to automatic segmentation (Figs. 6(f-j)). The images in Figs. 6(a-c) are examples of erroneous gold standard segmentation, while Fig. 6(j) is an example of erroneous automatic segmentation. Their corresponding images (Figs. 6(e-h)) are examples of accurate segmentation. Finally, Figs. 6(d, i) show an undetermined case where it is unclear which segmentation is correct.
All 23 images used in the study and their corresponding manual and automatic segmentation data are available at http://www.duke.edu/~sf59/Chiu_BOE_2012_dataset.htm. The “single click” algorithm was coded in MATLAB (The MathWorks, Natick, MA) and resulted in an average computation time of 2.95 minutes per image (with an average of 964 cells per image) on a desktop computer with a 64-bit operating system, Core i7 2600K CPU (Intel, Mountain View, CA), solid state hard drive, and 16 GB of RAM. This time included the overhead required for reading and writing operations.
Section 2 presents a generalized segmentation framework for closed-contour structures, and Sections 3 and 4.1 discuss the implementation and results for RPE cell segmentation. We have additionally implemented the algorithm to segment other ophthalmic features. These include intra-retinal cysts seen on SDOCT images (Spectralis; Heidelberg Engineering, Heidelberg, Germany) of patients with DME (Fig. 7 ). We have also included automatic segmentation of cone photoreceptors on AOSO images  of the retina (Fig. 8 ). To segment the cysts and photoreceptors, the same basic framework as described in Section 3 was utilized. Changes required to adapt the algorithm to different image types included a modification of image filters and pilot estimation techniques, changing the edge weights, and different refinement steps to exclude erroneous segmentation (e.g. removing segmentations of hypo-reflective blood vessels in SDOCT images).
To the best of our knowledge, no other technique has been reported for fully automatically segmenting confocal fluorescence images of RPE cells. The qualitative results shown in Fig. 4 demonstrate accurate automatic cell segmentation despite the presence of image artifacts, low signal-to-noise ratio, or inaccurate estimations. Furthermore, the cell shape and size were not constraining factors for the presented automatic algorithm.
Quantitative results show that our fully automatic algorithm accurately segmented RPE cells on confocal images of flat-mounted mouse retina with AMD, yielding cell count and mean cell area measurements with an average error of less than 1.5%. Automatic segmentation errors similar to Fig. 6(j) occurred due to the presence of bright imaging artifacts, which erroneously altered values in the first cost function of Eq. (6). However, as shown in Figs. 6(a-c), even the gold standard segmentation reported in PNAS contained errors as well. Such errors will naturally occur in most manual or semi-automatic segmentation tasks required for large-scale studies, where multiple human experts subjectively grade images at different dates. Furthermore, since manual correction is extremely time consuming, cells with inconsequential errors (Fig. 6(b)) may not have been a priority to correct. As a result, a 0% error from the gold standard does not imply perfect segmentation.
To further reduce the effect of pilot estimation and achieve more accurate segmentation, one may employ an iterative approach, in which the result of segmentation for one iteration is used as the pilot estimate for the next iteration. In addition, application of modern adaptive denoising methods to raw images, which remove imaging artifacts without altering the anatomic structures, may further improve the automated segmentation performance [48,49].
Our motivation for using cell count and area as validation metrics stems from the need to quantify the extent of cell damage present in different groups of diseased mice . While error in the absolute position of each automatically segmented cell may have been a more direct validation metric, the presence of errors in the gold standard made this unfeasible. Furthermore, the cell boundary is several pixels thick, making it difficult to assign a “true” boundary position. While the gold standard divided cells often through the middle of the cell border, this was not always the case. As a result, we deemed cell count and area to be a more meaningful and viable validation metric.
Finally, we showed preliminary results demonstrating the segmentation of cysts seen in SDOCT images of eyes with DME (Fig. 7) and cone photoreceptors seen in AOSO images of the retina (Fig. 8). Prior studies on retinal fluid segmentation have focused on the detection of larger regions of intra-retinal fluid  or fluid constrained by retinal layer segmentation . However, our algorithm is sensitive to the small changes in intra-retinal fluid buildup. With regards to cone photoreceptor segmentation, many previous works assumed a given photoreceptor size, shape, or packing [32–35,37]. Our method relaxes these constraints, allowing for a future version of the algorithm capable of segmenting both photoreceptor rods and cones in a single image. Lastly, while we plan to validate these algorithms in future studies, we believe these preliminary results are encouraging, as they show the flexibility of our algorithm to be extended to different structures, diseases, and imaging modalities.
In summary, we demonstrated the extension of our automatic GTDP framework to segment not only layered, but also closed-contour structures. Implementation of the algorithm for RPE cell segmentation in confocal fluorescence images of flat-mounted AMD mouse retina resulted in an accurate extraction of cell count and average cell size. We believe that such a tool will be extremely useful for future studies, which use cell morphology as a biomarker for the onset and progression of disease. While validation has yet to be shown for other ophthalmic applications, we have demonstrated preliminary results for intra-retinal cyst segmentation in a SDOCT image as well as cone photoreceptors segmentation in an AOSO image. This is highly encouraging for reducing the time and manpower required for segmenting closed-contour features in ophthalmic studies.
We would like to thank Dr. Jindong Ding for his supervision on the semi-automatic segmentation of RPE cells, Dr. Priyatham Mettu for his expertise on DME cyst segmentation, Prof. Glenn Jaffe for providing the image in Fig. 7(a), and Linkun (Sam) Xi for manual correction of the RPE cells. This work was supported in part by The American Health Association Foundation, NIH (1R21EY021321-01A1), Research to Prevent Blindness (Duke’s 2011 Unrestricted Grand Award), the North Carolina Biotechnology Center, and the John T. Chambers Scholarship.