O Ther

C Colin

M H Berger

L Gerbaud

A Dourfaye

Reactive Imbibition of WC-Co Substrate for PDC Cutters Used in Oil and Gas and Mining Drilling

Keywords: Cemented Carbides, Functionally Graded Materials (FGM), Imbibition

Cemented carbides are used in rock drilling for mining tools and wear resistant parts. These composite materials possess an excellent compromise between hardness and toughness. Nowadays, the concept of graded structure is widely studied to improve these two properties simultaneously, and so to increase the service life of drilling tools.

A continuous composition gradient on several millimetres is generated in commercial WC-Co substrate for PDC cutters by using Reactive Imbibition method. The effects of this process are analysed in terms of microhardness, cobalt concentration and WC grain size. A continuous gradient of about 300HV on 8mm-height substrate is obtained in one-step by imbibition process into combination with a boron-rich coating deposed on its free surface. In part, this gradient of hardness and its shape are preserved after HPHT (high pressure-high temperature) process that is used for the diamond table deposition on the WC-Co substrate. Such gradient must significantly increase the cutter service life.

Introduction

The goal of this study is to develop a new process to enhance the service life of drilling tools cutters. Indeed, drilling conditions in deep and heterogeneous wells lead to premature wear and fracture of tools because of abrasive rocks and repeated shocks. Drilling tools are made of cemented carbides which are composed of WC hard phase embedded in a Cobased ductile metallic matrix.

In order to increase hardness and toughness, which are antagonistic properties, Functionally Graded Material (FGM) way is chosen. Gradation of cemented carbides is not a recent idea, its interest has been demonstrated in 1972 [START_REF] Mb | Gradients in composite materials[END_REF]. Since, a lot of processes have been developed to grade WC-Co: Infiltration [START_REF] Kieffer | Uber die Herstellung von Hartmetallen nach dem Tränkverfahren[END_REF], Dual Properties (DP) carbides [3], [4], [START_REF] Eso | Kinetics of cobalt gradient formation during the liquid phase sintering of functionally graded WC-Co[END_REF], Carbide Free Layer (CFL) [START_REF] Gustafson | Binder-phase enrichment by dissolution of cubic carbides[END_REF], [START_REF] Chen | Fundamentals of liquid phase sintering for modern cermets and functionally graded cemented carbonitrides (FGCC)[END_REF], [START_REF] Zackrisson | Microstructure and performance of a cermet material heat-treated in nitrogen[END_REF], [START_REF] Lengauer | Tailoring hardness and toughness gradients in functional gradient hardmetals (FGHMs)[END_REF], Ion implantation [START_REF] Ie | The life of WC-Co cutting tools by plasma immersion ion implantation[END_REF], Cr3C2 coating on cemented carbides before sintering [START_REF] Ms | Hard composite and method of making the same[END_REF] and multilayers or stacking of different cemented carbides by their cobalt content and carbide grain size [START_REF] Colin | Processing of functional-gradient WC-Co cermet by powder metallurgy[END_REF], [START_REF] Colin | Equilibrium distribution of liquid during sintering of assemblies of WC/Co cermets[END_REF]. These processes rely on different techniques, which have two common drawbacks:

the small distance on which is developed the cobalt redistribution, -the discontinuity of gradients of the hardness and cobalt content. But, there is another gradation process which overcomes these drawbacks: Imbibition from dense cemented carbides allows to continuously enrich the substrates in binder phase (between 10 and 20wt.% Co) over several millimetres (about 10 to 20mm). First results on this process were published in 1972 by Lisovsky and Babich [START_REF] Af | Redistribution of molten cobalt in powder metallurgical WC-Co hard alloys[END_REF].

Recently, Reactive Imbibition [START_REF] Sorlier | Développement d'un procédé de graduation des carbures cémentés WC-Co basé sur l'imbibition, amélioration de la durée de vie des taillants de forage[END_REF], patented by Varel Europe/Armines is a one step process combining imbibition from dense cemented carbides and a boron-rich coating on their free surfaces. This new process generates a continuous variation of cobalt through several millimetres and with an excellent gradient shape control to increase at the same time the surface hardness and core toughness.

In this article, this new process has been chosen to grade the substrate of commercial cutters. The underlying principles governing the cobalt phase migration and distribution are examined and based on thermodynamic properties and phase reactions during imbibition and HPHT treatments with the aim of gaining an understanding of the phenomena.

Reactive Imbibition process

A.

Imbibition treatment Imbibition consists in a migration of liquid arising from a reservoir (filler material) in a two-phase system (solid / liquid) which is raised to high temperature. This phenomenon is different from well-known diffusion mechanism, because the liquid migration has the same composition than liquid within sample. The imbibition phenomenon is also different from infiltration mechanism which implies a three-phase system (solid / liquid / gas, i.e. a porous media). According to Lisovsky [START_REF] Af | Redistribution of molten cobalt in powder metallurgical WC-Co hard alloys[END_REF], three conditions are required in a two-phase system to make imbibition possible: a good wetting of solid grains by the liquid matrix, a partial dissolution of solid grains in this liquid matrix and the absence of equilibrium dihedral angle defined by γSS/2γSL>1. So, if the ratio of the solid-solid to solid-liquid surface energy is high (>2), then the penetration of grain boundaries by liquid is possible because it is energetically more favourable to system evolution. Also, in the process of imbibition, the liquid phase in contact with a fully dense cemented carbide flows through the liquid channels formed by refractory particles (WC grains). A complete reconstruction of the refractory skeleton by particle rearrangement occurs during migration and Park [START_REF] Hh | Effect of dihedral angle on the morphology of grains in a matrix phase[END_REF] deals with this phenomenon by determination of the minimum interface energy configurations of solid / liquid system with varying dihedral angles and liquid contents. A similar approach is undertaken by Delannay and al [START_REF] Delannay | Equilibrium distribution of liquid during liquid sintering of composition gradient materials[END_REF] that expresses the variation of driving force for the absorption of liquid by solid/liquid system as a function of the liquid volume fraction, U, and the coordination number, nc. On the other hand, the Ostwald ripening with a modification of grain shape, without necessary grain growth is not a condition essential to liquid migration in imbibition process, according to Sorlier [START_REF] Sorlier | Développement d'un procédé de graduation des carbures cémentés WC-Co basé sur l'imbibition, amélioration de la durée de vie des taillants de forage[END_REF]. However, the driving force of liquid migration for the assemblies of fully dense cemented carbides is the reduction of the interfacial energy by the grain shape accommodation, according to Colin [START_REF] Colin | Migration d'un liquide eutectique dans des carbures cémentés WC-Co complètement denses de façon à réaliser des matériaux à gradient de composition, in: Les traitements des poudres et leurs conséquences[END_REF], [START_REF] Colin | Equilibrium distribution of liquid during sintering of assemblies of WC/Co cermets[END_REF]. The author puts forward the dwell time very long (several hours) used for these assemblies which allows the grain shape accommodation.

