G Quellec

M Lamard

G Cazuguel
email: guy.cazuguel@telecom-bretagne.eu

B Cochener

C Roux

Wavelet Optimization for Content-Based Image Retrieval in Medical Databases

Keywords: Content-Based Image Retrieval, custom wavelet, lifting scheme, diabetic retinopathy, mammography

may come

Wavelet optimization for content-based image retrieval in medical databases

Gwénolé Quellec, Mathieu Lamard, Guy Cazuguel, Béatrice Cochener,

INTRODUCTION

Images have always been used in medicine for teaching, diagnosis, and management purposes. Now medical imaging systems produce more and more digitized images in all medical fields: visible, ultrasound, X-ray tomography, MRI, nuclear imaging, etc. These images are very interesting for diagnostic purposes: they are directly related to the patient pathology and medical history, and may be used for diagnosis aid. However the amount of images we can access nowadays is so huge that database systems require efficient indexing to enable fast access to images in databases. Automatic image indexing using image digital content (Content-Based Image Retrieval -CBIR) is one of the possible and promising solutions to effectively manage image databases [START_REF] Nastar | Indexation d'images par le contenu: un etat de l'art[END_REF]. In recent years, CBIR has been applied to several applications such as industrial [START_REF] Tao | Which components are important for interactive image searching[END_REF], forensic [START_REF] Algarni | A novel technique for automatic shoeprint image retrieval[END_REF], remote sensing [START_REF] Hafiane | Region-based CBIR in GIS with local space filling curves to spatial representation[END_REF][START_REF] Xie | A rotation-and flip-invariant algorithm for representing spatial continuity information of geographic images in contentbased image retrieval[END_REF] and medical applications [START_REF] Müller | A review of content-based image retrieval systems in medical applications -clinical benefits and future directions[END_REF]. In a medical context, the potential applications of CBIR are categorization (i.e. retrieving images of the same anatomical regions) [START_REF] Rahman | A framework for medical image retrieval using machine learning and statistical similarity matching techniques with relevance feedback[END_REF][START_REF] Horsthemke | Task-oriented medical image retrieval[END_REF][START_REF] Greenspan | Medical image categorization and retrieval for PACS using the GMM-KL framework[END_REF], data mining [START_REF] Siadat | Content-based image database system for epilepsy[END_REF], retrieving visually similar lesions [START_REF] Kim | A new way for multidimensional medical data management: Volume of interest (VOI)-based retrieval of medical images with visual and functional features[END_REF][START_REF] El-Naqa | A similarity learning approach to content-based image retrieval: application to digital mammography[END_REF][START_REF] Doyle | Using manifold learning for content-based image retrieval of prostate histopathology[END_REF][START_REF] Balmashnova | Content-based image retrieval by means of scale-space top-points and differential invariants[END_REF], typically for educational purposes, and finally diagnosis aid [START_REF] Dy | Unsupervised feature selection applied to content-based retrieval of lung images[END_REF][START_REF] Pourghassem | Content-based medical image classification using a new hierarchical merging scheme[END_REF][START_REF] Antani | Evaluation of shape similarity measurement methods for spine X-ray images[END_REF], which is the application we will focus on in this article. The purpose is to let an expert carry out a search in the database, without formulating a semantic description of the image he is examining: the expert simply places the image as a query to a system, which retrieves the most similar images in the database. Finally, the expert chooses the most similar images, from his point of view, in the set of retrieved images, and accesses the associated information. To compare the query image to another image in the database, feature vectors characterizing their numerical content, called signatures, are extracted. Then, comparing these two images simply comes down to comparing their signatures. In recent CBIR systems, images are mostly characterized visually by their shape [START_REF] Pourghassem | Content-based medical image classification using a new hierarchical merging scheme[END_REF][START_REF] Dai | Unbalanced region matching based on two-level description for image retrieval[END_REF][START_REF] Zhang | Study and evaluation of different fourier methods for image retrieval[END_REF][START_REF] Antani | Evaluation of shape similarity measurement methods for spine X-ray images[END_REF], their edges [START_REF] Rahmani | Localized content-based image retrieval[END_REF][START_REF] Rahman | A framework for medical image retrieval using machine learning and statistical similarity matching techniques with relevance feedback[END_REF][START_REF] Banerjee | Edge based features for content based image retrieval[END_REF][START_REF] Han | A shape-based image retrieval method using salient edges[END_REF], textural features [START_REF] Rahman | A framework for medical image retrieval using machine learning and statistical similarity matching techniques with relevance feedback[END_REF][START_REF] Kokare | Rotation-invariant texture image retrieval using rotated complex wavelet filters[END_REF][START_REF] Chun | Image retrieval using BDIP and BVLC moments[END_REF][START_REF] Rallabandi | Rotation-invariant texture retrieval using wavelet-based hidden markov trees[END_REF][START_REF] Pourghassem | Content-based medical image classification using a new hierarchical merging scheme[END_REF][START_REF] Oliveira | Towards applying content-based image retrieval in the clinical routine[END_REF][START_REF] Muneeswaran | Texture image segmentation using combined features from spatial and spectral distribution[END_REF][START_REF] Yu | A three-object model for the similarity searches of chest CT images[END_REF] or color histograms [START_REF] Rahman | A framework for medical image retrieval using machine learning and statistical similarity matching techniques with relevance feedback[END_REF][START_REF] Yap | Content-based image retrieval using legendre chromaticity distribution moments[END_REF][START_REF] Khanh | Image retrieval based on regions of interest[END_REF][START_REF] Konstantinidis | Image retrieval based on fuzzy color histogram processing[END_REF]. Several features proposed in the literature are invariant to translation [START_REF] Bishnu | Euler vector for search and retrieval of gray-tone images[END_REF], rotation [START_REF] Kokare | Rotation-invariant texture image retrieval using rotated complex wavelet filters[END_REF][START_REF] Bishnu | Euler vector for search and retrieval of gray-tone images[END_REF][START_REF] Rallabandi | Rotation-invariant texture retrieval using wavelet-based hidden markov trees[END_REF][START_REF] Algarni | A novel technique for automatic shoeprint image retrieval[END_REF][START_REF] Sastry | A wavelet based multiresolution algorithm for rotation invariant feature extraction[END_REF][START_REF] Xie | A rotation-and flip-invariant algorithm for representing spatial continuity information of geographic images in contentbased image retrieval[END_REF], size [START_REF] Bishnu | Euler vector for search and retrieval of gray-tone images[END_REF][START_REF] Algarni | A novel technique for automatic shoeprint image retrieval[END_REF][START_REF] Lo | An implementation of the scale invariant feature transform in the 2.5d domain[END_REF] or illumination [START_REF] Greenspan | Medical image categorization and retrieval for PACS using the GMM-KL framework[END_REF]. In medical applications, rotation and size invariance are mostly necessary if we want to retrieve images similar in shape: typically, if we want to retrieve anatomical regions [START_REF] Rahman | A framework for medical image retrieval using machine learning and statistical similarity matching techniques with relevance feedback[END_REF][START_REF] Horsthemke | Task-oriented medical image retrieval[END_REF] or if we want to characterize lesions or patterns within a user-defined ROI [START_REF] Kim | A new way for multidimensional medical data management: Volume of interest (VOI)-based retrieval of medical images with visual and functional features[END_REF] or within an automatically segmented region [START_REF] Yu | A three-object model for the similarity searches of chest CT images[END_REF][START_REF] Siadat | Content-based image database system for epilepsy[END_REF][START_REF] Antani | Evaluation of shape similarity measurement methods for spine X-ray images[END_REF]. In this article, we propose a diagnosis aid system intended to retrieve images with the same pathology and at the same disease severity level, not simply representing the same anatomical region. For that purpose, we have to characterize images by the lesions they contain. And we do not want the user to define ROIs: medical images sometimes contain hundreds of lesions, it would be too time consuming and physicians would not use such a system. Also, we do not want to rely on segmentation to characterize the lesions, since it is hard to design a robust segmentation algorithm for each type of lesions. So we propose to characterize images by their textural content, at different scales: it allows us to characterize lesions of different sizes, in a generic fashion. Our approach is invariant to translation. For many medical imaging modalities, such as eye fundus imaging or mammography, the orientation and the scale of images are approximately normalized, due to the acquisition protocol. If for a particular imaging modality they are not, then they should be normalized prior to image retrieval; this possible issue is not addressed in the article. The outstanding challenges of CBIR in medical applications are the following. First, how to bridge the gap between the extracted low-level features and the high-level perception of similarity of the physicians [START_REF] Rahman | A framework for medical image retrieval using machine learning and statistical similarity matching techniques with relevance feedback[END_REF]? We addressed this issue by automatically adapting the system parameters to a high-level criterion, using a supervised learning procedure. Then, how to integrate the proposed retrieval systems with the compression and archiving standards [START_REF] Greenspan | Medical image categorization and retrieval for PACS using the GMM-KL framework[END_REF]? We addressed this issue by defining a wavelet-based image signature that may be extracted as images are compressed by a wavelet-based compression standard, such as JPEG 2000 [START_REF] Taubman | JPEG2000: Image Compression Fundamentals, Standards and Practice[END_REF]. Then, how to combine visual features with medical metadata [START_REF] Greenspan | Medical image categorization and retrieval for PACS using the GMM-KL framework[END_REF][START_REF] Antani | A biomedical information system for combined content-based retrieval of spine x-ray images and associated text information[END_REF][START_REF] Le Bozec | Towards content-based image retrieval in a HIS-integrated PACS[END_REF][START_REF] Shao | Medical image retrieval based on visual contents and text information[END_REF]? And finally, how to efficiently retrieve images in large databases [START_REF] Oliveira | Towards applying content-based image retrieval in the clinical routine[END_REF][START_REF] Dy | Unsupervised feature selection applied to content-based retrieval of lung images[END_REF][START_REF] Iakovidis | A pattern similarity scheme for medical image retrieval[END_REF]? We addressed the last two issues in a previously published article [START_REF] Quellec | Recherche de cas médicaux multimodaux à l'aide d'arbres de décision[END_REF]: we used knowledge discovery trees to retrieve complete patient files consisting of possibly several images with metadata, thanks to wavelet-based image signatures, which enables a fast retrieval and an increased retrieval performance. Wavelets have been widely used in image processing since the publication of the JPEG 2000 standard, including in the CBIR field [START_REF] Xiong | Wavelet-based texture features can be extracted efficiently from compressed-domain for JPEG2000 coded images[END_REF][START_REF] Do | Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF]. Do and Vetterli, in particular, proposed a texture retrieval method based on the global distribution of wavelet coefficients in a texture image. We explore in this article two approaches to define image signatures and their associated distance measures from the global distribution of wavelet coefficients in images. Both approaches are based on a specific model for wavelet coefficient distributions: the first one is based on histograms and, like Do's approach, the second one is based on generalized Gaussian functions. In fact, this article generalizes Do's image representation, designed for texture images, to make it general, in particular to make it applicable to real-world medical images, and automatically tunable, so that it can be adapted to the studied images and pathologies. To adapt this image representation, we propose to use a custom wavelet basis and a custom decomposition scheme. The two approaches are evaluated on two medical image databases: one for diabetic retinopathy follow up and one for screening mammography. To show that the method is robust to illumination and slight orientation variations, it was also evaluated on a face database. The setup of the article is as follows. The signatures and distance measures are described in section 2. Section 3 introduces our approach to search for an optimal wavelet basis. We discuss the use of an optimal wavelet transform for both image indexing and compression in section 4. Sections 5 and 6 explain how our method is evaluated and calibrated. And the method is applied to the three databases in section 7. We end with a discussion and conclusion in section 8.

