Nicolas Coudray
email: nicolas.coudray@uha.fr

Argyro Karathanou

Sylvie Chambon
email: chambon@lcpc.fr

MULTI-RESOLUTION APPROACH FOR FINE STRUCTURE EXTRACTION Application and validation on road images

Keywords: Segmentation, multi-resolution, road images, cracks

In the context of fine structure extraction, this paper presents a new method based on multi-resolution segmentation applied for the detection of road cracks. A method already developed to detected low-contrasted biological membranes has been adapted to detect cracks on images: crack features are defined as heterogeneities rather than transitions of closed regions characterizing the membranes. This new methodology is quantitatively validated on reference segmentations and compared to an adapted filtering and Markovian modelling algorithm.

INTRODUCTION

Detecting fine structures is very helpful in a lot of domains: to extract ceramic damages [START_REF] Elbehiery | Surface defects detection for ceramic tiles using image processing and morphological techniques[END_REF], to find cracks in underground pipes [START_REF] Iyer | A robust approach for automatic detection and segmentation of cracks in underground pipeline images[END_REF], to detect road network in satellite images [START_REF] Geman | An active testing model for tracking roads in satellite images[END_REF], to follow vessels in medical images [START_REF] Frangi | Muliscale vessel enhancement filtering[END_REF]. Since 1990, a lot of algorithms have been proposed in the domain of crack detection on road pavement surface. In fact, every country needs to evaluate, periodically, the quality of roads and most of this work is done manually, which is expensive, non reproducible, dangerous and not very efficient. In consequence, a lot of efforts have been made in the field of research on automatic or semi-automatic procedures for detecting deterioration on the roads and, in particular detecting cracks. In 2003, the report of Schmidt [START_REF] Schmidt | Automated pavement cracking assessment equipment -state of the art[END_REF] gives a good summary of the technologies developed in this field. The main difficulty of this task is that the default is not well contrasted compared to the texture of the road, and that it represents a very small part of the images (about 1.5% of the image).

The most recent algorithms proposed to tackle this problematic are based on multi-resolution approaches. Even if they are the most efficient, they have shown their limits by giving a sparse detection. In the field of detection of biological membranes, recent methods based on multi-resolution watershed segmentation provided great performances, even if the objects of interest are low-contrasted (Coudray et al., 2007). By analogy, the problem of biological membranes detection is similar to crack detection. The goal of this paper is to adapt this watershed method to road crack detection, and compare its performances with a Markovian modeling-based algorithm developed specifically for this purpose [START_REF] Chambon | Road crack extraction with adapted filtering and markov model-based segmentation[END_REF].

Firstly, a brief state of the art of fine structure extraction is exposed. Secondly, the studied crack algorithms are presented: the Markovian modeling algorithm and the multi-resolution segmentation. Thirdly, the protocol used to analyze the performances is given. And finally, before concluding, results are summarized and discussed.

FINE STRUCTURE EXTRACTION

We distinguish four methods: those based on a THRESHOLD [START_REF] Koutsopoulos | Primitive-based classification of pavement cracking images[END_REF], methods combining a threshold with mathematical morphology tools (MORPHOLOGY) [START_REF] Tanaka | A crack detection method in road surface images using morphology[END_REF], NEURAL NETWORK-based methods [START_REF] Bray | A neural nework based technique for automatic classification of road cracks[END_REF], and MULTI-SCALE filter-ing [START_REF] Subirats | Automation of pavement surface crack detection with a matched filtering to define the mother wavelet function used[END_REF]. A detailed state of the art can be found in [START_REF] Chambon | Road crack extraction with adapted filtering and markov model-based segmentation[END_REF]. THRESHOLD methods are popular and simple but inefficient (the amount of false positives and of false negatives is important). Methods based on MOR-PHOLOGY reduce the number of false positives but they are strongly dependent on the choice of the parameters. NEURAL NETWORK techniques are efficient but they need a learning step that is quite expensive and non-automatic. Finally, MULTI-SCALE methods seem to be the most efficient: unlike other methods, they lead to less false positives, but more false negatives. In conclusion, we can use methods that detect most of the crack but give a lot of errors or methods that obtain only a part of the crack with low errors.

In the field of microscopy images, a watershed segmentation has been proposed to extract biological membranes (Coudray et al., 2007). The problem of this task is comparable to the problem of road crack detection: edges of the membranes (here the default) are low and heterogeneously contrasted compared to the background (the road). Consequently, we propose to adapt this method to the road application.

CRACK DETECTION ALGORITHMS

In this part, we briefly present a recent crack detection algorithm based on an adapted filter combined with Markovian modeling segmentation and, then, we introduce a new method, based on multi-scale analysis with watershed and we describe all the contributions in order to adapt this method to a road crack detection.

Adapted Filtering and Markovian Modeling (AFMM)

The method proposed in [START_REF] Chambon | Road crack extraction with adapted filtering and markov model-based segmentation[END_REF]) uses these two hypotheses, well known in the field of road crack detection: (H 1 ) Pixels of a crack are darker than the background (the road); (H 2 ) A crack is a set of connected segments with different orientations. Consequently, the principle of the method is, first, to binarize the image (by adapted filtering) and, second, to refine the binarization by a segmentation based on Markovian modeling. At the beginning, we have to choose the number of scales, which depends on the resolution of the image (the size of the pixels): if we suppose that we have 1 mm/pixel, 5 scales are sufficient (we can detect cracks from 2 mm to 1 cm).

Then, we also have to choose the number of directions for the adapted filtering and, considering the aspect of the road cracks, it is realistic to use these directions:

[0, π 4 , π 2 , 3π 4 ].
Interested readers can find more details about this work in [START_REF] Chambon | Road crack extraction with adapted filtering and markov model-based segmentation[END_REF]. The goal of this paper is to compare this method with the new one presented below.

Multi-Resolution Detection of Heterogeneities (MRDH)

This method is inspired from the 3-step segmentation method developed for the segmentation of low and heterogeneously contrasted biological membranes (Coudray et al., 2007). The aim was to partition the image by detecting the closed contours of membranes. In this initial algorithm, the edges are firstly thresholded at different scales, and, secondly, the resulting binary images are combined in a Reconstructed Gradient-Like (RGL) image. Finally, the watershed algorithm is applied on the RGL image to obtain the searched partition. For the cracks, the third step has been adapted and a fourth step has been added to discard spurious detections. The steps are detailed below.

A multi-resolution transform is first applied to the initial image I(x, y, 1) using a pyramidal transform: images I(x, y, s) are obtained at scale s using:

I(x, y, s) = +k ∑ m=-k +k ∑ n=-k G(m, n, s)I(sx + m, sy + n, 1), (1)
where G(m, n, s) is the smoothing average filter. Dyadic pyramids reduce the size of the image by a factor of 2 between each scale. Here, non-dyadic pyramids are used to better represent and identify the features, with s ∈ V s = {1, 2, . . . , 10} (at scales coarser than 10, the images become too small to be analyzed). Heterogeneities are detected and thresholded at each scale, leading to binary images B(x, y, s). In this study, features to be identified, i.e. the cracks, are better described as local heterogeneities than local transitions (e.g. Figure 1), making a local standard deviation filter more appropriate than a gradient filter like Sobel, for instance. In the next step, binary images are combined to form the RGL image I RGL (x, y):

I RGL (x, y) = max s∈V s ((11 -s) × B(x, y, s)) .
(2)

With the (11s) weighting, the gray-level on the RGL image corresponds to the finer scale at which the feature has been identified. Details identified with a better precision have therefore a higher value (e.g. Figure 2). The next step aims to thin the RGL image features to obtain 1-pixel wide objects, and has been adapted from our initial work to this road application. It needed to be modified since, contrary to membrane edges, cracks do no form closed regions: therefore, we can either apply the watershed algorithm after considering the background of the RGL image as seeds (i.e. pixels which value is 0 and which touch non -0 value pixels), or by applying a gray-scale thinning algorithm [START_REF] Redding | The autoscaling of oblique ionograms[END_REF]. For its simplicity and lower computational time, the watershed approach has been used in this study. The resulting image I pc (e.g. Figure 3) points to the potential crack pixels (PCP). A last step has been added to validate I pc and reduce the false detections: a PCP is validated if its contrast is above a given threshold C Tr . To make the measure more robust to the noise, instead of analyzing each PCP independently, the contrast is averaged on segments (a segment is a set of contiguous pixels having at most two neighboring PCP, see Figure 4 for an illustration of this vocabulary). The average contrast C seg of a segment is the difference between the average gray-level of the segment (GL seg ), and the average gray-level of the neighboring background (GL Bkg ):

C seg = GL Bkg -GL seg .
(3) GL Bkg (e.g. Figure 5) is measured using a large filter averaging the neighboring values of the background (pixels not labeled as PCP). The filter kernel size did not seem critical, but it should be large enough to reduce the influence of noise, and small enough to consider the irregularities of the road: a 25 × 25 pixels kernel was experimentally chosen. A segment is considered as targeting a crack if C seg > C Tr (e.g. Figure 6). Measures have been realized to set the main parameters of the algorithms, i.e. the kernel size of the standard deviation filter and the threshold C Tr ; we identified that the best results were obtained with a 5 × 5 pixels kernel and C Tr = 50. We will describe the analysis methodology, and we will restrain the presentation of the influence of C Tr on the results.

ASSESSMENT METHODOLOGY

To evaluate and compare the crack detection methods, we had to choose: a/ the tested images, b/ how to determine the "ground truth" segmentation or reference segmentation and c/ the criteria in order to quantitatively evaluate the results.

Images

The algorithms were evaluated on 3 sets of test images: on 14 synthetic crack images, on 10 real images, and on 10 pre-processed real images. For synthetic images, it is easy and reliable to give a ground truth segmentation. For real images, estimate a "ground truth" or a reference segmentation is more complicated but the images are more realistic than the synthetic images.

The 14 synthetic images were built using different kinds of backgrounds: cracks were added on 8 images built with a random dot texture (sample S 1 ), and on 6 images of road containing no cracks, 2 of them acquired with a static camera (sample S 2 ), and 4 of them dynamically acquired with a camera embedded on a vehicle (sample S 3 ). For the last four ones, controlled lights were added. For all the 14 images, the cracks were randomly added, with a random shape and a random gray-level (Figure 7).

