Julien Dehos
email: julien.dehos@gmail.com

Éric Z Éghers
email: eric.zeghers@u-clermont1.fr

Laurent Sarry
email: laurent.sarry@u-clermont1.fr

Franc ¸ois ROUSSELLE Christophe Renaud
email: renaud@lisic.univ-littoral.fr

Practical photoquantity measurement using a camera

An image output by a camera is generally not a faithful representation of the real scene, because it undergoes a series of radiometric disturbances during the imaging process. This paper proposes a method for obtaining a more accurate measure of the light seen by a camera. Our method requires no specific calibration apparatus and only minimal supervision. Nevertheless, it is quite comprehensive, since it accounts for response function, exposure, vignetting, spatial non-uniformity of the sensor and colour balancing. Our method works in two steps. First, the camera is calibrated off-line, in a photoquantity sense. Then, the photoquantity of any scene can be estimated in-line. Our method is therefore geared to a wide range of computer vision applications where a camera is expected to give a measurement of the visible light. The paper starts by presenting a photoquantity model of the camera imaging process. It then describes the key steps of calibration and correction method. Finally, results are given and analyzed to evaluate the relevance of our approach.

Introduction

Computer vision handles various kinds of data, from videos to 3D medical images. This paper solely deals with 2D stills typically obtained from a computer-controlled colour camera. Such images can be used, for example, to reconstruct scene geometry [START_REF] Horn | Shape From Shading[END_REF][START_REF] Zhang | Shape from Shading: A Survey[END_REF], estimate the material reflectances [START_REF] Yu | Inverse global illumination: recovering reflectance models of real scenes from photographs[END_REF] or compute a wide-field view by image mosaicing [START_REF] Mann | Video Orbit of the Projective group: A New Perspective on Image Mosaicing[END_REF]. To achieve these tasks accurately, input images generally have to be a quite faithful measure of the light seen by the camera. However, this is not the case since the imaging process introduces geometric and radiometric distortions.

In this paper, we only deal with radiometric distortions. Our objective is to deploy a comprehensive yet practical method which can reduce radiometric distortions and thus give a quite faithful measure of the light seen by the camera.

We make several assumptions about the camera:

• electronic gain and optical aperture are constant,

• exposure time is known and controlled,

• optics can be temporarily removed from the camera.

Our method comprises two steps. First, the camera is calibrated. This step has to be done a single time, off-line. The correction step then computes the photoquantity of any static scene, in-line, based on the characterization of the camera given by the calibration step.

Our method accounts for:

• non-linear response functions (one response for each colour channel),

• colour balancing (global adjustment of the colour channels),

• spatial non-uniformity of the sensor (due to imperfections in the building process),

• vignetting (spatial non-uniformity due to the internal geometry of the camera).

Furthermore, since the method performs high dynamic range imaging (HDRI) through multiple exposures, it can be used with highly contrasted scenes and possesses some inherent robustness to noise.

The main contribution of this paper is a practical method to characterize a camera and then measure photoquantity. It also proposes an extension of the Goldman-Chen vignetting model and a calibration scheme based on this model. The remainder of the paper is organized as follows. Section 2 reviews some related work. Section 3 presents a photoquantity model of the imaging process. Section 4 details the model inversion, which gives a photoquantity measurement method. Section 5 details the calibration step. Section 6 presents and discusses some of our results and Section 7 concludes.

Related work

There is a wealth of literature on radiometric camera correction. However, work is generally limited to a specific kind of radiometric distortion (such as, vignetting or camera response), application (mosaicing, HDRI) or technology (CCD, CMOS).

In [START_REF] Tarel | Calibration Radiométrique de Caméra[END_REF], Tarel gives a fairly comprehensive study of radiometric camera distortion, but only considers CCD cameras leading to a specific characterization. Based on the study, he also explains how to radiometrically correct camera acquisitions. However, no practical method is given, as it was not the purpose of his paper.

High dynamic range imaging is a related field which has gained a lot of attention. In [START_REF] Mann | On Being 'undigital' With Digital Cameras: Extending Dynamic Range By Combining Differently Exposed Pictures[END_REF], Mann and Picard propose a method to estimate camera response from multiple exposures. Their method works as follows. Let us consider a point in the scene. This point emits a constant light flux. Therefore, by comparing the exposures and the corresponding values output by the camera, we can characterize the camera response for the measured light. Thus, by considering several points in the scene, the full response can be measured.

In [START_REF] Debevec | Recovering high dynamic range radiance maps from photographs[END_REF], Debevec and Malik propose another method to estimate the camera response from multiple exposures. Their method consists in considering the values of several points under different exposures. This makes it possible to define a minimization problem. The resolution of this problem gives the camera response.

In [START_REF] Robertson | Estimation-theoretic approach to dynamic range enhancement using multiple exposures[END_REF], Robertson et al. propose a similar method to Debevec-Malik. However, they use a different formalization, starting out from a statistical approach. They obtain a minimization problem, and propose a relaxation method to solve it.

Finally, in [START_REF] Mitsunaga | Radiometric Self Calibration IEEE Conference on Computer Vision and Pattern Recognition (CVPR)[END_REF], Mitsunaga and Nayar propose a polynomial model of camera response. Using this model, they formulate the response characterization as an optimization problem and propose an iterative algorithm to solve it.

There has been a recent surge in work on vignetting correction. The methods proposed are generally based on a vignetting model. For example, in [START_REF] Tu | Practical anti-vignetting methods for digital cameras[END_REF], Yu proposes a hypercosine model of vignetting and an antivignetting method based on wavelet denoising and decimation. In [START_REF] Zheng | Single-Image Vignetting Correction[END_REF], Zheng et al. use the Kang-Weiss model [START_REF] Kang | Can We Calibrate a Camera Using an Image of a Flat, Textureless Lambertian Surface?[END_REF] to propose a single-image vignetting correction method.

Vignetting correction is often associated with HDRI and image mosaics. In [START_REF] Goldman | Vignette and Exposure Calibration and Compensation[END_REF], Goldman and Chen propose a method to estimate vignetting, response function and exposures from an image mosaic. They propose a 6 th -order polynomial vignetting model and use the empirical model of response given in [START_REF] Goldman | Vignette and Exposure Calibration and Compensation[END_REF] to formulate the problem as a non-linear optimization.

In [START_REF] Kim | Robust Radiometric Calibration and Vignetting Correction[END_REF], Kim and Pollefeys use the same vignetting and response models as Goldman and Chen. However, they decouple vignetting from the response function and propose a method to successively estimate response function and vignetting by solving two linear least-square problems.

Finally, in [START_REF] Litvinov | Addressing Radiometric Nonidealities: A Unified Framework[END_REF][START_REF] Litvinov | Radiometric framework for image mosaicking[END_REF], Litvinov and Schechner propose a method to estimate vignetting, response function and exposures from a frame sequence. Their method is not based on a model. Instead, they use the fact that regions are represented on several frames to state and resolve a non-parametric linear least-square problem.

In contrast to previous methods, ours is both general and practical. First, it accounts for the main radiometric distortions occurring throughout the camera imaging process. Second, it is not limited to a specific kind of camera or application. Finally, it is practical in the sense that it requires no specific calibration apparatus and only little supervision.

Photoquantity model

In this section, we start with a recap on how a camera creates an image. This enables us to explain the main radiometric distortions, and thus to establish a photoquantity model of the imaging process.

Image formation

The digital imaging process is illustrated Fig. 1. First, the light travels from the real scene onto the optics of the camera. It then hits the cells of the light sensor where it is converted into electric signals. These signals are then amplified, digitized and possibly post-processed. This gives the output image.

In the case of colour imaging, the sensor is composed of several kinds of cell. Each cell type has its own spectral response. The colour of a pixel is given by associating cells of different kinds.

Radiometric distortions

The imaging process introduces radiometric distortions in the output image. The main distortion sources are described below.

Spatial non-uniformity

Spatial non-uniformity means that for the same input light, the camera can output different values according to the position reached in the image. This phenomenon is sometimes called vignetting, but has several forms and origins [START_REF] Goldman | Vignette and Exposure Calibration and Compensation[END_REF]. In this paper, we distinguish "sensor non-uniformity" and "vignetting".

We refer to "sensor non-uniformity" as the variation in sensitivity between cells of the sensor. This stems from the fact that the industrial sensor building process is not perfect, which means that two different cells exposed to the same irradiance can output different values.

We refer "vignetting" as the variation in light reaching a cell according to the position of this cell. Vignetting is due to the geometry of the camera optics, which are generally aligned to the center of the sensor. Thus, at the sensor center, the aperture of the optical system is maximal, enabling a lot of light to be collected, whereas in the corners of the sensor, the aperture is reduced, and only allows a little light to be collected, hence the well-known corner darkening.

Exposure and response

Another important radiometric distortion is the response of the imaging system. This means that the value of a pixel is not linear in relation to the light collected by the cells of this pixel. There are two reasons for this phenomenon.

Firstly, sensor cells have a limited dynamic range. If the input light is too low, a cell detects nothing while if the light is too high, the cell outputs the maximum value. Thus, the input light can only be quantified if it maps between these two thresholds.

Secondly, the cells may not give a linear response. This means that there is no linear function that maps the input light to the output value of a cell.

Colour balancing

As stated earlier, a sensor can associate different kinds of cells (colour channels) to quantify colour. Consider a pixel. We call l the spectral distribution function of the light reaching this pixel and r λ the spectral response of the colour channel λ of the camera. Thus, the value l λ output by the camera for this colour channel is :

l λ = r λ (x)l(x)dx
where x denotes wavelength.

This means that a camera uses its own, generally 3-dimensional, colour system. Therefore, the values output by a camera are not radiometric values because we cannot deduce the light spectrum from the values output by the camera. Furthermore, since the spectral responses of the camera generally do not exactly match the spectral responses of the CIE standard observer, the camera values cannot be expressed faithfully in a standard CIE colour system. For this reason, in [START_REF] Mann | Comparametric equations with practical applications in quantigraphic image processing[END_REF], camera values are called photoquantities.

The problem is that we cannot exactly transform the photoquantities output by a camera in a standard device-independent colour system. As a result, the colours in the output image may be different to the colours label we can really see.

Noise

Finally, as with any measurement devices, cameras are subject to noise. Noise can be seen as a random distortion of the output values. There are roughly two kinds of noise in a camera imaging system. The first is thermal noise, which depends on heat.

It can be seen as a biased additive noise and is therefore relatively important at low output levels. The second is measurement noise, which depends on the output level.

Model

Here is the photoquantity model we use. First, light comes from the scene to the camera. Then, the camera optics produce vignetting and the non-uniform sensitivity of the sensor introduces high-frequency spatial variations. Finally, the sensor responses and exposure introduce non-linearity between the incoming light and the output image. As it converts light into photoquantities, the sensor introduces colour balancing in the output image.

Let l p be the radiance coming from the scene to a pixel p in the output image. The camera measures this radiance using the corresponding sensor cells. Let r λ be the spectral response of the cells for the colour channel λ . The matching photoquantity l pλ is therefore:

l pλ = r λ (x)l p (x)dx
Let c pkλ n be the value given by the camera for the pixel p, the exposure k of exposure time t k and during the acquisition n. This gives:

c pkλ n = a λ f λ (u pλ v pλ t k l pλ ) + b pkλ n (1)
where v pλ is the vignetting at p, u pλ the sensor sensitivity at p, f λ the response of the colour channel λ , a λ the colour balancing factor for λ , and b pkλ n the noise.

The radiometric distortions are modeled as follows:

• the sensor responses fit the non-parametric model of Debevec and Malik,

• colour balancing is modeled by a colour-channel-dependent factor,

• sensor non-uniformity is modeled by a look-up table (LUT) giving the sensitivity attenuation of each pixel,

• vignetting follows an extension of the Goldman-Chen polynomial model.

Model inversion

In the previous section, we gave a photoquantity model of the imaging process. We now give a camera correction method by inverting the photoquantity model. The problem is to compute the photoquantity l pλ assuming the camera output c pkλ n is known.

Denoising

First, we have to remove noise (b pkλ n ). There is a wealth of research into how to denoise images [START_REF] Motwani | A Survey of Image Denoising Techniques[END_REF]. Here, we exploit the fact that we can easily control the camera. Indeed, we do not explicitly know the b pkλ n , but if we assume that thermal noise is negligible and that measurement noise is unbiased, then the denoised value of a pixel can be estimated by computing its mean value. Thus, we can compute the denoised value cpkλ from N values c pkλ n of the same scene (with constant exposure time and view direction) with:

cpkλ = 1 N N ∑ n=1 c pkλ n (2)
Note that a greater N value translates better noise reduction. However, this requires more image acquisitions. In practice, we can use N = 1 to focus on speed (during the correction step for example) and N > 1 to focus on accuracy (during the calibration step for example). We can thus rewrite Equation 1 as:

cpkλ = a λ f λ (u pλ v pλ t k l pλ ) (3)
In the following, we consider cpkλ = c pkλ and just write c pkλ .

Linearization

Now, we have to cancel the response functions and colour balancing. Our statement is:

c pλ = u pλ v pλ l pλ (4)
Therefore, Equation 3 is equivalent to:

c pkλ = a λ f λ (c pλ t k )
We cannot consider that any sensor is linear but we can consider that the more input light there is, the higher the output value of the sensor will be. We can thus assume that f λ is monotonous and write:

f -1 λ c pkλ a λ = c pλ t k i.e. ln f -1 λ c pkλ a λ = ln(c pλ ) + ln(t k ) (5)
We use the method proposed by Debevec and Malik in [START_REF] Debevec | Recovering high dynamic range radiance maps from photographs[END_REF] to robustly compute c pλ :

c pλ = exp   ∑ k∈K w(c pkλ ) ln f -1 λ c pkλ a λ -lnt k ∑ k∈K w(c pkλ )   ( 6
)
where w is a weighting function that discards under-exposed and over-exposed values:

w(c) = c -c min if c ≤ 1 2 (c min + c max ) c max -c
otherwise c min and c max are respectively the minimum and the maximum values the camera is able to provide. We may notice that the photoquantity c pλ is computed from several camera values c pkλ . Thus, c pλ has values in R where as c pkλ has values in [c min , c max ]. In other words, the images output by the camera are low dynamic range (LDR) images where as the images obtained after linearization are high dynamic range (HDR) images.

Uniformization

Finally, we have to cancel spatial non-uniformity. From Equation 4, we have:

l pλ =
c pλ u pλ v pλ which gives us the photoquantity l pλ we are looking for.

Calibration

In the previous section, we described a method to compute photoquantity based on the photoquantity model of the camera. In this section, we explain how to estimate the model parameters for a given camera. Our calibration method consists in acquiring a set of images from which response functions, colour balancing, sensor spatial non-uniformity and vignetting are then successively estimated.

Image acquisition

Acquisition scenes

To characterize a camera using our model, we have to take images of several scenes:

1. a scene to estimate the response functions, 2. a scene containing a white chart (such as a simple sheet of paper) to estimate colour balancing, 3. a scene captured from a number V of slightly different view directions to estimate vignetting. 4. a "irradiance scene" captured by the sensor, with no optics, to estimate sensor spatial non-uniformity, In practice, we can place a white chart in the real environment and take images from the V view directions. This enables us to estimate response functions, vignetting and colour balancing1 (items 1, 2 and 3). Finally, we have to use one more image from the scene, with no optics, to estimate sensor spatial non-uniformity (item 4). This means a total of V + 1 scenes is required.

Exposures

To estimate the camera response functions and measure the photoquantity of a scene, we have to capture this scene using several exposure times. Here, we consider that a scene is captured K times and we call t k the exposure time of the k th capture. In fact, K and the t k do not have to be the same for all the scenes considered.

Denoising

Since the calibration step is done once and its result is used many times, it can be profitable to take some time to improve calibration accuracy by reducing noise. This can be done by taking each image N times (with N > 1) and computing Equation 2.

Response functions

As stated earlier, we use the method proposed by Debevec and Malik in [START_REF] Debevec | Recovering high dynamic range radiance maps from photographs[END_REF] to estimate the response functions f λ . Knowing the response functions enables us to linearize the pixel values and then to estimate the other factors of radiometric distortion. Thus, we first consider each colour channel independently, which means that we can set a λ = 1. We set g λ (c) = ln f -1 λ (c) and rewrite Equation 5 with:

g λ (c pkλ ) = ln c pλ + lnt k (7)
From this equation, Debevec and Malik propose to minimize the objective function:

O λ = ∑ p∈P - k∈K (w(c pkλ )[g λ (c pkλ ) -ln c pλ -lnt k ]) 2 + β c max -1 ∑ c=c min +1 [w(c)g ′′ λ (c)] 2
where P -is a subset of the set of image points and g ′′ λ is the second derivative of g λ defined by:

g ′′ λ (c) = g λ (c -1) -2g λ (c) + g λ (c + 1)
The objective function O λ can be expressed as an over-determined quadratic system which can be solved by singular value decomposition. The first term of O λ forces the unknowns g λ and c pλ to match the measured values c pkλ and t k in Equation 7. The second term of O λ smoothes the computed response function. Finally, the parameter β defines the weight of the smoothing.

Colour balancing

As explained before, the colours given by the camera are influenced by the spectral responses of the different kinds of cell and by the spectral density of the light. In our method, we adjust colour channels according to the reference white point given by the white chart we have placed in one of the captured scenes. The corresponding pixel should have approximately the same values for all colour channels. If this is not the case, we select a channel as a reference (for example, the green channel) and apply a specific correcting factor to each other channel. Let [R w G w B w ] T be the linearized colour (i.e. the colour given by Equation 6) of the white point. We define the green-based colour balancing factors a λ by:

a r = G w /R w a g = 1 a b = G w /B w
Note here that the colour balancing factors can be included into the camera response functions. Thus, inverting response functions also performs colour correction.

Sensor spatial non-uniformity

Sensor spatial non-uniformity can be calibrated using the irradiance scene captured with no optics. Indeed, in these conditions, sensor cells are exposed to an overall irradiance, which remains roughly constant over the sensor.

Thus, to estimate sensor spatial non-uniformity, we can compute linearized colour-corrected colours in the scene captured with no optics. After normalization, this gives us the u pλ , which we can store in a LUT.

Vignetting

The last part of the model we have to calibrate is vignetting. We assume the vignetting v pλ of a pixel p = (p x , p y ) is:

v pλ = 1 + m 1λ R λ (p) + m 2λ R λ (p) 2 + m 3λ R λ (p) 3 where R λ (p) = m 4λ (p x -m 5λ ) 2 + (p y -m 6λ ) 2
This vignetting model is an extension of the polynomial model proposed by Goldman and Chen in [START_REF] Goldman | Vignette and Exposure Calibration and Compensation[END_REF]. This model presents an axis-aligned elliptical symmetry. The parameters m 1λ , m 2λ and m 3λ are the polynomial factors. The parameters m 5λ and m 6λ define the center of the symmetry. Finally, m 4λ is the scale factor of the elliptical symmetry. The vignetting model parameters can be determined by considering a scene from different view directions. Let s be a point of the scene and l isλ the photoquantity going from s to the camera in the view direction i. Since we know all the components of the photoquantity model of the camera except vignetting, we have:

c isλ = v isλ l isλ
where c isλ is the colour coming from s, seen by the camera in the view direction i, after linearization, colour correction and sensor uniformization. The corresponding vignetting factor is v isλ .

Similarily, in another view direction j, we have:

c jsλ = v jsλ l jsλ
Since the V scenes we use to estimate vignetting are seen from the same camera viewpoint but from different camera directions, we can consider that l ksλ is constant for all directions k. Therefore, we have l isλ = l jsλ and:

c isλ v jsλ -c jsλ v isλ = 0
We cannot estimate the vignetting model from this equation alone. However, if we consider several points and view directions, we have the following optimization problem:

arg min m 1λ ,m 2λ ,m 3λ ,m 4λ ,m 5λ ,m 6λ     ∑ i< j s∈S i ∩S j c isλ v jsλ -c jsλ v isλ    
where S k is a set of points visible from the camera in the direction k. To solve this optimization problem we use the Levenberg-Marquardt method [START_REF] Marquardt | An Algorithm for Least-Squares Estimation of Nonlinear Parameters[END_REF].

Note that after calibration, the vignetting model can be evaluated for all the pixels and included in the spatial non-uniformity LUT to optimize the correction step.

Note also that this vignetting calibration method requires us to identify the position of a scene point in the different images (taken from different view directions). This can be done automatically (using feature detection or geometric calibration) or manually.

Results and analysis

In this section, we present and analyze some of our results. We first consider the calibration step, then the correction step, and finally, we analyze our approach from different standpoints. We also give the computation times we obtained with our non-optimized Matlab implementation. We used a 1024x768 CCD colour camera (MATRIX VISION mvBlueFox-121C).

Calibration

Image acquisition

We recall that the calibration method requires us to consider V + 1 scenes. Each scene is captured using K different exposure times. For one scene and one exposure, N images are taken to reduce noise. Thus, the calibration method implies taking NK(V + 1) images. In the experiment presented here, we used exposure times ranging from 1 ms to 200 ms and took a total of 320 images (N = 10, K = 8 and V = 3). This may seem a lot, but with a computer-controlled camera, the calibration (image acquisition and model computation) can be performed in less than 5 minutes. Furthermore, with N = 4 and K = 4, we only need 64 images while still getting good results.

Fig. 2 presents selected denoised images of the calibration scene under different exposure times and view directions, and some denoised images of the irradiance scene (taken with no optics). With N = 10, denoising all our images takes less than 1 s.

Response Functions

Using K exposure times and one view direction, the calibration method then computes the camera response functions (10 s using 500 reference points) (Fig. 3) and the corresponding linearized high dynamic range image (10 s per image) (Fig. 4). Computing the response functions and the HDR images therefore takes less than 1 minute.

Note that a CCD sensor has a quasi-linear response function, which is not the case for a CMOS sensor.

Colour balancing

Using the white chart (the sheet of paper to the left of the rolls of adhesive tape in Fig. 4), the method then computes colour factors (1 ms). These factors are then included in the response functions (Fig. 3) to give the colour-corrected images (Fig. 4) (20 ms per image). Computing the colour-corrected images therefore takes less than 1 s.

