Paolo Piro

Sandrine Anthoine

Eric Debreuve

Michel Barlaud

Sparse Multiscale Patches (SMP) for image categorization

In this paper we address the task of image categorization using a new similarity measure on the space of Sparse Multiscale Patches (SMP ). SMP s are based on a multiscale transform of the image and provide a global representation of its content. At each scale, the probability density function (pdf ) of the SMP s is used as a description of the relevant information. The closeness between two images is defined as a combination of Kullback-Leibler divergences between the pdfs of their SMP s. In the context of image categorization, we represent semantic categories by prototype images, which are defined as the centroids of the training clusters. Therefore any unlabeled image is classified by giving it the same label as the nearest prototype. Results obtained on ten categories from the Corel collection show the categorization accuracy of the SMP method.

Introduction

Image categorization is still one of the most challenging tasks in computer vision. It consists in labeling an image according to its semantic category. The main difficulty lies in using low-level information provided by digital images to retrieve semantic-level classes, which are generally characterized by high intra-variability.

Most approaches address the task of image categorization as a supervised learning problem. They use a set of annotated images to learn the categories and then assign to an unlabeled image one of these categories. These methods rely on extracting visual descriptors, which are to be highly specific, i.e. able to highlight visual patterns that characterize a category. Providing suitable sets of visual descriptors has been a topic of active research in the recent years and several approaches have been proposed.

Related works

Visual descriptors for image categorization generally consist of either global or local features. The former ones represent global information of images and are based on global image statistics such as color histograms [START_REF] Szummer | Indoor-outdoor image classification[END_REF] or edge directions histograms [START_REF] Vailaya | Image classification for contentbased indexing[END_REF]. Global feature-based methods were mostly designed to separate very general classes of images, such as indoor vs. outdoor scenes or city vs. landscape images. On the contrary, local descriptors extract information at specific image locations that are relevant to characterize the visual content. Local approaches are more adapted to cope with the intra-class variability of real image categories. Indeed these techniques are able to emphasize local patterns, which images of the same category are expected to share.

Early approaches based on local features work on image blocks. In [START_REF] Li | Automatic linguistic indexing of pictures by a statistical modeling approach[END_REF] color and texture features are extracted from image blocks to train a statistical model; this model takes into account spatial relations among blocks and across image resolutions. Then several region-based approaches have been proposed, which require segmentation of images into relevant regions. E.g. in [START_REF] Bi | A sparse support vector machine approach to regionbased image categorization[END_REF], an algorithm for learning region prototypes is proposed as well as a classification of regions based on Support Vector Machines (SVMs). More recent techniques have successfully used bags-of-features, which collect local features into variable length vectors. In [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF], the bags-of-features representing an image are spatial pyramid aggregating statistics of local features (e.g. "SIFT" descriptors); this approach takes into account approximate global geometric correspondences between local features. Other methods using bags-of-features are based on explicitly modeling the distribution of these vector sets. In fact, measuring the similarity between the bags-of-features' distributions is the main difficulty for this kind of categorization methods. For example, Gaussian Mixture Models (GMMs) have been used to model the distribution of bags of low-level features [START_REF] Liu | A similarity measure between unordered vector sets with application to image categorization[END_REF]. This approach requires both to estimate the model parameters and to compute a similarity measure to match the distributions.

Apart from what descriptors are used, computing a similarity measure between feature sets is crucial for most approaches. Measuring similarities is particularly adapted to categorization when one prototype of each category is defined in the feature space. In this context the similarity measure is used to find the prototype that best matches an unlabeled image. As pointed out in [START_REF] Zhang | Svm-knn: Discriminative nearest neighbor classification for visual category recognition[END_REF], this framework has provided the best results in image categorization. The main reason is that categorization based on similarity corresponds well to the way human beings recognize visual classes. Indeed human visual categories are mostly defined by similarity to prototype examples, as it results from research on cognitive psychology [START_REF] Rosch | Basic objects in natural categories[END_REF]. Taking advantage of these results, we propose a new categorization technique that is based on category prototypes and uses a similarity measure between feature set distributions.

Proposed statistical approach

The categorization method that we propose consists of two steps: the training step and the classification step. The training consists in selecting one prototype per category among a set of labeled images (training set). The classification step assigns to an unlabeled image (query image) the label of its closest prototype. Both training and query images are represented by their sets of Sparse Multiscale Patches (SMP s).

We have designed the SMP descriptors in order to exploit local multiscale properties of images, which generally convey relevant information about their category. Indeed SMP s describe local spatial structures of images at different scales by taking into account dependencies between multiscale coefficients across scale, space and color channels. An image is represented by the set of its SMP s at each scale and we measure the similarity between two images as the closeness between their SMP probability density functions (pdf ). Namely we combine their Kullback-Leibler divergence, which has already shown good performances in the context of image retrieval [START_REF] Puzicha | Empirical evaluation of dissimilarity measures for color and texture[END_REF].

We use this SMP similarity measure in both steps of the categorization, as depicted in the block diagram of Figure 1. During the training step we compute all pairwise SMP -based "distances" between images of the same training category C, thus obtaining a distance matrix for this category (see left column of Fig. 1). The entry s i,j of this matrix represents the distance between the image I i (query image) and the image I j (reference image); we denote this distance as S(I i |I j ). The distance matrix is used to select the prototype J C , which we define as the image minimizing the distance to all other images of the same category:

J C = argmin j|Ij ∈C i|Ii∈C s i,j (1)
Applying the same method to all training categories yields one prototype for each category in terms of SMP feature sets. 1In the second phase, which is the categorization, we compute the SMP similarity of an unlabeled image Q to all category prototypes (see right column of Fig. 1). The query image is given the same label as the nearest prototype.

Organization of the paper

In the rest of this paper we explain in more details how the SMP similarity measure is defined and give some experimental results of categorization. The SMP similarity measure is described in two steps. Firstly we define the proposed feature set of Sparse Multiscale Patches in Section 2. Secondly we define the similarity between SMP probability densities in Section 3. We also propose a method for estimating this measure non-parametrically, thus avoiding to model the underlying pdf s. Finally, in Section 4, we present some results of performing our categorization method on a subset of the Corel database.

Feature space: Sparse Multiscale Patches

Les us now define our feature space, which is based on a multiresolution decomposition of the images. Throughout the paper, we will denote by w(I) j,k the coefficient for image I at scale j and location in space k for a general multiresolution decomposition.

Multiscale patches

The wavelet transform enjoys several properties that have made it quite successful in signal processing and it naturally yields good candidates for image descriptors. Indeed, it provides a sparse representation of images, meaning that it concentrates the informational content of an image into few coefficients of large amplitude while the rest of the coefficients are small. Classical wavelet methods focus on these large coefficients and treat them separately, relying on their decorrelation, to efficiently process images. However, wavelet coefficients are not independent and these dependencies are the signature of structures present in the image. These dependencies have then been exploited in image enhancement (e.g [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF][START_REF] Romberg | Bayesian tree-structured image modeling using wavelet-domain hidden markov models[END_REF]). In particular, the authors of [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF] introduced the concept of patches of wavelet coefficients (called "neighborhoods of wavelet coefficients") to represent efficiently fine spatial structures in images.

Following these ideas, we define a feature space based on a sparse description of the image content by a multiresolution decomposition. More precisely, we group the Laplacian pyramid coefficients of the three color channels of image I into coherent sets called patches. Here the coherence is sought by grouping coefficients linked to a particular scale j and location k in the image.

In fact, the most significant dependencies are seen between a coefficient w(I) j,k and its closest neighbors in space: w(I) j,k±(0,1) , w(I) j,k±(1,0) and in scale: w(I) j-1,k , where scale j -1 is coarser than scale j. Grouping the closest neighbors in scale and space of the coefficient w(I) j,k in a vector, we obtain the patch → w(I) j,k (see Fig. 2):

→ w(I) j,k = w(I) j,k , w(I) j,k±(1,0) , w(I) j,k±(0,1) , w(I) j-1,k (2)
which describes the structure of the grayscale image I at scale j and location k.

The probability density functions of such patches at each scale j has proved to characterize fine spatial structures in grayscale images [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF][START_REF] Pierpaoli | Reconstructing sunyaev-zeldovich clusters in future cmb experiments[END_REF]. We consider color images in the luminance/chrominance space: I=(I Y ,I U ,I V ). Since the coefficients are correlated through channels, we aggregate in the patch the coefficients of the three channels:

w(I) j,k = → w(I Y ) j,k , → w(I U ) j,k , → w(I V ) j,k (3)
The low-frequency approximation that results from the Laplacian pyramid is also used to build additional feature vectors. The 3 × 3 pixel neighborhoods along the three channels are joined to form patches of dimension 27 (whereas patches from the higher-frequency subbands defined in Eq.( 3) are of dimension 18). The union of the higher-frequency and low-frequency patches forms our feature space. For convenience, the patches of this augmented feature space will still be denoted by w(I) j,k .

Sparse Multiscale Patches

The coefficients are obtained by a Laplacian pyramid decomposition [START_REF] Burt | The Laplacian pyramid as a compact image code[END_REF]. Indeed, critically sampled tensor wavelet transforms lack rotation and translation invariance and so would the patches made of such coefficients. Hence we prefer to use the Laplacian pyramid which shares the sparsity and inter/intrascale dependency properties with the wavelet transform while being more robust to rotations.