On the 10 real images, 4 were acquired using a static camera and 6 using the dynamic system (Figure 7).

In the third set of test images, the 10 real images were pre-processed using:

1. Threshold -In order to reduce the light halo in some images (the last six ones presented in Figure 9), each pixel over a given threshold is replaced by the local average gray levels.

2. Smoothing -A mean filter of size 3×3 is applied.

3. Erosion -An erosion with a square structuring element of size 3 × 3 is applied.

4.

Restoration -This last pre-processing tries to combine the advantages of all the previous methods in three steps: histogram equalization, thresholding (like Threshold), and erosion (like Erosion).

Reference Segmentation

For real images, we briefly explain how the manual segmentation is validated. Four experts manually segmented the images with the same tools and in the same conditions. Then, the four segmentations were merged, following these rules:

1. A pixel labeled as crack by almost two experts was kept;

2. Every pixel near to a pixel kept by step 1 was also kept.

The second rule is iterative and stops when no pixel is added.

Synthetic image Ground truth

Real image + simulated default Ground truth

Real image manually segmented Reference Then, the result is dilated with a structuring element of size 3 × 3. Results of the manual segmentations are presented in Figure 9 (for more details, the study of the reliability of these reference segmentations can be found in [START_REF] Chambon | Road crack extraction with adapted filtering and markov model-based segmentation[END_REF]).

Quantitative Assessment

In Figure 8, the evaluation criteria are presented. They are the most used in the literature, except "accepted" pixels. In consequence, we have included accepted pixels in the computation of the similarity coefficient. For estimating accepted pixels, a threshold is applied on the distance between the detected pixel and the reference segmentation. This threshold is estimated using the mean distance between each of the four segmentations used for establishing the reference.

The DICE values lie on the [0; 1] interval, but to simplify the visualization of the results, we display DICE×100 in all figures and tables of section 5. the edges, while the AFMM method is more adapted to keep the proportion of FN low. By comparing DICE values, the multi-resolution method appears to be the most adapted. Table 1 details the results for the different types of images (samples S 1 to S 3 presented in section 4.1). On fully synthetic images, both algorithms perform very well. It appears that the images acquired with the dynamic system are the most difficult to process, whatever the method used.

RESULTS AND DISCUSSION

Synthetic Images

Real Images

Figure 11 displays the results obtained on the real images. We see that the amount of FN is higher than with synthetic images, but we note that these values are probably over-estimated, the human-reference images tending to select only the major cracks. The AFMM algorithm identified almost half of the cracks selected by the human, and a bit more than half were identified with the multi-resolution method. Considering the DICE measure, the multi-resolution tests gave the best overall results; though, on camera acquired images, DICE values can be considered as equivalent for both methods.

Real Pre-Processed Images

Here, we compare the AFMM method with the MRDH method (only the one with the more efficient parameter, i.e. C Tr = 50, is discussed here) when the real images are pre-processed. In particular, we evaluate the influence of the 4 pre-processing described in section 4.1. Table 2 shows the DICE measures obtained. It can be seen that the four pre-processings enhance the overall performances.

When pre-processed with the threshold, smoothing, erosion and restoration, the AFMM's algorithm is better, with the best results being obtained after a smoothing. Also, smoothing and threshold have an even more beneficial effect when static images are used.

According to the DICE measure, the performances of MRDH applied after pre-treatment are quite disparate. In fact, it illustrates how this new method is globally quite robust against the acquisition conditions. Images do not need to be pre-processed in order to increase the performances and this is an important superiority compared to the AFMM method which is highly dependent on the acquisition conditions. Analyzed in details, we notice that for static images, DICE values seem to be enhanced by each pre-processing, while DICE values seem to decrease for "dynamic" images.

CONCLUSIONS

A methodology for adapting a multi-resolution-based segmentation to the field of road crack detection has been introduced. The experimental results demonstrate the efficient results obtained on 14 synthetic images and 10 significant real images. Moreover, it shows how it outperforms a previous method based on adapted filtering, in particular with the most difficult images acquired dynamically and that present illumination defaults.

The next step of this work will be to validate this work on more data (32 more images will be available soon). It has been shown that MRDH improves the percentage of true detections (this percentage is higher than with AFMM), but the AFMM presents less false negatives. The next improvement will be to reduce this phenomenon by introducing more constraints on the final step, like for example, an active contour constraint instead of using only the constraint on the contrast.

Figure 1 :

 1 Figure 1: Initial image example (I(x, y, 1)) -The crack corresponds to a local low-contrasted heterogeneity.

Figure 2 :

 2 Figure 2: RGL image of Figure 1 (I RGL (x, y)) -The brighter the pixels, the higher the precision of the detection.

Figure 3 :

 3 Figure 3: Potential crack pixels from Figure 2 (I pc displayed in blue on the original image).

Figure 4 :

 4 Figure 4: A potential crack composed of 5 segments (each having a different color).

Figure 5 :

 5 Figure 5: Estimated background of Figure 1.

Figure 6 :

 6 Figure 6: Final crack detected on Figure 1 (using C Tr = 40).

Figure 7 :

 7 Figure 7: Tested images.

Figure 10 Figure 8 :Figure 9 :

 1089 Figure 10 compares results for the different methods with synthetic images. It shows that the multiresolution with C Tr = 50 is the best to identify most of

Figure 10 :

 10 Figure 10: Summary of the algorithm performances on synthetic images.

Figure 11 :

 11 Figure 11: Summary of the algorithm performances on real images.

Table 1 :

 1 Results on Synthetic Images -The 3 image categories are: S 1 , fully synthetic images, S 2 , background acquired with a camera and S 3 , background acquired with the dynamic system. Results in bold indicate best values and it illustrates how the MRDH approach outperforms the AFMM, in particular, with the most difficult images (sample S 3 ) where DICE is correct (greater than 50) whereas it is very low with the other methods. S 2 S 3 S 1 S 2 S 3 S 1 S 2 S 3 S 1 S 2 S 3

	AFMM MRDH,	MRDH,	MRDH,
	C Tr = 30	C Tr = 40	C Tr = 50
	S 1 TP+Acc 87 46 15 68 41 30 86 53 43 91 66 59
	FN 1.5 3.7 2.6 2.8 4.4 17 5.3 6.9 18 6.7 6.9 18
	DICE 92 63 24 81 59 44 92 69 57 94 80 70

Table 2

 2

	: DICE of pre-treated images: the three numbers
	show the average values: with the 10 images (with the 4
	images acquired in static way, with the 6 images acquired
	in a dynamic manner).
		AFMM	MRDH, C Tr = 50
	None	46 (49,44)	55 (52,57)
	Threshold	64 (68,61)	55 (60,50)
	Smoothing 66 (77,59)	52 (63,42)
	Erosion	59 (58,60)	57 (64,51)
	Restoration 52 (49,55)	57 (59,55)


Romuald Jolivot

Hermawan Nugroho

Pierre Vabres

M H Ahmad Fadzil

Franck Marzani

Validation of a 2D multispectral camera: application to dermatology/cosmetology on a population covering five skin phototypes

Keywords: Skin imaging system, multispectral image, hyperspectral cube, spectral reflectance, spectral reconstruction, validation

This paper presents the validation of a new multispectral camera specifically developed for dermatological application based on healthy participants from five different Skin PhotoTypes (SPT). The multispectral system provides images of the skin reflectance at different spectral bands, coupled with a neural network-based algorithm that reconstructs a hyperspectral cube of cutaneous data from a multispectral image. The flexibility of neural network based algorithm allows reconstruction at different wave ranges. The hyperspectral cube provides both high spectral and spatial information. The study population involves 150 healthy participants. The participants are classified based on their skin phototype according to the Fitzpatrick Scale and population covers five of the six types. The acquisition of a participant is performed at three body locations: two skin areas exposed to the sun (hand, face) and one area non exposed to the sun (lower back) and each is reconstructed at 3 different wave ranges. The validation is performed by comparing data acquired from a commercial spectrophotometer with the reconstructed spectrum obtained from averaging the hyperspectral cube. The comparison is calculated between 430 to 740 nm due to the limit of the spectrophotometer used. The results reveal that the multispectral camera is able to reconstruct hyperspectral cube with a goodness of fit coefficient superior to 0,997 for the average of all SPT for each location. The study reveals that the multispectral camera provides accurate reconstruction of hyperspectral cube which can be used for analysis of skin reflectance spectrum.

INTRODUCTION

Skin lesion evaluations are mainly based on naked eyes/RGB camera analysis. However such analysis are limited due to poor discrimination of the human eye. It is linked to the fact that colour sensed by the human eye is a combination of three different bands called trichromatic model. Intensive researches in the non-invasive skin imaging area bring new imaging techniques that work as an adjunct help system and for some even surpass the human eye. The motivation is to go beyond the limited red, green, blue color imaging as it reflects the human light perception. The human eyes and RGB cameras are affected by a phenomenon called metamerism: two colour samples match under the same illuminant but may appear totally different when using a different illuminant; similarly, for skin lesions observed under different condition, the analysis made by a dermatologist can be affected by this phenomenon. Imaging techniques offer capacity to get extra information about the skin to the dermatologist for diagnosis assessment. Those efforts try to bring objective measurement of skin lesions. It can be done either by a spectroscopic measurement [START_REF] Tominaga | Principal component analysis-based reflectance analysis/synthesis of cosmetic foundation[END_REF][START_REF] Prince | Analysis of diffuse reflectance spectra of various skin conditions by principal component method[END_REF] or by imaging techniques [START_REF] Nugroho | Determination of skin repigmentation progression[END_REF][START_REF] Cotton | A noninvasive skin imaging system[END_REF] . Combining advantages of both spectrophotometer (spectral resolution) and digital camera (spatial resolution), MultiSpectral Imaging (MSI) systems overcome their respective limitations (lack of spatial variation and lack of spectral information).

