P Brown

D W Cunningham

V Interrante

J Mccormack

T Hurtut

Y Gousseau

F Cheriet

F Schmitt

Pictorial analysis of line-drawings

Keywords: (according to ACM CCS): I, 2, 10 [Vision and Scene Understanding]: Perceptual reasoning, shape

We present in this paper an approach to the analysis of the pictorial content of artistic line-drawings. The pictorial content is the combination of the stylistic content and of the visual features of the represented subject. This paper focuses on the pictorial content hold by line strokes in line-drawings. To this aim, we propose a parameterfree method to detect the hierarchical set of stroke contours. This structure allows to estimate the radius of the drawing tool that has been used. This information then efficiently tunes several methods to extract strokes curvature information, endpoints, stroke junctions and corners. The efficiency of the proposed methods is illustrated with several experiments.

Introduction

In this paper, we call pictorial content the combination of the artistic style and of the visual features of the represented subject contained in an artwork. This content contributes to the visual impression delivered to a human observer. "In an abstract painting, ideas, emotions, and visual sensations are communicated solely through lines, shapes, colors, and textures that have no representational significance" [START_REF] Owen | Painting. Encyclopaedia Britannica[END_REF]. The psychological effect delivered by the pictorial content, or pictorial effect, is one part of the socalled aesthetic effect [START_REF] Silverman R | Learning about Art: A Practical Approach[END_REF]. However, since computational methods do not succeed in describing the representational content, the aesthetic effect cannot be fully extracted. The semantic gap [START_REF] Sebe | Contentbased multimedia information retrieval: State of the art and challenges[END_REF] is still a topical research problem in computer vision. This gap is even wider in the visual arts field since pictures are often not perfectly realistic. We therefore limit our study to the pictorial level of analysis. The analysis of the pictorial content is useful in many applications for identifying and retrieving similar artworks in cultural heritage databases [START_REF] Stanchev | Some Issues in the Art Image Database Systems[END_REF]. It is also useful for non-photorealistic rendering (NPR) of scences in computer graphics that produces the same pictorial effect as some real artworks [JEGPO02, [START_REF] Bae | Statistical analysis and transfer of pictorial styles[END_REF][START_REF] Grabli | Programmable style for NPR line drawing[END_REF].

The pictorial content is partially related to the artistic style. The style detection problem has been tackled in many references, see e.g. [START_REF] Van Den Herik H | Discovering the visual signature of painters[END_REF][START_REF] Wang Z | Studying digital imagery of ancient paintings by mixtures of stochastic models[END_REF][START_REF] Yan | Indexing and Retrieving Oil Paintings Using Style Information[END_REF]. These ap-proaches are usually based on a list of low level features and a few definitions of some modern art movements such as cubism, impressionism, etc. The main problem is that artistic style inherits from many definitions. According to the American Dictionary [P * 00], style is "the combination of distinctive features of artistic expression, execution or performance characterizing a particular person, group, school or era". Unfortunately style depiction is often based on the same visual effects as subject depiction. Style recognition is thus a very difficult task, requiring the knowledge of numerous art historians and experts, which still cannot be fully carried out with computer vision techniques. Line-drawings have long been considered as a very important study in artist preliminary works [START_REF] Kandinsky | Point and Line to Plane[END_REF][START_REF] Paul | cours du Bauhaus[END_REF]. It is also often a standard format for illustration in print publications using black ink on white paper. Line content in artworks have been addressed in several perceptual studies. In its inspiring study [START_REF] Willats | Art and representation[END_REF], Willats proposed several depiction systems that let one person describe the artistic rules chosen by an artist. These depiction systems are well adapted to NPR line-drawings [START_REF] Grabli | Programmable style for NPR line drawing[END_REF]. Leyton studies the visual tension induced by the geometrical shapes in a painting [START_REF] Leyton | The structure of paintings[END_REF]. Leyton theory proposes to see shapes as memory storages. Practically, his theory is based on curvature extrema information and local neighborhood of these extrema. This approach is closely related to cognitive studies such as the seminal work of Attneave for instance [START_REF] Attneave | Some informational aspects of visual perception[END_REF].

Studying the line content in artworks has been rarely investigated in computer vision. Berezhnoy et al. [START_REF] Berezhnoy | Authentic: Computerized Brushstroke Analysis[END_REF] proposed a method using polynomial interpolation of painting brush-strokes on Van Gogh artworks. This approach is adapted to very small strokes such as Van Gogh's. Li and Wang studied ancient black and white ink drawings using wavelets and hidden Markov models [START_REF] Wang Z | Studying digital imagery of ancient paintings by mixtures of stochastic models[END_REF]. Their method is yet close to texture analysis in paintings. Onkarappa and Guru studied the spatial mutual arrangement of strokes in line-drawing images to achieve similarity retrieval [START_REF] Onkarappa | Modified 9DLT Matrix for Similarity Retrieval of Line-Drawing Images[END_REF]. This approach aims at describing the artistic composition of an image.

In this paper, we first propose in Section 2 a parameterfree approach to detect closed contours in line-drawings. Next, we introduce in Section 3 several methods for extracting and characterizing some important visual features: curvature, corners, junctions, endpoints. In Section 4, some similarity experiments on a database of artistic line-drawings will be used to demonstrate the efficiency of the features.

Detection of closed stroke contours

We thereafter call stroke cluster any set of connected strokes. In this section, we explain how to associate closed contours to the stroke clusters present in a line-drawing. This is a crucial step since the geometrical analysis introduced in the next section is based on these closed contours. Among the numerous segmentation approaches proposed in the literature, we choose to rely on level lines and draw our inspiration from the method presented by Desolneux et al. in [START_REF] Desolneux | Edge detection by helmohltz principle[END_REF]. Level lines hold several interesting properties. They are represented by closed curves in a hierarchical structure, they do not require any initialization, and the selection method proposed in [START_REF] Desolneux | Edge detection by helmohltz principle[END_REF] relies on only one parameter which is fixed for every image.

Level lines are defined as the connected components of the topological boundary of level sets. There are both positive level lines (corresponding to upper level sets) and negative ones (corresponding to lower level sets). These lines are organized in an inclusion tree, called the topographic map, that can be efficiently computed on digital images, as detailed in [START_REF] Monasse | Fast computation of a contrast-invariant image representation[END_REF]. Level lines provide a fairly complete but largely redundant representation of the geometry of an image. Desolneux et al. have proposed an image segmentation method relying on a filtering of the topographic. A meaningfulness measure called number of false alarms (NFA) is affected to each line. This measure depends on both the length of the line and its contrast. According to the NFA measure, selecting meaningful level lines practically yields lines that are well localized along contours but are still largely redundant. In [START_REF] Desolneux | Edge detection by helmohltz principle[END_REF], it is suggested to further filter the meaningful lines by using a maximality principle. The goal of this step is to associate exactly one line to each geometrical entity present in the image. The maximality principle relies on the fact that the set of meaningful level lines inherits the hierarchical structure of the topographic map (i.e. two meaningful level lines are either disjoint or linked by inclusion). For each sequence of meaningful level lines that have each exactly one child and are of the same type (positive or negative), only the more meaningful is kept. While giving good results on natural images, this maximality principle is not well adapted to the segmentation of line-drawings. This is mainly due to contrast variations along the strokes.

We therefore propose an alternative way to filter the tree of meaningful level lines that is specifically suited to linedrawings and aims at associating a unique contour to each stroke. The ideal structure we are aiming at is explained in Figure 1 for some manually selected level lines on a linedrawing. In order to achieve such a segmentation, we take advantage of the fact that line-drawings have a constrained structure, and in particular that a stroke is made of a single shape inside which there is no contrast inversion. We call Maximal Monotone Tree (MMT) a subtree of the set of meaningful level lines than contain only positive or only negative nodes and which cannot be included in another monotone subtree. For each MMT of the tree of meaningful level lines, we first look for the most meaningful level line L according to the NFA measure. We then remove all level lines from the MMT which are children or parents of L. These two steps are repeated until the MMT is empty. This principle allows us to keep only the most meaningful level lines and it insures in most situations that each contour of a stroke will be represented by a unique level line, without having to tune any parameter. The contours obtained using this MMT-principle on the tree of meaningful level lines in Figure 1a is exactly Figure 1b. More examples of stroke contours detection are to be seen in Figure 6.

a) b) c) Figure 1: A simple line drawing of two leaves (left). Four level lines are manually chosen to best represent the strokes contours (middle). On the tree structure of these four level lines (right), a positive or negative level line is represented by a white or black node respectively. The unique positive level line b describes the inner contour of the biggest leaf.

