1.

Introduction

Geometric deformation of the optic disc (OD) region is one of the signs of many pathologies affecting the eye. One of the most prevalent pathology that affects the OD is glaucoma. It is characterized by nerve fiber defects caused by an abnormal increase in intraocular pressure.¹ Glaucoma has a global prevalence of 3.5% in the population aged 40 years and over, the risk increasing with age.² It is the third cause of blindness after diabetic retinopathy and age-related macular degeneration.³ Evaluation of OD appearance in the fundus image (see Fig. 1, left/middle) is one of the diagnostic tests used clinically along with visual and intraocular pressure (IOP) tests.⁴ Visual field abnormalities and IOP higher than 22 mm Hg are two clinical findings of glaucoma.⁴ OD appearance can be visualized by standard color fundus photography. Changes can be very subtle, especially in early stages, and clinicians observe different abnormalities in the fundus images to assess early signs of glaucoma. These include generalized or localized enlargement of the optic cup, nerve fiber loss, asymmetry of the cups between the two eyes, narrowing of the neuroretinal rim, and hemorrhages in the OD. The enlargement of the cup is evaluated with respect to the size of the disc. Secondary signs include nasal displacement of the blood vessels, peripapillary atrophy, translucency of the neuroretinal rim, development of vessels overpass, and vessels kinks.⁴ The final diagnosis of glaucoma is supported when at least two or more findings are present, especially in the presence of other risk factors, such as age, sex, family history, identified genes associated with glaucoma, African descent, or myopia.⁴ If detected early, the progression of glaucoma can be slowed down by preventive treatments.⁵

Fig. 1

Anatomy of the OD [(a), (b)]; ISNT rule and CDR (c).

Patients with these symptoms need follow-up exams to check if there is an abnormal progression, but the early signs are difficult to detect on fundus image due to the variability of the OD shapes among patients.

The most commonly used index is the cup-to-disc ratio (CDR), which measures the ratio between the cup diameter and the OD diameter in the vertical direction (VCDR⁶; see Fig. 1). The index can also be defined as the square root of the ratio of the cup area over the disc area called the linear CDR (LCDR).⁷ These CDR measures are used in automatic glaucoma detection and achieve acceptable accuracy.^7–11 Using this indicator requires an adequate segmentation of the disc and the cup. Thus, much research effort has been devoted to segmenting the OD and cup in fundus images.^6,7,10,12,13 Other indices, such as the ISNT rule (in a normal retina the neuroretinal rim width follows the decreasing order: inferior $>$ superior $>$ nasal $>$ temporal¹⁴) and the measurement of the nasal displacement of vessels¹⁵ and vessels kinks,¹³ are used to determine the presence of glaucoma (see Fig. 1, right). While the CDR evaluates general enlargement of the cup, the ISNT rule evaluates its enlargement in the four main directions (temporal, nasal, inferior, and superior).

But the usefulness and accuracy of these indicators are debated among ophthalmologists. Indeed, an OD with a CDR as low as 0.3 can be glaucomatous, while a CDR of 0.8 does not necessarily indicate glaucoma,¹⁶ which reduces the diagnostic accuracy of this measurement. In a recent study, the ISNT rule, originally assessed on 457 normal fundus in 1988,¹⁷ was proven to have limited usefulness in assessing glaucomatous damage to the optic nerve.¹⁸

Simple geometric indicators extracted from fundus images, like the CDR and ISNT rule, are not sufficient to diagnose early geometric abnormalities, which are signs of potential glaucomateous condition of the OD.

The early signs of glaucoma are very subtle, whereas the CDR and ISNT measures describe only advanced stage of glaucoma. Furthermore, the geometric deformations of the vessels, cup, and disc are correlated because their geometric changes occur simultaneously. For example, a large disc size artificially increases the CDR, but it does not indicate pathological condition as nerve fibers need the same space across patients. In addition, when the cup is enlarged, it can induce a nasal rejection of the vessels. The CDR does not evaluate neither vessels deformation nor local deformations of the cup. The ISNT rule does not account for vessels deformation either and evaluates only simple local deformations of the cup.