Multiscale coefficients provide a sparse representation of images and, similarly, patches of multiscale coefficients of large overall energy (sum of the square of all coefficients in a patch) also concentrate the information. Since the total number of patches in an image decomposition is 4/3N with N the number of pixels in the image, the number of samples we have in the feature space is quite large as far as measuring a similarity is concerned. The possibility of selecting a small number of patches which represent the whole set well is therefore highly desirable. In practice, we selected a fixed proportion of patches at each scale of the decomposition and proved that the resulting similarity measure (defined in Section 3) remains consistent (see [START_REF] Piro | Image retrieval via kullbackleibler divergence of patches of multiscale coefficients in the knn framework[END_REF] for details). This is exploited to speed up our computations. We now define a similarity on this feature space.

3 Similarity measure

Definition

Our goal is to define a similarity measure between two images I 1 and I 2 from their feature space i.e. from their respective set of patches {w(I 1 ) j,k } j,k and {w(I 2 ) j,k } j,k . When images are clearly similar (e.g. different views of the same scene, images containing similar objects...), their patches w(I 1 ) j l ,k l and w(I 2 ) j l ,k l do not necessarily correspond. Hence a measure comparing geometrically corresponding patches would not be robust to geometric transformations. Thus, we propose to compare the pdfs of these patches. Specifically, for an image I, we consider for each scale j the pdf p j (I) of the set of patches {w(I) j,k } k .

To compare two pdfs, we use the Kullback-Leibler (KL) divergence which derives from the Shannon differential entropy (quantifies the amount of information in a random variable through its pdf). The KL divergence (D kl ) is the quantity [START_REF] Puzicha | Empirical evaluation of dissimilarity measures for color and texture[END_REF]:

D kl (p 1 ||p 2 ) = p 1 log(p 1 /p 2 ). ( 4
)
This divergence has been successfully used for other applications in image processing in the pixel domain [START_REF] Boltz | High-dimensional kullback-leibler distance for region-of-interest tracking: Application to combining a soft geometric constraint with radiometry[END_REF][START_REF] Angelino | Image restoration using a knn-variant of the mean-shift[END_REF], as well as for evaluating the similarity between images using the marginal pdf of the wavelet coefficients [START_REF] Do | Wavelet based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF][START_REF] Wang | Qualityaware images[END_REF]. We propose to measure the similarity S(I 1 |I 2 ) between two images I 1 and I 2 by summing over scales the divergences between the pdfs p j (I 1 ) and p j (I 2 ) (with weights α j that may normalize the contribution of the different scales):

S(I 1 |I 2 ) = j α j D kl (p j (I 1 )||p j (I 2 )) (5)

Parametric approach to the estimation: pros and cons

The estimation of the similarity measure S consists of the evaluation of divergences between pdfs p j (I i ) of high dimension. This raises two problems. Firstly, estimating the KL divergence, even with a good estimate of the pdfs, is hard because this is an integral in high dimension involving unstable logarithm terms. Secondly, the accurate estimation of a pdf itself is difficult due to the lack of samples in high dimension (curse of dimensionality). The two problems should be embraced together to avoid cumulating both kinds of errors.

Parametrizing the shape of the pdf is not ideal. The KL divergence is easy to compute when it is an analytic function of the pdf parameters. This is the case for generalized Gaussians models, which fit well the marginal pdf of multiscale coefficients [START_REF] Do | Wavelet based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF][START_REF] Wang | Qualityaware images[END_REF]. However, this model cannot be extended to account for the correlations we want to exploit in multiscale patches. Mixture of Gaussians on contrary are efficient multidimensional models accounting for these correlations [START_REF] Pierpaoli | Reconstructing sunyaev-zeldovich clusters in future cmb experiments[END_REF] but the KL divergence is not an analytic function of the pdf parameters.

Thus, we propose to make no hypothesis on the pdf at hand. We therefore spare the cost of fitting model parameters but we have to estimate the divergences in this non-parametric context. Conceptually, we combine the Ahmad-Lin approximation of the entropies necessary to compute the divergences with "balloon estimate" of the pdfs using the kNN approach.

Proposed non-parametric estimation of the similarity measure

The KL divergence can be written as the difference between a cross-entropy H x and an entropy H (see Eq.( 4)):

H x (p 1 , p 2 ) = -p 1 log p 2 , H(p 1 ) = -p 1 log p 1 (6)
Let us explain how to estimate these terms from i.i.d sample sets

W i = {w i 1 , w i 2 , .., w i
Ni } of p i for i = 1 or 2 (The samples are in R d .) Assuming we have estimates p 1 , p 2 of the pdfs p 1 , p 2 , we use the Ahmad-Lin entropy estimators [START_REF] Ahmad | A nonparametricestimationof the entropy absolutely continuousdistributions[END_REF]:

H al x ( p 1 , p 2 ) = -1 N1 N1 n=1 log[ p 2 (w 1 n )], H al ( p 1 ) = -1 N1 N1 n=1 log[ p 1 (w 1 n )] (7)
to obtain a first estimator of the KL divergence:

D al kl ( p 1 , p 2 ) = 1 N1 N1 n=1 log[ p 1 (w 1 n )] -log[ p 2 (w 1 n )]. (8)
General non-parametric pdf estimators from samples can be written as a sum of kernels K with possibly varying bandwidth h (see [START_REF] Terrell | Variable kernel density estimation[END_REF] for a review):

p 1 (x) = 1 N1 N1 n=1 K h(W 1 ,x) (x -w 1 n ) (9)
We use a balloon estimator i.e. the bandwidth h(W 1 , x) = h W 1 (x) adapts to the point of estimation x given the sample set W 1 . The kernel is binary and the bandwidth is computed in the k-th nearest neighbor (kNN) framework [START_REF] Terrell | Variable kernel density estimation[END_REF]:

K h(W 1 ,x) (x -w 1 n ) = 1 v d ρ d k,W 1 (x) δ ||x -w 1 n || < ρ k,W 1 (x) (10)
with v d the volume of the unit sphere in R d and ρ k,W (x) the distance of x to its k-th nearest neighbor in W. This is a dual approach to the fixed size kernel methods and was firstly proposed in [START_REF] Loftsgaarden | A nonparametric estimate of a multivariate density function[END_REF]: the bandwidth adapts to the local sample density by letting the kernel contain exactly k neighbors of x among a given sample set.

Although this is a biased pdf estimator (it does not integrate to one), it has proved to be efficient for high-dimensional data [START_REF] Terrell | Variable kernel density estimation[END_REF]. Plugging Eq. [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF] in Eq.( 8), we obtain the following estimator of the KL divergence, which is valid in any dimension d and robust to the choice of k:

2 D kl (p 1 ||p 2 ) = log N2 N1-1 + d N1 N1 n=1 log[ρ k,W 2 (w 1 n )] -d N1 N1 n=1 log[ρ k,W 1 (w 1 n )] (11)
4 Experiments

Database and parameter settings

An experimental evaluation of the SMP method has been made on a subset of the Corel database. It includes 1,000 images of size 384 × 256 or 256 × 384 which are classified in ten semantic categories (Africa, Beach, Buildings, Buses, Dinosaurs, Flowers, Elephants, Horses, Food, Mountains). This dataset is well known in the domain of content-based image retrieval and particularly it has been widely used to evaluate several methods of image categorization, like in [START_REF] Chen | Image categorization by learning and reasoning with regions[END_REF], [START_REF] Li | A semantic region descriptor for local feature based image categorization[END_REF], [START_REF] Bi | A sparse support vector machine approach to regionbased image categorization[END_REF] and [START_REF] Zhu | Content-based image categorization and retrieval using neural networks[END_REF].

The SMP descriptors of images were extracted as described in Section 2. In particular, to build the patches, the Laplacian pyramid was computed for each channel of the image (in the YUV color space) with a 5-point binomial filter w5 = [1 4 6 4 1]/16, which is a computationally efficient approximation of the Gaussian filter classically used to build Laplacian pyramids. Two high-frequency subbands and the low-frequency approximation were used. In the following experiments, 1/16 (resp. 1/8 and all) of the patches were selected in the first high-frequency (resp. second high-frequency and low-frequency) subband to describe an image (see Section 2.1). At each scale, the KL divergence was estimated in the kNN framework, with k = 10. The contributions to the similarity measure of the divergences in all subbands were equally weighted (α j = 1 in Eq. ( 5)).

Since the computation of KL divergences is a time-consuming task, we developed a parallel implementation of the kNN search on a Graphic Processing Unit (GPU) [START_REF] Garcia | Fast k nearest neighbor search using gpu[END_REF]. This implementation is based on a brute-force approach and was run on a NVIDIA GeForce 8800 GTX GPU (with 768 MB of internal memory). The computation of one similarity measure between two images required 0.1 s on average.

Categorization results

In order to allow the prototype learning, we have randomly split the whole image set into a training set and a test set. We have tested our method for different sizes of the training set (15, 30 and 50 images per category), while using all the remaining images as a test set. For each value of the training set size we have repeated the same experiment ten times; hence all reported results are averaged over ten experiments with identical parameter values.

These results are shown in Fig. 3 in terms of the average recognition rate per category. This rate corresponds to the average proportion of correct classification as a function of the number of training images. We can see that the accuracy of categorization is quite stable with the size of the training set. In the same figure we can see an example of which images are selected as prototypes from the training set.

