Itebeddine Ghorbel

Florence Rossant

Isabelle Bloch

Michel Paques

MODELING A PARALLELISM CONSTRAINT IN ACTIVE CONTOURS. APPLICATION TO THE SEGMENTATION OF EYE VESSELS AND RETINAL LAYERS

Keywords: Active contours, Parallelism, Retinal vessel segmentation, Optical Coherence Tomography, Segmentation, Quantitative Evaluation

Parametric deformable models are an important technique for image segmentation. In order to improve the robustness of the model, it may be interesting to incorporate a priori information about the shape of the objects to be segmented. In this paper, we propose to add a parallelism constraint. Such a model is relevant in many applications where elongated structures have to be detected. One main advantage of our formulation is that it only needs few parameters to be adjusted in addition to those of traditional snakes. The proposed model has been applied for the segmentation of OCT images of the retina and for the segmentation of retinal vessels. Experimental results, obtained on 25 OCT images and 30 eye fundus images, demonstrated the robustness, flexibility and large potential applicability of this new formulation. The accuracy of the method has been assessed by comparing manual segmentations, made by experts, with the automatic ones.

INTRODUCTION

The snake concept, also called "active contour model", was originally introduced by Kass et al. [START_REF] Kass | Snakes: Active contour models[END_REF] in order to detect the boundary of an object in an image. It consists of a parametric curve V , whose evolution is driven by the minimization of an energy function E, which includes terms controlling the smoothness of the curve (internal energy EInt) and terms related to information image (image energy EImage):

E(V ) = EInt(V ) + EImage(V ) + EExt(V ) (1)
The image energy attracts the curve toward the features to be detected, such as edges or lines. The external energy (EExt) can express additional contextual information. Thus, the curve is able to move dynamically from an initial position towards the desired boundary, without any further intervention. However, some problems still remain. Especially, the snake is sensitive to the initialization and it is unable to converge to the correct solution in the presence of high levels of noise. Many improvements of the original formulation have been proposed in the literature. Here, we focus on methods dedicated to the detection of several lines or contours. The authors in [START_REF] Gunn | A robust snake implementation; a dual active contour[END_REF][START_REF] Velasco | Robust parametric active contours: the sandwich snakes[END_REF] present the simultaneous use of two snakes, which are connected in a one-to-one correspondence ("sandwich snake", "dual snakes"), both looking for the same contour. For closed curves in the "sandwich snake" approach, one snake p is initially positioned inside the contour, and the other q outside. The energy equation includes a new term related to a mutual attraction force (spring force), applied in the direction of the line connecting each pair of points (pi, qi). For the dual snakes, both curves move towards a mean position, driven also by image forces and spring forces:

EExt(V ) = 1 0 λ 2 (V (s) -mean(s)) 2 ds, mean(s) = 1 2 (p(s) + q(s)) (2)
where V (s) = (x(s), y(s)) is parameterized by arc length s ∈ [0, 1] and λ controls the spring forces. However, as in the original formulation, the external force applied on one point is only related to very local image information, the greylevel or the gradient at this point. To gain in robustness, it may be very interesting to integrate higher level information, related to the shape of the object, when this a priori information is available. Following this idea, the "dual snake" model [START_REF] Gunn | A robust snake implementation; a dual active contour[END_REF] allows controlling the shape of the snake, by replacing the mean contour by a shape model in the energy term. Parallelism is also an interesting feature that has been modelled in several works [START_REF] Kerschner | Twin snakes for determining seam lines in orthoimage mosaicking[END_REF][START_REF] Laptev | Automatic extraction of roads from aerial images based on scale space and snakes[END_REF][START_REF] Al-Diri | An active contour model for segmenting and measuring retinal vessels[END_REF]. The "twin snake" algorithm [START_REF] Kerschner | Twin snakes for determining seam lines in orthoimage mosaicking[END_REF] is an extension to the traditional snake, designed for the detection of two parallel contours simultaneously, in high-resolution images. Traditional snakes cannot detect accurately close edges, and the curve often jumps from one to the other. To overcome this problem, a mutual attraction force (spring force) is introduced, which is expressed at point s as:

Eext(V (s)) = δ(s)(d(s) -d0) 2 (3)
where δ(s) is a weight factor, d(s) the actual distance between the point s and its twin partner and d0 the desired distance. This latter parameter has to be initialized beforehand, which constitutes a strong limitation of this model. Moreover, this formulation assumes that the distance is constant all the way (strict parallelism), which is generally not the case in real images.

The "ribbon snake" [START_REF] Laptev | Automatic extraction of roads from aerial images based on scale space and snakes[END_REF] is modeled by a centerline V (s) associated with a half width w(s), defining the left and right borders of the ribbon, VL(s) and VR(s) respectively. It has been initially applied for automatic road extraction from aerial images. The image energy is defined at point s as:

EImage (V (s)) = (∇I(VL(s)) -∇I(VR(s))) .n(s) (4)
where n(s) is the normal to the centerline and ∇I denotes the image gradient. This model does not allow explicitly controlling the variation of the ribbon width. Finally, the "twins ribbon" method [START_REF] Al-Diri | An active contour model for segmenting and measuring retinal vessels[END_REF] takes back the ideas developed in the "twin snake" and in the "sandwich snake" models. Four curves are evolving simultaneously: two of them inside the contour, the two others outside. The two outside contours are connected by pull forces to the corresponding inside contours (similarly to the sandwich model), while the inside contours are connected by push forces to each other (twin snake). This model is robust to high levels of noise and can locate boundaries under difficult conditions such as lack of gradient or occlusion. However, it requires a precise initialization, the prior knowledge of the ribbon width, and it is more computationally expensive than other models.

In the context of clinical studies related to ophthalmology, we need extracting thin, elongated and parallel structures that can be detected in several types of retina images, such as retinal layers in OCT images or retinal vessels in eye fundus images. These structures are difficult to segment since they present varying thicknesses and are generally weakly contrasted. Moreover, in some cases, diseases can affect the width or the grey levels of the searched structures.

In this paper, we propose a method to detect pairs of contours that are almost parallel. The advantages of our method are as follows:

(i) The distance between both contours has not to be known accurately as prerequisite. On the contrary, it is adjusted during the evolution process and can vary along the boundaries.

(ii) The parallelism constraint can be strengthened or relaxed locally by changing a weighting parameter.

(iii) All these interesting features are obtained without increasing the computational cost of the algorithm. Consequently, this model is well adapted to several applications, such as the ones mentioned above.