Usually, WC particles have an anisotropic solid / liquid interfacial energy. The system reaches its minimum energy for a certain amount of binder phase which allows adopting its equilibrium configuration, corresponding to lower energy level. Until the critical amount of binder phase has not been reached, a migration pressure exists in the system. The phenomenologic expression of the liquid migration pressure Pmig (MPa), established by Lisovsky is given by expression corrected by Sorlier [START_REF] Sorlier | Développement d'un procédé de graduation des carbures cémentés WC-Co basé sur l'imbibition, amélioration de la durée de vie des taillants de forage[END_REF]:

where U is the volume fraction of liquid and dWC (µm) is the mean WC grain size. It is thus predicted that the liquid migration pressure increases when the WC grain size and the Co content are smaller. This expression also points out that this migration pressure reaches zero for a critical volume fraction of liquid (Um=0.61) which corresponds to 32 wt.% Co. Thus, a WC-Co containing more than 30 wt.% Co cannot be enriched in binder phase.

B.

Reactive coating mechanism A boron nitride coating is uniformly deposited on WC-Co part by spraying (Fig. 1a). Just before the migration of liquid from a reservoir composed of green part with eutectic composition, the BN coating is destabilized by liquid phase, after the solidus of the cemented carbide is reached. According the mechanism proposed by Sorlier [START_REF] Sorlier | Développement d'un procédé de graduation des carbures cémentés WC-Co basé sur l'imbibition, amélioration de la durée de vie des taillants de forage[END_REF], this destabilization leads to a diffusion of boron inside the substrate through its liquid phase (Fig. 1b). The solubility of boron into the liquid phase is small and so, its saturation is quickly reached. Also, borides, in particular WCoB, precipitate under the substrate surface. The formation of this ternary compound involves a decrease of W and Co contents into the liquid phase. So, the carbon in excess in the liquid can diffuse towards the bulk (Fig. 1c), causing the diffusion of the cobalt in the same direction, as for the DP carbides. Indeed, it is well-known that when there is a carbon gradient into the liquid, cobalt appears to migrate in the direction of carbon diffusion, resulting in a redistribution of cobalt in the cemented carbide [START_REF] Zz | Liquid phase sintering of functionally graded WC-Co composites[END_REF]. However, the microstructure shows no free-carbon nor η-Co3W3C phase. Materials and sample preparation The WC-Co substrates used in this study are commercial ones. The cylinder-shape substrates for cutters have a 14.40mm diameter and a height ranging from 7.35 to 8.10mm. Their initial Co content is 12.25 wt.%, their grain size measured by linear intercept ranges from 0.5 to 4µm and they are fully dense (Porosity< 0.02% according to ASTM B276). The hardness of these substrates is measured under a 2kg load and reaches 1245 +/-30 HV.

The green part composed of the eutectic composition (WC-65wt.% Co) is prepared with powders supplied by Ceratizit. Once the powders have been blended in a mixer for 24h, 10mm diameter compacts are cold pressed at 200MPa by uniaxial compaction. During this time, the commercial substrates are cleaned in an ultrasonic alcohol bath for 5min. Coating is applied to those which need it.

Once the green part and the substrate have been prepared, the BN-coated substrate is superimposed on the green part before being put in a graphite furnace. The sample is processed at 1450°C under an Ar-5vol.% H 2 atmosphere in order to approximate industrial conditions. Then, the processed substrate has undergone a classic High Pressure / High Temperature (HPHT) step in order to assembly the diamond table on the WC-Co substrate. This treatment consisting in a quick heating of the substrate with diamond powder in a refractory metallic cap is carried out for a few minutes under a stable temperature / pressure condition adapted to diamond. In order to clarify the running of the experiment, Fig. 2 presents a schematic representation of the major steps of the process.

After the reactive imbibition treatment or the HPHT process, depending on the sample, the cutters are electric-discharge machined in order to cut them into two equal parts, ready to be polished and analysed. The samples are then ground and polished to a 1µm finish for microstructural examinations.

The diamond tables are processed under two Pressure/Temperature combinations, with the second combination at higher Pressure/Temperature. Two grain size distributions are used for the diamond (coarse and fine). The table 1 summarizes processed samples.

X B X C X Coarse 1 D X Coarse 1 E X Fine 1 F X Fine 1 G X Coarse 2 H X Coarse 2 I X Fine 2 J X Fine 2 B.
Characterization of cutters To characterize the microstructure, a Zeiss-DSM982 Scanning Electron Microscope (SEM) in secondary electron (SE) and back-scattered electron (BSE) modes under a 10 kV voltage is used. After microstructure binarization, linear intercept analysis is performed on nine SEM micrographie at x5000 magnification for each position. These nine micrographie are performed around hardness indentations in order to control the WC grain size along the revolution axis. This grain size measurement allows to verify that the two successive thermal treatments (imbibition and HPHT) have not generated normal and abnormal grain growth (Fig. 3.a).

In order to quantify the liquid migration in the sample, three experimental techniques based upon different principles have been compared. This first one is founded on the hardness measurement, since it is related to Co content for a constant grain size. A Buehler microhardness tester with a Vickers diamond tip is used to determine the Vickers hardness with a load of 2kg and an indentation time of 10s, according to NF A 03-154. Indentations are performed along the revolution axis and on the right part of the sample, because the Co migration is symmetrical across this axis. Three hardness measurements are performed for each pre-determined position. The second method to measure the Co content is based upon image analysis. For an isotropic microstructure, surface fraction is equivalent to volume fraction of each phase in a two-phase structure. Then, the volume fraction of cobalt (VCo) is converted into mass content (XCo) from the relation as follows: Fig. 2. Major steps of the process where ρCo and ρWC are the density of cobalt phase and tungsten carbide, respectively. In this case, the β-phase (solid solution with a large solubility of W) is assimilated to pure Co.