Proposed image signatures and distance measures

We propose to build image signatures from their wavelet transform. The wavelet transform is indeed convenient since it is highly tunable: in particular we can use a custom basis of functions, called wavelets, to decompose images. This wavelet basis will be tuned in order to improve the retrieval performance of the system. To implement the wavelet transform, we use the lifting scheme [START_REF] Sweldens | The lifting scheme: a custom-design design construction of biorthogonal wavelets[END_REF], which is particularly suited to design a custom wavelet basis (see section 2.2). It is also computationally efficient, hence its use in compression, for instance. Two different image signatures are then extracted from these decompositions in section 2.4. The associated distance measures are described in section 2.5. But first, let us set out the basics of multiresolution analysis.

Multiresolution analysis

A multiresolution approximation is a nested sequence of linear spaces (V j ) j∈Z for approximating functions f ∈ L 2 (R) such that:

• V j ⊂ V j+1 , ∀j ∈ Z • +∞ j=-∞ V j is dense in L 2 (R) and +∞ j=-∞ V j = {0} • f (x) ∈ V j ⇔ f (2x) ∈ V j+1 , ∀j ∈ Z • f (x) ∈ V j ⇒ f (x -2 -k ) ∈ V j , ∀j, k ∈ Z
The approximation of a function f ∈ L 2 (R) at a resolution 2 j is defined as the projection of f on V j . There is a function Φ ∈ L 2 (R), called scaling function, such that for any j ∈ Z,

(x → Φ jk (x) = √ 2 j Φ(2 j x -k)) k∈Z is a basis of V j . The scaling function Φ satisfies the refinement relation of equation 1. Φ(x) = 2 ∞ k=-∞ h k Φ(2x -k) (1)
The additional information available in the approximation of f at the resolution 2 j+1 as compared with the resolution 2 j , is given by the projection of f on the orthogonal complement of V j in V j+1 , noted W j . There is a function Ψ ∈ L 2 (R), called wavelet function, such that for any j ∈ Z,

(x → Ψ jk (x) = √ 2 j Ψ(2 j x -k))
k∈Z is a basis of W j . The wavelet function Ψ satisfies the refinement relation of equation 2.

Ψ(x) = 2 ∞ k=-∞ g k Φ(2x -k) (2)
h k and g k are the coefficients defining the family of the wavelet functions.

The wavelet and the scaling filters are noted G = (g k ) k∈Z and H = (h k ) k∈Z , respectively. {V j } j∈Z and {W j } j∈Z constitute what we call a multiresolution analysis framework.

Biorthogonal wavelets and the lifting scheme

To reconstruct a signal decomposed by the procedure above, we define a dual multiresolution analysis framework ({ Ṽj } j∈Z ,{ Wj } j∈Z ) associated with a dual scaling function Φ and a dual wavelet function Ψ, that are biorthogonal to Φ and Ψ (see equation 3).

                   < Φjk , Φ jk >= δ kk < Ψjk , Ψ jk >= δ kk < Ψjk , Φ jk >= 0 < Φjk , Ψ jk >= 0 (3)
Φ and Ψ satisfy refinement relations such as equations 1 and 2, with coefficients H = ( hk ) k∈Z and G = (g k ) k∈Z respectively. A function f ∈ L 2 (R) can be written as follows:

x → f (x) = j,k∈Z < f, Ψj,k > Ψ j,k (x) = k∈Z < f, ΦJ,k > Φ J,k (x) + j≤J,k∈Z < f, Ψj,k > Ψ j,k (x), J ∈ Z (4)
The lifting scheme allows the construction of compactly supported biorthogonal wavelets with compactly supported duals, i.e. such that the refinement filters H, G, H and G are finite filters (a finite number of coefficients are non-zero). It relies on a simple relationship between all multiresolution analysis frameworks that share the same scaling function: let Φ 0 , Φ0 , Ψ 0 , Ψ0 be an initial set of biorthogonal scaling and wavelet functions, then a new set Φ, Φ, Ψ, Ψ can be found according to equation 5 [START_REF] Sweldens | The lifting scheme: a custom-design design construction of biorthogonal wavelets[END_REF].

                   Φ(x) = Φ 0 (x) Ψ(x) = Ψ 0 (x) -k s k Φ 0 (x -k) Φ(x) = 2 k h0 k Φ(2x -k) + k s -k Ψ(x -k) Ψ(x) = 2 k gk Φ(2x -k) (5)
Coefficients s k of equation 5 can be freely chosen. Hence we can start with an initial set of biorthogonal filters and generate new biorthogonal filters by the previous relation. As a consequence, the lifting scheme is particularly suited to design a custom wavelet basis: we can generate several biorthogonal filters, using an optimization procedure, in order to find the wavelet basis maximizing a criterion (the retrieval performance of the system, in our application). [START_REF] Claypoole | Adaptive wavelet transforms via lifting[END_REF] proposed a filter bank adapted from equation 5. It uses the lazy wavelet [START_REF] Sweldens | The lifting scheme: a custom-design design construction of biorthogonal wavelets[END_REF] as initial biorthogonal filter and two filters P = (p i ) i∈N and U = (u i ) i∈N are defined from the s k coefficients. One step in the lazy wavelet transform consists in splitting the signal into its odd and even indexed samples. One step of the wavelet decomposition is given as follows:

• split: the signal x[n] is split into its odd x o [n] and even x e [n] coefficients. The filter bank is depicted in figure 1. Let N p and N u be the length of the linear filters P and U . [START_REF] Claypoole | Adaptive wavelet transforms via lifting[END_REF] demonstrated that, to form a biorthogonal wavelet, filters P and U only have to satisfy the following conditions:

             Np i=1 p i = 1 Nu i=1 u i = 1 2 (6)
There are then N p +N u -2 undetermined coefficients. Wavelets designed by the lifting scheme have a support length equal to s/t given in equation 7, where s is the support of the low-pass filter and t is the support of the high-pass filter.

     s = 2(N p + N u ) -3 t = 2N p -1 (7)
For instance, the following wavelets have been evaluated in the proposed system:

• the Le Gall 5/3 biorthogonal wavelet, which is used in the JPEG 2000 standard part I (support=5/3) (le Gall and Tabatabai, 1988) • the Daubechies 9/7 biorthogonal wavelet, also in the JPEG 2000 standard (support=9/7) • the orthogonal Haar wavelet (support=2/2) • Cubic B-spline wavelet (support=4/2) [START_REF] Cohen | Bi-orthogonal bases of compactly supported wavelets[END_REF] • the Daubechies 4-tap orthogonal wavelet (support=4/4) [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF] • the Daubechies 6-tap orthogonal wavelet (support=6/6) [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF] • Adapted wavelets (see section 3)

Decomposition of images

In the case of images, a separable wavelet decomposition is commonly used; it means that the one-dimensional filtering operations described above are performed on both rows and columns separately (the same filters are used along both directions). Three subimages of coefficients are obtained at each level l: lHL, lLH and lHH, and a unique approximation subimage at the last level N l : N l LL. These coefficients are obtained from the decomposition of (l -1)LL (0LL is the raw image):

• lHL contains high frequency coefficients along columns and low frequency coefficients along lines, • lLH contains low frequency coefficients along columns and high frequency coefficients along lines, • lHH contains high frequency coefficients along columns and along lines, • lLL contains low frequency coefficients along columns and along lines, it is decomposed if l < N l .

Subimages of coefficients are called subbands. An image decomposed on N l levels thus consists of 3×N l +1 subbands. An example of image and its wavelet transform is shown in figure 2. Typically, the interest of using a different decomposition scheme is to find a basis relative to which the given image or collection of images has the lowest information cost, for compression purposes. Also, some subbands may emerge in which it is easier to discriminate between images of different classes, for classification or image retrieval purposes.

Image signatures

In order to characterize the image texture at different scales, we characterize the distribution of the wavelet coefficients in each subband of such a decomposition.

Histogram Representation

The first approach we consider to define an image signature is to build a histogram of the transformed image coefficients in each subband. If images are decomposed with a pyramidal scheme on N l levels, they consist of 3 × N l + 1 subbands: there are 3 subbands of details at each scale l ≤ N l (lHH, lHL and lLH) plus an approximation (N l LL), 3 × N l + 1 histograms are thus built. The signature is a vector formed by the concatenation of these histograms. The size of the resulting signature is relatively large. Indeed, let N B be the number of bins in each histogram, if N l =4 and N B = 32, the signature consists of N B ×(3×N l +1) = 416 coefficients. For this reason, we propose an alternative method, in the next paragraph, that alleviates the problem.

Generalized Gaussian Representation

To characterize textures, [START_REF] Wouwer | Statistical texture characterization from discrete wavelet representations[END_REF] proposed to model the distribution of transformed coefficients in each subband with zero-mean generalized Gaussian functions. A generalized Gaussian function may be expressed as in equation 8:

         p(x; α, β) = β 2αΓ 1 β e -( |x| α ) β Γ(z) = ∞ 0 e -t t z-1 dt, z > 0 (8)
We propose to discover to what extent this representation holds for any image dataset. The goodness of fit of generalized Gaussians to the actual distribution of wavelet coefficients is assessed by Kolmogoroff-Smirnov (K-S) tests (von Mises, 1964). For a reference wavelet basis (the Daubechies 9/7 wavelet basis), a K-S test is performed for each subband of the wavelet decomposition of each image in the dataset. The K-S tests measure the maximal absolute difference D between the cumulative distribution of the data, i.e. the wavelet coefficients in a given subband of a given image, and the cumulative distribution function (CDF) of the model, i.e. the generalized Gaussian function that best fits the distribution of the data. The CDF of a zero-mean generalized Gaussian is given below:

             cdf (x; α, β) = 1 2 + sign(x) γ 1 β , |x| α β 2Γ 1 β γ(s, z) = z 0 e -t t s-1 dt, z > 0 (9)
We propose to build signatures consisting of the maximum likelihood estimator (α, β) of a generalized Gaussian distribution in each subband where the generalized Gaussian representation holds, and an histogram of the transformed image coefficients in the others. Intuitively, the last subband has no reason to have a specific distribution, since it is mostly an approximation of the image.

For the scenario presented in the last paragraph, assuming the representation holds for each high frequency subband but not for the last subband, the signature consists of 2 × (3 × N l ) + N B = 56 coefficients. To compute the maximum likelihood estimator of a distribution, we applied the method proposed by [START_REF] Do | Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF], extending the work of [START_REF] Varanasi | Parametric generalized Gaussian density estimation[END_REF]. This method relies on the Newton-Raphson root finding algorithm [START_REF] Kay | Fundamentals of statistical signal processing: estimation theory[END_REF], that sometimes diverges on our images. This phenomenon is more or less important depending on the wavelet basis used to decompose images. As a consequence, we applied a robust root finding method proposed in (Press et al., 1992c): a hybrid method combining Newton-Raphson and the bisection algorithm.