Further author information: (Send correspondence to Romuald Jolivot) : E-mail: romuald.jolivot@u-bourgogne.fr, Telephone: +33 (0)3 80 39 59 34

The aim of this work is to prove the accurate reconstruction of hyperspectral cube of cutaneous data from multispectral images. In order to validate this, a comparison study between a commercial spectrophotometer and the ASCLEPIOS 5 camera developed in our team has been performed on a study population covering different skin tones.

MATERIALS

The acquisition device used in the study is a multispectral camera covering visible and near infra red spectrum, called ASCLEPIOS (see figure 1). It is an extended version of a previous multispectral camera limited to the visible [START_REF] Jolivot | Reconstruction of hyperspectral cutaneous data from an artificial neural network-based multispectral imaging system[END_REF] . The MSI is composed of a Xenon light source covering wavelength range from 340 nm to infrared light. The light compartment houses a rotating wheel holding a set of interference filters which gradually rotates to position each filter in front of the light source. The spectral band of the set of ten interference filters covers from 430 nm to 970 nm. An opaque custom built hand-held device encloses a liquid light guide, which transmits the illumination, a camera and a remote wireless trigger. The camera consists of a monochromatic CMOS digital camera with spectral response from 370 to 1100 nm and resolution of 1312 x 1082 pixels. The hand-held device, in contact with the skin, keeps a constant focus distance between the camera and the skin area. A magnifying lens is mounted in front of the camera which provides an useful area of 32 x 38 mm size with a depth of field of 5 mm, yielding a spatial resolution of 33 pixels.mm -1 . The acquisition is performed in less that 2 seconds by illuminating the skin at different wavebands after the automatic rotation of the wheel for each filter. The system acquires one multispectral image composed of ten monobands. The originality of the ASCLEPIOS system is its capacity to reconstruct a Hyperspectral Cube from a multispectral image. This reconstruction is performed by an Artificial Neural Network-based (ANN) algorithm using hetero-associative memory. [START_REF] Mansouri | Neural networks in two cascade algorithms for spectral reflectance reconstruction[END_REF] The reconstruction aims to retrieve skin reflectance spectra in each pixel from the camera signal. The reconstructed hyperspectral cube of cutaneous data provides a 3-dimensional volume (x,y,z) where x and y are for spatial dimensions and z for spectral one. The ANN is divided into two parts, first a learning process based on GretagMacBeth ColorChecker and then a reconstruction process. Due to the use of hetero-associative memory, the algorithm is able to reconstruct hyperspectral cube with different wave ranges within the system capability (430 to 970 nm) and step up to 1 nm..

METHODS

The goal of hyperspectral cube reconstruction is to retrieve skin information from the reflectance spectrum leading to an increase of the amount of information available to the dermatologist for diagnosis purpose. To validate the accuracy of the reconstruction of hyperspectral cube of cutaneous data, a study on healthy participants from different skin tones has been performed. The classification of the different skin tones is based on the Fitzpatrick Skin PhotoType 8 (SPT) which classify a person's complexion and tolerance to sunlight. According SPT Type SPT II SPT III SPT IV SPT V SPT VI Nb 18 21 56 41 14

Table 1. SPT distribution to the Fitzpatrick scale, human skin can be divided into six categories, SPT I for very white skin (never tans and always burns) to SPT VI for black skin (never burns). The different SPT aims to show accurate reconstruction for a wide range of spectrum. Three acquisitions (by a commercial spectrophotometer CS 2600d and by the multispectral camera) were collected from each participant of the study: two data were acquired from facultative skin colour and one from constitutive skin colour. Facultative skin colour is skin area affected by the environment (sun, hormones) and it impacts the original appearance of the skin colour over time whereas constitutive skin colour is determined by genetics. For facultative colour, one sample from the back of the hand and one sample from the face are acquired. For constitutive skin colour, one sample from the lower back is acquired. The purpose of different body locations is to record a difference of skin chromophore composition between these samples. This is due to the fact that the face and the back of the hand are supposed to be more exposed to the sun than the lower back. Melanin production is linked to sun exposure (UV), these two body locations are expected to have higher melanin concentration and therefore they should have a different skin reflectance spectrum. The aim is to observe a difference within each participant between exposed (facultative skin colour) and non-exposed (constitutive skin colour) area.

In this study, using the wave-range flexibility offered by the ANN, three reconstruction configurations are selected, one reconstruction from 430 to 780 nm, another one from 430 to 850 nm and a third one from 430 to 950 nm with a step of 10 nm. The objective of testing different wave-range is to select the best configuration for accurate reconstruction for visible light range.

The study involves 150 recruited healthy (with no known skin disease) participants of different skin tones. The acquisition was conducted at University Teknologi PETRONAS (Malaysia) on international and local undergraduate and postgraduate students. The population covers five out of the six SPT groups, no participant of type I were available for the study. The distribution within the different SPT is shown in table 1. The average age of the study population is 22.6 years with a standard deviation of 4.2 years.

To evaluate the performance of the reconstruction of hyperspectral cube, three different metric scales are employed, the Root Mean Squared Error (RMSE), the Goodness of Fit Coefficient (GFC) and the Reconstruction Percentage (RecP). These metric scale use the spectrophotometer data as reference.

The RMSE calculates the mean of euclidean distance values between two spectra. Reconstruction is good when RMSE is close to zero. RMSE is calculated using the following formula:

RM SE = 1 N N j=1 R m (λ j ) -R r (λ j ) 2 (1)
where R m (λ j ) is the spectrum value measured with the spectrophotometer at wavelength λ j and R r (λ j ) is the reconstructed spectra value at wavelength λ j and N the number of samples.

The GFC is based on the Schwartz inequality and is calculated using the following formula:

GF C = j R m (λ j )R r (λ j ) j [R m (λ j )] 2 1/2 j [R r (λ j )] 2 1/2 . ( 2
)
This criterion has the advantage to be bounded between 0 and 1, providing an easy interpretation.

The RecP is a criterion which evaluates the distance between two spectra. It requires a higher precision for the spectrum part having low amplitude. A good reconstruction will be close to 100%. It is defined by:

RecP = 1 - N j=1 R m (λ j ) -R r (λ j ) 2 N j=1 R r (λ j ) 2 .
(3)

RESULTS

The data obtained from the spectrophotometer is a spectrum from 360 nm to 740 nm with a step of 10 nm. The hyperspectral cube is reconstructed using the three configuration previously mentioned (430-780 nm, 430-850 nm, 430-950 nm) all with a step of 10 nm.

The results are divided into six groups (see table 2). The first five sets of rows concern the average of each SPT according to their locations and reconstruction wave-ranges. The last row of table 2 is the average of the all SPT data for each location and reconstruction wave-range.

The average of all the data reveal that the GFC, RMSE and RecP for the hand, face and lower back is higher for reconstruction between 430 to 780 nm. It is also important to notice that the values of other reconstructions (430 to 850 nm and 430 to 950) are all superior to 0.9967 for the GFC. According to Hernandez et al., [START_REF] Hernández-Andrés | Colorimetric and spectroradiometric characteristics of narrow-field-of-view clear skylight in granada, spain[END_REF] reconstruction with GFC higher than 0.99 is considered good.

The GFC result for each specific SPT group shows that reconstruction for SPT II and III appears to be better in the reconstruction range of 430 to 850 whereas reconstruction for SPT IV, V and VI obtains better results in the range of 430 to 780. Similar conclusion arises for the PRec values and the RMSE. The results can be visually assessed with figure 2 and figure 3 which respectively show the average data per SPT acquired from the back of the hand for the spectrophotometer and Asclepios.

The shape of the curve from Asclepios is in good accordance with the curve obtained from the spectrophotometer.

Difference between constitutive and facultative skin colours is obtained by both system. Figures 4 and5 shows respectively the average SPT III data acquired from the lower back (constitutive skin colour) and the back of the hand (facultative). The constitutive skin colour have higher spectral reflectance compare to the facultative one. The difference of spectrum shape between area exposed to the sun against area non exposed to the sun can also be visually assessed. This revelation means that quantification of skin chromophores is possible with the two systems. However, with a spectrophotometer, several acquisitions are required whereas with Asclepios and its hyperspectral cube, only one acquisition is necessary.

Both numerical and visual analysis of the results reveal that the hyperspectral cube is accurately reconstruct.

CONCLUSION

This study validates the capacity of ASCLEPIOS system to accurately reconstruct hyperspectral cube of cutaneous data, which provides more information than conventional spectrophotometer, which provides only the average spectrum of an areas. Such hyperspectral cube provides both spectral and spatial information. The capacity of the system to discriminate area with different melanin composition, shows that within a same hyperspectral cube, the spectral information can be obtained offering a wide range of analysis to be performed (spectral analysis, morphological analysis) like source separation to extract within a single acquisition or follow up of wounds healing using area segmentation on a specific monoband plan.

Figure 1 .

 1 Figure 1. Asclepios system

Figure 2 .

 2 Figure 2. Spectrophotometer data from back of the hand.

Figure 3 .

 3 Figure 3. Asclepios data from back of the hand.

Figure 4 .

 4 Figure 4. Average SPT III at two different locations acquired by CS2600d.

Figure 5 .

 5 Figure 5. Average SPT III at two different locations acquired by Asclepios.

