Elodie Roullot
email: roullot@esme.fr

Jean-Eric Autegarden

Patrick Devriendt

Francisque Leynadier

Segmentation of Erythema from Skin Photographs for Assisted Diagnosis in Allergology

More than 2 people out of 10 suffer from allergies, which can take various forms, from eczema to anaphylactic reactions with possible lethal consequences. Diagnosis is achieved through so-called prick-tests or IDR (intradermo-reaction): the injection of a small quantity of substances suspected to cause the allergic manifestation induces an erythema, the size of which is a useful indicator for the diagnosis. The manual surface measurement is timeconsuming and inaccurate. This article presents a method for the semi-automatic measurement of the erythema from a photograph of the skin, taken in such conditions that lighting problems are minimized. The method is based on region growing and takes advantage of the most significant color spaces; the Lab space appears to be the best suited. It was tested on nearly 100 images, taken by various operators, on patients with various skin pigmentations; it gave promising results and proved to be robust.

Introduction

Allergies are classified as the fourth most important public health problem in the world. Numerous epidemiologic studies show that 20 to 25% of people suffer from allergic troubles. Allergic manifestations take various forms, among which eczema, rhinitis, asthma, or even Quincke oedema or anaphylactic reaction, with possible lethal consequences.

In order to prevent from the most severe anaphylactic manifestations, it is necessary to have an acute knowledge of the substances that are responsible for the allergic reaction, called the allergens . This diagnosis is achieved by the allergologist through cutaneous tests: prick-tests or IDR (intra dermo reaction). In prick-tests, a drop of solution containing the allergen to be tested is dropped on the forearm of the patient, which is then punctured in order to make the allergen penetrate into the skin ; in IDR s, the allergen is injected into the skin with a syringe.

An allergic reaction locally induces a dilatation of blood vessels and the release of chemical substances due to the axon reflex [START_REF] Nischik | Analysis of Skin Erythema Using True-Color Images[END_REF], leading after 15 to 20 minutes to the apparition of an erythema, characterized by a local reddening of the skin and an in-crease in its temperature; at the center of the erythema a papula is formed, looking like a mosquito bite. The surface of the erythema is an indicator for the severity of the allergy.

The mean diameter of each erythema is measured with a scale; each value, corresponding to a given allergen, is compared to 2 reference values: the positive and negative standards. The positive standard indicates the standard degree of reaction of the patient and is achieved by the injection of codein phosphate which involves for everybody a positive reaction. The negative standard is achieved by the injection of physiological serum which never involves a positive reaction (only the puncture is involved in an eventual change in the skin appearance). The reaction induced by a given allergen is considered to be positive if it is equivalent to or larger than the positive standard; if it is comparable to the negative standard, the reaction is considered to be negative.

This technique appears very fastidious especially due to the huge number of allergens to be tested (up to 20); moreover, the measures are corrupted by errors due to the inter-and intra-operators variability, thus making a follow-up of the patients more difficult. Moreover, the diameter of the erythema does not give a good estimation of the surface, the measure of which would be preferred by the allergologists.

All these reasons make an automatic measure of the surface very challenging; the expected improvements are: less time spent for measuring, less intra-and interoperators variations in the measurement, more accuracy thanks to the surface measurement instead of the diameter measurement, and the storage possibility, making a follow-up of the patients easier. However, the clinical motivation of this work is not to provide a fully automated diagnosis system, but rather to provide assistance to the allergologist to make his diagnosis. Thus, it should be focused on the segmentation problem but it should not be dealt with an automatic diagnosis whether the patient is allergic or not.

Several research works emerged during the last years, dealing with the automation of erythema detection. They make use either of spectroscopy, or of echography, or of visible imaging. Several teams including Kopola et al. [START_REF] Kopola | Two-channel fiber optic skin erythema meter[END_REF], Dawson et al. [START_REF] Dawson | A theoretical and experimental study of light absorption and scattering by in-vivo skin[END_REF], Ferguson-Pell et al. [START_REF] Ferguson-Pell | An empirical technique to compensate for melanin when monitoring skin microcirculation using reflectance spectrophotometry[END_REF], Hajizadeh-Saffar [START_REF] Hajizadeh-Saffar | An investigation of factors affecting the accuracy of in vivo measurement of skin pigments by reflectance spectrophotometry[END_REF] or Diffey [START_REF] Diffey | Quantitative aspects of ultraviolet erythema[END_REF] suggest calculating an erythema index from reflectance measures in different carefully chosen spectral bands, by means of a spectroscope. Several methods among these were evaluated by Riordan et al. [START_REF] Riordan | Testing the validity of erythema detection algorithms[END_REF]. Some of them allow detecting erythemas even in highly pigmented skins; however, all of them only provide an index (generally an erythema index and a melanin index) that characterizes the global change in color, whereas routinely the allergologist is interested in the size of the erythema.

Seidenari [START_REF] Seidenari | The Quantification of Patch Test Responses: A Comparison Between Echographic and Colorimetric Methods[END_REF] et al., for their part, use echography to quantify erythemas induced by patch-tests; once again, this provides values representing the intensity of the reaction whereas the allergologist is interested in the surface.

The surface can only be obtained automatically from images of the skin; although the problem of segmentation is dealt with very often in the context of melanoma detection for example, this problem has been addressed very rarely in the context of erythema detection. Unfortunately, the methods used for segmenting melanomas can not be transposed easily for erythema detection, since melanoma images are much more contrasted than erythema images. Nischik et al. addressed the problem of erythema segmentation from visible light images and proposed a method for segmenting the erythema, based on the CIELab color space and on the comparison between two images acquired before and after prick-testing [START_REF] Nischik | Analysis of Skin Erythema Using True-Color Images[END_REF]. This method presents a major drawback: it implies a motion analysis in order to suppress artifacts due to the unavoidable motion between both images, since the reaction can be read 15 to 20 minutes after prick-testing.

We therefore propose an original approach to segment the erythema and provide a surface measurement of it, using only one image taken after reaction. In Section 2 we briefly describe the image capture system and the image database; then we address in Section 3 a colorimetric study in order to extract the best color space to use. In Section 4 we describe the method, and the results are presented and discussed in Section 5.

Image Database

As said before, the detection is performed from digital photographs of the arm and forearm. These digital photographs are taken using a standard color CCD camera with a spatial resolution of 752 x 582 pixels. The camera is connected to a personal computer via a digitalization card (Ellips Riowin) allowing real-time visualization of the acquired image and eventual tuning of some acquisition parameters, as well as saving a color image in bitmap format (BMP).

The illumination system is composed of two neon tubes placed behind unpolished and translucent Plexiglas plates, in order to provide an homogeneous illumination and thus to avoid reflection artifacts. 90 images of the arm or forearm were acquired at anonymous hospital ( town , country ) from January 2003 to June 2004 on both male and female adult patients who agreed to take part in these tests. Most patients had lightly pigmented skins (Caucasian for most of them, and Metis or Mongoloid), and only one patient had a very dark skin (Negroid). Ten randomly chosen images formed the learning database and the others were used for evaluation purposes.

Colorimetric Study

A first examination of the learning database revealed that the classical RGB (red, green, blue) color space is not the best suited for the segmentation of erythemas. In this section, we study several color spaces to extract the best suited for discriminating the reddened skin from the normal skin. The CMY color space is the complementary color space to the RGB one: it is based on subtractive synthesis (as in printers) while the RGB space is based on additive synthesis (as in monitors). HSV (hue, saturation, value) and HSL (hue, saturation, luminance) systems correspond to a psychovisual description of color, based on the human eye function. The LMS color space is also based on the human eye: each stimulus L (long), M (middle) and S (short) represents the spectral sensitivity of the photoreceptive cells of the eye. The CIELab color space is derived from the XYZ standard color space, defined by the CIE in 1931 by a homographic transformation of the RGB space. L represents the luminance whereas a and b respectively represent the red-green chrominance and the blue chrominance.

All these values (R, G, B; C, M, Y; H, S, V; H, S, L; L, M, S; X, Y, Z; L, a, b) were computed for each image of the learning database. An SNR study revealed that the chrominance a provides the best discrimination of erythematous skin versus normal skin. These results confirm those obtained by Fullerton et al. [START_REF] Fullerton | Guidelines for the Measurement of Skin Color and Erythema[END_REF] and Weatherall et al. [START_REF] Weatherall | Skin Color Measurements in Terms of CIELAB Color Space Values[END_REF]; they recommend the use of the Lab color space, which takes into account the non linearity of color perception by the human eye.