One way to implicitly exploit all the geometric information contained in the fundus images is to use statistics on deformation fields by computing an atlas of the anatomy of interest. Retinal atlases found in the literature allow detection of abnormalities, such as imaging artifacts¹⁹ or exudates,²⁰ in the fundus image. The registration method used to build the atlas in these methods is an affine registration followed by a thin-plate spline warping guided by the main vessel arches.¹⁹ Thus, inside the OD, the geometric deformation is not fully described and is limited to affine registration, which restricts the statistical characterization of glaucoma that can be achieved using such an atlas. Finally, other researchers use “topological component analysis”²¹ and “principal orthogonal decomposition”²² to characterize the topological changes between one reference exam and its follow-up exams and to quantify the progression of glaucoma. These techniques were not devised to take into account the significant changes between images of different subjects and therefore are not applicable for interpatients statistical analysis.

In this article, we propose to build a statistical atlas of the OD region to analyze the variability of this part of the retina within a healthy population. We build this OD atlas from a population of healthy retinae using the most recent methods of nonlinear registration^23,24 and statistical analysis.^25,26 One major contribution of this work is a representation of the OD region (cup and rim), including the vessels, resulting in a well-defined atlas. The key to success in constructing our atlas is the use of log-demons registration combined with this local representation of the OD region. The proposed representation captures all the structures including the blood vessels that are likely to deform when the patient is subject to pathological condition. To construct this representation, we reject those vessels that cannot be registered between images and thereby form a probabilistic model containing only the statistically significant vessels. In this way, the atlas construction is not corrupted by the large residuals that would otherwise ensue from keeping all the vessels. An atlas-based shape descriptor (ASD) is then derived from the statistical atlas. This descriptor expresses the shape variability of the whole OD region. The main contribution of this work is to overcome the limitations of simple clinical indicators, such as the CDR by providing a descriptor to identify abnormalities by characterizing any significant deviation from the normal geometric variability. The paper is organized as follows. In Sec. 2, we present the methods for constructing the atlas, including registration, statistical analysis, and defining the ASD. We also introduce a specific representation of the OD region used to build the atlas. In Sec. 3, we present the results of the atlas construction from a healthy population, along with some statistics on the ASD showing its ability to characterize OD normality versus abnormalities. Section 4 concludes this article and presents ongoing and future work.

2.

Methodology

The construction of the ASD involves five major steps: the extraction of an image patch centered on the OD, the segmentation of the anatomical structures leading to a suitable representation for registration, the construction of an average model through nonlinear registration, the use of statistical tools to represent the atlas variability, and finally the construction of ASD.

Figure 2 illustrates the whole methodology to obtain the ASD, detailed in the following subsections.

Fig. 2

Methodology to form the ASD: each OD is represented as a four-cluster image. The OD region atlas is constructed by iteratively registering the four-cluster representations to the current reference.

2.2.

Four-Clusters Representation of the Optic Disc Region

We propose a representation of the OD region that contains the main structures inside that area of the retina: the neuroretinal rim, the cup, and the blood vessels. Generalized or localized enlargement of the optic cup, narrowing of the neuroretinal rim, and nasal displacement of the blood vessels are geometric deformation occurring on these structures for a pathological OD. Thus, the proposed representation takes into account these pathological signs.

2.2.1.

Optic disc and cup segmentation

Segmentation of the OD is achieved via a local K-means clustering applied to color coordinates in polar space followed by a polynomial fitting regularization step.²⁷ This segmentation method achieves a very competitive overlapping ratio of 0.9 on the MESSIDOR database.²⁸ The cup is located using an ellipse fitting method. This is achieved by locating the high gradients in polar space; the resulting ellipse is the one that maximizes gradients along its borders (see Fig. 3).

Fig. 3

OD and cup segmentation: (a) OD region patch with disc segmentation from Ref. 27; (b) cup segmentation by ellipse fitting (the fitted ellipse is marked in red and green, high polar gradients are marked in black); and (c) cup and neuroretinal rim clusters.