Sensor spatial non-uniformity

Using the colour-corrected linearized image of the irradiance scene, the method computes a LUT (Fig. 3) (0.13 s) for the sensor spatial non-uniformity correction (Fig. 4) (20 ms per image). Computing the colour-corrected images therefore takes less than 1 s.

Vignetting

Finally, using the different view directions of the calibration scene, the method estimates the vignetting model for the camera and computes a correction LUT (Fig. 3) which is combined with the sensor spatial non-uniformity LUT (40 s with 10 reference points). This enables us to compute the photoquantity of the calibration scene (Fig. 4).

Figure 5: Images of a test scene under different exposure times (columns 1 to 3). Linearly tone-mapped photoquantity of the test scene obtained with our method (column 4).

Correction

Image acquisition

First, we acquire some images of the scene. Since we are considering a single view direction, we need NK images only (Fig. 5). We may notice that N, K and the exposure times can be different than the ones used during the calibration step. We may even state N = 1 while having good results because the correction method inherently performs noise reduction (this is detailed below).

Image correction

The previouly estimated response functions enable us to compute a linearized colour-corrected HDR image. Then, using the estimated LUT, we can uniformize this image. This results in the expected photoquantity (Fig. 5).

Analysis

In this section, we analyze the method from several standpoints: the accuracy of the linearization and colour-correction processes, the influence of noise and the accuracy of the vignetting modeling.

We call e p the relative error for a pixel p of an image P:

e p = ||c p -ĉp || 2 ||c p || 2 (8)
where c p is the reference colour and ĉp the computed colour. Let e max be the maximal value of the relative error, e mean the mean value of the relative error and e std the standard deviation of the relative error:

e max = max

Linearization and colour-correction

In order to analyze the linearization and colour-correction accuracy of the correction method, we use a scene containing a ColorChecker. A ColorChecker has a set of patches, reflectances of which are calibrated. The manufacturer gives the CIE-RGB values for the reflected light under the D65 CIE Standard Illuminant. Our scene is exposed to daylight conditions, which corresponds to this illuminant. Thus, the measured values should match the given values (ground-truth). However, the values measured by a camera rely on exposure time, sensor response. . . so that we cannot expect absolute equality but an affine relationship.

Thus, for each channel of each ColorChecker patch, we have the measured value and the ground-truth (i.e. a 2D point). We report all these points on a graph and compute the best affine function between the measured values and ground-truth (by linear regression). The measured values are accurate if the points are aligned (linear correlation).

Fig. 6 (a) shows the uncorrected LDR image, directly output by the camera (exposure time: 20 ms), and Fig. 7 (a), its correlation to ground-truth. We may notice that below 230, the measured values are quite well correlated to ground-truth. This can be explained by the linear response of CCD sensors (CMOS sensors give less correlated values because of their non-linear response). However, above 230, the sensor is over-exposed and the measured values do not match ground-truth (outliers). The graph also reveals colour shifting of the measured values (blue channel data points are below the regression line).

Fig. 6 (b) shows the corrected HDR image (linearly tone-mapped in a displayable range), given by our correction method, and Fig. 7 (a), its correlation to ground-truth. Here, the measured values are quite well aligned with no significant outliers (more generaly, the method also linearizes data from CMOS sensors). Moreover, the points of each colour channel are distributed on both sides of the regression line, which means that the method corrects colour shifting.

Noise

We first investigate how many images (i.e. the value of N) are required to significantly reduce the noise. Theoritically, the standard deviation of cpkλ drops with 1/ √ N (see Appendix A). In our experiments, setting N = 4 was sufficient to obtain e mean < 0.5% (e max = 2.4%, e mean = 0.4% and e std = 0.2%).

We also consider another question about noise: is the correction method inherently robust enough to use N = 1? We compared correction results obtained with N = 10 (c p in Equation 8) and with N = 1 ( ĉp ). The relative error e p is shown Fig. 8 (e max = 4.9%, e mean = 0.8% and e std = 0.4%).

Vignetting

Finally, we evaluate the accuracy of the vignetting modeling. Fig. 9 (a) shows the vignetting model that we find in our experiment (m 1 = 0.48, m 2 = 0.66, m 3 = -1.09, m 4 = 1.48, m 5 = 0.5 and m 6 = 0.69). This figure gives the quantity of light u p that reaches the pixel (p x , p y ) of the sensor after vignetting. We compare this result with the vignetting ground-truth (Fig. 9 (b)), obtained in experimental conditions (using an industrial uniform 20 × 30 cm light source). We find e max = 9%, e mean = 1.6% and e std = 1% (Fig. 9 (c)).

Other results

We apply the proposed method using different kinds of camera (CCD, CMOS) and lighting conditions (outdoor scene lit by the sun, indoor scene lit by a halogen lamp). To evaluate the method, we consider some reference points and compare the uncorrected values, the values corrected by the proposed method and the values corrected using Hugin2 [START_REF] Angelo | Radiometric alignment and vignetting calibration[END_REF]. The corrected values are defined in a high dynamic range therefore we compute correlation between output values and ground-truth (Table 1).

We may notice that our test scenes are not very constrasted which makes the uncorrected values quite representative of ground-truth, excepted when lighting conditions are not neutral (indoor scene). Hugin has many parameters and should give better results using a finer tuning. Finally, the method proposed in this paper gives good results, independently of lighting conditions, while requiring only minimal supervision. This means that standard deviation of the denoised value is σ pkλ / √ N and therefore drops with 1/ √ N.

Figure 1 :

 1 Figure 1: Camera acquision radiometric process.

Figure 2 :

 2 Figure 2: The calibration scene under different view directions (columns 1 to 3) and the irradiance scene (column 4), under different exposure times (first row: 2 ms; second row: 5 ms; third row: 20 ms).

Figure 3 :Figure 4 :

 34 Figure 3: Calibration of the mvBlueFox-121C camera.

p∈P (e p ) , e mean = 1 #P ∑ p∈P e p , e std = 1 #P - 1

 111 ∑ p∈P (e pe mean ) 2

Figure 6 :

 6 Figure 6: Study of the linearization and colour-correction accuracy of the method: uncorrected LDR image (a) and corrected HDR image (b).

Figure 7 :

 7 Figure 7: Correlation to ground-truth: uncorrected LDR image (linear regression: y = 0.7x -29.47, residual norm: 215.35) (a) and corrected (linearly tone-mapped) HDR image (linear regression: y = 1.09x -11.68, residual norm: 126.34) (b). Each point corresponds to one colour channel (indicated by the colour of the point) of one patch of the ColorChecker.

Table 1 :

 1 Correlation between output values and ground-truth (the higher the better).

		CCD outdoor CCD indoor CMOS outdoor CMOS indoor
	uncorrected values	0.85	0.63	0.84	0.67
	corrected values	0.87	0.84	0.9	0.86
	Hugin values	0.84	0.76	0.9	0.75

The light reflected by the white chart depends on the material and on the incident light. The chart should therefore be placed at a position where the light can be considered as a reference.

Image stitching software accounting for sensor response, vignetting and colour balancing.

Acknowledgements

This work is funded by the Audiovisual and Multimedia Program of the French Agency for Research (ANR-06-RIAM-001).

Conclusion

This paper described a practical method for measuring photoquantity using a camera. First, we proposed a photoquantity model of the imaging process. This model accounts for non-linear response functions, sensor exposure, colour balancing, sensor spatial non-uniformity and vignetting. Then, by inverting this model, we gave a photoquantity correction method which uses a small set of images together with the model parameters of the camera. We then gave a method to estimate these parameters using the camera. Finally, we presented and discussed some of our results. Our method is general as well as practical. It accounts for the main radiometric distortions occurring during the end-to-end camera imaging process. Furthermore, the method is not limited to a specific kind of camera or application. Finally, it is practical in the sense that it requires no specific calibration apparatus and only minimal supervision.

A Analysis of the denoising technique

We model the noise terms b n (Equation 1) as independent random variables with expected value 0 and standard deviation σ pkλ . Therefore, the c pkλ n are also independent random variables, with expected value c pkλ and standard deviation σ pkλ .

In Equation 2, we define cpkλ as:

Since the expected value is linear, we have:

This means that the denoised value effectively converges to the desired value. Finally, the variance of cpkλ is:


Youssef Stitou
email: youssef.stitou@ims-bordeaux.fr

Nour-Eddine Lasmar
email: nour-eddine.lasmar@ims-bordeaux.fr

Yannick Berthoumieu
email: yannick.berthoumieu@ims-bordeaux.fr

COPULAS BASED MULTIVARIATE GAMMA MODELING FOR TEXTURE CLASSIFICATION

Keywords: Image texture analysis, Information retrieval, wavelet decomposition, Gamma distribution, Gaussian copula

This paper deals with texture modeling for classification or retrieval systems using multivariate statistical features. The proposed features are defined by the hyperparameters of a copula-based multivariate distribution characterizing the coefficients provided by image decomposition in scale and orientation. As it belongs to the multivariate stochastic models, the copulas are useful to describe pairwise non-linear association in the case of multivariate non-Gaussian density. In this paper, we propose the d-variate Gaussian copula associated to univariate Gamma densities for modeling the texture. Experiments were conducted on the VisTex database aiming to compare the recognition rates of the proposed model with the univariate generalized Gaussian model, the univariate Gamma model, and the generalized Gaussian copula-based multivariate model.

INTRODUCTION

For numerous application domains, texture analysis is used in the workflow of image filtering, classification, segmentation, indexing, and/or synthesis. The main issue is to provide a unified and well-founded model of homogenous textures. When dealing with statistical texture modeling, many authors proposed to make use real or complex wavelet transform in a multiorientation and multiscale scheme [START_REF] Freeman | The design and use of steerable filters[END_REF], [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF]. This transformation which consists in decomposing an image into a set of oriented and scaled subbands captures the directionality, the structuredness, and coarseness of a texture. Consequently, approaches based on orientation and scale features imply to work in high-dimensional space. The topic concerning the selection of a stochastic model to characterize this augmented data-space is the main issue discussed in this paper.

Several works use marginal densities to characterize separately each subband. Subband coefficients can be described by univariate Gaussian (Gaud), generalized Gaussian (gGaud) or Gamma (Gamd) densities [START_REF] Do | Wavelet-based texture retrieval using generalized gaussian and Kullback-Leibler distance[END_REF], [START_REF] Mathiassen | Texture Similarity Measure Using Kullback-Leibler Divergence between Gamma Distributions[END_REF], [START_REF] Ves | A New Wavelet-Based Texture Descriptor for Image Retrieval[END_REF].

This kind of representation leads to a simple and attractive approach, defined by a limited set of hyperparameters, but this does not provide a complete statistical description of the texture images. Indeed, to take into account the interscale and intrascale statistical dependencies between subband coefficients it is necessary to develop a statistical multivariate framework. The main problem is that the direct multivariate extension of the above model, e.g. gGaud or Gamd, are generally not analytically defined [START_REF] Boubchir | Multivariate Statistical Modeling of Images with The Curvelet Transform[END_REF], [START_REF] Cho | Multivariate statistical approach for image denoising[END_REF], [START_REF] Chatelain | Change Detection in Multisensor SAR Images Using Bivariate Gamma Distributions[END_REF]. Thus, alternative models such as sub-Gaussian model, alpha-stable, Gaussian scale mixture or other non-Gaussian stochastic fields for joint and marginal statistics has been studied in previous works [START_REF] Tzagkarakis | Rotation-invariant texture retrieval with gaussianized steerable pyramids[END_REF] [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF]. However, none of these models is very reliable due to the high computational complexity for estimating hyperparameters and due to the very high number of associated hyperparameters.