Color-based Segmentation

Choice of the Method

As said before, the segmentation of erythema is a difficult task, because of the colorimetric heterogeneity of the skin (both the normal and the reddened skin present a marmor aspect), as well as because of the low contrast between erythema and normal skin. This problem belongs to the family of problems dealing with the segmentation of low SNR images. For such images, classical segmentation methods such as contour detection approaches are not well suited, that is why we concentrate our attention on constrained region segmentation approaches. Among these methods, structural ones are well suited to the problem since they allow constraining more or less the pattern of the extracted region [START_REF] Cocquerez | Analyse d images : filtrage et segmentation[END_REF].

Moreover, each image acquired represents the full arm or forearm, and can thus contain several erythemas corresponding to different allergens. It does not appear judicious to develop a system that would be able to recognize automatically to each kind of allergen each area of the arm corresponds, since it will induce many technical constraints for the operators performing the test, whereas it is not very constraining for him/her to have a minimum level of interaction with the system. On the contrary, the medical specialists generally prefer keeping a given degree of control on the automatic diagnosis process rather than obtaining results from a fully automatic system in which they have low confidence. Thus, it appears as a reasonable compromise to require the allergologist to click on the corresponding erythema when prompted with the name of an allergen, for example, as follows: ? launch image acquisition ? for each allergen:

? if an erythema is present: ? mouse click within the reaction area ? detection is thus performed automatically from the given seed ? view, validate and store the results, both quantitative (image) and qualitative (surface) ? else, store the result (absence of reaction) and skip to next allergen.

Taking into account both the considerations about the segmentation problem and the integration of the algorithm within a future computer-aided diagnosis software, a method based on region growing with manual input of the seeds appears consistent. In the next subsection we describe our algorithm based on pixel aggregation, with a strong constraint on the connectivity of the region and a weak constraint on its pattern.

Segmentation Algorithm

The algorithm is based on region growing [START_REF] Cocquerez | Analyse d images : filtrage et segmentation[END_REF][START_REF] Maître | Le traitement des images[END_REF] where the initial regions have a unitary size and are chosen manually by mouse clicking. Aggregation of new pixels is iterative and depends on two criteria: adjacency and colorimetric similarity. Each iteration begins with the constitution of a candidate set, based on adjacency; all candidates satisfying the colorimetric criteria are aggregated, and the new region is then used for the following iteration. The algorithm stops when no more candidates are eligible for aggregation.

The different steps of the algorithm are detailed below: ? Initialization: manual choice of N seeds chosen within the erythema, forming the N initial regions R i . ? For each region i:

? selection of the candidate set C i , i.e. the set of pixels adjacent to the current region, obtained with a morphological dilatation:

C i = (R i ⊕ E S ) R i (1)
where ⊕ represents the m orphological dilatation with a 3x3 square structural element E s .