Finally we group these selected meaningful level lines to model the contours of every stroke cluster. These groups of c The Eurographics Association 2008. level lines, thereafter called coalescence sets, are obtained by grouping every retrieved negative level line with its positive children. More details can be found in [START_REF] Hurtut | Pictorial content analysis of line-drawings using geometrical shape information[END_REF].

Pictorial content analysis

This section proposes several methods to analyze the coalescence sets introduced in the previous section. The estimation of the tool radius (Section 3.1) tunes several geometrical parameters in the following steps. The proposed description of line-drawing is based on the curvature distribution (Section 3.2) and on the detection of several visual cues such as stroke junctions, corners and endpoints (Section 3.3). Eleven scalar features are computed based on these cues (Section 3.4).

Tool radius estimation

We make the following hypothesis. Each stroke cluster Φ is done with tools having the same radius R Φ . This may sound restrictive since for example calligraphic artworks are discarded, but line-drawings are often made all over with one single tool. This hypothesis let us easily estimate the tool radius. Following a constant width stroke model made with a circular tool, for each stroke cluster Φ associated with a coalescence set made of n level lines L i , the tool radius R Φ is given by:

R Φ ≈ sur f ace perimeter = -∑ n i=1 sign(L i )S i ∑ n i=1 P i , (1)
where S i (resp. P i ) is the polygonal surface (resp. perimeter) of the level line L i .

Curvature computation

Curvature being quite sensible to noise, level lines of the coalescence sets needs to be smoothed. Besides, the hierarchical structure inherited from the topographic map must be preserved. We thus rely on affine plane curve smoothing which is morphologically invariant [START_REF] Moisan | Affine plane curve evolution: a fully consistent scheme[END_REF]. This smoothing is applied directly to the polygonal models of level lines. A fixed sampling step at the curvilinear abscissa is used for every image: ∆s = 0.5 pixels. Let us stress that only one smoothing is operated here at a very fine scale. Our approach differs from a multi-scale approach such as [START_REF] Mokhtarian | Silhouette-based isolated object recognition through curvature scale space[END_REF].

For each curvilinear abscissa s i along a level line L, the curvature κ(s i ) is estimated as

κ(s i ) = θ ′ (s i ) ≈ θs i si+1 -θs i-1si (s i s i+1 + s i-1 s i )/2 , ( 2
)
where θs isi+1 is the local orientation of the segment s i s i+1 . This curvature estimation is less sensitive to noise than the ones based on first and second derivatives. It gives a curvature value relative to the pixel size.

Curvature is a signed information depending on the direction of shifting along a closed curve (clockwise or anticlockwise) and the local geometry (convex or concave). We choose the following convention. Negative level lines will be travelled clockwise. This induces negative curvature values at line convexities and positive values at line concavities. On the opposite, positive level lines will be travelled anticlockwise to invert the curvature sign. This convention first ensures us that pictorial elements such as stroke endpoints or junctions have the same local curvature sign wherever they are located in a coalescence set (Figure 2). Endpoints have always a negative curvature values, and regions where strokes create a non-reflex angle (e.g. at junctions) hold positive values.

anticlockwise). On this example (left), one negative level line delimits the outer contour, and one positive level line delimits the inner contour (right). Arrows notify the travel direction. This convention induces negative curvature values (red dots) and positive values (blue dots) wherever they represent an inner or an outer stroke cluster contour. This convention ensures us that pictorial elements such as stroke endpoints or junctions have the same local curvature sign.

Curvature zero-crossings contribute to artistic visual impression and object recognition [START_REF] Lowe | Perceptual Organization and Visual Recognition[END_REF]. They are estimated directly on the curvature signal. An example of a drawing with detected curvature zero-crossing is to be seen in

Detection of endpoints, junctions and corners

This section proposes to detect three different strong visual cues: stroke endpoints, stroke junctions and stroke corners, following to the taxonomy presented in Figure 4. In Section 3.3.1, we first select a set of candidates to that taxonomy.

Selection of candidates to a taxonomy of high curvature points

We first extract a large set of extrema that are candidates to the extrema point taxonomy shown in Figure 4. We iteratively consider every continuous portion of the curvature signal where |κ(i)| > κt with κt = 1/(kcR Φ ). We will discuss the choice of kc in Section 4.1. Zero-crossings of the derivative signal of each of these portion are considered as candidates to the taxonomy. Depending on the curvature sign, an extremum can be either positive or negative. An example of extrema detection is shown in Figure 3b. Thanks to the chosen convention on curvature presented in Section 3.2, stroke endpoints coincide with negative maxima (Figure 4b), and junctions coincide with several positive maxima (Figure 4d).

Corners have one positive maximum on the concave side, and one negative maximum on the convex side (Figure 4c). This information will be useful in the next three sections to analyze stroke endpoints, junctions and corners.

Stroke junctions

Stroke junctions are important geometrical characteristics of the 1D content of a drawing. Each junction indicate a possible occlusion denoting an actual level of perspective in the depicted scene [START_REF] Willats | Art and representation[END_REF]. A positive curvature maximum indicate possible junction nearby its location point pm in the drawing. Recall that R Φ is the tool radius. To state for a junction we center a disk of radius k j R Φ on pm. If there are three or more pieces of level lines covered by the disk, a stroke junction is detected and pm corresponds to one of the nonreflex angle of this junction. If there are only two pieces, this maximum is a line corner which will be characterized further in Section 3.3.3. Once all positive curvature maxima in a coalescence set have been analyzed, the ones that have been detected as belonging to a stroke junction are finally merged if their mutual distances are less than k j R Φ .

Stroke corners

A positive maximum that has not been classified as stroke junctions is considered as a stroke corner. Such maxima points have a very strong pictorial impact [START_REF] Leyton | The structure of paintings[END_REF]. Among the characteristics that have been previously proposed in the perception literature to measure the visual strength of a corner, we use the relative surface as a corner strength description [START_REF] Winter | Perceptual saliency of points along the contour of everyday objects: A largescale study[END_REF]. Considering a maximum pm, we iteratively consider the polygon p m-i . . . pm . . . p m+i for i > 0. This polygon is expanded while it does not contain any other point belonging to the coalescence set than the points {p m-j , . . . , p m+i }. The strength value is computed as the relative surface of this maximal polygon normalized by the total surface of the image.

If the underlying stroke describes a form that is perceptually visible such as an object contour, a line corner can be convex or concave and induces a different pictorial effect. This type effect has been recently studied in [START_REF] Fantoni | Contour curvature polarity and surface interpolation[END_REF][START_REF] Feldman | Information along contours and object boundaries[END_REF]. Recognizing the type of a corner, i.e. stating if it convex or concave, can be an ill-defined problem in a line-drawing. Indeed, strokes do not always describe a part of an object contour and objects described by an open contour also create some ambiguities. We propose a well-defined way to precise the visual orientation of a corner at stake, that we call sign of a corner. We say that a stroke corner is negative if it points toward the barycenter of the whole coalescence set to which it belongs, and positive if it points toward the opposite direction. The corner orientation is estimated with the orientation of the geometrical vector -----→ p m-iΦ pm + -----→ p m+iΦ pm, where i Φ is such that i Φ ∆s ≈ R Φ . This works in most of the situations. However a typical case where a corner sign and a corner type do not correspond is when the form describes a U-bend around the corner.

Stroke endpoints

To detect stroke endpoints, the set of negative curvature maxima is analyzed. Stroke endpoints are important components of the pictorial content. A drawing made of dotted lines delivers a different pictorial effect than the same drawing made with continuous strokes. For each negative maximum point pm, we consider the two points p i , p j where i < m < j along the level line at a keR Φ distance from the maximum. If the Euclidean distance between p i and p j is less than 3R Φ , and if there does not exist any point pq inside the triangle p i pm p j such that q < i and q > j, then pm is considered as an endpoint. This distance threshold tolerates a 50% precision error on R Φ . If the maximum is an endpoint, the distance between p i and p j should be indeed close to the width 2R Φ . This method is illustrated on Figure 5.

Indexing features

To illustrate how the proposed features can be used for the indexing of line-drawings, a similarity retrieval framework is considered in the experimental section. Based on the pictorial content analysis presented in previous sections, eleven scalar features are computed, summarized in Table 1. The three first features are computed on the distribution of curvature values κ(i) satisfying -κt ≤ κ(i) ≤ 0. Negative curvature values let us select one half of the curvature points along the strokes (see Figure 2). We thus indirectly consider points that are related to the underlying drawing stroke. We only consider values that respect κ(i) ≥ -κt to discard the strong curvature values that belong to stroke endpoints. These three features give an indication of how straight are the flat parts of the strokes independently of the rest of the geometrical information. Curvature values are normalized by the image diagonal. These three features thus become invariant to the drawing scale. This is required if we want that two images of a same drawing with two different resolutions have the same values.