In this paper, we focus on the multivariate statistical modeling using the copulas theory [START_REF] Sklar | Fonctions de répartition à n dimensions et leurs marges[END_REF], [START_REF] Joe | Multivariate models and dependence concepts[END_REF], [START_REF] Nelsen | An Introduction to Copulas[END_REF]. Copula models have become increasingly popular for multivariate modeling in many fields where the multivariate dependence is of great interest. Copulas are useful especially with non-Gaussian random variables and play an important role in developing a unified likelihood framework to analyze discrete or continuous stochastic processes. To develop a tractable multidimensional statistical models based on Gamd or gGaud, we address copula-based multivariate model which is able to describe the previous dependencies and the local structure inside each subband.

The paper is organized as follows. The next section provides the proposed statistical multivariate models for homogenous texture. A brief review and main proprieties of copulas is discussed followed by a description of the texture features corresponding to the proposed model. In section 3, experimental results are given to evaluate the retrieval performance.

PARAMETRIC TEXTURE MODELLING IN THE WAVELET TRANSFORM DOMAIN 2.1 A brief review of copulas

A copula is a multivariate cumulative distribution function (cdf) defined on the d-dimensional unite cube , , , ,

d d 1 1 d 1 x X x X P x x F ≤ ≤ = L L , n R x ∈ ∀ (1)
Then there exists a copula C such that

( ) ( ) ( ) (
).

, , , ,

d d 1 1 d 1 x F x F C x x F L L = (2)
Conversely, if C is a copula function and

d 1 F F , ,L
are cdf, then the function defined in ( 2) is a joint cdf with margins

d 1 F F , ,L
. This theorem demonstrates that joint distribution of a random vector of variables and the associated marginal distribution are necessary liked by a copula. Moreover if the function C is continuous and differentiable, then the copula density is given by

( ) ( ) . , , , , d 1 d 1 d n 1 u u u u C u u c ∂ ∂ ∂ = L L L (3)
Given that distribution F and copula C are absolutely continuous, the joint density function f is given as follows

( ) ( ) ( ) ( ) ( ) ∏ = = d 1 i i i d d 1 1 d 1 x f x F x F c x x f , , , , L L (4)
where i f is the density of i X and c is the density of the copula C . More details about copula theory can be consulted in [START_REF] Nelsen | An Introduction to Copulas[END_REF].

Proposed copula-based multivariate models

In this subsection we attempt to introduce a copula for characterizing simultaneously the marginal distributions, and the local structure of each subband. For this purpose, from each i th subband noted

( ) { } n m w i ,
we consider a ( ) ( )

1 q 1 p 2 + + 2 x
window. We concatenate neighbors' pixel in a column vector i W of size ( )( )

1 q 2 1 p 2 N + + = as follows ( ) ( ) T i i i q n p m w q n p m w ] , , , , [ + + - - = L W (5)
Under the spatial homogeneity assumption of each subband, observations of i W can be obtain by moving the window across the subband in an overlapping manner. To describe also the statistical interband dependencies, we consider the following d-dimensional random vector

T T B T 2 T 1 ] , , , [ W W W W L =
where B is the total number of bands, and each component i W is the N -dimensional vector described in [START_REF] Ves | A New Wavelet-Based Texture Descriptor for Image Retrieval[END_REF]. Thus, thanks to copula approach the multivariate density of the -

d dimensional vector W with NB d = is expressed as in (4). ( ) ( ) ( ) ∏ = = d 1 i i i d 1 d 1 f u u c f ω ω ω , , , , L L W (6) with ( ) ( ) d d d 1 1 1 F u F u ω ω = = , ,L , i F and i f
are respectively the univariate marginal pdf and cdf of W . Among a wide variety of copulas, we propose a Gaussian copula density defined as follows:

( ) ( )       - - = ∈ ∀ - 2 u c 1 0 u 1 T 2 1 d u I Σ u Σ Σ ẽxp , , ] , [ (7)
with u ~ being a vector of normal scores such that ( )

i 1 i u u - = Φ
~, and Φ is the cdf of the normalized Gaussian distribution. The matrix I implies the d dimensional matrix identity and Σ is the correlation matrix.. Note that the multivariate Gaussian pdf is a special case of ( 6) when all margins are univariate Gaussian. In this case the matrix Σ is the Pearson correlation matrix. If the margins are non-Gaussian, the ( )

k j, th element of Σ represents the linear correlation of two normal scores ( ) ( ) ( ) ( ) ( ) [ ] k k 1 j j 1 W F W F corr l j r - - = Φ Φ , , (8)
where Φ is the cdf of normal distribution ( )

1 0,

Ν

, j W and k W are respectively j th and the k th univariate component of W . According to the structure of the vector W , the correlation matrix Σ is a block matrix structured as follows

              = - B 1 B B 1 B 2 21 B 1 11 1 Σ Σ Σ Σ Σ Σ Σ Σ Σ ) ( L O O M M O L (9)
where i Σ is the correlation matrix associated to i th subband and ij Σ reflects the spatial crosscorrelation matrix between the i th and j th subbands. Moreover, since each band is homogenous, the diagonal matrices i Σ are structured as a ( ) ( )

1 p 1 p + × +
block Toeplitz matrices with Toeplitz blocks as follows

              = - - - i 0 i 1 i p i 1 i 0 i 1 i p i 1 i 0 i Σ Σ Σ Σ Σ Σ Σ Σ Σ Σ L O O M M O L (10) Each submatrix i k Σ is a ( ) ( ) 1 q 1 q + × + Toeplitz one as ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )             - - - = 0 k r 1 k r q k r 1 k r 0 k r 1 k r q k r 1 k r 0 k r i i i i i i i i i i k , , , , , , , , , L
O O M M O L Σ (11)
where the term ( ) l k r i , is the spatial correlation of the ith subband at the lags ( ) l k, . Each matrix ij Σ has the same structure as i Σ and contains the spatial crosscorelation between the i th end j th subbands. As mentioned in the previous section, it is well-known that a univariate gGaud or Gamd would be a reasonable statistical representation in the multi-resolution approach. These univariate models are simple to fit and capable of characterizing the first order statistics properties of the band coefficients and its magnitudes respectively. allows us to cover Using copula theory, multivariate versions of these models can be obtained as follows

Multivariate Gamma model

To model all bands coefficients magnitudes by a multivariate Gamma model we use equation ( 4) with a Gaussian density copula in [START_REF] Cho | Multivariate statistical approach for image denoising[END_REF] and the univariate Gamd defined by

( ) ( ) ( ) a b e x x f R x a b 1 a Γ θ γ - - + = ∈ ∀ , , , (12)
where 0 a > , is a shape parameter, and 0 b > is related to the scale distribution, and ( ) . Γ is the Gamma function, defined as follows

( ) ∫ = ∞ - - 0 1 z t dt t e z Γ , 0 z > ∀ ( 13
)

Multivariate Generalized Gaussian Model

The multivariate dGaud can be obtained using the Gaussian copula and the univariate dGaud as marginal. The univariate dGaud is given by

( ) ( ) ( ) β α β Γ α β β α x e 1 2 x f R x - = ∈ ∀ , , , (14)
where α is a shape parameter, and β is related to the scale of the distribution. This marginal modeling allows us to cover super gaussian

) ( 2 < β , Gaussian ) ( 2 = β and a sub gaussian ) ( 2 > β densities.

Parameters estimation and Feature extraction

In a parametric approach, the feature extraction step consists of estimating the model parameters. For multivariate copula-based model, the maximum likelihood criterion for the model in (4) requires two stages of estimation methods. The first stage consists in estimating the parameters of marginal distribution by the maximum likelihood estimator. In the second stage, these parameters are used to transform the data in the unit hypercube and involve a maximum likelihood of copula matrix. If we consider the coefficients of different subbands as independent, the correlation matrix of the Gaussian copula is a block diagonal matrix as follows

            = B 2 1 0 0 0 0 0 0 Σ Σ Σ Σ L O O M M O L (15)
Furthermore, the the joint pdf associated to the ith subband is given by

( ) ( ) ( ) ( ) ( )
∏ = = n 1 k i k i i n i 1 i i i f F F c f i i ε Σ ε Σ W W , , , , , , ω ω ω ω
L ( 16
)
where i W c is a Gaussian copula density parameterized by the matrix i Σ , i f is the univariate marginal pdf of the ith subband characterized by the parameters i ε and i F is the related cdf. As a result, the signatures at each image constitute a vector Θ which contains the marginal parameters i ε and a subset of elements from the matrix i Σ of each band.

EXPREMENTAL RESULTS

To compare the performances of the presented texture models, we now use the benchmark proposed by Do et al. in the framework of texture retrieval [START_REF] Do | Wavelet-based texture retrieval using generalized gaussian and Kullback-Leibler distance[END_REF]. The retrieval scheme is applied on a set of texture images obtained from the MIT Vision Texture (VisTex [START_REF] Nelsen | An Introduction to Copulas[END_REF]) database. From each of these texture images of size 512x512 pixels, 16 subimages of 128x128 pixels are created. A test database of 640 texture images is thus obtained. The objective of our experiments is to compare the retrieval performances of the following models:

-MGamdGC: Multivariate Gamd with Gaussian copula.

-MgGaudGC: Multivariate gGaud with Gaussian copula.

-UGamd: Univariate Gamd -UgGaud: Univariate gGaud For this purpose we use the steerable pyramid decomposition proposed in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF], [START_REF] Do | Wavelet-based texture retrieval using generalized gaussian and Kullback-Leibler distance[END_REF] with 2 Nsc = levels and 6 Nor = orientations. The univariate and multivariate generalized Gaussian models are used to fit only the real part of the steerable pyramid components. For the univariate approaches UGamd and UGaud, the marginal parameters of each band are estimated using the ML technique described respectively in [START_REF] Do | Wavelet-based texture retrieval using generalized gaussian and Kullback-Leibler distance[END_REF] and [START_REF] Mathiassen | Texture Similarity Measure Using Kullback-Leibler Divergence between Gamma Distributions[END_REF]. For each model the estimated parameters are concatenated to form a feature vector

[ ] B B 1 1 UGamd b a b a , , , , L = Θ (17) [ ] B B 1 1 UgGaud β α β α ; , ; ; L = Θ (18) where Nsc Nor B × =
is the total number of bands. For the two multivariate models the feature extraction is carried out using the autocorrelation samples within each matrix i Σ . In this experiments, the parameters p and q are fixed such as 3 q p = = . Thus, the matrix i Σ contains only the following