? among the candidates c i k , selection of those satisfying the colorimetric similarity criterion:

)

) (a(R T ) a(c i k i > (2)
where a represents the chrominance a in the Lab color space, T is an automatic threshold computed from statistics in the region R i .

? if no candidate was selected at the previous step, algorithm ends. ? else, creation of the new region R i+1 as the union of region R i and selected pixels c i k .

Automatic Threshold Tuning

A statistical study of the chrominance a was performed on the 10 images of the learning database: a coarse thresholding helped select the arm and suppress the background, and then the histogram of the chrominance a was computed. This study revealed that the mean and variance of the chrominance varied sensibly from an image to another. Therefore the threshold was chosen as a function of the mean and variance of the chrominance a . For the 10 images of the learning database, the erythema sur-face was computed for various thresholds and compared to the reference surface (which is described further in Section 5). The threshold T varies as follows:

) ( ) ( I I a T a λσ + = (3)
where a is the mean chrominance over the image I, a σ is the variance of a over the image I, and the parameter λ varies from 0.6 to 1.4. The value of λ giving the exact diameter is interpolated; over the learning database, the mean? λ equals 0.97 meaning that the best threshold is

) ( 97 . 0 ) ( I I a a σ + .
The robustness of the parameter λ can also be evaluated from these tests: a variation of 10% of the standard deviation in the threshold induces a variation of about 5% in the diameter. We can therefore affirm that the threshold is quite robust.

Results and Discussion

Qualitative Study

The algorithm was tested on the 90 erythema images of the database. Among the 90 images, 89 concern patients with lightly pigmented skins (Caucasian, Mongoloid, Metis) and one patient has a highly pigmented skin (Negroid). Some results are shown on Figure 1 for the clear skins and on Figure 2 for the dark skin. On these images, blue pencil marks can be noticed, corresponding to the marks used by the physician; note that these marks do not appear in the a chrominance image and therefore do not disturb the algorithm (this would not be the case if a red pencil was used, as it was the case for some images that were therefore useless). When the system will be used instead of the manual measurement, these marks will not be needed anymore and will thus disappear. At this step of development, they are useful to compare the contour provided by the algorithm with the diameter estimated by the physician.

On Figure 1, it can be noticed that the lighting conditions differ considerably from one image to another; this does not seem to affect the algorithm, and confirms its robustness. Concerning the dark skin on Figure 2, the erythema of the unique negroid patient taking part in the tests could be detected in a satisfactory way.

In 5 images, the algorithm was not able to detect any contour. Possible reasons for these failures are: the low contrast of the image, due either to the image capture or to the skin reactivity of the patient; artifacts on the patient skin; defaults in the lighting systems providing an image with reflection artifacts that disturb the segmentation process.

The detection of an erythema lasts a few seconds, so that for one patient tested for 10 allergens, the total time including the mouse clicking, the automatic detection and storage of the results amounts to far less than one minute: less than the time for the patient to get dressed after the measurements. Routinely, the operators needs to draw the limits of the erythema before measuring (when measuring the diameter, the operator has to flatten the skin, which makes the skin reddening disappear), then measure, and finally report the value in a paper or electronic report; all these steps require sev-eral minutes and are much more time-consuming than the mouse-clicking and detection. Moreover, as said before, the measures are corrupted by an important variability between different operators and even between different measures taken by the same operators. The automatic detection provides a reproducible measure and proved to be nearly independent of the manual input of the seed (if the seed is chosen within the erythema and not on the border). Fig. 1. Examples of results obtained on various skin types among Caucasian, Mongoloid, and Métis patients. Left, the regions of interest in the original images; Right, the results superimposed on the original regions of interest. Blue pencil marks can be noted, which correspond to the separations between the allergens being tested, and to the limits used by the physician for the diameter measurement. Fig. 2. Example of result obtained on a skin of type Negroid. On this unique example it seems that the algorithm does not encounter any particular problem detecting the contour of the erythema, although the erythema is hardly visible with the naked eye.

Quantitative Evaluation

Two options can be considered for evaluating the algorithm; the first one consists in comparing the results to the clinical measurements performed manually, and the second one consists in comparing to a ground-truth surface that remains to be defined.

Actually, the measurements conducted routinely are diameter measurements performed by the nurses with a scale. These measurements suffer from a very important intra-and inter-operators variability, mainly due to two reasons. Firstly, the use of the scale implies a deformation of the natural curve of the arm and pushing on the skin implies color changes that can affect the measurement. Secondly, the erythema contour is not necessarily a circle; the nurses take as diameter measurement the biggest distance between the erythema borders. This is not so easy to appreciate with the human eye and is another factor for the variability of the measurement. Therefore, a quantitative evaluation based on these measurements does not appear to be judicious, and it appears necessary to perform evaluation by comparing to a ground-truth value.

To this end, for each image of the database an expert was asked to draw the erythema contour as precisely as possible; actually he pointed as many points on the contour as he needed and the corresponding polygon (with many vertices) was considered to be the ground-truth boundary. A goodness measure was defined as the ratio between the provided surface and the ground-truth surface. Figure 3 shows the histogram of this ratio computed over the database.

Figure 4-(a) presents a good result of evaluation with the polygon and the detection result superimposed on the image. On this example, the evaluation metric (the surface ratio) is 95%; the area of the polygon is 6411 pixels² and that of the automatic contour is 6139 pixels². A fictive diameter was computed as π area 2 , giving 90 pixels for the ground-truth, 88 pixels for the automatic detection. These values can be compared to the diameter measured manually (corresponding to the pencil parks) to 111 pixels. The error of 2 pixels is very low compared to the 21 pixels difference between the manual measure and the ground-truth. The difference between the manual measurement and the ground-truth provides a surface ratio of 152%, which corresponds to the worst cases of the automatic detection, as shown on Figure 3 (left).

Figure 4-(b) and 4-(c) shows two of the worse results obtained with the automatic detection. For the left image, the surface ratio is 76%; for the right image, the surface ratio is 155%. When comparing the results of Figure 4, it can be noticed that even in the worse cases, the contour is not so bad and is anyway at least as good as the one that could have been drawn manually. This illustrates the fact that surface ratios of 80% or 150% for example are considered as acceptable by the experts, and also contribute to the improvement of the measurements when compared to the manual ones. Fig. 3. Left: Histogram of the ratio between the surface provided by the algorithm and the ground-truth surface. This surface comparison leads to ratios that go very fast far away from the ideal ratio of 100%; this is due to the fact that the variation of surface is quadratic. Right: the same histogram after correction of the quadratic effect due to the surface measurement. As a consequence the ratios are homogeneous with diameter ratios and are thus easier to understand intuitively.

Conclusion

In this article, we addressed the problem of the automatic detection and quantification of reactions to allergic tests. This original problem was not addressed often in the literature, although it derives from an actual query from the allergologists.

In this study we focused on the segmentation problem; the automatic decision process was not addressed. Indeed the segmentation problem is the most difficult to solve because of the low contrast of erythema images and of the heterogeneity of the skin color. To this end a method was developed that combines color image analysis with a region growing algorithm; it provides a robust and reproducible contour detection from the manual input of a seed point within the erythema.

The method was applied on a database of nearly 100 images acquired during the consultations of allergologists at anonymous hospital. A rigorous evaluation allowed proving that the method has a satisfying detection rate and is quite robust. The physicians judged the software potentially helpful and very promising. In the future, it will be necessary to perform more tests especially on dark skins, since they were not well represented in the present database.

It is also planed to work on the acquisition and lighting systems: the system described in Section 2 is a provisory one and work is already being performed on a best suited system, such as those used by some dermatologists consisting of a digital camera together with a ring light.

Fig. 4 .

 4 Fig. 4. Comparison between the automatic contour (solid line), the ground-truth (dashed line) and the manual diameter (on image (a): dark pencil marks looking like parenthesis) in a good case (a) and in two bad cases ((b): the surface ratio is 76%; (c): the surface ratio is 155%).

Proceedings of ICAPR 05, Bath, UK, August 2005.


S Msolli
email: smsolli@enit.fr

O Dalverny

J Alexis

M Karama

Implicit integration scheme for porous viscoplastic potentialbased constitutive equations

Keywords: Viscoplastic, Porous, Potential, Integration scheme, FEM

This paper deals with a viscoplastic potential-based model allowing thermomechanical damage behavior modeling of porous materials. The model describes rate dependent effects, hardening, creep as well as defects coalescence and propagation. Kinematic and isotropic hardening effects are taken into account by a set of internal state variables. The integration and implementation of the model into the FE code using a fully implicit integration scheme is exposed. Finally, its used to predict mechanical behaviour degradation of solder layers used in power electronic packaging. Stress-strain behaviour and the evolution of volumic fraction of voids for the material under cyclic loading are presented.

Introduction

Solder joints are basic parts in power electronic modules which are subjected to severe thermomechanical cyclic loading, involving fatigue and creep at high temperatures [START_REF] Tien | Creep-fatigue interactions in solders[END_REF][START_REF] Micol | An investigation into the reliability of power modules considering base plate solders thermal fatigue in aeronautical applications[END_REF]. Many factors could have influence on the solder behaviour. Among these factors, we note the metallurgical aspect and the initial state of the solder material, the amount of porosity in the layer, the temperature history during the manufacturing stages. All these factors induced a set of defects of solder matrix, and as resulting effects, they modify their thermomechanical behaviour leading to a fall in solder lifetime and so of the whole packaging assembly. Volumic fraction of porosity may be considerable within the total volume of the solder layer and thus leads to weaken the whole thermomechanical behaviour of the junction [START_REF] Dudek | Three-dimensional (3D) visualization of reflow porosity and modeling of deformation in Pb-free solder joints[END_REF]. The fact that lifetime is influenced by defects coming before and during service, proves the need of pertinent theoretical models which are able to reproduce entire material behaviour during thermal-mechanical loading cycles and to accurately predict material degradation. The models have to take into account hardening variables and porous damage for viscoplastic materials which still a challenging task. Hardening rate controlling parameter

Nomenclature

Formulation of the viscoplastic constitutive equations

We introduce an internal state variable based unified viscoplastic model combined with a porous damage as first step for the viscoplastic-damage predictions. The damage model is based on Haghi and Anand framework [START_REF] Haghi | A constitutive model for isotropic, porous, elastic-viscoplastic metals[END_REF], which considers a Duva and Hutchinson approaches [START_REF] Duva | Constitutive potentials for dilutely voided nonlinear materials[END_REF], for the evolution of the volumic fraction of voids. As for the existing unified models, the constitutive equations are composed of a flow law for the inelastic deformation and a set of evolution equations for the state variables. From the thermodynamic viewpoint, the inelastic flow is derived from the free Helmholtz energy. This free energy depends on the volumic fractions of voids f and the internal state variables. Once the free energy expressed, then the flow law and the state variables can be directly derived from it. In our case, the kinematic and isotropic hardening rules are expressed respectively as:

* + * + ( ) ( ) ' ( ) ( ) ( ) 0 ' 1 in k k k ij k ij k k k k ij ij ij ij k k k i i k k k C b C p f C b T C b b T C T b T c g c c o e c c c e y sc o e Ã Ô / - Ä Õ Ç Ä Õ ? / • • • È Ù Ä Õ - - / Y • • • Ä Õ È Ù É Ú Å Ö (1)
* + * + * + 1 f p T T e oe e e
o e e oe s oe

Ã Ô • ? / / - / Y Ä Õ • Å Ö (2)
Where

2 : 3 in in ij ij p g g ?
is the effective inelastic strain rate. Thus, the total backstress is expressed as the sum of the ( )

k ij c such as ( ) k ij ij k c c
? Â . Note that the kinematic hardening formulation could represent hardening, dynamic recovery and static thermal recovery effects. Description of thermal variation effects is also improved by introducing the energy activation term in the static thermal recovery term.

Lets remember that for virgin materials, the inelastic flow rate is derived from a viscoplastic potential f, which could be expressed in our case as follows:

1 1 exp 1 1 m D A B B m z - Ç H ? S È Ù Ã Ô É Ú - Ä Õ Å Ö (3)
Where is expressed as:

( ) 0 k ij ij k J R R D u c z Ã Ô / / / Ä Õ Å Ö ?  (4)
and, e is a thermal diffusivity parameter expressed as an Arrhenius function.

Í Í Ì Í Í Ë Ê Õ Õ Ö Ô Ä Ä Å Ã Ý Ü Û Ì Ë Ê - Õ Ö Ô Ä Å Ã / ? S Õ Ö Ô Ä Å Ã / ? S 2 for 1 2 ln 2 exp 2 for exp m m m m T T T T kT Q T T kT Q (5)
<g> is defined as g if g oe 0 else g = 0. By introducing the volumic fraction of voids f, the viscoplastic potential for porous materials is calculated using Duva and Hutchinson approaches based on the expression: f p = F×f (6) F is given by Haghi and Anand as function of the effective stress triaxiality X and takes the form:

* + (1/ ) 1 /2 2 2/(1/ ) 1 3 1 2 m m f F C mX g m - - Ç Ê Û ? - - È Ù Ë Ü Ì Ý È Ù É Ú (7)
Where

( ) ( ) k mm mm k k ij ij k X J u c u c / ? Ã Ô / Ä Õ Ä Õ Å Ö Â Â and * + 1/ 1 1 m m C f ? /
. C is used here for voids interaction effects.

Usually, g is assumed to be constant and takes the value of 1.3 for a corresponding value of 0.161 for m. From Zavaliangos and Anand [START_REF] Zavaliangos | Thermo-elasto-viscoplasticity of isotropic porous materials[END_REF], an alternative form of g(f,m) is applicable for a wide range of m and is written as:

* + * + * + 2/ 1 1/ 0 , 1 m m g f m F f - Ç ? / È Ù É Ú (8)
Where

* + 1/ 2 2 0 2 1 3,34 0, 25 1 m f f F f Ã Ô - - ? Ä Õ / Å Ö (9)
The evolution of the volumic fraction of voids is governed by the next equation taking into account void nucleation and neglecting elastic compressibility.

* +

1 in G N mm N f f f f A p g ? - ? / - (10)
where

2 1 exp 2 2 N N N N N f p A s s g r Ã Ô Ã Ô / Ä Õ ? / Ä Õ Ä Õ Å Ö Å Ö
means that nucleation is governed by a Gauss distribution.

Finally, the normality rules permit to determine the inelastic flow law by deriving the porous viscoplastic potential expression in order to have,

p in ij ij g u •H ? • (11)
The built model in its actual form makes it possible to take into account kinematic and isotropic hardening, time-dependent and time-independent effects, and degradation due to the volumic fraction of voids.

Development of the integration algorithm

The model is composed of a set of differential constitutive equations. Weber and Anand developed an integration scheme for isothermal deformations of fully dense materials [START_REF] Weber | Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelasticviscoplastic solids[END_REF]. Zavaliangos and Anand presented an extension of this scheme for porous materials based on Zavaliangos framework [START_REF] Zavaliangos | Towards a Capability for Predicting the Formation of Defects During Bulk Deformation Processing[END_REF]. The aim of this fully implicit, fully coupled integration scheme is to find the solution of the following equation.

* +

in in ij ij ij h t g g u •H F ? F / F • ( 12
)
If we assume that * +

1 n in ij g - F
is the solution at iteration n+1, then we can write: :

* + * + * + * + * + * + * + 1 1 : 0 n n n n in in in in ij ij ij ij in kl h h h g g g g g - - • F ? F - F / F ? F (
n n n in in in ij ij ij in kl h h g g g g / - Ç • F ? / F -F È Ù F É Ú (14)
The value of * +

1 n in ij g - F
is then calculated by a Newton-Raphson iteration scheme.

As a result, if the relative error * + * + * +

1 1 n n in in ij ij tol n in ij g g g g - - F
/ F F , the solution is retained, else the calculation is repeated for the next iteration step n + 1. The integration algorithm is resumed in the following steps below: | -Assume an inelastic strain increment tensor ~ -Compute the volumic fraction of voids using equation (10) ¡ -Compute the variables of Haghi and Anand damage model using ( 7), [START_REF] Zavaliangos | Towards a Capability for Predicting the Formation of Defects During Bulk Deformation Processing[END_REF], and (9) ¢ -Compute the trial stress tensor using Hookes law. £ -Compute the effective stress triaxiality X ⁄ -Compute the internal state variables in the equations ( 1) and ( 2)

ƒ -Compute residual function ij h
a -If the residual function is less than the fixed tolerance then retain computed variables and stress as solution and go to next increment. b -Else update inelastic strain increment using (14) and go to step ~.

Numerical results

The exposed algorithm was implemented in the AbaqusAE FE code. Pure mechanical tests on Sn96Ag3Cu1 solder were used to prove the capability of the model. In this way, we first simulate uniaxial tensile and relaxation tests with many strain rate conditions. The simulations results shown in figure 1 and 2, proved the accuracy of the algorithm.

Conclusion

Porous viscoplasticity in SnAgCu solder is studied within a porous viscoplastic potential-based model. The model assumes an initial uniform distribution of voids in the whole volume and follows the phenomenological concept of internal state variables is able to describe hardening, creep and static thermal recovery among others. The implicit integration scheme proves to be stable and successfully implemented into finite element code AbaqusAE. Effect of nucleation is observed clearly after the first few cycles and increases under a nearly saturated level of about 6.×10 -2 .

  13)That allows determining the solution increment by transforming equation (13) into,

Fig. 1 .Fig. 2 .

 12 Figure1(a)shows the ability of the model to reproduce recovery effects. In the case of SnAgCu, kinematic hardening is dominant, and the saturation of stress is observable especially at high temperature. We denote a reduction of the final stress from 60 MPa at 35°C to nearly 30 MPa at 95°C. On figure1(b) relaxation tests at various temperatures demonstrate the effect of the thermal or static recovery on the stress response. In a second way, weve simulated cyclic tests to illustrate the damage behaviour coupled with the viscoplastic law. In the ten first cycles of the tensile loading, the mean volumic fraction of voids is constant so it means that there is no nucleation of voids but only regular voids growth due to !openclose" cycles (Figure2(a) stage A). This fact is improved also by taken into account voids shrinkage.


Nicole El Zoghby
email: nicole.el-zoghby@hds.utc.fr

Véronique Cherfaoui
email: veronique.cherfaoui@hds.utc.fr

Bertrand Ducourthial
email: bertrand.ducourthial@hds.utc.fr

Thierry Denoeux
email: thierry.denoeux@hds.utc.fr

Distributed Data fusion for detecting Sybil attacks in VANETs

come

Introduction

Exchanging data in a Mobile Ad hoc NETwork (MANET) in a safe manner becomes an important issue. These networks are vulnerable to different attacks such as intrusion. The need for security requires the introduction of the notion of confidence, as each node should have confidence in other nodes or in the received data before using the exchanged information in different applications. By broadcasting messages, nodes will discover their neighborhood. These neighbors can be fake or real nodes, they can also be attackers. Different research papers have been dedicated to find a solution to these problems. Many recent works deal with reputation mechanisms ( [START_REF] Zacharia | Trust management through reputation mechanisms[END_REF], [START_REF] Singh | An evidential model of distributed reputation management[END_REF], [START_REF] Liu | Enhanced reputation mechanism for mobile ad hoc networks[END_REF]) and trust evaluation ( [START_REF] Theodorakopoulos | Trust evaluation in ad-hoc networks[END_REF], [START_REF] Wang | A new evidential trust model for open communities[END_REF]) to manage the confidence in the source of information. Others were interested in data aggregation without taking into account the source [START_REF] Chen | Dempster-shafer theory for intrusion detection in ad hoc networks[END_REF][3][10] [START_REF] Raya | On data-centric trust establishment in ephemeral ad hoc networks[END_REF].

We propose a method to fuse data in a distributed system in order to build confidence over the network. Nodes broadcast their opinions, which are then used at the reception to evaluate other nodes. Since local opinion is uncertain and incomplete, the use of belief functions to evaluate the received messages seems appropriate. The fusion of a node's local knowledge with all the received messages is done by Demp-ster's rule. The network can suffer from cycles of data dissemination where the same information can be combined many times as it is coming from independent sources [START_REF] Evans | Bayesian network model for data incest in a distributed sensor network[END_REF], [START_REF] Mitchell | Multisensor Data Fusion: An introduction[END_REF]. To avoid that, we use the cautious rule of combination [START_REF] Denoeux | Conjunctive and disjunctive combination of belief functions induced by nondistinct bodies of evidence[END_REF].

We are interested in studying the confidence in a node for the purpose of detecting sybil attacks in VANETs (Vehicular Ad Hoc NETworks). The sybil attack is the case where a single faulty entity, called a malicious node, can present multiple identities [START_REF] Douceur | The sybil attack[END_REF] known as sybil nodes or fake nodes. This attack can affect the functionality of the network for the benefit of the attacker. Several techniques have been developed to detect misbehaving or fake nodes in VANETs. Gole et al [START_REF] Golle | Detecting and correcting malicious data in vanets[END_REF] represented an adversial parsimony that means finding the explanation for corrupted data. Vehicles can distinguish their neighbors by using cameras or exchanging messages in infrared light spectrum. The technique described by Xiao et al. is based on statistic signal strength analysis with the help of roadside infrastructure to detect sybil nodes [START_REF] Xiao | Detection and localization of sybil nodes in vanets[END_REF]. Yan et al. [START_REF] Yan | Providing vanet security through active position detection[END_REF] used an on-board radar to detect neighbors and to confirm their announced position. Piro et al [START_REF] Piro | Detecting the sybil attack in mobile ad hoc networks[END_REF] showed that the sybil attack can be detected passively through single or multiple observers. Due to the dynamics of the vehicular networks, of the number of vehicles and of the lack of permanent infrastructure access, deploying a Public Key Infrastructure in vehicular network (Vehicular PKI) is a very challenging task. As shown in [START_REF] Guette | On the sybil attack detection in vanet[END_REF], by simply comparing the received signal strength, half of the vehicles can detect the Sybil nodes and it is expected that cooperative techniques would decrease the number of cheated vehicles. Our work proposes such a cooperative algorithm between vehicles, based on the theory of belief functions, and could allow to avoid cryptographic schemes.

In this paper, we develop a distributed fusion technique based on the theory of belief functions. We first describe the system and how we represent the confidence using mass functions. We present the distributed data fusion approach and the proposed algorithm. We validate our approach by simulations and finally we conclude.

Distributed data fusion approach

We consider a network composed by nodes exchanging messages. It can be modeled by a directed graph G = (V, E), where V represents the set of nodes V = {v 1 , v 2 , ..., v n } and E represents the set of edges. The neighbors of each node are represented by Γ (v) = {v j ∈ V, {v i , v j } ∈ E}. For the sake of simplicity, we suppose that each node knows n = |V |. Figure 1 shows an example of network configuration. Each node periodically sends regular messages composed of its true identity and geographical position. Moreover, one of the node sends both its regular messages and fake messages composed of a forged identity and a forged position. By receiving the fake messages, other nodes are cheated and consider a non existing node, called fake node or Sybil node. We consider a single malicious node, which creates several Sybil nodes. All nodes use the same transmission system (same antenna, same transmission power). The topology of the network is given by the transmission radio range of the nodes (unit disk graph). We propose a data fusion methodology to combine data exchanged in a mobile ad hoc network, with the aim of quantifying the confidence in the nodes of the network. For this purpose, the messages contain also the sender's confidence in the nodes of the network.

Representing the confidence by mass function: Each node is able to assign a confidence to each other node of the network. This confidence is represented by a basic belief assignment (bba) denoted by m, defined from the frame of discernment Ω = {0, 1} where 0 represents FakeNode and 1 represents RealNode.

We denoted by m i j the corresponding bba that represents the opinion of node v i about node v j . The bba m i j is defined in Ω by:

m i j ( / 0) = 0 ; m i j (0) = p i j ; m i j (1) = q i j ; m i j (Ω ) = 1 -p i j -q i j . (1)
Principle of the approach: Node v k sends a message to v i containing its identity, its coordinates and its opinion about the network. When node v i receives the message, it calculates, after analyzing the signal strength, what we call a direct confidence. It is a mass vector denoted by m d ik . This direct confidence is saved in a local memory called local knowledge or private knowledge.

Note that each node has two bodies of knowledge: local knowledge and public knowledge. Local knowledge represents what each node can collect from its neighborhood. It is combined with the public knowledge of other nodes in order to update the public knowledge and rebroadcast it through the network. We thus have a distributed system. Local knowledge depends only on the signal strength of the messages and not on their content: consequently, it cannot be cheated. In contrast, public knowledge is based on the combination of the content of the messages and can be cheated by fake messages. This is why we separate local and public knowledge. The internal memory of each node is thus represented by two mass vectors (arrays of |V | cells initialized at m(Ω ) if i = j and m(1) if i = j):

K private i (t) = [m (t) l i j ] ; K public i (t) = [m (t) p i j ]. (2)
Distributed fusion Algorithm: The processing steps performed at the reception are presented in Algorithm 1, and explained hereafter.

Distributed Fusion: When node v i receives a message, it computes the direct confidence m d ik . This confidence is independent of previous messages and it is not the result of any other combination. So we use it to update the receiver's local knowledge about the transmitter by Dempster's rule [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF]. The function UpdateLocalKnowledge (m

(t-1) l ik , m (t)
d ik ) is calculated as: Algorithm 1: Received Message Processing on node v i Require: message from v k to v i , the signal strength P, message contains m p k j ∀ j Ensure:

K private i =[m (t) l i j ] and K public i =[m (t) p i j ] ∀ j ∈ V m (t) d ik ← DirectCon f idence(message, P) m (t) l ik ← U pdateLocalKnowledge(m (t-1) l ik , m (t) d ik ) m (t) p ik ← U pdatePublicKnowledge(m (t-1) p ik , m (t) l ik ) α ← DiscountingFactor(m l ik ) for each node j ∈ V such as j = i , j = k do α m (t) p k j ← DiscountTransmitterKnowledge(α, m pk j (t) , m (t)
Ω ) m (t)

p i j ← U pdatePublicKnowledge(m (t-1) p i j , α m (t) p k j ) m (t) l ik = m (t-1) l ik ⊕ m (t) d ik , (3)
where ⊕ denotes Dempster's rule. Since fake nodes might falsify the opinion of each node, the knowledge of other nodes is needed. To this end we use a distributed fusion to collect other opinions. As we consider that the transmitter is not totally reliable, we discount its opinion before combining it with the receiver's knowledge. The discounting factor α = 1m l ik (1) is defined as the plausibility that the transmitter is unreliable. The transmitter's opinion is discounted with the function DiscountTransmitterKnowledge(α, m

Ω ) as follows:

α m (t)

p k j = (1 -α).m (t) p k j + α.m (t) Ω . (4)
To update the receiver's public knowledge, we use the cautious rule [START_REF] Denoeux | Conjunctive and disjunctive combination of belief functions induced by nondistinct bodies of evidence[END_REF]. In a distributed system, the same information can be received and treated many times. While combining the knowledge, it is useful to use an idempotent rule to avoid counting the same information several times (data incest) as if it is provided by different independent sources. So the function UpdatePublicKnowledge(m (t-1)

p i j , α m (t) p k j )
allows us to combine the receiver's public knowledge with the transmitter's discounted knowledge about its neighbors as follows:

m (t) p i j = m (t-1) p i j ∧ α m (t) p k j , (5)
where ∧ denotes the cautious rule.

Direct confidence: Different methods can be used to compute the direct confidence m d ik . We propose a method that allows us to convert a real measure into a mass function. The real measure is based on signal strength analysis. Each receiver can analyze the signal strength to detect if the announced position is the real one [START_REF] Guette | On the sybil attack detection in vanet[END_REF]. It measures the strength of the received signal and calculates a theoretical value in terms of the node's coordinates. The estimated value of the signal strength is calculated by the Friis formula as µ = P e .G SR /d 2 ik , where

• P e is the transmitted signal power, depending on the transmitter antenna;

• G SR = G t .G r .λ 2 16.π 2 is the antenna gain, G t and G r are the gains of the transmit antenna and the receive antenna, respectively, and λ is the wavelength; • d ik is the distance between the transmitter node v k and the receiver node v i .

The comparison between the estimated power and the theoretical one allows the detection of a misbehavior. We propose to compute the plausibility that the received signal power P is equal to x, given that the transmitting node is a true node (ω = 1) as follows:

pl(P = x/ω = 1) = f (x/ω=1) sup x ∈R ( f (x /ω=1) , (6)
where f (x/ω = 1) is the normal density function with mean µ and standard deviation σ depending on the receiver antenna.

The plausibility pl(P = x/w = 0) is defined as shown in Figure 2: if the estimated and the theoretical powers are equal, we leave the possibility that the transmitter can be a fake node. Indeed, if the transmitter is a fake node but its position is near the malicious node therefore the estimated position will be approximately equal to the measured position. This result can influence the detection of the fake node. The direct confidence is computed using the Generalized Bayes theorem [START_REF] Ph | Belief functions: the disjunctive rule of combination and the generalized Bayesian theorem[END_REF]. It is obtained by combining the prior knowledge about the transmitter m Ω 0 with the plausibility that the node is a fake node knowing that it is a real {0} pl(x/w=1) and the plausibility that the node is a real node knowing that it is a fake {1} pl(x/w=0) :

m (t) d ik = m Ω (./x) = m Ω 0 ∩ {0} pl(x/w=1) ∩ {1} pl(x/w=0) , (7)
where ∩ denotes the unnormalized Dempster's rule.

Results

In order to validate our approach, Algorithm 1 has been implemented in Matlab. Simulations were performed on static and dynamic network. For simplicity of analysis, we first assumed all nodes in the network to be static. We performed simulations on different random network configurations. Next we tested our approach on a dynamic network, where nodes were moving in the same direction following a highway scenario. Implementation: In this part we will represent an example of a network composed from six true nodes, one of them is a malicious node that creates three fake nodes. The transmitted signal Power P e is about 600 mW and the antenna ranges is in order of 400 m. We consider that each transmitter sends its id, its position and its public knowledge. The receiver uses these informations to perform all the calculations and to verify if the node is true or fake. Simulations are performed until the convergence of the algorithm. We consider that the algorithm has converged when |m

(t-1) i j -m (t) i j | < ε,
where ε is a defined small threshold. The results of the simulation will be represented by gray scale matrices.

Static network: We present in Figure 3 an example of a network configuration where the nodes are static (left figure) and the result of the simulation (right figure). The white color in the right figure corresponds to a mass equal to 1 representing true nodes. The black color correspond to a mass equal to 0 representing fake nodes. The malicious node 3 will try to convince other nodes that the fake nodes [START_REF] Golle | Detecting and correcting malicious data in vanets[END_REF][START_REF] Guette | On the sybil attack detection in vanet[END_REF][START_REF] Liu | Enhanced reputation mechanism for mobile ad hoc networks[END_REF] are true nodes. The fake nodes have the same opinion as the malicious node. The first part of the rightmost figure represents the private knowledge. Each node has only information about its neighbors. The second part represents the public knowledge. We see that m p i j ({1}) = 0 for i = {1, 2, 4, 5, 6} and j = {7, 8, 9}, which means that the true nodes have correctly identified nodes {7, 8, 9} as untrustworthy. To verify the convergence of the algorithm, we performed simulations on different random network configurations by changing the number of the fakes nodes. Table 1 represents the result with different proportions of fake nodes. Each iteration represents the simulation of the process of a message. It needs more time to converge when the proportion of the fake nodes is greater. Our approach can detect sybil nodes with different static configurations.

Dynamic Network: Static configurations have some limits, especially when a malicious node is not in the neighborhood of the true nodes: in that case, fake nodes cannot be detected. So, we simulated a dynamic scenario in which nodes move in the same direction as on a highway. While moving, the neighborhood of each node changes. It influences the private knowledge because it depends on the neighborhood. Thanks to public knowledge, each node can get information about the whole network and can quantify its confidence. Table 2 shows results for different dynamic network configurations. The number of iterations until convergence changes at each simulation, because the node motions and neighborhoods are random. These preliminary results suggest that true nodes can successfully detect fake nodes in the network while moving on a highway.

Conclusion

A distributed data fusion approach based on belief functions for detecting sybil attacks in VANETs has been developed. The method uses both Dempster's rule and cautious rule to combine information and to compute a distributed confidence over the network. The results are promising and demonstrate that we can determine the reliability of nodes and detect fake nodes in a VANET. More realistic scenarios are currently being studied using an ad hoc network simulator. The method presented in this paper computes the confidence in the nodes without taking into account the contents of the messages exchanged in the network. The joint analysis of information and node reliability is currently being investigated. Results along these lines will be reported in future publications.

Fig. 1

 1 Fig. 1 Network Configuration

Fig. 2

 2 Fig. 2 Plausibility of received power values for true (ω = 1) and fake (ω = 0) nodes

Fig. 3

 3 Fig. 3 Network configuration and simulation results.

Table 1

 1 Results with different nodes configurations These results represent the average of 20 simulations.

	Nodes Configurations	Average of the number of iterations a Standard deviation
	True Nodes=6 Fake Nodes =3	207.05	7.86
	True Nodes=6 Fake Nodes =4	227.55	6.89
	True Nodes=6 Fake Nodes =5	255.8	6.33
	True Nodes=6 Fake Nodes =6	304.7	7.55
	a

Table 2

 2 Results for dynamic networks with different node configurations These results represent the average of 10 simulations.

	Nodes Configurations	Average of the number of iterations a Standard deviation
	True Nodes=6 Fake Nodes =3	119.3	45.88
	True Nodes=6 Fake Nodes =4	274.4	40.96
	True Nodes=6 Fake Nodes =5	361.1	54.23
	True Nodes=6 Fake Nodes =6	376.3	32.05
	a


Mickaël Garnier

Thomas Hurtut

Laurent Wendling

Object Description Based on Spatial Relations between Level-Sets

Object recognition methods usually rely on either structural or statistical description. These methods aim at describing different types of information such as the outer contour, the inner structure or texture effects. Comparing two objects then comes down to averaging different data representations which may be a tricky issue. In this paper, we introduce an object descriptor based on the spatial relations that structures object content. This descriptor integrates in a single homogeneous representation both shape information and relative spatial information about the object under consideration. We use this description in the context of image retrieval and show results on a butterfly image database compared with both GFD and dense SIFT descriptors. These results show that our method is more efficient to distinguish the objects where the spatial organization is a discriminative feature.

I. INTRODUCTION

Many computer vision applications rely on the automatic description of an object image. For instance, object recognition and image classification usually use features that endeavour to describe different types of information such as the outer contour, inner structure or texture effects. These different image informations often lead to different types of data which can be tricky to combine and may lead to inhomogeneous mixing of data [START_REF] Jolion | Feature similarity[END_REF].

A. Our Contribution

In this paper, we introduce a new image representation that capture these heterogeneous information in a single homogeneous representation. Given a decomposition of an object image into several disjoint layers of pixels, representing the different patterns presented by the object, the key idea of our method is to encode the pairwise spatial relations between all these layers. When applied to each layer with itself, the spatial self-relations encode first order shape information whereas for two different layers, the spatial relations encode relative structure and texture aspects. We show that a simple image decomposition such as the quantized level sets of the image gives interesting results, thus preventing from considering complex segmentation techniques.

B. Related Works 1) Spatial Relations:

A core aspect of our method is the encoding of the pairwise spatial relations. Literature in this domain can be structured in two main categories : qualitative and quantitative approaches. Qualitative approaches use symbolic relations such as positioning relations (left, right, below, above, etc.) and topological relations (inside, outside, etc.), see for instance [START_REF] Freeman | The modelling of spatial relations[END_REF] [START_REF] Egenhafer | Categorizing binary topological relations between regions, lines, and points in geographic databases[END_REF]. In this paper, we seek to capture a precise description of possibly complex objects and to characterize both large-scale and low-scale directional relations. Depending on the content meaning, the object patterns may also contained unconnected subsets of pixels. Therefore, in our context the spatial relations cannot be summarized in a symbolic manner.

Quantitative approaches gather methods that precisely describe the relative positions between two binary objects. Fuzzy quantitative methods are popular in different application domains such as spatial reasoning in medical images [START_REF] Bloch | Fuzzy spatial relationships for image processing and interpretation : a review[END_REF] and handwritten symbol recognition [START_REF] Delaye | Fuzzy relative positioning templates for symbol recognition[END_REF]. These methods produce a fuzzy landscape per considered potential direction. Combining these landscapes in order to capture the omnidirectional spatial organisation of possibly sparse object is not obvious. In this paper, we build on a quantitative model called force histogram [START_REF] Matsakis | A new way to represent the relative position between areal objects[END_REF], thereafter noted F-histograms. This model straightforwardly handles sparse objects and summarize their relative position in every directions in a single histogram. Basically, a F-histogram between two objects is a circular distribution measuring the relative attraction between these objects along every desired directions.

