Fabrice Devaux
email: fabrice.devaux@univ-fcomte.fr

Jean-Pierre Huignard
email: jphuignard@free.fr

François Ramaz
email: francois.ramaz@espci.fr

Modelization and optimized speckle detection scheme in photorefractive self-referenced acousto-optic imaging

Keywords: Photorefractive optics, Imaging through turbid media, Speckle interferometry, Acousto-optics

A photorefractive BSO single crystal can be used for axially resolved acousto-optic imaging of thick scattering media in absence of a reference beam. This configuration renders the experimental setup easier to realize for imaging through thick scattering media with an improved optical etendue.

We present here a model and simulations that explains these results. It is based on the spatial heterogeneity of the speckle pattern incident on the crystal. Optimization of the detector position and of the speckle grain size is confirmed by the model.

INTRODUCTION

Optical imaging through thick biological tissues of several centimeters requires indirect methods in order to retrieve a millimeter resolution over the entire volume. The major drawback comes from scattering that generates multiple optical paths, so that the origin of the photons that propagate through the medium is lost. Many attempts have been made over the past few decades based on flux techniques such as diffuse optical tomography [START_REF] Arridge | [END_REF]2], but the extraction of the signal is not straightforward, since it requires solving an inverse scattering problem. A more recent concept couples light with an ultrasonic beam within the medium [3]. Dealing with the acousto-optic (AO) effect, the ultrasound (US) locally modulates the optical path, and thus the phase of light crossing the US field. This phase modulation mainly creates two weak optical sidebands (< 1%), shifted from the career frequency ω L of the US frequency +/ -ω U S . The magnitude of these so-called tagged photons depends on the optical properties of the medium, such as absorption and scattering. A 3D image of the optical properties of the medium is thus performed with a scan of the US beam. Mainly two approaches are implemented to detect the optical sidebands. The first one consists in a highly spectrally resolved flux measurement with a confocal Fabry-Perot interferometer [4], or in spectral holeburning filtering at liquid helium temperature [5]. A second approach uses a coherent detection of the tagged photons, initially performed by digital off-axis holography with camera detection (CCD, CMOS) [6,7]. This coherent detection requires accurate sampling of the speckle and the use of a fast camera, since the speckle decorrelates in living tissues on a millisecond time scale or less [8]. For these reasons, wavefront adaptive setups have been implemented with photorefractive (PR) crystals and large-area single detectors, operating much faster than a camera [9]. In dynamic PR holography, a refractive index grating is recorded within the volume of the crystal, resulting from the interference between the tagged photons and a plane wave off-axis reference beam. Consequently, the reference diffracts a spatial replica of the tagged photons' wavefront, which is time averaged according to the relative response time τ P R of the crystal (1-20ms), depending on light irradiance (a few 100mW/cm 2 ) [10]. The use of a reference beam to record the hologram generally brings an important coherent background noise (beam fanning), which limits the signal-to-noise ratio (SNR) of detection, since the AO signal is less than two orders of magnitude compared to the DC-flux components when dynamic axial resolution along the ultrasonic beam is considered. In this paper, we will show that AO imaging with PR detection can also be performed without applying a reference beam onto the crystal. This type of interaction was demonstrated by Kamshilin et al. [11] for application of PR correlation filtering to the detection of vibrating rough surfaces or vibrating multimode fibers. The proof of principle for AO imaging is illustrated by Benoit à la Guillaume et al. [12] with a dynamic 2D image of ink inclusions embedded within 40mm of a multiply scattering phantom. The PR crystal used is the sillenite type Bi 12 SiO 20 (BSO), known for its high photosensitivity in conventional holographic experiments in the blue-green spectral region. We outline that this configuration is much simpler to implement after having optimized the light-collection efficiency and the speckle size incident on the crystal. The experimental setup detailed in Ref. [12] is shown in Fig. 1. We apply an excitation burst (carrier at 2.23 M Hz) on the US transducer to get a good axial resolution of 1.35 mm while the lateral resolution, given by the tranverse profile of the US, is 2 mm. The used gel phantom has a reduced scattering coefficient µ ′ s of about 10 cm -1 with two included china ink cylinders to simulate local absorbers. The detection of the AO signal is performed with a large area photodiode of 3.6×3.6mm 2 placed behind the optical system (L 3 L 4 ) that improves the collection of the speckle light. As analyzed in [12], the self referenced interaction can be interpreted as follows: in absence of US the incident speckled wavefront creates in the PR BSO a photoinduced space charge field which modulates the refractive index through the diffusion of carriers at no applied electric field. This complex phase volume hologram results from the interference of the incident speckle field with self-diffracted beam fanning noise fields. The speckle pattern photo-induces micrometer size random refractive index dots within the BSO volume. The operating conditions of the self-referenced PR holography for acousto-optic imaging are explained as follows: when the US transducer is off, the incident speckle field is perfectly phase matched with the recorded phase volume hologram. It thus results an optimum coherent diffraction on the detector. When the US transducer is now operating over short pulse duration the incident speckle structure is changed and it diffracts a wavefront which does not perfectly match with the complex phase hologram recorded in the volume of the BSO. In such conditions the diffracted intensity on the detector is reduced or enhanced as shown in [12] and permits an acousto-optic imaging of a phantom embedded in the scattering medium. It was also shown experimentally in [12] that in the case of self referenced holographic configuration an optimum signal detection is obtained when the two following conditions are fulfilled: first, the speckle grain size recorded in the PR crystal must be of the order of the Debye screening length. Second, the position of the detector in the transverse plane must be optimized due to the spatial frequency response of the photorefractive effect when speckle is involved.

The objective of the paper is to develop hereafter an analytic model and numerical simulations which confirm these optimum conditions for signal detection in the conditions of self referenced speckle beam holography with a PR crystal operation at zero electric field.

THEORETICAL MODEL AND NUMERICAL SIMULATIONS

According to the experiment described in [12], we first define the incident speckle field (i.e. the reference speckle field) and its intensity inside the BSO crystal respectively as: S r (r) and I r (r) = |S r (r)| 2 , where r is the position vector inside the PR medium. This random distribution of coherent light intensity I r photoinduces an internal space charge electric field in the whole volume of the PR crystal noted E sc (r). When the classical approximations are made in the PR model where electrons are considered as the single specy of charge carriers [START_REF] Yeh | Introduction to Photorefractive Nonlinear Optics[END_REF], the Ohm's law governing the current density and the rate equations of the PR charge carriers lead to the following differential equation [START_REF] Daniel | [END_REF]15] :

E sc (r) - 1 k 2 0 ∇ (∇ • E sc (r)) = - k B T e ∇ ln [I r (r) + I d ] , (1)
where ∇ is the nabla operator, k B is the Boltzmann constant, T is the temperature, e is the elementary charge, k 0 is the Debye screening wave number (typically for BSO k 0 ≃ 4.65µm -1 [16]) and I d is a uniform background illumination responsible of dark current photoconductivity. Using the Fourier transform properties, Eq. ( 1) becomes:

FT {E sc (r)} = - k B T e ik 1 + |k| 2 k 2 0 FT {ln [I r (r) + I d ]} , (2)
where FT denotes the 3D Fourier transform and k is the wavevector. Assuming that the variation of the PR variables along the propagation axis can be neglected and a twodimensional illumination where I r (x, y) ≪ I d , Eq. ( 2) may be rewritten as:

       Ẽscx (k t ) = -iαkx 1+ |k t | 2 k 2 0 Ĩr (k t ) Ẽscy (k t ) = -iαky 1+ |k t | 2 k 2 0 Ĩr (k t ) (3)