Features 4 to 8 are linear densities of curvature zerocrossings, stroke endpoints, stroke junctions, positive and negative stroke corners. The linear density of each of these geometrical elements is defined as its total number in the image divided by the total length of the level lines of the coalescence sets in the image. Each linear density is normalized by the image diagonal. Features 9 and 10 are the sum of the strength values of convex and concave stroke corners. These two features are normalized by the image surface and therefore represent the relative image surface corresponding to each signed stroke corners. The last feature is the tool radius normalized by to the image diagonal.

Let us notice that linear densities of endpoints and junctions have been previously proposed by Julesz for perceptive preattentive discrimination of textures [START_REF] Julesz | Texton gradients: The texton theory revisited[END_REF]. We conducted several experiments on Julesz textons using the methods proposed in this paper. They led to the same perceptive results as in Julesz studies.

Table 1: Features vector used for similarity retrieval

Experiments and results

Section 4.1 studies the setting of the three parameters that have been used to extract the characteristic points: kc, ke and k j . Some similarity retrieval results on a database of artistic line-drawings are discussed in Section 4.3.

Setting of parameter kc, ke and k j

Ideally, κt is set in such a way that the three first features of Table 1 do not depend on the presence of endpoints. This is not true when kc is too low since endpoints regions may be integrated in the distribution. On the opposite, a too high value of kc raises the risk of not discriminating between a drawing holding strictly straight linear parts and another made of slightly curvy strokes. A good tradeoff has been empirically set with kc = 5, using synthetic images made with different tool radius and real line-drawings.

Endpoints are detected using a triangle of scale keR Φ around a maximum candidate (Section 3.3.4). Increasing ke raises the risk to miss some endpoints that are very close to a junction. Lowering ke raises the risk to detect some elongated junctions of two strokes as an endpoints. Stroke junctions are detected using a disk region whose radius is k j R Φ . Increasing k j raises the risk to consider too large regions than what is necessary to cover a stroke junction, and thus provoke some false detection. Lowering k j raises the risk to miss some stroke junction. Empirical experiments with real and synthetic drawings led us to choose k j = ke = 5. It is not surprising that k j and ke hold the same value since they are related to the same local scale of extremum neighborhood. In the remaining part of this paper, we will thus use those two values.

Detection of endpoints, junctions and corners

Three examples of analyzed drawings are to be seen in

Similarity results on a line-drawings database

In this section, we rely on a database made of 105 linedrawings. Most of this database comes from an exhibition catalogue of drawings made by Henri Matisse and Ellsworth Kelly [START_REF] Remi | Henri Matisse and Ellsworth Kelly -Dessins de plantes[END_REF]. Images have been scanned with an HP Scanjet 8200 with 4800 ppi resolution. 21 images come from several Picasso sketch books released on DVD-Rom in 2006 as a 800 × 600 digital facsimile edition [START_REF] Rmn) R | Les carnets de picasso[END_REF]. Finally, 31 computer drawings from [START_REF] Explosion | Hallogram[END_REF] have also been included.

Computational times for the complete features extraction procedure takes around 15 seconds for a 800 × 600 image using a PC Pentium IV running at 4.3gHz. Most of this time is dedicated to the stroke contours detection. Thanks to the Euclidean distance and the size of the feature vector, a query runs almost instantaneously.

Each drawing is described by its eleven scalar features (Table 1). Once each feature is computed for every drawing in a database, it is centered by the mean and normalized by the standard deviation of the feature distribution over the database. Let V i (k), j ∈ [1, 11] be the eleven normalized features of a drawing, then the similarity measure d(V 1 ,V 2 ) between two drawings is the following Euclidean distance:

d(V 1 ,V 2 ) = (∑ 11 k=1 (V i 1 (k) -V i 2 (k))
2 ) 1/2 . Retrieval results on this database are shown in Figure 7. On these figures, the query image is on the top left corner, followed by the closest results according to the Euclidean distance. First example is based on the query image of Figure 6a. All results hold a high linear density of endpoints. This example gives a good insight of what we call picto-rial effect, results being semantically different, yet delivering a very similar pictorial expression. Same remark can be done for examples 2, 4, 7 and 8. Examples 3, 6, 5 and 9 are very similar both from the pictorial and semantical point of view. The visual features of the represented subject are indeed very similar and discriminative compared to the rest of the database content. The two last results of the 10 th example for instance are yet so different from the query and the three closest images that they look unsatisfying, knowing that more similar rendered birds are in the database. This points out a limit of the proposed description: in some cases, the semantic impact delivered by the represented subject might overwhelm the pictorial effect. In the image retrieval context, this limit could be moved away by integrating some subject metadata when available, if someone seeks images having the same semantical content and being rendered with the same style as the query. The whole database and all retrieval results are to be seen online on a website † . This database has also been classified by a professional artist into 13 classes. Discounted cumulative gain (DCG) and Bull-eye percentage (BEP) [START_REF] Shilane | The Princeton Shape Benchmark[END_REF] are 93.8% and 88.2% respectively using these classes.

Conclusion and future work

In this paper we introduced a method for the pictorial content analysis of line drawings using the geometrical information of stroke contours. Computer Graphics and Non-Photorealistic Rendering (NPR) could benefit from methods that are able to analyze the artistic content. This would enable style transferring for the rendering of drawings or line textures. We believe that the proposed methods are a promising step toward this objective. This work opens several perspectives. First we are presently collaborating with several artists to evaluate the efficiency of the proposed features. These artists will be asked to classify a database on which we will rely to quantify the performance of our approach. Another perspective concerns more sophisticated similarity measures using some kernel optimization for non-uniform weighting. We are also investigating the possibility to describe separately every coalescence set in a line-drawing and to model the comparison of two groups of coalescence sets through a global optimization procedure, e.g. by relying on the Earth Mover Distance [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF]. This paper is limited to line-drawings, that is drawings fully based on 1D content.

Our goal was to observe the importance of 1D content independently of other types of primitive. It would be interesting to further investigate some methods to extract 1D content in any type of artwork as for example in paintings.

Figure 2 :

 2 Figure 2: During curvature computation, convention is chosen so that negative (resp. positive) level lines are travelled clockwise (resp.anticlockwise). On this example (left), one negative level line delimits the outer contour, and one positive level line delimits the inner contour (right). Arrows notify the travel direction. This convention induces negative curvature values (red dots) and positive values (blue dots)wherever they represent an inner or an outer stroke cluster contour. This convention ensures us that pictorial elements such as stroke endpoints or junctions have the same local curvature sign.

Figure 3 :

 3 Figure 3: Left, maximal meaningful level lines from a drawing, and its superposed curvature zero-crossings. Right, curvature maxima detected with the method presented in Section 3.3.1. Red (resp. green) diamonds are positive (resp. negative) maxima.

Figure 4 :

 4 Figure 4: Taxonomy of high curvature values along stroke contours: a) a simple high curvature point along the gesture trajectory, b) an endpoint, c) a stroke corner, d) two types of stroke junctions: X and T junctions. Plus and minus sign refers to the curvature sign on the stroke contour. Strokes are grey lines, and stroke contours are blue surrounding lines.

c

  The Eurographics Association 2008.

Figure 5 :

 5 Figure 5: Left: negative maximum is classified as a stroke endpoint. Right: some points pq which are not belonging to the level line portion [p i , p j ], are included in the triangle. This maximum is thus not classified as a stroke endpoint.

Features 1 .

 1 mean of the curvature distribution of κ(i) ∈ [-κt , 0] 2. standard deviation of the curvature distribution of κ(i) ∈ [-κt , 0] 3. kurtosis of the curvature distribution of κ(i) ∈ [-κt , 0] 4. linear density of curvature zero-crossings 5. linear density of stroke endpoints 6. linear density of stroke junctions 7. linear density of positive stroke corners 8. linear density of negative stroke corners 9. sum of the depths of positive stroke corners 10. sum of the depths of negative stroke corners 11. tool radius normalized by the image diagonal

Figure 6 :

 6 Figure 6: Examples of analyzed line-drawings. Green dots are detected stroke endpoints, black dots are detected stroke junctions, and green (resp. red) stars are positive (resp. negative) stroke corners. Coalescence sets have different colors. Left, middle and right images are the queries for the first, second and third retrieval example respectively in Figure 7.

Figure 7 :

 7 Figure 7: Similarity retrieval results in 105 line-drawings database. For each example, image query is in the upper-left corner followed by its five closest images.

c The Eurographics Association 2008.

† URL: http://www.hurtut.net/geometry/ c The Eurographics Association 2008.


Thomas Hurtut

Yann Gousseau

Francis Schmitt