In the third method, Co content is evaluated using a microprobe (CAMECA-SX100) to make composition profiles along the symmetrical axis of the sample. Measurements are undertaken under 10kV voltage and 80nA intensity. The beam is defocused with a 20µm spot size. For each measurement, five spots around a hardness indentation are used. Contrary to image analysis method, it is directly obtained the Co content in weight percent by microprobe.

An experimental relation between hardness (HV) and Co mass content (XCo) obtained by image analysis is established from several samples (as-received, imbibed and HPHT processed samples).

where HV (kg/mm 2 ) is the Vickers Hardness (between 1000 and 1600HV) and XCo (wt.%) is Co mass content. These values of cobalt content deduced from microhardness are compared to values measured by image analysis and microprobe experiments. Fig. 3.b compares the results obtained with the three methods. These methods are well correlated, and so the microhardness is considered for Co gradient characterization. A.

Gradient before HPHT Microhardness measurements are used to quantify the gradients generated by the reactive imbibition process. The microhardness profiles (Fig. 4) are produced along the axis of symmetry (r=0) and the axis oblique (r=z). Fig. 5. show Co distribution for sample A and B along the axis of symmetry. The hardness maps on the right part of the samples are presented in Fig. 6. On the maps, each black spot represents a hardness test mark.

Sample A (as-received substrate) is an unprocessed sample. No significant gradient is present on the hardness maps. The hardness variation is mainly in the range from 1215 to1275HV with an average value around 1245HV.

Sample B is a processed substrate by Reactive Imbibition. On this sample, hardness varies from a minimum of 1080HV in the lower part to a maximum of 1350HV in the upper part. The decrease of hardness in the lower part is undoubtedly due to the migration of the Co-rich liquid phase from the green part. The hardness increase in the upper part is due to the borides formation. Gradient after HPHT: effect of substrate gradation Whatever the effect studied (substrate gradation, HPHT cycle and diamond grain size), the HPHT process induces a decrease of Co content in the cemented carbide with a hardening of the substrate.

Fig. 6 underlines that the samples C, E, G and I are harder than sample A, in the case of an ungraded substrate. For HPHT processed samples from as-received substrates, the difference of hardness varies from 1275 to 1485HV against 1215 to 1275HV for unprocessed sample by HPHT treatment. So, HPHT process clearly develops a gradation from homogeneous substrates with a amplitude that can reach about 200HV on 6mm-height. These hardness gradients, more particularly, the iso-Borides values of hardness are parallel to the lower surface.

On the same way, it is observed a similar phenomenon with the graded substrates. The samples D, F, H and J are harder than sample B. In this case, the maximum variation of hardness for HPHT processed samples from a graded substrate ranges from 1160 to 1485HV against 1080 to 1350HV for unprocessed sample by HPHT treatment. The HPHT treatment significantly increases the difference of hardness (∆HV=325 against ∆HV=270) but it always displays a dome-shaped gradient, generated by reactive coating enriched in boron. This dome-shaped gradient allows associating a hard shell with a softer bulk. The global or mean hardness of these new processed cutters (D, F, H and J) is lower than that of the commercial cutters (C, E, G and I).

C. Gradient after HPHT: effect of HPHT treatment at "low" and 'high" Pressure/Temperature It can be observed the effect of HPHT treatment at "low" and "high" Pressure/Temperature on Fig. 8 and Fig. 9, respectively. A slight decrease of the hardness appears between 3 and 5mm from the bottom of each ungraded substrate (C and E). This phenomenon doesn't occur for HPHT process at "high" Pressure/Temperature from the same substrate ungraded (G and I).

On graded substrates, the "low" Pressure/Temperature more preserves the dome-shaped gradient and its amplitude is lowered (D and F), compared to "high" Pressure/Temperature (H and J). This difference of hardness can reach about 150HV on 6mm for samples C/E and D/F against about 200HV to 300 HV on 6mm for samples G/I and H/J on the revolution axis. The same statement can be done on the r=z axis. So, for "High" HPHT treatment, the amplitude of gradient can be strongly increased (on an average of 50 to 150HV). All these results support the conclusion established earlier by several authors [START_REF] Miess | Rai Fracture toughness and thermal resistance of polycrystalline diamond compacts[END_REF], [START_REF] Cabral | Cobalt diffusion in different microstructured WC-Co substrates during diamond chemical vapor deposition[END_REF] that the binder phase of the WC-Co substrate migrates, more exactly infiltrates into the diamond compact during the HPHT process in order to act as catalyst for the sintering of diamond powder (neck formation between diamond grains).

D.

Gradient after HPHT: effect of diamond grain size distribution Fig. 7 and Fig. 8 summarize hardness profiles obtained after HPHT treatment for coarse and fine diamond grain size distributions, respectively. Hardness gradients seem to be very close and then independent of the mean diamond grain size, as shown in Fig. 9 and Fig. 10 These results are surprising and disagree with previous researches [START_REF] Uehara | High pressure sintering of diamond by cobalt infiltration[END_REF].

Axis r=0

Axis r=z Fig. 7. Hardness profiles with a coarse diamond grain size distribution

Axis r=0

Axis r=z Fig. 8. Hardness profiles with a fine diamond grain size distribution

Axis r=0

Axis r=z Fig. 9. Hardness profiles at "low" Pressure/Temperature for different diamond grain size distributions

Axis r=0