We have also measured the accuracy of our method by using a criterion for rejecting unreliable categorization results, so as to reduce the probability of misclassifications. Namely we have empirically fixed a threshold value for the query-to-prototype distance. Hence, whenever the SMP -based "distance" between the query image and its nearest prototype is larger than the threshold, we reject the result, thus giving no label to the image. Table 1 contains the confusion matrix which summarizes the categorization results for the case of 50 training images per category. The generic element a i,j of this matrix represents the average percentage of the test images belonging to category i which have been classified into category j. Therefore, the entries on the main diagonal show the categorization accuracy for each category. On the contrary, off-diagonal entries show classification errors. Interestingly note that big values of some of these entries are related to a certain semantic "intersection" between two different categories. This is particularly clear for misclassifications between categories 2 (Beaches) and 9 (Mountains), since several images belonging to these categories share similar visual patterns (e.g. rocky seaside pictures are visually similar to mountain sceneries containing lakes or rivers).

Our method shows good overall categorization performances, typically of the same order of magnitude as those of state-of-the-art methods developed by Chen and colleagues [START_REF] Chen | Image categorization by learning and reasoning with regions[END_REF][START_REF] Bi | A sparse support vector machine approach to regionbased image categorization[END_REF], although not as accurate on some categories". This can be explained by the difficulty of choosing a representative prototype amongst the training images. Indeed, for diverse categories, the classification rate may vary a lot according to the randomly chosen training set which yields to choosing different prototypes. For example, in the category "Africa", some training sets lead to a prototype being a picture of a face (as displayed in Fig. 3)

Fig. 1 .

 1 Fig. 1. Overview of the proposed method for image categorization.

Fig. 2 .

 2 Fig. 2. Building a patch of multiscale coefficients, for a single color channel image.

Note that the proximity matrix is not symmetric here si,j = sj,i.

Note that in the log term, N1 has been replaced by N1 -1, which corresponds to ommitting the current sample w 1 i in the set W 1 when estimating the entropy (Eq. (7) and (9)).

while others lead to village picture prototypes. Different methods are under study to tackle the problem of prototype selection, one of which being to create a poll of representative SMP s to define a pdf prototype. The patch selection strategy presented here is also basic: one retains a predefined percentage of the patches at each scale. The selection does not adapt to the specificity of the image considered and thus probably results in taking into account outliers of the underlying pdf s. Therefore we also study adaptive patch selection strategies. Considering the large prototype variations and the simplistic patch selection process used, the reasonably good performances obtained here show that the SMP similarity measure is promising regarding the classification step. Entry on the row i and column j is the average percentage (over the 10 experiments) of images belonging to the category i which have been classified into the category j (see Fig. 3 for category names). The last column lists the percentage of non-classified images for each category.

Conclusion

In this paper we tackled the task of image categorization by using a new image similarity framework. It is based on high-dimensional probability distributions of patches of multiscale coefficients which we call Sparse Multiscale Patches or SMP s. Image signatures are represented by sets of patches of subband coefficients, that take into account their intrascale, interscale and interchannel dependencies. The similarity between images is defined as the "closeness" between the distributions of their signatures, measured by the Kullback-Leibler divergence. The latter is estimated in a non-parametric framework, via a k-th nearest neighbor or kNN approach. Our approach to image classification is to represent each category by an image prototype. The latter is defined as the image minimizing the SMP -based measure with all other images in the category's training set.

The experiments we made on a subset of the Corel collection show that the SMP similarity measure is a promising tool for the categorization problem. The prototype selection as well as the patch selection are to be improved and are among the subjects of ongoing work.


Paolo Piro

Sandrine Anthoine

Eric Debreuve

Michel Barlaud

Sparse Multiscale Patches for Image Processing

Keywords: multiscale transform, sparsity, patches, Kullback-Leibler divergence, k-th nearest neighbor

This paper presents a framework to define an objective measure of the similarity (or dissimilarity) between two images for image processing. The problem is twofold: 1) define a set of features that capture the information contained in the image relevant for the given task and 2) define a similarity measure in this feature space. In this paper, we propose a feature space as well as a statistical measure on this space. Our feature space is based on a global descriptor of the image in a multiscale transformed domain. After decomposition into a Laplacian pyramid, the coefficients are arranged in intrascale/interscale/interchannel patches which reflect the dependencies of neighboring coefficients in presence of specific structures or textures. At each scale, the probability density function (pdf) of these patches is used as a descriptor of the relevant information. Because of the sparsity of the multiscale transform, the most significant patches, called Sparse Multiscale Patches (SMP), describe efficiently these pdfs. We propose a statistical measure (the Kullback-Leibler divergence) based on the comparison of these probability density functions. Interestingly, this measure is estimated via the nonparametric, k-th nearest neighbor framework without explicitly building the pdfs. This framework is applied to a query-by-example image retrieval method. Experiments on two publicly available databases showed the potential of our SMP approach for this task. In particular, it performed comparably to a SIFT -based retrieval method and two versions of a fuzzy segmentation-based method (the UFM and CLUE methods), and it exhibited some robustness to different geometric and radiometric deformations of the images.

Introduction

Similarity in image processing

Defining an objective measure of the similarity (or dissimilarity) between two images (or parts of them) is a recurrent question in image processing that is dealt with in quite different ways. When dealing with inverse problems such as denoising or deconvolution of images, a similarity measure is needed to evaluate how well the estimate explains the observations. However, for these problems, efforts have been concentrated in the conditioning of the inverse operator as well as the spatial properties of the estimated images. The measure of fitness to the data has been less studied and is usually a simple euclidean norm in pixel space such as:

d(I 1 , I 2 ) = i∈{pixel} |I 1 (i) -I 2 (i)| 2 .
At the other end of the spectrum, for some applications, the similarity measure is at the core of the problem and has received much more attention. This is the case for applications such as tracking or image retrieval, where the task is to rank the images of a database according to their visual similarity to the given query image. In any case, defining a similarity measure is a two steps process:

1. Define a set of properties that capture the information contained in the image relevant for the given task. This step defines the so-called feature space. 2. Define a similarity measure in the feature space. This measure is often (but not always) a distance.

Number of possibilities have been explored for the feature space itself. Some spaces involve a transform domain (e.g. wavelet transforms), some are based on various descriptors. A variety of descriptors (see [START_REF] Deselaers | Features for image retrieval: An experimental comparison[END_REF] for a review) has been proposed in the literature. Local descriptors (e.g. salient points [START_REF] Loupias | Wavelet-based salient points for image retrieval[END_REF]) are based on a subset of the pixels of the image while global descriptors give information about the image as a whole (e.g. color histograms [START_REF] Swain | Color indexing[END_REF]). Local descriptors exploit the information given by a limited of number of points of interest together with their spatial neighborhood. Hence much information in the image is not used in these methods (see [START_REF] Mikolajczyk | A performance evaluation of local descriptors[END_REF] for an extensive comparison and performance evaluation of most local descriptors). On the contrary, global descriptors include information of the whole image (e.g. histograms of intensity values). Global descriptors may be defined at the pixel level (e.g. color histograms [START_REF] Swain | Color indexing[END_REF]) and include no notion of spatial correlation or at the patch level including spatial and/or scale correlations. The concept of global patch descriptors is supported by statistical studies on images [START_REF] Huang | Statistics of natural images and models[END_REF]. Here, we propose a new descriptor of this kind. The similarity measure can range from simple euclidean norm to more sophisticated measures: robust estimators have been used for optical flow [START_REF] Black | A framework for the robust estimation of optical flow[END_REF][START_REF] Black | The robust estimation of multiple motions: parametric and piecewise-smooth flow fields[END_REF], Bhattacharya's distance for tracking [START_REF] Comaniciu | Real-time tracking of non-rigid objects using mean shift[END_REF], entropic measure such as entropy, Kullback-Leibler divergence or mutual information for registration [START_REF] Viola | Alignment by maximization of mutual information[END_REF][START_REF] Bansal | Entropy-based, multiple-portal-to-3dct registration for prostate radiotherapy using iteratively estimated segmentation[END_REF]. A general requirement for the similarity measure is the visual relevance, i.e. a strong correlation with human perception of similarity itself. Research in vision science has already brought some perspectives on how to do so [START_REF] Marques | On the potential of incorporating knowledge of human visual attention into cbir systems[END_REF]. Nevertheless, designing systems purely based on the perceptual characteristics of the human visual system is a difficult task. Therefore, once meaningful features have been selected, we prefer to employ metrics that have a mathematical foundation and can be easily implemented. For this purpose, several distance metrics have been used to compare feature vectors for various tasks of image processing. The authors of [START_REF] Puzicha | Empirical evaluation of dissimilarity measures for color and texture[END_REF] give a variety of such measures and empirically show how the selection of a metric affects the performances of a retrieval system.

In this paper we propose a feature space as well as a statistical measure on this space. Our feature space is based on a global descriptor in a transformed domain that we call Sparse Multiscale Patches. The measure we propose on this space is statistical: it compares the probability density function (pdf) of these patches.

Proposed feature space and measure

We propose a new descriptor based on Sparse Multiscale Patches. In short, we integrate using probability distributions the local information brought by the SMP. The key aspects of these descriptors are the following:

-A multiscale representation of the images; -Inter/intrascale patches that describe locally the structure of the image at a given scale; -A sparse repartition: most of the energy is concentrated in a few patches.

Note that the occurrence in different parts of an image of similar patches of spatially neighboring pixels has been exploited in image processing [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF][START_REF] Awate | Unsupervised, information-theoretic, adaptive image filtering for image restoration[END_REF][START_REF] Angelino | Image restoration using a knn-variant of the mean-shift[END_REF].

Here the concept is used for multiscale coefficients as proposed in [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF].

The visual content of images is represented by patches of multiresolution coefficients. The extracted feature vectors are viewed as samples from an unknown multidimensional distribution. The multiresolution transform of an image being sparse, a reduced number of patches yields a good approximation of the distribution. We estimate the similarity between images by a pseudo-distance (or measure) between these multidimensional probability density functions.