Adaptive image retrieval based on the spatial organization of colors

Keywords: Spatial organization of colors, color image retrieval, a contrario method, earth mover's distance, optimal transport, dead leaves model, dead leaves model, image distance

come

Introduction

Image retrieval is the process of finding an image or a group of images of interest in an image database. It usually relies on some distance or dissimilarity measure between image features that takes information about color, shape or texture into account. Since the seminal work of Swain and Ballard [START_REF] Michael | Color indexing[END_REF], many image retrieval systems have been strongly oriented toward representing the color content of images by using histograms. This is probably due to the ability of color histograms to concisely represent the global aspect of natural images and to the striking impression of global coherence given by series of images having similar color content. Tutorials about such methods can for instance be found in the monographs [START_REF] Bimbo | Visual information retrieval[END_REF][START_REF] Lew | Principles of Visual Information Retrieval[END_REF].

More recently, many works took interest in the spatial organization of colors and in dependence relationships between pixels. A first type of approach to this problem has been to select some spatial statistics between pixels color, such as correlograms or some filter responses, see [START_REF] Pass | Comparing images using joint histograms[END_REF][START_REF] Huang | Spatial color indexing and applications[END_REF][START_REF] Ciocca | Image indexing and retrieval using spatial chromatic histograms and signatures[END_REF][START_REF] Siggelkow | Improvement of histogram-based image retrieval and classification[END_REF][START_REF] Lambert | Image retrieval using spatial chromatic histograms[END_REF]. A related category of methods makes use of points of interest within images in a way similar to classical object recognition methods in computer vision, see [START_REF] Heidemann | Combining spatial and colour information for content based image retrieval[END_REF][START_REF] Van De Weijer | Coloring local feature extraction[END_REF]. In both cases, the features of color images that are accounted for are various types of color transitions, usually using various scales and orientations.

A second type of approach aims at describing and comparing images with respect to their global spatial content, that is the coarse organization of colored regions in the image. Region-based approaches first perform a segmentation of images which is then used for image comparison. The Blobworld image representation, for instance, uses multi-dimensional Gaussian mixtures to represent the image [START_REF] Carson | Blobworld: image segmentation using EM and its application to image querying[END_REF]. Many other unsupervised segmentation methods have been used to index images, see e.g. [START_REF] Hsu | An integrated color-spatial approach to color image retrieval[END_REF][START_REF] Smith | Tools and techniques for color image retrieval[END_REF][START_REF] Smith | Image classification and querying using composite region templates[END_REF][START_REF] Wang | Simplicity: semantic-sensitive integrated matching for picture libraries[END_REF][START_REF] Rugna | Color coarse segmentation and regions selection for similar images retrieval[END_REF][START_REF] Dvir | Context-based image modelling[END_REF][START_REF] Prasad | Regionbased image retrieval using intergrated color, shape, and location index[END_REF][START_REF] Jing | An effective region-based image retrieval framework[END_REF][START_REF] Liu | Region-based image retrieval with high-level semantic color names[END_REF]. Other methods rely on a small number of predefined image regions, see [START_REF] Stricker | Color indexing with weak spatial constraints[END_REF][START_REF] Lipson | Configuration based scene classification and image indexing[END_REF]. Several works specifically study the mutual arrangement of regions in images. Representation and comparison of these relationships can be symbolic [START_REF] Sk Chang | Iconic indexing by 2-D strings[END_REF] or nonsymbolic [START_REF] Berretti | Weighted walkthroughs between extended entities for retrieval by spatial arrangement[END_REF]. The goal of a global description is to be able to detect images that have a similar color layout, such as "a green stripe on the top of the image, above a red blob that is itself in between two blue blobs". The present work is concerned with the development of a method for image querying relying on the similarity of such a global, coarse organization of colors. In what follows, this aspect of image content will be denoted by spatial organization of colors.

The spatial organization of colors is an important aspect of the structure of an image and definitively one of the first to be perceived. It is an intermediate feature between low-level content such as color histograms and image semantics. In the context of image retrieval, various taxonomies of both image content and user-needs have recognized its importance. It is part of the visual relationships stressed by Burford et al. [START_REF] Burford | A Taxonomy of the Image: On the Classification of Content for Image Retrieval[END_REF], which is characteristic of the spatial arrangement of the scene in two dimensions and of the color content. In the cultural heritage domain, it also concerns the cultural abstraction content since it models the artistic composition [START_REF] Eakins | Image retrieval interfaces: A user perspective[END_REF][START_REF] Chen | An analysis of image retrieval tasks in the field of art history[END_REF]. Such a feature can be used for illustration, generation of ideas, aesthetic decoration, and learning [START_REF] Conniss | Information Seeking Behaviour in Image Retrieval VISOR I Final Report[END_REF]. Within one semantic category reached with metadada for instance, a journalist may pursue its illustration needs according to composition clues. Image retrieval with such features may be a source of ideas for finesarts students, professors, artists and curators. Image retrieval using the spatial organization of colors can also be serendipitous, meaning one can accidentally discover a fortunate association of semantically different images delivering the same visual impact, which is practically impossible using metadata or ontology-based description.

Our aim in this paper is twofold. First we introduce a robust method enabling the comparison of the spatial organization of color between images. This method does not require any segmentation step. By making use of the well-known Earth Mover's Distance (EMD), first introduced in [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF] in the context of image querying, we show that it is unnecessary and sometime misleading to segment the image into regions. The comparison of the spatial organization of colors between two images is addressed as a global optimization procedure. To achieve reasonable computational costs, the EMD is computed between pixels of coarse thumbnail representations of images. Second we address the matching step of our querying system. In short, given a query image, we want to decide automatically which images from a database should be paired with the query. In particular, this yields a retrieval system in which the number of result images may vary from one query to another. This is what we refer to as adaptive image retrieval. This question has rarely been addressed by the image querying literature. It is however of primary concern since databases tend to become huge and knowing the number of images that are similar to a query gives a quick insight on the content of a database. To reach a decision, we follow an a contrario approach similar to the one in [START_REF] Musé | An a contrario decision method for shape element recognition[END_REF]. Roughly speaking, the system returns images for which the proximity to a query image is very unlikely to be due to chance. More precisely the method relies on a model of random generic color layouts to estimate the probability that two images are at a small distance apart by chance. To the best of our knowledge, this is the first proposition of an automatic matching threshold in the context of color image retrieval.

This article extends our conference paper [START_REF] Hurtut | Spatial color image retrieval without segmentation using thumbnails and the earth mover's distance[END_REF], providing more detailed description and discussion. The plan of the paper is as follows. Section 2 introduces the distance used to compare the spatial organization of colors. In Section 3, an automatic matching threshold is derived in the context of image querying by example. Finally, the efficiency of the method is illustrated in Section 4 by querying two databases : a database of 1500 high quality illuminated manuscripts and a database of 15000 non-specialized digital photographs. This section includes an analysis of the behavior of the matching threshold depending on the database and the query. We give perspectives and conclude in Section 5.

Global comparison of the spatial organization of colors

In this paper, the problem of comparing the spatial organization of colors between two images is viewed as a global optimization procedure, in which the best transform matching an image to the other one is sought. The best transform refers to a set of pixels displacements that optimize some cost function relying on the spatial and color differences between pixels.

Earth Mover's Distance

Intuitively, if one is allowed to move the pixels of a first image and to change their color, one seeks the minimal amount of work needed to match a second image. The natural mathematical modeling of this approach relies on the theory of Optimal Transport, initiated by G. Monge at the end of the 18th century. The application of this theory to the problem of image matching and image querying is not new, and started with the work of Rubner et al. [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF].

In [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF], the Earth Mover's Distance (EMD) is introduced as the smallest amount of work needed to match a set of weighted features (f i , w i ) i=1,...,n to another one. The f i 's are features belonging to some space E on which a metric (i.e. a distance) d e is defined, and the w i ≥ 0 are the weights. To illustrate this notion, the authors suggest an analogy between EMD and the minimal amount of work needed to put some mass of earth spread in space (the first set of features) into a collection of holes (the second one), which is precisely the problem addressed by G. Monge more than two centuries ago. In this context, the work corresponds to the quantity of earth that is displaced times the length of the displacement. Formally, the EMD between two sets of weighted features

f 1 = (f 1 i , w 1 i ) i=1,...,n1 and f 2 = (f 2 i , w 2 i ) i=1,.
..,n2 having same total weight1 i w 1 i = j w 2 j is defined in [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF] as

d(f 1 , f 2 ) = min fij n1 i=1 n2 j=1 d e (f 1 i , f 2 j )x i,j , (1)
where the flow x i,j is constrained by