Axis r=z Fig. 10. Hardness profiles at "high" Pressure/Temperature for different diamond grain size distributions Discussion The goal of this article was to grade a commercial substrate for PDC cutters in order to improve their shock resistance without reducing their abrasion resistance. Nowadays, it is widely accepted that a component presenting a tough bulk and a hard surface shows superior mechanical properties. The most important gradient obtained in this study reaches 300 HV on 6mm-height, due to a variation of the binder phase into the graded substrate, after HPHT treatment at "high" Pressure/Temperature. Furthermore, the special shape of the gradient, obtained thanks to the use of the reactive coating, is truly innovative and interesting for mining applications. This particular shape is due to two phenomena before HPHT process: a competition between the liquid migration from the green part and the reaction of the BN coating with the cutter surface, leading to the formation of hardening phases and to a diffusion of the carbon and cobalt towards the bulk i.e. in the opposite direction to liquid migration by imbibition. To confirm the borides formation mechanism proposed by Sorlier [START_REF] Sorlier | Développement d'un procédé de graduation des carbures cémentés WC-Co basé sur l'imbibition, amélioration de la durée de vie des taillants de forage[END_REF], some TEM observations were performed in order to validate the presence of the WCoB phase (Fig. 11). This phase is a ternary orthorhombic compound known to enhance the wear resistance of carbide alloys in metal cutting due to its extremely high hardness (about 4300HV).The hardening of the cutter surface under the coating providing the dome shape is thus explained by the formation of a finely dispersed WCoB phase and by the decrease in the amount of the binder phase. The most interesting result is the conservation or even sometimes improvement of the gradient amplitude after HPHT treatment. This process at high temperature, involves the formation of a liquid phase which could give rise to a complete homogenization of the binder phase into the substrate. In fact, the gradients are modified but preserved after HPHT treatment. According to Sorlier [START_REF] Sorlier | Développement d'un procédé de graduation des carbures cémentés WC-Co basé sur l'imbibition, amélioration de la durée de vie des taillants de forage[END_REF], the liquid phase quickly infiltrates by capillarity into the polycrystalline diamond powder compact before sintering takes place. As long as porosity is concerned in powder compact, this infiltration generates a liquid flow under pressure in substrate similar to mass transport during creep tests at high temperatures in cemented carbides [START_REF] Favrot | Etude des différentes étapes de l'élaboration de structures graduelles en WC-Co : compactage à froid, frittage, homogénéisation[END_REF].

A surprising result is that the initial diamond grain size has no obvious influence on the redistribution of cobalt in the cemented carbides. However, it is widely known that the diamond grain size distribution would have a great influence on the porosity of powder compact before infiltration. Uehara [START_REF] Uehara | High pressure sintering of diamond by cobalt infiltration[END_REF] has studied the relation between the volume fraction of cobalt and the initial diamond grain size. The Co content increases in diamond table as the mean diamond grain size decreases. We can assume the different diamond grain size distributions ensure the same porosity in powder compact just before infiltration by liquid phase arising from substrate. This point would be further verified.

Another surprising result comes from the hardness decrease in the center of the ungraded substrate which occurs after a "low" HPHT process. According to the shape of the gradient formed, we can assume that the time at liquid state during "low" HPHT process is longer than "high" HPHT process in order to ensure a good sintering of the diamond table. In this case carbon dissolution from the diamond grains into binder phase would be more important and a carbon gradient takes place into the liquid. Then, a Co migration occurs from the interface towards the lower part of the substrate in the direction of carbon diffusion, as for the DP carbides [3], [4], [START_REF] Eso | Kinetics of cobalt gradient formation during the liquid phase sintering of functionally graded WC-Co[END_REF].

Conclusion

Reactive imbibition allows obtaining in one-step a dome-shaped gradient of about 300HV on 6mm without WC grain growth. After HPHT process, the magnitude of the gradient is preserved and sometimes improved. But the mean hardness of the substrate increases due to the migration of Co phase from the substrate into the diamond compact during the HPHT treatment.

The first mechanical tests show that the service life in abrasion of a graded PDC cutter is improved by 50% in comparison with a standard PDC cutter, tested in the same conditions. Its shock resistance exhibits an increase in performance by 50%.

Fig. 1 .

 1 Reactive coating mechanism proposed by Sorlier [14] ExperimentationA.

Fig. 3 .

 3 WC grain size distribution after each process step (a); Comparison of Co content measured along the revolution axis by three methods from a reactive imbibed substrate (b) Results

  (a) r= 0 (b) r=z Fig. 4. Microhardness profiles along the axis of symmetry (a); along the oblique axis r=z (b) for as-received (sample A) and reactive imbibed substrates (sample B) Sample A (z = 4mm) Sample B (z = 0.25mm) Sample B (z = 4mm) Sample B (z = 7.75mm) Fig. 5. SEM micrographs showing Co distribution in samples A and B on r =

11 .

 11 TEM observation of a WCoB phase (lighter phase) in a WC-Co substrate (a), diffraction of the WCoB phase on [110] (b), diffraction of the WCoB phase on [120] (c)

Table 1 .

 1 Summary of the samples

	WC-Co Cutters	HPHT Conditions
	As-received	Reactive Imbibition	Diamond	Pressure/Temperature
	A

Acknowledgements

The authors wish to thank PROFOR (Nouveau PROcédé d'élaboration d'outils à gradient de propriétés pour le FORage de roches abrasives en conditions sévères) partners for participating and the Agence Nationale de la Recherche (ARN-09-MAPR-441690) for funding this research.


Arthur Leclaire
email: arthur.leclaire@parisdescartes.fr

Lionel Moisan
email: lionel.moisan@parisdescartes.fr

Blind Deblurring Using a Simplified Sharpness Index

Keywords: global phase coherence, sharpness, blind deconvolution, noreference image quality assessment, oracle deconvolution filter

It was shown recently that the phase of the Fourier Transform of an image could lead to interesting no-reference image quality measures. The Global Phase Coherence, and its recent Gaussian variant called Sharpness Index, rate the sharpness of an image in contrast not only with blur, but also noise, ringing, etc. In this work, we introduce a new variant of these indices, that can be computed with one Fourier Transform only, hence four times quicker than the Sharpness Index. We use this new index S to build an image restoration algorithm that, in a stochastic framework, selects a radial-unimodal deconvolution kernel for which the S-value of the restored image is optimal. Experiments are discussed, and comparison is made with a radial oracle deconvolution filter and the recent blind deconvolution algorithm of Levin et al.

Introduction

No-reference image quality assessment consists in designing algorithms to evaluate the quality of an image (in particular in relation with its level of blur and noise) without requiring either an ideal version of this image (full-reference) or features extracted from this ideal image (reduced-reference). Finding good image quality (and sharpness) metrics has several applications, like, e.g., parameter selection, image restoration [START_REF] Zhu | Automatic Parameter Selection for Denoising Algorithms Using a No-Reference Measure of Image Content[END_REF], benchmarking, or depth estimation [START_REF] Calderero | Evaluation of Sharpness Measures and Proposal of a Stop Criterion for Reverse Diffusion in the Context of Image Deblurring[END_REF].