1 4 × vector ( ) ( ) ( ) ( ) [ ] 1 1 r 1 1 r 0 1 r 1 0 r i i i i i , ,
[ ] UGamd d 1 MgGamdGC Θ R R Θ , , ,L = (20) [ ] UGaud d 1 MgGaudGC Θ R R Θ , , ,L = (21)
We note that if we use a small window of size 3 3 × i.e.

1 q p = =
, the number of the extracted features for multivariate models is given by ( )

2 q 2 p Nor Nsc N + = . . Θ (22)
In order to determine similarity between two images in the database, we use the normalized L1 distance between the feature vectors 1 Θ and 2 Θ , i.e.,

∑ = - = Θ σ Θ Θ N 1 k 2 1 2 1 k k k d Θ Θ , ( ) ( ) ( ) ( )
where ( ) k σ refers to the standard deviation of the component k of all the feature vectors in the database. In the retrieval stage a query image is any one in the database. The relevant images for each query are the other 15 images obtained from the same original 512x512 image. The performance of each model is evaluated in the average percentage of retrieving relevant images. In Table 1, we provide a comparison of performances when the signatures are computed by the four previous methods. The maximum and minimum percentage rates were 73.4% and 83.4%. From this table we first observe that the proposed multivariate models improve the performances of the retrieval system than univariate models. Second, from these results the copula based multivarite Gamma model has about 3% improved retrieval rate than generalized Gaussian multivariate model. .

UgGaud

CONCLUSION

In this paper we have addressed the problem of texture image modeling in a multiorientation and multiscale scheme. A multivariate Gamma model using Gaussian copula is proposed and its performances are evaluated for texture classification. The experiments results show that the multivariate model driven by copula presents performances in the retrieval system better than that of marginal modeling. However, the Gaussian copula used in this work is just a solution among others. As a perceptive we intend in future work to use other copulas with more flexible dependence.

Table 1 :

 1 Average retrieval rates (%) comparison

			UgGamd	MgGauGC	MgGamGC
	Average retrieval rate	73.4375	77.9785	80.3625	83.4543
	(%)
	Θ N	24	24	72	72


Nour-Eddine Lasmar
email: nour-eddine.lasmar@ims-bordeaux.fr

Youssef Stitou
email: youssef.stitou@ims-bordeaux.fr

Yannick Berthoumieu
email: yannick.berthoumieu@ims-bordeaux.fr

MULTISCALE SKEWED HEAVY TAILED MODEL FOR TEXTURE ANALYSIS

Keywords: Image texture analysis, Asymmetric Generalized Gaussian density, Kullback-Leibler Divergence, Dual-Tree Complex Wavelet Transform, Texture Retrieval

This paper deals with texture analysis based on multiscale stochastic modeling. In contrast to common approaches using symmetric marginal probability density functions of subband coefficients, experimental manipulations show that the symmetric shape assumption is violated for several texture classes. From this fact, we propose in this paper to exploit this shape property to improve texture characterization. We present Asymmetric Generalized Gaussian density as a model to represent detail subbands resulting from multiscale decomposition. A fast estimation method is presented and closed-form of Kullback-Leibler divergence is provided in order to validate the model into a retrieval scheme. The experimental results indicate that this model achieves higher recognition rates than the conventional approach of using the Generalized Gaussian model where asymmetry was not considered.

INTRODUCTION

The accurate characterization of texture is fundamental in various image processing applications, ranging from retrieval in large image databases to segmentation and texture synthesis.

Over the last decade, numerous works have been proposed for texture modeling using multiscale image representations [START_REF] Van De Wouwer | Statistical texture characterization from discrete wavelet representations[END_REF] [2] [START_REF] Do | Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF]. The conventional scheme of multiscale texture analysis consists of extracting features from subband coefficients and uses these features as a signature for a specific texture class. Portilla and Simoncelli presented a statistical model based on joint statistics of steerable pyramid coefficients [START_REF] Portilla | A Parametric Texture Model Based on Joint Statistics of Complex Wavelet Coefficients[END_REF]. In their work, efficient algorithm of texture synthesis was developed giving increased synthesis quality. However this model is not tractable for classification applications due to the largeness of the signature. Some others conventional approaches consist in representing detail subband coefficients by their probability density functions (PDFs). In [START_REF] Van De Wouwer | Statistical texture characterization from discrete wavelet representations[END_REF] [3] Generalized Gaussian (GG) density was used as a parametric model for subband marginal PDFs where orthogonal wavelets are employed as a multiscale decomposition. In a related recent work [START_REF] Kwitt | Image similarity measurement by Kullback-Leibler divergences between complex wavelet subband statistics for texture retrieval[END_REF], Weibull distribution was proposed to model the detail subbands magnitudes employing Dual-Tree Complex Wavelet Transform (DT-CWT).

The GG distribution was introduced to model the detail subband coefficients because of the almost heavy tailed aspect of the empirical PDFs of these coefficients. It is adequate for Super Gaussian (leptokurtic), Gaussian (mesokurtic) and sub Gaussian (platykurtic) behavior of the marginal distribution of subband coefficients. This tail aspect is directly linked to the Kurtosis which shouldn't be the only statistic taking into account. However, Skewness is also important for modeling textures [START_REF] Portilla | A Parametric Texture Model Based on Joint Statistics of Complex Wavelet Coefficients[END_REF] and GG density presents limitations for fitting asymmetric subband distribution because its Skewness is null by definition. In order to take into account simultaneously the two statistics (Kurtosis and Skewness) we present a new model using a skewed heavy tailed distribution in the framework of texture analysis: the Asymmetric Generalized Gaussian density (AGG) [START_REF] Tesei | The Asymmetric Generalized Gaussian Function: A New HOS-Based Model for Generic Noise Pdfs[END_REF], [START_REF] Tesei | HOS-based generalized noise pdf models for signal detection optimization[END_REF], [START_REF] Lee | Maximum likelihood parameter estimation of the asymmetric generalised Gaussian family of distributions[END_REF].

In this paper, we work with the Dual-Tree Complex Wavelet Transform (DT-CWT) proposed by Kingsbury [START_REF] Kingsbury | Complex Wavelets for Shift Invariant Analysis and Filtering of Signals[END_REF]. However others multiscale decompositions can be used (e.g. steerable pyramid [START_REF] Portilla | A Parametric Texture Model Based on Joint Statistics of Complex Wavelet Coefficients[END_REF], orthogonal wavelets [1] [3]).

The organization of this paper is as follows: we review in section 2 the AGG density model and introduce a new fast deterministic method to estimate the marginal distribution parameters. In section 3 the Kullback-Leibler divergence between AGG densities is derived. Finally, experimental results indicate the significant improvement in retrieval rates using the new model.

ASSYMETRY MODELING

PDFs proposed in previous works rely on symmetric assumption for the wavelet coefficient histogram [1] [3]. In our knowledge, no parametric models explicitly incorporate possible non null texture Skewness. Nevertheless, Van de Wouwer et al [START_REF] Van De Wouver | Multiscale asymmetry signatures for texture analysis[END_REF] defined a multiscale asymmetry signature consisting in an extra texture descriptor defined by: . Experimental results [START_REF] Van De Wouver | Multiscale asymmetry signatures for texture analysis[END_REF] confirmed that adding asymmetry signatures resulted in an improved texture characterization. However, this non-parametric model has some limitations; it is bins-dependant and signature estimation uses histogram computing and integral approximation which need more computation.

du u u h u h A ni ni ni ∫ +∞ - - = 0 ) ( ) ( (1
Overall, Skewness is a universal measure of the asymmetry of probability distribution. We present a model that explicitly incorporates possible texture asymmetry. The AGG model is a generalization of the GG one and can track not only the forth high order statistic (HOS) but also the third HOS: the Skewenss. Figure 1 shows how AGG density is more suitable than GG density to respect joint (Skewness,Kurtosis) information. We represent a scatterplot of the two statistics calculated from 1level DT-CWT coefficients of a textured images set. The GG domain is represented by the vertical dashed straight line. When GG can only fit zero-Skewness subbands, AGG have by definition a huger domain including the GG one and can de facto better model the skewed subbands.

Skewed GG density: the AGG

The PDF of the AGG density with zero mode is given by [START_REF] Tesei | The Asymmetric Generalized Gaussian Function: A New HOS-Based Model for Generic Noise Pdfs[END_REF] [6]:

( ) ( ) ( ) ( ) ( )        ≥         - Γ +         - - Γ + = 0 exp 1 0 exp 1 , , ; x x x x x f r r l l r l r l α α β α β β α β α β β α β β α (2)
where ( ) The Skewness of this distribution can take any real value (negative, null or positive):

∫ +∞ - - = Γ 0 1 dt e t x
( ) ( ) ( ) ( ) ( ) ( ) 3 3 3 3 4 4 3 1 4 α β β α β β α β β Γ + Γ + Γ - = l r l r l r s (3)
The AGG distribution was introduced by Tesei and Regazzoni [START_REF] Tesei | The Asymmetric Generalized Gaussian Function: A New HOS-Based Model for Generic Noise Pdfs[END_REF] in order to optimize signal detection in non-Gaussian environments. It is able to describe many kinds of wavelet coefficients; symmetric, asymmetric, with variable sharpness of marginal PDFs (Figure 2).

Second Order Statistics (SOS) estimation

To estimate AGG density parameters Maximum Likelihood (ML) method was proposed [START_REF] Lee | Maximum likelihood parameter estimation of the asymmetric generalised Gaussian family of distributions[END_REF]. It uses several steps to find the best value of the shape parameter α . This method provides a better convergence rate for reaching parameter value but does not guarantee convergence for bad initial value. So, ML estimation is inherently complex and needs a lot of computations. We now develop a simple method that allows us to find model parameters given the one-order absolute moment and the second-order moment: SOS estimation. We note that the proposed estimate can be used to initialize the ML approach.

We first calculate the k-order moments k µ and k-order absolute moments k m in the case of AGG density with parameters ) , , (

r l β β α ; ) ( k k X E = µ and ) ( k k X E m =
. By using the integral calculation [START_REF] Abramowitz | Handbook of Mathematical Tables[END_REF]:

∫ +∞ - - -       Γ = 0 1 1 γ ν µ γ γ ν µ ν γ x e x ( 4
)
we establish that

( ) ( ) ( ) α α β β β β µ 1 ) 1 ( 1 1 1 Γ + Γ × + - + = + + k l r k l k k r k ( ) ( ) α α β β β β 1 ) 1 ( 1 1 Γ + Γ × + + = + + k m l r k l k r k for 2 , 1 = k
We use the second order statistic 2 µ with absolute moment 1 m to calculate the ratio . We find:

( ) ( )( ) ( ) α ρ γ γ γ × + + + = 1 1 1 3 2 2 r ( 5
)
where

r l β β γ = and ( ) ( ) ( ) ( ) α α α α ρ 3 1 2 2 Γ Γ Γ =
is the generalized Gaussian ratio function.

The value of γ can be estimated by [6][7]:

∑ ∑ ≥ = = - - = r k l k N x k k r N x k k l x N x N 0 , 1 2 0 , 1 2 1 1 1 1 ˆ γ (6)
where l N ( r N ) is the number of k x samples 0 ( 0 ≥ ).

An unbiased estimate of r is:

( ) ∑ ∑ = 2 2 ˆk k x x r (7)
A very good approximation of ( ) α ρ and its inverse was proposed in [START_REF] Domínguez-Molina | A practical procedure to estimate the shape parameter in the generalized Gaussian distribution[END_REF], that leads to a fast and deterministic shape parameter estimation. So, we propose to estimate AGG parameters using the following algorithms:

1-we calculate R ˆusing γˆand r ˆestimation above, ( 6) and ( 7):

( )( ) ( ) 2 2 3 1 ˆ1 1 ˆ+ + + = γ γ γ r R (8)
2-According to R ˆ value we estimate α using the approximation of the inverse generalized Gaussian ratio [START_REF] Domínguez-Molina | A practical procedure to estimate the shape parameter in the generalized Gaussian distribution[END_REF]:

( ) R 1 - = ρ α

SIMILARITY MEASURMENTS BETWEEN AGG DISTRIBUTIONS

The conventional scheme of multiscale texture analysis consist in extracting features from subband coefficients and use them as a signature associated with an appropriate distance to distinguish several texture classes. Even that several probabilistic distances can be used, Do and Vetterli justified the use of Kullback-Leibler divergence (KLD) to have relevant results in retrieval application [START_REF] Do | Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF].

Considering two AGG PDFs, ) ; ,; (

1 1 1 1 r l x agg β β α and ) ; , ; ( 2 2 2 2 r l x agg β β α
, KLD is defined by:

( ) ( ) ( ) ( )         = ∫ r l r l r l x agg x agg x agg agg agg KLD 2 2 2 2 1 1 1 1 1 1 1 1 2 1 , , ; , , ; ln , , ; , β β α
β β α β β α
We calculate a closed-form expression of KLD using the integral expression (4), we find:

( ) ( ) ( ) ( ) 1 1 1 2 1 1 2 2 1 1 2 1 1 2 1 1 1 2 2 2 2 1 2 1 1 1 1 1 1 ln , 2 2 2 2 2 2 α α α α β β β β β β β β α β β α α β β α α α α α α α -         Γ         + Γ × + + +                       Γ +         Γ + = + + r l r l r l l r r l r l agg agg KLD

EXPERIMENTAL RESULTS

The experiments give an evaluation of the proposed model in the framework of texture retrieval. We use almost the same experimental setup presented in [3] [4]. We work with 40 texture classes from Brodatz album [12] (Figure 3). From each of these texture images of size 640x640 pixels, 16 subimages of 160x160 are created. A test database of 640 texture images is thus obtained. A query image is any one of these images in the database. The relevant images for each query are the other 15 images obtained from the same original 640x640 image. We must note that this dataset was not selected in any way to contain textures that showed specific asymmetry. Statistical analysis of all computed asymmetry parameter mentioned above [START_REF] Van De Wouwer | Statistical texture characterization from discrete wavelet representations[END_REF] showed that the 90 th quantile is equal to 9.68 10 -5 . The Kingsbury's DT-CWT is employed to have a multiscale image representation with six complex detail subbands at each level. Advantages of the DT-CWT can be resumed in its nearly shift invariance and directional selectivity. We use Kingsbury's Q-Shift (14,14)tap filters in combination with [START_REF] Kingsbury | Complex Wavelets for Shift Invariant Analysis and Filtering of Signals[END_REF]19)-tap near-orthogonal filters [START_REF] Kingsbury | Complex Wavelets for Shift Invariant Analysis and Filtering of Signals[END_REF]. First, we tested the impact of using a skewed model by comparing the performance of retrieving. The conventional criterion on recall/precision shows the improvement obtained with the skewed model AGG (Figure 4 (a)).