x i,j ≥ 0, ∀i = 1, . . . , n 1 , n2 j=1 x i,j = w 1 i , ∀j = 1, . . . , n 2 , n1 i=1 x i,j = w 2 j .
With the same analogy as before, the first equation means that all available earth is displaced, while the second one means that holes are completely filled.

In the context of color image retrieval, features f i either correspond to some color in a 3-dimensional space equipped with a perceptually adapted metric or to some point in a 5-dimensional space accounting for both color and spatial information. In the first case, the EMD compares the global color content of images. It is shown in [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF] that this method is more satisfying than approaches relying on bin to bin distances on histograms that are very sensitive to the chosen quantization and to small color changes. In the second case (features in a 5-dimensional space), the (f i , w i ) are clusters representing important color regions of the image. In [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF], the clustering relies on the use of recursive splitting and k-d trees. This same approach can rely on any region-based segmentation method by assigning a position and various features to each region. This has been recently investigated by various authors [START_REF] Dvir | Context-Based Image Modelling[END_REF][START_REF] Jing | An effective region-based image retrieval framework[END_REF][START_REF] Liu | Region-based image retrieval with high-level semantic color names[END_REF]. The region integrated method developed in [START_REF] Wang | Simplicity: semantic-sensitive integrated matching for picture libraries[END_REF] can also be considered as a simplification of the same idea.

Applying EMD to pixels

Instead of relying on a preliminary segmentation to perform region-based image comparison, we suggest using EMD directly at the smallest possible detail level. We thereafter explain why this approach is more robust than using segmentations and allows the same degree of efficiency. An image is associated with a set of features (f i , 1/n) i=1,...,n where n is the number of pixels in the image and f i the color and position of the ith pixel. Each feature is given the same weight because pixels correspond to a fixed area in the original image. In this case (see [START_REF] Ford | Solving the transportation problem[END_REF][START_REF] Levina | The earth movers distance is the Mallows distance: Some insights from statistics[END_REF]), the computation of the EMD is reduced to an assignment problem and is therefore greatly simplified:

if

f 1 = (f 1 i , 1/n) i=1,...,n and f 2 = (f 2 i , 1/n) i=1,.
..,n are two sets of identically weighted features in a metric space (E, d e ), then

d(f 1 , f 2 ) = min φ n i=1 d e f 1 φ(i) , f 2 i , (2)
where φ belongs to the set of permutations of {1, . . . , n}. This means that when features all have the same weight, they are not split apart by EMD: each feature of the first image is assigned to exactly one feature of the second one. Each x i,j in Equation ( 1) is either one or zero. Computation of the EMD is therefore simplified and can be achieved by solving an optimal assignment problem [START_REF] Burkard | Linear assignment problems and extensions[END_REF][START_REF] Burkard | Assignment Problems[END_REF]. In the absence of specific structure for the cost matrices involved in the computation of Equation (2) (in particular they are not sparse), we chose to solve the assignment problem using the classical Hungarian method, also called the Kuhn-Munkres algorithm [START_REF] Kuhn | The hungarian method for the assignment problem[END_REF], which has a complexity upper bound in O(n 3 ). An alternative to this choice would be to bound the displacements of pixels, yielding sparse cost matrices and enabling faster computations, see [START_REF] Burkard | Linear assignment problems and extensions[END_REF][START_REF] Burkard | Assignment Problems[END_REF].

Using thumbnails

Of course, even with this simplification, the spatial resolution of even a small image entails excessive computational time when querying large databases. Therefore, the method is applied to coarsely sampled thumbnails, see Figures 1 and2. This is not a strong limitation, as we are interested in the global color content of images, and not in their small scale structure. This point will be further discussed in the next paragraph. The total number of pixels of these thumbnails is fixed to n, depending on the database content and on computing constraints. The influence of this important parameter of the method will be studied in Section 4.1. The typical order of magnitude of n is several hundreds. Image sub-sampling is achieved by averaging the image pixels. This is done first by gamma correcting the numerical RGB image components. The resulting tristimulus values can then be summed according to the Grassman's law of the colorimetry [START_REF] Wyszecki | Color Science: Concepts and Methods, Quantitative Data and Formulae[END_REF]. Assuming the image colorimetric space to be sRGB space, thumbnails are then converted into the psychometric CIELab color space [START_REF] Anderson | Proposal for a Standard Default Color Space for the Internet-sRGB[END_REF] in order to enable the use of the Euclidean distance to compare color between their pixels. The signature for an image is therefore composed of the n thumbnail pixel features

f i = {L i , a i , b i , X i , Y i },
where L, a, b are the color coordinates and X, Y the spatial positions.

Thumbnails versus segmentation

When using EMD with region based methods, each region is usually weighted according to its area. In this case, EMD breaks apart the weights of the features when optimizing the cost of the mapping. Therefore, the usefulness of a segmentation step before EMD is questionable since pixels that have been grouped into regions can be mapped separately to different target regions by EMD. Indeed no constraints are included in EMD that take into account the fact that regions represent groups of pixels. On the other hand, a segmentation method may fail and produce an under-segmentation which can yield regions that are not well represented by for instance their mean color, thus affecting the performances of the matching step. In short, using large regions with EMD is useless in favorable cases and misleading in unfavorable cases. Our approach is simpler, avoids the complex tunings of segmentation methods, is more robust and takes full advantage of EMD: we oversegment the image into small and regularly spaced regions (the thumbnail pixels), and then let the EMD distance resolve the matching problem. These regions are chosen to be of fixed size (using thumbnails) to reduce the computation of EMD to an assignment problem. Of course, using a fixed number of elementary regions (the thumbnail pixels) can be seen as a disadvantage. However, as stressed in the previous discussion, it is only necessary to use a number of pixels that is large enough to account for coarse regions in the image. The number of regions that are considered is only limited by computation constraints. We claim that it is generally enough to use an order of magnitude of several hundred pixels to achieve this. This fact is corroborated by the number of regions usually considered by Region Based Image Retrieval systems. The influence of the size of thumbnails will be further studied in the experimental Section 4.1.

Distance

The last tool that is needed in order to compare two images is a distance between two pixel features f 1 i and f 2 j . We choose to use a weighted exponential distance:

d e (f 1 i , f 2 j ) = α *   1 -e - ((L 1 i -L 2 j ) 2 +(a 1 i -a 2 j ) 2 +(b 1 i -b 2 j ) 2 ) 1/2 δc   +(1 -α) *   1 -e - ((X 1 i -X 2 j ) 2 +(Y 1 i -Y 2 j ) 2 ) 1/2 δs   . (3)
The parameter α balances color and spatial contributions. Parameters δ c and δ s are chosen according to the color and spatial dynamics. Typical values used in this paper are α = 0.5, δ c = 15 and δ s equal to a quarter of the thumbnails diagonal. The weighted exponential distance supports the idea that beyond a certain distance, two different images are just considered as not similar. This nonlinear distance is adapted to an image retrieval task, where a target image is close to a query, in which case the similarity is measured, or is simply different and penalized. All of them have the same weight. In this context, EMD is equivalent to an assignment problem. Each pixel of the first thumbnail (top middle) is assigned to exactly one pixel of the second one (bottom middle). To illustrate the EMD operation, the thumbnail on the bottom right represents the optimal permutation of the pixels from the first thumbnail to match the pixels of the second one. Pixels of the red region are moved towards their location in the second image. The white column of paper is also correctly moved from right to left.

3. Automatic matching threshold using an a contrario approach.

In this section, we consider the task of adaptive image retrieval. Our framework is the classical query by example. Even though this framework is known to have some limitation (see for instance [START_REF] Fang | An interactive system for mental face retrieval[END_REF]), it is nevertheless usually one of the building blocks of a querying system and it enables us to address the problem we are concerned with : how to decide which answers are meaningful when querying a very large amount of data. More precisely, we tackle the following problem: which images from a database should be matched with a query image?

A contrario method

Our approach is similar to the one in [START_REF] Musé | An a contrario decision method for shape element recognition[END_REF][START_REF] Gousseau | Comparaison de la composition de deux images, et application à la recherche automatique[END_REF] : we match two images as soon as their proximity is unlikely to be due to chance. This is an application of the general principle of a contrario methods: we do perceive events that are unlikely in a white noise situation. Such methods have been successfully applied to several tasks in computer vision following the work of [START_REF] Desolneux | Meaningful alignments[END_REF]. A recent exposition can be found in the monograph [START_REF] Desolneux | From Gestalt Theory to Image Analysis: A Probabilistic Approach[END_REF]. In the field of shape matching, related ideas can be found in [START_REF] Lindenbaum | An integrated model for evaluating the amount of data required for reliable recognition[END_REF][START_REF] Grimson | On the verification of hypothesized matches in modelbasedrecognition[END_REF]. In short, the method works as follows. We first introduce a probabilistic model of generic images, called a background model. Given a query image, we then match it to images that are at a distance that is too small to be plausible under the background model. Such matches will be called meaningful matches. We now detail how this principle can be applied to the retrieval of color images.