A way to address the notion of image quality is to think in terms of precision of its geometric elements (contours, alignments, etc.). Since the pioneering work of Oppenheim and Lim [START_REF] Oppenheim | The Importance of Phase in Signals[END_REF], it is well known that the geometry of an image is mainly encoded in the phase of its Fourier Transform. And yet, the phase information itself is still very difficult to understand. A first definition of local phase coherence was given in [START_REF] Morrone | Feature Detection in Human Vision: a Phase-Dependent Energy Model[END_REF], [START_REF] Kovesi | Phase Congruency: a Low-level Image Invariant[END_REF] and used for edge detection. Later, it was used to design a local sharpness measure in [START_REF] Wang | Local Phase Coherence and the Perception of Blur[END_REF] and [START_REF] Hassen | No-Reference Image Sharpness Assessment Based on Local Phase Coherence Measurement[END_REF]. In 2008, the authors of [START_REF] Blanchet | Measuring the Global Phase Coherence of an Image[END_REF] defined a notion of Global Phase Coherence (GPC), which rates the sharpness of an image depending on how the regularity of the image is destroyed as its phase information is lost. Very recently in [START_REF] Blanchet | An Explicit Sharpness Index Related to Global Phase Coherence[END_REF], a variant of GPC called Sharpness Index (SI) was introduced. It has the advantage of being described by an explicit closed-form formula, without needing computationally expensive Monte-Carlo Simulations like GPC.

In the present paper, we show that the SI metric can be further simplified to yield a new metric S that can be computed with only one Discrete Fourier Transform (versus four for the SI metric), while being an excellent approximation of the latter (Section 2). The behavior of this new metric is analyzed, in particular in comparison with the Q metric proposed by Zhu and Milanfar [START_REF] Zhu | Automatic Parameter Selection for Denoising Algorithms Using a No-Reference Measure of Image Content[END_REF] (Section 3). Then, a blind deblurring algorithm is built in Section 4, that looks for the linear filter that maximizes the S-value of the restored image while imposing the Fourier Transform of the convolution kernel to be radial, unimodal and smooth. The results of this algorithm are discussed, and compared with the corresponding linear oracle and with the blind deconvolution algorithm recently proposed by Levin et al. [START_REF] Levin | Efficient Marginal Likelihood Optimization in Blind Deconvolution[END_REF].

Global Phase Coherence and Derived Sharpness Metrics

Let us first introduce some useful notations. In all the following, we consider gray-level images u : Ω → R defined on a discrete M × N rectangular domain

Ω = Z 2 ∩ - M 2 , M 2 × - N 2 , N 2 .
The discrete Fourier transform (DFT) of u is the complex function û defined by The Ω-periodization of u is the image u :

∀ξ ∈ Z 2 , û(ξ) = x∈Ω u(x)e -i ξ,x , (1)
Z 2 → R that extends u to Z 2 by u(x) = u(x ′ ), where x ′ is the unique element of Ω such that x ′ -x ∈ M Z × N Z.
The gradient of u is defined by

∀(x, y) ∈ Z 2 , ∇ u(x, y) = ∂ x u(x, y) ∂ y u(x, y) = u(x + 1, y) -u(x, y) u(x, y + 1) -u(x, y) , (2)
and the (periodic and anisotropic) Total Variation of u is

TV(u) = ∂ x u 1 + ∂ y u 1 = x∈Ω |∂ x u(x)| + |∂ y u(x)| . (3)
The autocorrelation of ∇ u is the function Γ : Ω → R 2×2 defined by

Γ (z) = Γ xx (z) Γ xy (z) Γ xy (z) Γ yy (z) = y∈Ω (∇ u(y)) (∇ u(y + z)) T . (4)

Global Phase Coherence

As mentioned earlier, the phase of an image u encodes a great part of the geometry of u: if one reproduces the famous experiment of [START_REF] Oppenheim | The Importance of Phase in Signals[END_REF] consisting in imposing the phase of an image u to another image v, one can see on the result that several edges from u have appeared and all the geometric content of v has disappeared. Indeed, phase coefficients need strong alignment constraints in order to produce sharp edges and clean flat regions in an image. The Global Phase Coherence metric introduced in [START_REF] Blanchet | Measuring the Global Phase Coherence of an Image[END_REF] quantifies how the loss of this phase coherence affects the image regularity, measured by its Total Variation (3). More precisely, the phase of an image u is randomized to produce the Random Phase Noise image U ψ defined in Fourier Domain by

∀ξ ∈ Ω, U ψ (ξ) = |û(ξ)|e iψ(ξ) , (5)
where ψ : Ω → R is a uniform random phase (the coefficients of ψ are independent and uniformly distributed in (0, 2π), modulo the relation ψ(-ξ) = -ψ(ξ) ensuring that U ψ is real-valued, see [START_REF] Galerne | Random Phase Textures : Theory and Synthesis[END_REF]), which leads to Definition 1 (Blanchet, Moisan, Rougé, 2008 [START_REF] Blanchet | Measuring the Global Phase Coherence of an Image[END_REF]). The Global Phase Coherence (GPC) of u is the number

GPC(u) = -log 10 P(TV(U ψ ) ≤ TV(u)) . (6)
For an image u with sharp edges and clean uniform zones, the Total Variation is expected to be low amongst the ones of its phase randomizations. Therefore, for such an image, the probability of the event {TV(U ψ ) ≤ TV(u)} will be very small, and the value of GPC(u) will be large. That is why this phase coherence index (and the variants that follow) is expected to behave like an image quality measure. Note that without the logarithm in [START_REF] Desolneux | A Compact Representation of Random Phase and Gaussian Textures[END_REF], the values of GPC(u) would often cause a numerical underflow (a value like, e.g., 10 -1000 cannot be represented in most computer environments).

The main issue with [START_REF] Desolneux | A Compact Representation of Random Phase and Gaussian Textures[END_REF] is that no closed-form formula has been found so far to compute GPC(u), so that a computationally expensive Monte-Carlo simulation (coupled with a Gaussian approximation of the random variable TV(U ψ )) is proposed in [START_REF] Blanchet | Measuring the Global Phase Coherence of an Image[END_REF], which limits the potential application of the GPC metric.

Sharpness Index

Hopefully, a closed-form variant of GPC was recently found. In [START_REF] Blanchet | An Explicit Sharpness Index Related to Global Phase Coherence[END_REF], the periodic convolution of u with a conveniently normalized Gaussian white noise W is considered instead of U ψ , and the first two moments of TV(u * W ) are explicitly computed in function of u. Note that u * W is nothing but the natural Gaussian approximation of U Ψ , and these two random images only differ by a convolution with the texton of a white noise, which is close to a Dirac distribution [START_REF] Desolneux | A Compact Representation of Random Phase and Gaussian Textures[END_REF] (in Fourier Domain, u * W and U Ψ differ by a multiplicative Rayleigh Noise). Even if the exact law of TV(u * W ) seems difficult to compute, it is expected to be approximately Gaussian, which leads to Definition 2 (Blanchet, Moisan, 2012 [2]). The Sharpness Index of u is