2.2.2.

Vessels segmentation

Ultimately, we want to be able to analyze not only the optic cup’s deformation but also the vessels’ displacement. One of the signs of pathological OD is the displacement of blood vessels. To segment the vessels inside the OD, a K-means clustering is performed in RGB space, considering that there are three classes inside the OD (see Fig. 4). The class whose average color is closest to dark red is then labeled as the vessels.

Fig. 4

Initial vessels segmentation: (a) OD region patch; (b) K-means clustering with $K=3$; and (c) selection of the vessels cluster.

However, this initial K-means method segments all the vessels that are visible in the OD region but some smaller vessels will be present in only few of the images. It would not be relevant to consider their displacement in our model. Moreover, due to topological changes that can affect these vessels, it is not guaranteed that we will be able to register them without large residuals between two OD representations, which would corrupt the atlas statistics. Therefore, to find the statistically significant vessels that the registration method will handle correctly, a vessels average model is constructed from the vessels and disc representation obtained previously [see Fig. 4(c)], using the method described in Sec. 2.3.2 below. The resulting vessels average model is used to generate a simplified vessels model that contains only the statistically significant vessels, i.e., those that are common to all the OD images in the atlas population and that we will be able to register without residuals. This simplified vessels model is obtained by thresholding the average model (see Fig. 5), so that at each nonzero locations of the simplified model, the probability to have a vessel after registration onto the vessels average model is over 75%. This value is a trade-off that is set to have the most simple vessels model while keeping the most significant vessels.

Fig. 5

Construction of simplified vessels model: (a) initial vessels segmentations via K-means; (b) vessels average model; and (c) simplified vessels model obtained by thresholding the average model.

We can observe in Fig. 5 that this simplified model is actually related to the main branches of the vessels in the OD. The final vessel segmentation is obtained by registering this simplified vessel model onto the initial K-means vessel segmentation using the log-demons algorithm,²³ presented in Sec. 2.3.1 below [see Fig. 6(b)]. We end up with an image comprising four clusters representing the background, the vessels, the cup, and the rim [Fig. 6(c)]. The OD region atlas, constructed from this simplified four-clusters representation, is noted $R$.

Fig. 6

Construction of four-cluster representation: (a) initial vessels segmentation via K-means; (b) simplified vessels model registered onto K-means vessels segmentation; and (c) final four-cluster representation.

3.

Experiments and Discussion

We first validated the construction of the statistical atlas by comparing the registration residuals resulting from using different representations and different registration methods. Then, the OD statistical atlas was further validated by analyzing the average model and its variance field and principal variation modes. Finally, we examined whether the ASD was able to distinguish between healthy cases and pathological conditions, such as glaucomateous or prone to develop abnormalities ODs.

3.1.

Construction of the Optic Disc Atlas

We selected 60 healthy fundus images from the MESSIDOR database,²⁸ each having a wholly visible OD and cup and therefore a high confidence in the segmentation accuracy, to construct the statistical atlas. These images were captured using a Topcon TRC NW6 nonmydriatic retinal camera with a 45-deg field of view. The image sizes were $1440\times \mathrm{960,2240}\times 1488$ or $2304\times 1536\text{pixels}$.

A glaucoma specialist was asked to eventually correct the automatic results of the segmentation by manually adjusting the optic cup and disc boundaries. This step was required to avoid any bias in the atlas construction method and the proposed descriptor.

To validate the construction of the OD atlas, including the choice of representation and the registration method, we quantitatively evaluated the residuals from the registration phase. The residual is defined as follows, with $R_{\infty }^{\mathrm{ref}}$ the average model:

Eq. (21)

$${\text{Residual}}^i={\Vert R_{\infty }^{\mathrm{ref}}-R^i\circ \exp (v^i)\Vert }^2.$$

The regularization parameters of the log-demons algorithm, ${\sigma }_{\mathrm{fluid}}$ and ${\sigma }_{\mathrm{diff}}$, were both set to 1.5 to constrain the smoothness of the transformation. These values were set by minimizing the registration residuals on a random set of OD representation from MESSIDOR dataset.