Table 2 .

 2 GFC, PRec and RMSE for each location and reconstruction range of the average data of each SPT and average of all SPT.

	SPT II	Hand 780 Hand 850 Hand 950 Face 780 Face 850 Face 950 Back 780 Back 850 Back 950
	GFC	0,9977	0,9977	0,9970	0,9963	0,9964	0,9958	0,9963	0,9965	0,9957
	PRec	0,9842	0,9846	0,9838	0,9776	0,9777	0,9694	0,9793	0,9809	0,9724
	RMSE	18,87	18,30	18,91	32,23	32,11	44,42	35,15	32,67	47,93
	SPT III	Hand 780 Hand 850 Hand 950 Face 780 Face 850 Face 950 Back 780 Back 850 Back 950
	GFC	0,9980	0,9980	0,9974	0,9972	0,9973	0,9966	0,9974	0,9975	0,9967
	PRec	0,9852	0,9856	0,9821	0,9797	0,9674	0,9677	0,9847	0,9848	0,9745
	RMSE	12,78	12,37	20,02	22,08	37,31	38,22	20,88	20,74	34,54
	SPT IV	Hand 780 Hand 850 Hand 950 Face 780 Face 850 Face 950 Back 780 Back 850 Back 950
	GFC	0,9976	0,9975	0,9970	0,9972	0,9972	0,9966	0,9975	0,9976	0,9969
	PRec	0,9792	0,9788	0,9749	0,9802	0,9797	0,9727	0,9781	0,9775	0,9612
	RMSE	16,30	16,67	17,55	20,84	21,32	27,74	21,06	21,53	37,15
	SPT V	Hand 780 Hand 850 Hand 950 Face 780 Face 850 Face 950 Back 780 Back 850 Back 950
	GFC	0,9971	0,9968	0,9963	0,9971	0,9970	0,9965	0,9980	0,9979	0,9972
	PRec	0,9814	0,9793	0,9807	0,9836	0,9839	0,9790	0,9765	0,9767	0,9548
	RMSE	15,61	16,94	15,67	13,22	13,05	15,71	15,53	15,31	29,45
	SPT VI	Hand 780 Hand 850 Hand 950 Face 780 Face 850 Face 950 Back 780 Back 850 Back 950
	GFC	0,9972	0,9961	0,9962	0,9977	0,9969	0,9967	0,9975	0,9962	0,9961
	PRec	0,9842	0,9797	0,9790	0,9860	0,9726	0,9775	0,9631	0,9691	0,9251
	RMSE	4,89	5,92	6,35	4,44	8,51	6,80	7,81	7,03	16,20
	ALL SPT Hand 780 Hand 850 Hand 950 Face 780 Face 850 Face 950 Back 780 Back 850 Back 950
	GFC	0,9975	0,9973	0,9968	0,9971	0,9970	0,9965	0,9975	0,9974	0,9968
	PRec	0,9817	0,9807	0,9790	0,9813	0,9782	0,9738	0,9773	0,9778	0,9590
	RMSE	14,86	15,33	16,50	18,77	21,40	25,97	19,77	19,42	33,71

ACKNOWLEDGMENTS

The authors would like to thank the financial support provided by Conseil Régional de Bourgogne, France, laboratoire Bioderma, France and Cancéropôle du Grand Est, France.


Marwa Hamza

Abdelhamid Hassairi
email: abdelhamid.hassairi@fss.rnu.tn

New characterizations of the Letac-Mora class of real cubic natural exponential families

Keywords: Natural exponential family, variance function, cumulant function, Monge-Ampère equation, prior distribution

 and satisfied by the Wishart and the simple quadratic NEF's. We also show that in the real case, each of these properties provides a new characterization of the Letac-Mora class of real cubic NEF's.

Introduction and preliminaries

For the convenance of the reader, we first introduce some notations and recall some facts concerning the natural exponential families and their variance functions, our notations are the ones used in [START_REF] Letac | Lectures on natural exponential families and their variance function[END_REF]. Let E be a linear vector space with finite dimension n, denote by E * its dual, and let E * × E → IR : (θ, x) → θ, x be the duality bracket. If µ is a positive radon measure on E, then

L µ (θ) = E exp( θ, x )µ(dx) (1.1)
denotes its Laplace transform. We also denote by M(E) the set of measures µ such that the set

Θ(µ) = interior{θ ∈ E * ; L µ (θ) < +∞} (1.2)
is non empty and µ is not concentrated on an affine hyperplane of E. The cumulant function of an element µ of M(E) is the function defined for θ in Θ(µ) by k µ (θ) = log L µ (θ).

To each µ in M(E) and θ in Θ(µ), we associate the probability distribution on E P (θ, µ)(dx) = exp ( θ, xk µ (θ)) µ(dx). The set

F = F (µ) = {P (θ, µ); θ ∈ Θ(µ)}
is called the natural exponential family (NEF) generated by µ. We also say that µ is a basis of F .

The function k µ is strictly convex and real analytic. Its first derivative k ′ µ defines a diffeomorphism between Θ(µ) and its image M F . Since k ′ µ (θ) = the probability of F with mean m, we have

F = {P (m, F ); m ∈ M F } ,
which is the parametrization of F by the mean. Now the covariance operator of P (m, F ) is denoted by V F (m) and the map

M F -→ L s (E * , E); m -→ V F (m) = k ′′ µ (ψ µ (m))
is called the variance function of the NEF F . It is easy proved that for all m ∈ M F ,

V F (m) = (ψ ′ µ (m)) -1 ,
and an important feature of V F is that it characterizes F in the following sense: If F and F ′ are two NEFs such that V F (m) and

V F ′ (m) coincide on a nonempty open set of M F ∩ M F ′ , then F = F ′ .
Now, let us examine the influence of an affine transformation and a power convolution on a NEF

F = F (µ). If ϕ(x) = a(x) + b, where a ∈ GL(E) and b ∈ E, is an affine transformation of E, then ϕ(F (µ)) = F (ϕ(µ)), M ϕ(F ) = ϕ(M F ), and
V ϕ(F ) (m) = a V F (ϕ -1 (m)) a * ,
where a * is the transpose of a. On the other hand the set

Λ(µ) = {λ > 0; ∃ µ λ ∈ M(E) such that L µ λ (θ) = (L µ (θ)) λ for all θ ∈ Θ(µ)}
is called the Jorgensen set of µ and the measure µ λ is the λ-power of convolution of µ.

For λ in Λ(µ), we have that

M F (µ λ ) = λM F , and V F λ (m) = λ V F ( m λ ).
A very interesting fact is that the most common real and multivariate probability distributions belong to the natural exponential families such that the variance function is a polynomial of degree less then or equal to three in the mean m. For instance, up to affine transformations and power of convolution (up to the type), the Gaussian, Poisson, gamma, binomial, negative binomial and hyperbolic cosine distributions form the class of all real NEF's whose variance function is a polynomial of degree less than or equal to 2 characterized by Morris [START_REF] Morris | Natrual exponential families with quadratic variance-function[END_REF]. Letac and Mora [START_REF] Letac | Natural real exponential familiies with cubic variance functions[END_REF] have added six others types of distributions, namely, the inverse Gaussian, Ressel, Abel, Tackàs, strict arcsine and large arcsine, to get the class of real cubic NEF's, that is the class of NEF's such that variance function is a polynomial of degree less than or equal to three. The classification of NEF's with polynomial variance function have been extended to the multivariate NEF's. The multivariate version of the Morris class, called the class of simple quadratic NEF's, has been completely described by Casalis [START_REF] Casalis | The 2d+4 simple natrual exponential families on IR d[END_REF], it contains 2n + 4 types. Hassairi [6] has defined and characterized the so-called class of multivariate simple cubic NEF's which is the natural extension of the class of real cubic NEF's. It is worth mentioning here that the simple quadratic NEF's are not the only families which have quadratic variance functions, the Wishart families on symmetric matrices have also quadratic variance functions. The classifications of NEF's by the form of the variance function provide an important tool in the study of distributions. In fact, in many important cases, the variance function is very simple and is easier to use than the distribution itself or the Laplace transform. Moreover, the fact that the variance function is quadratic or cubic, is not only a question of form, but the form corresponds to some very interesting analytical characteristic properties. In this respect, let us mention that for the Morris class of real quadratic NEF's, we have the Meixner characterization based on some families of orthogonal polynomials which generate exactly the Morris class (see [START_REF] Meixner | Orthogonal Polynomsysteme mit einer besonderen Gestalt der erzengenden function[END_REF]). Another characterization due to Feinsilver [START_REF] Feinsilver | Some classes of orthogonal polynomials associated with martingales[END_REF] states that a certain class of polynomials naturally associated to a NEF is orthogonal if and only if the family is in the Morris class. This characterization has been extended to the Casalis class of simple quadratic NEF's by Labeye-Voisin, and Pommeret [START_REF] Labeye-Voisin | Polynmes orthogonaux associs aux familles exponentielles de[END_REF]. Concerning the cubic NEF's, Hassairi and Zarai [START_REF] Hassairi | Characterization of the cubic exponential families by orthogonality of polynomials[END_REF] introduced a notion of 2-orthogonality for a sequence of polynomials to give an extended version of the Meixner and Feinsilver characterization which subsume the Letac-Mora class of real cubic NEF's. Hassairi and Zarai [START_REF] Hassairi | Characterization of the simple cubic multivariate exponential families[END_REF] have also introduced a notion of trans-orthogonality for a sequence of multivariate polynomials to extend their characterization result to the class of multivariate simple cubic NEF's. Besides these characterizations based essentially on different notions of orthogonality of polynomials, it is stated in Casalis [START_REF] Casalis | The 2d+4 simple natrual exponential families on IR d[END_REF] that the simple quadratic NEF satisfies a property based on two conjugates families of prior distributions related to the NEF. For a NEF F = F (µ), consider the family of prior distributions Π introduced by Diaconis and Ylvisaker [START_REF] Diaconis | Conjugate priors for exponential families[END_REF] and defined by

Π = {π t,m 0 (dθ) = C t,m 0 exp t( m 0 , θ -k µ (θ)) 1 Θ(µ) (θ)dθ, t > 0, m 0 ∈ M F } (1.3)
where C t,m 0 is a normalizing constant. Consider also the family Π * introduced by Consonni et al [START_REF] Consonni | Reference priors for exponential families with simple quadratic variance function[END_REF], see also [START_REF] Gutierrez-Pena | Reference priors for exponential families[END_REF] and defined by