The remainder of the article is organized as follows: the novel Parallel Double Snake method is presented in Section 2, and experiments, including quantitative measurements, are discussed in Section 3. Let us denote by n(s) = nx(s) ny(s) the normal vector to the curve V (s) (Figure 1). Then, points of the curves V1 and V2 are given by:

PARALLEL DOUBLE SNAKE MODEL

V1,2(s) = V (s) ± b(s)n(s) (5)
Our aim is to find the centerline V (s) and the associated halfthickness b(s), so that the curves V1 and V2 match the searched image contours, which are assumed to be almost parallel. For this, we introduce a new energy defined as:

E(V, b) = EImage(V1)+EImage(V2)+EInt(V )+R(V1, V2) (6)
where EImage and EInt are defined as in "classical" snakes, in order to attract the curves V1 and V2 towards the high intensity gradient values while preserving regularity properties of V :

EImage(Vi) = 1 0 P (Vi(s))ds = - 1 0 |∇I(Vi(s))| 2 ds, i = 1, 2 (7)
EInt(V ) = 1 2 α(s) ∂V (s, t) ∂s 2 + β(s) ∂ 2 V (s, t) ∂s 2 2 ds (8)
nx ∂P ∂x V1(s) - ∂P ∂x V2(s) +ny ∂P ∂y V1(s) - ∂P ∂y V2(s) -2ϕb (s) = 0 (11) n, F (V2(s)) -F (V1(s)) -2ϕb (s) = 0 (12)
This equation allows us to adjust the local distance between V1 and V2. Let us suppose that the centerline is at the correct position, but that b(s) is locally too high (case 1) or too low (case 2). The local forces F (V2(s)) and F (V1(s)) point globally towards opposite directions (Figure 2), and the projection of their difference on the normal is negative (case 1) or positive (case 2). Thus, b(s) increases or decreases so that V1 and V2 move toward the correct position. The last term of Equation 11, related to b (s), imposes the curves V1 and V2 to be regular.

To solve these equations, we classically introduce the time variable t and discretize Equations 10 and 11, by using classical numerical approximations of derivatives (central difference in space, and backward difference in time). The global minimization is realized by updating sequentially V (s) and b(s). Thus, the algorithm can be summarized as follows:

step1: Initialize the curve V Repeat step 2: Calculate the normal vector n step 3: Update the curves V1 and V2 (Equation 5) step 4: Deform the curve V as a function of internal forces and image forces (Equation 10) step 5: Update the thickness parameter b (Equation 11)

until (Number of iterations)

Algorithm 1: Parallel Double Snake algorithm process.

The initialization of the curve V and the parameter setting are specific for each type of application.

EXPERIMENTS AND RESULTS

In this section, we present two series of experiments, performed on two eye imaging modalities. First, the Parallel Double Snake model has been applied on OCT images for the segmentation of some retinal layers. In our previous work [START_REF] Ghorbel | Automated Segmentation of Macular Layers in OCT Images and Quantitative Evaluation of Performances[END_REF], a complete segmentation method was proposed for healthy subjects. Unfortunately this method may fail in some pathological cases, where the destruction of photoreceptors induces a loss of definition of the inner nuclear layer (INL), and a loss of contrast of the hyper-reflective complex (HRC). Therefore, the INL and HRC are difficult to segment, especially the outer interface of each one. Here, we propose to first filter the original image using a nonlinear diffusion filter [START_REF] Weickert | Efficient and reliable schemes for nonlinear diffusion filtering[END_REF]. Then, the location of the inner side of the HRC (I. S. HRC) and of the INL (I. S. INL) is obtained from the k-means algorithm, followed by a Markovian regularization [START_REF] Ghorbel | Automated Segmentation of Macular Layers in OCT Images and Quantitative Evaluation of Performances[END_REF] (respectively k = 4 and k = 2 classes). Indeed, the INL appears as the cluster corresponding to the lowest mean intensity. So, the inner side of the INL boundary can be easily deduced from the labeled image. A similar process is applied for the localization of the inner side of the HRC. This provides, in both cases, an initialization of the centerline V (s), chosen as a curve parallel to the interface found previously with an initial layer radius b, which is set to a constant value, adapted to each layer to be segmented. The proposed model is applied from this initialization, in order to get an accurate segmentation of the considered layers, as illustrated in Figure 3 We compared the result of our method with the one provided by the "classical" snake algorithm in order to detect the HRC area. We used the same initial parameters, in both cases. As shown in Figure 4, we can see that the proposed method successfully segments the HRC zone (Figure 4(b)), while the "classical" snake method fails in the region outside fovea (Figure 4(a)). The parallelism constraint prevents the snake to escape towards high gradient present in the choroid (under the HRC). OCT imaging was performed on 25 images, from 8 retinitis pigmentosa subjects, with various grades of pathology evolution. A vertical and a horizontal cross-section were generally acquired for each eye. This study was conducted in accordance with French ethics regulation and all subjects gave informed consent to participate. The same parameters were used for all images. The results were evaluated by comparing the obtained segmentations with manual segmentations performed by a highly experimented medical doctor, on 5 OCT images. Intra-expert variability was assessed too. The quantitative results are summarized in Table 1, showing a very good agreement between automated and manual results (and within the range of intra-expert variability). As shown in Table 1, the intra-expert errors are on average of the same order of magnitude as the automatic/manual errors. Moreover, for the outer side of the HRC and the ONL/IS interface, the MSE between our method and the manual segmentation is lower than the intra-expert variability. The index of similarity (S) was also computed for all regions. This index is sensitive to variations in shape, size and position and a value of S > 0.7 indicates a strong agreement [START_REF] Zijdenbos | Morphometric Analysis of White Matter Lesions in MR Images: Method and Validation[END_REF]. We found that S is in all cases greater than 0.7 and that it ranges from 0.81 for the IS layer to 0.91 for the HRC layer. These evaluations confirm that the Parallel Double Snake model allows us to get accurate and robust segmentations, which are comparable to the expert's ones. In the second experiment, the proposed model has been applied on eye fundus images for the segmentation of retinal vessels in the vicinity of the optical disk. The vessels are enhanced, based on morphological operations and linear filtering, and segmented. A graph is then extracted from the skeleton of the binary image, providing a high-level representation of the vascular tree [START_REF] Rossant | A morphological approach for vessel segmentation in eye fundus images, with quantitative evaluation[END_REF]. A region of interest is then defined as the ring delimited by the circles of radius 2R and 3R (Figure 5), where R is the optical disk radius. The proposed method is then used to refine the segmentation of the branches in this region. The centerline V (s) is initialized as the skeleton branch.