where Ẽscx (k t ), Ẽscy (k t ) and Ĩr (k t ) are respectively the 2D Fourier transforms, calculated in the (x, y) plane, of the transverse components of the space charge field and of the speckle intensity. k t = k x + k y is the transverse component of the wavevector and α = k B T I d e is a constant. These equations show that each internal space charge field component corresponds to the product of an odd function, presenting extreme values for ±k 0 , with Ĩr (k t ) that corresponds to the autocorrelation function of the speckle intensity.

Following the analytical expression of the internal space charge field given by Eq. 3, numerical calculations are performed to illustrate the electric field components. Figures 2(a) and 2(b) represent respectively the intensity of a typical speckle pattern in the (x, y) spatial domain and its spectrum in the (k x , k y ) spatial frequency domain. Figure 2(c) depicts the modulus of the Fourier transform of the speckle intensity | Ĩr (k t )|. The calculations are performed at the wavelength λ = 0.532µm and with the sampling steps ∆x = ∆y = λ 4 in the (x, y) plane. The speckles are simulated using a complex field amplitude with a real part and an imaginary part that obey to a gaussian random distribution. The 1 e 2 half width σ s of the speckle grains is adjusted by filtering, in the Fourier domain, the spectrum of the speckle with a gaussian filter. We verified that the statistical property of the intensity fluctuations are in agreement with the thermal statistic of speckles [START_REF] Goodman | Laser Speckle and Related Phenomena[END_REF] where the probability 3, we assume that the pedestal of the autocorrelation function of the speckle intensity has a gaussian shape with a 1 e 2 half width σ k related to the size σ s of the speckle grains. We can easily demonstrate that the amplitude of the space charge field components exibit maxima for the spatial frequencies k x = ±σ k and k x = ±k 0 when the size of the speckle grains is respectively greater and smaller (or equal) than the Debye screening length.

The crystallographic orientation of the BSO crystal used in [12] is shown in Fig. 1. The light is incident on the plane (x, y) where the x axis is perpendicular to the plane (0,0,1)

and the y axis is perpendicular to the plane (1,-1,0). In [12], measurements of the angular dependency of the PR signal were performed with a photodiode placed in the far field of the PR crystal and that scanned the PR signal along the k x spatial frequencies where k y ≈ 0.

For this spatial frequency range, the internal space charge field has only a component along the x axis (i.e. the Ẽscx (k t ) component) and two birefringence axes appear, one parallel to the laboratory x axis, and one parallel to the y axis. The refractive indices found by the light polarized along these axes and propagating in the directions such that k y = 0, are n 0 + ∆n x and n 0 + ∆n y where n 0 is the linear refractive index of the BSO crystal and [13]

   ∆n x (k x ) = 0 ∆n y (k x ) = - 1 2 r 41 n 3 0 E scx (k x ). (4)
As scattered light is strongly depolarized, we emphasize that only the component of the speckle field polarized along the y axis reads the modulation of the refractive index ∆n y induced by the component Ẽscx of the internal space charge field. These results suggest that the x component of the speckle field provides no contribution to the PR signal. Then, in the experimental setup, the use of a polarizer oriented along the y axis at the output of the BSO crystal could improve the signal to noise ratio of the PR signal. For the y axis polarization considered in this paper, BSO exhibits both photoinduced birefringence and large optical activity at the green wavelength. In such conditions wave progation of an incident depolarized speckle beam is complex and its detailed analysis may be out of the scope of this paper. When taking account of the optical activity for linearly polarized plane waves, it leads to a reduction of the BSO holographic efficiency [16].

Because the modulation of the refractive index is weak (10 -5 to 10 -6 ), the PR phase hologram induced by the speckle and read by the y polarized light can be expressed by :

H y (x, y) = exp -i 2π∆n y L λ ≃ 1 -i 2π∆n y L λ ( 5
)
where L is the thickness of the hologram. Then, the y component of the output speckle field is given by : S out ry (x, y) = S in ry (x, y)H y (x, y).

Using the Fourier transform expression of the space charge field component given by the Eq. ( 3) and the expression of the phase function of the hologram given by the Eq. 5, the (k x , k y = 0) spatial frequency spectrum of the output speckle is given by:

Sout ry (k x ) = Sin ry (k x ) ⊗ δ(k x ) + βkx 1+ k 2 x k 2 0 Ĩr (k x ) = Sin ry (k x ) + Sin ry (k x ) ⊗ βkx 1+ k 2 x k 2 0 Ĩr (k x ), (7)
where ⊗ denotes the convolution product and δ(k x ) is the Dirac function centered at the zero spatial frequency and β a new constant including all the constants of the previous equations. We can see that the amplitude of the output spectrum of the speckle diffracted by the photorefractive hologram is made of the incident speckle field (i.e. the zero order of diffraction), and a second term corresponding to the diffraction of the incident speckle by its own hologram. From an algebric point of view, this term is of the same sign as k x and explains the origin of the reduction or the enhancement of the PR signal with respect to the position of the detector along the k x axis. We stress out that this term is null in the mean direction of the incident light (i.e. the zero spatial frequency) so that no PR signal can be recorded in this direction. This property can be physically explained if considering the speckle field within the crystal as a superposition of plane waves whose incidences are included within the aperture of the lenses that shine the cristal : a multi two-wave mixing process is performed along the mean propagation axis z, leading by symmetry to a photorefractive gain which can be positive or negative, depending on the sign of the angle of the plane wave with respect to z. As a consequence, the superposition of all these contributions cancels the effect along z.

Because of the small phase modulation induced by the hologram, the diffracted speckle is very similar to the reference speckle and the contribution of the PR signal in the output speckle is very weak. Then, the PR signal is obtained by subtracting the background level of light such that:

I P R (k x ) = Sout ry (k x ) 2 -Sin ry (k x ) 2 . (8)
In order to study the diffraction efficiency of the PR hologram, numerical simulations are performed with different sizes of the speckle grains varying in the range [0.1µm -3µm] and speckles with a mean intensity kept constant. Figures 4(a) and 4(b) show, for different sizes of the speckle grains, the variation of the modulus of the x component of the space charge field and of the PR signal with respect to the normalized spatial frequency kx k 0 . In agreement with Eq. 3, the amplitude and the position of the maximum of the space charge field depend on the size of speckle grains. The best efficiency is obtained at the spatial frequencies k x = ±k 0 when the size of the speckle grains corresponds approximately to the Debye screening length (σ s ≈ 1 k 0 = 0.21µm). The experimental results reported in Fig. 2 of [12] show the angular dependency of the PR signal measured with a large area photodiode of 3.6×3.6 mm 2 placed at a distance of 25mm of the PR crystal (in this case the optical system L 3 L 4 in the experimental setup is removed). In order to fit with these experimental conditions, the PR signal curves presented in fig. 4(b) are the result of the convolution product between the PR signal calculated with Eq. 8 and a square function, with a width ∆k D = 1.7 µm -1 , that simulates the 0.14 rad angular aperture of the photodiode. Moreover, to smooth out the intensity fluctuations of the PR signal due to the small number of speckle grains in the area of the simulated photodiode, we average the calculated PR signal over more than 10 numerical realizations to obtain these curves. Finally, the PR signal is plotted 9 simultaneously with respect to the normalized spatial frequency kx k 0 and with θ = asin kxλ 2π , the external angle (outside the crystal). In fig. 4(b), the shape of the PR signal is in a good agreement with the results presented in [12] where amplitude of the PR signal was maximized for a position of the detector corresponding to an external angle θ = ±17 • when the size of speckle grains was on the order of l D . These results also suggest that experimental parameters used in [12] were not optimized. Indeed, fig. 4(b) shows that the PR signal is maximized for larger speckle grains (0.4µm ≤ σ s ≤ 0.9µm) while the amplitude of the space charge field is maximized when the size of the speckle grains fits with the Debye screening length. Secondly, the efficiency of the PR signal is optimized when the detection is performed in a direction corresponding to a spatial frequency smaller than k 0 . Physically, this result differs from the two-wave mixing experiments where the photorefractive gain is optimized when the grating wave number K equals k 0 [16]. Indeed, in this case the terms related to the speckle in Eq. 7 will be replaced by Dirac functions δ(k x ±K) that leads to a diffraction efficiency which is governed only by the impulse response of the PR crystal βkx