Π * = {π * t,m 0 , t > 0, m 0 ∈ M F } (1.4) where π * t,m 0 (dm) = C * t,m 0 exp t( m 0 , ψ µ (m) -k µ (ψ µ (m)) ) 1 M F (m)dm,
and the constant C * t,m 0 is a normalizing constant. Then, when F (µ) is a Wishart or a simple quadratic NEF, we have that k ′ µ (Π) = Π * . It is also shown that this property is equivalent to two other properties expressed in terms of some differential equations satisfied by the cumulant function k µ . In the real case, the property characterizes the Morris class of real quadratic NEF's, that is k ′ µ (Π) = Π * if and only if the NEF F is in the Morris class. In the present paper, we extend these results to the class of multivariate simple cubic NEF's. We construct two families of prior distributions related to a multivariate NEF, and we show that these families coincide when the NEF is simple cubic. We then show that this property is equivalent to the fact that the cumulant function is a solution of some Monge-Ampère equation and also equivalent to the fact that the variance function satisfies a differential equation. As a corollary, we obtain three new characterizations of the Letac-Mora class of real cubic NEF's.

Some equivalent properties

Throughout this section, we suppose that F = F (µ) is a NEF on a linear vector space E with dimension n. Besides the family Π of prior distributions defined in (1.3), we introduce another family Π of prior distributions. Let β be in E * such that the set

Θ = {θ ∈ Θ(µ); 1 + β, k ′ µ (θ) > 0}
is nonempty, and denote M = k ′ µ ( Θ). Consider the family of prior distributions

Π = { π t,m 0 ; t ∈ IR * + , m 0 ∈ M F }, (2.5)
where

π t,m 0 (dm) = C t,m 0 (1 + β, m ) -n-2 exp t{ m 0 , ψ µ (m) -k µ (ψ µ (m))} 1 M (m)dm.
With these notations, we next state and prove our first main result.

Theorem 2.1 The three following properties are equivalent (1) There exists (a, b, c)

∈ E × IR 2 such that for all m in M F , det V F (m) = (1 + β, m ) n+2 exp{ a, ψ µ (m) + bk µ (ψ µ (m)) + c}.
(2) There exists (a, b) ∈ E × IR such that for all m in M F and any basis, (e i ) n i=1 of E, with dual basis (e * i ) n i=1 , we have

n i=1 [V ′ F (m)(e i )]e * i = n + 2 1 + β, m V F (m)(β) + a + bm. ( 2

.6)

(3) There exists an open subset

Ω of IR * + × M F such that k ′ µ (Π) = Π Ω = { π t,m 0 ; (t, m 0 ) ∈ Ω}.
Note that (1) may be stated in terms of the cumulant function as there exists (a, b, c) ∈ E × IR 2 such that for all θ in Θ(µ),

det k ′′ µ (θ) = (1 + β, k ′ µ (θ) ) n+2 exp{ a, θ + bk µ (θ) + c},
that is the cumulant function is solution of some Monge-Ampère equation (see [START_REF] Zuily | Local existence and regularity of the Dirichlet problem for the Monge-Ampère equation, Journes "équations aux Drives Partielles[END_REF]). Proof We will show that (1) ⇔ (2) and ( 1) ⇔ (3).

(1) ⇒ (2) Suppose that V F (m) satisfies ( 1), then we have

log det V F (m) = (n + 2) log(1 + β, m ) + { ψ µ (m), a + b k µ (ψ µ (m)) + c}.
Taking the derivative, we get

trace(V -1 F (m)V ′ F (m)(.)) = (n + 2) β, . 1 + β, m + ψ ′ µ (m)(.), a + b m, ψ ′ µ (m)(.) , which is equivalent to n i=1 [V ′ F (m)(.)V -1 F (m)(e i )](e * i ) = (n + 2) β, . 1 + β, m + ψ ′ µ (m)(.), a + b m, ψ ′ µ (m)(.) .
Replacing (.) by V F (m)(.), and using the condition of symmetry [START_REF] Letac | Lectures on natural exponential families and their variance function[END_REF], page 103), we obtain