We propose to use the Kullback-Leibler (KL) divergence as a similarity measure that quantifies the closeness between two probability density functions. Such measure has already shown good performances in the context of image retrieval [START_REF] Puzicha | Empirical evaluation of dissimilarity measures for color and texture[END_REF]. It has already been used for the simple case of parametrized marginal distributions of wavelet coefficients [START_REF] Do | Wavelet based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF][START_REF] Wang | Qualityaware images[END_REF], assuming independence of the coefficients. In contrast, we define multidimensional feature vectors (patches), that capture interscale and intrascale dependencies among subband coefficients. These are better adapted to the descriptor of local image structures and texture. In addition, for color images, we take into account the dependencies among the three color channels; hence patches of coefficients are also interchannel. This approach implies to estimate distributions in a high-dimensional statistical space, where fixed size kernel options to estimate distributions or divergences fail. Alternatively, we propose to estimate the KL divergence directly from the samples by using the k-th nearest neighbor (kNN) approach, i.e. adapting to the local sample density.

Organization of the paper

This paper is organized as follows. In Section 2, we define our feature space which consists of inter/intrascale and interchannel patches of Laplacian pyramid coefficients for color images, called Sparse Multiscale Patches. We then define the global similarity on this feature space in Section 3 by combining similarities between the probability density functions of these patches at different scales. The comparison between pdfs is measured by the KL divergence. We also explain how to estimate this quantity. Finally, in the last section we illustrate the use of the proposed measure in a particular application: image retrieval.

2 Feature space: Sparse Multiscale Patches Throughout this paper, we will denote the input image by I, the scale of the multiresolution decomposition by j, and the location in the 2D image space by k.

Multiscale coefficients: strengths and limits

The wavelet transform enjoys several properties that have made it quite successful in signal processing and that are relevant for the definition of similarity between images. Indeed, it provides a sparse representation of images, meaning that it concentrates the informational content of an image into few coefficients of large amplitude while the rest of the coefficients are small. This combined with a fast transform is what makes wavelet thresholding methods so powerful: in fact just identifying large coefficients is sufficient to extract where the information lies in the image. For example, denoising can be done very efficiently by simply thresholding wavelet coefficients as proved in [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF]. Such simple coefficient-wise treatments provide results of excellent quality at a reduced computational cost.

In fact, these classical wavelet methods treat each coefficient separately, relying on the fact that they are decorrelated. However, the wavelet coefficients are not independent and these dependencies are the signature of structures present in the image. For example, a discontinuity between smooth regions at point k 0 will give large coefficients at this point at all scales j (w(I) j,k0 large for all j). Classical methods using coefficient-wise treatments may destroy these dependencies between coefficients and hence alter the local structure of images. Therefore models using the dependencies between coefficients have been proposed and used in image enhancement (e.g. [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF][START_REF] Romberg | Bayesian tree-structured image modeling using wavelet-domain hidden markov models[END_REF]). In particular, the authors of [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF] introduced the concept of patches of wavelet coefficients (which they called "neighborhoods of wavelet coefficient") to represent efficiently fine spatial structures in images.

Multiscale patches for color images

Following these ideas, we define a feature space based on a sparse descriptor of the image content by a multiresolution decomposition. More precisely, we group the Laplacian pyramid coefficients of the three color channels of the image I into coherent sets called patches. Here the coherence is sought by grouping coefficients linked to a particular scale j and location k in the image.

In fact, the most significant dependencies are seen between a coefficient w(I) j,k and its closest neighbors in space: w(I) j,k±(0,1) , w(I) j,k±(1,0) and in scale: w(I) j-1,k , where j -1 represents the scale a step coarser than the scale j. Grouping the closest neighbors in scale and space of the coefficient w(I) j,k in a vector, we obtain the patch → w(I) j,k (see Fig. 1): → w(I) j,k = w(I) j,k , w(I) j,k±(1,0) , w(I) j,k±(0,1) , w(I) j-1,k

which describes the structure of the grayscale image I at scale j and location k.

The probability density functions of such patches at each scale j have proved to characterize fine spatial structures in grayscale images [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF][START_REF] Pierpaoli | Reconstructing sunyaev-zeldovich clusters in future cmb experiments[END_REF]. Such patches are therefore relevant features for our problem as will be seen in Section 4.3.

We consider color images in the luminance/chrominance space: I = (I Y , I U , I V ). Since the coefficients are correlated through channels, we aggregate in the patch the coefficients of the three channels:

w(I) j,k = → w(I Y ) j,k , → w(I U ) j,k , → w(I V ) j,k (2)
with

→ w(I Y ) j,k , → w(I U ) j,k and → w(I V ) j,k
given by Eq.( 1). The low-frequency approximation that results from the Laplacian pyramid is also used to build additional feature vectors. Namely, 3 × 3 pixel neighborhoods along all three channels are joined together to form patches of dimension 27 (whereas patches from the higher-frequency subbands are of dimension 18, as defined in Eq.( 2)). The union of the higher-frequency and low-frequency patches forms our feature space. The patches of this augmented feature space will be denoted by w(I) j,k .

Multiscale transform

The coefficients are obtained by a Laplacian pyramid decomposition [START_REF] Burt | The Laplacian pyramid as a compact image code[END_REF]. Indeed, critically sampled tensor wavelet transforms lack rotation and translation invariance and so would the patches made of such coefficients. Hence we prefer to use the Laplacian pyramid which shares the sparsity and inter/intrascale dependency properties with the wavelet transform while being more robust to rotations. Moreover, the Laplacian pyramid is at the basis of the SVC standard and thus our approach will be compatible with SVC.

Sparsity of the multiscale patches

As we have seen earlier, multiscale coefficients provide a sparse representation of images by concentrating the information into a few coefficients of large amplitude and this sparsity is exploited on the raw coefficients in thresholding methods. As illustrated in Fig. 2, our experiments show that the patches of multiscale coefficients of large overall energy (sum of the square of all coefficients in a patch) also concentrate the information. More precisely, we selected a fixed proportion of patches at each scale of the decomposition and proved that the resulting similarity measure (defined in section 3) remains consistent (see [START_REF] Piro | Image retrieval via kullbackleibler divergence of patches of multiscale coefficients in the knn framework[END_REF] for details). Since the total number of patches is 4/3N with N the number of pixel in an image, the number of samples we have in the feature space is quite large as far as measuring a similarity is concerned. The possibility of selecting a small number of patches which represent the whole set well is highly desirable and we will exploit it to speed up our computations.

Note that other selecting procedures may be investigated such as using the energy of the central coefficient, using the sum of absolute differences in the patches or thresholding based on the variance of the patches.

Let us now explain how we define a similarity in this feature space.

3 Similarity measure

Definition

Our goal is to define a similarity measure between two images I 1 and I 2 from their feature space i.e. from their respective set of patches {w(I 1 ) j,k } j,k and {w(I 2 ) j,k } j,k . When images are clearly similar (e.g. different views of the same scene, images containing similar objects...), their patches w(I 1 ) jo,ko and w(I 2 ) jo,ko do not necessarily correspond. Hence a measure comparing geometrically corresponding patches would not be robust to geometric transformations. Thus, we propose to compare the pdfs of these patches. Specifically, for an image I, we consider for each scale j the pdf p j (I) of the set of patches {w(I) j,k } k .

To compare two pdfs, we place ourselves in the framework of Bregman divergences, which allows to generate a class of metrics that generalize the classical squared Euclidean distance. These divergences do not necessarily satisfy the triangle inequality nor are symmetric (they are not metrics) but share similar properties. A Bregman divergence is derived from a convex function. For example, the square Euclidean distance stems from the square function f (x) = x 2 , while the Kullback-Leibler divergence derives from the function f (x) = x log x [START_REF] Nielsen | On bregman voronoi diagrams[END_REF][START_REF] Nielsen | On the smallest enclosing information disk[END_REF]. In this paper, we use the Kullback-Leibler divergence because it is a Bregman divergence that derives from the Shannon differential entropy (quantifies the amount of information in a random variable through its pdf). The Kullback-Leibler divergence (D kl ) is the following quantity [START_REF] Puzicha | Empirical evaluation of dissimilarity measures for color and texture[END_REF]:

D kl (p 1 ||p 2 ) = p 1 log(p 1 /p 2 ), (3)
This divergence has been successfully used for other applications in image processing in the pixel domain [START_REF] Boltz | High-dimensional kullback-leibler distance for region-of-interest tracking: Application to combining a soft geometric constraint with radiometry[END_REF][START_REF] Angelino | Image restoration using a knn-variant of the mean-shift[END_REF], as well as for evaluating the similarity between images using the marginal pdf of the wavelet coefficients [START_REF] Do | Wavelet based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF][START_REF] Wang | Qualityaware images[END_REF]. In this paper, we propose to measure the similarity S(I 1 , I 2 ) between two images I 1 and I 2 by summing over scales the divergences between the pdfs p j (I 1 ) and p j (I 2 ):

S(I 1 , I 2 ) = j α j D kl (p j (I 1 )||p j (I 2 )) (4)
where α j is a positive weight that may normalize the contribution of the different scales.

Limits of the parametric approaches to the estimation

The estimation of the similarity measure S consists of the evaluation of divergences between pdfs p j (I i ) of high dimension. This raises two problems. Firstly, estimating the KL divergence, even with a good estimate of the pdfs, is hard because this is an integral in high dimension involving unstable logarithm terms. Secondly, the accurate estimation of a pdf itself is difficult due to the lack of samples in high dimension (curse of dimensionality). The two problems should be embraced together to avoid cumulating both kinds of errors.