Setting

We write Q for a color query image, and B = {T 1 , . . . , T mB } for a database composed of m B target color images. For two images Q and T , their EMD distance (as detailed in the previous section) is denoted as d(Q, T ). The goal here is to automatically compute a threshold for these distances that is able to distinguish between matches and non-matches. Using a classical statistical setting, we want, for each image T , to develop a test for the hypothesis H 1 :

H 1 : " T is similar to Q ",
relying on the quantity d(Q, T ).

Meaningful matches

We make use of an a contrario approach similar to the one of [START_REF] Musé | An a contrario decision method for shape element recognition[END_REF]. The test relies on a background model2 M of generic images. The threshold on distances is fixed in such a way that it is very unlikely that an image generated by M would match the query. In this sense, matches are rare events with respect to M. More precisely, we fix the rate of false alarms when testing H 1 against H 0 : H 0 : " T follows the background model ".

The main idea is to fix the matching thresholds by bounding the probability that two images are matched by chance. Assume that we are able to compute the probability P r(d(Q, T ) ≤ c), for any c ≥ 0, where T is an image distributed according to M (we write P r H0 for this probability). Recall that m B is the number of images in the database B. For ǫ > 0, we say that

T is an ǫ-meaningful match of Q if d(Q, T ) ≤ c Q , where c Q is defined as c Q = sup c : P r H0 (d(Q, T ) ≤ c) ≤ ǫ m B . (4)

Control of the number of false matches

The upper bound ǫ/m B in the previous definition is chosen because it ensures that: if all images T i from the database B follow the model M, then the expectation of the number of ǫmeaningful matches is less than ǫ.

The proof of this fact is quite simple and relies on the linearity of the expectation. Indeed, if T 1 , . . . , T mB are distributed according to M, then, by writing E for the mathematical expectation, we obtain,

E(there exists i such that d(Q, T i ) ≤ c Q ) ≤ 1≤i≤mB P r H0 (d(Q, T i ) ≤ c Q ) ≤ m B ǫ/m B = ǫ.
In short, this means that ǫ-matches are defined so that, on the average, there are fewer than ǫ casual matches when querying a generic database having the same size as B. This is possible because the matching threshold depends upon the size of the database. In practice, the value ǫ = 1 is always considered in this paper, that is, the expected number of false matches is fixed at one. In what follows, the term meaningful match is used instead of 1-meaningful match.

Background model

Of course, the background model M now has to be defined. The main idea when applying a contrario methods is to use a background model in which the basic elements (in our case the spatial color content) are independently distributed. Drawing our inspiration from [START_REF] Desolneux | Edge detection by Helmohltz principle[END_REF][START_REF] Musé | An a contrario decision method for shape element recognition[END_REF], we choose a model leading to images that have similar color marginal distributions as images from the database B. Since the distance d only takes thumbnails into account, we seek a background model for thumbnails.

The first idea was to use a uniform white noise model, in which all pixels of a thumbnail are drawn independently and identically in the CIELab space following the color distribution observed on B. This was not satisfactory, since too many images were then meaningful matches of the query. Our interpretation of this fact is that, even when they are perceptually different, images share some common structures (e.g. the presence of homogeneous regions). In short, the high number of matchings follows from the fact that two random natural images are often closer one to another than they are to realizations of white noise.

A more structured model is needed, with geometric features similar to the ones encountered in real images, such as homogeneous regions and edges. The model we consider is the dead leaves model ( [START_REF] Matheron | Modèle séquentiel de partition aléatoire[END_REF][START_REF] Bordenave | The dead leaves model: a general tessellation modeling occlusion[END_REF], consisting in the superimposition of random objects, which provides a simple way to generate random images with uniformly colored regions. The random objects are assumed to be independently and identically distributed (i.i.d.). This model has been experimentally proved to be well adapted to the structure of natural images when using power law distributions for the size of objects [START_REF] Alvarez | The size of objects in natural and artificial images[END_REF][START_REF] Lee | Occlusion models for natural images: A statistical study of a scale invariant dead leaves model[END_REF]. For the sake of simplicity, we restrict ourselves to an isotropic model where objects are simply disks. The radius of these disks are i.i.d. random variables with density f (r) ∝ r -γ . The model is then characterized by three parameters : the scaling parameter γ, and r 0 , r 1 the minimal and maximal sizes of objects 3 . We chose γ = 3, a typical value for natural images, r 0 = 1 and r 1 of the same magnitude as the dimensions of thumbnails. Each object is colored independently of its size and independently of other objects. Observe that this implies that the marginal color distribution of the model is the same as for objects. The distribution of the color of objects (and therefore also the color marginals of the model) is then learned from the color distribution of pixels of the thumbnails from the database. Samples of the model (using color marginals from the database of illuminated manuscripts used in the experimental section) are displayed in Figure 3. The last step, in order to be able to compute meaningful matches of a given query image, is to compute the matching threshold c Q such that P r(d(Q, T ) ≤ c Q ) ≤ ǫ/m B where T follows M. 3 One can get rid of r 0 by considering a limit model yielding the same small scale structure as in natural images [START_REF] Gousseau | Modeling occlusion and scaling in natural images[END_REF], but this model is unnecessarily complex for the modeling of coarse representations of images.

This quantity relies on geometrical properties of the dead leaves model for which no explicit formula are known. We therefore rely on numerical simulations. The probability is approximated by a frequency over a large number of samples. The choice of the number of samples will be addressed in Section 4.4.

Discussion

A more classical approach to the statistical testing of hypothesis H 1 (an image from the database is similar to the query) would be to use a probabilistic model for the images that are similar to the query and, for example, to perform a Bayesian test. This would imply the building up of specific models for each category of color images that we deal with, a task that is not compatible with a generic contentbased querying system. An alternative is provided by a contrario methods, as explained above.

A major advantage of the proposed approach to the matching problem is that the threshold automatically adapts to both the query and the database. First, the threshold is expected to become more conservative as the database gets larger, which can be crucial when dealing with very large databases. Second, since the threshold is computed using color marginal learned on the database B, it adapts to the scarceness or the commonness of the query image. This can be especially useful when dealing with specialized bases. More generally, we will observe in the experimental section that for a fixed ǫ the matching threshold value c Q varies very much from one query image to another. In fact, this is one of the main advantages of the proposed method : the threshold is only driven by the parameter ǫ, which intuitively corresponds to a number of false matches when submitting a query and is therefore easier to choose than an upper bound on distances. In all experiments to be performed in this paper, we choose ǫ = 1.

Experiments

Tests have been performed on two databases. The first one is an illuminated manuscript database provided by the Institut de Recherche et d'Histoire des Textes (IRHT), thereafter referred to as the IRHT database. This database contains 1500 high quality and color-calibrated illuminations. The second database, the CLIC [START_REF] Moellic | Evaluating content based image retrieval techniques with the one million images CLIC testbed[END_REF] database, consists in 15000 non-specialized color photographs and has been gathered by researchers from the French Com-missariat à l'Energie Atomique (CEA). The first database is especially adapted to the purpose of this paper since it contains a large variety of colorful arrangement of patterns. As stressed in the introduction, we believe our approach to be particularly adapted to the querying of cultural heritage databases since the spatial organization of colors summarizes the artistic composition. The CLIC database is less varied in terms of spatial organization despite its larger size. Yet it is quite useful for studying the behavior of the proposed matching threshold. Indeed, it contains very large series of images with similar compositions, such as landscapes in which the sky more or less covers the top half of the picture.

Besides section 4.1 where the choice of thumbnail size n is discussed, all experiments to be presented here use n = 150 pixels (10 × 15 or 15 × 10). Except where mentioned, tests are performed using parameters α = 0.5, δ c = 15 and δ s = 5 for the distance d e (see Equation ( 3)). The parameters used for the background model are γ = 3, r 0 = 1 and r 1 = 30, as justified in Section 3. The expected number of false detection is fixed at one, that is ǫ = 1.

A point that is rarely discussed and that concerns all spatial color methods is how to manage different image ratios. Various choices are available to handle this problem. In the experiments, we chose to separate portrait and landscape orientations. However EMD has no constraints on signature length and total weight (see [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF]), and therefore we could also have compared images with different ratios. Yet, this would transforms the assignment problem into a transportation problem leading to a bigger computational load since transportation problems are usually solved using the simplex method.

In Section 4.1, we illustrate the compromise between thumbnail size, result accuracy and computational time. In Section 4.2, we illustrate the effect of the parameter α (see Equation ( 3)) that balances color and spatial content, and in Section4.3, we compare our results with a region-based method relying on a segmentation followed by the use of EMD, as discussed in Section 2. Both experiments are performed on queries from the IRHT database. Next, we explain in Sectionsec:variability how to numerically compute the matching threshold and we study its variability according to both the query and the database. In Section4.6, we present several query examples from the IRHT and CLIC databases. Before proceeding, let us mention that our study does not include any quantitative evaluation of the results, for instance through the use of a pre-classified database. Yet the definition of the matching threshold would have enabled the comparison between predicted and obtained false matchings. Moreover, it would have been interesting to evaluate the efficiency of the threshold using ROC curves, for instance. However, to the best of our knowledge, none of the existing predefined classification of image databases relies on the spatial organization of colors. Therefore, the efficiency of the retrieval method and the precision of the matching threshold are essentially illustrated by examples, comparisons and a study of the variability of the matching threshold.

Size of thumbnails

In this section we illustrate the effect of the size of thumbnails on some query examples. As explained in Section 2, solving the assignment problem (2) for images having a standard resolution is not computationally feasible. Two query examples have been computed with the following thumbnails size n = 150, 600, 2400 (corresponding to 10×15, 20×30 and 40 × 60 for horizontal thumbnails). Results are displayed in Figure 4 and5. One can see on these two examples that increasing the thumbnail size has a low impact on results. We observed similar results on other queries from the two databases used in this experimental section. This illustrates that an order of magnitude of several hundreds pixels is enough to achieve coarse description of the spatial organization of colors. Using a PC under Linux with Pentium IV running at 4.3gHz, and C programming code, computational time for querying a database of 1000 images is 18 seconds for n = 150, compared to about 10 minutes and 8 hours for n = 600 and 2400 respectively. This is coherent with the O(n 3 ) complexity upper bound of the Hungarian method used to solve the assignment problem. In the rest of the paper, thumbnails of n = 150 pixels are always used.

Color-spatial weigthing

As detailed in Section 2, we use EMD to measure the distance between signatures of images. The weighting parameter α ∈ [0, 1] biases the distance d e between two features (see Equation ( 3)) towards either the spatial or the color constraint. If the user wishes to be severe on spatial organization (pixels cannot move freely on the image), α must be set at a low value. The higher the value, the more the objects can move and be split apart. The extreme case where α = 1 corresponds to a color histogram method, with EMD as a similarity measure. Two results are displayed in Figure 6, using different values for α.

Comparison with the use of a segmentation step before computing EMD

In this section, we address the usefulness of a segmentation step before applying EMD, as proposed for instance in [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF][START_REF] Dvir | Context-Based Image Modelling[END_REF][START_REF] Jing | An effective region-based image retrieval framework[END_REF][START_REF] Liu | Region-based image retrieval with high-level semantic color names[END_REF]. In our experiments, we choose the Blobworld segmentation 4 [10] that represents regions by using Gaussian mixtures. The original method uses texture and location to extract blobs. We only use color and location information for the results to be directly comparable with our method. We nevertheless observe similar retrieval a) b) Fig. 6. Color spatial weighting using parameter α. Two query images (framed in red on the left) are followed by their three nearest matches with α = 0.5 in Figure a) and α = 1 in Figure b). This second case corresponds to color histogram retrieval methods using EMD. We see in Figure b) the loss in the spatial organization.