1+ k 2 x k 2 0
and maximized when K = k 0 . With speckles, convolution product in Eq. 7 implies that diffraction efficiency depends on the impulse response of the PR crystal as well as the spatial frequency width of the speckle σ k = 1 σs . Finally, finite size of the detector is also a factor to be taken into account to optimize measurement of the PR signal.

Figure 4c shows the amplitude of the PR signal convolved with a square function with a width ∆k D and calculated as a function of the normalized spatial frequency kx k 0 and the normalized half width of the gaussian envelop of the speckle spectrum σ k k 0 . The black dotted curve depicts the position of the detector that gives the higher PR signal as a fonction of σ k . Vertical and horizontal black dotted lines make visible the couple of parameters (σ k = 0.20k 0 , k x = 0.35k 0 ) optimizing the PR signal. We can stress out that this optimized couple of parameters depends on the angular aperture of the detector. Now, in order to modelize the ultrasound modulated optical tomography by selfreferenced photorefractive holography reported in [12], let us consider the case where an incident speckle tagged by a short US pulse, noted S in t , is diffracted by the hologram induced by the reference speckle S r . Because the interaction between the short US pulse with the light propagating in the scattering medium leads to changes of the speckle pattern, the tagged speckle is transiently not perfectly matched with the complex phase hologram recorded in the whole volume of the BSO crystal. As a consequence, the level of the PR signal varies during the time of the US pulse.

In our simulations, the tagged speckle is modelised by the coherent superposition of the reference speckle S r with a "noise" noted S n (i.e. a fully decorrelated speckle) so that the amplitude of the tagged speckle may be written in the (x, y) domain as: where γ represents the ratio of the incident light interacting with the US pulse in the scatering medium. Then, the intensity of the tagged speckle is given by:

S in t = √ γS n + 1 -γS r , (9)
I in t = γI n + (1 -γ)I r + 2 γ(1 -γ)Re{S n S * r }. (10)
Because the noise speckle S n is fully decorrelated with the reference speckle S r , the mean value of the interference term between the two speckles is null. Figure 5 shows the typical result of the interference between two speckles fully decorrelated. Therefore, this term can be neglected when γ ≪ 1 so that the intensity of the tagged speckle is assumed to be the incoherent sum of the intensities of the two speckles. According to this assumption, the detected intensity of the tagged speckle that reads the PR hologram after US interaction can be written, in the spatial frequency domain, as : Experimentally, the position of the photodiode is adjusted along the k x dimension to maximize the PR signal. Then, the measured AO signal is directly the intensity of the spectrum of the tagged speckle at this spatial frequency and the DC background corresponds to the intensity of the untagged speckle measured at the same spatial frequency. Hence, we will define the relative AO signal as :

| Sout t | 2 ≃ γ| Sout n | 2 + (1 -γ)| Sout r | 2 . (11)
R AO = ∆k D | Sout t | 2 d 2 k -∆k D | Sout r | 2 d 2 k ∆k D | Sout r | 2 d 2 k (12)
where ∆k D denotes the integration of the intensity of the spectrum of the diffracted speckle field over the spatial frequency bandwidth ∆k D of the photodiode. For numerical calculation, the mean radius of the speckle grains in the (x, y) is adjusted to 0.21 µm (i.e. the Debye screening length). Figure 6(a) represents the typical shape of the intensity of the spatial frequency spectrum of a tagged speckle and the two yellow squares symbolise the optimized position and the area of the detector centered on the spatial frequencies k x = ± k 0 2 (see Fig. 4b). For these two positions of the photodiode and using Eq. 12, we calculate the variation of the relative AO signals as a function of the level γ of the "noise" in the tagged speckle.

Assuming that the scattered light totally filled the volume of the scattering medium, the maximum value taken by γ is determined very approximately to 10 -3 by calculating the ratio between the focusing volume of the US pulse (∼ 5mm 3 ) and the total volume of the scattering medium (< 40 3 mm 3 ). Figure 6(b) gives the variation of the relative AO signals, respectively noted R AO+ and R AO-when the detector is centered on + k 0 2 or -k 0 2 , as a function of the level of noise in the tagged speckle. Because of the small number of speckle grains in the area of the simulated detector, the curves correspond to the AO signals averaged over few tens of numerical realizations. In agreement with the shape of the PR signal given by the fig. 4(b) and the Eq. 11, the relative AO signal measured around + k 0 2 and -k 0 2 respectively decreases and increases linearly, when the level of noise increases. All these results give a satisfactory quantitative interpretation of the experimental results reported in [12]. However, the developed model mainly outlines the optimal values of the speckle grain size and of the detector position and, for the moment, does not permit to quantify the holographic efficiency and the AO signal dynamic range. In future work, several experimental parameters, such as fraction of tagged photons, beams polarization, light collection efficiency on the detector need to be carefully considered in order to refine amplitude consistency between the model and experiment.

CONCLUSION