We compare the residuals from the atlas construction to evaluate the contribution of the log-demons registration method over simpler affine registration used in this work and in Lee et al.¹⁹ Likewise, the proposed four-clusters representation of the OD region was compared to the grayscale representation, obtained by averaging the R,G, and B intensities. The impact of a large registration residual is that when we average the registered image intensities to form the average model, the edges will be blurry and hence not well preserved. Furthermore, the variability will not be correctly calculated as the ODs will not be registered correctly onto the average model.

Figure 7 shows how the residuals evolve as a function of the number of iterations in the average model construction process. Clearly, using the four-cluster representation is better than the grayscale representation and log-demons registration performs better than affine registration. We achieve the best average residual using log-demons and the four-cluster representation. As expected, the affine registration residual does not decrease after the first iteration because the average model is directly calculated from the affine transformation. From these results, we can first see that the smooth diffeomorphic transformation, which has more degrees of freedom than an affine transformation, and reduces registration errors. Second, the four-cluster representation is less complex to register for the log-demons than the grayscale representation. The grayscale representation leads to incorrect registrations because the gradients and the disc and cup edges that contribute to the velocity update are noisy and the small vessels are difficult to register. The strength of our representation is that we extract only the statistically significant vessels inside the OD region and the registration is not corrupted by vessels that are not present in the whole population. For example, the cilioretinal artery coming from the temporal side of the OD is present only in 30% of retinae.¹

Fig. 7

(a) Registration residual during average model construction for different registration methods and OD region representations. (b) Norm of the average of velocity fields $\Vert \frac{1}{N}\sum _{i=0}^N{v}_k^i\Vert$ during average model construction.

We evaluate that the resulting average model is well defined by calculating the norm of the average of velocity fields $\Vert \frac{1}{N}\sum _{i=0}^N{v}_k^i\Vert$ and observe how this quantity evolves throughout the minimization (see Fig. 7). The norm of the average of velocity fields converges after four iterations and its final value corresponds to a change of less than 0.5 pixels on average.

Next, we conducted a test to determine the influence of the choice of initial reference image $R_{k=0}^{\mathrm{ref}}$; five different images, chosen randomly, were tested. We can observe, in Fig. 8, the initial references, their four-cluster representations, and the corresponding final average models. The resulting average models are very similar. This is confirmed by the fact that after four iterations, the average and standard deviation of the root-mean square error of intensities between the five average models is very low (see Fig. 9).

Fig. 8

Influence of initial reference selection: (a) initial image chosen as reference; (b) corresponding four-cluster representation $R_{k=0}^{\mathrm{ref}}$; (c) resulting final average model $R_{\infty }^{\mathrm{ref}}$; and (d) corresponding average model $I_{\infty }^{\mathrm{ref}}$ in RGB space.

Fig. 9

Root mean square error between the five average models constructed from five different initial reference images.

These results indicate that the average model construction behaves well with the representation we have chosen and that it is robust with respect to the choice of the first reference image.

3.2.

Optic Disc Atlas Analysis

We can now present the results of the statistical analysis of the average model and its variability measures (local variance and principal components). First, we analyzed the average model qualitatively. In Fig. 10, we can observe that the average model constructed with only affine registration is very blurry. By contrast, using the log-demons registration results in well-preserved edges of the principal anatomical structures, namely the cup, the disc, and the vessels, while keeping a smooth and invertible transformation. The variance field is shown in Fig. 10(c). This map emphasizes that the variability is high for the vessels inside the OD, whereas it is moderate on the boundaries of the OD and cup and low inside the neuroretinal rim.

Fig. 10

Atlas analysis: (a) average model calculated using only affine registration; (b) average model using log-demons algorithm; and (c) variance map of model in (b).