Distance Measures

Each of the signatures presented above are associated with a particular distance measure, as described in the following paragraphs. The proposed distance measures are designed in order to profit from the characteristic of their associated signature.

Distance Between Histogram Based Signatures

The distance used to compare two images Im 1 and Im 2 by the first approach is base on the l1-distance between histograms. [START_REF] Varela | Indexation et recherche d'images par le contenu, utilisant des informations de compression d'images : application aux images médicales[END_REF] compared other distance measures between histograms, on several medical databases, and this distance measure turned out to give the highest precision. The distance measure is given in equations 10 and 11:

d(Im 1 , Im 2 ) = 3N l +1 i=1 λ i (H 1 i -H 2 i ) (10)
H 1 i -H 2 i = N B j=1 H 1 i (j) -H 2 i (j) (11)
where H n i (j) is the value of the j th bin of the i th normalized histogram of image n and (λ i ) i=1..3N l +1 is a set of tunable weights.

Distance Between Generalized Gaussian Based Signatures

Do and Vetterli (2002) used the Kullback-Leibler divergence to compare two generalized Gaussian distributions. It measures the difference from a query probability distribution q to an arbitrary probability distribution i, defined by parameters θ q and θ i , respectively. Its general form is given in equation 12:

D(p(X; θ q )||p(X; θ i )) = p(X; θ q ) log p(X; θ q ) p(X; θ i ) dx (12)
Injecting equation 8 in 12, we obtain equation 13:

D(p(X; α 1 , β 1 )||p(X; α 2 , β 2 )) = log β 1 α 2 Γ( 1 β 2 ) β 2 α 1 Γ( 1 β 1 ) + α 1 α 2 β 2 Γ( β 2 +1 β 1 ) Γ( 1 β 1 ) -1 β 1 (13)
The proposed distance measure between two images Im 1 and Im 2 is a combination of equations 13 and 11 (it is not actually a distance measure since it is not symmetric):

d(Im 1 , Im 2 ) = 3N l i=1 λ i D(p(X; αi,1 , βi,1 )||p(X; αi,2 , βi,2 )) +(H 1 3N l +1 -H 2 3N l +1 ) (14)
where (α i,n , βi,n ) is the maximum likelihood estimator for the i th subband of image n.

The two couples (signature,distance measure) presented in this section will be compared experimentally in section 7.4. In order to understand whether the potential performance difference between the two couples come from the signature or the distance, the couple (generalized Gaussian representation, l1-distance) will also be assessed: the area between two probability distribution functions, with parameters (α i,1 , βi,1 ) and (α i,2 , βi,2 ), is estimated according to the trapezoidal rule on an interval adapted to both distributions

([-20α i,1 ; 20α i,1 ] ∪ [-20α i,2 ; 20α i,2 ]
). But first, we will see in the following section how the wavelet transform can be tuned in order to adapt these generic image representations to a given database.

Wavelet Adaptation

A major interest of wavelet decompositions is that we can choose the wavelet basis used to decompose images. Thus it can be adapted to the images and pathologies studied. However, the design of matched wavelet bases is usually a difficult task, as several constraints have to be satisfied simultaneously, in particular perfect reconstruction and vanishing moments. Moreover, these properties might not be sufficient to process specific signals or images. Several studies have thus been carried out to adapt a wavelet to a signal of reference, with different constraints on the wavelet properties and on the adaptation criteria. Usually, the goal is to minimize the L 2 norm between the signal and an approximation derived from its wavelet decomposition [START_REF] Tewfik | On the optimal choice of a wavelet for signal representation[END_REF][START_REF] Gupta | A new method of estimating wavelet with desired features from a given signal[END_REF], for compression purposes, for instance. These methods have been applied to medical classification problems, in order to match waveform signals [START_REF] Maitrot | Design of wavelets adapted to signals and application[END_REF]de Sobral Cintra et al., 2004). However, in our case, we do not intend to match a signal of reference: we want to find the wavelet basis that maximizes the retrieval performance of our system.

Adapting the wavelet basis with a high level criterion