Second, we compare the performances of the proposed skewed model with those obtained using Weibull distribution to model DT-CWT coefficients magnitude [START_REF] Kwitt | Image similarity measurement by Kullback-Leibler divergences between complex wavelet subband statistics for texture retrieval[END_REF]. As we can see in Figure 4 (b), AGG model outperforms Weibull model. But we must note that Weibull model has three times less parameters to be estimated. It can be seen from the Figure 5 that, compared to the two others models, the AGG converge faster. For example, we retrieve 94% of the relevant images with 20 retrieved images when we require 24 ones to have the same percentage if we employ the two others models (GG and Weibull). A second series of experiments was conducted using orthogonal wavelets with Daubechies' D4 filters. The results of average retrieval rates are summarized in Table 1. As expected, the use of AGG model improves the performance of retrieval and in consequence texture characterization.

Retrieval

CONCLUSION

In this work we have addressed the issue of modeling textures in multiscale scheme. The probability density functions (PDFs) of detail subband coefficients are skewed and heavy tailed. So, Asymmetric Generalized Gaussian model is proposed in order to take into account third and fourth Higher Order Statistics. The model is validated in the retrieval system and achieves better recognition rates compared to GG-based approach. Future work includes the AGG parameter estimation by Negentropy Matching and the model extension to multivariate joint PDFs.

Fig. 1 .

 1 Fig.1. Scatterplot of subbands' (Skewness,Kurtosis) couple

  left-scale and rightscale parameters respectively. AGG density extends GG when β β β = = r l .

7 Fig. 2 .

 72 Fig.2. PDFs differentiated by Skewness and with same shape parameter

  we estimate left and right scale parameters:

Fig. 3 .

 3 Fig.3. Selection of images from Bordatz album: from top to bottom and left to right: D6, D8, D11, D16, D17, D18, D20, D21, D24, D26, D28, D32, D33, D34, D47, D49, D51, D53, D55, D56, D57, D64, D65, D66, D75, D76, D77, D78, D82, D83, D84, D85, D87, D92, D95, D101, D102, D103, D104, D110.

Fig. 4 .

 4 Fig.4. Recall-Precision curves showing the impact improvement of skewed model, with DT-CWT and 1 scale

Fig. 5 .

 5 Fig.5. models convergence comparison

Table 1 :

 1 Average retrieval rate (%) comparison

			Rate
	Number of scales	AGG	GG	Weibull
	1	87.9980	87.2656	87.4023
	2	93.7988	91.8945	92.4316
	3	98.3301	94.9805	94.5898


Nour-Eddine Lasmar
email: nour-eddine.lasmar@ims-bordeaux.fr

Yannick Berthoumieu
email: yannick.berthoumieu@ims-bordeaux.fr

MULTIVARIATE STATISTICAL MODELING FOR TEXTURE ANALYSIS USING WAVELET TRANSFORMS

Keywords: image texture analysis, multivariate model, Kullback-Leibler Divergence, wavelet transforms, information retrieval

In the framework of wavelet-based analysis, this paper deals with texture modeling for classification or retrieval systems using non-Gaussian multivariate statistical features. We propose a stochastic model based on Spherically Invariant Random Vectors (SIRVs) joint density function with Weibull assumption to characterize the dependences between wavelet coefficients. For measuring similarity between two texture images, the Kullback-Leibler divergence (KLD) between the corresponding joint distributions is provided. The evaluation of model performance is carried out in the framework of retrieval system in terms of recognition rate. A comparative study between the proposed model and conventional models such as univariate Generalized Gaussian distribution and Multivariate Bessel K forms (MBKF) is conducted.

INTRODUCTION

The accurate characterization of texture is fundamental in various image processing applications, ranging from retrieval in large image databases to segmentation and texture synthesis.

Recent works are shown that wavelet-based statistical characterization is a very efficient solution to describe textures for texture application. The conventional scheme of multiscale texture analysis in wavelet domain consists in modeling subband coefficient distributions by parametric density functions as a signature for a specific texture class.

Many univariate prior models such as the Generalized Gaussian distribution (GGD) [START_REF] Do | Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF] and the Student tdistribution model [START_REF] Huang | Statistics of natural images and models[END_REF] have been used to successfully characterize the marginal distribution of subband coefficients. In order to model magnitude coefficients, Gamma and Weibull distributions have been introduced to achieve better retrieval rates [3] [4]. This kind of representations leads to a simple and tractable approach but univariate modeling does not provide a complete statistical description of the textured images.

Latterly, researchers started to study the joint statistics of the wavelet coefficients of both textured and natural images. Some models and methods were formulated to explore statistical dependencies existing across scale, orientation and position. Portilla and Simoncelli presented a statistical model based on joint statistics of steerable pyramid coefficients [START_REF] Portilla | A Parametric Texture Model Based on Joint Statistics of Complex Wavelet Coefficients[END_REF]. In their work, efficient algorithm of texture synthesis was developed giving increased synthesis quality. However this model is not tractable for classification applications due to the largeness of the signature. Tzagkarakis et al. [START_REF] Tzagkarakis | Rotation-invariant texture retrieval with gaussianized steerable pyramids[END_REF] proposed a computationally complex Gaussianization procedure of the filter banks output in order to model wavelet coefficients with a multivariate normal distribution. Powerful statistical algorithms have been developed for image denoising using a Multivariate Generalized Gaussian Distribution (MGGD) [START_REF] Cho | Multivariate statistical modeling for image denoising using wavelet transforms[END_REF] and Elliptically Contoured Distribution (ECD) [START_REF] Tan | Multivariate Statistical Models for Image Denoising in the Wavelet Domain[END_REF], but no closed expression exists for the KLD between these joint distributions to measure similarity in a retrieval or classification context. Recently, Boubchir et al. [START_REF] Boubchir | Multivariate Statistical Modeling of Images in Wavelet and Curvelet Domain using the Bessel K Form Densities[END_REF] introduced Multivariate Bessel K Forms density (MBKF) to characterize joint statistics in wavelet domain. Oppositely to MGGD and ECD, the MBKF has a closed-form for the KLD, so the MBKF model seems a good candidate for modeling texture for retrieval or classification issues.

In this paper, a joint probability distribution of wavelet coefficients based on SIRV model is proposed. A SIRV process is a non-homogeneous Gaussian model defined by the product between the independent zero mean Gaussian vector and the root of a positive random variable. We provide a closed form solution to the KLD to measure similarity.

The remainder of this paper is as follows. We review in section 2 the statistical retrieval framework. In section 3 SIRV based multivariate model is introduced and the related KLD is calculated. Finally, in section 4 experimental results are given to evaluate retrieval performance.

PROBABILISTIC IMAGE RETRIEVAL

We establish the formal framework of probabilistic image retrieval. Consider an image database with M images

M i I i ≤ ≤ 1 ,
. Each The retrieval task is to search the N most similar images to a given query image q I . It is natural to select the most similar image r I to q I as the one whose parameter s θ leads to a maximization of the log-likelihood function,

i.e.

( ) (

)

∑ = = n j i ij i x p n s 1 ; log 1 max arg θ (1)
Using the weak law of large number, we have