V ′ (m)(V (m)(α))(β) = V ′ (m)(V (m)(β))(α) ∀ α, β ∈ E * (see
n i=1 [V ′ F (m)(e i )(.)](e * i ) = (n + 2) β, V F (m)(.) 1 + β, m + a, (.) + b m, (.) . (2.7) As V F (m) is symmetric, we get n i=1 [V ′ F (m)(e i )](e * i ) = (n + 2) 1 + β, m V F (m)(β) + a + bm.
(2) ⇒ (1) Suppose that (2) holds. Then, we easily get (2.7). Replacing, in (2.7), (.) by V -1 F (m)(.), one obtains

n i=1 V ′ F (m)(e i )V -1 F (m)(.)(e * i ) = (n + 2) β, (.) 1 + β, m + a, ψ ′ µ (m)(.) + b m, ψ ′ µ (m)(.) .
This is equivalent to

trace(V -1 F (m)V ′ (m)(.)) = (n + 2) β, (.) 1 + β, m + a, ψ ′ µ (m)(.) + b m, ψ ′ µ (m)(.) .
Integrating, we deduce that there exists c in IR such that

log det(V F (m)) = (n + 2) log(1 + β, m ) + { ψ µ (m), a + b k µ (ψ µ (m)) + c},
and the result follows.

(1) ⇒ (3) Suppose that (1) holds, and define

Ω = {(t, m 0 ) ∈ IR * + × M F ; t > b and m 0 ∈ (1 - b t )M F - a t }.
Take (t, m 0 ) in Ω and denote ν the image of π t,m 0 by ψ µ . Then it is easy to verify that

ν(dθ) = C t,m 0 e c exp{ tm 0 + a, θ -(t -b)k µ (θ)} 1 Θ (θ)dθ. Since (t, m 0 ) is in Ω, we have that t -b > 0 and tm 0 + a t -b ∈ M F .
Thus taking t 1 = tb and m 1 = tm 0 + a tb , we obtain that

ν(dθ) = C t 1 ,m 1 exp{t( m 1 , θ -k µ (θ))} 1 Θ (θ)dθ.
Hence ψ µ ( Π Ω ) ⊂ Π, and it follows that Π Ω ⊂ k ′ µ (Π). Conversely, if π t,m 0 is an element of Π, then its image σ by k ′ µ is given by

σ(dm) = C t,m 0 e -c (1 + β, m ) -n-2 exp{ tm 0 -a, ψ µ (m) -(t + b)k µ (ψ µ (m))} 1 M (m)dm.
Taking

t 1 = t + b and m 1 = tm 0 -a t + b . Then (t 1 , m 1
) is in Ω, and we have

σ(dm) = C t 1 ,m 1 (1 + β, m ) -n-2 exp t 1 { m 1 , ψ µ (m) -k µ (ψ µ (m))}1 M (m)dm ,
which is an element of Π Ω .

(3) ⇒ (1) Suppose that k ′ µ (Π) = Π Ω . Then, for an element π t,m 0 of Π, we have on the one hand,

k ′ µ (π t,m 0 )(dm) = (det V F (m)) -1 C t,m 0 exp t{ m 0 , ψ µ (m) -k µ (ψ µ (m))} 1 M (m)dm.
On the other hand, since

k ′ µ (π t,m 0 ) is in Π Ω , there exists (t 1 , m 1 ) in Ω such that k ′ µ (π t,m 0 )(dm) = C t 1 ,m 1 (1 + β, m ) -n-2 exp t 1 { m 1 , ψ µ (m) -k µ (ψ µ (m))} 1 M (m)dm. Comparing these two expressions of k ′ µ (π t,m 0 ) gives det V (m) = (1 + β, m ) n+2 exp{ a, ψ µ (m) + bk µ (ψ µ (m)) + c}, where a = tm 0 -t 1 m 1 , b = t 1 -t, and c = log( c t,m 0 c t 1 ,m 1
).

3 Characterizations of the Letac-Mora class of real cubic NEFs

In this section, we prove that a multivariate simple cubic NEF satisfies the properties in Theorem (2.1), and that the real versions of these properties characterize the real cubic NEFs. Recall that a simple cubic NEF is obtained form a simple quadratic NEF by the so-called action of the linear group GL(IR × E) on the NEFs of E. For more details, we refer the reader to [START_REF] Hassairi | La classification des familles exponentielles naturelles sur IR n par l'action du groupe linaire de IR n+1[END_REF], where a complete description of this class is given. This action is in fact an extension of the way in which the Letac-Mora class of real cubic NEFs is obtained from the Morris class of real quadratic NEF's. For our purposes here, we need only to mention that, up to affine transformations and power of convolution, a simple cubic variance function is of the form

V (m) = (1 + β, m ) (I + m ⊗ β) V 1 ( m 1 + β, m ) (I + β ⊗ m), (3.8)
where V 1 is the variance function of a simple quadratic NEF F 1 , and m is in (M F 1 ) β , where

(M F 1 ) β = {m ∈ M F 1 ; 1 + β, m > 0 and m 1 + β, m ∈ M F 1 }.
The relation (3.8) is invertible and conversely, we have

V 1 (M ) = (1 -β, M ) (I -M ⊗ β) V ( M 1 -β, M ) (I -β ⊗ M ), (3.9)
where M is in (M F ) -β . We also mention that the relation between a simple cubic NEF F (µ) and a simple quadratic NEF F (ν) may also be expressed in terms of the cumulant functions by

     k µ (λ) = k ν (θ) -k 0 λ = -βk ν (θ) + θ -λ 0 (3.10)
or equivalently by

     k ν (θ) = k µ (λ) -k 1 θ = βk µ (λ) + λ -θ 1 (3.11)
where (k 0 , λ 0 ) and (k 1 , θ 1 ) are constants in IR × E. Note that if β = 0 in (3.8) we obtain the simple quadratic class. Then for more accuracy we exclude this case and we keep only β in E * \ {0}.

We now prove that the multivariate simple cubic NEF's satisfy the properties in Theorem 2.1.

Proposition 3.1 Let F = F (µ) be a simple cubic NEF on E, then there exists (a, b, c) in E * × IR 2 such that det(V F (m)) = (1 + β, m ) n+2 exp{ ψ µ (m), a + b k µ (ψ µ (m)) + c}.
Proof Given that the family F is simple cubic, then there exist β in E * and F 1 = F (ν) a simple quadratic NEF such that

V F (m) = (1 + β, m ) (I + m ⊗ β) V F 1 ( m 1 + β, m ) (I + β ⊗ m), see (3.8). As det(I + m ⊗ β) = 1 + β, m , we obtain det(V F (m)) = (1 + β, m ) n+2 det(V F 1 ( m 1 + β, m
)).

We now use the fact for a simple quadratic NEF F 1 (see [START_REF] Casalis | The 2d+4 simple natrual exponential families on IR d[END_REF]), there exist a ′ in E * and b ′ , c ′ in IR such that, for all As the Letac-Mora class of real cubic NEFs is nothing but the simple cubic class, when the dimension n is equal to 1, this class satisfies the real version of the properties in Theorem 2.1. We will show that, in this case, these properties are characteristic. Solving this differential equation by standard methods gives

M in M F 1 , det(V F 1 (M )) = exp{ a ′ , ψ ν (M ) + b ′ k ν (ψ ν (M )) + c ′ }.
V F (m) = λ (1 + βm) 3 - b β 2 (1 + βm) 2 + b -βa 2β 2 (1 + βm),
which is a polynomial of degree less then or equal to 3. 2

E

  xP (θ, µ)(dx), M F is called the domain of the means of F . The inverse function of k ′ µ is denoted by ψ µ and setting P (m, F ) = P (ψ µ (m), µ)

  It follows that det(V F (m)) = (1 + β, m ) n+2 exp{ ψ ν ( m 1 + β, m ), a ′ + b ′ k ν (ψ ν ( m 1 + β, m)) + c ′ }.From (3.10), putting λ = ψ µ (m) and θ = ψ ν ( m 1 + β, m), we getk ν (ψ ν ( m 1 + β, m )) = k µ (ψ µ (m)) + k 0 , ψ ν ( m 1 + β, m ) = ψ µ (m) + βk µ (ψ µ (m)) + βk 0 + λ 0 . Then det V F (m) = (1 + β, m ) n+2 exp{ ψ µ (m), a ′ + (b ′ + a ′ , β )k µ (ψ µ (m)) +(b ′ k 0 + a ′ , βk 0 + λ 0 + c ′ )}.Setting a = a ′ , b = b ′ + a ′ , β and c = b ′ k 0 + a ′ , βk 0 + λ 0 + c ′ , we obtain the desired result. 2

Theorem 3 . 2

 32 Let F = F (µ) be a NEF on the real line, then F is cubic if and only ifk ′ µ (Π) = Π.Proof Suppose that k ′ µ (Π) = Π. Then according to Theorem(2.1), the variance function V F (m) satisfies the differential equation(1 + βm)V ′ F (m) -3β V F (m) = (a + bm)(1 + βm).


Sandro Bertolino
email: sandro.bertolino@unito.it

Barbara Ingegno

Marco Girardello

Modelling the habitat requirements of invasive Eastern Cottontail (Sylvilagus floridanus) introduced to Italy

Keywords: Introduced species, Lagomorph, Management, Species-habitat relationship

We analysed the spatial distribution of Eastern Cottontail (Sylvilagus floridanus) and its relationship with the habitat attributes of the agro-ecosystems in north-western Italy. We used generalised linear models to analyse the relationship between the occurrence of the cottontail and a series of landscape-level variables. Hierarchical partitioning was used to estimate the independent contribution of these variables. Our analysis showed that cottontail is typically a lowland species, in that it selected municipalities below 400 m a.s.l. and avoided higher elevations. Cottontails were more likely to occur where the edges of meadows and crop fields were higher. This might indicate the importance of these herbaceous habitats in cottontail feeding behaviour and also the significance of edge vegetation in providing cover. In addition, the extent of the hydrographic network positively affected cottontail presence, probably due to the permanent vegetation cover associated with the banks of rivers, streams and channels of the otherwise intensively cultivated and densely inhabited Po Plain. On the other hand, a prevalence of woodlands, tree plantations and large pastures negatively affected cottontail presence. This is surprising, because these habitats usually encourage cottontail presence in their native range. In general, however, cottontails prefer small woodlots and utilise only the edges of forests, while in Italy forests and pastures are large and typically located in hilly and mountain areas generally avoided by this species. In conclusion, our results confirm that the cottontail is an edge species whose presence is encouraged by crop margins, natural herbaceous habitats and extensive hydrographic networks with the associated riverside vegetation. These characteristics are common in the Po Plain and are likely to facilitate a further expansion of cottontail.

Introduction

The introduction of species is changing ecosystem composition and function [START_REF] Ricciardi | Are modern biological invasions an unprecedented form of global change?[END_REF]. For this reason, the management of non-native species has become an important factor in wildlife and habitat conservation, as well as a priority for preserving human activities [START_REF] Wittenberg | Invasive alien species: a toolkit of best prevention and management practices[END_REF]. The promotion of appropriate management strategies requires a deep understanding of the factors that either affect species distribution or abundance. This kind of information is generally derived either from intensive field work [START_REF] Gurnell | The effects of food availability and winter weather on the dynamics of a grey squirrel population in southern England[END_REF][START_REF] Bonesi | Evaluation of sign surveys as a way to estimate the relative abundance of mink[END_REF] or from habitat models [START_REF] Barreto | The role of habitat and mink predation in determining the status and distribution of water voles in England[END_REF][START_REF] Guisan | Predictive habitat distribution models in ecology[END_REF][START_REF] Tattoni | Modelling the expansion of grey squirrels (Sciurus carolinensis) in Lombardy, Northern Italy: implications for squirrel control[END_REF][START_REF] Bertolino | Predicting the spread of the American grey squirrel (Sciurus carolinensis) in Europe: a call for a co-ordinated European approach[END_REF]. Model-based analyses of species-habitat relationships can help to clarify which factors influence the establishment and spread of non-native species [START_REF] Bertolino | Modelling the distribution of an introduced species: the coypu Myocastor coypus (Mammalia, Rodentia) in Piedmont region, NW Italy[END_REF][START_REF] Roura-Pascual | From introduction to equilibrium: reconstructing the invasive pathways of the Argentine ant in a Mediterranean region[END_REF]. In turn, knowledge of the ecological requirements of these species is a prerequisite for the development of appropriate mitigation strategies [START_REF] Braysher | Managing vertebrate pests: principles and strategies[END_REF][START_REF] Bertolino | Effectiviness of coypu control in small Italian wetland areas[END_REF].

The Eastern Cottontail (Sylvilagus floridanus) is a lagomorph native to America and more precisely to southern Canada, USA, central and north-western South America [START_REF] Chapman | Sylvilagus floridanus[END_REF]. The species was introduced to other areas of the USA [START_REF] Chapman | Sylvilagus floridanus[END_REF]) and to several European countries (DAISIE 2009) for hunting purposes. As to date, wild populations in Europe that have survived are only located in Italy where the earliest introduction event dates back to 1966 [START_REF] Mussa | Il Silvilago in Provincia di Torino[END_REF]DAISIE 2009). The species expanded its range eastward from the original area of introduction, across the Po Plain, mainly through the network of rivers which still retain some vegetation cover [START_REF] Silvano | Distribution of the Eastern Cottontail Sylvilagus floridanus) in the province of Alessandria[END_REF][START_REF] Vidus-Rosin | Introduced Lagomorphs as a threat to "native" Lagomorphs: The case of the Eastern Cottontail (Sylvilagus floridanus) in northern Italy[END_REF]. Nowadays, the species is present in a large area of north-western Italy that includes most of Piedmont and the western portion of Lombardy, while smaller populations are located in north-eastern and central regions of Italy [START_REF] Angelici | Mammalia II: Erinaceomorpha, Soricomorpha, Lagomorpha, Rodentia[END_REF].

The cottontail is a possible vector of diseases that can be transmitted to native lagomorphs [START_REF] Jacobson | Disease and physiologic characteristics of two cottontail populations in Virginia[END_REF][START_REF] Tizzani | Presence of infectious agents and parasites in wild populations of cottontail (Sylvilagus floridanus) and consideration on its role in the diffusion of pathogens infecting hares[END_REF]. In particular, it can be a carrier of pseudotuberculosis and myxomatosis. The former is transmissible to hares, which are prone to it, and the latter is lethal to the European rabbit (Oryctolagus cuniculus). Imported cottontails introduced to Italy have also introduced seven species of Eimeria, a genus of protozoan intestinal parasites of North American origins [START_REF] Bertolino | Richness, origin and structure of an Eimeria community in a population of Eastern Cottontail (Sylvilagus floridanus) introduced into Italy[END_REF]. Studies conducted in France and Italy have highlighted the risks connected to crop damage as well as those related to possible competition for food and space with other lagomorphs [START_REF] Arthur | L'introduction de Sylvilagus floridanus en France: historique, dangers et expérimentation en cours[END_REF][START_REF] Chapuis | Etude de Sylvilagus floridanus en region mediterraneenne dans des conditions de semi-liberte. Regime alimentaire au cours d'un cycle annuel par l'analyse micrographique des feces[END_REF][START_REF] Vidus-Rosin | Introduced Lagomorphs as a threat to "native" Lagomorphs: The case of the Eastern Cottontail (Sylvilagus floridanus) in northern Italy[END_REF], 2009).

Cottontail habitat requirements in Italy have been investigated by means of field studies that included the observation of the animals [START_REF] Vidus-Rosin | Introduced Lagomorphs as a threat to "native" Lagomorphs: The case of the Eastern Cottontail (Sylvilagus floridanus) in northern Italy[END_REF], 2009), radiotracking [START_REF] Bertolino | Progetto di ricerca sull'ecologia della Lepre comune (Lepus europaeus) e del Silvilago (Sylvilagus floridanus) in un'area di sovrapposizione spaziale nella Riserva Naturale del Torrente Orba[END_REF]) and pellet-plot counts [START_REF] Vidus-Rosin | Introduced Lagomorphs as a threat to "native" Lagomorphs: The case of the Eastern Cottontail (Sylvilagus floridanus) in northern Italy[END_REF]. The main results of these local scale studies have shown a strong association of the species with herbaceous habitats and fallow lands. These two broad habitat categories comprise hedgerows and other shrubby habitats that provide permanent cover and small patches of cultivations that provide food supply. The aim of the present study was to analyse the spatial distribution of cottontail in relation to habitat attributes of agro-ecosystems at a large (regional) scale in an area where the cottontail was first introduced and that include most of the present range of the species.

