Scholarly article on topic 'Fully automatic evaluation of the corneal endothelium from in vivo confocal microscopy'

Fully automatic evaluation of the corneal endothelium from in vivo confocal microscopy Academic research paper on "Medical engineering"

Share paper
Academic journal
BMC Medical Imaging
OECD Field of science

Academic research paper on topic "Fully automatic evaluation of the corneal endothelium from in vivo confocal microscopy"

Selig et al. BMC Medical Imaging (2015) 15:13 DOI 10.1186/s12880-015-0054-3


^BMC Medical Imaging

Open Access

Fully automatic evaluation of the corneal endothelium from in vivo confocal microscopy

Bettina Selig1, Koenraad A Vermeer2, Bernd Rieger3, Toine Hillenaar4 and Cris L Luengo Hendriks1*


Background: Manual and semi-automatic analyses of images, acquired in vivo by confocal microscopy, are often used to determine the quality of corneal endothelium in the human eye. These procedures are highly time consuming. Here, we present two fully automatic methods to analyze and quantify corneal endothelium imaged by in vivo white light slit-scanning confocal microscopy.

Methods: In the first approach, endothelial cell density is estimated with the help of spatial frequency analysis. We evaluate published methods, and propose a new, parameter-free method. In the second approach, based on the stochastic watershed, cells are automatically segmented and the result is used to estimate cell density, polymegathism (cell size variability) and pleomorphism (cell shape variation). We show how to determine optimal values for the three parameters of this algorithm, and compare its results to a semi-automatic delineation by a trained observer. Results: The frequency analysis method proposed here is more precise than any published method. The segmentation method outperforms the fully automatic method in the NAVIS software (NidekTechnologies Srl, Padova, Italy), which significantly overestimates the number of cells for cell densities below approximately 1200 mm-2, as well as previously published methods.

Conclusions: The methods presented here provide a significant improvement over the state of the art, and make in vivo, automated assessment of corneal endothelium more accessible. The segmentation method proposed paves the way to many possible new morphometric parameters, which can quickly and precisely be determined from the segmented image.

Keywords: Frequency analysis, Stochastic watershed segmentation, Image-based quantification, In vivo imaging, Corneal endothelium


Corneal endothelium

The cornea is the frontal, transparent layer of the eye covering the pupil and iris. Its endothelium is composed of a monolayer of hexagonal cells, and plays a pivotal role in the homeostasis of the cornea. Ion exchangers at the basolateral border of the endothelial cells create an osmotic gradient that passively drives corneal fluid towards the anterior chamber. By maintaining the water fraction at 78%, the electrolyte balance ensures corneal transparency. Distortion of this equilibrium caused by damage of endothelial cells leads to swelling of the


1 Centre for Image Analysis, Uppsala University, Box 337, 75105 Uppsala,


Full list of author information is available at the end of the article

hydrophilic stroma and a consequent visual loss. In the normal cornea, endothelial cell density slowly declines with age [1]. Several factors such as intraocular surgery or inflammation, however, may speed up the endothelial cell loss. Since the endothelium has very limited regeneration capacity, a decrease in endothelial cell density is compensated by thinning and elongation of the remaining cells. At the same time, the ion exchangers are upregulated to uphold the fluid pumping mechanism [2]. Corneal transparency is affected when the decrease in endothelial cells becomes too large and the pumping mechanism fails. A critical endothelial cell density of 300-500 mm-2 seems needed to prevent corneal decompensation [3]. In ophthalmic practice, this threshold is often used to determine whether or not a patient requires corneal transplantation. The donor cornea, in turn, should have a density of at least

O© 2015 Selig et al.; licensee BioMed Central. This is an Open Access article distributed under the terms of the Creative Commons Central Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the originalwork is properly credited. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

2000 mm-2 to be eligible for keratoplasty. Endothelial cell quantification is also used as a screening method to elect patients for corneal refractive surgery and to determine the rate of endothelial cell loss after iris claw lens implantation.

Even if cell density is the single most important measure when judging the condition of the corneal endothelium, other measures such as variation of cell size and shape are, in principle, also interesting. These measures are, however, not often used in the clinic because of the difficulty in obtaining reliable results. The frequent clinical use of endothelial cell density measurement and its important consequences warrant the development of an objective, fully automatic quantification method. In this paper we introduce two such methods, using images from in vivo white light slit-scanning confocal microscopy. The first simply obtains an estimate of the cell density. The second method produces a delineation of all cells in the field of view. From such a delineation, it is possible to obtain many different morphometric parameters. These methods will enable the use of other morphometric parameters in clinical practice.

State-of-the-art analysis tools

It is common practice to use commercially available semiautomatic software tools (such as NAVIS by Nidek, IMAGEnet by Topcon, Bambi, Tomey, or EAT by Rhine-Tec) to delineate endothelial cell images and determine morpho-metric quantities. Even though the procedure is faster and less subjective than manual segmentation, semi-automatic analysis is still tedious and time consuming. The automatic approaches that are commercially available (such as the automatic segmentation mode in the NAVIS software) do not give satisfactory results [4,5] (compare Figures 1b-d).

Methods that automatically estimate the cell density using frequency analysis seem promising [6-9], but the commercial implementations give poor results for both low and high cell densities [10].

Several automatic segmentation methods for corneal endothelium have been proposed recently [11-13]. Of these, only Gavet and Pinoli [11] attempts to solve the problem for images obtained in vivo. These authors use a specular microscope, which produces similar images to those of the confocal microscope used in this paper. The other two papers concern images of in vitro specular microscopy [12] and in vitro inverse contrast phase microscopy [13].

Frequency analysis

Due to the fairly regular formation of the corneal endothelium (Figure 1), frequency analysis is commonly used to obtain an average cell size. In the late 1980's and early 1990's, Masters and colleagues [14-16] observed a

concentric ring in the Fourier transform of an image of corneal endothelium. The ring's radius was related to the average cell size. Some years later, this relation was further explored through the analysis of synthetic images [8,17,18].

The procedure of the most recent approaches to estimate endothelial cell density [6,7,9] is to project the two-dimensional frequency spectrum onto a one-dimensional function. One salient spatial frequency on this function is then found, often a peak. This frequency, which we will here call the characteristic frequency f*, is related to the average cell size, and thus to cell density. We review these methods, and propose a new one, under the heading "Frequency analysis" in the "Methods" section.