( ) ( ) [ ] i X i n x p s q θ θ ; log max arg Ε = ∞ → (2) ( ) ( ) ( ) ∫ ∫ = x d x p x p i X q X i θ θ ; log ; max arg (3)
where the term q θ Ε denotes the expectation with respect to ) ; ( q X x p θ . Equation ( 3) can be rewritten as the following minimizing problem [START_REF] Portilla | A Parametric Texture Model Based on Joint Statistics of Complex Wavelet Coefficients[END_REF] This can be seen as equivalent to minimizing the KLD between

( ) ( ) ( ) x d x p x p s i X q X i ∫ ∫ - = θ θ ; log ; min arg (4) ( ) ( ) ( ) x d x p x p x p i X q X q X i ∫ ∫         = θ θ θ ; ; log ; min arg
) ; ( q X x p θ and ) ; ( i X x p θ noted ) ( i q p p KLD
. So, to select the N top matches to the query image q I we retrieve

the set of images { } N k k k I I I , , , 2 1
such as:

) ( ) ( ) ( 2 1 N k q k q k q p p KLD p p KLD p p KLD ≤ ≤ ≤ (6)

SIRV BASED MODELING

Spherically Invariant Random Vectors (SIRVs) has been appropriately used in modeling non-Gaussian problems. This is for instance, the case for radar clutter returns [START_REF] Conte | Recursive estimation of the covariance matrix of a compound-Gaussian process and its application to adaptive CFAR detection[END_REF], radio fading analysis [START_REF] Yao | Unified theory on wireless communication fading statistics based on SIRP[END_REF], or sonar interferences [START_REF] Barnard | Statistical normalization of spherically invariant non-Gaussian clutter[END_REF]. The joint statistics of wavelet coefficients also exhibit the obvious non-Gaussianity and SIRV model is suitable to characterize these statistics.

Let us consider an image decomposed into oriented subbands at multiple scales. We denote as . We denote as x a neighborhood of coefficients clustered around this reference coefficient. We assume the coefficients within each local neighborhood around a reference coefficient of a subband are characterized by a SIRV model. Formally, a random vector x is a Spherically Invariant Random Vector [START_REF] Yao | A representation theorem and its applications to spherically-invariant random processes[END_REF] if it is the product of the square root of positive random variable τ called the texture and a d -dimensional independent zero-mean Gaussian vector g with covariance

) ( g g M t Ε = verifying d M Tr = ) ( : g x × = τ (7)
The joint density of vector x is determined by the covariance matrix M and the mixing density

) (τ τ p : ( ) ∫ +∞ -         - = 0 1 2 1 2 2 exp 2 ) ( ) ( τ τ τ π τ τ d x M x M p x p t d X (8)
To complete the model, we need to specify the probability density

) (τ τ p
. We propose Weibull distribution as an appropriate description of the texture τ , given by:

0 , exp ) , ; ( 1 τ τ τ τ τ               -       = - a a b b b a b a p (9)
where 0 a is the shape parameter and 0 b is the scale parameter. By inserting Eq. ( 9) to Eq. ( 8), the joint density which models the vector of wavelet neighbors is: [START_REF] Conte | Recursive estimation of the covariance matrix of a compound-Gaussian process and its application to adaptive CFAR detection[END_REF] In this case, the hyperparameters of the correspondent joint distribution

( ) ∫ +∞ - -                       + - = 0 1 2 1 2 1 2 exp 2 ) ( τ τ τ τ π τ d b x M x M b a x p a t d a a X
) (x p X are noted ) , , ( b a M
.

To estimate these parameters we will first use a fixed point (FP) estimate for M . We note

N i x i , , 1 , =
the realizations of the d -dimensional vector x (in our case, d is the size of the neighborhood). In their recent work [START_REF] Pascal | Covariance Structure Maximum-Likelihood Estimates in Compound Gaussian Noise: Existence and Algorithm Analysis[END_REF], Pascal et al provide the proof of existence and uniqueness of M ˆas the solution of the following: ( )

M f M = (11) where f is given by ( ) ∑ = - = N i i t i t i i x M x x x N d M f 1 1 ˆ (12)
Equation ( 11) is solved using an iterative procedure with the initial guess from the Maximum Likelihood (ML) estimate

∑ = = N i t i i x x N M 1 0 1 ˆ
Experiments show that typically only around five iteration steps are required to obtain convergence. We note that FP estimation of the covariance matrix M does not depend on the texture τ but only on the vectors i x . Then, the ML estimate of texture τ is given by: Once we estimate

N i d x M x i t i i , , 1 , ˆ1 = = - τ (13)
N i i , , 1 , ˆ = τ
, the Weibull parameters ) , ( b a can be estimated using ML which we can found in [START_REF] Krishnamoorthy | Handbook of Statistical Distributions with Applications[END_REF]:

( ) ( ) a N i a i N i i a i N i i a i N b a 1 1 1 1 ˆ1 ; ln ln 1         = - = ∑ ∑ ∑ = = = τ τ τ τ τ (14)
Using SIRV modeling, joint distribution of wavelet coefficient is represented by a univariate Weibull distribution and a multivariate Gaussian distribution. Taking a neighborhood of dimension equal to 2, this is illustrated in Fig. 1 (b) that shows the good fitting of Weibull density with the normalized histogram of the multiplier τ estimated from a detail subband obtained by steerable pyramid decomposition; the empirical density of the correspondent bivariate Gaussian vectors g is presented in Fig. 1 (c).

In the best of our knowledge, the proposed joint distribution in Eq. ( 10) doesn't have a closed analytical form. However, we can derive a closed form solution for the KLD thanks to SIRV representation. Consider two joint distributions ) , ,; (

1 1 1 1 b a M x f and ) , , ; ( 2 2 2 2 b a M x f
. In a SIRV representation, the texture τ and the Gaussian vector g are independent, so the KLD between the two joint distributions is the sum of the KLD between the KLD of the two Weibull distributions and the KLD between the two multivariate Gaussian densities:

) ) ; ( ) ; ( ( ) ) , ; ( ) , ; ( ( ) ) , , ; ( ) , , ; ( ( 2
2 1 1 2 2 2 1 1 1 2 2 2 1 1 1 1 1 M g p M g p KLD b a p b a p KLD b a M x f b a M x f KLD Gaussian Weibull + = τ τ (15)
In other hand ) and (.) Γ is the gamma function.

1 ) )( ln( ) ln( ) ln( ) )( 1 ( ) ) , ; ( ) , ; ( ( 1 2 2 1 1 2 2 1 1 2 1 1 2 2 2 2 1 1 1 2 1 2 - - + - + - + + Γ = - - γ γ τ τ
and the KLD for the d -dimension Gaussian case is:

        -         + = - d M M M M tr M g p M g p KLD Gaussian 1 2 1 1 2 2 2 1 1 ln ) ( 5 . 0 ) ) ; ( ) ; ( (

EXPERIMENTAL RESULTS

The experiments give an evaluation of the proposed model in the framework of texture retrieval. We use the same experimental setup presented in [START_REF] Do | Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF] and [START_REF] Kwitt | Image similarity measurement by Kullback-Leibler divergences between complex wavelet subband statistics for texture retrieval[END_REF]. We work with 40 texture classes from VisTex database [START_REF]MIT Vision and modeling group[END_REF]. From each of these texture images of size 640x640 pixels, 16 subimages of 160x160 are created. A test database of 640 texture images is thus obtained. A query image is any one of these images in the database. The relevant images for each query are the other 15 images obtained from the same original 640x640 image.

We employ the steerable pyramid decomposition proposed in [START_REF] Portilla | A Parametric Texture Model Based on Joint Statistics of Complex Wavelet Coefficients[END_REF] ( we note Nsc the number of scales and Nor the number of orientations). However, orthonormal or biorthogonal wavelet representations can be used. We use the conventional criterion of precision/ recall to compare the performance of the proposed SIRV model with the retrieval approach using GGD presented in [START_REF] Do | Wavelet-based texture retrieval using generalized Gaussian density and Kullback-Leibler distance[END_REF] and with the MBKF distribution [START_REF] Boubchir | Multivariate Statistical Modeling of Images in Wavelet and Curvelet Domain using the Bessel K Form Densities[END_REF]. Ψ is the digamma function.

number
The recall/precision curves presented in Fig. 2 show the improvement obtained by using SIRV model. The results of average retrieval rates according to different scales and orientations are summarized in Table 1. We can see that SIRV modeling significantly improves recognition rates, e.g from 73% to 79%, compared with GGD and MBKF. Furthermore, we observe also that the proposed method converge faster than the two others. For example we retrieve 90% of the relevant images with a query of size 45 when we must consider a query of size 70 to retrieve the same percentage if we employ the two others models.

CONCLUSION

In this work we have shown that image retrieval improve considerably when wavelet coefficients are jointly modeled using SIRV model. In a statistical retrieval framework, we have proposed a closed form for the Kullback-Leibler Divergence as a similarity measure. The model is validated in the retrieval system and achieves better recognition rates compared to GGD and MBKF distributions.

  scale s , orientation o and centered at spatial

1 .

 1 multiplier τ histogram fitted with Weibull density (c) empirical joint density of the related bivariate Gaussian vector g Fig. Example texture and SIRV representation in the case of bivariate modeling for the subband 11 B resulting from steerable pyramid decomposition ( so B denotes the subband at scale s and orientation o )

Fig. 3 .

 3 Fig.3. models convergence comparisonWe use KLD as image similarity for the tree compared model. For MBKF distributions, developing the solution leads to:

Table 1 :

 1 Average retrieval rate (%) comparison

	recall	=		of relevant retrieved images
				number of relevant images
	precision	=	number of relevant retrieved images
				number of retrieved images
	For multivariate modeling, the neighborhood may
	include wavelet coefficients from other subbands (i.e.
	corresponding to nearby scales and orientations) as well as
	from the same subband. In our experiments, we used a
	neighborhood drawn from the same subband.


Ahmed Drissi

El Maliani

Nouredine Lasmar
email: nour-eddine.lasmar@ims-bordeaux.fr

Mohammed El Hassouni
email: mohamed.elhassouni@gmail.com

Yannick Berthoumieu
email: yannick.berthoumieu@ims-bordeaux.fr

Texture Classification Based On The Generalized Gamma Distribution And The Dual Tree Complex Wavelet Transform

This paper deals with stochastic texture modeling for classification issue. A generic stochastic model based on three-parameter Generalized Gamma (GG) distribution function is proposed. The GG modeling offers more flexibility parameterization than other kinds of heavy-tailed density devoted to wavelet empirical histograms characterization. Moreover, Kullback-leibler divergence is chosen as similarity measure between textures. Experiments carried out on Vistex texture database show that the proposed approach achieves good classification rates.

I. INTRODUCTION

Coming up with accurate texture representation is central in computer vision applications such as covering texture retrieval, texture synthesis and texture classification. Considering the literature, recent works point out that the marginal stochastic modeling approach is well suited for texture retrieval. It exhibits demonstrated effectiveness properties and its formulation simplicity contributes to its popularity for characterizing empirical wavelet histograms. The Generalized Gaussian density (GGD), proposed in [START_REF] Do | Wavelet-Based Texture Retrieval Using Generalized Gaussian Density and Kullback-Leibler Distance[END_REF], has a considerable notoriety. Do et al. in their pioneering work proposed the GGD model for characterizing histograms of subbands issued from the Discrete Wavelet Transform (DWT). Recently, other kinds of decomposition appeared such as Dual Tree Complex Wavelet Transform (DTCWT) [START_REF] Kingsbury | The Dual-Tree Complex Wavelet Transform: A new Technique for Shift-Invariance and Directional Filters[END_REF]. This latter decomposition has many advantages in terms of shift invariance, rotation invariance and directional selectivity offering six orientations in each decomposition level. In the case of complex coefficients characterization, Weibull and Gamma distributions showed a good ability to modeling complex wavelet coefficients [3] [4]. The combination of these models with the Kullbak-Leibler divergence outperforms the GGD based approach in the context of texture retrieval. In this paper, our aim is to propose a more generic model than the previous ones devoted to the complex coefficient magnitudes distribution modeling. We propose a three-parameter Generalized Gamma (GG) distribution, as a flexible model including many distributions, especially, the Weibull and Gamma distributions. The interpretation and use of the GG modeling are discussed, and we present a texture classifier using both the GGD model and the GG one. To measure the similarity between training and testing textures, in the classification framework, we propose the Kullback-Leibler divergence, which measures the discrepancy according to the model. The outline of this paper is as follows, in the next section, we give a definition to the Generalized Gamma density and associated Kullback-Leibler divergence.

In the third section we present the texture classification framework, while experimental results, are presented in section 4, which is followed by conclusions.

II. GENRALIZED GAMMA DENSITY FOR TEXTURE MODELING

As we have already mentioned, the Generalized Gamma distribution is a generic model, which is suitable for duration analysis, and thus, very useful in different fields [START_REF] Stacy | A Generalization of the Gamma Distribution[END_REF]. The probability density function of the Generalized Gamma distribution is:

f (y; α, τ, λ) = τ λ ατ Γ(α) y ατ -1 exp -( y λ ) τ , y ≥ 0, α, τ, λ > 0 (1)
Where α and τ are shape parameters, λ is the scale parameter, and Γ(.) is the gamma function. The Generalized Gamma family encompasses a variety of well-known models as limiting distributions. It provides the Gamma density for τ = 1, Weibull for α = 1, and Exponential for α = τ = 1.

We use the maximum likelihood method for estimating the GG parameters, supposing y = (y 1 , y 2 , ..., y M ), a set of M independent coefficients. Let us define the maximum likelihood function of the sample

y : L(y; α, τ, λ) = log M i=1 f (y; α, τ, λ) ∂L(y; α, τ, λ) ∂α = -M (τ log λ -ψ(α)) + M i=1 τ log y i = 0. ( 2
)
Where

ψ(z) = Γ ′ (z)/Γ(z), is the digamma function. ∂L(y; α, τ, λ) ∂τ = M ( 1 τ -α log λ)+ M i=1 α log y i -( y i λ ) τ log y i λ = 0. ( 3
)
∂L(y; α, τ, λ) ∂λ = - M ατ λ + τ λ -τ λ M i=1 y i = 0. (4)
Thus, the parameters are deduced by solving a system of three equations:

λ = [ 1 M α M i=1 y i τ ] 1 τ . ( 5
) α = 1 τ [ M i=1 y i τ log y i M i=1 y i τ -log y i ] -1 . ( 6
) log M α( M i=1 y i ) τ M M i=1 y i τ -ψ( α) = 0. (7)
To make up for the high non-linearity of the system, we use the numerical algorithm of Cohen et al. [START_REF] Cohen | Parameter Estimation in Reliability and Life Span Models[END_REF]. Figure 1, shows the fitting of our GG pdf to DTCWT subband coefficient histogram, compared to Weibull and Gamma fitting. The Similarity Measurement is a crucial step in the classification process. Thus, we need an efficient distance between the images of our database. We adopt the Kullback Leibler divergence which is fully defined by the model parameters of the images. The Kullback Leibler divergence between two PDFs p(X; θ q ) and p(X; θ i ) is:

KLD(p(X; θ q ) p(X; θ i )) = p(X; θ q ) log p(X; θ q ) p(X; θ i ) dx (8)
By inserting Eq. ( 1) into Eq. ( 8), we derive the KLD between two GG PDFs:

KLD(f 1 (.; α 1 , τ 1 , λ 1 ), f 2 (.; α 2 , τ 2 , λ 2 )) = log( τ 1 λ α2τ2 2 Γ(α 2 ) τ 2 λ α2τ2 1 Γ(α 1 ) )+ +( λ 1 λ 2 ) Γ(α 1 + τ2 τ1 ) Γ(α 1 ) + ψ(α 1 ) τ 1 (α 1 τ 1 -α 2 τ 2 ) -α 1 (9)
So, two GG PDFs are similar to each other when the KLD between them is small. The mean idea here is to calculate the complex transform of signals using two filter's trees A and B (two separate DWT decompositions), that independently generate the real and imaginary coefficients (Figure 2). Contrary to the classic DWT, DTCWT don't suffer of the problems of shift invariance and directional selectivity, providing six orientations at each decomposition level.

A. Feature extraction

Each texture of the database is decomposed via the Dual Tree Complex Wavelet Transform introduced by N.Kingsbury [START_REF] Kingsbury | The Dual-Tree Complex Wavelet Transform: A new Technique for Shift-Invariance and Directional Filters[END_REF]. We choose two scale decomposition Nsc=2, and six orientations Nor=6 in each scale. So we may have twelve subbands for each texture. The histograms of the coefficients magnitudes are then modeled by the Generalized Gamma distribution, and the signature (feature vector) of a texture is composed by concatenation of the estimated parameters: v GG = {α 1 , τ 1 , λ 1 , ..., α 12 , τ 12 , λ 12 }

B. Classifier

We adopt the K-nearest-neighbor (KNN) classifier, as a simple and good reference in texture classification [START_REF] Fukunaga | Introduction to Statistical Pattern Recognition[END_REF]. KNN is a type of Instance-based classifiers. The mean idea of this approach is that an unknown instance can be done by relating the unknown according to some similarity measure. In our case the Kullback-Leibler divergence is chosen as a similarity measure. The training samples are the feature vectors obtained after the feature extraction step, associated to a class label vectors. In the classification (recognition) phase, an unknown sample is classified by assigning the label of the class to which the majority of its K nearest neighbor belongs; it's a sort of majority voting.

V. EXPERIMENTAL RESULTS

This section aims to show the effectiveness of the Generalized Gamma and DTCWT based approach in the framework of texture classification. We conduct our experiments using the setup presented in [START_REF] Van De Wouver | Statistical texture characterization from discrete wavelet representation[END_REF], considering 30 textures of size 512 × 512 from the Vistex database [START_REF]MIT vision and modeling group[END_REF], as it's shown in the Figure 3. Each image, is splited into sixteen Sub-images of size 128×128, so then we obtain a database of 480 samples. Before anything else, we construct a training and testing sets for the classification process. We will, randomly, choose Q samples from each of the 30 texture classes, this will be considered as the training set. Testing set is the rest of the samples. The classification accuracy is obtained by varying Q from 2 to 8, repeating the KNN test 100 times and then returning the average performance. Firstly, let's point out the flexibility of the Generalized Gamma model in relation to his limiting cases, in the context of texture classification. Table 1 shows the average classification rate vs. the number of training samples for Generalized Gamma, Weibull and Gamma models. We can see the superiority of Generalized Gamma model's classification rates, which is due to its high flexibility.

Secondly, to demonstrate the efficiency of our approach, we compare it to the one based on Generalized Gaussian Density (GGD). This approach was proposed by Woower et al. [START_REF] Van De Wouver | Statistical texture characterization from discrete wavelet representation[END_REF], they used the same KNN classifier using an L2-norm distance as similarity measure. Thirty texture images: From left-right and top bottom: Bark.0009, Bark.0004, Bark.0006, Bark.0008, Bark.0009, Brick.0001, Brick.0004, Brick.0005, Fabric.0000, Fabric.0004, Fabric.0007, Fabric.0009, Fabric.0011, Fabric.0013, Fabric.0016, Fabric.0017, Fabric.001, Food.0000, Food.0002, Food.0005, Food.0008, Grass.0001, Sand.0000, Stone.0004, Tile.0001, Tile.0003, Tile.0007, Water.0006, Wood.0001, Wood.0002 The results, of the 100 random experiments (Table 2), show that our approach achieves the highest average classification rate 96.33% for Q=8 (nb. of training samples from each class), while the GGD-based approach obtain 92.41% . Finally, we will focus on the suitability of the DT-CWT decomposition with the Generalized Gamma based approach. For this purpose, we shall compare it to the steerable pyramid decomposition [GG+steerpyr] performances [START_REF] Simoncelli | The Steerable Pyramid: A Flexible Architecture for Multi-Scale Derivative Computation[END_REF]. Table 3 points out the choice of the decomposition; we compare the same model (GG) in combination with the two decompositions. We observe a higher performance of the DTCWT based approach, in term of higher average classification rate.

VI. CONCLUSION

This paper presents, a statistical approach for texture classification, using the Generalized Gamma Density (GG), and the Dual Tree Complex Wavelet Transform decomposition. The Generalized Gamma Density has the advantage to be a more flexible and a generic model. Each, texture of the database, was decomposed via the DTCWT decomposition, and the histograms of the magnitudes of the subbands coefficients are then modeled using the GG density, which the parameters serve as a features for the texture classification process. A KNN classifier is designed for classifying the unknown textures. The experiments, on a thirty class of the VisTex database, show that the combination of the Generalized Gamma Density and the DTCWT gives good performances over other models and decompositions. In future work, we study the performance of the Generalized Gamma Density in combination with the Dual Tree M-band Wavelet Transform.

Figure 1 .

 1 Figure 1. Wavelet subband coefficient magnitude histogram fitted with Generalized Gamma, Weibull and Gamma densities. Example for the Bark.0000.05 subimage. The estimated parameters are: α=3.1, τ =1.46 and λ=0.77 for GG, α=7.45 and β=1.44 for Weibull, α=1.95 and β=3.44 for Gamma

Figure 2 .

 2 Figure 2. Dual trees for complex wavelet transform

Figure 3 .

 3 Figure 3.Thirty texture images: From left-right and top bottom: Bark.0009, Bark.0004, Bark.0006, Bark.0008, Bark.0009, Brick.0001, Brick.0004, Brick.0005, Fabric.0000, Fabric.0004, Fabric.0007, Fabric.0009, Fabric.0011, Fabric.0013, Fabric.0016, Fabric.0017, Fabric.001, Food.0000, Food.0002, Food.0005, Food.0008, Grass.0001, Sand.0000, Stone.0004, Tile.0001, Tile.0003, Tile.0007, Water.0006, Wood.0001, Wood.0002

Table I AVERAGE

 I CLASSIFICATION RATE VS. NUMBER OF TRAINING SAMPLES FOR GENERALIZED GAMMA, WEIBULL AND GAMMA MODELS

	Q	2	3	4	5	6	7	8
	GG	85.63 88.65 91.15 92.43 93.34 94.85 96.33
	Weibull 84.88 87.29 89.56 91.25 92.44 93.16 95.08
	Gamma 84.19 86.73 88.96 91.33 91.84 92.76 94.42
				Table II
	AVERAGE CLASSIFICATION RATE VS. NUMBER OF TRAINING SAMPLES
		FOR GENERALIZED GAMMA AND GGD
	Q	2	3	4	5	6	7	8
	GG 85.63 88.65 91.15 92.43 93.34 94.85 96.33
	GGD 83.9 85.37 88.31 89.26 91.12 91.96 92.41

This work was supported by franco-morrocan CNRS-CNRST project STIC 05/10