The vessel radius b is estimated by computing the Euclidean distance transform to the skeleton branch and retaining the maximum value found on the vessel borders. The Parallel Double Snake model is run from this initialization to get an accurate location of the vessel boundaries, as illustrated in Figure 5. The parameters are set to α = 20, β = 0, ϕ = 50. The method was applied on a database of 30 images acquired in color or grey levels. Compared to the initial segmentation [START_REF] Rossant | A morphological approach for vessel segmentation in eye fundus images, with quantitative evaluation[END_REF], the vessels are much more precisely delineated, whatever their thicknesses and the image contrast. These results, visually assessed by the expert, show again the robustness and the accuracy of the Double Parallel Snake method. From the obtained segmentation results, measurements of arteria and vein diameters, and their combination into central arteriolar and venular equivalents (A/V ratio) can be then derived. From a clinical point of view, such measures are used for the follow up of various retinal vascular abnormalities, such as arteriolar narrowing, or arteriovenous (A/V) crossing changes.

CONCLUSION

In this paper, a new parametric deformable model has been proposed. It incorporates a priori information by adding a parallelism constraint compared to other models. The Parallel Double Snake algorithm reduces the complexity of evolution as well as the number of parameters, and it is well adapted to several applications. The quantitative evaluation shows a very good agreement and high correlation between the manual and the automated segmentations. Such segmentation results will help ophthalmologists to improve the accuracy of diagnosis of blinding retinal diseases. A more detailed evaluation will be performed in future work to assess the effectiveness of the proposed method compared to well known methods.

Parallel

  Double Snakes are a new model for the simultaneous extraction of two almost parallel curves. The proposed model consists of a centerline, V (s) = (x(s), y(s)), parameterized by s and a local half-thickness b(s). It implicitly defines two parallel curves, V1 and V2.

Fig. 1 .( 10 )

 110 Fig. 1. Parametric representation of the Parallel Double Snakes.

Fig. 2 .

 2 Fig. 2. Variations of the half-thickness b(s) given by the second Euler-Lagrange Equation (11): (a) the obtained forces on each contour (green and blue arrows), (b) the variation of half-thickness.

  . The parameters are set to b = 30 pixels, α = 20, β = 20, ϕ = 20 for the outer side of the HRC (O. S HRC), b = 10 pixels, α = 100, β = 0, ϕ = 20 for the outer side of the INL (O. S. INL). In the same way, we can then apply our method to find the outer nuclear layer/inner segments of the photoreceptors (ONL/IS) interface, which is parallel to the inner side of the HRC previously detected. The parameters are set to b = 15 pixels, α = 50, β = 10, ϕ = 20.

Fig. 3 .

 3 Fig. 3. Segmentation of some retinal layers in OCT images.

Fig. 4 .

 4 Fig. 4. A comparison between (a) "classical" snake model and (b) the Parallel Double Snake model.

Fig. 5 .

 5 Fig. 5. Segmentation of retina vessels near the optical disk in eye fundus images.

Table 1 .

 1 Mean square error (MSE) between manual and automated segmentations (in pixels).

	Layer	Automatic vs. manual Intra-expert variability
	O. S. HRC	1.06 ± 0.22	1.28 ± 0.19
	ONL/IS	1.23 ± 0.44	1.30 ± 0.18
	O. S. INL	1.35 ± 0.51	1.14 ± 0.29

Acknowledgments: This work was partially funded by a CIFRE grant from Fovea Pharmaceuticals and the French Ministry of Research (Contract CIFRE 897-2008).


Florence Rossant
email: florence.rossant@isep.fr

Itebeddine Ghorbel

Isabelle Bloch

Michel Paques

Sarah Tick

AUTOMATED SEGMENTATION OF RETINAL LAYERS IN OCT IMAGING AND DERIVED OPHTHALMIC MEASURES

Keywords: OCT retinal imaging, Automated segmentation, Quantitative evaluation

This paper proposes an automated method for the segmentation of eight retinal layers in high resolution OCT images. It has been evaluated based on comparison with manual segmentation performed by five different experts. The method has been successfully applied on a database of 72 images. Quantitative measures are then derived as an aid to ophthalmic diagnosis. A good agreement with measures derived from manual segmentation is obtained which allows us to use the proposed method for retinal variability studies.

INTRODUCTION

Recently, optical coherence tomography (OCT) has revolutionized the clinical imaging of the retina, the photosensitive tissue of the eye. The spectral domain OCT technology now reaches a 5-10 µm resolution while the high acquisition speed reduces motion artifacts, allowing reliable image averaging [START_REF] Topcon | Topcon Europe medical BV[END_REF]. These new technologies have contributed to improve the OCT image quality, and innovative segmentation approaches can now be proposed, for reaching a better accuracy and a better reliability. Then, measures derived from the segmented retinal layers (Fig. 1) will help to better understand the anatomy of the human fovea. Fig. 1. Cross-sectional OCT image of the retina, acquired with the 3D OCT-1000, Topcon [START_REF] Topcon | Topcon Europe medical BV[END_REF], and retinal layer definition.