The average OD model is slightly elliptical with the vertical axis as the major axis, which is consistent with normal anatomy.⁴ The LCDR of the average model is 0.49 and the VCDR is 0.51, which are close to the average of the CDR measures in the atlas population (0.51 for both the LCDR and VCDR). The cup’s location is slightly temporal with respect to the center of the OD and the ISNT rule ($\mathrm{I}>\mathrm{S}>\mathrm{N}>\mathrm{T}$) is respected (see Fig. 11).

Fig. 11

Anatomical validity of average model: VCDR measure and conformity with ISNT rule (lengths indicated in pixels.)

The PCA applied to the OD atlas is another way to visualize the variability. Rather than expressing the variability locally, the eigenvectors represent velocity fields and therefore global displacement. Each eigenvector explains a proportion of the total variability equal to its eigenvalue divided by the sum of all eigenvalues.

We observe the first eigenvector, explaining 29% of the variability, in Fig. 12. This figure shows the principal deformation covering the range of variability from $-2\sigma$ to $2\sigma$ and the corresponding closest OD image from the atlas projected onto the eigenvector. The first mode corresponds to a simultaneous enlargement of the disc and the cup. The size of the cup is more variable horizontally, thus a vertical enlargement of the cup would more likely be abnormal than a horizontal one.

Fig. 12

First (a) and second (b) mode of variation: (first row) $\exp ({(VE)}_0)$ applied to RGB average model; (second row) corresponding four-clusters average model; (last two rows) closest element of the database to the corresponding $\exp (\alpha {(VE)}_{\mathrm{1,2}})$, with $\alpha \in [-2;-1;0;1;2]$.

The second eigenvector, explaining 22% of the variability, corresponds to a nasal displacement of the vessels as they are pushed by the enlargement of the cup (see Fig. 12 left). Interestingly, this mode of variation characterizes a type of deformation that is one of the less-specific clinical signs of glaucoma.⁴

The two following eigenvectors, explaining, respectively, 10% and 9% of the variability, are close to small rotations with different centers (see Fig. 13). These two eigenvectors are certainly related to the eye placement during the image acquisition. The fifth eigenvector (accounting for 6% of the variability) is a vertical stretching of the OD combined with a vertical stretching of the cup (see Fig. 13). These first five eigenvectors explain 76% of the variability in the atlas population.

Fig. 13

Third, fourth, and fifth modes of variation applied to RGB average model [$\exp ({(VE)}_i)$] with $i\in [2;3;4]$).

The sixth eigenvector, shown in Fig. 14, explains 4% of the variability and corresponds to a horizontal enlargement of the cup while the disc boundary does not move. The seventh eigenvector, illustrated in Fig. 14, explains 3% of the variability and corresponds mainly to a vertical enlargement of the cup. It is important to note that these two modes characterize specific deformations that can help to detect glaucoma.

Fig. 14

(a) Seventh mode of variation applied to RGB and four-clusters average model $[\exp ({(VE)}_6)]$. (b) Sixth mode of variation applied to RGB and four-clusters average model $[\exp ({(VE)}_5)]$.

The first 10 eigenvectors explain 90% of the variability, whereas the first 20 account for 96.4% of the variability (see Fig. 15).

Fig. 15

Cumulative percentage of variability explained by the eigenvectors.

3.3.

Assessing Abnormal Conditions Using the Atlas-Based Shape Descriptor

We want to evaluate the ability of the ASD to capture early changes in the geometry of OD and characterize pathological cases. Experiments are conducted on two datasets. The first one contains 16 OD-centered fundus images from Ibn-El-Haythem, Algers, Algeria. These fundus images have been identified as pathological by two ophthalmologists. Pathological case present symptomatic signs and needs, according to the opthalmologists, close follow-up examination to control the potential evolution of the pathology, which is important to detect as early as possible. For the second dataset, 255 fundus images were acquired from the CARA screening platform (Diagnos Inc., Montreal, Canada). These images were captured using digital retina cameras from different manufacturers (Centervue, Canon, Topcon and Zeiss) with a 45-deg field of view and are of size $1620\times 1444$, $2196\times 1958$, $2592\times 1944$, or $3872\times 2592\text{pixels}$. Among the 255 images, 126 fundus images were identified as glaucomatous cases by two ophthalmologists. The other 129 fundus images were not suspected to have glaucoma and considered to be healthy.