We propose in this article a method to find, within the lifting scheme framework, the wavelet basis that maximizes a given high level criterion, such as the retrieval performance of our system. The lifting scheme is convenient since it generates biorthogonal wavelet bases, and it is known that biorthogonality implies perfect reconstruction of decomposed images [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. Indeed, any couple of filters P /U (of size N p and N u , respectively) satisfying the simple relations of equation 6 defines a biorthogonal wavelet filter bank.

First solution: adapting biorthogonal wavelets

A first solution to adapt the wavelet filter, within the lifting scheme framework, is thus to browse the space of biorthogonal filter banks, by an optimization procedure, and to select the filter bank that maximizes the adaptation criterion. The first N p -1 (resp. N u -1) coefficients of P (resp. U ) can be freely chosen, and the last coefficient of both vectors are computed in order to satisfy equation 6; as a consequence, we have N p + N u -2 degrees of freedom. The optimization procedure we used is described in appendix A: it is based on a genetic algorithm and Powell's direction set method.

Second solution: adapting biorthogonal wavelets with vanishing moments

In an image compression application, the N p + N u -2 degrees of freedom may be tuned in order to minimize the prediction error. To achieve this goal, [START_REF] Claypoole | Adaptive wavelet transforms via lifting[END_REF] proposed a solution to eliminate the first N p moments of the primal wavelet Ψ and the first N u moments of the dual wavelet Ψ.

In our application, we wish to find a filter bank that would imply a high retrieval performance. However, it can be useful to design wavelets with several vanishing moments: it makes it possible to reduce the dimension of our optimization problem. Another interest of designing wavelets with a desired number of vanishing moments is discussed in section 4: it allows the design of wavelets suitable for both indexing and compression. So, we propose a solution to find a trade-off between these two goals: maximizing the retrieval performance and eliminating the lowest order moments of the wavelet. First, we eliminate the first M p (resp. M u ) moments of the primal (resp. dual) wavelet, where M p ≤ N p and M u ≤ N u : it sets M p and M u degrees of freedom to design filters P and U , respectively. And the remaining N p -M p and N u -M u degrees of freedom, respectively, are set in order to maximize our high level-criterion. Note that the first solution described above is a particular case of this second solution, where M p = M u = 1; the same optimization procedure (described in appendix A) is used to find the (N p -M p ) + (N u -M u ) free coefficients left after eliminating the first wavelet moments. The way we adapt the predict filter (P ) and the update filter (U ) is detailed in the next sections.

Adapting the predict filter (P)

In order to control the regularity of the primal wavelet, we will eliminate its first M p ≤ N p moments:

t n Ψ(t)dt = 0, n = 0..M p -1 (15)
For that purpose, we define the following matrix V , of size M p × N g , where N g = 2N p + 1 is the size of the wavelet filter G (see equation 7):

V (k, l) = (l + 1 -N p ) k , k = 0..M p -1, l = 0..N g -1 (16) V.G = 0 (17)
If the wavelet filter G satisfies equation 17, then the first M p moments of the primal wavelet vanish (each line in V eliminates one moment) [START_REF] Claypoole | Adaptive wavelet transforms via lifting[END_REF].

We would like to transform this system in order to determine filter P . Comparing the two equivalent filter banks depicted in figure 1 Let the central coefficient in G be indexed by 0; we notice that the even coefficients in G are equal to 0, except the central coefficient which is equal to 1. As a consequence, we can reduce the dimension of the system above. Let C be the matrix containing the even coefficients in V , equation 17 becomes:

C. p 0 . . . p Np-1 t = 1 0 . . . 0 t ( 18
)
The right-hand side of this equation comes from the product of the central column of V , i.e. (0 0 = 1, 0 1 = 0, 0 2 = 0, . . .) t , by the 1 in the central coefficient of G when solving equation 17 (the convention 0 0 = 1 was proposed by [START_REF] Claypoole | Adaptive wavelet transforms via lifting[END_REF]). Notice that the constraint associated with the first line in C is equivalent to the first biorthogonality relation (see equation 6).

So we have M p ≤ N p equations and N p unknowns. We propose to generate coefficients (p i ) i=0..Np-Mp-1 by the optimization module, in order to maximize the retrieval performance of the system (see appendix A) and to determine coefficients (p i ) i=Np-Mp..Np-1 thanks to equation 18. As a consequence, we have to solve a linear system A.x = b of dimension M p , where: We propose to solve this system by the LU method (Press et al., 1992a).

• A is a square matrix of size M p ,

Adapting the update filter (U)

In order to control the regularity of the dual wavelet, we will eliminate its first M u ≤ N u moments:

t n Ψ(t)dt = 0, n = 0..M u -1 (19)
As in the previous section, we design a matrix V of size M u × N g to eliminate the M u moments of G. Since, we work with biorthogonal wavelets, filters G and H are related by equation 20 [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] and

N g = N h = 2(N p +N u )-3 (see equation 7). gk = (-1) k h k , ∀k (20)
V is thus defined by equation 21 and G must satisfy equation 22.

V (k, l) = (l + 2 -N p -N u) k , k = 0..M u -1, l = 0..N h -1 (21) V . G = 0 (22)
We can express this system as a function of filters P and U . Indeed, c[n], the low-pass filtered version of the input signal x[n], satisfies the relations

c[n] = H(x[n]) = x e [n] + U (d[n]) = x e [n] + U (x o [n] -P (x e [n]
)): as a consequence, if for instance P = (p 0 , p 1 ) and U = (u 0 , u 1 ), then H = (-p 0 u 0 , u 0 , 1p 1 u 0 -p 0 u 1 , u 1 , -p 1 u 1 ) (see figure 4 (b)) and according to equation 20, G = (-p 0 u 0 , -u 0 , 1 -p 1 u 0 -p 0 u 1 , -u 1 , -p 1 u 1 ). Suppose that filter P has already been determined by the method proposed in section 3.3.1, then we can express equation 22 as a function of U . We propose a convenient solution below.

Let C be the matrix of size M u × N u defined in equation 23, in which gj i is the factor of u j in the expression of gi (for instance, if gi = 1 -p 1 u 0 -p 0 u 1 , then g0 i = -p 1 , g1 i = -p 0 , g2 i = 0, etc.).

C (k, l) = u V (k, u)g l u (23)
Then, equation 22 becomes:

C .u = -1 0 . . . 0 t ( 24
)
The right-hand side of this equation comes from the product of the central column of V , i.e. (0 0 = 1, 0 1 = 0, 0 2 = 0, . . .) t , by the 1 in the central coefficient of G. So we have M u ≤ N u equations and N u unknowns. We propose to generate coefficients (u i ) i=0..Nu-Mu-1 by the optimization module, in order to maximize the retrieval performance of the system (see appendix A) and to determine coefficients (u i ) i=Nu-Mu..Nu-1 thanks to equation 24. As a consequence, we have to solve a linear system A .x = b of dimension M u , where:

• A is a square matrix of size M u , which consists of the last M u columns in C . • Let B be the matrix containing the first N u -M u columns in C . The right-hand side b is equal to -1 0 . . . 0 t -B . u 0 . . . u Nu-Mu-1 t .
We also propose to solve this system by the LU method. This concludes the description of the method we propose to adapt the wavelet basis with some high level criterion. We will now see how to adapt it in order to maximize the performance of a retrieval system.

Image indexing and compression

We envisage to extract the image signatures described above as images are compressed to be stored in a database: the same decomposition would be used for both indexing and compressing images. With this in mind, the use of the lifting scheme is particularly convenient, since it is used in compression standards such as JPEG 2000. In our application, we wish to find a filter bank that would imply a high retrieval performance and that may be used to compress images efficiently before they are stored in the database. That is why the second solution we proposed to adapt the wavelet (see section 3.3) is particularly suited: we build a wavelet basis with a desired number of vanishing moments, which ensures a good compression rate [START_REF] Claypoole | Adaptive wavelet transforms via lifting[END_REF]. In other words, the custom wavelet we build can be used to store images in a database efficiently.

Retrieval Performance

To tune the proposed retrieval system for a given database, we search for the system settings that maximize the retrieval performance. To evaluate the performance of a medical retrieval system, it is not possible to ask physicians to examine every image in the dataset and decide if it is relevant for each other image in this set. So, the evaluation of such systems simply relies on a classification of the images in the database, defined a priori by the physicians. However, images belonging to different predefined classes may be considered relevant by the physicians: for instance, two images from patients with different diseases, and consequently assigned to different classes, may contain lesions of the same type. So the retrieval performance we measure is a lower bound on the actual performance. Let N be the number of images in the database and V n the relevance of an image n for a query: V n = 1 if n and the query image belong to the same class, V n = 0 otherwise. The usual way to evaluate a retrieval system performance is the following [START_REF] Müller | Performance evaluation in content-based image retrieval: Overview and proposals[END_REF]: 1) each image in the database is presented to the system as a query, 2) the system ranks the N -1 remaining images by increasing order of the distance measure, 3) for each integer k = 1..N -1, the number of relevant images selected is computed:

A k = k-1 n=0 V n .
Then, the following standard performance values are derived:

• the recall: the fraction of relevant images that are selected R k = A k N -2 n=0 V n • the precision: the fraction of selected images that are relevant

P k = A k k
To evaluate the system, precision-recall curves are usually plotted: they consist of linked points (P k ,R k ), k = 1..N -1. However, to optimize the system settings, a numerical value is needed: we measured the retrieval performance of the system by the mean precision evaluated at a cut-off rank k of interest, which is called mean precision at k. Note that, to retrieve images in a large database, it may not be relevant to compare the query image to each image in the database: to speed up the retrieval process, some efficient search algorithms may be applied [START_REF] Quellec | Recherche de cas médicaux multimodaux à l'aide d'arbres de décision[END_REF][START_REF] Oliveira | Towards applying content-based image retrieval in the clinical routine[END_REF][START_REF] Dy | Unsupervised feature selection applied to content-based retrieval of lung images[END_REF][START_REF] Iakovidis | A pattern similarity scheme for medical image retrieval[END_REF].

Retrieval system Calibration

In order to design an optimal wavelet family within the lifting scheme framework, for our retrieval system, we first tune the important parameters jointly with usual wavelet families, and the optimal value of these parameters are then used in the wavelet optimization step. The important parameters to be tuned are:

• P 1 : the model used to represent the distribution of coefficients in each subband (either histograms or generalized Gaussians) and the distance between models (either the l1-distance or the Kullback-Leibler divergence), • P 2 : the number of decomposition levels N l , • P 3 : the decomposition scheme (either pyramidal, packet or SPACL decomposition).

Let P 4 be the wavelet family. a priori, the weights in the distance measure have to be tuned for each set of important parameters (as described in appendix A). Each parameter is evaluated independently, the other three are set to a default value; the following default values are used:

• for P 1 : generalized Gaussian model and Kullback-Leibler divergence,

• for P 2 : N l = 3,

• for P 3 : pyramidal decomposition,

• for P 4 : Daubechies 9/7 wavelet.