Several retina segmentation methods have already been presented in the literature. They all start with a pre-processing step that allows reducing the speckle noise: mean filter [START_REF] Ishikawa | Macular Segmentation with Optical Coherence Tomography[END_REF], median filter [START_REF] Koozekanani | Retinal Thickness Measurements from Optical Coherence Tomography Using a Markov Boundary Model[END_REF][START_REF] Baroni | Towards Quantitative Analysis of Retinal Features in Optical Coherence Tomography[END_REF], directional filters [START_REF] Bagci | Thickness Profile of Retinal Layers by Optical Coherence Tomography Image Segmentation[END_REF] or diffusion filters [START_REF] Fernandez | Automated Detection of Retinal Layer Structures on Optical Coherence Tomography images[END_REF][START_REF] Garvin | Intraretinal Layer Segmentation of Macular Optical Coherence Tomography Images Using Optimal 3-D Graph Search[END_REF]. Then, most authors make extensive use of local boundary detection methods, as local profile analysis or local gradient analysis [START_REF] Ishikawa | Macular Segmentation with Optical Coherence Tomography[END_REF][START_REF] Koozekanani | Retinal Thickness Measurements from Optical Coherence Tomography Using a Markov Boundary Model[END_REF][START_REF] Baroni | Towards Quantitative Analysis of Retinal Features in Optical Coherence Tomography[END_REF][START_REF] Bagci | Thickness Profile of Retinal Layers by Optical Coherence Tomography Image Segmentation[END_REF][START_REF] Fernandez | Automated Detection of Retinal Layer Structures on Optical Coherence Tomography images[END_REF], which are not very robust to noise and varying contrast, and lack global consistency. Consequently, tests are required to select and label the significant edges: continuity check [START_REF] Ishikawa | Macular Segmentation with Optical Coherence Tomography[END_REF][START_REF] Koozekanani | Retinal Thickness Measurements from Optical Coherence Tomography Using a Markov Boundary Model[END_REF][START_REF] Fernandez | Automated Detection of Retinal Layer Structures on Optical Coherence Tomography images[END_REF][START_REF] Garvin | Intraretinal Layer Segmentation of Macular Optical Coherence Tomography Images Using Optimal 3-D Graph Search[END_REF], comparison with retina models [START_REF] Koozekanani | Retinal Thickness Measurements from Optical Coherence Tomography Using a Markov Boundary Model[END_REF][START_REF] Bagci | Thickness Profile of Retinal Layers by Optical Coherence Tomography Image Segmentation[END_REF], thus limiting the practical use of these techniques in a clinical environment. These approaches require a lot of parameters, empirically tuned (e.g. [START_REF] Ishikawa | Macular Segmentation with Optical Coherence Tomography[END_REF][START_REF] Garvin | Intraretinal Layer Segmentation of Macular Optical Coherence Tomography Images Using Optimal 3-D Graph Search[END_REF]) or set through a large training set [START_REF] Koozekanani | Retinal Thickness Measurements from Optical Coherence Tomography Using a Markov Boundary Model[END_REF] which is again a limitation for clinical routine. Global optimizations are sometimes used to refine the first segmentation result [START_REF] Baroni | Towards Quantitative Analysis of Retinal Features in Optical Coherence Tomography[END_REF][START_REF] Garvin | Intraretinal Layer Segmentation of Macular Optical Coherence Tomography Images Using Optimal 3-D Graph Search[END_REF]. But in all cases, not all retinal layers are detected (from two [START_REF] Koozekanani | Retinal Thickness Measurements from Optical Coherence Tomography Using a Markov Boundary Model[END_REF] to six [START_REF] Bagci | Thickness Profile of Retinal Layers by Optical Coherence Tomography Image Segmentation[END_REF][START_REF] Fernandez | Automated Detection of Retinal Layer Structures on Optical Coherence Tomography images[END_REF]), and the segmentation is not always made around the foveola [START_REF] Fernandez | Automated Detection of Retinal Layer Structures on Optical Coherence Tomography images[END_REF][START_REF] Garvin | Intraretinal Layer Segmentation of Macular Optical Coherence Tomography Images Using Optimal 3-D Graph Search[END_REF].

To overcome these limits, we propose a new method, based on more global segmentation algorithms, such as active contours, k-means and Random Markov Fields. Another originality of our approach is the modeling of the approximate parallelism between layers, based on a Kalman filter. As a result, eight retinal layers can be determined, including the inner and outer segments of the photoreceptors, which, to our knowledge, was not yet addressed in existing methods.

The remainder of the article is organized as follows: the segmentation method is presented in Section 2, experiments, including quantitative measurements, are discussed in Section 3.

SEGMENTATION

The images were acquired with the 3D OCT-1000, Topcon [START_REF] Topcon | Topcon Europe medical BV[END_REF], to provide vertical and horizontal retina cross-sections of 6mm length. The images were stored with two resolutions. We use the knowledge we have about the relative position of the retinal layers and how they appear in the OCT images, to detect and localize sequentially the retinal layers. Figure 2 presents the main steps that will be described in the next sub-sections: Fig. 2. Flowchart of the methodology.

Preprocessing

A non linear diffusion filter [START_REF] Weickert | Efficient and Reliable Schemes for Nonlinear Diffusion Filtering[END_REF] is applied to the image (Fig. 3). This filter performs better than median or Gaussian filters since edges are better preserved. The output image is then normalized between 0 and 1. In what follows, we denote by W (typically 600 or 780 pixels, depending on the image resolution) and H the image width and height. The coordinates are defined by the origin at the top left corner, the vertical x-axis and the horizontal y-axis.

Localization of the inner limiting membrane (ILM)

The contrast between the inner retina and the vitreous is high enough to detect the ILM by an edge-tracking algorithm, based on the maximization of the local mean gradient. The result is then refined and regularized by applying an active contour algorithm. The energy functional is defined by internal forces, the rigidity and the tension of the contour, and external forces, the GVF field derived from the gradient image [START_REF] Xu | Snakes, Shapes, and Gradient Vector Flow[END_REF].

Then, the foveola F(x F ,y F ) is defined as the ILM curve pixel whose x-coordinate is maximal in the central area (Fig. 4).

Detection of the hyper-reflective complex (HRC)