A first idea consists in parametrizing the shape of the pdf. The marginal pdf of multiscale coefficients is well modeled by generalized Gaussians. In this case, the KL divergence is an analytic function of the pdf parameters. This technique has been used in [START_REF] Do | Wavelet based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF][START_REF] Wang | Qualityaware images[END_REF] to compare images on the basis of the marginal pdf of their wavelet coefficients. To our knowledge, the generalized Gaussian model cannot be extended to account for correlations in higher dimension. Mixture of Gaussians on the other hand are efficient multidimensional models accounting for correlations that fit well the pdf of wavelet coefficients patches [START_REF] Pierpaoli | Reconstructing sunyaev-zeldovich clusters in future cmb experiments[END_REF]. However the KL divergence is not an analytic function of the model parameters.

Thus, we propose to make no hypothesis on the pdf at hand. We therefore spare the cost of fitting model parameters but we have to estimate the divergences in this non-parametric context. Conceptually, we combine the Ahmad-Lin approximation of the entropies necessary to compute the divergences with "balloon estimates" of the pdfs using the kNN approach.

Non-parametric estimation of the similarity measure

The KL divergence can be written as the difference between a cross-entropy H x and an entropy H (see Eq.( 3)):

H x (p 1 , p 2 ) = -p 1 log p 2 , H(p 1 ) = -p 1 log p 1 . (5)
Let us explain how to estimate these terms from an i.i.d sample set

W 1 = {w 1 1 , w 1 2 , .., w 1 N1 } of p 1 and an i.i.d sample set W 2 = {w 2 1 , w 2 2 , .., w 2 N2 } of p 2 . (The samples are in R d .)
Assuming we have estimates p 1 , p 2 of the pdfs p 1 , p 2 , we use the Ahmad-Lin entropy estimators [START_REF] Ahmad | A nonparametricestimationof the entropy absolutely continuousdistributions[END_REF]:

H al x ( p 1 , p 2 ) = -1 N1 N1 n=1 log[ p 2 (w 1 n )], H al ( p 1 ) = -1 N1 N1 n=1 log[ p 1 (w 1 n )]. (6)
General non-parametric pdf estimators from samples can be written as a sum of kernels K with (possibly varying) bandwidth h (see [START_REF] Terrell | Variable kernel density estimation[END_REF] for a review):

p 1 (x) = 1 N1 N1 n=1 K h(W 1 ,x) (x -w 1 n ). (7)
-Parzen estimators h(W 1 , x) = h: the bandwidth is constant. They perform very well with samples in one dimension but become unadapted in high dimension due to the sparsity of the samples: the trade-off between a bandwidth large enough to perform well in low local sample density (which may oversmooth the estimator) and a bandwidth small enough to preserve local statistical variability (which may result in an unstable estimator) cannot always be achieved. To cope with this problem, kernel estimators using adaptive bandwidth have been proposed;

-Sample point estimators h(W 1 , x) = h W 1 (w 1 i ), i ∈ {1
, N 1 }: the bandwidth adapts to each sample w 1 i given the sample set W 1 ; -Balloon estimators h(W 1 , x) = h W 1 (x): the bandwidth adapts to the point of estimation x given the sample set W 1 .

We use a balloon estimator with a binary kernel and a bandwidth computed in the k-th nearest neighbor (kNN) framework [START_REF] Terrell | Variable kernel density estimation[END_REF]. This is a dual approach to the fixed size kernel methods and was firstly proposed in [START_REF] Loftsgaarden | A nonparametric estimate of a multivariate density function[END_REF]: the bandwidth adapts to the local sample density by letting the kernel contain exactly k neighbors of x among a given sample set:

K h(x) (x -w 1 n ) = 1 v d ρ d k,W 1 (x) δ ||x -w 1 n || < ρ k,W 1 (x) (8)
with v d the volume of the unit sphere in R d and ρ k,W (x) the distance of x to its k-th nearest neighbor in W. Although this is a biased pdf estimator (it does not integrate to one), it has proved to be efficient for high-dimensional data [START_REF] Terrell | Variable kernel density estimation[END_REF]. Plugging Eq.( 8) in Eq.( 6), we obtain the following estimators of the (cross-)entropy:

H knn ( p 1 ) = log[N 1 v d ] -log(k) + d N 1 N1 n=1 log[ρ k,W 1 (w 1 n )] , (9)
H knn x ( p 1 , p 2 ) = log[N 2 v d ] -log(k) + d N 1 N1 n=1 log[ρ k,W 2 (w 1 n )] . (10)
As previously stated, these estimates are biased. A correction of the bias has been derived in [START_REF] Goria | A new class of random vector entropy estimators and its applications in testing statistical hypotheses[END_REF] in a different context. In the non-biased estimators of the (cross)-entropy the digamma function ψ(k) replaces the log(k) term:

H knn ( p 1 ) = log[(N 1 -1)v d ] -ψ(k) + d N 1 N1 n=1 log[ρ k,W 1 (w 1 n )] , (11)
H knn x ( p 1 , p 2 ) = log[N 2 v d ] -ψ(k) + d N 1 N1 n=1 log[ρ k,W 2 (w 1 n )] . (12)
And hence the KL divergence reads:

D kl (p 1 ||p 2 ) = log N2 N1-1 + d N1 N1 n=1 log[ρ k,W 2 (w 1 n )]-d N1 N1 n=1 log[ρ k,W 1 (w 1 n )]. ( 13
)
This estimator is valid in any dimension d and robust to the choice of k.

4 Application: Image Retrieval

Content-based image retrieval

With the rapid growing of general-purpose image collections, performing efficiently a search on such large datasets becomes a more and more critical task.

Content-based image retrieval (CBIR) systems tackle this task by analyzing the content of images in order to provide meaningful signatures of them. Automatic search of the target images is made possible by defining a similarity measure on the underlying signature space which has a reduced dimension. As a result, content based image retrieval mainly relies on describing the image content in a relevant way (the feature space) and defining a quantitative measure on this space (the similarity measure): the retrieval task is then accomplished by ranking images in increasing order of the pseudo-distance between their feature vector and the one of a given query image.

As seen in the introduction, a variety of descriptors and similarity measures have been proposed to represent image content. In this paper, we will compare our SMP approach to three different approaches to image retrieval, two of them relying on the same similarity measure. The first approach is based on SIFT descriptors [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], which are considered state-of-the-art amongst local descriptors. Salient points are extracted by detecting the highest coefficients in the wavelet transform of the image and SIFT features are then represented by histograms of the gradient orientation in regions of interest. Matching the SIFT features obtained in two images allows then to quantify their similarity. The other methods to which we compared ours use a segmentation-based fuzzy logic approach called UFM for Unified Feature Matching [START_REF] Chen | A region-based fuzzy feature matching approach to contentbased image retrieval[END_REF]. The descriptors are fuzzy features (called fuzzy sets) reflecting the color, texture, and shape of each segmented region. The UFM measure then integrates the fuzzy properties of all the regions to quantify the similarity between two images. Using this measure, the authors proposed two image retrieval algorithms. The first one is a strictly content-based approach (similarly to ours): it consists in ranking the database images based solely on their UFM distance to the query. We refer to it as the UFM approach. The retrieval accuracy is improved by a second method called CLUE : the UFM distances between target images themselves are used to obtain a clustering of the data from which the ranking is obtained. Consequently, this method involves additional information compared to strict content-based systems such as our approach.

Database and parameter settings

Databases

Numerical experiments were performed on two different databases. The first one contains small categories and allows to evaluate specific performances of a retrieval system such as its robustness to deformations; while the second database, with larger categories, allows to test global retrieval performances.

One of these databases contains 1,000 images of the Nister Recognition Benchmark collection [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF]. The images of size 640x480 pixels are grouped by sets of four images showing the same scene or object. Their content is quite various, from indoor scenes with a single object to outdoor scenes. Images belonging to the same group are related by geometric deformations (rotation, translation, zoom and perspective) as well as radiometric deformations (changes of bright-ness and contrast). The ground-truth for any query image is clear: exactly the three other images of the same group are relevant.

The SMP retrieval method was also tested on a general-purpose image database from COREL that has been widely used for CBIR evaluation purposes. In particular, results presented in [START_REF] Chen | Clue: Cluster-based retrieval of images by unsupervised learning[END_REF] can be considered as a reference. We used the same subset of the COREL database as in [START_REF] Chen | Clue: Cluster-based retrieval of images by unsupervised learning[END_REF]. It includes 1,000 images of size 384 × 256 or 256 × 384 which are classified in 10 semantic categories (Africa, Beach, Buildings, Buses, Dinosaurs, Flowers, Elephants, Horses, Food, Mountains). In some categories, the visual similarity between two given images is not always obvious since the grouping has been made in a semantic sense (e.g., category "Africa").

Parameter settings

To build the patches as defined in section 2.2, the Laplacian pyramid was computed for each channel of the image (in the YUV color space) with a 5-point binomial filter w 5 = [1 4 6 4 1]/16, which is a computationally efficient approximation of the Gaussian filter classically used to build Laplacian pyramids. Three high-frequency subbands plus the low-frequency approximation were used.

In the following experiments, 1/64 (resp. 1/32, 1/16 and all) of the patches were selected in the first high-frequency (resp. second, third and low-frequency) subband to describe an image (see 2.4). At each scale, the KL divergence was estimated in the kNN framework, with k = 10. The contributions to the similarity measure from the divergences of all subbands were equally weighted (α j = 1 in Eq.( 4)).