Methods

Study area

The study area is Piedmont, a region in north-western Italy (Fig. 1). Piedmont can be divided into four main geographical areas: the Alps (with peaks above 1,500 m a.s.l.), the pre-Alps (with peaks below 1,500 m a.s.l.), inland low hills and Apennines (the hills bordering the River Po, Monferrato, Langhe and Apennines uplands) and lowlands (including wide Alpine valley floors below 600 m a.s.l.). The total surface area is 25,383 km 2 , with 29% of lowlands and 71% of mountains and hills; 54% of the total is below 1,000 m of elevation [START_REF] Gottero | I boschi del Piemonte, conoscenze ed indirizzi gestionali[END_REF]. At the time of the study, most of the region was covered by agricultural lands (37%) and rangelands (13%). Woodlands were mainly spread across the hilly areas, Alps and Apennines and almost absent from lowlands. They covered 34% of the surface. Urban areas comprised 1,209 municipalities (mean area 2,100±2,299 ha), facilities and road infrastructures and covered 6% of the total surface [START_REF] Gottero | I boschi del Piemonte, conoscenze ed indirizzi gestionali[END_REF]).

Habitat analysis

Measurements of landscape attributes were obtained from a CORINE (third level, scale 1:100,000) digital database with For each municipality considered in the analysis, we used 17 variables to describe the structure and composition of the landscape, thus focusing on those features that are relevant to the ecological requirements of cottontail (Table 1). We used ArcView 3.1 GIS software (ESRI-Environmental Systems Research Institute, California) to assess the total surface area, the mean elevation and the length of the hydrographic network (elevation and hydrographic network source: Regione Piemonte website at http://www.regione.piemonte.it/repertorio). The landscape composition was assessed in terms of the relative proportions of the municipality's area covered by woodlands and tree plantations, open habitats, wetlands, crops and urban areas. The pattern and heterogeneity of the landscape mosaic were assessed using two indexes of patch characteristics, describing mean patch size and edge density (amount of edges in proportion to the relative area). These were calculated for each municipality and habitat category using the Patch Analyst extension of Arc View 3.1 (Table 1).

Data analysis

The distribution of the cottontails in 2000-2002 was mapped by compiling information provided by the Osservatorio Faunistico regionale (regional Wildlife Service) and the Wildlife Services of the provinces. The Osservatorio Faunistico regionale receives the data regarding the damage caused by wildlife within the region and is responsible for the collection of information on species presence. Further information was collected and made available by the staff of parks and protected areas, local hunting units, wildlife managers and scientists. Records of cottontails were collected at a municipality level, and we used only information derived from direct animal observations. The cottontail range in 2002 is reported in Fig. 1. The map was produced joining the 467 out of 1,209 municipalities where the species was present. The species distribution was modelled on a municipality scale.

The relationship between the presence of the cottontail and the independent variables was analysed by the means of a generalised linear model (GLM). Because there were no records of sites from which cottontails were surely absent, 467 further municipalities were selected at random throughout Piedmont following [START_REF] Stockwell | Controlling Bias in Biodiversity Data[END_REF]. However, it is possible that some of the presumed absences are sites with cottontail, i.e. false negatives. In order to assess the impact of false negatives on the analysis, a false record permutation test (FRPT) was carried out [START_REF] Greaves | Predicting species occurrence using information-theoretic approaches and significance testing: an example of dormouse distribution in Cumbria, UK[END_REF]. The model was run with varying numbers of false absences and false presences; all sites and parameters were included. Initially, 437 of the presence sites were converted at random to false absences. This was carried out 50 times. This process was repeated, sequentially decreasing the number of false absences by ten, then increasing the number of false presences by ten until there were 437 false presences. Deviance explained (D 2 ) values were used to assess the impact of this on model output. The Akaike information criterion (AIC) was used to select the minimum adequate model (MAM; [START_REF] Manel | Evaluating presenceabsence models in ecology: the need to account for prevalence[END_REF]. Model selection was carried out using the stepAIC function from R package MASS [START_REF] Venables | Modern applied statistics with S, 4th edn[END_REF] function is a stepwise algorithm that sequentially searches through all possible models for the one that minimises the AIC. The stepwise search was performed both ways (i.e. backward and forward), as regular one-step approximations are considered less robust and less profitable [START_REF] Venables | Modern applied statistics with S, 4th edn[END_REF]. We finally estimated the significance of the parameters of our MAM using z tests and evaluated the predicted performance of the final model using receiver operating characteristic (ROC) plot [START_REF] Fielding | A review of methods for the assessment of prediction errors in conservation presence: absence models[END_REF]. Because no independent data were available for model validation, we used a tenfold cross-validation approach.

A cross-validation is made with subsets of the training dataset, where each subset contains an equal number of randomly selected data points. Each subset is then dropped from the model, the model is recalculated and predictions were made for the omitted data points. A ROC plot was calculated for the output of the cross-validation. The discrimination ability of each of the ten random portions of data was assessed using the area under the curve (AUC; [START_REF] Fielding | A review of methods for the assessment of prediction errors in conservation presence: absence models[END_REF]. The AUC values range from 0.5 for models with no discrimination ability to 1 for models with perfect discrimination [START_REF] Swets | Measuring the accuracy of diagnostic systems[END_REF]).

We used hierarchical partitioning (HP) analyses to calculate the independent contribution of each predictor selected by the MAM. In HP all possible models for the distribution of a species are considered in a hierarchical regression setting. HP involves measuring the increase in the goodness-of-fit of all models with a particular variable compared with the equivalent model without that variable (Mac Nally and Horrocks 2002;[START_REF] Luoto | Determinants of biogeographical distribution of butterflies in boreal regions[END_REF][START_REF] Radford | The relative importance of landscape properties for woodland birds in agricultural environments[END_REF]. The improvement in fit is then averaged across all possible models in which that variable occurs to provide a measure of its independent effects. HP is considered superior to other multiple regression techniques for inferring probable causality in multivariate data sets [START_REF] Watson | Determinants of diversity in a naturally fragmented landscape: humid montane forest avifaunas of Mesoamerica[END_REF]. We specified a logistic model with log-likelihood as the goodness-of-fit measure.

To assess the degree of utilization of different elevation ranges by cottontail, we compared the observed frequency of the species at each elevation range (proportion of municipalities where cottontails were recorded, at each elevation range) with the expected values (proportion of municipalities at each elevation range) by means of the chisquare goodness-of-fit test and Bonferroni's confidence interval analysis [START_REF] Neu | A technique for analysis of utilization-availability data[END_REF][START_REF] Manly | Resource selection by animals[END_REF].

All of the analyses were carried out in R version 2.9.2 (R Development Core Team 2005).

Results

The 467 municipalities (38.6% of the total) where the cottontail was recorded covered a surface of 9,421 km 2 (37.1% of the total surface). The proportion of municipalities where cottontails were recorded differed according to the province (χ 2 =227, df=7, P<0.001): in the provinces of Asti, Vercelli, Alessandria, and Biella, the presence of cottontail was large, while in Verbania, the species was absent (Fig. 2). The species used municipalities below 400 m a.s.l. more than expected and avoided municipalities at higher elevations (Fig. 3).

The FRPT analysis carried out to assess the impact of using pseudo-absences in the model showed that the highest D 2 value was obtained when no false presences or false absences were included in the analysis. It was therefore justified to use the pseudo-absences in the analysis. The D 2 values for all combinations of false presences and absences included in the analysis are shown in Fig. 4.

Ten of seventeen environmental variables were retained in the minimum adequate model. The probability of cottontail occurrence was positively related to the amount of edge in crops and meadows, the hydrography extension, the mean size of meadow and pasture patches, and the mean size of all habitat patches contained within each municipality (Table 2). The variables related to the mean size of patches were selected in the model with coefficients close to zero; therefore, they probably did not really affect cottontail presence. Instead, cottontail presence was negatively related to the extension of municipalities covered by meadows, pastures, wetlands and woodlands. The crossvalidation exercise produced an AUC score of 0.92±0.02 (mean ± SD; Fig. 5).

The hierarchical partitioning analysis indicates that woodlands and hydrography brought the highest independent contribution in explaining cottontail presence (Fig. 6). A second group of variables were related to pastures cover, crops edges and pastures and all habitats sizes.

Provinces Proportion with cottontail

Surface % Municipalities % Fig. 2 Percentage of municipalities (number) and proportion of surface area with cottontail presence for each province. As for the proportion of surface, the areas of the municipalities in the province where cottontail was reported were summed and compared to the total surface area of the province. Symbols are explained in Fig. 1 Discussion The present distribution of the Eastern Cottontail in Piedmont was modelled as a function of landscape attributes. The probability to find cottontails increased in municipalities where crops and meadows patches had a high perimeter/area ratio, i.e. a higher amount of edges, and an extensive hydrographic network. The mean size of meadows and pastures patches, as well as of all habitats patches, were positively related with the cottontail occurrence, but their effect was probably low. On the other hand, increasing the surface of municipalities covered by meadows and pastures, wetlands and woodlands, including tree plantations, had a negative effect.

The "edge" variables of open habitats such as crops and meadows that were retained in the minimum adequate model showed a positive effect on the species occurrence. This means that the extension of ecotones is more important than the overall surface covered by these habitats. An increase in the total surface of meadows in a municipality had a negative effect, while the probability of cottontail occurrence increased when meadows had a more convoluted shape (i.e. higher edge extension). Besides providing food, crops and meadows are often surrounded by residual of natural vegetation that may assure cover. According to [START_REF] Swihart | Habitat feature influencing use of farmstead shelterbelts by the Eastern Cottontail Sylvilagus floridanus[END_REF], small patches of cultivations that are characterised by the growth of edges in proportion to the relative area may provide a source of food supply. In their native range, cottontail search areas with dense shrubs or other escape cover that is found near open foraging areas such as grasslands and pastures [START_REF] Trent | Home range and survival of cottontail rabbits in south-western Wisconsin[END_REF][START_REF] Althoff | Daytime habitat selection by cottontails in central Pennsylvania[END_REF][START_REF] Mankin | A regional model of the Eastern Cottontail and land-use changes in Illinois[END_REF][START_REF] Bond | Habitat use by cottontail rabbits across multiple spatial scales in Mississippi[END_REF]. A similar behaviour was observed also in Italy [START_REF] Vidus-Rosin | Introduced Lagomorphs as a threat to "native" Lagomorphs: The case of the Eastern Cottontail (Sylvilagus floridanus) in northern Italy[END_REF]; in Piedmont, during night censuses, cottontails were mainly found in edge environments, between woody vegetation and open lands [START_REF] Silvano | Distribution of the Eastern Cottontail Sylvilagus floridanus) in the province of Alessandria[END_REF]. The size of the hydrographic network positively affected the cottontail presence, while an increase in the surface covered by wetlands had a negative effect. The cottontail is not an aquatic species, but the extension of rivers and other basins resulted important for the species. This can be arguably explained by the presence of residual vegetation cover on the banks of the rivers of the heavily inhabited plain of Piedmont whose natural vegetation was drastically reduced by intensive cultivation. Other authors suggested that cottontails spread from their area of introduction by progressing through the river network of the Po Plain [START_REF] Silvano | Distribution of the Eastern Cottontail Sylvilagus floridanus) in the province of Alessandria[END_REF][START_REF] Vidus-Rosin | Introduced Lagomorphs as a threat to "native" Lagomorphs: The case of the Eastern Cottontail (Sylvilagus floridanus) in northern Italy[END_REF].

The increase of pastures and woodlands negatively affected the likelihood of cottontail occurrence. Woodland cover in Piedmont is fragmented into small forests and tree plantations at low elevations while forms large forests at higher elevations, i.e. across hills and mountains. In Illinois, [START_REF] Mankin | A regional model of the Eastern Cottontail and land-use changes in Illinois[END_REF] reported the extent of the decline of cottontail to be lesser in areas with more pastures, hay and woodlands. However, most of the woodlands were composed of fragmented woodlots and narrow riparian borders along rivers, interspersed with agricultural crops, thus suggesting that the cottontail is not so much a forest as an edge species [START_REF] Mankin | A regional model of the Eastern Cottontail and land-use changes in Illinois[END_REF]. [START_REF] Beckwith | Ecological succession on abandoned farmlands and its relationship to wildlife management[END_REF] observed that in abandoned fields, the number of cottontail decreases as succession proceeds towards a habitat with an increase in tree and canopy cover that usually accounts for a decrease in shrubby ground cover. In general, it seems that cottontail only enter either small woodlots or external portions of large forests and avoids their core [START_REF] Beckwith | Ecological succession on abandoned farmlands and its relationship to wildlife management[END_REF][START_REF] Mankin | A regional model of the Eastern Cottontail and land-use changes in Illinois[END_REF]. On the contrary, pastures were found to be an important feeding habitat in North America [START_REF] Trent | Home range and survival of cottontail rabbits in south-western Wisconsin[END_REF][START_REF] Mankin | A regional model of the Eastern Cottontail and land-use changes in Illinois[END_REF]. The negative effect of pastures extension on the likelihood of cottontail occurrence in Piedmont may be explained by the fact that this habitat is mainly present in mountain areas, which are generally avoided by the cottontail. The cottontail proved to be a lowland species, in that it selected municipalities below 400 m a.s.l. and avoided higher elevations. On the other hand, pastures covered only a mean (±SE) of 1.4 (0.1)% of the surface area of the municipalities below 600 m a.s.l., and a remarkable 14.3(11.3)% at higher elevations (unequal variance t test: t=-24.44, p<0.001).

Overall, cottontail resulted to be a lowland species, typical of open herbaceous habitats near vegetation that provides cover and the edges between cultivations and natural vegetation. Its presence is favoured by crops and meadows with high ecotones extension and wide hydrographic networks associated with riverside vegetation. Our results are in agreement with those reported by [START_REF] Vidus-Rosin | Introduced Lagomorphs as a threat to "native" Lagomorphs: The case of the Eastern Cottontail (Sylvilagus floridanus) in northern Italy[END_REF] in their study at a local scale. These authors find cottontails in woody and row habitats, including field margins and hedgerows with herbs and woody layers along fields and streams. Thus, in Italy the cottontail has a pattern of habitat use that is similar to that reported in the agriculture-dominated portions of their native range, selecting habitats with dense and permanent cover near crops, meadows and other herbaceous habitats used for feeding [START_REF] Chapman | Sylvilagus floridanus[END_REF][START_REF] Swihart | Habitat feature influencing use of farmstead shelterbelts by the Eastern Cottontail Sylvilagus floridanus[END_REF]. These characteristics are common in the Po Plain and may account for the increase in abundance distribution of the species during the last few decades [START_REF] Angelici | Mammalia II: Erinaceomorpha, Soricomorpha, Lagomorpha, Rodentia[END_REF][START_REF] Vidus-Rosin | Introduced Lagomorphs as a threat to "native" Lagomorphs: The case of the Eastern Cottontail (Sylvilagus floridanus) in northern Italy[END_REF]. In 20 years, the cottontail has expanded its range into north-western Italy [START_REF] Vidus-Rosin | Introduced Lagomorphs as a threat to "native" Lagomorphs: The case of the Eastern Cottontail (Sylvilagus floridanus) in northern Italy[END_REF]. Populations in the province of Alessandria, for instance, have increased from a mean of 4.3 animals per 100 ha with peaks of 25-27 animals per 100 ha [START_REF] Silvano | Distribution of the Eastern Cottontail Sylvilagus floridanus) in the province of Alessandria[END_REF] to 10-30 animals per 100 ha and a peak of 110 animals per 100 ha (Bertolino et al. unpublished data). Given the remarkable adaptability of this species, a further expansion eastwards in North Italy is likely to occur in the years to come and densities could probably increase further [START_REF] Trent | Home range and survival of cottontail rabbits in south-western Wisconsin[END_REF][START_REF] Chapman | Sylvilagus floridanus[END_REF].