For the purpose of training, the database is divided into a training and an evaluating set. The weights are set in order to maximize the mean precision at k on the training set. Then, we compute the mean precision at k on the evaluation dataset, using the weights found on the training set. Finally, to design an optimal wavelet basis within the lifting scheme framework, we use the optimal values found for parameters P 1 , P 2 and P 3 , i.e. the values for these parameters that maximize the mean precision at k of the system, while the other parameters are set to their default values. Note that, to find the optimal wavelet basis, we have to conduct a two-stage optimization: each time a wavelet basis is evaluated, we search for the optimal weight vector, as illustrated in figure 5. The wavelet basis is also designed in order to maximize the mean precision at k on the training dataset. The retrieval method proposed below is evaluated on two medical databases.

It is also evaluated on a face database, to show the robustness of the method to illumination and slight orientation variations. The databases are first described and the results we obtained are given in section 7.4.

Databases

Diabetic Retinopathy Database (DRD)

This database has been built at the LaTIM laboratory 1 (Inserm U650) for research on diabetic retinopathy follow up. Diabetes is a metabolic disorder characterized by sustained inappropriate high blood sugar levels. This progressively affects blood vessels in many organs, which may lead to serious renal, cardiovascular, cerebral and also retinal complications. In the latter case, the pathology, namely diabetic retinopathy, can cause blindness. The database consists of 78 patient files: one for each diabetic patient seen at Brest University Hospital from June 2003 to April 2007 for an eye test. It contains 1261 photographs altogether, with associated anonymous information on the pathology. Images have a definition of 1280 pixels/line for 1008 lines/image. There are two possible orientations (right eye or left eye). Images have the same resolution. There are intra-image illumination variations but few interimage illumination variations. They were acquired by experts using a Topcon Retinal Digital Camera (TRC-50IA) connected to a computer and are lossless compressed. Four types of images are obtained: color, red free, blue-light images and the angiographic sequences. An image series is given in figure 6 as an illustration. The disease severity level, ranging from 0 to 5, according to ICDRS classification [START_REF] Wilkinson | Proposed international clinical diabetic retinopathy and diabetic macular edema disease severity scales[END_REF], was determined by experts for each patient.

Digital Database for Screening Mammography (DDSM)

The DDSM project [START_REF] Heath | Current status of the digital database for screening mammography[END_REF], involving the Massachusetts General Hospital, the University of South Florida and the Sandia National laboratories, has built a mammographic image database for research on breast cancer screening. It consists of 2500 patient files. Each one includes two images of each breast, along with some associated patient information. The database consists of 10000 images. These images have a varying definition, of about 2000 pixels/line for 5000 lines/image. There are four possible orientations (see figure 7). Images have the same resolution for a given acquisition device. There are intra-image illumination variations but few inter-image illumination variations. An example of image sequence is given in figure 7. Each patient file has been graded by a physician. Patients are then assigned one of these labels: normal, benign and cancer.

Face Database (FD)

The database consists of 400 images: ten photographs of 40 distinct subjects.

For some of them, the images were taken at different times, with different lightings, facial expressions (open / closed eyes, smiling / not smiling) and facial details (glasses / no glasses). All the images were taken against a dark homogeneous background with the subjects in an approximately frontal position.

Images have a definition of 92 pixels/line for 112 lines/image. The database can be viewed on the website2 . Figure 8 shows the 10 face images of a subject. Even though it is not our main purpose, this database can be classified easily. Images belong to the same class if and only if they represent the same subject.

Class Definition

Concerning the two medical databases, we are mainly interested in the medical grade of each image. As a consequence, we associate each possible grade with a class: six classes are then defined in DRD and three in DDSM. In the face database, we associate each person with a class, so we have 40 classes of ten elements. The goal of the system is thus to retrieve images belonging to the same class.

Objectives of the retrieval system

In a medical application, the retrieval system is intended to be used to assist the diagnosis of physicians. Such a system would be used as follows: if a physician has a doubt on his diagnosis, he could send the acquired images to the retrieval system, that will display the k most similar images along with their medical interpretations. The physician would then analyze the retrieved images and, reasoning by analogy, these images may help him confirm or invalidate his diagnosis for the new case. In order to limit the time required by physician to analyze the retrieved cases, ophthalmologists in Brest University Hospital want the system to return only the five most similar images for a query. For comparison purposes, the same number of images is displayed for the other two databases. As a consequence, we want the five retrieved cases to be as relevant to the query as possible: we will tune the method parameters in order to maximize the mean precision when k=5 images are returned by the system, i.e. the mean precision at k=5.

Results

The goodness of fit of the generalized Gaussian representation to the wavelet coefficient distribution, in each decomposition subband, is first reported for each of these databases: table 1 reports the average and standard deviation of the K-S statistic (D) measured for each subband across all the images in a dataset. Figure 9 illustrates the goodness of fit of the representation for one image. The relevance of the generalized Gaussian representation is high for high frequencies and decreases as the frequency decreases.

To assess the performance of the method, each database is divided into a training and an evaluation dataset. The sizes of these datasets are chosen according to the number of images available for each database. For DRD and FD, the training dataset consists of 50% of the whole database; on DDSM, it consists of 5% of the database. For both medical databases, the training and the evaluation datasets have the same proportion of images from each class. For these two databases, images from the same person are always assigned to the same dataset, either the training or the evaluation dataset. For FD, all images in a given sequence are assigned to the same set. The images are assigned to the two datasets at random, with respect to the constraints above.

The same training and evaluation datasets are used for each experiment. For wavelet adaptation, the size of filters P and U , respectively N p and N u , Mean precision at 5 of the retrieval method on the face database, as a function of the number adapted wavelet moments forced to vanish. Note that for M p = N p = M u = N u = 4, there is no wavelet adaptation. As an illustration, figure 11 depicts the distribution of images from DDSM in the feature space. The proposed method is compared to image retrieval based on color (or gray level) histograms and to [START_REF] Do | Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF]'s method, which are both particular cases of the proposed method. The color histogram based method is obtained by the proposed method without decomposition (N l = 0): the score obtained with this method is thus the first line of table 2. Do's method uses Daubechies 4-tap wavelets and a pyramidal transform; the same weight is assigned to each subband and the low-frequency subband is discarded. The proposed method is also compared to a retrieval method based on Zernike moments [START_REF] Khotanzad | Invariant image recognition by zernike moments[END_REF], in which the influence of each moment in the distance measure is weighted using an SVM [START_REF] Mladenić | Feature selection using linear classifier weights: interaction with classification models[END_REF].

M u =1 M u =2 M u =3 M u =4 M p =1
The results are reported on figure 12. For information, the precision-recall curves obtained for each evaluated value of the system parameters (P 1 , P 2 , P 3 and P 4 -see section 6) are given in figure 13. On medical databases, we observe that the different parameters only affect the precision-recall curves for small recall values. The reason is that the learning procedure focuses on small recall values: it is supervised by the mean precision at five. Note that we obtain low precisions when the recall is high, because there are actually several clusters within each class (see figure 11 for instance), so two images from the same class can be far from each other in the feature space. Finally, we study the distribution of the images selected for a query, among the predefined classes.

DISCUSSION AND CONCLUSION

In this article, two kinds of generic wavelet-based image signatures, with associated distance measures, have been evaluated in a CBIR system. They take advantage of the flexibility inherent in the wavelet transform framework to adapt the system to any specialized database. In particular, a way to adapt the wavelet transform to a high level criterion, within the lifting scheme framework, is proposed in this paper. It makes it possible to generate any wavelet transform, respecting the biorthogonality relations, and with a desired number of vanishing moments. A controlled random search, based on a genetic algorithm, is then performed in the predict and update filter space in order to find the optimal wavelet transform, and a similar search is performed in the as table 4 shows, a high mean precision score can still be achieved while eliminating a desired number of moments, which ensures a good compression rate.