SI(u) = -log 10 Φ µ -TV(u) σ (7)
where

Φ(t) = 1 √ 2π +∞ t exp - x 2 2 dx , µ = (α x + α y ) 2 π √ M N , (8)
σ 2 = 2 π z∈Ω α 2 x • ω Γ xx (z) α 2 x + 2α x α y • ω Γ xy (z) α x α y + α 2 y • ω Γ yy (z) α 2 y , (9)
α x = ∂ x u 2 = x∈Ω |∂ x u(x)| 2 1 2 , α y = ∂ y u 2 ,
Γ is the autocorrelation of ∇ u given in (4), and ω is the function defined by

∀t ∈ [-1, 1], ω(t) = t • Arcsin(t) + 1 -t 2 -1 . ( 10
)
In practice, the numerical computation of SI(u) requires the computation of T V (u), α x and α y (linear time), plus the three different components of Γ that can be computed with four Fast Fourier Transforms (a direct FFT of u, and 3 inverse FFTs for the cross correlations of ∂ x u and ∂ y u). The overall cost is hence dominated by these four M × N FFT computations, which represents a complexity of O(M N log(M N )) for well-suited image dimensions. Note that the function Φ is available in most mathematical libraries through the complementary error function (often written erfc), but when t is greater than say, 20, the following (almost exact) approximation is systematically used to avoid numerical underflow:

-log 10 Φ(t) ≃ t 2 + log(2πt 2 ) 2 log(10) . (11)

A Simplified Version of SI

We now introduce a new index S which is analytically close to SI and faster to compute. For that, let us observe that ω(t) is equivalent to t 2 /2 when t → 0. Therefore, this approximation can be used to replace σ 2 = Var(TV(u * W )) by

σ 2 a = 1 π z∈Ω α 2 x • Γ xx (z) α 2 x 2 + 2α x α y • Γ xy (z) α x α y 2 + α 2 y • Γ yy (z) α 2 y 2 ,
which after simplification yields the following Definition 3. The S-metric of an image u is

S(u) = -log 10 Φ µ -TV(u) σ a , (12)
where

σ 2 a = 1 π Γ xx 2 2 α 2 x + 2 • Γ xy 2 2 α x α y + Γ yy 2 2 α 2 y , (13)
and α x , α y , µ, and Γ are as in Definition 2.

Whereas SI needed all the coefficients of the gradient auto-correlation matrix, the S metric only depends on the overall energy of the three components Γ xx , Γ xy and Γ yy . Thanks to Parseval's formula, they can be computed in Fourier domain, recalling that

Γ xx (ξ) = | ∂ x u(ξ)| 2 = 2 sin 2 πξ 1 M |û(ξ)| 2 , Γ yy (ξ) = 2 sin 2 πξ 2 N |û(ξ)| 2 ,
and

| Γ xy (ξ)| = | ∂ x u(ξ)|| ∂ y u(ξ)| = 2 sin πξ 1 M sin πξ 2 N |û(ξ)| 2 .
The computation of S only involves the l 1 and l 2 norms of the gradient, and the FFT of u. Thus, the overall dominant cost is only one FFT (compared to 4 for SI), while the approximation of SI by S is very good as stated by

Proposition 1. We have 0 ≤ σ 2 -σ 2 a σ 2 a ≤ π -3 ≈ 0.142.
Proof. With the expressions of σ 2 and σ 2 a , one can write

σ 2 -σ 2 a = 2 π x∈Ω α 2 x ω Γ xx (x) α 2 x - 1 2
Γ xx (x) α 2 x 2 +2α x α y ω Γ xy (x) α x α y - 1 2 Γ xy (x) α x α y 2 + α 2 y ω Γ yy (x) α 2 y - 1 2 Γ yy (x) α 2 y 2 .
Besides, Taylor Formula applied to ω yields

∀t ∈ [-1, 1], 0 ≤ ω(t) - 1 2 t 2 ≤ ct 4 ≤ ct 2 , (14)
with c = ω(1) -1 = π-3 2 , and thus

0 ≤ σ 2 -σ 2 a ≤ 2c π x∈Ω α 2 x Γ xx (x) α 2 x 2 +2α x α y Γ xy (x) α x α y 2 +α 2 y Γ yy (x) α 2 y 2 = 2cσ 2 a .
3 Validation of S as a Quality Measure

As in [START_REF] Blanchet | Measuring the Global Phase Coherence of an Image[END_REF], we shall systematically apply two simple image transforms to an image u before computing S(u) with [START_REF] Moisan | Periodic Plus Smooth Image Decomposition[END_REF], in order to avoid periodization and quantization biases. First, as the S metric is defined (like GPC and SI) through a periodic setting, the periodic component of u (see [START_REF] Moisan | Periodic Plus Smooth Image Decomposition[END_REF]) is first extracted to avoid discontinuities across the image frame border. Then, a simple dequantization procedure (a (1/2, 1/2) sub-pixel translation with Fourier interpolation, see [START_REF] Desolneux | Dequantizing image orientation[END_REF]) is applied to ensure that the quantization of the original image (generally with 256 gray levels) does not artificially decreases its Total Variation. On Fig. 1, we give a first empirical evidence that S behaves as an image quality measure. Indeed, this blur-noise diagram shows that the value of S decreases as the the level of blur or noise increases, with a correspondence between noise and blur which is similar to GPC and SI (see [START_REF] Blanchet | Measuring the Global Phase Coherence of an Image[END_REF][START_REF] Blanchet | An Explicit Sharpness Index Related to Global Phase Coherence[END_REF]). The S metric not only decreases with respect to blur and noise, but unlike the Q metric of Zhu and Milanfar [START_REF] Zhu | Automatic Parameter Selection for Denoising Algorithms Using a No-Reference Measure of Image Content[END_REF], it also decreases when ringing artifacts (that may result from excessive deblurring) appear, as shown on Fig. 2. This suggests that S could be used in a parametric or even non-parametric blind deblurring algorithm, as will be done in Section 4. Note however, that the Q metric performs slightly better than S for parameter selection in the SKR denoising method presented in [START_REF] Zhu | Automatic Parameter Selection for Denoising Algorithms Using a No-Reference Measure of Image Content[END_REF]. In fact, the reason is logical considering the origin of S: the noise left by SKR in uniform zones is structured and the coherence of its phase makes S prefer less denoised images than Q.