results on the IRHT database when also using texture. The number of blobs per image is an important parameter. It is set by an MDL method (see [START_REF] Carson | Blobworld: image segmentation using EM and its application to image querying[END_REF] for details). For the sake of coherence in our comparison, we also use EMD to measure the distance between two blob signatures, the weight of each blob being equal to its size.

Two comparative results are displayed in Figures 7 and 8. In all our tests, the blobs+EMD method performs well on queries where simple objects stand on a homogeneous background. In this case both methods perform similarly. As explained before, since the weights of features (regions) are not equal, EMD breaks them apart when optimizing the cost of the mapping. Such a situation is displayed in Figure 7. On the other hand, our method performs better where segmentation methods fail. Indeed, wrong segmentations yield regions that are not well represented by a single feature, highly affecting the performances of the matching step.

Numerical computation of the matching threshold

Recall that in Section 3.3, we introduced a method to compute a matching threshold when submitting a query image Q to a database B. This threshold, c Q , depends on both a positive parameter ǫ and the query Q. It enables to define the meaningful matches of Q as images T for which d(Q, T ) ≤ c Q . Recall also that c Q is defined by Equation (4). As explained at the end of Section 3.5, the practical computation of c Q relies on numerical simulations. Namely, it is estimated as an empirical frequency. We first sample n bg ×m B realizations M j of the back-a) blobs+EMD b) our method ground model, where m B is the number of images in the database B and n bg is an integer parameter. We then approximate c Q by the empirical quantile of order ǫ/m B of distances from Q to M j 's. That is, denoting by

B Q the n th bg closest image from Q among {M 1 , . . . , M n bg ×mB }, c Q is obtained as c Q = d(Q, B Q ).
The next point is to investigate the number of samples that is needed in order to get correct approximations of c Q . In order to do so, we perform the following experiment. Background databases of n bg × m B realizations of the background model are synthesized, for n bg = 1, . . . , 9. Then, for a query image Q, we count how many images from B are meaningful matches. Let r Q be this number of matches.

Figures 9 and10 show the standard deviation of r Q as a function of n bg , when considering ten background databases for each n bg . Standard deviations are divided by the average value of r Q for n bg = 1. Figure 9 deals with the IRHT database. This experiment suggests that, despite the high variability of realizations of the dead leaves model, n bg = 3 is enough to obtain consistent results with different background databases, and is therefore the numerical value retained in this paper. Figure 10 deals with CLIC database and suggests that n bg = 3 is not related to the specific IRHT content, and is also sufficient in the case of databases which are 10 times 10 bigger.

Variability of the matching threshold

In the following experiment, we computed all the 1500 possible IRHT queries Q with n bg = 3. The bidimensional histogram of (c Q , r Q ) is displayed in Figure 11 (recall that for a query Q, c Q is the matching threshold and r Q the number of matches). This experiment gives several insights into the proposed adaptive retrieval. The range of possible matching thresholds c Q is large, as well as the range of the number of meaningful matches r Q , which varies from 0 to 91 images. This shows that the matching threshold is quite different from a bound on distances or ranks and that it adapts to the query. This experiment also shows that r Q and c Q are not strongly correlated. A large r Q indicates series of images in the database like the one appearing in Figure 12. We see on the histogram that this can happen with low and high values of c Q . This example corresponds to the upper point in Figure 11, near (0.3, 90). In Figure 13, the same histogram is shown with the CLIC database, on 2000 different queries. In this case, there is a stronger correlation between r Q and c Q . Simultaneous high values for r Q and low values for c Q reveal the presence of very large series in this database which are mostly landscapes series.

Examples of queries

In the experiments, we display the query image, followed by meaningful matches from the database, in decreasing order up to the matching threshold indicated by a grey thumbnail followed by the first five images that do not match the query.

Retrieval results with the IRHT illuminated manuscripts database are displayed in Figures 14 to 22. In these examples, the number of meaningful matches r Q varies from 0 (second result of Figure 17) to 23 (Figure 16) according to the query and the the database. The matching threshold c Q sometimes stops after or before where an observer would stop the list of results, but results are most of the time very relevant. The matching threshold is a good indication of where the results should stop and it also informs the user on the quantity of relevant images present in the database.

In Figure 22, a retrieval is performed with a sketch query made with a simple painting tool. If a user roughly remembers a color composition, querying can also be performed this way.

Retrieval results with the CLIC photo database are shown in Figures 23 to 28. The original CLIC database is divided into approximately a hundred categories such as Landscape, City, Architecture. Many images from different categories can have the same spatial organization of colors and this appears on many results. For instance, in Figure 26, the query belongs to the category Architecture/Detail. The first six results belong to the category City/Building.

An example which illustrates how EMD works on thumbnails is shown in Figure 27 with a synthetic image from the CLIC database. The blue region at the center of the query moves a little on the closest result. On the last result, the blue region breaks apart into four small regions. Even though its composition does not respect the query, its highly similar color distribution allows larger pixel movements.

Increasing the α parameter would discard this image.

Even though the CLIC database is 10 times bigger than the IRHT database, the matching threshold still enables us to identify rare spatial organizations of colors. In Figure 28, from a query that is indeed quite unusual in the CLIC database, no match is returned. On the opposite, large series of similar composition can yield a high number of matches as in Figure 29, where a query has 72 meaningful matches. Only the first four matches and last five matches are visible here, as well as the first five images beyond the matching threshold.