Note that the use of the Jensen-Shannon divergence, which is a symmetrized version of the KL divergence, has also been studied. We found that the performances of this symmetric measure are less good than with the KL divergence, and so until further understanding of this phenomenum, we report here only the results with the KL divergence.

Numerical experiments

This section presents an experimental analysis of the SMP method; the patchbased retrieval algorithm is evaluated in terms of its ability to retrieve similar images in a query-by-example context. Images belonging to the Nister database were used to evaluate the robustness of the method to different geometric transformations. A set of artificially-degraded images of this database was also used to evaluate the retrieval performances with respect to radiometric deformations (JPEG2000 compression noise). The global retrieval performances on the Nister database were evaluated by ROC (Receiver Operating Characteristic) curves and our method was compared to a reference SIFT -based retrieval algorithm. For the COREL database, the global retrieval performances were evaluated by precision curves and our method was compared with the fuzzy, segmentationbased UFM approach. Note that for all the following experiments, the given distance between images is S (Eq. ( 4)), hence the smaller the given distance is the more similar the two considered images are.

Robustness to geometric deformations

The robustness of a retrieval system to geometric deformations is its ability to find relevant images in spite of some transformations of the query, such as changes of viewpoint, rotations, zoom. This is an important requirement in image retrieval, e.g. for finding a given object in different images independently of the viewpoint. Because of its structure, the Nister database allows to evaluate the robustness of the proposed method to geometric deformations. Indeed, the database is composed of groups of four images containing the same object or scene under different viewpoints and/or lightening conditions.

Examples of retrieval for five query images taken from the database are presented in Fig. 3. In this figure, each row displays the retrieval result for the query image shown in the leftmost column. From the second column on, one can see the first 4 retrieved images ranked in increasing order of their distance to the query. Hence the second leftmost image is the most similar one, excluding the query image which is always ranked first with a distance of zero. The first retrieved images are generally relevant for the query, in spite of rotations (row 2), changes of viewpoint (rows 1, 3, 5) and zooms (rows 2, 4). This shows that the proposed descriptors and similarity measure are robust in terms of geometric deformations for the retrieval problem.

Robustness to JPEG2000 compression

Another important requirement for content-based retrieval systems is the robustness to radiometric deformations. Transmission on heterogeneous networks requires compression. This process induces a loss of quality that can be significant especially in critical transmission conditions. A retrieval system is expected to be robust to compression quality. To test the proposed method on this specific point, groups of images from the Nister database were expanded. Namely, three highly-compressed versions of one image were added to each group. They were obtained by setting three different quality levels of JPEG2000 compression.

Queries were launched on this dataset with both original and compressed images. An example of the results is shown in Fig. 4, where a non-compressed image is used as a query. The distance from the three compressed versions to the query image being quite small, the system ranked them first and before any geometrically deformed version of the query. This behavior is general and still holds when compressed images are used as queries, confirming the reliability of the proposed similarity measure in terms of its robustness to compression. Moreover, the distance to the query increases as the compression level increases. This is shown in Fig. 4, where images A, B, C are compressed versions of the query image in decreasing order of quality, the PSNR being respectively of 31.8, 29.7 and 29.3 dB.

Image retrieval performances (I): ROC curves and comparison with a SIFT -based method

The overall performances of the SMP retrieval method were evaluated by ana- performance is still satisfactory for a larger dataset. In any case, the best tradeoff between precision and recall was reached when we retrieved three images, i.e. when the cut-off value matches exactly the number of relevant images; as a result, there is a high probability that the retrieved images are all and only the relevant ones.

Finally, the results for our SMP and the SIFT -based approach are shown in Fig. 6. The latter were obtained by running a publicly available Matlab implementation of the SIFT algorithm [START_REF] Lowe | Sift keypoint detector[END_REF]. Because of the long processing time of the SIFT implementation (4.8 s on average for each comparison between two images), performing a query with each image of the database could not be done in a reasonable time. In consequence, a comparison was made by querying a subset of 100 images. In light of the ROC curves, the performances of our SMP method and the SIFT -based algorithm are comparable for this experiment.

Image retrieval performances (II): precision curve and comparison with the UFM method

The SMP retrieval method was also tested on a subset of the COREL database and compared to the UFM and CLUE methods [START_REF] Chen | Clue: Cluster-based retrieval of images by unsupervised learning[END_REF]. This database is made of a small number of categories [START_REF] Viola | Alignment by maximization of mutual information[END_REF] containing a large number of images per category (100). Hence, ROC curves are not adapted to evaluate the global performances of a retrieval system in this case. Instead, we used the Average Precision to evaluate the retrieval performances for each category (the precision values for a cut-off equal to 100 were averaged over all images of the category) as in [START_REF] Chen | Clue: Cluster-based retrieval of images by unsupervised learning[END_REF].

Examples of our retrieval results are shown in Fig. 7 and the Average Precision is given for each category in Fig. 8 (dark blue bars). The results of the UFM and CLUE approaches are also displayed in this latter figure for comparison. Fig. 7 illustrates the fact that the most of time, the first four retrieved images belong to the query's category (row 1, 4, and 5). This figure also illustrates well the difficulties encountered in this task: since the categories are quite large and diverse, images belonging to different categories may have very similar visual properties that are picked by our method. For example, the elephant and building (row 2 of Fig. 7) have dominating vertical structures and same dominant colors. Likewise, images belonging the "mountains" or "beaches" are freqently mismatched (row 3 of Fig. 7). These retrieval errors are common to all methods comparing images solely on the basis of the image content (i.e. introducing no semantics) and explain the fluctuation of the results displayed in Fig. 8 for all three methods. Our method compares well with the two established methods displayed here: it is more accurate than UFM (gray bars) for six categories out of ten; the accuracy is also better than or comparable to CLUE (white bars) for five categories out of ten. On average, our method performs better than the UFM approach and slightly less well than the CLUE one. As pointed out in Section 4.1, the SMP and UFM approaches are strictly content-based approaches. The CLUE method, while performing better, uses additional image distances and is therefore much more time-consuming. Thus, the performances of our method seem quite promising for three reasons: -It performs slightly better than the UFM approach which relies on the same information. -The results are not far from those of the more advanced CLUE approach which relies on more information. -A similar clustering processing as the one applied with the UFM measure in CLUE may be applied to improve the SMP approach.

In conclusion, in its current state of development, the proposed SMP measure does not outperform the state-of-the-art methods selected as benchmark here. However, it does bring a novel approach to tackle the problem of image retrieval.

Computational speed-up(s)

The evaluation of our SMP similarity requires the computation of several KL divergences in a non-parametric framework. Since this is a time-consuming task, we propose two ways to speed-up the computations. The first one is based on a GPU implementation of the algorithm, the second on a preselection of the relevant images in the database.

GPU implementation

When computing the similarity between two images with the SMP approach, most of the time is devoted to the search of the k-th nearest neighbors in the evaluation of the KL divergences. Indeed, finding a k-th nearest neighbor requires to compute and sort distances between features (here the patches). The "brute force" algorithm has a complexity of order O(N 2 ) for N samples in the feature For each retrieved image, the SMP similarity measure to the query is also shown.

set. Smarter algorithms with a lower complexity (typically of order O(N log N )) such as the classical KD-tree (ANN) algorithm [START_REF] Arya | An optimal algorithm for approximate nearest neighbor searching fixed dimensions[END_REF] have been designed. Nevertheless, in practice, the computation time of a similarity between two images with the SMP approach remains large even with this low-complexity algorithm: on average 2.2s on a Pentium 4 3.4 GHz (2GB of DDR memory) with the ANN algorithm.

To speed up the computation time, we developed a parallel implementation of the kNN search on a Graphic Processing Unit (GPU) [START_REF] Garcia | Fast k nearest neighbor search using gpu[END_REF] using CUDA. This implementation is based on a brute force approach since recursive algorithms (the preferred strategy when using trees such as in ANN) are not parallelizable. It was implemented on an NVIDIA GeForce 8800 GTX card with 768 MB of memory. The computation time for one similarity measure between two images required 0.2s on average (i.e., 10 times less than with the CPU implementation of ANN). As of today, the brute force algorithm parallelized on GPU is by far the fastest implementation of our method. Developing smart algorithms (such as the KD-tree one), which may not be parallelizable but have a very low complexity, is a topic of active research, as is the development of GPU for computational purposes. Hence both types of methods should be kept in mind for efficient implementations in the near future.

Preselection of the relevant images

The computational speed can be improved by splitting the retrieval task into two steps:

1. Only the low frequency contribution to the similarity measure defined in Eq. ( 4) is computed for all images in the database. This "partial" similarity measure produces a first ranking of the database images from which the first n images are selected for the next step. 2. The complete similarity measure is computed between the query and the n selected images.

This procedure saves computation time as it computes the whole similarity measure only for a reduced number of images (computing only part of it for images that are unlikely to be relevant to the query). The smaller the size of the preselected subset, the greater the improvement in terms of computation time. For example, when a query on the Nister database is processed following the described two-step procedure with a selected subset of 50 images, the average computation time per image drops from 0.2s to about 0.06s with the GPU implementation (and with similar retrieval performances). It is clear however that the number of preselected images cannot be arbitrarily small without seriously affecting retrieval performances. It should be large enough compared to the number of images in the database as well as the number of relevant images for the query.

Conclusion