Application to Blind Deblurring

Removing blur from a single image is a difficult task. If the blur is linear and spatially uniform, it can be modeled as a convolution. Several algorithms (see, e.g., the recent efficient scheme for T V -L 2 deblurring proposed in [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF]), have been proposed to invert the effect of this convolution when the blur kernel is known. Addressing the problem of blind deconvolution, i.e. when the kernel is not known, is even more difficult, and several solutions have been proposed in the last decades. In the present paper, we shall consider in particular (for comparison purposes) the very recent work of Levin et al. [START_REF] Levin | Efficient Marginal Likelihood Optimization in Blind Deconvolution[END_REF].

Here, rather than trying to reverse the effect of a convolution, we shall try to improve the image directly by convolving the blurry image u with a unitmass (that is, average-preserving) filter k that maximizes F u (k) = S(k * u). Indeed, Fig. 2 shows that, in a parametric Wiener deconvolution, the S metric is able to select the blur parameter. In [START_REF] Calderero | Evaluation of Sharpness Measures and Proposal of a Stop Criterion for Reverse Diffusion in the Context of Image Deblurring[END_REF], Calderero and Moreno made a similar observation for the SI metric in a context of reverse diffusion. In this section, we will show that S can be used for non-parametric blind deblurring. [START_REF] Zhu | Automatic Parameter Selection for Denoising Algorithms Using a No-Reference Measure of Image Content[END_REF] obtained when applying (in Fourier domain) a H 1 -regularized deconvolution filter to a natural image (Yale). The abscissa (ρ) corresponds to the standard deviation of the supposedly Gaussian convolution kernel. Over the results obtained for the whole range of ρ, one can see that contrary to Q, S attains a maximum for ρ = 0.55, which roughly corresponds to the value beyond which ringing artifacts begin to appear. This suggests that S is able to discriminate ringing, whereras Q does not.

Kernels with Compact Support

Algorithm 1 below can be used to optimize F u on particular sets of kernels k, as soon as the number of coefficients that define k remains small enough. For example, one can optimize F u on the set of symmetric 5 × 5 kernels (12 coefficients), or on the set of separable symmetric 21 × 21 kernels (20 coefficients, or 10 if the same one-dimensional kernel is used for each coordinate). In general, the results obtained with Algorithm 1 are good, but some images lead to interesting failure cases, in particular when regions with highly structured textures or dominant orientations are present (see Fig. 3). Since the functional F u is not concave, it does not necessarily have a unique local maximum, and a reason could be that Algorithm 1 does not manage to converge to the actual global maximum of the objective function F u . However, experiments suggests that the failure is more likely to be due to an inadequate set of kernels (in particular the set of separable kernels). To avoid such degenerated cases, and get rid of the small-kernel-support constraint, we consider in Section 4.2 other sets of kernels for which constraints are considered in Fourier domain.

Algorithm 1

-Begin with k = δ 0 (discrete Dirac kernel) -Repeat n times ⊲ Define k ′ from a random perturbation of k ⊲ If S(k ′ * u) > S(k * u) then k ← k ′ -Renormalize k to a unit-

Optimization of a Radial-Unimodal Kernel

Instead of imposing that the convolution kernel has a fixed compact support, we here consider the set of kernels that are radial in Fourier domain, with a unimodal profile, which is a plausible assumption for a deconvolution kernel. More precisely, we assume that the DFT of the restoration kernel k r is given by

∀ξ ∈ Ω, k r (ξ 1 , ξ 2 ) = L r 2(d -1) ξ 1 M 2 + ξ 2 N 2 , (15)
where

L r : [0, d -1] → R is the piecewise affine interpolate on [0, d -1] of the finite sequence r(0) = 1, r(1), r (2)
, . . . , r(d -2), r(d -1) = 0. This sequence is supposed to be unimodal, which means that there exists a value i m (mode index) that satisfies ∀i < i m , r(i + 1) ≥ r(i) , and ∀i ≥ i m , r(i + 1) ≤ r(i) .

One possible perturbation strategy for this set of kernels consists in the addition of a uniform random value to a randomly chosen coefficient of r, followed by a projection on the set U of unimodal sequences (this projection can be computed in O(n 2 ) operations with the Pool Adjacent Violators algorithm [START_REF] Frisen | Unimodal Regression[END_REF]). We observed that with this strategy, moving the mode position was difficult, so we relaxed the unimodality constraint and incorporated in the objective function the l 2 -distance d(r, U ) between r and the set U of unimodal sequences. We also found useful to add the possibility to increase the regularity of the radial profile r by incorporating a term depending on

r 2 H 1 = d-2 i=0 (r(i + 1) -r(i)) 2 . (16)
Finally, the objective function (to be maximized) is

F u (r) = S(k r * u) -λ d(r, U ) -µ r H 1 , (17)
where λ and µ are weighting parameters.

In order to maximize F u , we used Algorithm 2 below. We observed that n = 10000 was sufficient to ensure convergence on r (the relative changes after 10000 iterations were less than 10 -3 ); moreover, we checked that several realizations of this stochastic technique led to the same local maximum. The other parameters were set to d = 20, initial i m = 5, a = 0.1, λ = 10000, and µ varying from 0 to 100. Let us comment the choice of i m . As soon as several local maxima are present, the result of an optimization technique may depend on the initialization, and the natural solution would be to apply the algorithm for all possible values of the initial mode index. But in a wide majority of cases, we observed that the result of this algorithm was not depending on this initial value. In a few cases, two different local maxima could be found, but the higher value of the objective function F u was always obtained for i m ∈ [d/4, 3d/4]. This is why it seems empirically sufficient to run the algorithm only once with the initial value of i m in that range.

Algorithm 2 -Initialize r with the piecewise-linear profile such that r(0) = 1, r(i m ) = 2, and r(d

-1) = 0. -Repeat n times ⊲ Pick a random index i ∈ [1, d -2] ⊲ Draw a random value ε ∈ [-a/2, a/2] with uniform distribution ⊲ Set r ′ ← r, then r ′ (i) ← r(i) + ε ⊲ If F u (r ′ ) > F u (r) then r ← r ′ -Return r, k r and k r * u