The self referenced wavefront acousto-optic imaging technique has been analysed in the condition where a single speckle field is projected on a PR crystal. In these experiments, no reference beam is applied and the incident speckle creates a complex space charge field inducing a phase volume hologram in the BSO. The model leads to determine the parameters controling the efficiency of the hologram. We show that the PR signal efficiency is maximized for a spatial frequency that depends on the crystal Debye screening wave number, the width of the spatial frequency spectrum of the speckle but also on the angular aperture of the detector. It results that there is an optimum position of the photodetector and an optimum size of the speckle grains to detect the speckle pattern in presence of the tagged photons in an acousto-optic imaging experiment. The numerical simulations agree well with the experimental results previously reported. The method is robust and easy to implement in comparison with conventional holography. It can be applied with other types of nonlinear media in particular operating in the near IR region for acousto-optic imaging in highly scattering biological tissues and in turbid media and vibrometry. This work was supported by the LABEX WIFI (Laboratory of Excellence within the

FIG. 1 .

 1 Photodiode

FIG. 2 .

 2 FIG. 2. (a) Intensity of a typical speckle pattern in the (x, y) spatial domain where the 1 e 2 half width of the speckle grains is σ s = 0.21µm. (b) Intensity of the spatial frequency spectrum of the speckle field. (c) Modulus of the Fourier transform of the speckle intensity. In the spatial frequency domain, the axes are normalized by k 0 , the Debye screening wave number.

FIG. 3 .

 3 FIG. 3. (a) and (c) Modulus of the transverse components of the space charge electric field in the (k x , k y ) spatial frequency domain. (b) and (d) Modulus the transverse components of the space charge electric field in the (x, y) spatial domain.

FIG. 4 .

 4 FIG. 4. Variation of the modulus of the x component of the space charge field (a) and of the PR signal (b) as a function of the normalized spatial frequency kx k 0 (k y = 0) and of the mean size σ s of the speckle grains. (c) Normalized PR signal as a function of the normalized spatial frequency kx k 0 and the normalized half width of the gaussian envelop of the speckle spectrum σ k k 0 . The black dotted curve depicts the optimized position of the detector that gives the higher PR signal as a fonction of σ k .

10 FIG. 5 .

 105 FIG. 5. (a) and (b) intensities in the (x, y) domain of the reference speckle and the fully decorrelated noise speckle. (c) Real part of the interference term between the reference and the noise speckles.

4 FIG. 6 .

 46 FIG. 6. (a) Intensity of the spatial frequency spectrum of a tagged speckle. The two yellow squares represent the detector area in the (k x , k y ) domain when it is centered respectively on the spatial frequencies k x = ± k 0 2 . (b) Amplitude of the AO signals as a function of the level γ of noise in the tagged speckle for the two positions of the detector.


M A Veganzones

J Frontera-Pons

J Chanussot

J P Ovarlez

ON THE USE OF THE HOTELLING'S T 2 STATISTIC FOR THE HIERARCHICAL CLUSTERING OF HYPERSPECTRAL DATA

Keywords: hypespectral imaging, hierarchical clustering, Fixed Point estimates

In this work we propose a hierarchical clustering methodology for hyperspectral data based on the Hotelling's T 2 statistic. For each hypespectral sample data, the statistical sample mean is calculated using a window-based neighborhood. Then, the pairwise similarities between any two hyperspectral samples are computed based on the Hotelling's T 2 statistic. This statistic assumes a Gaussian distribution of the data while hyperspectral data have been observed to be long tailed distributed. In order to improve the statistic robustness we use the Fixed Point estimates, and compare them to the classical sample mean estimator. The similarities are then used to hierarchically cluster the hyperspectral data. We give some preliminary qualitative results of the proposed approach over the Indian Pines hyperspectral scene. Results show that the use of the Fixed Point estimator does not significantly affect the clustering results. Further work will be focused on the use of the robust Hotelling statistic.

INTRODUCTION

Hyperspectral data have been observed not to be multivariate normal but long tailed distributed [START_REF] Manolakis | Detection algorithms for hyperspectral imaging applications[END_REF]. In order to take into account these features, the class of elliptical contoured distributions (ECD) is considered to describe clutter statistical behavior. It provides a multivariate location-scatter family of distributions that primarily serves as heavy tailed alternative to the multivariate normal model. A m-dimensional random vector y = [y 1 y 2 ...y m ] T with mean µ and scatter matrix Σ has an elliptical distribution if its probability density function (PDF) has the following form:

f y (y) = |Σ| -1 h m ((y -µ) T Σ -1 (y -µ)), (1)
where T denotes the transpose operator and h m (.) is any function such as [START_REF] Manolakis | Detection algorithms for hyperspectral imaging applications[END_REF] defines a PDF in R m .

This work was supported by the Delegation Generale de l'Armement (DGA) under contract PDE12C2009.

The ECD class includes a large number of distributions, notably the Gaussian distribution, multivariate t distribution, K-distribution or multivariate Cauchy. Thus it allows for heterogeneity of the background power with the texture. Although non-Gaussian models are assumed for background modelling and test design, the parameters estimation is still performed using classical Gaussian based estimators; as for the covariance matrix, it is generally determined according to the Sample Covariance Matrix (SCM), MSCM = N i=1 y i y T i . We use here the Fixed Point (FP) estimators [START_REF] Tyler | A distribution-free m-estimator of multivariate scatter[END_REF] as a robust alternative for background statistical characterization.

In [START_REF] Formont | Statistical classification for heterogeneous polarimetric sar images[END_REF], authors propose a general approach for highresolution polarimetric SAR (POLSAR) data classification in heterogeneous clutter, based on a statistical test of equality of covariance matrices. In hyperspectral domain the information provided by the mean is of high relevance. Thus, we propose here a similar approach for hyperspectral data using the Hotelling's T 2 statistic as a statistical test of equality of means [START_REF] Bilodeau | Theory of multivariate statistics[END_REF]. We define a pairwise similarity between two pixels by using the output of the Hotelling's test for the means of the pixels calculated using a sliding window. As a preliminary test, we segment the hyperspectral scene by means of a hierarchical clustering algorithm [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF] whose input is given by the previously calculated pairwise distances. We test the proposed approach in the Indian Pines hyperspectral scene.

THE FIXED POINT ESTIMATES

The Fixed Point estimates proposed by Tyler in [START_REF] Tyler | A distribution-free m-estimator of multivariate scatter[END_REF] belong to the wider class of robust M-estimates [START_REF] Huber | Robust estimation of a location parameter[END_REF]. According to the FP approach, the joint estimation of the scatter matrix, M, and the mean vector, µ, leads to [START_REF] Tyler | A distribution-free m-estimator of multivariate scatter[END_REF]:

MFP = m N N n=1 (y n -μ)(y n -μ) H (y n -μ) H M-1 FP (y n -μ) (2)
and

μ = N n=1 y n (y n -μ) T M-1 FP (y n -μ) 1/2 N n=1 1 (y n -μ) T M-1 FP (y n -μ) 1/2 (3)
where N denotes the number of secondary data and y n the vector under observation. For the matrix estimate, existence and uniqueness have been established in [START_REF] Pascal | Covariance structure maximum-likelihood estimates in compound gaussian noise: Existence and algorithm analysis[END_REF]. Although the proof for simultaneous scatter and location estimates is still an open question, they have been found to be useful and reliable for the estimation of ECD parameters because of its easy implementation. They are specified by implicit equations and can be easily computed using a recursive algorithm. We refer to [START_REF] Pascal | Performance analysis of covariance matrix estimates in impulsive noise[END_REF] for a detailed performance analysis of the FP covariance matrix estimate. The main results of the statistical properties of the MFP are summarized: MFP is a consistent and unbiased estimate of M; its asymptotic distribution is Gaussian and is the same as the asymptotic distribution of a Wishart matrix with mN/(m + 1) degrees of freedom. They provide a robust alternative to the classical estimators and its definition on hyperspectral detection framework was introduced in [9].

HOTELLING'S T 2 HYPOTHESIS TESTS

Let X ∼ ECD(µ 1 , Σ 1 , h m,1 ) and Y ∼ ECD(µ 2 , Σ 2 , h m,2 )
be two independent random vectors, elliptically distributed. We intend to decide if they belong to the same class comparing their mean vectors µ 1 and µ 2 . The classification problem can be formulated as a binary hypothesis test. Symbolically, we aim to distinguish between:

H 0 : µ 1 = µ 2 H 1 : µ 1 = µ 2
where the decision on the hypothesis test is made according to some statistical based criteria. We propose a similarity function based on the parameters estimation.

Hotelling's T 2 statistic. Assuming x 1 , ..., x N1 and y 1 , ..., y N2 are i.i.d. samples from X ∼ N m (µ 1 , Σ) and Y ∼ N m (µ 2 , Σ) respectively. And the mean vector computed as, μx =1 N1 N1 i=1 x i and μy = 1 N2 N2 i=1 y i and the covariance matrix estimated according to:

W = N 1 i=1 (xi -μx )(xi -μx ) T + N 2 i=1 (yi -μy )(yi -μy ) T N1 + N2 -2
which correspond to the unbiased pool estimate. Then the Hotelling T 2 statistic test is defined as:

t 2 = N 1 N 2 N 1 + N 2 ( μx -μy ) T W -1 ( μx -μy ). ( 4
)
This expression is easily derived from the likelihood function of a multivariate normal model when H 0 is assumed, see e.g. [START_REF] Bilodeau | Theory of multivariate statistics[END_REF]. The distribution of the test under the null hypothesis can be related to the F-distribution according to

N 1 + N 2 -m -1 (N 1 + N 2 -2)m t 2 ∼ F (m, N 1 + N 2 -1 -m).
A different situation arises when samples do not fit on a multivariate gaussian distribution. We can improve the robustness of the test replacing the mean and matrix estimates for some robust alternative (e.g. the FP) if its asymptotic distribution is gaussian. The equality of the covariances matrices is a simplification that will certainly modify the test, as the resulting matrix estimate will not correspond to a Wishart distribution

(1/N 1 Σ 1 +1/N 2 Σ 2
). The results found on statistics literature for robust Hotelling and for both matrix and mean comparison will be further investigated.

STATISTICAL HIERARCHICAL CLUSTERING

A hierarchical clustering algorithm returns a hierarchy of clusters built by merging smaller components into bigger clusters (agglomerative clustering) or by splitting the whole image into smaller regions (divisive clustering). In order to do that, a pairwise similarity function is necessary, i.e. the angular distance, to compare each pair of components. Then, a linkage function indicates which of any possible pair of components is merged to (or split from) a bigger cluster. A common linkage is the Unweighted Pair Group Method with Arithmetic Mean (UPGMA) method which consists in averaging all distances between the samples in each node, that is, the mean distance between elements of each cluster. Given two nodes A and B, the UPGMA linkage distance between the nodes, denoted l (A, B) is given by:

l (A, B) = 1 |A| |B| x∈A y∈B d (x, y) (5)
where |•| denotes the cardinality of a set and d (x, y) is any pairwise distance. In this work we propose to use the Hotelling's T 2 statistic (4) as the pairwise distance. The result is a hierarchy of clusters that can be represented as a dendrogram. A common criterion to stop merging (splitting) is to a priori set the number c of clusters one is looking for.

EXPERIMENTS AND RESULTS

Indian Pines scene

The Indian Pines scene 1 was gathered by airborne AVIRIS sensor over North-western Indiana and consists of 145 × 145 pixels and 224 spectral reflectance bands in the wavelength range 0.4-2.5 µm . We have reduced the number of bands to 200 by removing bands covering the region of water absorption: [104-108], [150-163], 220. Since the scene is taken in June some of the crops present, corn, soybeans, are in early stages of growth with less than 5% coverage. The ground Fig. 1. Reduced groundtruth for the Indian Pines scene truth available is designated into sixteen classes with variable number of samples for each class. In order to simplify the problem in this preliminary stage, we have grouped some similar classes (see Fig. 1). The grouped classes are Corn (corn, corn-notill and corn-mintill), Grass (grass-pasture, grasspasture-mowed and grass-trees) and Soybean (soybean-notill, soybean-mintill and soybean-clean).

Experimental design and results

We compare the use of the classic estimator for the sample mean and SCM to the FP estimator for the clustering of the Indian Pines scene using the Hotelling's T 2 statistic (4) and the hierarchical clustering with UPGMA linkage function [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF]. For each pixel in the image we calculate its mean and covariance matrix by the classic estimator and the FP estimator using a 15 × 15 sliding window in order to ensure enough samples for the estimation stage (N > 2m). Those pixels on the borders of the image for which the window is incomplete are discarded resulting in a sample matrix of size 131 × 131.

We have performed two clustering experiments, a first one using the entire dataset (17161 samples) and a second one using only those data samples for which there is available groundtruth (9343 samples). For the latter, the statistics have been estimated using the entire dataset, but the linkage given by UPGMA algorithm have been built using only the samples from the groundtruth subset. Figure 2 shows the clustering results for the former and figure 3 does it for the latter. In order to visualize the results we have assigned to each cluster the RGB color resulting of averaging the groundtruth RGB color (see Figure 1) of the pixels of each cluster. Thus, a cluster with a color similar to any groundtruth color represents that the cluster is almost contained in a single groundtruth class. Clusters grouping several different groundtruth classes will present a mixed color.

The obtained results show that the use of the Hotelling's T 2 statistic yields to meaningful clusters, specially if we allow some overclustering (see figures 2 and 3 e-j). Moreover, it is possible to visually appreciate a better preformance of the FP estimator for a small number of clusters (see figures 2 and 3 a-d). The use of robust Hotelling's test will be further investigated to deal with the ECD assumption.

CONCLUSIONS

We have presented an approach to statistically characterize the hyperspectral data and then use the Hotelling's T 2 statistic to define a pairwise similarity function. The definition of a pairwise similarity allows to use a hierarchical clustering algorithm to segment the hyperspectral image. We have shown that even when the Hotelling's statistic Gaussian assumption differs from the actual ECD, the clustering results are meaningful. The use of a robust FP estimator shows a better performance than the use of the classic estimator. The use of a robust Hotelling's estimator capable of dealing with non-Gaussian assumption will be further investigated.

Fig. 2 .

 2 Fig. 2. Results of the hierarchical clustering over the entire dataset using the Hotelling T 2 test with the classic estimator (left column) and the Fixed Point estimator (right column). Each row identifies with a different value of the number of clusters, c.

Fig. 3 .

 3 Fig. 3. Results of the hierarchical clustering over the reduced dataset using the Hotelling T 2 test with the classic estimator (left column) and the Fixed Point estimator (right column). Each row identifies with a different value of the number of clusters, c.

https://engineering.purdue.edu/ biehl/MultiSpec/

Acknowledgements M.A. Veganzones was supported by the Delegation Generale de l'Armement (DGA) under contract PDE2012600079.


M A Veganzones

J Frontera-Pons

F Pascal

J.-P Ovarlez

J Chanussot

BINARY PARTITION TREES-BASED ROBUST ADAPTIVE HYPERSPECTRAL RX ANOMALY DETECTION

Keywords: Anomaly detection, RX AD, binary partition trees

The Reed-Xiaoli (RX) is considered as the benchmark algorithm in multidimensional anomaly detection (AD). However, the RX detector performance decreases when the statistical parameters estimation is poor. This could happen when the background is non-homogeneous or the noise independence assumption is not fulfilled. For a better performance, the statistical parameters are estimated locally using a sliding window approach. In this approach, called adaptive RX, a window is centered over the pixel under the test (PUT), so the background mean and covariance statistics are estimated using the data samples lying inside the window's spatial support, named the secondary data. Sometimes, a smaller guard window prevents those pixels close to the PUT to be used, in order to avoid the presence of outliers in the statistical estimation. The size of the window is chosen large enough to ensure the invertibility of the covariance matrix and small enough to justify both spatial and spectral homogeneity. We present here an alternative methodology to select the secondary data for a PUT by means of a binary partition tree (BPT) representation of the image. We test the proposed BPT-based adaptive hyperspectral RX AD algorithm using a real dataset provided by the Target Detection Blind Test project.

INTRODUCTION

Target detection (TD) and anomaly detection (AD) of multidimensional signals have proved to be valuable techniques in a wide range of applications, including search-and-rescue, surveillance, rare minerals detection, mine detection, etc [START_REF] Nasrabadi | Hyperspectral target detection : An overview of current and future challenges[END_REF][START_REF] Manolakis | Detection algorithms in hyperspectral imaging systems: An overview of practical algorithms[END_REF]. TD aims to discover the presence of a specific signal of interest (the target) among a set of signals. TD is based on the Neyman-Pearson (NP) criterion, which maximizes the probability of detection for a given probability of false alarm. AD is a special case of TD in which no a-priori target is provided. Hence, the goal of AD is to detect signals that are anomalous with respect to the background [START_REF] Matteoli | A tutorial overview of anomaly detection in hyperspectral images[END_REF].

The Reed-Xiaoli (RX) AD algorithm [START_REF] Irving | Adaptive multipleband cfar detection of an optical pattern with unknown spectral distribution[END_REF] is considered as the benchmark algorithm in multidimensional AD. However, the RX detector performance strongly relies on the statistical parameters estimation. Hence, when the background is non-homogeneous or the noise independence assumption is not fulfilled, the detector performance is deteriorated [START_REF] Acito | On the cfar property of the rx algorithm in the presence of signaldependent noise in hyperspectral images[END_REF][START_REF] Matteoli | A locally adaptive background density estimator: An evolution for rx-based anomaly detectors[END_REF][START_REF] Frontera-Pons | Robust anomaly detection in hyperspectral imaging[END_REF]. The statistical parameters can be estimated globally or locally. The local version is called adaptive RX, since it adapts to the local properties of the data. The conventional approach to develop an adaptive RX algorithm is by using sliding windows. For a given pixel under test (PUT), a small window with the size of the expected anomalies maximum size, named the guard window, is centred on the PUT. A second large window, usually named outer window, is also centred on the PUT. The pixels inside the outer window, except those lying inside the guard window, compose the secondary data used to estimate the statistical parameters for the given PUT (see Fig. 1).

In this work we propose to use a binary partition tree (BPT)-based approach to define the secondary data. The BPT is a hierarchical tree representation of the data that exploits the spatial and spectral information contained on the image [START_REF] Salembier | Binary partition tree as an efficient representation for image processing, segmentation, and information retrieval[END_REF]. By means of the BPT representation of the image, we propose an alternative definition of the guard and outer windows that serve to obtain the secondary data of a given PUT. This BPT-based definition of secondary data allows to overcome the limitations of the conventional sliding window approach. On one hand it helps to define more homogeneous background regions and on the other hand it is flexible in order to fit adequately the geometry of the background.

We assessed the performance of the proposed BPT-based approach compared to the conventional sliding windowsbased approach on a real hyperspectral dataset provided by the Target Detection Blind Test project 1 .

The remainder of the paper is as follows: in Sec. 2 and Sec. 3 we overwiew the hyperspectral RX AD algorithm and BPT representation respectively. In Sec. 4 we introduce the proposed BPT-based approach. In Sec. 5 we provide the experimental validation and finally, in Sec. 6 we give some concluding remarks.

RX ADAPTIVE ANOMALY DETECTION

The RX algorithm was derived from the Generalized Likelihood Ratio Test (GLRT) assuming Gaussian hypothesis [START_REF] Irving | Adaptive multipleband cfar detection of an optical pattern with unknown spectral distribution[END_REF]:

H 0 : y = b H 1 : y = s + b , (1)
where b represents the background and s denotes the presence of an anomalous signal. The adaptive detector is obtained by replacing the unknown parameters by their estimates. For example, an estimate may be obtained from the range cells surrounding the cell under test. The size of the cell has to be chosen large enough to ensure the invertibility of the covariance matrix and small enough to justify both spectral homogeneity (stationarity) and spatial homogeneity. The use of a sliding mask (adaptive RX) provides a more realistic scenario than when estimating the parameters using all the pixels in the image. Thus, the mean vector µ and the background covariance matrix Σ, are estimated from N signal free secondary data, y i = 1, . . . , N . The resulting GLRT decision rule is the following:

t RX (y) = (y -μSMV ) T Σ-1 SCM (y -μSMV ) H1 ≷ H0 λ, (2)
where μSMV and ΣSCM denote the sample mean vector (SMV) and sample covariance matrix (SCM) estimators respectively:

μSMV = 1 N N i=1 y i , (3)
1 http://dirsapps.cis.rit.edu/blindtest Fig. 2. Construction of the Binary Partition Tree (BPT).

ΣSCM = 1 N N i=1 y i y T i . ( 4
)
For a Gaussian distribution, the quadratic form (yµ)

T Σ -1 (yµ) follows a χ 2 distribution when µ and Σ are known. This quadratic form is usually known as the Mahalanobis distance [START_REF] Mahalanobis | On the generalised distance in statistics[END_REF]. When the parameters µ and Σ are replaced by their maximum likelihood estimate parameters, μSMV and ΣSCM , the distribution of the quadratic form:

N -m + 1 N m (y -μSMV ) T Σ-1 SCM (y -μSMV ) ∼ F m,N -m+1 , (5)
can be approximated by the non-central F -distribution with m and Nm + 1 degrees of freedom [START_REF] Weisstein | CRC Concise Encyclopedia of Mathematics[END_REF]. For high values of N (i.e. N > 10m), the distribution can be approximated by the χ 2 distribution.

BINARY PARTITION TREES

In the BPT representation, the leaf nodes correspond to an initial partition of the image, which can be the individual pixels, or a coarser segmentation map [START_REF] Salembier | Binary partition tree as an efficient representation for image processing, segmentation, and information retrieval[END_REF]. From this initial partition, an iterative bottom-up region merging algorithm is applied until only one region remains. This last region represents the whole image and corresponds to the root node. All the nodes between the leaves and the root result from the merging of two adjacent children regions. An example of BPT is displayed in Fig. 2. If the initial partition contains n leaf nodes, the final BPT contains 2n -1 nodes.

Two notions are of prime importance when defining a BPT representation:

1. The region model M R which specifies how a region R is modelled.

2. The merging criterion O(M Rα , M R β ), which is a distance measure between the region models of any two regions R α and R β .

Each merging iteration involves the search of the two adjacent regions which achieve the lowest pair-wise dissimilarity among all the pairs of adjacent regions in the current segmentation map. Those two regions are consequently merged. Given a hyperspectral region R, with N R hyperspectral samples y j ∈ R q , j ∈ 1 . . . N R , the first-order parametric model M R is defined by the sample mean vector of the hyperspectral samples μR [START_REF] Valero | Comparison of merging orders and pruning strategies for binary partition tree in hyperspectral data[END_REF][START_REF] Valero | Hyperspectral image representation and processing with binary partition trees[END_REF]:

M R : μR = 1 N R N R j=1 y j . (6)
Using the first-order parametric model ( 6), a merging criterion is defined as the spectral angle distance, d SAM , between the sample mean vectors of any two adjacent regions [START_REF] Veganzones | Hyperspectral image segmentation using a new spectral mixture-based binary partition tree representation[END_REF]:

O M Rα , M R β : d SAM μRα , μR β , (7)
where

d SAM (a, b) = arccos ab a b
is the angular distance between any two vectors a and b.

The building of a BPT may suffer from small and meaningless regions resulting in a spatially unbalanced tree. To overcome this limitation, a priority term is included in the merging criterion that forces those regions smaller than a given percentage of the average region size to be merged first [START_REF] Tochon | Binary partition tree as a hyperspectral segmentation tool for tropical rainforests[END_REF].

BPT-BASED ADAPTIVE RX

Given a BPT representation of a hyperspectral image, let a branch of a PUT, x, denoted as B (x), define the sequence of nodes ascending on the BPT representation from the leaf containing the PUT up to the root node. For instance, given the BPT representation depicted in Fig. 2, the branch of a PUT in the leaf number 5 is: B (x) = {5 → 7 → 9}. This branch definition is a sorted list of nodes starting in the leaf node and ending in the root node. Then, the guard and outer nodes as defined as follows:

• The guard node of a PUT, G (x), is the first node in the PUT's branch, B (x), containing the PUT and at least, a given number of guard pixels.

• The outer node of a PUT, O (x), is the first node in the PUT's branch, B (x), containing the PUT, the guard pixels and at least a given number of secondary pixels.

Thus, in order to define the guard and outer nodes it is necessary to set the number of guard and secondary pixels, which work in a similar fashion to the guard and outer windows size respectively. The secondary data of a PUT, S (x), is then defined by the pixels contained in the outer node once the pixels in the guard node have been removed:

S (x) = O (x) \G (x) , (8)
where A\B denotes the complement operation between sets A and B.

EXPERIMENTAL METHODOLOGY AND RESULTS

Methodology

In order to test the proposed approach, we used a hyperspectral dataset provided by the Target Detection Blind Test project. The dataset includes a high-resolution hyperspectral image, spectral libraries of targets in the scene, and the location of targets in the scene. The hyperspectral imagery was collected by the HyMap sensor operated by HyVista 2 . The georegistered HyMap imagery has approximately 3 meter ground resolution and it is available both in calibrated spectral radiance as well as in spectral reflectance after atmospheric compensation [START_REF] Snyder | Development of a web-based application to evaluate target finding algorithms[END_REF]. Fig. 3 depicts a false color image of the dataset. Three civilian vehicles and four small (1m -3m) fabric panels were used as targets, for a total of 129 target pixels in the image. Fig. 4 shows the location of the targets in the map. The output of the RX adaptive detector, either using the sliding window or the BPT-based approaches, has been calculated using equation [START_REF] Acito | On the cfar property of the rx algorithm in the presence of signaldependent noise in hyperspectral images[END_REF]. For the sliding window approach, the guard window is of size 3 × 3 and the outer window is of size 21 × 21. For the BPT-based approach, the guard node should have at least 9 pixels while the outer node should have at least 400 pixels. For each of the two approaches, the probability of detection (P D ) has been calculated for different values of the probability of false alarm (P FA ), using as ground-truth the location of the 129 target pixels given by the dataset provider. hyperspectral RX algorithm over the test dataset for probabilities of false alarm set to 0.05 and 0.10 respectively, using the conventional windowing approach and the proposed BPT-based approach. It is possible to notice that the detected anomalies are more homogeneous and less noisy on the proposed BPT approach, what it suggests that this general approach to identify the secondary data could successfully be extended to other anomaly and target detection techniques.

Results

CONCLUSIONS

The proposed BPT-based approach to estimate the statistical parameters need for the RX adaptive hyperspectral AD show an improvement over the results obtained by the conventional sliding window-based approach on the Target Detection Blind Test's hyperspectral dataset. These preliminary results encourage the authors to extend this approach to other adaptive anomaly and target detector algorithms. Further experiments will introduce robustness of the proposed approach in terms of the CFAR property of the detection algorithms.

  This work was supported by ERC AdG-2013-320594 DECODA and by the Delegation Generale de lArmement (DGA) under contract PDE2012600079.

Fig. 1 .

 1 Fig. 1. Sliding window over the pixel under test (green circle): (red) 3 × 3 guard window and (blue) 15 × 15 outer window.

Fig. 3 .

 3 Fig. 3. Hyperspectral dataset provided by the Target Detection Blind Test project.

Fig. 4 .

 4 Fig. 4. Location of the targets on the image.

Fig. 5 .

 5 Fig. 5. P FA -P D plot comparing the sliding window-based and the BPT-based RX adaptive hyperspectral anomaly detectors over the Target Detection Blind Test project's hyperspectral dataset.

Fig. 5 Fig. 6 .

 56 Fig. 5 shows the P FA -P D plot (ROC curves) comparing the window-based and the BPT-based RX adaptive hyperspectral anomaly detectors over the dataset. The plot shows that the BPT-based approach outperforms the window-based approach in almost all the probability of false alarm range. This improvement is due to the selection of more homogeneous regions to perform the estimation stage. These preliminary results support the use of BPT-representations for an improved estimation of the statistical parameters need in RX adaptive hyperspectral anomaly detection.

Figs. 6 Fig. 7 .

 67 Figs. 6 and 7 show the detections obtained by the adaptive


Isabelle Bloch

Alain Bretto

come

Mathematical morphology on hypergraphs, application

to similarity and positive kernel.


J Frontera-Pons

M A Veganzones

S Velasco-Forero

F Pascal

J P Ovarlez

J Chanussot

ROBUST ANOMALY DETECTION IN HYPERSPECTRAL IMAGING

Keywords: hypespectral imaging, anomaly detection, elliptical distributions, M-estimators

Anomaly Detection methods are used when there is not enough information about the target to detect. These methods search for pixels in the image with spectral characteristics that differ from the background. The most widespread detection test, the RX-detector, is based on the Mahalanobis distance and on the background statistical characterization through the mean vector and the covariance matrix. Although non-Gaussian distributions have already been introduced for background modeling in Hyperspectral Imaging, the parameters estimation is still performed using the Maximum Likelihood Estimates for Gaussian distribution. This paper describes robust estimation procedures more suitable for non-Gaussian environment. Therefore, they can be used as plug-in estimators for the RX-detector leading to some great improvement in the detection process. This theoretical improvement has been evidenced over two real hyperspectral images.

INTRODUCTION

Target detection (TD) and anomaly detection (AD) of multidimensional signals have proved to be valuable techniques in a wide range of applications, including search-and-rescue, surveillance, rare mineral and land mines detection, etc. TD aims to discover the presence of a specific signal of interest (the target) among a set of signals. Statistical TD is based on the Neyman-Pearson (NP) criterion, which maximizes the probability of detection for a given probability of false alarm. AD is a special case of TD in which no a-priori target is provided. Hence, the goal of AD is to detect signals that are anomalous respect to the background. The Reed-Xiaoli (RX) AD algorithm [START_REF] Reed | Adaptive multiple-band cfar detection of an optical pattern with unknown spectral distribution[END_REF] is considered as the benchmark algorithm in multidimensional AD. However, the RX detector performance strongly relies on the statistical parameters estimation. Accordingly, when the background is non-homogeneous or the noise independence assumption is not fulfilled, the detector performance can be deteriorated. Here, we highlight a third drawback in the estimation problems: the presence of outliers in the secondary data used for the parameters estimation. In hyperspectral imaging, the actual distribution of the background pixels differs from the theoretically predicted under Gaussian hypothesis. In fact, as stated in [START_REF] Manolakis | Non gaussian models for hyperspectral algorithm design and assessment[END_REF], the empirical distribution usually has heavier tails compared to the Gaussian distribution, and these tails strongly influence the observed false-alarm rate of the detector. One of the most general and acknowledged model for background statistics characterization is the family of Elliptically-Contoured Distributions (ECD). They account for non-Gaussianity providing a long tailed alternative to multivariate normal model. They are proven to represent a more accurate characterization of HSI than models based on Gaussian assumption [START_REF] Manolakis | Non gaussian models for hyperspectral algorithm design and assessment[END_REF]. Although non-Gaussian distributions have already been assumed for background modeling, the parameters estimation is still performed using classical Gaussian based estimators; as in the case of covariance matrix, generally determined by the Sample Covariance Matrix (SCM) and the mean vector with the Sample Mean Vector (SMV). These classical estimators correspond to the Maximum Likelihood Estimators (MLE) for Gaussian assumption. However, they lead to suboptimal detection schemes when the noise is a non-Gaussian process. When working on ECD framework the model can be used to assess the robustness of statistical procedures and to derive alternative robust estimators of the parameters, the mean vector and the covariance matrix [START_REF] Gini | Covariance matrix estimation for CFAR detection in correlated heavy tailed clutter[END_REF][START_REF] Frontera-Pons | A class of robust estimates for detection in hyperspectral images using elliptical distributions background[END_REF]. These can then be used as plug-in estimators in place of the unknown mean vector or/and of the covariance matrix. This is a simple but often efficient method to obtain robust properties for signal processors derived under the Gaussian assumption.

ELLIPTICALLY CONTOURED DISTRIBUTIONS

Hyperspectral data have been proven not to be multivariate normal but long tailed distributed. In order to take into account these features, the class of elliptically-contoured distributions is considered to describe clutter statistical behavior. It provides a multivariate location-scatter family of distributions that primarily serves as heavy tailed alternative to the multivariate normal model. An m-dimensional random complex vector y = [y 1 , y 2 , . . . , y m ] T with mean µ and scatter matrix Σ has an elliptical distribution if its probability density function (PDF) has the form [START_REF] Kelker | Distribution theory of spherical distributions and a location-scale parameter generalization[END_REF]:

f y (y) = |Σ| -1 h m ((y -µ) H Σ -1 (y -µ)) (1)
where H denotes the conjugate transpose operator and h m (.) is any function such as (1) defines a PDF. If the second-order moment exists, then Σ reflects the structure of the covariance matrix of the elliptically distributed random vector y, i.e. the covariance matrix is equal to the scatter matrix up to a scalar constant. It serves to characterize the correlation structure existing within the spectral bands. It is worth pointing out that the ECD class includes a large number of distributions, notably the Gaussian distribution, multivariate t distribution, K-distribution or multivariate Cauchy. Thus, it allows for heterogeneity of the background power with the texture.

ROBUST PARAMETERS ESTIMATION

Along with their well-known properties and their simplicity of analysis, the SCM and the SMV are the most extended estimates since they are the MLEs for Gaussian case.

μSMV = 1 N N i=1 y i ΣSCM = 1 N N i=1 (y i -μ)(y i -μ) H
(2) where N denotes the number of secondary data. However, such widespread techniques are suboptimal when the noise is a non-Gaussian stochastic process. This article reviews some robust procedures particularly suited for estimating the covariance matrix and the mean vector of elliptical populations.

The Fixed Point estimators, according to the definition proposed by Tyler in [START_REF] Tyler | A distribution-free M -estimator of multivariate scatter[END_REF], satisfy the following equations:

μF P = N i=1 x i (x i -μF P ) H Σ-1 F P (x i -μF P ) 1/2 N i=1 1 (x i -μF P ) H Σ-1 F P (x i -μF P ) 1/2
(3)

ΣF P = m N N i=1 (x i -μF P ) (x i -μF P ) H ((x i -μF P ) H Σ-1 F P (x i -μF P )) (4)
The Fixed Point estimates have been widely investigated in statistics and signal processing literature. We refer to [START_REF] Pascal | Performance analysis of covariance matrix estimates in impulsive noise[END_REF] for a detailed performance analysis. It is worth pointing out that ΣSCM and ΣF P have the same asymptotic Gaussian distribution which differs on their second order moment by a factor m+1 m N , i.e. for N sufficiently large, ΣF P behaves as a Wishart matrix with m m+1 degrees of freedom.

RX ADAPTIVE ANOMALY DETECTION

The RX algorithm was derived from the Likelihood Ratio assuming Gaussian hypothesis:

H 0 : y = b H 1 : y = s + b , (5)
where s denotes the presence of an anomalous signal. The adaptive detector is obtained by replacing the unknown parameters by their estimates. For example, an estimate may be obtained from the range cells surrounding the cell under test.

The size of the cell has to be chosen large enough to ensure the invertibility of the covariance matrix and small enough to justify both spectral homogeneity (stationarity) and spatial homogeneity. The use of a sliding mask provides a more realistic scenario than when estimating the parameters using all the pixels in the image. Thus, the mean vector µ and the background covariance matrix, Σ are estimated from N signal free secondary data surrounding the pixel under test, y i , i = 1, . . . , N . The resulting GLRT decision rule is the following:

t RX (y) = (y -μSMV ) H Σ-1 SCM (y -μSMV ) H1 ≷ H0 λ. ( 6
)
and λ is a given threshold. When Gaussian assumption is valid, the quadratic form (yµ) H Σ -1 (yµ) follows a χ 2 distribution for Σ and µ perfectly known. This quadratic form is usually known as the Mahalanobis distance [START_REF] Chandra | On the generalized distance in statistics[END_REF]. When the parameters and under Gaussian assumptions Σ and µ are replaced by their MLE parameters and under Gaussian assumptions (2), the distribution of the quadratic form

N -m + 1 m N (z-μSMV ) H Σ-1 SCM (z-μSMV ) ∼ F m,N -m+1
follows a Hotelling T 2 distribution F m,N -m+1 which is the non-central F -distribution with m and Nm + 1 degrees of freedom [START_REF] Weisstein | CRC concise encyclopedia of mathematics[END_REF]. For high values of N, (N > 10 m), the distribution can be approximated by the χ 2 distribution. However, real hyperspectral scenes can not be described only with Gaussian distribution, as mentioned above. In this work we explore the use of Fixed Point estimators in the classical RX detector :

t RX-F P (y) = (y -μF P ) H Σ-1 F P (y -μF P ) H1 ≷ H0 λ. ( 7
)
It is important to highlight that the distribution of this detector is still an open question, as far as the authors are aware.

RESULTS

The experiments were conducted firstly on a real hyperspectral image where artificial targets with known spectral signature were introduced as anomalies in the background, see Fig. 1. The original data set consist in 50 × 50 pixels with 126 bands. Most of the theory on covariance matrix estimation have been recently extended to complex value signals [START_REF] Mahot | Asymptotic properties of robust complex covariance matrix estimates[END_REF].

Since hyperspectral data are real and positive, we proposed to use a Hilbert filter in order to render them complex. However, it is important to note that the real component after Hilbert transform is still the original signal. To avoid the well-known problem due to high dimensionality we have chosen sequencially eleven bands in the complex representation. In this approach, both covariance matrix and mean vector are estimated using a sliding window of size 9 × 9, having N = 80 secondary data. The results for this image are shown on the Fig. The results obtained with [START_REF] Frontera-Pons | A class of robust estimates for detection in hyperspectral images using elliptical distributions background[END_REF] show that the robust detector t RX-F P is capable of locate all the artificial targets and present a lower number of false alarms. This improvement is due to the fact that Fixed Point estimators treat the outliers and impulsive samples in order for them to have a smaller contribution to the background characterization process, while the SMV-SCM estimates suffer from the presence of strong reflectance pixels in the secondary data.

The algorithm has also been applied for galaxy detection on the MUSE data cube. The Multi Unit Spectroscopic Explorer (MUSE) project (see [START_REF]Official website of the muse project[END_REF]) aims to provide astronomers with a new generation of optical instrument, capable of simultaneously imaging the sky (in 2D) and measuring the optical spectra of the light received at a given position on the sky. MUSE was installed on the VLT telescope and operational in 2013, and its performances are expected to allow observation of far galaxies up to 100 times fainter than those presently detectable. MUSE will deliver a 3D data-cube made of a stack of images recorded at 3578 different wavelengths over the range 465-930 nm. Each monochromatic image repre- sents a field of view of 60 × 60 arcsec, recorded with a spatial sampling of 0.2 arcsec. Each record results in a data cube of size 1570 MB encoding 3578 images of 300 × 300 pixels, possibly containing thousands of objects (galaxies) existing over different subsets of wavelengths! An example of MUSE data cube image is displayed in Fig. 3, from the 3578 available bands, we have chosen one band of each 100 after Hilbert transformation. The results for anomaly detection are presented in Fig. 3 These examples illustrate the robust behavior of Fixed Point estimators in non-Gaussian environments or for close targets detection problems.

CONCLUSIONS

The family of elliptical distributions is considered for impulsive background characterization in hyperspectral imaging. In this context, robust estimation methods for mean vector and covariance matrix are used to overcome the non-Gaussianity of the background and the presence of outliers or strong scatters in the secondary data. Moreover, the use of the robust Fixed Point estimators for anomaly detection purposes has been discussed and compared to the classical SMV-SCM Gaussian estimators. The theoretical improvement provided by the robustness of the estimators is borne out through two real hyperspectral images.

2

 2

  .

Fig. 1 .

 1 Fig. 1. (a) Original background image with artificial anomalies, (b) Endmember used in the experiment.

  (a) Classical RX detector (b)RX detector built with the FP estimates

Fig. 2 .

 2 Fig. 2.

  (a) MUSE data cube (a) Classical RX detector (b) RX detector built the FP estimates

Fig. 3 .

 3 Fig. 3.