In this paper, we proposed a new image similarity framework based on highdimensional probability distributions of patches of multiscale coefficients which we call Sparse Multiscale Patches. Feature sets are represented by these patches of subband coefficients that take into account intrascale, interscale and interchannel dependencies. The similarity between two images was defined as a linear combination of the "closeness" between the distributions of their features at each scale measured by the Kullback-Leibler divergence. The Kullback-Leibler divergences are estimated in a non-parametric framework, via a kNN approach. The proposed similarity measure seems to be stable when selecting a reduced number of patches, proving that a few significant patches are enough to represent the image features. This is a consequence of the sparsity of the multiscale transform.

We applied this framework to image retrieval. The proposed approach takes advantage of the properties of its global multiscale descriptors. In particular, it is robust to JPEG2000 compression (i.e. it matches the visual similarity between images with different amounts of blur or compression noise). Retrieval experiments were conducted on two publicly available datasets of real world images (Nister Recognition Benchmark and the COREL database) to evaluate the average performances of the method. In particular, the Nister dataset was used to benchmark the robustness to several geometric image deformations, such as change of viewpoint, rotation and zoom. Our results showed the reliability of the SMP approach with respect to these deformations. In addition, although our method is new, its performances tested on two databases are very close to those of several established retrieval methods: a reference retrieval method based on (local) SIFT descriptors and two versions of a fuzzy, segmentation-based UFM approach: UFM and CLUE . This indicates that the SMP approach adapts to quite different retrieval tasks, from the object level (on the Nister database) to the level of general categories (on the COREL database). Finally, our Sparse Multiscale Patches approach follows the same multiscale philosophy as the new compression standard SVC [START_REF]Scalable video coding -joint draft 6[END_REF], presuming nearly straightforward use of lowlevel bitstream components for retrieval purposes. This framework can also be used for other image processing tasks such as tracking or denoising.

Fig. 1 .

 1 Fig. 1. Building a patch of multiscale coefficients, for a single color channel image.

Fig. 2 .

 2 Fig. 2. White indicates the location of patches of largest energy (1/8 of the patches is selected for each subband).

Fig. 4 .

 4 Fig. 4. Evaluation of the robustness to JPEG compression for one query image. Displayed distances are from the query to the 6 relevant images -3 compressed (A, B, C) and 3 geometrically transformed versions of the query -and to the first 2 non-relevant images. PSNR of the compressed versions: A: 31.8dB, B: 29.7dB and C: 29.3dB.

Fig. 5 .

 5 Fig. 5. Retrieval performance of the SMP method for different subset sizes of the Nister database; the roc curves were obtained for cut-off values ranging from 1 to 9.

Fig. 6 .

 6 Fig. 6. Comparison of the retrieval performances of the SMP approach and the SIFTbased algorithm; the roc curves were obtained for cut-off values going from 1 to 9.

  .4 Dist = 12.6 Dist = 12.8 Dist = 12.9 Buildings #87 Dist = 14.4 Dist = 15.7 Dist = 17.7 Dist = 20.0 Food #56 Dist = 16.3 Dist = 18.1 Dist = 18.5 Dist = 19.1

Fig. 7 .

 7 Fig. 7. Retrieval results for 5 images of the Corel database. For each row, left to right: query image; first 4 ranked images of the database (excluding the query image).For each retrieved image, the SMP similarity measure to the query is also shown.

Fig. 8 .

 8 Fig. 8. Average Precision for each category of the Corel database. Dark blue bars: SMP approach; gray bars: UFM approach; white bars: CLUE approach.

lyzing retrieval results on the Nister dataset; namely, each of the 1,000 images was used as a query and the similarity measure to all other images was computed. The same experiment was conducted by using a state-of-the-art retrieval method based on (local) SIFT descriptors [START_REF] Lowe | Sift keypoint detector[END_REF]. For this method, the similarity measure is defined as the number of points of interest that can be matched between two images. The results of both methods were quantitatively compared by means of ROC curves. These are recall versus 1 -precision curves 1 averaged over all queries. The larger the precision and recall values, the better the retrieval performances (this corresponds to the top left side of the plot of an ROC curve).

The results of our SMP retrieval method are shown in Fig. 5 for different subset sizes of the database. Namely, average results on the first 100, 200 or 500 images are compared to those on the whole dataset (1000 images). Although the probability of retrieval errors increases with the size of the database, global


Sandrine Anthoine
email: anthoine@i3s.unice.fr

Eric Debreuve
email: debreuve@i3s.unice.fr

Paolo Piro
email: piro@i3s.unice.fr

Michel Barlaud
email: barlaud@i3s.unice.fr

Using neighborhood distributions of wavelet coefficients for on-the-fly, multiscale-based image retrieval

In this paper, we define a similarity measure to compare images in the context of (indexing and) retrieval. We use the Kullback-Leibler (KL) divergence to compare sparse multiscale image descriptions in a wavelet domain. The KL divergence between wavelet coefficient distributions has already been used as a similarity measure between images. The novelty here is twofold. Firstly, we consider the dependencies between the coefficients by means of distributions of mixed intra/interscale neighborhoods. Secondly, to cope with the high-dimensionality of the resulting description space, we estimate the KL divergences in the kth nearest neighbor framework, instead of using classical fixed size kernel methods. Query-by-example experiments are presented.

Introduction

Comparing two images in the context of (indexing and) retrieval often relies on global descriptions such as dominant colors or color distribution, or on extracted information such as salient points/regions together with local features or segmentation along with region arrangement [START_REF] Mezaris | Objectbased mpeg-2 video indexing and retrieval in a collaborative environment[END_REF][START_REF] Zhang | Optimizing metrics combining low-level visual descriptors for image annotation and retrieval[END_REF]. The philosophy here is to use a synthetic, multiscale image description based on the sparse representation in a wavelet domain. Such experiments have been conducted using the marginal distributions of the wavelet coefficients at different scales associated with the Kullback-Leibler (KL) divergence as a similarity measure between distributions [START_REF] Do | Wavelet based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF]. Nevertheless, independence between the coefficients was assumed, preventing from taking into account local image structures such as texture. In contrast, we propose to consider dependency by means of distributions of mixed intra/interscale neighborhoods of coefficients. However, this approach implies to deal with a high-dimensional statistical description space. The number of samples being too small to reasonably fill this space, fixed size kernel options to estimate distributions or divergences fail. Alternatively, we propose to estimate the KL divergence in the k-th nearest neighbor (kNN) framework [START_REF] Boltz | High-dimensional kullback-leibler distance for region-of-interest tracking: Application to combining a soft geometric constraint with radiometry[END_REF], i.e., adapting to the local sample density and directly from the samples.

Similarity between images

A central question in content-based image indexing is to define a similarity measure between images that matchesor at least is close enough to -our perception of the similarity of images. Once this is done, the images in the database can be simply ranked in increasing order of their similarity to the reference (or example) image for a query-by-example task. Perceptual studies to understand how human perceive the similarity between images are still a topic of ongoing research. Therefore, content-based image indexing system relying on such studies may be subjective and hard to implement. Here, we focus on developing an objective and mathematically defined measure that will be easily implementable.

Neighborhoods of wavelet coefficients

Let us denote by w(I) j,k the wavelet coefficient of image I at scale j and location k. I.e. it is the scalar product w(I) j,k = ψ j,k , I of I with ψ j,k , the mother wavelet ψ translated at location k and dilated at scale j.

The wavelet transform enjoys several properties that have made it quite successful in signal processing and that are relevant for the definition of similarity between images. Indeed, it provides a sparse representation of images, meaning that it concentrates the informational content of an image into few coefficients of large amplitude while the rest of the coefficients are small. This combined with a fast transform is what makes wavelet thresholding methods so powerful: in fact just identifying large coefficients is sufficient to extract where the information lies in the image. Thus it seems natural to define the feature space in the wavelet domain.

Initial thresholding wavelet methods treated each coefficient separately relying on the decorrelation of these coefficients. However, they are not independent and these dependencies are the signature of structures present in the image. For example, a discontinuity between smooth regions at point k 0 will give large coefficients at this point at all scales j (w(I) j,k0 large for all j). The most significant dependencies are seen between a wavelet coefficient w(I) j,k and its closest neighbors in scale (w(I) j-1,k ) or space (w(I) j,k±(0,1) , w(I) j,k±(1,0) ). Several models using these dependencies have been proposed and used in image enhancement [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF][START_REF] Romberg | Bayesian tree-structured image modeling using wavelet-domain hidden markov models[END_REF]. Here we use the concept of wavelet neighborhoods introduced in [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF]; these are vectors of wavelets coefficients of the form:

w(I) j,k = w(I) j,k , w(I) j-1,k , w(I) j,k±(1,0) , w(I) j,k±(0,1) . (1)
It was shown that the probability density function (pdf) of such neighborhoods allow to characterize and estimate fine spatial structures in images [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF][START_REF] Pierpaoli | Reconstructing sunyaev-zeldovich clusters in future cmb experiments[END_REF]. Hence we will define our feature space on the set of neighborhoods of wavelet coefficients of the form of (1). Critically sampled tensor wavelet transforms lack of translation and rotation invariance and so would the neighborhoods made of such coefficients. Since it is desirable to find rotated and translated versions of an image to be similar to the original one, we prefer to use a slightly more redundant transform, namely the steerable pyramid [START_REF] Portilla | Image denoising using a scale mixture of Gaussians in the wavelet domain[END_REF]. This is the decomposition on the set of dilated, translated and (Fourier-)rotated version of a mother wavelet. The image I is then represented by a set of wavelet coefficients of the form: {w(I) j,o,k } j∈Z,k∈Z 2 ,o=1..No , o indexing the orientations. Subsampling is not performed on the first level of decomposition (but is done subsequently) thus allowing that all subbands are aliasing-free (translation invariance in each scale). Moreover the different orientations allow to get some rotation invariance. Although this transform is redundant (with a factor 4N o /3), it is fast and enforces sparsity of image decomposition as do the critically sampled wavelet transforms, but it also enjoys more invariance properties than the latter.