The cottontail is already widespread in northern Italy and any eradication programme should be presently considered unfeasible. Within this context, possible interactions with native hares should be carefully evaluated. Although first results suggest that the two species can coexist [START_REF] Vidus-Rosin | Density and habitat requirements of sympatric hares and cottontails in northern Italy[END_REF], 2010;Bertolino et al. unpublished data), cottontail habitat requirement is partially similar to hares. Thus, further research is necessary to exclude a potential competitive advantage of the cottontail species in respect to the hare, under different environmental conditions.

Hare populations in agricultural landscapes are often to be found at low-density levels, and managers should consider management alternatives to improve habitat quality and increase population size. In this context, Fig. 6 Independent contribution for each of the ten variables selected in the MAM. Independent contribution was calculated using hierarchical partitioning (see text) and is expressed as the percentage of the independent contribution of each variable in explaining the total variance conservation measures have to improve native populations of hares while avoiding the introduction or expansion of non-native species. Our results may help to better plan landscape and habitat management for hare conservation. For instance, large meadows that usually favour hares [START_REF] Hutchings | The current status of the brown hare (Lepus europaeus) in Britain[END_REF] turned out to have no or negative effects on cottontail presence. These habitats should therefore be the subject to improve farmland biodiversity in the areas where cottontails are present.

Fig. 1

 1 Fig. 1 Location of the Piedmont Region in Italy and its partition into eight provinces (AL Alessandria, AT Asti, BI Biella, CN Cuneo, NO Novara, TO Torino, VC Vercelli, VB Verbania). The range of cottontail in 2002 is represented in grey. The map was produced joining the 467 municipalities where the species was present

Fig. 5 Fig. 3

 53 Fig. 5 ROC plot for the final model after the tenfold cross-validation. The dashed lines show the ROC curves calculated for each random partition of the cross-validation. The thick bold line shows the average ROC curve calculated across all the random partitions

  . This

	Table 1 Summary of explanatory variables (and their	Variables	Description
	abbreviations) assessed for
	every municipality inside the	Hydrography	Hydrographic network size per surface unit (m/ha)
	study area	Urban %	% of the surface covered by urban areas
		Wetlands %	% of the surface covered by wetlands
		Woodlands %	% of the surface covered by woodlands, orchards and tree plantations
		Meadows %	% of the surface covered by meadows
		Pastures %	% of the surface covered by pastures
		Crops %	% of the surface covered by crop fields (rice fields excluded)
		PASTsize	Mean size of pastures (ha)
		PASTedge	Amount of edge (perimeter of patches) of pastures relative to the area (m/ha)
		MEADsize	Mean size of meadows (ha)
		MEADedge	Amount of edge of meadows relative to the area (m/ha)
		CROPsize	Mean size of crop (rice excluded) patches (ha)
		CROPedge	Amount of edge of crops (rice excluded) relative to the area (m/ha)
		WOODsize	Mean size of woodland patches (ha)
		WOODedge	Amount of edge of woodlands relative to the area (m/ha)
		TOTsize	Mean size of all patches (ha)
		TOTedge	Amount of edge of all patches relative to the area (m)

Table 2

 2 Minimum adequate model (MAM), based on stepwise AIC selection (AIC=647.41, D 2 =0.52)

	Variable	Coefficient	df	SE	z	P
	CROPedge	6.79	1	2.95	2.30	0.02
	Hydrography	0.14	1	0.01	10.53	0.00
	MEADedge	2.99	1	1.45	2.06	0.04
	Meadows	-8.10	1	2.76	-2.94	0.00
	MEADsize	0.01	1	0.00	2.09	0.04
	Pastures	-29.95	1	9.40	-3.19	0.00
	PASTsize	0.01	1	0.01	1.78	0.07
	TOTsize	0.01	1	0.00	2.96	0.00
	Wetlands	-5.74	1	2.72	-2.11	0.04
	Woodlands	-5.04	1	0.56	-8.97	0.00

Acknowledgements We are grateful to the Osservatorio Faunistico del Piemonte, the provincial Wildlife Services and all the institutions and people that gave us information on the distribution of the cottontail in Piedmont. Two anonymous referees greatly improved a previous draft of the manuscript.