Results

We first used Algorithm 2 with µ = 0 on Room image, and checked that the failure of Algorithm 1 was avoided (Fig. 3, right).

Then, to produce the results shown in Fig. 4, we took two classical images (Capitol and Parrots) and corrupted them with a Gaussian blur kernel (width 1 pixel) and an additive white Gaussian noise (variance 1). We then applied several deblurring algorithms (detailed below) and evaluated their performances by computing their respective PSNR values with respect to the original clean image. Notice, however, that for blind deblurring tasks the PSNR value is not very reliable (in particular because even the original clean image is necessarily, in some sense, blurry and noisy), and visual inspection is often preferable to compare the algorithms.

First, we used Algorithm 2 with µ = 0 and µ = 10. We observed that the results were stable, and that the restoration resulted in a significant sharpness increase. However, for µ = 0 some low-frequency noise is still visible on uniform zones. Increasing the value of µ to µ = 10 reduces the residual noise (because it reduces the amplitude of L r , that is, the amplification of noisy Fourier coefficients), but also attenuate some details in textured zones.

We then compared these results with the state-of-the-art blind deconvolution algorithm of Levin et al. [START_REF] Levin | Efficient Marginal Likelihood Optimization in Blind Deconvolution[END_REF]. One can see in Fig. 4 that Algorithm 2 is more precise on the fine details of the image, but it also keeps much more noise than this method. In fact, the result of [START_REF] Levin | Efficient Marginal Likelihood Optimization in Blind Deconvolution[END_REF] is really "clean" and has a small Total Variation, which explains incidentally why its S value is significantly larger compared to Algorithm 2. Notice also that the method [START_REF] Levin | Efficient Marginal Likelihood Optimization in Blind Deconvolution[END_REF] is more general, and has been shown to perform particularly well in the case of a motion blur, while the radial constraint of Algorithm 2 cannot handle such motion blurs.

Another interesting experiment consists in computing the optimal kernel k o that maximizes the expected distance between the reconstructed image and the clean image u 0 , knowing the parameters (kernel κ and noise level β) of the degradation process (this is an oracle since neither these parameters nor the clean image are supposed to be known). One has

k o = Arg min k E u 0 -k * (κ * u 0 + βW ) 2 (18)
and if all kernels k were considered, the solution would be given by

∀ξ ∈ Ω, k o (ξ) = κ(ξ) | û0 (ξ)| 2 |κ(ξ)| 2 | û0 (ξ)| 2 + σ 2 M N . (19)
Now, since we only consider kernels whose DFT is built from a radial profile (linearly interpolated on d points), one can show that the optimal radial profile is the minimum of a quadratic function and thus can be obtained by solving a small linear system. In Fig. 5 we can see that the profile of the oracle radial kernel is unimodal, with a mode at a position which is close to the one estimated by Algorithm 2. The restored images obtained with this oracle filter are also displayed in Fig. 4: they are a little more precise on the details, but present a significant amount of structured noise; indeed, such a noise is not very costly for a l 2 risk function. This is somehow reassuring: this shows that the structured noise also appearing with Algorithm 2 (in particular with µ = 0) is truly a limit of techniques based on linear filtering.

Conclusion

We introduced a new variant of the GPC and SI image quality metrics, that can be computed four times faster than SI. This new index S provides a sharpness measure that can be used in a stochastic optimization framework to achieve blind deblurring through linear convolution with a radial unimodal kernel. Though suffering from the limits of linear filtering, the obtained results are convincing, and visually similar to the best possible ones (oracle) obtained by such an approach. The extension to motion blur kernels could be an interesting generalization, as well as the use of more sophisticated (non-linear) restoration techniques to make the best usage of the selection performances of the S metric.

  where ξ, x = 2π x1ξ1 M + x2ξ2 N with ξ = (ξ 1 , ξ 2 ) and x = (x 1 , x 2 ). The function |û| will be called the modulus of u. A phase function for u is any function ϕ : Z 2 → R such that for all ξ ∈ Z 2 , one has û(ξ) = |û(ξ)|e iϕ(ξ) .

Fig. 1 .

 1 Fig.1. Some level lines of the function (r, β) → S(gr * u + βW ) where gr is the 2-D Gaussian convolution kernel with standard deviation r, and W is a white noise image with unit variance in each pixel. The absolute values of S and the exact shape of the level lines depend on the image considered (here, Barbara), but the overall shape remains similar.

Fig. 2 .

 2 Fig.2. These diagrams plot the proposed metric S (left) and the metric Q (right) of Zhu and Milanfar[START_REF] Zhu | Automatic Parameter Selection for Denoising Algorithms Using a No-Reference Measure of Image Content[END_REF] obtained when applying (in Fourier domain) a H 1 -regularized deconvolution filter to a natural image (Yale). The abscissa (ρ) corresponds to the standard deviation of the supposedly Gaussian convolution kernel. Over the results obtained for the whole range of ρ, one can see that contrary to Q, S attains a maximum for ρ = 0.55, which roughly corresponds to the value beyond which ringing artifacts begin to appear. This suggests that S is able to discriminate ringing, whereras Q does not.

Fig. 3 .

 3 Fig.3. Blind deblurring of the original Room image (left). Algorithm 1 is applied for two different sets of kernels: 5 × 5 symmetric kernels (second column) and 21 × 21 symmetric separable kernels (third column). The right image is the result obtained with the method proposed in Section 4.2 (Algorithm 2, µ = 0). We observe that Algorithm 1 fails in both cases (probably in reason of the large striped texture), while the radial-unimodal constraint imposed in Algorithm 2 yields a nice-looking result.

OriginalFig. 4 .

 4 Fig. 4. Blind deblurring of a degraded version of Capitol and Parrots images (Gaussian blur of width 1 pixel plus Gaussian noise with variance 1). We present in each case the original image, the blurred and noisy input, the result of Levin et al. algorithm [11], the oracle output (best possible result obtained by a radial unimodal convolution filter) and the results of Algorithm 2 with µ = 0 and µ = 10. The PSNR values are computed in each case with respect to the original image.

Fig. 5 .

 5 Fig. 5. Different radial profiles (oracle, and Algorithm 2 with µ = 0 and µ = 10) obtained on images Capitol (left) and Parrots (right).

Acknowledgments. This work has been supported by the French National Research Agency under grant ANR-09-BLAN-0029-01.