All published approaches make the same, implicit assumption A = if*)-2 about the relationship between f* and the average cell area A. This is, in general, not correct, and under the heading "Density estimation" in the "Methods" section we show how to establish the correct relationship.

The stochastic watershed algorithm

For more complex analyses of the endothelial layer, a full segmentation of the cells is needed. In confocal images of the endothelial layer, the cells are visualized as bright areas separated by dark lines (Figure 1). An excellent method to segment such images is the watershed algorithm [19,20]. This algorithm segments the image into regions with one local minimum each. The segmentation lines run along the brightest line in between each pair of local minima. Thus, the confocal images of the endothelium must be inverted first, such that the cell boundaries are bright lines. Unfortunately, due to the noise in the image, each cell typically has many local minima, and therefore will be segmented into many regions by the watershed algorithm. One solution is to apply the H-minima transform before the watershed algorithm [21]. This transform suppresses all local minima that have a depth less than the parameter h. However, it is often impossible to find a value h such that each cell has exactly one local minimum.

An alternative to avoid oversegmentation is to apply a seeded watershed [22]. One seed point needs to be placed in each cell; a minimum imposition algorithm will transform the image such that these seeds will be the only local minima in the image. Consequently, the watershed algorithm will generate a single region for each seed. The logic required to automatically place one seed point in each cell is far from trivial, and thus this method is often used with some manual interaction.

In the stochastic watershed approach, the seeded watershed is repeatedly applied to the image with randomly placed seeds [23]. The segmentation results are added together to build a map that shows the likelihood of each

Figure 1 Different endothelial cell count methods performed with NAVIS software on the same confocal image of corneal endothelium after keratoplasty. (a) Region of interest. Scale bar = 100 um. (b) Automatic cell count. Endothelial cell density: 2213 mm-2; Polymegathism: 79.0%; Pleomorphism: 33.3%. User interaction time: <1 min. (c) Polygonal area selection with manual cell count. Endothelial cell density: 687 mm-2; Polymegathism: could not be determined; Pleomorphism: could not be determined. User interaction time: 4 min. (d) Manual selection of individual cell borders with manual cell count. Endothelial cell density: 669 mm-2; Polymegathism: 21.5%; Pleomorphism: 52.4%. User interaction time: 20 min.

pixel belonging to a boundary. This map is referred to as probability density function, PDF. More salient boundaries will have a larger value in the PDF. A segmentation is obtained from this map by applying a smoothing, the H-minima transform, and finally the classical watershed algorithm.

Bernander et al. [24] proposed two modifications to the stochastic watershed algorithm to greatly increase its ability to distinguish salient boundaries and reduce over-segmentation. The segmentation method proposed under the heading "Segmentation" in the "Methods" section is based on this version of the stochastic watershed.



We used a set of 52 confocal images of corneal endothelium. These images were acquired in 23 patients with Fuchs' endothelial corneal dystrophy using a white light slit-scanning confocal microscope (Confoscan 4; Nidek Technologies Srl, Padova, Italy). All patients were examined in the first year after Descemet stripping automated endothelial keratoplasty.

Because the diseased endothelium was replaced by donor endothelium, the images varied considerably concerning the morphometric parameters, but did not contain corneal guttae, characteristic of Fuchs' endothelial corneal dystrophy. Before the examination, the eyes were

anesthetized with one drop of 0.4% oxybuprocaine (Ceban BV; Breda, The Netherlands). We used a 40x objective lens together with a z-ring adapter to minimize the eye movements during the 12 seconds of image acquisition. A coupling gel (Vidisic; Dr. Mann Pharma, Berlin, Germany) was applied between the objective lens and the corneal surface. All scans were performed using fixed device settings: full-thickness mode, 72% light intensity, and a 6 |im scan step. The images were captured with a charge-coupled device camera, producing 768x576 pixels and 8-bit integer pixel values; we manually cropped each image to remove the dark areas, leaving only the well-illuminated central part. The cropped images had a mean area of 219146 px2 (range [124740,400575]). According to the manufacturer, the lateral sampling density was 0.557 |im per pixel.

In each of these images, an expert marked a set of cells using the NAVIS software (Nidek Technologies Srl, Padova, Italy). He placed a point in the middle of each cell and drew an outline around the set of cells (polygonal frame method [25]). Thus, in the outlined region (evaluation area Aeval) is a known number of complete cells Keval. We consider the morphometric quantities that were then determined by the NAVIS software as ground truth. We also used the fully automatic segmentation mode of the NAVIS software to obtain a second set of measurements for comparison.

Data was collected in accordance with the tenets of the Declaration of Helsinki and informed consent was obtained from all participants. Approval was obtained from the Medical Ethical Committee of the Erasmus Medical Center, Rotterdam, The Netherlands (METC-2007-131). All data is freely available at http://rod-rep. com.

Frequency analysis

Below, we review published methods to derive the characteristic frequency from the frequency spectrum of the image, and propose a new one.

Radial mean of the frequency spectrum

The standard method to obtain an orientation-independent representation of the frequency content of a two-dimensional image is to compute the radial mean of the magnitude of its Fourier transform,