It is particularly suited to JPEG 2000 standard part II [START_REF] Jpeg | Coding of still pictures -JPEG 2000 part2 ISO/IEC 15444-2[END_REF], which allows the user to define a custom wavelet basis and a custom decomposition scheme.

The system efficiency, measured with the mean precision when five images are returned by the system, i.e. the mean precision at five, is satisfactory: it reaches 56.50% on DRD, 70.91% on DDSM and 96.10% on FD. On DRD for instance, it means that on average two or three images out of the five retrieved by the system belong to the class of the query; and most images that do not belong to the right class are from nearby classes (see table 8). It is important to keep in mind that the scores of mean precision at five reported in this paper for medical datasets (56.50% on DRD and 70.91% on DDSM) are lower bounds on the actual performance of the system; since the system is evaluated through the use of predefined classes, it does not take into account the fact that images from different classes might be considered relevant by physicians (images at different severity scales can contain similar lesions for instance).

Compared with the reference methods, the improvement is notable (see figure 12): for instance, performance increases by approximately a third on DRD. In particular, the proposed image representation is more efficient than Zernike moments, which are known to be rotation-, scale-, and translation-invariant; however, these invariances are not really important, since we are mainly interested in textures as well as lesions, which are details of these images. It emerges from the precision-recall curves and the mean precision at five that the most influential parameters are the wavelet basis and the number of decomposition levels, which justifies our interest in finding an optimal basis within the lifting scheme framework. The optimization procedure always leads to a "better" wavelet basis than the standard bases evaluated. As for the decomposition scheme, it does not significantly influence the precision; as a consequence, we did not explore this trail any further. However, one could use the best basis algorithm [START_REF] Coifman | Remarques sur l'analyse de Fourier à fenêtre (French. English summary) [Remarks on windowed Fourier analysis[END_REF] to find an optimal decomposition scheme.

A higher performance was obtained on the three databases when wavelet distributions were compared using the Kullback-Leibler (KL) divergence rather than using the l1-distance (see table 2: GG / KL versus GG / l1). Using a generalized Gaussian (GG) representation is required to compute the Kullback-Leibler divergence: although this divergence was also defined for histograms, it is not suited for wavelet coefficient histograms that usually contain empty bins. However, we observe that representing the wavelet coefficient distribution by histograms does not lead to a high decrease of performance (see table 2: histograms / l1 versus GG / l1). From a computational point of view, both representations have their advantage. On one hand, computing the signatures is almost 20 times faster on average with the histogram based approach, although parameter estimation for generalized Gaussians could be approximated quickly from an adaptive histogram of each subband. But on the other hand, it is twice as fast to compute the distance measure using generalized Gaussians. As a consequence, the larger the database size is, the more adequate the generalized Gaussian based approach is; and the larger the query image definition is, the more adequate the histogram based approach is. Figure 11 illustrates another advantage of the generalized Gaussian representation: because the dimension of the feature space is much smaller than that of the histogram-based feature space, it is easier to vizualise it and extract new knowledge. In particular, in the feature space of the mammography dataset, depicted in this figure, we can see that images from three different devices are clustered, and within these clusters, semantically similar images are close to each other. It suggests that if a dataset is large enough, we may be able to find relevant images for query images of different resolutions and from different modalities. On the other hand, we can see that normal images from the top-right cluster are close to the abnormal images from the bottom-right cluster. It suggests that if the cardinality of the dataset was smaller, then only images from the same modality should be looked for, in order to avoid retrieving semantically irrelevant images.

Although the face database is not a medical database, it is very interesting.

Because we get a very high retrieval performance on this dataset, consisting of photographs of faces acquired with different illumination conditions and with slight orientation variations, we show that our method is robust to both these factors.

The adaptive framework proposed in this method turned out to be useful for the problems it was applied to. For instance, in DRD, more emphasis was given to the scales where the lesions are most visible. Moreover, in DDSM, the presence of lesions affects the distribution along specific directions: more emphasis was given to these discriminant directions (see table 5).

As mentioned in the introduction, the problem that remains to be solved, for such a CBIR system to be applied to any medical imaging modalities, is the normalization of the orientation and of the scale of images, when the variations are very large. However, the orientation normalization is mainly necessary in the presence of anisotropic textures or lesions, and the resolution of images is usually available in the metadata. Fortunately, if necessary, robust image registration techniques are available. As a conclusion, because it returns at least two or three relevant images out of five in the worst case (DRD), we believe this system is a useful tool that can help physicians reducing the search space when looking for images in a database clinically similar to an image acquired during an examination.

Fig. 1 .

 1 Fig.1. The two-channel filter bank for the lifting scheme, compared with the classic scheme. Both filters are equivalent.

Fig. 2 .

 2 Fig. 2. Wavelet transform of an image -Pyramidal decomposition schemeSuch a decomposition scheme is said to be pyramidal. Other decomposition schemes are conceivable, indeed it is possible to decompose other subbands of level (l -1) in addition to (l -1)LL, and the choice of the subbands to decompose can be different at each scale. JPEG 2000 standard part II[START_REF] Jpeg | Coding of still pictures -JPEG 2000 part2 ISO/IEC 15444-2[END_REF], for instance, lets us specify any decomposition scheme. In particular, the following four schemes are specified in the standard, they are illustrated in figure3:

  , we see that filter G can be expressed as a function of P . Indeed, d[n], the high-pass filtered version of the input signal x[n], satisfies the relations d[n] = G(x[n]) = x o [n] -P (x e [n]): as a consequence, if for instance P = (p 0 , p 1 , p 2 , p 3 ), then G = (-p 0 , 0, -p 1 , 1, -p 2 , 0, p 2 ) (see figure 4 (a)).

  Fig. 4. Equivalence between the classic and the lifting scheme based filter bank.

  which consists of the last M p columns in C. • Let B be the matrix containing the first N p -M p columns in C. The righthand side b is equal to 1 0 . . . 0 t -B. p 0 . . . p Np-Mp-1 t .

  Fig. 5. Optimization procedure

Fig. 6 .

 6 Fig. 6. Photograph series of a patient eye. Images (a), (b) and (c) are photographs obtained by applying different color filters on the camera lens. Images (d) to (j) form a temporal angiographic series: a contrast agent is injected and photographs are taken at different stages (early (e), intermediate (d), (f), (g), (h), (j) and late (i)). For the intermediate stage, photographs from the periphery of the retina are available.

Fig. 7 .

 7 Fig. 7. Mammographic image sequence of the same patient. (a) and (b) are images of the left breast, (c) and (d) are images of the right one.

Fig. 8 .

 8 Fig. 8. Image sequence of the same person's face

Fig. 9 .

 9 Fig. 9. Comparison of the generalized Gaussian model with the histogram model. The comparison was performed on a typical eye fundus red-free image (see figure 6 (b)).

Fig. 11 .

 11 Fig. 11. Distribution of images from DDSM in the feature space, for a particular subband (HH2), using the generalized Gaussian representation (image labels: green='normal', blue='benign', red='cancer'). Each cluster is associated with a particular acquisition device. Within these clusters, similar images are clustered (only approximately, because this figure depicts a projection of the feature space).

  • predict: we generate the wavelet coefficients d[n] as the error in predicting x o [n] from x e [n] using predictor operator P : d[n] = x o [n] -P (x e [n]). • update: combine x e [n] and d[n] to obtain scaling coefficients representing a coarse approximation to the original signal x[n]. This is done by using an update operator U : c[n] = x e [n] + U (d[n]).

		xe[n]		+	c[n]
		even/
	x[n]	odd	P	U
		split
			-		d[n]
		x o [n]