2) Object Recognition: Aside from spatial relations reasoning, object recognition is an attested issue with many approaches offering good results irrespective of the spatial relations of the object patterns. Among them, the Generic Fourier Descriptors [START_REF] Zhang | Shape based image retrieval using generic fourier descriptors[END_REF], hereafter named GFD, measure a shape descriptor. This descriptor belongs to the MPEG-7 standard and its key idea is to compute several polar Fourier transform of the image for several angular and radial frequencies, the descriptor being the normalized histogram composed by the values of these transforms.

Another classical method is the SIFT descriptor [START_REF] David | Object recognition from local scale-invariant features[END_REF] which consists in a 128-dimensional vector containing a set of gradient orientation histograms computed. The keypoints where histograms are computed are either extracted using the SIFT detection algorithm or consist in a matrix of points regularly spaced on the image (dense SIFT). Around every keypoints, a neighbourhood, divided in 4x4 smaller areas, is considered. On each areas histograms with 8 intervals are then computed. The final descriptor for each keypoint consists in a concatenation and a normalisation of the 16 histograms.

In this paper, we compare our proposed descriptor to these two generic and widely used methods. Our aim is to show Based on the set of S i , we then compute all F-histograms between every pairs (S i , S j ) for all 1 ≤ i ≤ j ≤ N inside a symmetric table (c). Each F-histogram in this table represents the description of the spatial relations between each considered pair of layers. Note that the diagonal part of the table represents self-relations, hence encoding shape information. The upper triangular of this symmetric structure embeds pairwise relations information between different layers.

the interest of our proposed model considering the inner spatial relations to efficiently recognize objects, compared to methods such as the GFD, that encodes exclusively shape information and as the SIFT, that models the organisation of local information.

C. Method Overview

Our goal in this paper is to demonstrate the descriptive strength of the inner pairwise spatial relations of an object content. We therefore assume that the objects we consider have been preprocessed using a background-foreground segmentation step. That is, the background pixels are not considered in our computations. As summarized in Figure 1b, our method first decomposes the object into N layers S i of its level sets. We then build a table formed by all the F-histograms computed between every pairs (S i , S j ), ∀i ∈ {1..N }, ∀j ∈ {i..N } (see Figure 1c). Consequently, this structure contains two types of information encoded in a single homogeneous representation. First the upper triangular encodes the pairwise spatial relations between every different layers. The diagonal encodes self spatial relations, which naturally models shape information of each considered layer.