Conclusion and future work

In this paper, we first introduced a method for the comparison of the spatial organization of colors between images, relying on a global optimization procedure through the use of the Earth Mover's Distance. The main claim here is that it is unnecessary to segment the image, a step which usually reduces the robustness of comparisons. Next, we apply the resulting distance to color image querying and address the decision stage of the system, that is the problem of deciding which images from a database match a query. This step is achieved by controlling the rate of false matchings using an a contrario approach. This work opens several perspectives. As mentioned earlier, we believe the method to be particularly adapted to the querying of cultural heritage databases. We are currently working with several artistic institutions to build large (typically over 100,000) specialized databases. We are also working on the installation of the system directly on web servers aiming at presenting the work of a collective artists studio [START_REF] Hurtut | A novel collaborative website and artworks database management system for artist-run centres[END_REF]. For such applications, it will be interesting to investigate further the implications of the various parameters involved in the definition of distance (3) or to develop distances specially adapted to, for example, black and white photographs. Next, several accelerations of EMD computations have been recently proposed [START_REF] Indyk | Fast image retrieval via embeddings[END_REF][START_REF] Assent | Approximation techniques for indexing the earth mover's distance in multimedia databases[END_REF][START_REF] Ling | An efficient earth mover's distance algorithm for robust histogram comparison[END_REF] and can be of interest for our approach. Another possible direction to speed up the query would be to impose constraints on the EMD between thumbnails, by forbidding large distances between pixels. This would yield sparse cost matrices and a faster computation of the assignment problem, see [START_REF] Burkard | Linear assignment problems and extensions[END_REF][START_REF] Burkard | Assignment Problems[END_REF]. Alternatively, one could adopt a coarse to fine strategy. For such a strategy, the solution of the optimization problem at a coarse scale would constrain the computation at a finer scale. We also started to test the approach on thumbnails having more information than just one color at each pixel, for instance a color array or texture features. In another direction, it is of interest to investigate the use of EMD to compare the color content of small local image patches, as it has been done recently [START_REF] Ling | An efficient earth mover's distance algorithm for robust histogram comparison[END_REF] for comparing the geometrical content of such local patches. 15. We see here that the matching threshold allows very variable number of matches according to the database content and the query. There are 23 meaningful matches. The database contains approximately a hundred white illuminations with a blue rectangle. On images that are not considered as matches, the blue rectangle has a different size and/or is not at the same location as in the query. Fig. 17. Same layout as in Figure 15. Two examples where the matching threshold provides very few meaningful matches. On the second example, it gives no result and indeed no such spatial organization can be encountered in the database.

Fig. 1 .

 1 Fig. 1. A coarsely sampled thumbnail. Left: source image (3200 × 2200 pixels). Right: the thumbnail (10 × 15 pixels).

Fig. 2 .

 2 Fig. 2. Two images (left) and the associated thumbnails (middle), made of n = 10 × 15 pixels. Each pixel of a thumbnail is associated to a feature in 5 dimensions (L, a, b, X, Y ).All of them have the same weight. In this context, EMD is equivalent to an assignment problem. Each pixel of the first thumbnail (top middle) is assigned to exactly one pixel of the second one (bottom middle). To illustrate the EMD operation, the thumbnail on the bottom right represents the optimal permutation of the pixels from the first thumbnail to match the pixels of the second one. Pixels of the red region are moved towards their location in the second image. The white column of paper is also correctly moved from right to left.

Fig. 3 .

 3 Fig. 3. Samples of the background model M using color marginals from the illuminated manuscript database used in the experimental section.

Fig. 4 .

 4 Fig. 4. Three query results with different thumbnail sizes (Left column, from top to bottom: n = 150, 600, 2400). The query image is framed in red and is followed by the three closest results. Results and their order slightly change according to the thumbnail size.

Fig. 5 .

 5 Fig. 5. Same layout as in Figure 4. Results and their order do not change in this case.

  Figure 8 illustrates this drawback of region-based methods. a) blobs+EMD b) our method

Fig. 7 .

 7 Fig. 7. Comparison between a regions + EMD method and our method. A query image (framed in red) followed by its three nearest matches is shown on the first lines of Figure a) and b), according to each method. The blob representations and the thumbnails (10 × 15) of images are shown on the second line of Figure a) and b) respectively. Here the blob segmentation is relatively satisfactory and both methods perform similarly.

Fig. 8 .

 8 Fig. 8. Same layout as in Figure 7. Here the blob segmentation fails to give a satisfactory representation of the color spatial organization of each image. Images are clearly under-segmented, yielding large regions with average colors of different patterns, leading to unsatisfactory matching.

Fig. 9 .

 9 Fig. 9. Standard deviations of the number of retrieved images for 9 queries on the IRHT database. Each colored line corresponds to a different query image. The retrieval is made on a database of m B images. The standard deviation of the rank of the last meaningful match (r Q ) is shown for background databases of size m B × n bg , for n bg = 1, . . . , 9. For each n bg , ten background databases have been generated in order to estimate this standard deviation. Standard deviations are divided by the average value of r Q for n bg = 1.

Fig. 10 .

 10 Fig. 10. Same layout as on Figure 9. Standard deviations of the number of retrieved images for 9 CLIC queries.

Fig. 11 .

 11 Fig. 11. Bidimensional histogram of (c Q , r Q ) over all the 1500 IRHT queries. c Q is the matching threshold and r Q the number of images closer to the query than c Q (meaningful matches).

Fig. 12 .

 12 Fig.[START_REF] Chen | An analysis of image retrieval tasks in the field of art history[END_REF]. The query and the first 9 matches (among a total of 90 matches) for the upper point in Figure11. The spatial organization of colors of the matches is very similar to the one of the query, denoting the presence of a large series of similar images in the database. The last thumbnail (framed in blue) is the background model thumbnail retained to compute the matching threshold.

Fig. 13 .

 13 Fig. 13. Bidimensional histogram of (c Q , r Q ) over 2000 CLIC queries. c Q is the EMD distance matching threshold, and r Q is the number of results closer to the query than c Q (meaningful matches).

Fig. 14 .

 14 Fig. 14. The query image is framed in red (top left) in figure a) and is followed (left from right and top to bottom) by its meaningful matches in decreasing order, up to the grey image which corresponds to the matching threshold. Beyond this limit, the next five nearest images are shown. The corresponding features (10 × 15 thumbnails) are shown in figure b). The last thumbnail (framed in blue) is the background model thumbnail retained to compute the matching threshold.

Fig. 15 .

 15 Fig. 15. The query image is framed in red (top left) in figure a) and is followed (left from right and top to bottom) by its meaningful matches in decreasing order, up to the grey image which corresponds to the matching threshold. Beyond this limit, the next five nearest images are shown.

Fig. 16 .

 16 Fig.[START_REF] Desolneux | From Gestalt Theory to Image Analysis: A Probabilistic Approach[END_REF]. Same layout as in Figure15. We see here that the matching threshold allows very variable number of matches according to the database content and the query. There are 23 meaningful matches. The database contains approximately a hundred white illuminations with a blue rectangle. On images that are not considered as matches, the blue rectangle has a different size and/or is not at the same location as in the query.

Fig. 18 .

 18 Fig. 18. Same layout as in Figure 15.

Fig. 19 .

 19 Fig. 19. Same layout as in Figure 15.

Fig. 20 .

 20 Fig. 20. Same layout as in Figure 15.

Fig. 21 .

 21 Fig. 21. Same layout as in Figure 15. Three examples with a small number of meaningful matches.

Fig. 22 .

 22 Fig.[START_REF] Ford | Solving the transportation problem[END_REF]. Same layout as in Figure15. Example of a retrieval from a sketch obtained with a simple painting software. In order for such a methodology to enable the retrieval of images we have a vague memory of, one can allow more flexibility on color through the parameter δ C of the distance de (see (3)).

Fig. 23 .

 23 Fig. 23. Same layout as in Figure 15 with a query from the CLIC database.

Fig. 24 .

 24 Fig. 24. Same layout as in Figure 15 with a query from the CLIC database.

It is possible to extend the definition to the case where total weights are not equal, see[START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF].

In this setting, the term background has no connection with background models used in object recognition and which refer to models of non-objects.

Code available at: http://elib.cs.berkeley.edu/blobworld/

Acknowledgments

We thank Gilles Kagan from IRHT for providing the IRHT database and for his helpful comments. We also thank Farida Cheriet from LIV4D for her support. This work was partially supported by a CNRS ACI supervised by Anne-Marie Edde and Dominique Poirel from IRHT, a grant from the LIV4D laboratory, and a FQRNT grant from the Quebec government.