Table 1

 1 Goodness of fit of the generalized Gaussian representation to the Daubechies 9/7 wavelet coefficient distribution in each decomposition subband. The distribution of the K-S statistics (D), measured for each subband across all the images of a dataset, is reported: a small D indicates a good fit. We emphasized the subbands where an average D greater than 0.1 was measured. The mean precision at k=5 is given in table 2. To study the influence of the training dataset cardinality, we repeat the experiment with a training dataset consisting of 25% of the whole database for DRD and FD, and of 2.5% for DDSM. The scores obtained for the freely adapted wavelets (i.e. without forcing the moments to vanish) are reported in table 3. It emerges that the performance of the method increases with the size of the training dataset, we can thus expect an improvement as the databases will grow. We show in table 4 that we can also successfully adapt the wavelet if we eliminate some moments (see section 3.3). The experiment has been carried out on FD. The optimal weight vectors and P /U filters are reported in tables 5 and 6. The wavelet and scaling filters, respectively G and H, corresponding to filters P and U , are plotted in figure10. The mean computation time to retrieve the 5 closest images is given in table 7 for the two types of signatures: histogram based (H) or generalized Gaussian based (GG). All experiments were conducted using an AMD Athlon 64-bit based computer running at 2 GHz, with the default values for P 2 , P 3 and P 4 , for comparison purposes.

	database	DRD	DDSM	FD
		mean	std	mean	std	mean	std
	1HL	0.017 0.0027 0.014 0.0072 0.024 0.0066
	1LH	0.009 0.0034 0.017 0.0095 0.032 0.0098
	1HH	0.023 0.0070 0.027 0.0063 0.030 0.0102
	2HL	0.007 0.0032 0.012 0.0062 0.038 0.0161
	2LH	0.018 0.0043 0.017 0.0110 0.043 0.0258
	2HH	0.038 0.0108 0.026 0.0062 0.047 0.0233
	3HL	0.024 0.0066 0.009 0.0063 0.061 0.0239
	3LH	0.026 0.1004 0.023 0.0156 0.074 0.0328
	3HH	0.038 0.0227 0.019 0.0073 0.073 0.0329
	3LL	0.141 0.0358 0.225 0.0456 0.112 0.0323
	are both set to 4. At first, these filters only verify the biorthogonality relations
	(see section 3.2).

Table 2

 2 Mean and standard deviation (std) of the precision at five (in %) obtained when varying the value of one important parameter of the method while setting the others to their default value. The default value for each parameter is underlined. GG and KL stand respectively for generalized Gaussians and Kullback-Leibler divergence.

		database	DRD		DDSM	FD
			mean	std	mean	std	mean	std
		0	42.88 19.73 63.97 24.96 94.20	9.07
	number	1	51.58 18.96 67.43 24.23 95.10 3.24
	of decom-2	52.42 18.93 67.72 25.77 95.50	1.99
	position 3	54.17 18.78 68.23 26.24 93.10	6.87
	levels	4	54.12 18.89 68.37 25.82 86.40	9.54
		5	53.22 19.74 67.99 25.82 78.70 12.23
		Le Gall 5/3	52.90 19.37 68.17 25.53 89.20	8.01
		Daubechies 9/7 54.17 18.78 68.23 26.24 93.10 6.87
	wavelet	Haar	49.95 19.20 64.15 24.21 72.00 15.67
	basis	Cubic B-spline	53.30 18.77 66.34 24.71 85.70 10.84
		Daubechies 4	56.10 19.54 66.10 24.96 78.50 15.56
		Daubechies 6	56.14 19.95 64.95 25.76 80.00 12.35
	signature histograms / l1 51.77 18.73 69.75 25.53 86.40 15.49
	/	GG / l1	50.85 18.93 68.01 25.82 85.70 14.31
	distance GG / KL	54.17 18.78 68.23 26.24 93.10 6.87
	decom-	pyramidal	54.17 18.78 68.23 26.24 93.10 6.87
	position wavelet packet	51.36 19.97 68.63 26.32 85.70 14.20
	scheme	SPACL	53.81 19.22 68.55 25.57 91.20 12.93
	freely adapted wavelets
	(i.e. without forcing the	56.50 19.26 70.91 25.10 95.50 4.66
	moments to vanish)

Table 3

 3 Influence of the cardinality of the training dataset on the precision at 5 of the retrieval method based on freely adapted wavelets. The large (resp.

	small) training

Table 4

 4

  The repartition of the retrieved images among the different classes, as a function of the class of the query, is reported in table 8 for DRD.

			100
			mean precision at five	0 20 40 60 80	color histograms	DRD Do Zernike	proposed method	DDSM color histograms proposed method Do Zernike	color histograms	FD Do Zernike	proposed method
				Fig. 12. comparison with other methods
	Table 8
	Confusion matrix -repartition of the retrieved images among the different classes,
	as a function of the class of the query
							class of the query
				0		1		2	3	4	5
	class of the results	0 53.03% 5.41% 1 8.48% 57.88% 12.08% 10.71% 8.94% 5.60% 8.33% 6.59% 2 10.91% 11.53% 55.56% 10.48% 8.00% 3 10.91% 8.94% 12.08% 54.05% 7.06% 4 8.79% 9.41% 6.76% 7.62% 58.82% 10.49% 5.63% 7.57% 9.51% 7.18%
		5	7.88%	6.82%	7.92%	8.81% 10.59% 59.61%

http://latim.univ-brest.fr

www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html

distance weight vector space in order to maximize the precision of the system. The adapted wavelet transforms could also be used to store images in a compressed format, such as JPEG 2000, as their signature is computed. Indeed, Appendix A -Optimization procedure When we search for a weight vector between subbands or the most adapted wavelet, we have no preconceived ideas about either the fitness function or the free coefficients. Then we start the search by a controlled random process. Genetic algorithms [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF] are amongst the quickest and most popular methods for finding rough estimations of the best local maxima within a given region of space. The general principle of genetic algorithms is to:

• generate a random population of solutions • combine previous good solutions (crossover) and slightly modify some of them (mutation) to make the population evolve • keep the fittest solutions at each generation In particular, we used the steady state algorithm (at each generation, new individuals are created and added to the former population, the worst individuals of the resulting population being removed), with the following parameters:

• population size = 50 • maximum number of generations = 30 • selection methods: tournament selector (2 individuals are selected for crossover with a probability proportional to their score, and only the fittest is actually kept) • crossover probability (CR) = 70% (CR: probability for two selected individuals to combine their genes and form two new individuals, otherwise they are cloned) • mutation probability (MU) = 60% (MU: probability for a gene to be swapped with another after a crossover)

This configuration was adopted for its fast convergence. Once approximations of the fitness function optima are found by a genetic algorithm, we use the best approximations as initial points for descents to reach the actual local minima. We used the Powell's Direction Set method (Press et al., 1992b). The interest of this method is that it does not require computing the function gradient, which is unknown. The descent typically leads to an improvement of 2 or 3% in precision over genetic algorithms alone, so this step is optional.


Thérèse Hèrvée Mayi

Christian Duhem

Corinne Copin

Mohamed Amine Bouhlel

Elena Rigamonti

François Pattou

Bart Staels

Giulia Chinetti-Gbaguidi

Visfatin is induced by peroxisome proliferator-activated receptor gamma

  Figure1