The sampling of orientations is rather coarse (usually N o = 4). Therefore dependencies between coefficients at different orientations are less significant than across scale or space. Thus we confine the neighborhoods to each orientation, i.e the neighborhood of w(I) j,o,k is:

w(I) j,o,k = w(I) j,o,k , w(I) j-1,k,o , w(I) j,o,k±(1,0) , w(I) j,o,k±(0,1) . (2)
Hence our feature space is the set of the neighborhoods as in ( 2) for all scales j, orientations o and locations k. Let us now turn to the measure of similarity on this space.

Similarity measure between images

Since geometrically modified or slightly degraded versions of the same image as well as images containing similar objects should be close, one cannot define a measure comparing directly the neighborhoods one by one, but rather their probability distribution. More specifically, we consider the pdf of the neighborhoods of (2) for each scale and orientation, i.e. we consider the pdf p wj,o(I) of the set neighborhoods {w(I) j,o,k } k for each fixed j and o.

The considered pdf are those of coefficients that carry the informational content of the signal. The natural way to compare such pdf is to use measures derived from information theory. Here we use the KL divergence between pdfs, an approach that has also been successfully taken for other applications [START_REF] Boltz | High-dimensional kullback-leibler distance for region-of-interest tracking: Application to combining a soft geometric constraint with radiometry[END_REF]. This was also done in [START_REF] Do | Wavelet based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF][START_REF] Wang | Quality-aware images[END_REF] in the context of evaluating the similarity between images using the marginal pdf of the wavelet coefficients. We propose to use this measure on the multidimensional pdf of the neighborhoods of coefficients: the similarity between images I 1 and I 2 is a weighted sum over orientations and scales of the KL divergences between the pdf p wj,o(I1) and p wj,o(I2) :

S(I 1 , I 2 ) = j,o α j D kl (p wj,o(I1) ||p wj,o(I2) ) (3)
with p wj,o(Ii) the pdf of the wavelet neighborhoods of image I i at scale j and orientation o and α j > 0 are weights (chosen according to the redundancy of the transform). Previous works on neighborhoods of wavelet coefficients or indexation using marginal pdf of these coefficients all assumed a parametric model for the pdf involved. In the marginal case, efficient models (e.g. generalized Gaussians [START_REF] Do | Wavelet based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF][START_REF] Wang | Quality-aware images[END_REF]) lead to an analytic expression of the KL divergence as a function of the model parameter; but they are not easily generalizable to the multidimensional correlated case of wavelet neighborhoods. On the other hand, efficients multidimensional models accounting for correlations (e.g. Gaussian mixtures [START_REF] Pierpaoli | Reconstructing sunyaev-zeldovich clusters in future cmb experiments[END_REF]) fit a wide variety of multidimensional pdf but impose to estimate the KL divergence after estimating the model parameters. Besides the heavy computational cost of the consecutive estimations, the numerical stability of such cascading estimates is difficult to obtain. We prefer to make no hypothesis on the pdf at hand, hence sparing the cost of fitting the model parameters but needing to estimate the KL divergences in this non-parametric case.

Estimation of the Kullback-Leibler divergences

Let us first remind the reader that the KL divergence between two continuous pdf p 1 and p 2 is:

D kl (p 1 ||p 2 ) = p 1 (x) log p 1 (x) p 2 (x) dx = H x (p 1 , p 2 ) -H(p 1 ) (4)
where H is the differential entropy and H x is the cross entropy.

The estimation of statistical measures in the multidimensional case is hard. In particular, kernel-based methods such as Parzen estimates become unadapted due to the sparsity of samples in high dimension (curse of dimensionality): the tradeoff between a kernel with a large bandwidth to perform well in low local sample density (which oversmoothes the estimator) and a kernel with a smaller bandwidth to preserve local statistical variabilities (which results in an unstable estimator) cannot always be achieved. We use instead the kth nearest neighbor (kNN) framework [START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF] to compute the KL divergence. Indeed it follows the dual approach to the above fixed size kernel: the bandwidth adapts to the local sample density by letting the kernel contain exactly k neighbors of a given sample. Moreover it allows direct estimation of the divergence without explicitly estimating the pdf.

Assume that ǫ is a set of N ǫ samples w 1 , w 2 ,..,w Nǫ of pdf p ǫ . Fix a non-zero integer k. Denote by v d the volume of the unit sphere in R d , and ψ the digamma function. Denote by µ ǫ (g) the mean of g over ǫ:

µ ǫ (g) = 1 Nǫ Nǫ n=1 g(w n ). (5)
ρ k,ǫ (s) is the distance for s ∈ R d to its kth nearest neighbor in ǫ -{s}. kNN balloon estimates are based on the principle that p ǫ (s) is inversely proportional to the volume of the sphere containing the k nearest neighbors of s in ǫ [START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization[END_REF]:

p ǫ (s) ∼ k v d ρ d k,ǫ (s) (6)
An unbiased estimator of the Ahmad-Lin approximation of entropy [START_REF] Ahmad | A nonparametric estimation of the entropy absolutely continuous distributions[END_REF]

H al (p ǫ ) = - 1 N ǫ Nǫ n=1 log p ǫ (w n ) = -µ ǫ log p ǫ (7)
in the kNN framework was proposed in [START_REF] Goria | A new class of random vector entropy estimators and its applications in testing statistical hypotheses[END_REF] by replacing log k by ψ(k):

H(p ǫ ) = log[(N ǫ -1)v d ] -ψ(k) + d µ ǫ (log[ρ k,ǫ ]) (8)
The cross entropy estimate is then [START_REF] Boltz | High-dimensional kullback-leibler distance for region-of-interest tracking: Application to combining a soft geometric constraint with radiometry[END_REF]:

H x (p ǫ1 , p ǫ2 ) = log[N ǫ2 v d ]-ψ(k)+dµ ǫ2 (log[ρ k,ǫ1 ]). (9)
And the KL divergence estimate is:

D kl (p ǫ1 ||p ǫ2 ) = log Nǫ 2 Nǫ 1 -1 + d µ ǫ2 (log[ρ k,ǫ1 ]) -d µ ǫ1 (log[ρ k,ǫ1 ]) (10)
This expression is valid in any dimension and it is robust to the choice of k.

Numerical experiments 3.1. Setting

The database used in our numerical experiments contains twenty five 128x128 color images from the VisTex database (available at [START_REF]Vision and modeling group[END_REF]). Given the small size of the images, only two levels of the decomposition with the steerable pyramid were computed. The number of orientations is fixed to four and the number of neighbors in the kNN procedure to ten.

So far, we have described the feature space and similarity measure considering implicitly single channel images (like gray level images). To extend them to the multichannel case, we consider the luminance/chrominance space (Y, Cb, Cr). Since the luminance and chrominance channels are fairly well decorrelated, one can in first approximation consider them independent. Hence, we simply sum the KL divergences obtained for each channel separately.

Retrieval results

The results for 5 of the images in the database are displayed in Fig. 1. In this figure, each row displays the retrieval result for the example (or reference) image shown on the leftmost column. From the second column on, one can see the first three images in the database ranked by our similarity distance (the leftmost, the most similar), excluding the example image (which is always at a distance of zero).

In general, our method seems to perform very well. In particular, images coming from the same scene (see rows 2 to 4 in Fig. 1) are ranked first. In this database such images are usually translated versions of one another. Hence this experiment shows that our method is robust to translation. Similar textured images such as trees, grass, and grids are also correctly classified (see the last row of Fig. 1).

Complexity and computation time

For one query image, the computational cost of the retrieval procedure is the number of image in the database times the cost of computing a similarity between two images. Denoting by N the number of pixel in an image, the similarity computation is made in three steps of complexity:

• O( √ N ) for the steerable transform, • O(N ) for the design of the wavelet neighborhoods,

• O(N log N ) for the evaluation of kNN distances (via a classical KD-tree implementation). Accordingly, most of the computational effort is put in the evaluation of the kNN distances. To improve this, we reduce our feature space by selecting a small proportion of the neighborhoods to evaluate the KL divergences. We select those with the largest central coefficient, thus exploiting the sparsity of wavelet representation. Fig. 2. shows how the computational time evolves then in O(M log M ) where M is the number of selected coefficients (green curve with circles) while the similarity measure remains consistent (the similarity between image 9 and its 3 closest matches are displayed). Selecting only 1/32 of the coefficients leaves us with results of the same accuracy as with all coefficients while greatly reducing the computation time.

Conclusion

In this paper, we proposed a similarity measure between images based on the KL divergence between multidimensional pdf of wavelet coefficients grouped in coherent sets called neighborhoods. The KL divergence is estimated nonparametrically via a kNN approach.

Experiments on small images show good performances of the proposed measure in the retrieval problem, particularly its robustness to simple geometric tranforms and to the sparsity of the feature space. Future works will focus on dataset with larger images as well as evaluation of the performances through recall-precision curves.

Figure 1 .

 1 Figure 1. Retrieval results. Left to right: reference image; first 3 ranked images.

Figure 2 .

 2 Figure 2. Evolution of similarity and computing time with proportion of coefficients used.