This area, including the ChCap, RPE and OS layers, appears as a high intensity image strip above the choroid (Fig. 4). A Gaussian 1-D filter is applied column-wise, to provide a vertically smoothed image S(x,y). The maximum response below x F corresponds to an inner pixel of the searched area. Starting from this point, the HRC median line is iteratively deduced, column by column, by looking for the maximum output of the following recursive low-pass filter:

) 1 , ( ) , ( ) 1 ( ) , ( ± + - = y x C y x S y x C α α
(1) (-from left to right, + from right to left) In this method, the S(x,y) coefficients are continuously integrated before making the decision, in order to model the continuity of the searched layer. Setting α=0.75 leads to a regular median line, despite the residual noise and the retinal blood vessel shades (Fig. 4).

The cumulative profile computed along the median line, centered on it, allows us to determine the approximate distance between the median line and the HRC edges. Thus, it is possible to initialize an active contour [START_REF] Xu | Snakes, Shapes, and Gradient Vector Flow[END_REF], which will converge accurately towards the boundaries (Fig. 4). This result will be refined in the next steps to separate the RPE+ChCap layer from the OS layer.

Segmentation of the photoreceptor segments (IS, OS)

The junction between the inner (IS) and the outer (OS) segments of the photoreceptors appears as a bright narrow band just above the RPE+ChCap layer, with a larger separation around the foveola ycoordinate. We apply a peak detector on an image area determined with respect to HRC. The maxima are detected in each column and labeled in order to form peak lines (Fig. 5a). The nearby peak lines are then merged, based on proximity criteria, and the gaps are filled using a linear interpolation (Fig. 5b). Then, the IS/OS junction is extracted by selecting the peak line that minimizes the distance to the inner side of the hyper-reflective complex. The interface between the ONL and the IS layers corresponds to a second fine line, parallel to the IS/OS junction, but less bright and more noisy. We propose to apply a Kalman filter [START_REF] Brown | The Discrete Kalman Filter, State-Space Modeling and Simulation[END_REF] to track column by column the ONL/IS boundary. The Kalman filter is defined by a state vector X, that characterizes the tracked pixels, and a theoretical model that allows predicting the evolution of X. In our application, the state vector includes the intensity value of the pixel and its distance to the IS/OS line, and both values are supposed to be constant along the searched curve. Note that the distance provides an information about the parallelism between lines, as an original feature of our approach. State noise and measurement noise variances complete the model. All the parameters are dynamically initialized from the longest peak line found just above the IS/OS junction, except the state noise which has been set for all images. The Kalman equations formalize the principle of prediction/verification, realized in 3 steps. Firstly, the state vector

( ) 1 ˆ-n n X
at the column n is predicted from the previous estimate ( )

1 1 ˆ--n n X
, according to the evolution model. Secondly, a measure is made in the image, by seeking the maximum intensity pixel around the predicted position. Thirdly, the error between the prediction and the measure is computed to correct the prediction and to provide the new state vector estimate ( ) n n X ˆ. The correction is weighted by the filter gain, which is also adjusted at each iteration. Nevertheless, when no image pixel fits the prediction, the predicted value is kept and no measure is integrated in the filter. This allows overcoming the problem of short interruptions. The process stops after several iterations without acceptable measures. The proposed method leads to an accurate localization of the ONL/IS boundary, despite the low contrast and the noise, as illustrated in Figure 6.

The foveola coordinates are then refined by minimizing the distance between the ILM and the IS/OS junction. A second Kalman filter is applied on the gradient image, to deduce the inner side of the RPE layer from the outer side of the RPE+ChCap layer, since both curves are almost parallel (Fig. 6).

This method leads to very good results, even when retinal blood vessel shades cross the layers, but fails for some blurred images, without any contrast between the RPE and the OS layers.

Alignment and clivus determination

The image is aligned against the outer side of the RPE+ChCap layer. The clivus is defined by the two highest points of the ILM found on both sides of the foveola, denoted by (x Cl ,y Cl ) and (x Cr ,y Cr ), as illustrated in Figure 7:

Inner layer segmentation (RNFL, ONL, GCL+IPL, INL)

The inner retinal layers are segmented using Markov random field (MRF) based techniques [START_REF] Geman | Stochastic Relaxation, Gibbs Distribution, and the Bayesian Restauration of Images[END_REF]. The Markovian hypothesis allows us to take into account spatial interactions between connected pixels. In our application, the classification is performed according to the Bayesian maximum a posteriori (MAP) criterion, i.e. by looking for the label configuration that maximizes the probability of the class field (the labels) conditionally to the observation field (the intensity image). This optimal configuration corresponds to a minimum state of an energy function, defined as follows.

Let us denote by ( )

i s s f P = ω
the probability distribution of the pixel intensities, conditional to the class i. Based on a statistical analysis of the inner retinal regions, we proved that the image noise is Gaussian. The probability distribution is then defined by two parameters, the mean µ i and the standard deviation σ i . We use the Potts model φ(ω s ,ω t ) to express the interactions between the 8connected sites s and t. Under these assumptions, the energy function is defined by:

( ) ( ) ( ) ( ) ∑ ∑ + ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ + - = s t s t s s f f U s s s ) , ( , 2 ln 2 2 2 ω ω ϕ β ω σ π ω σ ω μ ω (2)
The first term is related to the image data, while the second one is a regularization term. The parameter β, empirically set, weights the relative influence of both.

We first apply the k-means algorithm to initialize the labels and estimate automatically the noise parameters (µ i , σ i ) from the classification results. So, the segmentation process is fully unsupervised. Then, the energy function is minimized by running the simulated annealing (SA) algorithm. This segmentation method is applied several times, on the non pre-processed aligned image, the region of interest being adjusted each time. The contours are deduced from the labeled image, and regularized by a snake algorithm [START_REF] Xu | Snakes, Shapes, and Gradient Vector Flow[END_REF]. More details are given in what follows for each layer. The final results are shown in Figure 11 of Section 3.

RNFL: the retinal nerve fiber layer is segmented by applying the k-means algorithm on the top part of the aligned image, around the ILM membrane, with k=3 classes. The aligned image has been denoised beforehand by the non linear diffusion filter. As the contrast is very high, the classification results are very homogeneous (Fig. 8) and the MRF classification is not required. The clusters corresponding to the higher mean intensity and connected to the ILM membrane are selected. From this first result, it is easy to deduce the contours of the RNFL, on both sides of the foveola, and to regularize them by applying the snake algorithm. OPL/ONL: the k-means/MRF segmentation process is applied on the foveal region (Fig. 9), defined horizontally by the clivus, limited vertically by the ILM or the RNFL and the ONL/IS boundary, with k=3 and β=2. The ONL appears as the cluster corresponding to the lowest mean intensity (label 1) just above the ONL/IS boundary. So the OPL/ONL boundary can be easily deduced from the labeled image. It is then regularized using the active contour algorithm (Fig. 9). GCL+IPL/INL: the region of interest is now more restricted, knowing the location of the OPL/ONL interface (Fig. 10). The method is applied with k=2 classes and β=5, on the left and right parts independently, since the average contrast may be slightly different on both sides. Clusters of low mean intensity (label 1) and connected to the upper boundary of the region of interest are fused with the surrounding cluster (label 2). The location of the searched boundary is then estimated at 4 y-coordinates (y Cl , y Cl /2, y Cr /2, y Cr ) and forced at the foveola y F to pass between the ILM and the OPL/ONL interface. These 5 points are linearly interconnected to provide a first estimate of the searched boundary (Fig. 10b). Finally, the labeled image is explored column by column from top to down. Pixels corresponding to a transition between label 2 and label 1 and close to the first estimate are marked as boundary points. These points serve as initialization of a snake, for the purpose of filling in the gaps and regularizing the curve (Fig. 11). INL/OPL: a similar segmentation process is applied on the region limited by the two curves found previously.

EXPERIMENTS AND RESULTS

OCT imaging was performed on 72 images, from 25 normal subjects (either both eyes or only the right one was observed). A vertical and a horizontal cross-section were generally acquired for each eye. All subjects had a best corrected visual acuity of at least 80 points on the ETDRS (Early Trial Diabetic Retinopathy Study) charts, and a normal fundus. This study was conducted in accordance with French ethics regulation and all subjects gave informed consent to participate.

The automated detection of all layer interfaces was performed successfully on the 72 images, using the methods described in this paper. The mean processing time is around 2 minutes for each image (Matlab program running on a 2.4GHz Intel Core2 DUO CPU). Only less than 1% of the interfaces were not accurately localized. This occurs for the detection of the ONL/IS or the OS/RPE boundary in some especially blurred images. Additionally to a visual assessment by an expert, a quantitative evaluation was performed on 14 images. These images have been chosen so as to be representative of the different quality and noise levels that occur in the complete set of images. They have been segmented manually by 5 different experts, using a graphical user interface (GUI) developed with Matlab. The evaluation then consists in comparing measures (Fig. 11, Table 1) obtained from the manual segmentations and from the automated ones. Table 1. Some measures obtained from the manual and automated segmentations: maximal (D1 cl , D1 cr ), foveal (D1 f ), INL (D2), OPL (D3) and ONL (D4) thicknesses, RNLF surface (S5).

Table 1 indicates a very good agreement between the measures performed manually by the experts and those obtained from the automatic segmentation. Interestingly, computing the measures on the whole database provides similar mean values and slightly higher standard deviations (which was expected), which illustrates the representativeness of the chosen examples for the quantitative evaluation.

Finally, a distance between the manually traced curves and the ones obtained by the proposed method was calculated. The distance varies from 1.8 pixels for the RNFL/GCL+IPL interface, to 4 pixels for the OPL/ONL, which is again very accurate.

The proposed method was a basis for a preliminary study of variation of the morphology of foveal and perifoveal layers within a population of healthy subjects. We found that the profile of the fovea changed with axial length, with a relative increase of central ONL thickness which underlies a corresponding increase of foveal thickness. Independently of axial length, a large variability of the OPL/ONL complex was noted. On the opposite, the maximal retinal thickness was stable over a large range of axial length and had a low coefficient of variation. No gender-related differences were noted.

CONCLUSION

In this paper, a new automatic method for segmenting eight layers in OCT retinal images has been proposed. A quantitative evaluation has been performed on 14 images covering the usually observed variability, by comparing the results with those obtained manually by five experts. The quantitative results show a very good agreement and high correlation between the two sets of measures. Also the distance between the manually segmented interfaces and the automated ones shows the accuracy of our segmentation method, including for interfaces that were not addressed before in the literature, such as the inner and outer segments of the photoreceptors. Moreover a database of 72 images has been successfully segmented. A more detailed evaluation will be performed in future work.

Fig. 3 .

 3 Fig. 3. Pre-processed image.

Fig. 4 .

 4 Fig. 4. Detection of the ILM and the HRC. In green, the median line of the HRC detected by recursive low-pass filtering.

Fig. 5 .

 5 Fig. 5. Detection of the IS/OS junction by peak intensity analysis: peak lines before (a) and after (b) concatenation.

Fig. 6 .

 6 Fig. 6. ONL/IS and OS/RPE edges detected by Kalman filtering.

Fig. 7 .

 7 Fig. 7. Alignment and clivus determination.

Fig. 8 .

 8 Fig. 8. Segmentation of the RNFL layer by the k-means algorithm.

Fig. 9 .

 9 Fig. 9. Segmentation of the ONL layer: (a) k-means/MRF segmentation outputs (label 1 in red), (b) OPL/ONL localization.

Fig. 10 .

 10 Fig. 10. Localization of the GCL+IPL/INL boundary: (a) k-means/MRF segmentation outputs (label 1 in red, label 2 in green), (b) first estimate of the boundary GCL+IPL/INL.

Fig. 11 .

 11 Fig. 11. Segmentation of the retinal layers and derived measures.

ACKNOWLEDGMENTS

The authors would like to thank the medical team of the Clinical Investigation Center of the Quinze-Vingts hospital for contributing to the manual segmentations. This work was partially funded by a CIFRE grant from FOVEA and ANRT.


Florence Rossant
email: rossant@isep.fr

Isabelle Bloch
email: isabelle.bloch@enst.fr

OPTICAL MUSIC RECOGNITION BASED ON A FUZZY MODELING OF SYMBOL CLASSES AND MUSIC WRITING RULES

We propose an OMR method based on fuzzy modeling of the information extracted from the scanned score and of musical rules. The aim is to disambiguate the recognition hypotheses output by the individual symbol analysis process. Fuzzy modeling allows to account for imprecision in symbol detection, for typewriting variations, and for flexibility of rules. Tests conducted on a hundred of music sheets result in a global recognition rate of 98.55%, and show good performances compared to SmartScore.

INTRODUCTION

Optical Music Recognition aims at reading automatically scanned scores. Although it has been an active research field since the 80's, there are still a lot of unsolved problems. Difficulties result from segmentation problems, printing imprecision, variability of typewriting [START_REF] Blostein | A Critical Survey of Music Image Analysis[END_REF]. The need of structural approaches is now recognized. Research has been conducted in this direction at the symbol level [START_REF] Watkins | The Use of Fuzzy Graph Grammar for Recognising Noisy Two-Dimensional Images[END_REF] or at a higher level and involves mostly local graphical musical rules [START_REF] Coüasnon | Using a Grammar for a Reliable Full Score Recognition System[END_REF][START_REF] Fahmy | A graph-rewriting paradigm for discrete relaxation: application to sheet-music recognition[END_REF]. Due to the limited work at the higher level, we concentrate on the introduction of local or non-local music writing rules, either graphical or syntactic. Another contribution of this paper is to solve ambiguities by accounting for sources of imprecision and uncertainty. Existing approaches are based on statistical methods [START_REF] Stückelberg | An Architecture for Musical Score Recognition using High-Level Domain Knowledge[END_REF], or graph-rewriting rules [START_REF] Fahmy | A graph-rewriting paradigm for discrete relaxation: application to sheet-music recognition[END_REF]. Here we rely on the fuzzy sets and possibility framework, since it is well adapted to the integration of flexible rules [START_REF] Dubois | Fuzzy Sets and Systems: Theory and Applications[END_REF] and for dealing with spatial imprecision [START_REF] Bloch | Fusion of Image Information under Imprecision[END_REF].

We first present an overview of the proposed recognition method (Section 2). The next two sections [START_REF] Coüasnon | Using a Grammar for a Reliable Full Score Recognition System[END_REF][START_REF] Fahmy | A graph-rewriting paradigm for discrete relaxation: application to sheet-music recognition[END_REF] expose the individual symbol analysis step and the fuzzy modeling. This paper focuses on the fuzzy modeling of symbol classes and graphical rules, that is the main contribution with respect to our previous work. We still refer to [START_REF] Rossant | A fuzzy model for optical recognition of musical scores[END_REF] for the syntactic part dealing with metric and tonality rules. The global decision algorithm is summarized in Section 5, and we conclude with experimental results in Section 6.

SYSTEM OVERVIEW

The inputs of the program are the binary image of a monophonic music sheet scanned at a resolution of 300 dpi and some global information (clef, tonality, metric).

The processing flow (Fig. 1) is divided into three main parts. The individual analysis process performs the segmentation of the objects and the correlation with symbol models stored in a reference base [START_REF] Rossant | A Global method for music symbol recognition in typeset music sheet[END_REF]. It provides for each pattern s some recognition hypotheses. Then, the fuzzy modeling part provides for each classification hypothesis a possibility degree of membership to the class. It also introduces a fuzzy representation of the common music writing rules by expressing graphical and syntactic compatibility degrees between the symbols. Finally, the decision process evaluates bar per bar all the hypotheses combinations and chooses the most consistent one.

INDIVIDUAL SYMBOL ANALYSIS

Individual symbol analysis is mainly based on pattern matching. We compute in a small search area around each segmented object s the normalized correlation scores

y) (x, C max ) y , (x C ) s ( C k s ) y , x ( k k k s k = = (1)
These coordinates represent the position of the center of the musical symbol in the hypothesis of the class k. At most three recognition hypotheses are selected. A recognition hypothesis assigns the pattern s to a class k, if the associated correlation score is greater than the minimum threshold t m (0.3). We add also the possibility that there is no symbol ('-' in Table 1) when the highest correlation score obtained for class k is less than the decision threshold t d (k). These thresholds allow to take into account that some classes are more sensitive to typewriting variations (small value for t d (k)) or have a higher probability of false detection (large value).

Fig. 2 shows two bars, (a) and (b), and the recognition hypotheses superimposed on the original images. Table 1 indicates some of the associated correlation scores.

Table 1: Some hypotheses and correlation scores

We can see the limits of the individual analysis: the correlation scores may be very ambiguous, and the highest does not always correspond to the right hypothesis.

FUZZY MODEL

Ambiguities result from segmentation defects and variations in typewriting. This section aims at modeling the information extracted from the score during the first analysis and the knowledge about music writing. Fuzzy sets and possibility theory offer a good framework for this aim. It allows us to deal with uncertainty concerning the shape and the location of the symbols, and to express musical rules that have different degrees of flexibility.

Symbol classes

The correlation scores provide similarities between each analyzed symbol and models. Therefore we define the degree of possibility π k (s) that s belongs to class k as an increasing function of C k (s): The shape of the possibility distribution for class k (Fig. 3) is defined by two parameters: D (always 0.3), and S k that is learnt from the individual analysis results: The shape of the distribution π k does not need to be estimated very precisely [START_REF] Dubois | Fuzzy Sets and Systems: Theory and Applications[END_REF]. It has experimentally proved to be robust. The most important is that it is not a binary function (there is no crisp threshold) and that it is increasing: the higher is the correlation score, the higher is the degree of possibility.

( ) ) s ( C f ) s ( k k k = π (2)
1 ) k ( n ) k ( m ) k ( n ) k ( t S d k + + = ( 3
Table 2 shows the results obtained for the two bars. Compared to Table 1, the possibility degrees present less ambiguity. It should be noticed that the classification rank may change: see for example object 9 in bar (b), for which the possibility of a flat is now eliminated.

Graphical consistency

Musical writing rules impose some compatibility of position between two successive objects. Due to possible imprecision in the score printing and in the segmentation, as well as to symbol density variations, these rules cannot be used in a crisp way, and are rather a matter of degree.

Compatibility between accidental and note

An accidental should be placed before a note and at the same height. The possibility degree that the object s n is an accidental of class k n , and that a following nearby object s m (m>n) is a note of class k m , is a function of the compatibility degree

) s , s ( C m n k m k n p between both symbols: ( ) ( ) ( ) ( ) otherwise 0 h f and 0 l f if 0 h f l f ) s , s ( C h l h h l l k m k n p m n > >    + = ∆ ∆ ∆ α ∆ α (4)
where f l , and f h . are two functions defining the admissible values for ∆l and ∆h, the differences in horizontal and vertical position between s n and s m : This combination is a compromise between two criteria, excluding the cases where at least one is not satisfied at all. The chosen coefficients α l =0.2 and α h =0.8 express their relative importance. Using a degree between 0 and 1 instead of a crisp threshold on each criterion allows us not to discard completely an accidental which is not exactly at the theoretically expected position. For object 2 in bar (a), we obtain a compatibility degree with the next note of 1.0 in the hypothesis of a flat, and of 0.77 in the hypothesis of a sharp. These results rank correctly the flat before sharp and strengthen the right hypothesis H1.

Compatibility between notes and augmentation dot

During the first analysis process, augmentation dots are searched for in an area next to the note head [START_REF] Rossant | A Global method for music symbol recognition in typeset music sheet[END_REF]. The typical location is represented in light gray in Fig. 5a. But smaller or larger horizontal distances ∆l may occur. That is why the searching area has been extended to the black rectangle. A new compatibility degree has also been defined as a function of ∆l (Fig. 5b). It expresses that the typical distance is in the range of 1 to 2 staff spacing and that smaller or larger distances become less and less possible. Some confusion may occur if a staccato dot above a note head is also in a searching area for an augmentation dot (see Fig. 7b). Therefore we define a compatibility degree between the hypothesis that s n is an augmentation dot and the hypothesis that a following object s m is a note, as a function of the horizontal and vertical differences between them (Fig. 6, Eq. ( 5)):

( ) ( ) ( ) [ ] h f , l f Max s , s C h l k m k n p m n ∆ ∆ = (5)
According to [START_REF] Stückelberg | An Architecture for Musical Score Recognition using High-Level Domain Knowledge[END_REF], the compatibility degree is high as soon as one of both criteria is well satisfied (Fig. 7).

Compatibility between any two symbols

For any other pair of nearby symbols, we just express that they should not overlap. For that, we estimate the positions of the bounding box of each symbol (Fig. 8a), knowing the typical sizes of the two considered symbol classes. Let us denote by ∆l the difference between the right border of s n and the left border of s m , m>n. The admissible values for ∆l are defined by function f (Fig. 8b). It allows negative values for ∆l in order to deal with imprecision on the symbol size and position, and with high symbol density. In the vertical direction, we must consider different configurations. For example, if s n is above s m , we compute the difference between the bottom border of s n and the upper border of s m . The admissible values for ∆h are again defined by f. The global compatibility is expressed as the maximum of both criteria: just one has to be satisfied.

( )

( ) ( ) [ ] h f , l f Max s , s C m n k m k n p ∆ ∆ = (6)

Fusion

We have computed compatibility degrees between all pair of nearby symbols, according to their classes. In case of music sheets with high symbol density, there may be more than one nearby symbol before or after s n . Equation [START_REF] Bloch | Fusion of Image Information under Imprecision[END_REF] defines the global graphical compatibility degree for object s n classified in class k n :

( ) ( )                   = > < l n j n n k l k n p n l k j k n p n j k n p s , s C min . s , s C min s C (7)
It is the product of two terms representing respectively the compatibility with the previous nearby objects, and with the next nearby objects. The use of the product t-norm, instead of the minimum for instance, makes this rule more precise.

FUSION AND GLOBAL DECISION

The global decision is made bar per bar. We consider all possible combinations of symbol hypotheses. The ones that include at least one zero graphical compatibility degree are impossible and thus immediately discarded. For the others, we compute the average of all possibility degrees (2), graphical compatibility degrees [START_REF] Bloch | Fusion of Image Information under Imprecision[END_REF] and syntactic compatibility degrees [START_REF] Rossant | A fuzzy model for optical recognition of musical scores[END_REF]. The resulting coefficient conveys the symbol class consistency for the configuration. It is then multiplied by another possibility degree checking duration consistency [START_REF] Rossant | A fuzzy model for optical recognition of musical scores[END_REF], so that we get just one final possibility degree merging all criteria. The decision algorithm chooses finally the configuration that results in the highest final possibility degree, with priority to the ones that achieve the strict metric constraint (number of beats per bar).

RESULTS AND CONCLUSION

Tests of the proposed method have been conducted on a large database (100 music sheets, about 42500 symbols) coming from various composers and publishers. The average recognition rate is now 98.55% . It has been slightly increased with respect to the results in [START_REF] Rossant | A fuzzy model for optical recognition of musical scores[END_REF], although we have introduced more difficult scores, thanks to the improvements realized in the graphical modeling part. The latter has proved to be very efficient especially in case of high symbol density.

Comparison between the results provided by our program (Fig. 9a,10a) and those output by SmartScore 3.2 Pro Demo [START_REF]SmartScore 3.2 Pro Demo[END_REF](Fig. 9b,10b), one of the most efficient commercial software, shows that the proposed method is able to solve problems for which SmartScore fails. (b) Results with SmartScore: 9 errors For example, there are some confusions between staccato dots and duration dots with SmartScore (Fig. 9b), while we avoid this problem thanks to our graphical model. We can also point out that symbols touching each other are often not recognized with SmartScore (Fig. 9b), unlike our program that models such configurations. Symbol 9 in bar (b) is now correctly recognized as a natural, thanks to the fusion with the syntactic rules (Fig. 10a), while it is suppressed by SmartScore. Lastly, the models of symbol classes and syntactic consistency allow here to avoid the confusions made by SmartScore on the eight rest (Fig. 9b) and the half note (Fig. 10b). Examples were extracted from two music sheets. The global symbol recognition rate is 92.0% for SmartScore, 98.7%. for our program. 85.3% of quavers get the correct length with SmartScore, 99.3% with our method. These results show the interest of the proposed method, in particular the fuzzy modeling of symbol classes, that allows to deal with imprecision and variations of symbol shapes, and the fuzzy modeling of the graphical rules, expressed as compatibility degrees between symbols, that allows to deal with spatial imprecision and flexibility.

Figure 1 :

 1 Figure 1: The global processing flow

  the models M k of the reference base. Let us denote by C k (s) the highest value obtained at the (x k ,y k ) coordinates:

Figure 2 :

 2 Figure 2: Original images and recognition hypotheses (a) 2 3 8 (b) 2 9 11 H1or 0.66

Figure 3 :

 3 Figure 3: Possibility distribution of class kThe shape of the possibility distribution for class k (Fig.3) is defined by two parameters: D (always 0.3), and S k that is learnt from the individual analysis results:

  ) where n(k) is the number of objects having a correlation score with model M k larger than the threshold value t d (k), and m(k) is the average value of the scores. Because of typewriting variations, the reference model M k may not match closely the objects of class k in the processed score. In this case, m(k) is close to t d (k), and S k also close to t d (k); m(k) takes larger values otherwise, so that the possibility distribution is shifted to the right.

Figure 4 :

 4 Figure 4: Accidental / note graphical compatibilityThis combination is a compromise between two criteria, excluding the cases where at least one is not satisfied at

Figure 5 :

 5 Figure 5: (a) Search area for an augmentation dot, (b) Note / augmentation dot graphical compatibilitySome confusion may occur if a staccato dot above a note head is also in a searching area for an augmentation dot (see Fig.7b). Therefore we define a compatibility degree between the hypothesis that s n is an augmentation dot and the hypothesis that a following object s m is a note, as a function of the horizontal and vertical differences between them (Fig.6, Eq. (5)):

Figure 7 :

 7 Figure 6: Dot / note graphical compatibility

Figure 8 :

 8 Figure 8: Graphical compatibility for any two symbols

Figure 9 .

 9 Figure 9. (a) Results with our program: 0 error.(b) Results with SmartScore: 9 errors For example, there are some confusions between staccato dots and duration dots with SmartScore (Fig.9b), while we avoid this problem thanks to our graphical model. We can also point out that symbols touching each other are often not recognized with SmartScore (Fig.9b), unlike our program that models such configurations. Symbol 9 in

Figure 10 .

 10 Figure 10. (a) Results with our program: 0 error (b) Results with SmartScore: 7 errors.These results show the interest of the proposed method, in particular the fuzzy modeling of symbol classes, that allows to deal with imprecision and variations of symbol shapes, and the fuzzy modeling of the graphical rules, expressed as compatibility degrees between symbols, that allows to deal with spatial imprecision and flexibility.

Table 2 : Some hypotheses and possibility degrees

 2