The box plots of the distribution of several ASD components are shown in Fig. 16 for the four resulting populations (atlas, healthy, pathological, and glaucomateous). We can see that for all of these measures, the atlas and healthy groups’ distributions follow each other closely, whereas the pathological and glaucomateous group is always shifted with respect to the other two. The results of Mann–Whitney test show strong difference between healthy and pathological or glaucoma population with most of the time a $p$-value less than 0.0001 (see Table 1).

Fig. 16

(a) Distribution of measures between the four image sets, for $p_{32}$, $p_{47}$, and $p_{59}$. $P$-value of Mann–Whitney test of ***, $<0.0001$; **, $<0.001$; *, $<0.01$. (b) Distribution of measures between the four image sets, for $d_{\max }$, $p_5$, and $p_{27}$. $P$-value of Mann–Whitney test of ***, $<0.0001$; **, $<0.001$; *, $<0.01$.

Table 1

Significance difference between healthy and abnormal OD with CDR and ASD.

	Healthy	Pathological		Glaucoma
	Mean±Std	Mean±Std	p-valuea	Mean±Std	p-valuea
LCDR	$0.48\pm 0.05$	$0.54\pm 0.05$	$<0.0001$	$0.59\pm 0.07$	$<0.0001$
Max deviation	$0.22\pm 0.05$	$0.31\pm 0.17$	$<0.01$	$0.41\pm 0.19$	$<0.0001$
$p_5$	$-0.07\pm 0.12$	$0.05\pm 0.13$	$<0.001$	$0.08\pm 0.16$	$<0.0001$
$p_{27}$	$-0.06\pm 0.20$	$0.09\pm 0.18$	$<0.0001$	$0.11\pm 0.23$	$<0.0001$
$p_{32}$	$-0.05\pm 0.22$	$0.09\pm 0.20$	$<0.01$	$0.17\pm 0.24$	$<0.0001$
$p_{47}$	$0.10\pm 0.09$	$-0.23\pm 0.24$	$<0.001$	$-0.37\pm 0.38$	$<0.0001$
$p_{59}$	$0.05\pm 0.19$	$-0.28\pm 0.21$	$<0.0001$	$-0.49\pm 0.47$	$<0.0001$

^aFor Mann–Whitney test.

Scatter plot on the two variables $p_{47}$ and $p_{59}$ show that these two components of the ASD are uncorrelated (see Fig. 17).

Fig. 17

(a) Scatterplot of atlas, healthy, and pathological along two components of the ASD. (b) Scatterplot of atlas, healthy, pathological, and glaucomateous along two components of the ASD.

These components of the ASD correspond to the projection on a specific eigenvector of the atlas. We have already seen that $p_5$ looks like an enlargement of the cup horizontally. In Fig. 18, we show the projection corresponding to $p_{47}$ and $p_{59}$. Interestingly, this seems visually to correspond to a vertical and horizontal displacement of the cup, respectively.

Fig. 18

$p_{47}$ and $p_{59}$ corresponding mode of variation applied to the RGB average model.

We can observe in Fig. 19 that the LCDR is uncorrelated with $p_{47}$. The LCDR is the square root of the ratio of the cup area over the disc area. It shows that the ASD has at least one component with discriminative information not present in LCDR.

Fig. 19

Scatterplot showing that CDR is uncorrelated with the measure extracted from the ASD.

In summary, several components of the ASD show significant differences between glaucomatous and healthy ODs while the healthy ODs lie inside the expected variability of the atlas. While some components of the ASD have low discriminating power to characterize glaucoma by themselves, most of the components of the ASD are uncorrelated with each other, which is very promising for the prospect of combining several ASD components to improve their discriminating power and specialize them for a specific pathology affected by specific geometric deformations.

4.

Discussion and Conclusion