II. F-HISTOGRAM DECOMPOSITION

In this section, in order to keep the paper self-contained, we first recall the definition of a F-histogram between two binary objects. Then, we will introduce the proposed object image representation based on a table of F-histograms.

A. F-Histograms

Originally, force histograms were introduced to solve the problem of measuring the fuzzy relative direction between two objects [START_REF] Matsakis | A new way to represent the relative position between areal objects[END_REF]. Basically, these are circular histograms measured along the directions θ ∈ [0, 2π[. Looking for the principal mode of this histogram is somehow equivalent to find the best θ that support the proposition "the first object is in direction θ from the second one".

To compute such a F-histogram, the objects are immersed in a space where an attractive force ϕ r operates. The definition of this force can vary widely, depending on the feature searched inside the objects. In order to have an intuitive representation of the involved spatial relations, this attractive force ϕ r is typically defined by a gravitational force based on the pairwise point distance:

∀(x, y) ∈ R 2 × R 2 , ϕ r (x, y) = 1 (d xy ) r (1)
where d xy is the Euclidean distance between two points x and y.

The F-histogram value along a direction θ between two objects A and B corresponds to the global force exerted by A with regard to B in the direction θ. In other words, this force F AB r (θ) is the integral sum of the infinitesimal forces ϕ r (a, b) where (a, b) ∈ A × B and the vector ab is along direction θ. Due to computational considerations, this global force is calculated on a set C θ of θ-oriented longitudinal cuts of the objects. Each of these longitudinal cut is built upon a straight line δ θ along direction θ, and is composed of two sets c A and c B of possibly disjoint segments (see Figure 2):

c A = δ θ ∩ A c B = δ θ ∩ B
The force exerted along δ θ by every points of c A with regards to every points of c B can be written as:

f θ (δ θ ) = c A c B ϕ r (a -b) db da (2)
The F-histogram value between A and B along the direction θ is then:

F AB r (θ) = δ θ ∈C θ f θ (δ θ ) = C θ δ θ ∩A δ θ ∩B ϕ r (a -b) db da (3)
Fig. 2. The global attractive force between A and B along the direction θ is the integral sum of infinitesimal forces computed on longitudinal cuts

(c A , c B ).
According to the value given to r, the resulting F-histogram can vary widely, giving more or less importance to closer objects. For r = 0 the F-histogram gives the same influence to distant organized structures from any local point. In this specific case, a histogram of forces is somehow equivalent to the previously proposed histogram of angles [START_REF] Miyajima | Spatial organization in 2D images[END_REF]. Yet practically, it has been shown that histograms of angles are not isotropic, not robust to rasterization and more computationally expensive [START_REF] Matsakis | A new way to represent the relative position between areal objects[END_REF]. Moreover, histogram of angles do not handle overlapping objects which will appear in our situation when computing the histogram between the same layers of a decomposed object. Note that F-histograms are naturally translation invariant, symmetric and isotropic [START_REF] Matsakis | A new way to represent the relative position between areal objects[END_REF].

B. F-Histogram Decomposition

Unlike the classical mode searching use of a F-histogram [START_REF] Matsakis | A new way to represent the relative position between areal objects[END_REF], we propose to consider the whole F-histogram as a complete signature of the relative positions between two possibly sparse objects, and therefore benefit from the whole information contained in the F-histogram.

In order to capture the inner spatial relations that structure the object, we first have to break this object up into multiple parts. Object segmentation is still an open research issue and no generic robust algorithm exists yet. Besides, our goal is different here since we aim at cutting out the object into its subparts. Therefore instead of using a complex segmentation algorithm, we choose a more pragmatic approach by decomposing a greyscale object image Q with scalar values Q(x) ∈ [0; 1[, into its N intensity level sets S i (Figure 1b). That is:

∀i ∈ [1; N ], S i = x ∈ R 2 , i -1 N ≤ Q(x) < 1 N (4)
This underlying assumption is also supported by the fact that mathematical morphology has shown that image contours locally coincide with level-set borders [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF], [START_REF] Cao | A theory of shape identification[END_REF]. By doing so, we thus adopt an oversegmentation, yet following the inner contours and whose behaviour will not vary from an image to another.

Based on this level-sets decomposition, we then compute all F-Histograms of every pair (S i , S j ). These F-histograms encode first-order shape information for each layer when i = j and second-order spatial relations information when i = j, both in the same mathematical formalism (Figure 1c). Considering the force ϕ r used during these computation, one may have a twofold strategy. First, the self (i = j) F-histograms along the table diagonal are computed with a force ϕ 0 . The natural overlapping induced by a level set with itself indeed lead to infinite forces when using ϕ r =0 . Dealing with the spatial relations information contained in (i = j) F-histograms, the ϕ 2 force (gravitational case) can be used. This choice is preferable in this specific case, since it models relative spatial relations, as suggested in [START_REF] Matsakis | A new way to represent the relative position between areal objects[END_REF].

The FHD of an object image Q is thus defined as:

F HD(Q) = F SiSi 0 ∀i∈{1..N } ∪ F SiSj 2 ∀(i,j)∈{1..N } 2 ,j>i (5)
This decomposition sums up N (N + 1)/2 F-histograms, made of N shape descriptors elements (diagonal), and N (N -1)/2 relative spatial relations elements (upper triangular).

Let notice that due to the invariant properties of Fhistograms, FHD are naturally translation invariant and symmetric. Depending on the application requirements, one can make FHD scale invariant by normalizing the histogram surfaces by the object surface. Finally, F-Histograms being isotropic [START_REF] Matsakis | A new way to represent the relative position between areal objects[END_REF], rotation invariant can be pursued by estimating the principal mode of the FHD or by minimizing the distance between globally shifted FHDs.

C. FHD Matching

In order to test the FHD descriptor on recognition and retrieval tasks, a dissimilarity measure is needed. The Fhistograms are first normalized by their surface in order to give an equal potential contribution to each F-histogram of the FHD. The distance between a query image Q and a target image T is then defined as:

D(Q, T ) = α × D shape (Q, T ) +(1 -α) × D spatial (Q, T ) (6)
where :

D shape (Q, T ) = 1 N N i=1 d χ 2 F SiSi 0 (Q), F SiSi 0 (T ) (7) D spatial (Q, T ) = 2 N (N -1) × N i=1 N j=i+1 d χ 2 F SiSj 2 (Q), F SiSj 2 (T ) (8)
where N is the number of layers of the FHD and α is the weight level given to the shape information compared to the spatial relations information. We use a chi-square distance to measure the distance between two single F-histograms, defined as:

d χ 2 (a, b) = imax i=1 (a(i) -b(i)) 2 a(i) + b(i) (9)

III. EXPERIMENTAL RESULTS

Several experimentations have been conducted using different sets of parameters for the FHD compared to two classical recognition methods, the Generic Fourier Descriptors (GFD) and the dense SIFT descriptor (dSIFT).

A. Image Database

Our experiments have been conducted on a database made from a subset of the Peale Collection [1]. This database is composed of 318 greyscale butterfly images grouped in 28 classes along the butterfly species. The typical height of an image is 640 pixels. All these images being over an homogeneous background, we first preprocessed the database by easily segmented the image backgrounds using a simple magic wand thresholding technique. Samples of this database are shown in Figure 3 and the subset we use is available online 1 . Butterflies are a typical case wherein inner spatial relations are a distinguishing feature making the wings patterns a direct link with the species. Another interesting application domain one can think about is botany taxonomy (flowers, mushrooms, . . . ).

B. Method Settings

In order to assess the descriptive strength of our Force Histogram Decomposition (FHD), we compare them wih the Generic Fourier Descriptors (GFD) and the Dense Scale Invariant Feature Transform (dSFIT). We present in this section all the settings we use for these methods in our experiments. The cross validation is made using a leave-one-out method, that is, for every object image, the remaining of the dataset serves as the training data. This method has been favoured due to the size of our database, leaving other approaches less statistically significant.

1 http://www.math-info.univ-paris5.fr/~mgarnier/dicta2012/ 1) FHD settings: The FHD are tested with several parameter sets, evaluating the gain from multiple forces histograms and the weighting between shape information and spatial relations parts of the feature. The robustness of F-histograms to the directions quantization has been studied in [START_REF] Matsakis | A new way to represent the relative position between areal objects[END_REF]. In order to avoid any binning effects, all the FHD are computed along 180 directions, regularly spanning the [0, 2π] interval with a 2 degrees step. In Equation ( 6), we test α values in {0, 0.2, 0.4, 0.5, 0.6, 0.8, 1}. Let notice that for α = 0 the shape descriptor is suppressed to promote the spatial relations and vice versa for α = 1.

2) GFD Settings: The GFD are based on the Polar Fourier defined as:

P F (ρ, θ) = x y I(x, y) × e [2jπ( r(x,y) R ρ+v(x,y)θ)] (10)
where r(x, y) and v(x, y) are respectively the radius and angle of the polar coordinates of the point (x, y), I is the intensity function and the parameters ρ and θ are bounded: 0 ≤ ρ < R and 0 ≤ θ < T with R and T respectively the the radial and angular resolutions. Finally, the GFD is written:

GF D(m, n) = |P F (0, 0)| M 11 , |P F (0, 1)| |P F (0, 0)| , ..., |P F (m, n)| |P F (0, 0)| (11)
where m and n are the radial and angular frequencies and M 11 is the order 1 moment. In our experiments, the GFD are computed on the object images with ρ = 4 and θ = 9, thus giving a signature of 37 bins, as suggested in [START_REF] Zhang | Image retrieval based on shape[END_REF].

3) dSIFT settings: The dSIFT are extracted with a step of 16 pixels and at several scales, 4 and 8 giving both local information and a more global one. The recognition is then made using a classical pair-wise image matching. For every points of the query image, the matching algorithm searches for the best matching point in the target image, if the resulting match gives a good contrast, this point vote goes to the target image. The finally matched image is the one with the higher votes. Although this pairwise image matching protocol is time consuming, it has been chosen due to its similarity with the matching used for the GFD and the FHD. In our experiments, we use the VL_Feat library SIFT implementation [START_REF] Vedaldi | VLFeat: An open and portable library of computer vision algorithms[END_REF].

C. Results and Discussion

Mean computational times for the processing of one image using N = 4 on an Intel CPU Xeon 3.0 GHz are the following. The FHD computation is approximately 2.3 seconds using a C programming implementation. Querying on the butterflies database using an unoptimised Matlab implementation takes around 5.2 seconds. Let notice that the overall complexity is O(N 2 ).

1) Recognition Rates: The recognition rates of the FHD are shown in Table I. The Table II shows the recognitions rates for thebest set of parameters {α = 0.8, N = 4}, as opposed to the GFD and the dSIFT performance rates. The rates in Table I approximately double from N = 2 to N ≤ 4 which show the interest of encoding inner layers information. N = 1 is indeed equivalent to only considering the whole object shape and thus gives results close to the GFD, see Table II. A second observation is that the values' increase tend to slow while N increases and start to decrease for N ≤ 8. This fact is first related to our database. Butterflies present homogeneous patterns with a rather limited range of luminance: black, white and one to three grey values, pleading to adopt an optimal value of N = 4 or N = 8. Secondly, the χ 2 distance compare F-histograms from the same combination of indices (i, j) for both Q and T objects. This histogram-tohistogram distance is thus neither illumination invariant nor contrast invariant. More sophisticated distances such as interhistograms distances should be investigated in order to tackle this issue.

The columns α = 0 and α = 1 in Table I show respectively the results with only the spatial relations descriptor and with only the shape description of every layer. These two sets of results point out the interest of extracting the information from both shape and spatial relations.

2) Precision-Recall Tests: Since we used a classified database, for each query, we can also compute the precisionrecall curve, classically defined as follows. Consider a query image belonging to a class of size C. For a given number W of images returned by this query, W tp is defined as the number of returned images belonging to the same class as the query (true positives). The precision-recall curve is then obtained by plotting the ratio P = W tp /W (precision rate) as a function of R = W tp /C (recall rate) [START_REF] Del Bimbo | Visual information retrieval[END_REF]. Averaged precision-recall curves over the whole database using the best parameters of the different compared methods are shown in Figure 4.