1 (2n FRM(f) = — J \F (f,

0) | de

where F(f, 0) is the Fourier transform of the image in polar form, i.e. f is the radial frequency and 0 the angle. To compute the radial mean, we follow this procedure: First, we define the output to be sampled at frequencies f = {0, Nf, N •••, 5} px-1, given N the width of the image in pixels (assuming a square image). Next, we determine the frequency for each sample in the frequency spectrum, and round the value to the nearest frequency in the list above. Finally, we compute the mean of all values assigned to the same frequency.

The mode of FRM is the characteristic frequency/*, and is related to the characteristic length X of the structures in the image by X = 1/f*. We determine the position and strength of all local maxima in FRM by fitting a parabola to three consecutive points around each local maximum, excluding the one atf = 0. The position of the local maximum with the highest value is the characteristic frequency. We refer to this method as modeRM.

Bucht et al. [9] take the first peak in FRM as the characteristic frequency. We refer to this method as IstRM.

Radial maximum of the frequency spectrum

Ruggeri et al. [7] proposed to compute the radial maximum of the magnitude of the Fourier transform,

FRMAxf) = max | Ff,0)|

0 e[0,2n)

They observed two peaks in this function, and stated that the first peak was related to the slow variation in the image intensity and the second peak represented the frequency of the repetitive cell pattern. Therefore, they suggested to use the second peak of FRMAX as characteristic frequency. We compute the radial maximum in an analogous way to the radial mean, and determine the

precise location of the peak by fitting a parabola. We refer to this method as 2ndRMAX.

Radial mean of the power spectrum

Foracchia and Ruggeri [6] used the mode (modePS) of the radial mean of the power spectrum,

|F(f, 0) |2 de

They also discuss using the mean (meanPS) or median (medianPS) instead. We compute the radial projection of the power spectrum in an analogous way to the radial projections above. However, Foracchia and Ruggeri [6] computed the cumulative radial mean by integrating the frequency spectrum over increasingly large circles, then computed the derivative of this cumulative distribution. Depending on the method used to determine the derivative, this results in either the same or a very similar function to that obtained by our implementation. They also applied a smoothing to Fps before computing the mode, mean or median. The size parameter of this smoothing operator was not given in their paper, but presumably has a very strong effect on the result, especially in the case of the mode. We chose to not apply any smoothing. We compute the mode of FPS in the same way as in the method modeRM.

The methods meanPS and medianPS require the application of a high-pass filter to the image. We implement this by setting the first N entries of the function FPS to zero. We indicate this parameter by appending it to the method name: meanPS 3, medianPS 5, etc.

Enhancing the ring prior to computing the radial mean

This paragraph proposes a new method, based on modeRM.

On occasions when the modeRM method fails at determining the correct cell density, we observed that the peak of |F(f, 0)| at f = 0 decays very slowly, and drowns the ring that we want to detect. We therefore developed a simple procedure to completely remove the central peak from the frequency spectrum before computing the radial mean.

The procedure applies the dilation by reconstruction [26] of the pixel at f = 0 with |F| as mask. That is, the image

Hof, 0) =

|F (0,0)| for f = 0 0 otherwise

is iteratively dilated, until stability, under the constraint that the result be smaller than | F| after every iteration:

Ht = (Ht-i e b) a|F| .

The structuring element B used has the smallest possible shape: the origin and its four direct neighbors. In practice, this operation is very fast due to a priority-queue

algorithm that addresses each pixel exactly once [27]. The result of the operation, Hx, is the central peak with its long tail, and without the ring. The function H^(f, 0) is decreasing for increasingf and constant 0. This would not be the case if a larger structuring element B was used.

The image F' = |F | — Hx contains the ring, but not the central peak. Note that this process is equivalent to a highpass filter, very specifically constructed for the image at hand. Any other high-pass filter (such as a Gaussian highpass filter) would have a parameter to be tuned, as is the case in the meanPS and medianPS methods. See Figure 2 for an example of the ring enhancement and how it affects the radial mean projection.

The method proposed here will be referred to as modeRMrec; it computes the radial mean of F', and finds the position of the maximum by fitting a parabola.

Density estimation

Iff * is the characteristic frequency determined by one of the methods above, then l = 1/f* is related to the most common cell width. The average area A of one cell can now be calculated by

A = al =

practice, one would expect a < 1.15. Note also that there is no reason to expect a > 1, as irregular patterns are likely to have lower values.

Previous publications have not considered the effect of the cell shape on the relationship between characteristic frequency and cell density. Except mentioning that "the geometrical pattern [of corneal endothelial cells] is not regularly hexagonal" [6], no explicit assumption of the cell shape is made [7,9]. Therefore, it seems that all these authors made the implicit assumption a = 1.

Assuming that the cell density of the imaged section is representative for the whole endothelium, the endothelial cell density can be calculated with

where a is a factor that depends on the shape and regularity of the cells.

Figure 3 shows two synthetic, band-limited images of 350x350 pixels with a hexagonal and a square cell pattern; in both images, each cell has a side-to-side length of 25 pixels. The same figure also shows the magnitude of the frequency spectrum of these two images. For the square pattern, the characteristic frequencyf* = 25 px—1, and the area of each cell A = 252 px2, leading to a = Af*"2 = 1. However, for the hexagonal pattern, f* = ■^325 px—1,

A = ^23252 px2, and thus a = ■ ^ 1.15. Note that the hexagonal packing is the most dense possible. Thus, in


In this section, we propose a new method to fully automatically segment the corneal endothelium. The method is based on the stochastic watershed, as described in the introduction. This algorithm has two steps. First, the seeded watershed is applied m times to the inverted input image. For each of these repetitions, uniformly distributed noise in the range [0, u] is applied to the image, and nseeds seeds are placed randomly. The results of these m segmentations are added together, yielding a PDF. In the second step, the PDF is blurred with a Gaussian filter with parameter aPDF, local minima are removed by a H-minima transform with parameter h, and then the classical watershed algorithm is applied. Optionally, the resulting segmentation can be refined to better match the location of the cell boundaries.

The procedure is described in detail below, including the methodology to determine the optimal values for all parameters, given a collection of images with a known ground truth. Refer to Figure 4 for a schematic representation of the segmentation procedure.

Figure 2 Enhancement of the central ring. (a) Magnitude of the frequency spectrum |F | for a typical image; only the central part is shown. (b) The image F; the central peak has been removed, which enhances the ring. (c) Radial mean projections for |F| (top) and F (bottom); note that the peak corresponding to the characteristic frequency is much more salient in the bottom graph. The vertical scaling in the two graphs is not the same.

Figure 3 Relation between cell shape, characteristic frequency, and cell density. Synthetic, band-limited images (350 x 350 pixels) of (a) square and (c) hexagonal cell pattern with 25 pixel side-to-side length for each cell. (b, d) Magnitude of the central region of the respective frequency spectra. Characteristic frequencies are JA = 25px-1 and J5O) ^ ^TIs px-1, respectively.

Determining the number of seed points nseeds

Ideally, «seeds is equal to the number of regions in the segmented image. We can estimate this number using the cell density estimate we obtained through frequency analysis. As will be shown in the results section, the mod-eRMrec method is the most suitable one to determine the characteristic frequency/*. Given a total image area of Aj, the number of cells expected to be found in the image is given by

«seeds = Al Sf* = - Alf

The value of «seeds is typically an underestimation of the number of regions that needs to be segmented, because some of the cells are intersected by the image border. Nonetheless, the stochastic watershed as used here is not very sensitive to the number of seeds [24,28], so it is acceptable to use this value. We may assume a = 1at this stage for the same reason.

Computing the stochastic watershed

To compute the stochastic watershed, we apply m = 100 [24] repetitions of the seeded watershed with randomly placed seeds. To place the seeds, we take a hexagonal grid with edge length

3sß «seeds f*\ 3s/3

give it a random translation and rotation, and use the vertices as seed locations [24]. The given edge length t creates a seed density equivalent to the cell density.

At each repetition of the seeded watershed, we add uniformly distributed noise in the range [0, u] to the image. This noise randomizes the order in which the pixels are processed, avoiding the reinforcement of spurious segmentation lines [24]. The results of all seeded watershed segmentations are added together to create a PDF, an image where pixels have a value between 0 and m, indicating how often each pixel was chosen to be part of a region boundary. Pixels chosen more frequently are more likely to be part of an actual cell boundary.

Due to the varying noise and the width of the boundary between cells, it is likely that the same boundary was slightly shifted in various repetitions of the seeded watershed. To unify these boundaries, we apply a smoothing filter with a Gaussian kernel. The optimal size aPDF of this kernel depends on the cell size, so we define a parameter ka = aPDF /* that is independent of the cell size.

Next, we apply the H-minima transform with a parameter h = khm/aPDF, to discard small minima. Here, kh is the parameter to be optimized, the fraction makes this value independent of m and aPDF. Finally, the classical watershed algorithm segments the image into individual cells.

Input image 100 randomly seeded watersheds PDF Blur & H-minima Watershed Seeds Refined result

Figure4 Schematic representation ofthe segmentation procedure, as detailed in the "Methods" section.

As we will see later, the computation of the PDF is not greatly affected by a wrong estimation of the value off*. In contrast, the smoothing of the PDF is rather sensitive to changes in this estimate. To have a more consistent smoothing, it is possible to estimate f* again from the PDF. We will study this improvement under the heading "Influence of estimated cell density on the segmentation algorithm" of the results section

Correcting segmentation borders

The segmentation obtained by the stochastic watershed follows ridges in the image smoothed by a Gaussian kernel with size aPDF. Such a smoothing tends to shift edges and ridges. Hence, for large aPDF, the ridges found are displaced with respect to those in the input image. This final, optional step corrects the delineation and improves subsequent cell shape measurements.

We first shrink all segmented regions by 20%, producing a marker image

S(x) =

1 lf Dc(x)+Db(x)

0 otherwise

where Dc(x) and Db(x) are the Euclidean distance to the center of mass and to the boundary of the region, respectively. Next we slightly smooth the inverted, original image with a Gaussian kernel with a = 2. Finally, we compute the seeded watershed of this slightly smoothed image using the marker image S as seeds.

By shrinking the segmented regions only slightly, we avoid large shifts of the boundaries. The purpose of this last step is not to find new boundaries, only to better align them with the boundaries in the unsmoothed image.

Evaluating the segmentation

The ground truth for the data set used consists of a point (marker) near the middle of each cell in the evaluation region («eval markers). The evaluation region was manually outlined as described in the materials section. Because it is hand-drawn, this outline cannot perfectly match the delineation by the segmentation algorithm, and thus we need to allow for small deviations.

We consider a cell correctly segmented if it contains exactly one ground-truth marker, and lies for at least 85% inside the evaluation region. We ignore any segmented regions outside of the evaluation region, as well as any region that overlaps the evaluation region by less than a quarter of the average region size. Using these definitions, we can determine «total, the total number of regions found, and «corr, the number of correctly segmented regions.

We then compute the precision p = Kcorr/«total and the recall r = Kcorr/«eval, and combine them into the F -measure

F = 2-

Following this definition, oversegmented cells decrease the precision p and under-segmented cells decrease the recall r. Because F is the harmonic mean of p and r, lowering either of them lowers F.

Parameter training

The segmentation method as presented here depends on 3 parameters: u (noise range), ka (blur size), and kh (local minima to ignore).

To determine the best values for these parameters, given a set of training images with known ground truth, we apply the algorithm to each image in the set with different values for the three parameters. We then evaluate every resulting segmentation, determining the F-measure as described in the previous section, and compute the average (mean) F for each combination of parameters. The set of parameters that yields the highest average F is the optimal set.


The health condition of the endothelial layer can be determined using three morphometric quantities: cell density, polymegathism (cell size variability) and pleo-morphism (cell shape variation). Cell density is defined as the inverse of the mean cell size; polymegathism is defined as the coefficient of variation of the cell size, and typically expressed as a percentage; and pleomorphism is defined as the percentage of cells that are hexagonal, which we compute by counting the percentage of cells with six neighbors. To avoid measuring partially imaged objects, we discard all segmented regions touching the image border.

Method evaluation

We evaluate our methods using leave-one-out cross validation [29]. Given a total of M images for both training and testing, we train the method with M — 1 images, and test on the remaining image. This process is repeated M times, such that each image has been used once for testing. The M test results are then averaged.


Density estimation

We tested the various frequency analysis methods by computing the endothelial cell density and comparing that with the manually estimated values. We used leave-one-out cross validation to determine an appropriate a for each image and each method. Figure 5 plots the relative errors.

The methods meanPS and medianPS have one parameter N (for the high-pass filter). We illustrate results with N e {2,3} for meanPS and N e {3,4,5,6,7} for medi-anPS, larger or smaller values produced worse results. With smaller values of N, these methods produced values

b 100 o

s/f / *

J? ß'

Figure 5 Relative error in estimated cell density for the various frequency analysis methods, as standard box plots. The box indicates the interquartile range, the line inside the box indicates the median, the whiskers indicate the extrema, and the dots indicate outliers.

of/* very close to zero, which obviously does not indicate any relevant frequency in the image. Thus, the high-pass filtering is indispensable for these methods.

The narrowest distributions of relative errors are obtained by the methods modeRM, modeRMrec, modePS and medianPS 5. Of these four, modePS and modeRM have the widest distributions and the most outliers. modeRMrec has a slightly narrower distribution and fewer outliers than medianPS 5, but it has one extreme outlier. Furthermore, the method medianPS 5 has a parameter that might need to change if the imaging setup changes; note that we selected this parameter using the same data set as we used to test the method.

Values of a for the modeRMrec method were in the range 0.97-0.99.

Training parameters of the segmentation algorithm

The proposed segmentation algorithm has 3 parameters for which optimal values need to be determined: the noise range u, the smoothing size ka, and the local minima depth kh. To avoid influence of the chosen frequency analysis method in the parameter training, we use a number of seeds given by the ground truth,

, «eval /.. ^

«seeds = A^-- . (12)

We then ran the segmentation algorithm on all images in the data set for all values of u between 20 and 50 in steps of 10, values of ka between 0.10 and 0.22 in steps of 0.01, and values of kh between 0.000 and 0.010 in steps of 0.001. We calculated the F-measure for each result. We then determined parameters to be used for each image i as follows: we took all the images in our set excluding image i, computed the mean F over these images for each combination of parameters, and determined the combination of parameters that yielded the largest average F. We will refer to

these values as ui, ka,i and kh,i. Testing the algorithm using these parameters is analogous to leave-one-out cross validation. For all images i, the optimal parameters were the same: u = 30, ka,,• = 0.17 and khj = 0.002. The parameter space in all these cases was smooth and convex (see Figure 6). The mean F-measure over all images was 0.907.

We also tested the method by Gavet and Pinoli [11] (using code provided by the author), as well as a much older method by Vincent and Masters [30]. We had to smooth and correct illumination of the input image to the latter method, since it was designed for in vitro images of a much lower resolution, with less noise and fewer artefacts. These methods had four and three parameters, respectively, which were optimized in the same way as

0.1 k 0.2

Figure 6 Mean F-measure for all images in the data set, with u = 30 and different values of ka and kh. Black indicates F < 0.8, white is for F = 0.9142, and occurs at ka = 0.17 and kh = 0.002 (marked by the red circle); this is the global maximum. This same combination of parameters were found to be optimal in each of the cases during leave-one-out cross validation.

described above. Vincent and Masters [30] provided the better results, with a mean F-measure of 0.821. Gavet and Pinoli [11] yielded a mean F-measure of 0.722. We did not have direct access to the segmentation results of the automatic mode in the NAVIS software, so were not able to compute F-measures, but the computed cell densities lie somewhere in between those of the methods by Vincent and Masters [30] and Gavet and Pinoli [11].

Stochasticity of the segmentation

The segmentation algorithm proposed here is stochastic, meaning that it can produce a different result every time it is run on the same image. We studied the influence of the stochasticity on the result by looking at the variation in f-measure over 12 segmentations of each image, using the optimal parameters found in the previous section. The interquartile range (IQR) of F varied from image to image, between 0 and 0.074. This indicated that for some images most segmentations were identical, whereas for other images the quality varied somewhat. In general, the images that yielded a larger average F also had a smaller variation in F, with a clear inverse correlation between the median and the range of F (p = -0.78).

Note that increasing the number of repetitions m of the seeded watershed in the stochastic algorithm will reduce the variation of the result, but not improve the results. Using the four images in the set that had the median F closest to the median over all images, we repeated the segmentation with increasing m, and observed the reduction in the range of obtained F-measures. The average F for these four images did not change significantly as m increased.

Influence of estimated cell density on the segmentation algorithm

The cell density Sf* is used to determine the number of seeds for the segmentation algorithm. A rough estimate is obtained using the frequency analysis method. The following experiment supports our claim that the segmentation algorithm is very robust against changes in the number of seeds.

We applied the segmentation algorithm to each of the images i, using parameters ui, ka,i and khi as before, but varied the number of seeds between 25% and 400% of the optimal value as determined from the ground truth. The results are shown in Table 1 (Setup 1). Note how the mean F-measure for the segmentations does not change significantly as «seeds changes by a factor of 2 in either direction. This shows that it is not problematic to assume a = 1in the segmentation algorithm.

The characteristic frequency f * also influences the size aPDF of the smoothing filter. If we do not only change the number of seeds, but change the value of f * correspondingly according to Equation 8 (simulating what

Table 1 Effect of the density estimation on the proposed segmentation algorithm

"seeds Corresponding f* Setup 1 Setup 2 Setup 3 (fraction of (fraction of mean F mean F mean F optimum) optimum)

0.25 0.50 0.883 0.206 0.874

0.50 0.71 0.907 0.766 0.908

0.75 0.87 0.926 0.913 0.923

1.00 1.00 0.914 0.914 0.915

1.50 1.22 0.928 0.892 0.925

2.00 1.41 0.918 0.791 0.910

4.00 2.00 0.831 0.465 0.807

F-measures are averaged over all 52 images. Setup 1: only nseeds changes. Setup 2: f* changes, affecting both nseeds and opdf. Setup 3: As in Setup 2, but f* is re-calculated from the PDF to determine ctpdf.

would happen if the frequency analysis produced a wrong result), the smoothing filter will change as well, affecting the segmentation result significantly. This case is reported as Setup 2 in Table 1.

One approach to reduce the influence of the frequency analysis on the segmentation algorithm is to apply modeRMrec again on the PDF produced in the stochastic watershed, and use its result to determine aPDF. We tested this procedure and reported it as Setup 3 in Table 1. Note that the average F-measures are now much less influenced by the initial estimatef*.


Next, we again applied the segmentation algorithm to each of the images in the set, using the optimal parameters ui, ka,< and kh,i obtained through leave-one-out cross validation. We estimated f* through the mod-eRMrec method, and re-estimated it from the PDF before applying the smoothing filter, as suggested in the previous section. Some example segmentation results are shown in Figure 7. From the segmentations we determine the morphometric quantities (cell density, polymegath-ism and pleomorphism) as described in the "Methods" section. We compared these measurements with the ground truth, and plotted the results in Figure 8, 9 and 10. These figures also contain the results obtained by the automatic mode in the NAVIS software. Estimated values for each image can be found in Additional file 1: Table S1. We can observe an overall improvement compared to the NAVIS software, especially for images with a cell density below ~ 1200 mm-2. Figure 8 also includes cell density estimates obtained through frequency analysis (modeRMrec), and shows that the segmentation yielded a much improved result. Note that, especially for low cell densities, images have few cells, and the slightly different choice of cells measured between manual and automatic method can have a large effect on estimated quantities.

Figure 7 Example results of the fully automatic segmentation algorithm. (a) The image used in Figure 1, and (b) the segmentation result. This is an endothelium with a low cell density, for which the NAVIS automatic method failed (compare with Figure 1b and d).The result of the proposed method produced a reasonable result, with only two cells too many (white arrows) and a few misplaced cell boundaries (black arrows). (c) One of the images for which the segmentation (d) had a perfect score. (e) Atypical high-density endothelium, for which the segmentation (f) only had two oversegmented cells. Scale bars = 50 |im.


Density estimation

We introduced the factor a to determine the endothelial cell density based on frequency analysis. The value of

a depends on the shape of the repeating pattern in the image. For a perfect hexagonal pattern a & 1.15, as illustrated in Section "Density estimation". However, for the less uniform pattern of the typical endothelium, the value

£ '¡o

S 2000

a) ■a

■a &

V> 1000 a>

2000 rr

8 i— 0

density (mm"

density (mm

3000 J>

Figure 8 Cell density. Estimated cell density (left) and error in estimated cell density (right), compared with the manual ground truth, for the segmentation method proposed here (x) and the NAVIS software in fully automatic mode (•). Note that the NAVIS software has an error that depends on the cell density. To the right, standard box plots summarize the results. The error of the modeRMrec method is included for comparison.

for a is on average closer to 1. This is presumably the reason previous authors were able to make the implicit assumption of a = 1 [6-9].

We observed that the cell structure as seen in the images can have a strong directional preference (Figure 11). In these cases, the ring in the frequency spectrum is elliptic rather than circular. The elliptic ring yields a wider peak in the one-dimensional radial projection, and leads to a mode that is a rough approximation of the average radius of the ring. In extreme cases, the peak could become so wide as to not be salient enough for detection. Nonetheless, the results on the data set used in this paper did not warrant specifically addressing this issue.

The modeRMrec method proved to be the most precise of the frequency analysis methods, and does not

have any parameters to tune. However, for one image from the data set, the modeRMrec method completely failed to detect the ring (see the point close to -100% in Figure 5). Figure 12 shows this image and the magnitude of its frequency spectrum. When applying the dilation by reconstruction algorithm, most of the ring was removed, causing the peak for the correct f * to disappear. The estimated f* was 4 times smaller than it should be. A larger data set is needed to know how common this effect is, and whether software should be able to detect or correct for this situation.

A problem with two of the frequency analysis methods we studied, meanPS and medianPS, is that they employ a high-pass filter. We implemented the high-pass filter by discarding a set number of low-frequency bins from

E 60 ■o

• • ■ • x

• •

X •. " ■ x v,' ' * —x

X *r*xX

20 30 40

polymegatism (%)

£ 40 &


.¿'ft * ■ A*r,X-

2000 density (mm-2)

3000 £

Figure 9 Polymegathism. Estimated polymegathism (left) and error in estimated polymegathism (right), compared with the manual ground truth, for the segmentation method proposed here (x) and the NAVIS software in fully automatic mode (-).To the right, standard box plots summarize the results.

Q. k_ O

E o 0) Q. T3

x •./ .

* ^ * v x

X r ,XXx .

X xX xX' Xx, 'x xx • xx ■ 'X • • "

>4 • • ■ ..

40 50 60 70 pleomorphism (%)

® -20

X. X- . X X x X •

. X X >o< X X X

2000 density (mm-2)

3000 J>

Figure 10 Pleomorphism. Estimated pleomorphism (left) and error in estimated pleomorphism (right), compared with the manual ground truth, for the segmentation method proposed here (x) and the NAVIS software in fully automatic mode (-).To the right, standard box plots summarize the results.

the radial projection. The number of bins to discard, i.e. the parameter to the high-pass filter, needs to be determined from example images. In our tests, the medianPS approach gave very good results when discarding the first five bins (medianPS 5). However, the method is very sensitive to changes in this parameter, as can be seen from the results for medianPS 4 and medianPS 6 in Figure 5. Other ways of implementing a high-pass filter would allow a more fine-grained tuning of the filter's parameter.

The NAVIS software has problems segmenting corneal endothelium with low cell density, as seen in Figures 8, 9 and 10, and previously mentioned by [4]. In Figure 8, we can observe that the automatic method strongly overestimated the density in images with a cell density below ~1200 mm-2. In fact, the software never produced a density estimate below 1000 mm-2. One possible explanation for this phenomenon could be that the method by Foracchia and Ruggeri [6], which was implemented in the NAVIS software, does not use

frequencies below 0.05 px-1, which corresponds to a cell density of 939 mm-2 ([6], Figure 5). That is, any density below 939 mm-2 yields a peak in the frequency spectrum below 0.05 px-1, and would not be found by the NAVIS software. A peak at a higher frequency, but close to this cutoff, would be affected as well. Given that 300-500 mm-2 is a critical density for the endothelium to uphold its fluid pumping function, this cutoff is puzzling.


The proposed segmentation method uses f * to determine a suitable number of seed points Kseeds for the stochastic watershed, and to determine the parameter aPDF to the smoothing applied to the PDF. The experiments in Section "Influence of estimated cell density on the segmentation algorithm" show that the influence off* on Kseeds is small for the final segmentation, but its influence on aPDF is significant. To reduce the influence of a poor frequency analysis result on the final segmentation, we recommend

to repeat the frequency analysis step on the PDF produced in the stochastic watershed. This PDF is a clean representation of the cell boundaries, in which neither the non-uniform illumination nor the other details are visible. Consequently, the frequency analysis of the PDF yields a much more precise estimate off*, and is less likely to be influenced by imaging artifacts.

The two steps used to simplify the PDF reduce the number of local minima. However, they both do this in different ways. The H-minima transform (directed by the parameter kh) directly removes local minima, without affecting other areas of the PDF. In contrast, the smoothing (directed by the parameter ka) changes the appearance of the PDF by widening and straightening the lines. As can be seen in Figure 6, these two parameters are strongly related, as changing one of them changes the optimal choice for the other. Both these operations are needed, as the global optimal combination has non-zero values for both parameters; however, the H-minima transform is the least important, as it is possible to obtain an average F close to the optimal with kh = 0.

The segmentation method is fully automatic, meaning that no user intervention is required to segment a given image into individual cells. However, there are three parameters that need to be tuned. We derived optimal values for these parameters given data from a specific instrument and a manually drawn ground truth. This procedure might need to be repeated if the method is to be applied to images from a different instrument. However, once these values have been optimized, it is not necessary to tweak these values for individual images.

Manually-drawn ground truth is often problematic when training and testing segmentation algorithms, because they are never perfect. Different experts will draw outlines differently, and will not always agree on diagnoses. However, because of the way we define and use the ground truth in our experiments, we can be quite certain that a different expert's ground truth will not produce

different F-measures: the data we used does not pose difficulties for an expert determining the extent of cells, and we allow small variations in the outline drawn around the selected cells. The only difference we have seen with repeated creation of a ground truth for an image is in which cells are contained in the evaluation region.

The segmentation algorithm proposed here is stochastic in nature. That means that it can produce a different result every time it is applied to the same image. In principle, after m = repetitions, the result is deterministic. This is, of course, impractical, and so m = 100 has been chosen in line with the recommendation by Bernander et al. [24]. As was discussed in Section "Stochasticity of the segmentation', increasing m does not improve the average F, even if it reduces its range. This implies that the probability for both extremely poor as well as extremely favorable segmentations is reduced. That is, increasing m can increase the confidence in the segmentation result. However, m is directly proportional to the time it takes the algorithm to produce its result. Hence, a compromise has to be made between precision and computation time. m = 100 has been deemed as an acceptable compromise [24]. Malmberg and Luengo Hendriks [31] have recently shown that it is possible to efficiently compute the PDF of the stochastic watershed for m = repetitions. However, this algorithm can only work for the case u = 0 (i.e. without the noise added for every repetition), which produces very poor results in this application.


It is common procedure to judge the condition of the corneal endothelium based on cell density, polymegath-ism and pleomorphism. However, it is unclear how reliable these measurements are. The main issue is that only a small fraction of the endothelium is visible in one image, and there exists no test to confirm whether the region imaged is representative of the whole endothelium [13], meaning that the results have an unknown precision.

Furthermore, some authors advised caution when comparing measurements obtained by different methods [4,5,32], indicating that different methods introduce different biases. Having a fast, fully automatic method as presented in this paper would allow recording more images, to cover a larger fraction of the endothelium, and thereby guaranteeing more precise morphometric estimates. However, this does not solve the problems with possible biases.

Of the three morphometric quantities, pleomorphism is the one that is most affected by errors in the segmentation. This is because each under- or over-segmented cell causes several neighboring cells to have a wrong polygon edge count. Furthermore, both errors of over- or under-counting polygon edges cause the pleomorphism measure to increase; these errors do not compensate for each other. Given that this measure is so unreliable, it might be worth considering alternative measures for the regularity of the endothelium, such as cell compactness or elongation, which can be computed on each cell independently.

Polymegathism has larger errors for very low density (Figure 9), as a segmentation error in that case affects the estimate more. In contrast, pleomorphism has similar errors at all densities (Figure 10).

Figure 8 shows a slight bias in the density estimation from the segmentation. One possible cause could be the way that cells are selected for this estimate: we discard cells partially in the image, which causes a bias towards smaller cells because large cells are more likely to be cut by the image border. This source of bias is carefully studied in the field of stereology, and can be solved using a counting frame [33]. We have not applied this method because the field of view in these images is rather small, with few cells imaged. Using a counting frame would further reduce the number of cells counted.

Bucht et al. [8] analyzed the frequency spectrum of a large number of simulated cell patterns. They concluded that the ratio of amplitudes of the center and the first harmonic peaks is related to the variation of cell sizes. In our notation, the center peak amplitude corresponds to the D.C. value, Frm(0) (equal to the mean grey value of the image), and the first harmonic peak amplitude corresponds to FRM (f *). Unfortunately, this relationship only holds for their synthetic images. As Bucht et al. [8] noted, larger cell variation causes the peak at f* to be smeared out, making it wider. In the idealized synthetic images they used, this translated to a lower peak as well. However, in real-world images, several other factors also affect the height of this peak as well as that of the center peak, such as the contrast, uneven illumination, noise, etc. Furthermore, to measure the width of a peak is not trivial in these data, and we have not found a reliable method to estimate the cell size variation using frequency analysis only.


In this paper we have introduced an algorithm that can automatically analyze the corneal endothelium, as imaged using in vivo confocal microscopy, with a reasonable accuracy. In some images, the algorithm's result may need further refinement. A simple and quick correction could be performed by the operator, as typically only a few cells need to be merged or split. With the right user interface, such a correction could take less than half a minute, yet yield a wealth of accurate morphometric parameters of the endothelium; currently an operator needs four minutes to simply estimate the cell density. Morphometric parameters other than density have not been practical in a clinical setting up to now, though polymegathism and pleomorphism are used in some situations. The segmentation method presented in this paper paves the way for a broader use of such parameters in the clinic. Now that it is easy to obtain a segmentation of the cells in the corneal endothelium, it will be interesting to study which morphometric parameters are important to monitor in various diseases. Measures like compactness (a ratio of the perimeter and the area of each cell) and elongation are easily quantified, and might provide better quality measures.

Additional file

Additional file 1: Values of manual and automatic estimates of endothelial parameters for each of the images in our data set. We

include the automatic estimates from the NAVIS software and from the segmentation algorithm proposed in this paper.

Competing interests

The authors declare that they have no competing interests. Authors' contributions

BS carried out the research. BS and CLLH developed the methods. BS, CLLH and BR designed the experiments. KAV conceived the study. TH collected images and performed manual measurements. All authors wrote the manuscript. All authors read and approved the final manuscript.


TH was funded through Stichting Wetenschappelijk Onderzoek Oogziekenhuis -Prof. Dr. H.J. Flieringa. No funding body had any influence in the research presented here, nor in the writing of this paper.

Author details

1 Centre for Image Analysis, Uppsala University, Box 337, 75105 Uppsala, Sweden. 2Rotterdam Ophthalmic Institute, Rotterdam,The Netherlands. 3Quantitative Imaging Group, Department of Imaging Physics, Delft University of Technology, Delft, The Netherlands. 4Rotterdam Eye Hospital, Rotterdam, The Netherlands.

Received: 9 October 2014 Accepted: 7 April 2015 Published online: 26 April 2015


1. Bourne WM, Nelson LR, Hodge DO. Central corneal endothelial cell changes over a ten-year period. Invest Ophthalmol Vis Sci. 1997;38(3): 779-82.

2. YeeRW, Geroski DH, Matsuda M, Champeau EJ, Meyer LA, Edelhauser HF. Correlation of corneal endothelial pump site density, barrier function, and morphology in wound repair. Invest Ophthalmol Vis Sci. 1985;26(9): 1191-2101.

3. Waring III, GOBourneWM, Edelhauser HF, Kenyon KR. The corneal endothelium. Normal and pathologic structure and function. Ophthalmol. 1982;89(6):531-90.

4. Jonuscheit S, Doughty MJ, Ramaesh K. In vivo confocal microscopy of the corneal endothelium: comparison of three morphometry methods after corneal transplantation. Eye. 2011;25(9):1130-7.

5. Klais CMC, Bühren J, Kohnen T. Comparison of endothelial cell count using confocal and contact specular microscopy. Ophthalmol. 2003;217(2):99-103.

6. Foracchia M, Ruggeri A. Automatic estimation of endothelium cell density in donor corneas by means of Fourier analysis. Med Biol Eng Comput. 2004;42(5):725-31.

7. Ruggeri A, Grisan E, Jaroszewski J. A new system for the automatic estimation of endothelial cell density in donor corneas. Br J Ophthalmol. 2005;89(3):306-11.

8. Bucht C, Söderberg P, Manneberg G. A model for corneal endothelial morphometry by diffraction. In: Proc. SPIE 6138, Ophthalmic Technologies XVI. SPIE; 2006. p. 61380O.

9. Bucht C, Söderberg P, Manneberg G. Fully automated corneal endothelial morphometry of images captured by clinical specular microscopy. In: Proc. SPIE 7163, Ophthalmic Technologies XIX. SPIE; 2009. p. 716315.

10. HirneissC, Schumann RG, Grüterich M, Welge-Luessen UC, KampikA, Neubauer AS. Endothelial cell density in donor corneas: a comparison of automatic software programs with manual counting. Cornea. 2007;26(1): 80-3.

11. Gavet Y, Pinoli JC. Visual perception based automatic recognition of cell mosaics in human corneal endothelium microscopy images. Image Anal Stereology. 2011;27(1):53-61.

Tinena F, Sobrevilla P, Montseny E. On quality assessment of corneal endothelium and its possibility to be used for surgical corneal transplantation. In: IEEE International Conference on Fuzzy Systems. IEEE; 2009. p. 1326-31.

Ruggeri A, Scarpa F, De Luca M, MeltendorfC, Schroeter J. A system for the automatic estimation of morphometric parameters of corneal endothelium in alizarine red-stained images. Br J Ophthalmol. 2010;94(5): 643-47.

Masters BR. Characterization ofcorneal specular endothelial photomicrographs by their Fourier transforms. In: Proc. SPIE 0938, Digital and Optical Shape Representation and Pattern Recognition. SPIE; 1988. p.


Masters BR, LeeYK, Rhodes WT. Fourier transform method for statistical evaluation ofcorneal endothelial morphology. In: Noninvasive Diagnostic Techniques in Ophthalmology. Springer; 1990. p. 122-41.

16. Masters BR, LeeYK, Rhodes WT. Fourier transform method to determine human corneal endothelial morphology. In: Proc. SPIE 1429, Holography, Interferometry, and Optical Pattern Recognition in Biomedicine. SPIE; 1991. p. 82-90.

17. Fitzke FW, Masters BR, Buckley RJ, Speedwell L. Fourier transform analysis ofhuman corneal endothelial specular photomicrographs. Exp Eye Res. 1997;65(2):205-14.

18. Doughty MJ, Spiteri M, Dilts DM. Determination of the unit size of the corneal endothelial cell mosaic from Fourier component image analysis. Tissue Cell. 1997;29(2):229-38.

19. BeucherS, Lantuejoul C. Use of watersheds in contour detection. In: International Workshop on Image Processing: Real-time Edge and Motion Detection/Estimation. France: Rennes. p. 1979.

20. Vincent L, Soille P. Watersheds in digital spaces: an efficient algorithm based on immersion simulations. IEEE Trans Pattern Anal Mach Intell. 1991;13(6):583-98.

21. Salembier P, Serra J. Flat zones filtering, connected operators, and filters by reconstruction. IEEE Trans Image Process. 1995;4(8):1153-60.

22. Meyer F, Beucher S. Morphological segmentation. J Vis Commun Image Representation. 1990;1(1):21-46.

23. Angulo J, Jeulin D. Stochastic watershed segmentation. In: Proceedings of the 8th International Symposium on Mathematical Morphology. Säo José dos Campos Brazil: Instituto Nacional de Pesquisas Espaciais; 2007. p. 265-76.

24. Bernander KB, Gustavsson K, Selig B, Sintorn IM, Luengo Hendriks CL. Improving the stochastic watershed. Pattern Recognit Lett. 2013;34(9): 993-1000.

25. Jonuscheit S, Doughty MJ, Ramaesh K. Assessment of a variable frame (polygonal) method to estimate corneal endothelial cell counts after corneal transplantation. Eye. 2012;26:803-9.

26. Goetcherian V. From binary to grey tone image processing using fuzzy logic concepts. Pattern Recognit. 1980;12(1):7-15.

27. Robinson K, Whelan PF. Efficient morphological reconstruction: a downhill filter. Pattern Recognit Lett. 2004;25(15):1759-67.

28. Selig B, Luengo Hendriks CL. Stochastic watershed—an analysis. In: Proceedings of Swedish Society for Image Analysis, SSBA. Stockholm, Sweden: KTH Royal Institute of Technology; 2012. p. 2012.

29. Duda RO, Hart PE, Stork DG. Pattern classification. John Wiley & Sons; 2012.

30. Vincent LM, Masters BR. Morphological image processing and network analysis of cornea endothelial cell images. In: Proc. SPIE 1769, Image Algebra and Morphological Image Processing III. SPIE; 1992. p. 212-226.

31. Malmberg F, Luengo Hendriks CL. An efficient algorithm for exact evaluation of stochastic watersheds. Pattern Recognit Lett. 2014;47:80-4.

32. Doughty MJ, Aakre BM. Further analysis of assessments of the coefficient ofvariation ofcorneal endothelial cell areas from specular microscopic images. Cl Exp Optom. 2008;91(5):438-46.

33. Gundersen HJG. Notes on the estimation of the numerical density of arbitrary profiles: the edge effect. J Microsc. 1977;111(2):219-23.

Submit your next manuscript to BioMed Central and take full advantage of:

• Convenient online submission

• Thorough peer review

• No space constraints or color figure charges

• Immediate publication on acceptance

• Inclusion in PubMed, CAS, Scopus and Google Scholar

• Research which is freely available for redistribution

Submit your manuscript at

BioMed Central