We have proposed a method to construct a statistical ASD to represent the geometric deformation of the OD region in fundus images. It is the first time that a validated atlas of the OD region is constructed. The atlas construction methods have been validated, and the resulting components of the ASD have been analyzed. Finally, we have showed that these components were able to characterize discriminative deformation and therefore distinguish between healthy ODs and pathological ODs.

The CDR is a good measure for comparative exam but is not able to deal with the variability of OD shape among patients, for example, OD size. This is the strength of our new atlas approach that through the shape descriptor, obtained from a PCA of the atlas velocity fields, we are able to separate the OD size variability among healthy patients from the pathological variability. Therefore, the ASD is able to characterize abnormalities from only one fundus image. We have shown that the ASD is able to distinguish normality as well as abnormalities including glaucoma. Many components of the ASD are uncorrelated with CDR measures, which confirm that all the available geometric information is not being exploited when using only the CDR. Additionally, the CDR is specific to describing glaucoma, whereas the ASD can highlight any significant deviation from the average mode and can reveal more signs than the CDR.

One big advantage of the atlas approach is the easiness to interpret results and what have been learned in the atlas. We have already seen that each element of the ASD has a geometric meaning and that we can visualize the modes as a geometric deformation and therefore explain on which deformation mode a pathological OD is outside the healthy variability. Furthermore, the local deviation map $d_{\mathrm{new}}$ is very useful to track the local deformations that are outside the atlas variability (see Fig. 20). The ASD is able to assist the clinicians in the assessment of pathological ODs and supplement the use of clinical measures, such as the CDR. The local deviation map highlights the local deformation and can assist the clinicians to detect these early local deformations.

Fig. 20

Three pathological ODs with column (a) the original cropped OD; column (b) the OD with the velocity field toward the average model applied; and column (c) the local deviation map superimposed onto the original cropped OD. The local abnormalities lie on the temporal side of the cup and also on the vessels on the nasal side (row 1), on the superior side of the cup and on the nasal rejection of the inferior vessels (row 2), and on the inferior side of the cup (row 3).

To further assist the clinicians, the statistical atlas will help the clinicians in early glaucoma detection and stratification with the computation of a new glaucoma likelihood score calculated from the ASD. More generally, this statistical atlas will help to study the correlation of the ASD components with clinical assessment in pathologic cases. We will be able to analyze and quantify the morphological variability of OD region in healthy populations and better understand ethnic differences. As well, we will be able to study the asymmetry of the OD region abnormalities between right and left eyes.

In this work, the statistical atlas is constructed with a healthy OD region to model the variability of a healthy population regardless of other demographic data. With added information such as age, we would be able to construct a longitudinal statistical atlas to analyze the variability due to aging and those related to a pathology. Furthermore, the atlas construction framework could be extended to characterize 3-D deformations with optical coherence tomography exams of the OD. The 3-D shape of the OD carries useful information and the variability could be studied using the same framework by adapting the registration procedure to 3-D velocity fields.

However, the number of fundus images chosen to construct the atlas is currently limited. We chose to validate each segmentation by an ophthalmologist so that the results of the study are not biased by any segmentation errors. The similarity of the ASD components between the atlas population and the other 129 healthy population, however, indicates that the current atlas represents well an unknown healthy population. In addition, our method for constructing the average model is incremental so we can easily add new ODs to the average model. An easy way is to begin with the current average model as the representation and iterate toward the new average model. So, adding new patients to the atlas is fast, which is a real advantage of an atlas approach.

In this study, only one ophthalmologist has validated the segmentation. Thus, the interexpert variability in OD region segmentation should be studied. For future improvements, we could construct the atlas from the average segmentation of multiple experts when available.

As for now, the proposed method is semiautomatic. Full automatic cup segmentation could be achieved with the recent development of deep learning techniques. However, lack of annotated data requires realistic data augmentation, especially for medical images. Our proposed atlas would be able to procure realistic deformation of the OD region and therefore improve automatic segmentation of the OD region.

Ongoing work focus on gathering more pathological data and then proposing a method to construct specific scores derived from the ASD, specialized to detect different pathologies.