3) Qualitative Retrieval Results: Several comparative retrieval results are shown in Figure 5. These results show that the FHD are more efficient to distinguish the butterflies where the spatial organization is a discriminative feature, compared to the dSIFT and the GFD.

The GFD focuses on the global shape and is thus not able to correctly discriminate two butterflies having the same global contour, yet being from different species due to different wing patterns (see for instance second column, third row of Figure 5). On the first query of the same Figure, and contrary to the two other methods, the FHD manages to capture the spatial organisation of the sparse black stripes.

The dSIFT are actually not meant to be used with a classical matching but rather with a learning and a bag of word for instance. In this experiments, numerous matches are made between keypoints in homogeneous areas with very low distances overwhelming relevant matches thus giving descriptors with a very low variance all over the database hence some of the poor visual results. Despite this, the dense SIFT have been prefered to the classical SIFT keypoints detection. Since our database present objects after a background segmentation, most of the keypoints using classical SIFT remain located on the contrasted outer border. Such an approach lead to results similar to the GFD.

4) FHD Noise Robustness: The robustness of our descriptor is also evaluated on the same image database altered with noise. To do this, we choose the set of parameters giving the best recognition results and apply them with the same protocol but on a noised version of the database. We used two different kind of noise, that are speckle and Gaussian noise, with increasing variances, as shown in Figure 6. Speckle noise tests assess the robustness of the F-histograms to possibly low quality parsing of the level sets. It also shows that the FHD are more sensitive to Gaussian noise, as the recognition rates drops along with the variance increase.

IV. CONCLUSION

We proposed in this paper a feature based on spatial relations for grayscale object recognition. Based on a homogeneous stack of F-histograms, it naturally embeds both absolute and relative spatial informations about the considered object. It thus encodes information both on outer and inner contours and information on the spatial relations that organize the underlying grey levels of the object. We showed that a simple level-sets quantization is sufficient to capture enough information to discriminate highly structured images. Recognition rates on the butterfly database altered with either Gaussian or speckle noise.

The limitations of our approach are twofold. First we presently do not encore the color content. Since this information might be of prime importance in several applications, this future work is one of our very next goals. By using bin-to-bin distances, another aspect that has been eluded is the circularity of F-histograms. Circular distances such as CEMD [START_REF] Rabin | Circular Earth Mover's Distance for the comparison of local features[END_REF] could be investigated here although normalization constraints raise several issues.

Fig. 1 .

 1 Fig. 1. Overview of the proposed descriptor. A query image (a) is first decomposed into its N grey level sets S i (b), N = 3 on this Figure.Based on the set of S i , we then compute all F-histograms between every pairs (S i , S j ) for all 1 ≤ i ≤ j ≤ N inside a symmetric table (c). Each F-histogram in this table represents the description of the spatial relations between each considered pair of layers. Note that the diagonal part of the table represents self-relations, hence encoding shape information. The upper triangular of this symmetric structure embeds pairwise relations information between different layers.

Fig. 3 .

 3 Fig. 3. Sample images from our database (reproduced with the permission of The Academy of Natural Sciences, Philadelphia) [1]

Fig. 4 .

 4 Fig.4. Precision-recall curves over the whole database using the best sets of parameters for the FHD, the GFD and the dSIFT descriptors.

Fig. 5 .

 5 Fig.5. Eight retrieval results on a 318 butterflies public database, obtained with the different compared descriptors: the proposed F-Histogram Decomposition (FHD), the Generic Fourier Descriptor (GFD) and the dense SIFT descriptor (dSIFT). The FHD well succeeds to capture both the inner relative spatial organisation of the grey levels and absolute shape of the different patterns composing each butterfly.

Fig. 6 .

 6 Fig. 6.

TABLE I RECOGNITION

 I % FOR THE WHOLE FHD WITH THE MIXED FORCES ϕ 0 AND ϕ 2 .

	N \α	0	0.2	0.4	0.5	0.6	0.8	1
	1	3.46 29.9 29.9 29.9 29.9 29.9 29.9
	2	12.6 32.7 36.2 37.4 38.1 39.6 41.5
	4	45.9 53.5 53.5 53.1 53.8 56.3 44.3
	8	49.7 52.2 52.8 53.1 52.2 49.7 41.8
	16	41.8 45.3 46.9 47.5 48.4 51.9 47.2

TABLE II RECOGNITION

 II % FOR THE FHD, THE GFD AND THE DSIFT.

	descriptor	FHD	GFD dSIFT
	recognition %	56.3	28.6	43.4

Acknowledgments This work has been sponsored by the ANR SPIRIT #11-JCJC-008-01.