Alexandre Barachant
email: alexandre.barachant@gmail.com

Anton Andreev

M Congedo

The Riemannian Potato: an automatic and adaptive artifact detection method for online experiments using Riemannian geometry

Keywords: EEG, Artifact detection, Riemannian geometry, BCI

The Riemannian Potato: an automatic and adaptive artifact detection method for online experiments using

Riemannian geometry

Introduction

In this work we consider as artifacts any kind of EEG signal different enough as compared to the normal baseline signal. Based on this new definition, covariance matrices are used as descriptors of EEG signals and a Riemannian metric is employed to compare these covariance matrices with an average covariance matrix estimated on the signal baseline. This framework is not specific to a particular kind of artifacts and allows us to take into account the spatial properties of the artifacts. A practical implementation of this method will be described, and results of the online detection will be shown.

Methods

The goal of the detection algorithm is to determine if a portion of EEG signal

T N X × ℜ ∈
recorded during a time window of T samples over N electrodes contains artifacts. In order to achieve this, a trial X will be represented by its spatial covariance matrix

T XX T = C 1 1 -
and the criterion for the detection will be based on a Riemannian distance computation. The main idea is to estimate a reference covariance matrix C and reject every trial which is too far, in term of Riemannian distance, from this reference matrix. The Riemannian distance between C and C is defined by [Förstner and Moonen, 1999]:

( ) ( ) n N = n R = C C d λ 1 1 2 log , (1) with n λ the eigenvalues of 2 1 2 1 - - C C C
. The trial corresponding to C will be considered as an artifacts if R d is greater than a threshold th. Thus, the detection algorithm requires two parameters: C , the reference point in the Riemannian manifold and the threshold th for the detection. The estimation of those two parameters is the important part of the algorithm. The reference point could be estimated in an adaptive manner during the whole recording session according to the following equation:

[ ] 2 / 1 / 1 2 / 1 2 / 1 2 / 1 1 ) ( ) ( ) ( ) ( t t t t t C C C C C C α - - + = (2)
with t C the reference matrix from the previous iteration, C the current covariance matrix andα a coefficient which defines the speed of the adaptation. This adaptation is done only when clean signal is detected, i.e., the distance is lower than the threshold. The threshold th is estimated based on the mean µ and standard deviation σ of the distance to the reference matrix defined in Eq. 1 : σ µ

5 . 2 + = th (3)
These two parameters define a the Riemannian manifold. Sin metric is non-linear, this corresponds to a "potato" manifold. Fig. 1 shows the po 100 2x2 covariance matrices on Each point represents a covar manifold. The big black poin reference matrix and the grid r the potato where the Rieman reference point is equal to th.

Results

This algorithm was implem during a P300 experiment. EE electrodes. After a bandpass filt a step of 100ms). The parameter point is done at the beginning of In the example Fig. 2, the p artifacts) and eye movements (4

Discussion

The efficiency of this me considering the covariance matr spatial structure of the artifacts structure of the EEG channels. enough from the reference activ the initialization of the referen sensitivity and specificity strict other hand, because of the good standard deviations away from t method is based on two facts: first, we use multiva atrices as EEG signal descriptors. This allows us to tak acts. For this reason the artifact detection is sensitive ls. Second, by using a strategy where an artifact is ev ctivity, the algorithm is sensitive to many kinds of artif rence point is critical for the good functioning of this ictly depend upon its correct initialization and correct a ood sensitivity of the Riemannian metric, the artifacts u m the reference point, so the threshold estimation

Figure 2 :

 2 Figure 2: Results of the online det shown (Fp1,Fp2, AFz,Cz, Pz) .Grey

  is not c supported by ANR (Agence Nationale de la Recherché on Française contre les Myopathies) ROBIK. tric for covariance matrices. In Tech. Report of the Dpt of Geodesy o M, Maby E, Delannoy V, Bertrand O, Lécuyer A. OpenViBE: An rain-Computer Interfaces in Real and Virtual Environments. Presen 2010.

Figure 1 :

 1 Figure 1: Riemannian potato for covariance matrices on simulated

Acknowledgements

This research has been sup ROBIK and AFM (Association


Florent Lafarge

Xavier Descombes

Josiane Zerubia

Marc Pierrot-Deseilligny

Building reconstruction from a single DEM

We present a new approach for building reconstruction from a single Digital Elevation Model (DEM). It treats buildings as an assemblage of simple urban structures extracted from a library of 3D parametric blocks (like a LEGO R set). This method works on various data resolutions such as 0.7 m satellite and 0.1 m aerial DEMs and allows us to obtain 3D representations with various levels of detail. First, the 2D supports of the urban structures are extracted either interactively or automatically. Then, 3D blocks are placed on the 2D supports using a Gibbs model. A Bayesian decision finds the optimal configuration of 3D blocks using a RJMCMC sampler. Experimental results on complex buildings and dense urban areas are presented using data at various resolutions 1 .

Introduction

Three dimensional models of urban areas are very useful for many applications such as urban planning, radiowave reachability tests for wireless communications, disaster recovery or computer games.

Problem statement Many automatic methods have been proposed. It is difficult to compare these methods efficiently since they have been developed in different contexts (kinds of data, types of reconstructed buildings, level of user interactivity, etc) and use different evaluation criteria. Multiple view images are the most common inputs. Scholze et al. extracted 3D-lines and grouped them into faces which allow the building reconstruction through a semantic interpretation [12]. Rooftop hypotheses were generated from 3Dlines and junction information by Kim et al.[5]. Baillard et al. present a method based on planar facet hypothesis which can be generated from single 3D-lines [1]. Mueller et al. developed a procedural model based on a shape grammar [START_REF] Mueller | Procedural modeling of buildings[END_REF]. One of the most efficient methods uses a hierarchical approach which combines three different levels of detail [START_REF] Lee | Interactive 3D building modeling using a hierarchical representation[END_REF]. Laser scans are also popular inputs due to the measurement accuracy. Fruech et al. [3] uses laser scans to model buildings with a detailed reconstruction of the facades. These methods provide convincing 3D-models using aerial images, ground views or laser scans. However, most of them have been developed using a specific kind of data and cannot easily be adapted if the image characteristics change (resolution, SNR, etc). We propose a new method adapted to varying data resolution.

Global strategy There are two main families of approaches in 3D building reconstruction. Generic representations are theoretically able to reconstruct any shape of building through connected planar facets, but they demand high resolution data. Parametric representations are known to be robust with respect to data quality [START_REF] Weidner | Fowards Automatic Building Reconstruction from High Resolution Digital Elevation Models[END_REF], but these reconstructions are limited -most parametric representations consider a symmetric two-plan roof reconstruction. In this paper, we use an approach which is halfway between generic and parametric representations : the structural approach. It consists in reconstructing buildings by assembling simple urban structures extracted from a library of 3D parametric blocks, as a LEGO R game (see Figure 1). This approach is particularly interesting since it combines the advantages of generic and parametric representations: 1-the robustness of parametric approaches is preserved since the library objects are defined by parameter sets, 2-an efficient library permits modeling a large range of buildings. It is even possible to reconstruct buildings such as curved roof structures that some generic models cannot construct, 3-assembly of urban structures is particularly adapted to multiresolution data and allows us to obtain both basic roof shapes from 70 cm resolution satellite images or details including superstructures (chimneys, dormer-windows,etc) from 10 cm resolution aerial images. This approach is based on important prior knowledge concerning urban structures and their assembly. It is necessary to define the interactions between blocks correctly to have a convincing modeling without artefacts. A stochastic framework is especially well adapted to introduce such knowledge. Lafarge et al. proposed a similar approach [START_REF] Lafarge | An automatic building reconstruction method : A structural approach using high resolution images[END_REF]. However, that work was limited and suffered from many drawbacks: 1-generation of many artefacts, 2-lack of 3Dmodeling generality, 3-tuning of many parameters, 4-long computation times and 5-restriction to satellite data. The input was a single Digital Elevation Model (DEM) which is well adapted to global geometric descriptions. First, the 2D-supports of the urban structures were extracted either interactively or automatically using Ortner et al.'s method [START_REF] Ortner | Building outline extraction from Digital Elevation Models using marked point processes[END_REF]. Then, 3D blocks are positioned on the 2D supports using a Gibbs model in Section 3. A Bayesian decision finds the optimal configuration of 3D blocks using a RJMCMC sampler. Experimental results on complex buildings and dense urban areas are shown using data at various resolutions in Section 4.

Building extraction

The first step extracts the 2D-supports of the urban structures from a DEM. Each 2D-support is represented by an unspecified quadrilateral (or triangle) and is associated with a specific part of a building. The 2D-supports of a building correspond to a set of connected quadrilaterals (i.e. non-overlapping quadrilaterals with common edges). We propose two different ways for the user to extract the 2Dsupports: an interactive one and an automatic one.

Interactive extraction An operator controls interactive extraction: four clicks add a 2D-support (these clicks are the four points of the quadrilateral). This method allows accurate extraction as we can see on Figure 2, but takes considerable operator time.

Automatic extraction

The second method is automatic: Ortner et al. obtain 2D-supports using an object approach based on marked point processes [START_REF] Ortner | Building outline extraction from Digital Elevation Models using marked point processes[END_REF]. This method generates a set of rectangles which provide a rough description of the 2D-supports. Then, the rectangles are transformed into connected quadrilaterals by fusing neighboring rectan-gles. This process is automatic but the supports are not as accurate as those from interactive extraction -see Figure 2.

3D reconstruction

Once the 2D-supports have been extracted, the buildings are automatically reconstructed through a density formulation. The first step specifies the 3D objects.

Library of 3D-blocks

The content of the library is a key point: if it is too limited (e.g. Lafarge et al. [START_REF] Lafarge | An automatic building reconstruction method : A structural approach using high resolution images[END_REF] had only flat and gable roof forms), the method loses generality. The proposed library, denoted by M and presented in Figure 3, allows us to reconstruct a large range of buildings through a collection of blocks. Each block possesses both a roof form and a variant:

• The proposed roof forms (denoted by F and illustrated in Figure 3(top)) include monoplane (F 1x ), multiplane (F 2x ) and curved roofs (F 3x ). Each roof form has a specific set of parameters F (the number of parameters varies between 1 and 6).

• The variants (denoted by V and shown in Figure 3(bottom) for a gable roof type) are specific to a roof form. They correspond to types of structure ends (hipped or straight ends) or structure junctions ("-","L","T" or "+" junctions). The variants also specify the orientation of the roof with respect to the quadrilateral 2D-support. There are 1 to 3 parameters in the set of the variant, denoted by V .

To sum-up, each block of the library M is defined by a tuple (F, V) and an associated parameter set θ = (F, V ). Some blocks can topologically be degenerated in some situations2 : these cases are not allowed in the process in practice. Details concerning this library are available in [6].

Bayesian formulation

The notation for the Bayesian formulation is: • S, a set of sites and Λ = {Λ(s)/s ∈ S}, a given DEM where Λ(s) represents the intensity of the site, s.

• C, the quadrilateral configuration representing the building 2D-supports associated with Λ. N is the number of quadrilaterals (see Figure 2). • S i , the subset of S whose sites are inside the quadrilateral, i ∈ C. • D = {Λ(s)/s ∈ S i , i ∈ C}, the set of data. • x, an element of the state space, T , which corresponds to a configuration of 3D-parametric blocks knowing the 2D-supports C. x = (x i ) i∈C = (m i , θ i ) i∈C where each block, x i , is specified by both a model, m i of the library M and an associated set of parameters, θ i . In the following, x i = (m i , θ i ) and m i will be referred to as an object (or block) and a model, respectively. • S xi , the function from S i to R which associates the roof altitude of the object, x i , to each site of S i .

We consider the random variable, X, distributed in T , which follows an unnormalized density, h. h is the posterior density of a configuration, x, of objects, given D. In a Bayesian framework, this density can be obtained from:

h(x) = h(x/D) ∝ h p (x)L(D/x) (1)
We must be able to build both a prior density, h p (x), and a likelihood, L(D/x).

Likelihood

The likelihood represents the probability of observing the data, D, knowing the configuration, x. By considering the hypothesis of conditional independence, it can be expressed in terms the local likelihood of objects, L(D i /x i ):

L(D/x) = i∈C L(D i /x i ) ∝ i∈C exp -Γ α (i) (S xi , D i ) (2)
where Γ α (i) (., .) is the distance from R card(Si) × R card(Si) to R defined by :

Γ α (i) (S xi , D i ) = s∈Si |S xi (s) -Λ(s)| α 1 α
(3)

To sum-up, the likelihood corresponds to the Z-error of the L α norm between the DEM and the objects. In practice, α = 3 2 is a good compromise between sensitivity and robustness to DEM errors.

Prior

The prior density introduces interactions between neighboring objects. It allows us to both assemble objects in order to propose a realistic building and compensate for the lack of information contained in the DEM. A neighborhood relationship on C must be set up to define the interactions: two distinct quadrilaterals, i and j ∈ C, are said to be neighbors if they have a common edge. The neighborhood relationship, denoted by ⊲⊳ (i ⊲⊳ j represents the set of neighboring pairs in C). In previous work, too many interactions were set up [START_REF] Lafarge | An automatic building reconstruction method : A structural approach using high resolution images[END_REF]. The number must be minimal to preserve robustness and avoid problems in parameter setting. We propose a simple and efficient prior which is defined through a single interaction.

To do so, we define an assembly law which tests whether two objects can be assembled together. Two objects x i = (m i , θ i ) and x j = (m j , θ j ) are said "joinable" (denoted by x i ∼ a x j ), if:

1 -F i = F j and 2 -roof top orientations are compatible and 3 -the common edge of the quadrilateral 2D-supports, i and j, is not a roof height discontinuity.

The first condition checks that the two blocks have the same roof form. The second and third conditions test whether the roof tops of the two objects can be connected. The prior favors "joinable" objects, i.e. homogeneous structures, which means that heterogeneous structures can also be reconstructed. Moreover, in order to avoid artefacts, the parameters of two "joinable" objects are encouraged to have similar values. To do so, the unnormalized density, h p , is expressed through a Gibbs energy, U p (i.e. h p (x) = exp -U p (x)), defined by:

∀x ∈ T , U p (x) = β i⊲⊳j 1 {xi∼axj } g(x i , x j ) (4)
where 1 {.} is the characteristic function. The parameter, β ∈ R + , weights the importance of the prior density with respect to the likelihood. The density only depends on this parameter which is computed using the Maximum Likelihood estimator under regularization constraints detailed in [6]. The function, g, taking values in [-1, 0], measures the distance between the parameters of two "joinable" objects:

g(x i , x j ) = D(x i , x j ) D max -1 = k ω k | θ i,(k) -θ j,(k) | D max - 1
(5) θ i,(k) and θ j,(k) are the k th element of the parameter sets F i and F j of the objects x i and x j , respectively. D max = max xi,xj D(x i , x j ) is the maximum value of the distance. ω k are weights which are introduced in this distance in order to normalize the parameter values according to the metric system. These weights are computed from the X, Y and Z resolutions and the configuration of quadrilaterals C. Figure 4 illustrates this interaction. If the two blocks belong to different roof types (for example a mansard roof and a semi-elliptic roof on the top right) or if the two objects do not have compatible roof orientations (bottom right), they will not be "joinable" and the energy will be null. On the contrary, if the two objects are "joinable", the energy will be negative : these configurations are favored. The nearer the parameters of the two objects, the lower the energy. The left configuration is the best one.

Optimization

We now find the object configuration which maximizes the posterior density,h(.), i.e. the Maximum A Posteriori (MAP) estimator,x M AP . This is a non convex optimization problem in a high and variable dimension space,T , since the blocks in library, M, are defined by a different numbers of parameters.

RJMCMC sampler

The Reversible Jump Markov Chain Monte Carlo (RJM-CMC) algorithm [START_REF] Green | Reversible Jump Markov Chains Monte Carlo computation and Bayesian model determination[END_REF] is well adapted to our problem. Several papers have shown the efficiency of the RJMCMC sampler for the problem of multiple parametric object recognition. For example, Dick et al. use such a sampler to reconstruct architectural buildings from terrestrial images where the parametric object set includes structures such as columns, buttresses, entablatures or drainpipes [2]. The RJMCMC sampler simulates a discrete Markov Chain (X t ) t∈N on T having π as an invariant measure (specified by the posterior density, h(.)) which makes small jumps between spaces of variable dimension respecting the reversibility assumption of the chain. The jumps are proposed according to the three kernels specified in the following and detailed in [6]:

• Kernel Q 1 : uniform jumps This is the classic kernel proposing a new state according to uniform distributions. It is enough to ensure that the Markov chain can visit any configuration of the state space. However, using only this kernel requires long computation time [START_REF] Lafarge | An automatic building reconstruction method : A structural approach using high resolution images[END_REF]. Thus we propose two efficient additional kernels Q 2 and Q 3 .

• Kernel Q 2 : data-driven jumps This kernel cleverly explores the state space using a data-driven process [START_REF] Tu | Image Segmentation by Data-Driven Markov Chain Monte Carlo[END_REF]. To do so, the state, x, is proposed knowing the data, i.e. according to a probability, p(x|D). More precisely, it firstly estimates the gutter roof height H g and the roof top height H t of the object concerned by the jump and, secondly, chooses its height values according to the Gaussian distributions, N ( H g , σ) and N ( H t , σ), respectively (in practice, σ = 1 m).

• Kernel Q 3 : regularization jumps In our application, the visual aspect of the result is very important: we need a kernel which proposes well-regularized objects, i.e. objects which are perfectly aligned with their neighbors. The new object, x i , must be proposed knowing the neighboring objects {x j /j ⊲⊳ i}, i.e. according to p(x i |{x j /j ⊲⊳ i}). The model is uniformly selected according to the models of the neighboring objects. The parameter values are chosen according to Gaussian mixtures depending on the parameter values of the neighboring objects.

In summary, the RJMCMC sampler is: at iteration t, if X t = x: 1-Choose the kernel Q i (x, .) with probability q i . 2-According to Q i , propose a new state y. 3-Take x (t+1) = y with probability:

min π(dy) π(dx) Q i (y, dx) Q i (x, dy) , 1 (6)
and take x (t+1) = x otherwise.

Simulated annealing

Simulated annealing is used to ensure convergence: the density, h(.), is substituted with h(.)

1
T t , where T t is a sequence of temperatures which tends to zero as t tends to infinity. Simulated annealing theoretically ensures convergence to the global optimum for any initial configuration, x 0 , using a logarithmic temperature decrease. In practice, we used a geometric decrease which is faster and gives an approximate solution close to the optimal one. The initial and final temperatures are estimated through the variation of the energy following White [15]. The process has two stages. At the beginning, i.e. when the temperature is high (see Figure 5-(d) showing two simple examples of simulations), the process explores the density modes and favors configurations which have a high density. In this exploration stage, the data-driven kernel Q 2 is mainly used (q 1 = q 3 = 1 8 q 2 = 0.1)3 . At low temperature 4 , the configuration is close to the optimal solution and does not evolve very much : it involves a detailed adjustment of the 3D-block parameters. In this second stage, the regularization kernel Q 3 is mainly used (q 1 = q 2 = 1 8 q 3 = 0.1).

Experiments

The results were obtained from satellite DEMs (0.7 m resolution) and aerial DEMs (0.1 and 0.25 m resolution). The results show the reconstruction of complex buildings and dense urban areas whose level of detail depends on both the choice of the extraction process (automatic or interactive) and the kind of data. There is no comparison with other approaches for the reasons underlined in the introduction, except where the context is similar [START_REF] Lafarge | An automatic building reconstruction method : A structural approach using high resolution images[END_REF]. DEMs have been generated from 3-view images using a multi-resolution implementation of a Cox and Roy optimal flow matching image algorithm [11]. The 3D ground truths are raster images. Generic textures were applied to the objects of Figure 6 for visualization. Figures 6 and7 present various examples of reconstruction (showing different roof types, roof height discontinuities, closed structures or complex roof junctions). These results are convincing. The 3D-blocks are correctly assembled and few artefacts are generated which means the process adapts to buildings with complex roof junctions. Results shown in Figure 6 were obtained from satellite DEMs with automatic 2D extraction. Even if some details are omitted, the shapes of buildings compare well to the ground truth and the generalization level is satisfactory with respect to the context. The different roof types were correctly identified and the roof height discontinuities were accurately located. The ground errors for the automatic 2D extraction process were satisfactory. The over-detection rate (in term of surface) was 9.7%. This rate can be improved by adding a vegetation mask to prevent the detection of trees. The rate of missed detection was quite high (15%). However, it was mainly due to low flat buildings in inner courtyards (one floor height structures) that the 2D extraction process cannot detect since these buildings have low DEM discontinuities. Without taking into account these low flat buildings, this rate falls to 4.5%. The altimetric Root Mean Square Error (RMSE) in this context (satellite data / fullyautomatic process) is 2.3 m. This is better than the 3.2 m error obtained earlier [START_REF] Lafarge | An automatic building reconstruction method : A structural approach using high resolution images[END_REF] in the same context, but still remains high. It is mainly due to both a non optimal positioning of 2D-supports in the automatic extraction process (which engenders important local altimetric errors at some locations) and inaccuracies in the DEMs (which mainly correspond to matching problems of non Lambertian surfaces such as glass rooves). The second example underlines the limits of the automatic 2D extraction process: some footprints (especially curved footprints) cannot be modeled accurately by sets of quadrilaterals. The proposed kernels achieved ac- ceptable computation times. Less than one minute is necessary to obtain the buildings of Figure 6 using a 3Ghz processor (vs 5 minutes [START_REF] Lafarge | An automatic building reconstruction method : A structural approach using high resolution images[END_REF]). Figure 7 presents examples of buildings reconstructed by the interactive 2D extraction process. In the two first examples, the results are clearly better since the 2D-supports are accurately located. The reconstruction is more detailed and the altimetric RMSE is 1.1 m from the satellite data. Figure 8 shows results on two typical European down town areas. The computation time for the second example was 35 minutes (0.6km 2 -about 700 objects). The proposed method also allows modeling roof details such as chimneys or dormer-windows. In fact, the library of 3D blocks is general enough to reconstruct such roof superstructures. Figure 9 presents accurate results of building superstructure reconstruction from a 0.1 m resolution aerial DEM. The altimetric RMSE for these images is 0.6 m -a very good value.

Conclusion

This new method is an interesting alternative to generic and parametric approaches. It presents several important characteristics. First, it obtains very good results from a single DEM. Moreover, contrary to other methods, this approach works efficiently on various data resolutions: a global description of the buildings from 0.7 m resolution satellite data or a detailled building reconstruction including roof superstructures from 0.1 m aerial data. The user can also choose the level of automation of the process since the 2D-supports can be extracted either interactively or automatically. To conclude, it is an adaptive method since other 3D-block types can be added to the library depending on the context. In future work, it would be interesting to improve the optimization step to achieve both higher precision and shorter computing time. Adaptive cooling schedules could be used in the RJMCMC sampler or the Jump-Diffusion processes which are efficient for similar optimization problems. Moreover, we should evaluate the potential of this method on other kinds of cities such as typical North American urban areas.

Figure 1 .

 1 Figure 1. Principle of the structural approach.

Figure 2 .

 2 Figure 2. Extraction of 2D-supports: DEMs (left), interactive (center)and automatic extraction (right).

Figure 3 .

 3 Figure 3. Library of 3D-blocks -the roof forms (3D and profile views) (top) and the variants of a gable roof (bottom)

Figure 4 .

 4 Figure 4. Prior energy -local configurations of various energies.

Figure 5 .

 5 Figure 5. Two simple examples of the optimization process -evolution of the configuration as the temperature decreases (d), associated with the satellite images (a), the ground truths (b) and the final results (c).

Figure 6 .

 6 Figure 6. Reconstructed buildings with automatic 2D extraction from satellite data (4 th row), satellite images (1 st row), ground truths (2 nd row) and DEMs with 2D-supports (3 rd row).

Figure 7 .

 7 Figure 7. Reconstructed buildings with interactive 2D extraction from satellite data (3 first examples) and 0.25 m resolution aerial data (2 last examples) (2 nd row), DEMS with 2D-supports (1 st row).

Figure 8 .

 8 Figure 8. Urban area reconstruction without and with textures (3 rd /2 nd columns) from satellite data (1 st row) and aerial data (2 nd row), DEMs (1 st column).

The first author thanks the French Mapping Agency (IGN) and the French Space Agency (CNES) for partial financial support during his PhD. The authors thank the CNES for providing satellite images.

For example, a semi-elliptic roof on a triangular support

The q i correspond to the probability of choosing the kernel Q i in an iteration of the RJMCMC sampler.

In practice, the second stage is detected when the accepted proposition rate computed on 1000 iterations becomes lower than 0.05.


Florent Lafarge

Xavier Descombes

Josiane Zerubia

FOREST FIRE DETECTION BASED ON GAUSSIAN FIELD ANALYSIS

The Thermical InfraRed (TIR) channel contains wave lengths sensitive to the emission of heat. The forest fires can be characterized by intensity peaks on TIR images. We present a fully automatic method of forest fire detection from TIR satellite images based on the random field theory. First, preprocessing is used to model the image as a realization of a Gaussian field which presents interesting properties. The fire areas which are supposed to be a minority are considered as anomalies of that field. We present a statistical analysis of Gaussian field to determine a degree of belonging of a cluster to the image background (i.e. a realization of a Gaussian field). We then extend this application to the estimation of the fire propagation direction.

INTRODUCTION

The forest fires cause a lot of damage and participate to the deterioration of the Earth ecosystem, especially to the global warming. One of the main issues is to quickly locate them in an automatic way. With the recent progress in the spatial domain, this detection problem can nowadays be tackled through satellite images. In this domain, various methods have been proposed. Most of them are based on radiometric analysis from Thermical InfraRed (TIR) images as the work of Flannigan and Vonder Haar [START_REF] Flannigan | Forest fire monitoring using NOAA satellite AVHRR[END_REF]. Such algorithms allow to estimate the ground temperature and efficiently detect fires using thresholding techniques. Hybrid approaches have also been developed such as [START_REF] Li | A hybrid contextual approach to wildland fire detection using multispectral imagery[END_REF] where TIR data are coupled with varied kinds of information. Other works such as [START_REF] Lafarge | Textural kernel for SVM classification in remote sensing : Application to forest fire detection and urban area extraction[END_REF] try to detect the forest fires through the emitted smokes by using learning methods. However, most of these methods need either training sets or a priori knowledge concerning the used sensors or the observed areas. We present a new method which does not use this kind of information. The proposed approach allows to automatically detect the forest fires from TIR satellite images. TIR channel contains wave lengths sensitive to the emission of heat. It means that forest fires can be characterized by intensity peaks on this kind of images. However, all intensity peaks of a TIR image are not necessary forest fires. The main difficulty, as we can see on Figure 1, consists in making the distinction between peaks due to fire and other peaks of the image background which can correspond to industrial activities for example. An anomaly detection method is proposed to achieve this distinction. Such an approach is particularly efficient to extract image irregularities, as it can be seen in [START_REF] Carlotto | A cluster-based approach for detecting man-made objects and changes in imagery[END_REF]. The proposed method is especially adapted to the extraction of anomalies specified by intensity peaks in remote sens-ing data. The idea consists in modeling the image as a realization of a Gaussian field and determining a degree of belonging of the peaks to the field with respect to both radiometric and spatial characteristics. Such a principle has been used in functional MRI by Worsley, Friston and Poline. The proposed approach is based on their works, especially on [START_REF] Friston | Assessing the significance of focal activations using their spatial extent[END_REF][START_REF] Poline | combining spatial extent and peak intensity to test for activations in functional imaging[END_REF][START_REF] Worsley | A three-dimensional statistical analysis for CBF activation studies in human brain[END_REF]. First, definitions and interesting properties concerning the Gaussian fields are mentioned. Then, a preprocessing is set up in order to consider the image as a realization of a Gaussian field. A statistical analysis allowing to detect the forest fires is presented in a third part. We then propose an extension to the estimation of the fire propagation direction. Finally, results are shown. The proposed method is tested on TIR satellite images (BIRD satellite -300 meter resolution) representing Portugal in May 2003.

GAUSSIAN FIELDS

In this section, definitions and properties related to the Gaussian field theory are presented. More details are available in [START_REF] Adler | The Geometry of Random Fields[END_REF] and [START_REF] Lafarge | Détection de feux de forêt par analyse statistique de la radiométrie d'images satellitaires[END_REF].

Definition 1 Let (Ω, F , P) be a complete probability space and T ⊂ R N a topological space. A N-dimensional real valued random field is a measurable mapping X : Ω → R T . X is homogeneous if ∀τ,t 1 , ...,t k ∈ R N , the vector (X(t 1 ), ..., X(t k )) follows the same distribution that the vector (X(t 1 + τ), ..., X(t k + τ)).

A u (F, S) = {t ∈ S : F(t) ≥ u} (1)
In imagery, A u represents the set of clusters of an image F at the threshold u (see Figure 2). Adler [START_REF] Adler | The Geometry of Random Fields[END_REF] proposes an estima- tion of the number m u of clusters in the excursion set A u . For a 2-dimensional homogeneous Gaussian field (centered and reduced), the expectation of m u is given by:

E[m u ] = S(2π) 3 2 |Λ| 1 2 exp- u 2 2 (2)
where S represent the volume of the field (i.e the number of pixels in the image) and Λ, the covariance matrix of the field derivatives. This matrix allows to characterize the smoothness of the field. This equation is a key point in the proposed method. However, in order to be able to use it in the imagery domain, the image must verify some hypothesis. In particular, the image must be modeled as a realization of a homogeneous Gaussian field. To do so, we approximate a Gaussian field by considering its F d distributions as Gaussian. It represents an acceptable approximation since the covariance matrix Λ has a minor importance in the algorithm that we propose. In other words, the pixel distribution of the image must, by approximation, be Gaussian (i.e. the image histogram must be Gaussian).

PREPROCESSING

Preprocessing must be applied on the image in order to satisfy equation ( 2). The image must be modeled as a realization of an homogeneous Gaussian field, which imposes that at least the pixel distribution must be well approximated by a Gaussian law.

To do so, the first step consists in separating the various image modes (ground, clouds, sea,...) and in keeping the mode which have the highest intensity mean (i.e. the "ground" mode which is the closest to pixels representing forest fires). An iterated "2-means" algorithm, detailed in [START_REF] Lafarge | Détection de feux de forêt par analyse statistique de la radiométrie d'images satellitaires[END_REF], is used to extract the mode of interest. The obtained sub-image has a unimodal distribution which includes potential forest fires.

The forest fires are supposed to be rare events which means that their surface is negligible with respect to the surface associated with the considered mode. This hypothesis, which is necessary to obtain a correct separation of the modes, is acceptable since the complete TIR images represent vast areas of countries such as Portugal.

Then, the obtained class is regularized by using mathematical morphology. This regularization allows to delete the isolated elements of the class and "to fill in the holes". Finally, the last stage consists in normalizing the resulting distribution by a histogram specification. The obtained subimage has a pixel distribution which is Gaussian. So, This sub-image can be considered, by approximation, as a realization of an homogeneous Gaussian field. These three stages are illustrated in Figure 3 where both the histogram of the image and a piece of the image are shown for each preprocessing steps.

STATISTICAL ANALYSIS

The forest fires correspond to intensity peaks in the subimage obtained in the previous section. The idea consists in considering this sub-image as a background (which is represented by a realization of a Gaussian field) where some anomalies (i.e. intensity peaks which represent forest fires) can exist. In this context, results concerning the Gaussian field theory, especially the equation ( 2), can allow to estimate whether a cluster belongs to the realization of a Gaussian field (i.e. the background). To do so, both radiometric and spatial information related to the cluster are used to define probabilities of belonging.

Radiometric information

The belonging of a cluster to the background can be estimated in function of the maximal intensity of this cluster. Let us consider C

x 0 u , a cluster at the threshold u having a maximal intensity x 0 (see Figure 4 -(left)). In the following, a cluster established at the threshold u is called a u-cluster. The probability that C

x 0 u belongs to the realization of a Gaussian field (noted B as background) is given through the equation ( 2) by the expectation of the number of x 0 -clusters to the expectation of the number of u-clusters ratio:

P(C x 0 u ∈ B) = E[m x 0 ] E[m u ] = x o u exp u 2 -x 2 0 2 ( 3
)
This probability is noted P H . The behavior of P H in function of u is particularly interesting: it has an unique minima which allows to fix a reference threshold u 0 . In the following, this threshold allows to extract all the clusters which correspond to potential forest fires. More details about the computation of u 0 are available in [START_REF] Lafarge | Détection de feux de forêt par analyse statistique de la radiométrie d'images satellitaires[END_REF].

Spatial information

The spatial characteristics of a cluster and more precisely its surface can also allow to estimate whether the cluster belongs to the background. Let N u be the number of pixels with an intensity superior to u, m u be the number of u-clusters , and n u be the surface (i.e. the number of pixels) of a u-cluster. Nosko [START_REF] Nosko | Local structure of Gaussian random fields in the vicinity of high level shines[END_REF] has established that n u follows an exponential distribution. We aim at estimating the parameter of this exponential distribution, i.e. the inverse of the expectation of n u . These variables are linked by the following expression [START_REF] Friston | Assessing the significance of focal activations using their spatial extent[END_REF]:

E[n u ] = E[N u ] E[m u ] (4)
The pixel distribution is Gaussian, so the expectation of N u is given by:

E[N u ] = S ∞ u (2π) -1 2 e -x 2 2 dx = S Φ(-u) (5)
The equation ( 2) allows to estimate the expectation of m u . Then, the expectation of n u can be expressed by the following equation:

E[n u ] = E[N u ] E[m u ] = Φ(-u) (2π) 3 2 |Λ| 1 2 exp-u 2 2 (6)
Let us consider C n 0 u , a u-cluster having a surface n 0 (see Figure 4 -(left)). The probability that C n 0 u belongs to the background corresponds to the probability that n u is superior to n 0 :

P(C n 0 u ∈ B) = P(n u ≥ n 0 ) = e - n 0 E[nu] = e (2π) 3 2 |Λ| 1 2 n 0 exp-u 2 2 Φ(-u) (7)
This probability is noted P S . In practice, Λ is computed using the empirical estimators of the variances and co-variances of the field derivatives.

Algorithm

The probabilities P H and P S allow to estimate the belonging of a cluster to the background using two different characteristics of this cluster (i.e. its maximal intensity and its surface). We aim at defining a decision law taking into account both probabilities. This decision law consists in considering the cluster as foreigner to the background if:

min(P H , P S ) < α ( 8
)
where α is the limit probability, i.e. a confidence coefficient on the result which acts as a "P-value". In practice, α is close to 10 -2 . This decision law mainly relies on the probability P H (i.e. on the radiometric characteristic of the cluster). The probability P S is tested only if the test on P H is rejected. The decision law allows to locally analyze a cluster. We aim at defining a global algorithm dealing with the complete image. First, we extract all the clusters at the threshold u 0 (defined in section 4.1) which correspond to potential forest fires. Then, each detected cluster is extended to a family of sub-clusters with respect to the family of thresholds U = (u -k , ..., u 0 , ..., u k ). This extension allows to have a better precision concerning the location of forest fires. In practice, we take k = 2 (see Figure 4 -(right)). Finally, the decision law is applied to all these clusters. This algorithm is summed-up as follows:

Algorithm 1 .

1. Fix the limit probability α, 2. Extract the set of u 0 -clusters, 3. For each detected cluster, -Establish the family of clusters w.r.t. the threshold set U, -For each element of this family, compute P S and P H , 4. Keep the cluster accepted by the decision law (equation ( 8)).

ESTIMATION OF THE FIRE PROPAGATION DIRECTION

The detected areas can provide additional information by studying the topography of their intensity peaks. Figure 5 represents the typical topography of an intensity peak corresponding to a forest fire. This peak is composed of two parts. The first one represents the pixels of very high intensity (right side) which present high discontinuities with the neighboring pixels: it corresponds to the fire fronts. The second part (left side) represents lower intensity pixels which refers to the burned area. By considering the alignment of these two areas, the fire propagation direction can be estimated. It provides interesting information related to the evolution of the fires, which can be useful for firebrigades for instance. In order to estimate the fire propagation direction, we propose a very simple but efficient process which consists in comparing the intensity weighted center to the spatial center of the cluster. Let C be the spatial center of a cluster and M, its intensity weighted center (i.e. M = ∑s s×x s ∑s x s where the s and x s are respectively the pixels and the pixel intensities of the cluster). The fire propagation direction θ is estimated by the following expression:

θ = arg( -→ CM) ( 9
)
where arg is the argument of a 2D-vector. The value of the modulus of -→ CM can be seen as a confidence coefficient on the estimated orientation. Figure 7 shows results obtained on different clusters. We remark that neighboring clusters have similar orientations. It can be explained by the local influence of the wind in the propagation of the forest fires.

RESULTS

The proposed approach has been tested on TIR satellite images (BIRD satellite -300 meter resolution -1000 × 4000 pixels) representing Portugal acquired at various dates of May 2003. Figure 8 shows a result on a piece of image associated with the ground truth. On the ground truth (which has a lower quality), the white clusters (corresponding to forest fires) have been voluntarily enlarged in order to make a better distinction. 16 forest fires are present on the ground truth. We detect 15 of the 16 forest fires with a confidence coefficient α = 0.05. One false alarm is located. This result is convincing compared to results obtained by [START_REF] Flannigan | Forest fire monitoring using NOAA satellite AVHRR[END_REF]. However, the choice of α mainly depends on the application. For example, firemen will prefer a low confidence coefficient (i.e. detecting as many forest fires as possible even if many false alarms are located). The computing time is acceptable: two minutes are necessary to detect the forest fires using a 2Ghz processor on a 1000 × 4000 image (i.e. an area corresponding to the surface of Portugal).

CONCLUSION

The proposed method allows to automatically extract the forest fires from TIR satellite images. The results only depend on α which acts as a confidence coefficient. Both detection rate and false alarm rate provide convincing values. Interesting information, related to the evolution of the fires, can also be obtained through the estimation of fire propagation direction. However, the radiometric and spatial information, which allow to define the decision law, are set up independently. It could be interesting to introduce the probability P(C x 0 ,n 0 u ∈ B) (i.e. the probability that a u-cluster having both a maximal intensity x 0 and a surface n 0 belongs to the background) in the decision law in order to get a more robust method. In the future, we should work on this point.

Figure 1 :

 1 Figure 1: Height-map of a piece of TIR image

Definition 3

 3 Let F(t) : R N → R be a function. ∀u ∈ R and S ⊂ R N , the excursion set of F in S over the threshold u is defined by :

Figure 2 :

 2 Figure 2: A realization of a 2-dimensional random field (left) associated with 2 excursion sets at different thresholds (right)

Figure 3 :

 3 Figure 3: Preprocessing steps -(a): initial histogram of the complete image (top) associated with a piece of the image (bottom) -(b): extraction of the mode of interest -(c): regularization of the extracted mode (d): histogram specification by a Gaussian distribution

Figure 4 :

 4 Figure 4: Left: Scheme of a u-cluster having a maximal intensity x 0 and a surface n 0 (profile and top views of the intensity peak) -Right: Family of sub-clusters and sup-clusters related to the threshold family U

Figure 5 :

 5 Figure 5: Typical topography of an intensity peak corresponding to a forest fire

Figure 6 :

 6 Figure 6: (top): Detection rate (red) and false alarm rate (blue) vs the confidence coefficient α -(bottom): the corresponding ROC curve

Figure 7 :

 7 Figure 7: Various clusters detected as forest fires (red to yellow : intensity value) associated with the estimation of their propagation direction (black arrow) Generally speaking, Figure 6 presents the detection rate and false alarm rate (in term of number of clusters) in function of the confidence coefficient. This figure also shows the ROC (Receiver-Operator Characteristics) curve. The optimal performances are obtained for a confidence coefficient contained between 0.05 and 0.01. For instance, α = 0.01 allows to detect 85% of the forest fires with 15% of false alarms.

Figure 8 :

 8 Figure 8: (left): Extract of the result obtained on a TIR image (BIRD c DLR), white clusters represent the detected forest fires -(right): ground truth provided by Alcatel Alenia Space

  The collection of measures F t 1 ,...,t n defined byF t 1 ,...,t n (B) = P((X(t 1 ), ..., X(t n )) ∈ B) where B ∈ n ( is a Borelian of R) is called the family of the F d (finite dimensional) distributions.

	Definition 2 A random field is Gaussian if its F d distribu-
	tions are Gaussian vectors.

ACKNOWLEDGMENTS

The authors thank the German Aerospace Center (DLR) and Alcatel Alenia Space for respectively providing the data and the ground truth.


Wilfried Glastre
email: wilfried.glastre@ujf-grenoble.fr

Olivier Hugon

Olivier Jacquin

Hugues Guillet De Chatellus

Eric Lacot

Demonstration of a plenoptic microscope based on laser optical feedback imaging

A new kind of plenoptic imaging system based on Laser Optical Feedback Imaging (LOFI) is presented and is compared to another previously existing device based on microlens array. Improved photometric performances, resolution and depth of field are obtained at the price of a slow point by point scanning. Main properties of plenoptic microscopes such as numerical refocusing on any curved surface or aberrations compensation are both theoretically and experimentally demonstrated with a LOFI-based device.

Introduction

In conventional imaging setups, when impacting the detector, only the position of incoming photons is recorded while the information on the angle of arrival is lost. This results in a severe loss of information (stereoscopic information for instance). To handle that problem a first theoretical solution was introduced by the French Nobel Prize Gabriel Lippmann in 1908 under the name "photographie totale" [START_REF] Lippmann | Epreuves réversibles. Photographies intégrales[END_REF], later called "light field" or "plenoptic imaging". Contrary to classical imaging setups, plenoptic imaging is a technology where both spatial and angular information are recorded. This property leads to unusual features in imaging such as the possibility to numerically refocus a blurry photo after the picture has been taken. This results in an extended depth of field without reducing the aperture of the objective lens and thus there is no need for increasing the exposure time i.e. no movement blurring [START_REF] Ng | Digital light field photography[END_REF]. The opportunity to correct the objective's aberrations by digital post-processing has also been demonstrated [START_REF] Levoy | Recording and controlling the 4D light field in a microscope using microlens arrays[END_REF]. One of the first experimental device based on a microlens array was presented by Levoy et al. [START_REF] Ng | Digital light field photography[END_REF][START_REF] Levoy | Recording and controlling the 4D light field in a microscope using microlens arrays[END_REF] but suffers from several drawbacks like a low sensitivity and a trade-off between the resolution and the latitude of refocusing.

In the first part of this paper a new kind of plenoptic microscope based on Laser Optical Feedback Imaging (LOFI) is presented [START_REF] Glastre | Synthetic aperture laser optical feedback imaging using a translational scanning with galvanometric mirrors[END_REF] and compared to the previous conventional technology (using microlens array). In this setup, the latitude of refocusing no longer depends on the resolution. Next to that, this imaging setup also has the advantage of being shot noise limited (ultimate sensitivity) [START_REF] Lacot | Comparative study of autodyne and heterodyne laser interferometry for imaging[END_REF][START_REF] Glastre | Deep and optically resolved imaging through scattering media by space-reversed propagation[END_REF] and is thus suitable for optically resolved imaging through scattering media where a few ballistic photons are available. These benefits are at the price of a slow acquisition caused by point by point scanning. In the two last parts, numerical aberration compensation and refocusing on a curved surface are theoretically and experimentally demonstrated with LOFI.

Theoretical background: defocusing

LOFI microscope [START_REF] Glastre | Synthetic aperture laser optical feedback imaging using a translational scanning with galvanometric mirrors[END_REF] (see fig. 1) is an ultra-sensitive laser autodyne interferometer imaging technique combining the high accuracy of optical interferometry with the extreme sensitivity of class B lasers to optical feedback [START_REF] Lacot | Comparative study of autodyne and heterodyne laser interferometry for imaging[END_REF]. Photons are emitted, frequency shifted by a value F e , retro-reflected by the target and finally reinjected into the laser cavity. In this autodyne method, an optical beating between the reference wave and the signal wave (the light backreflected by the target) takes place inside the laser cavity. This beating at the frequency shift F e of the laser is detected by a photodiode and demodulated, leading to amplitude and phase information on the reinjected electric field. Images are obtained because of a point by point scanning by galvanometric mirrors [START_REF] Glastre | Synthetic aperture laser optical feedback imaging using a translational scanning with galvanometric mirrors[END_REF]. Because one sets the total round-trip frequency shift F e close to the relaxation frequency of the laser, a large amplification of the optical beating inside the cavity by the laser gain is achieved. More precisely, it was previously shown that the detection is shot noise limited even at low power (~mW) [START_REF] Lacot | Comparative study of autodyne and heterodyne laser interferometry for imaging[END_REF][START_REF] Glastre | Deep and optically resolved imaging through scattering media by space-reversed propagation[END_REF] which is well suited for the exploration of biological media (no damage and high optical attenuation of the tissues). In addition, because the laser plays the role of both the emitter and the detector it results as a self-aligned and easily transportable device.

                            L y x j r L y x y x L h R (1)
This equation contains information about the spatial position of photons retroreflected by the target (thus a raw blurred image is obtained). By taking the Fourier transform of this expression, one obtains ( and  are the spatial frequency coordinates along x and y directions):

) 2 ) ( exp( ) 2 exp( ) , , ( 2 2 2 2 2                       L j r L H R (2)
This equation gives an angular information about the retro-reflected photons. Indeed because of the translational scanning, each plane wave contained in the signal is associated with a unique Doppler frequency; each frequency coordinates (,) corresponds to a plane wave propagating with an angle of   / 2 in the x direction (  / 2 in the y direction) with the optical axis. Finally from Eq. 1 and 2, one can see that information on both the position and the direction of propagation of retroreflected photons is accessible (the wavefront contains a complete information).

On the contrary to classical imaging device (fig. 2(a This double spatial/angular information allows obtaining images with unusual properties such as the possibility to refocus a blurry image after raw acquisition. Indeed, in a defocused image, photons incoming from a punctual target impact different pixels on the detector instead of arriving at the same pixel. But because they arrive with different angles too (fig. 2), the image can be numerically refocused by summing these light rays incoming from the target.

In the case of the conventional plenoptic setup based on microlens array (fig. 2), as each pixel corresponds to one light rayi.e. position plus angle of a photon -, the refocusing is possible by simply summing the N appropriates pixels [START_REF] Ng | Digital light field photography[END_REF]. Ng showed that with this refocusing technique the equivalent depth of field is multiplied by a factor N: this is related to the angular resolution, the aperture of the objective is sampled by N elements. This device however suffers from a trade-off between the spatial and the angular resolution which depends on the choice of N. An increase of N increases the "plenoptic property" that is to say the angular resolution but reduces the final resolution (the size of the microlens is proportional to N).

In the case of LOFI setup (fig. 1), the defocus is visible in Eq. ( 2) in the second term and corresponds to the quadratic phase dependence of the plane waves in the signal. To numerically refocus raw images, this phase has to be cancelled by multiplying the signal in Fourier space (Eq. ( 2)) by the phase filter H filt (L,,): This filter corresponds to the free space retropropagation transfer function over a distance L/2 (factor 2 is due to round trip configuration of LOFI, see [START_REF] Glastre | Synthetic aperture laser optical feedback imaging using a translational scanning with galvanometric mirrors[END_REF]). After filtering and inverse Fourier transform in spatial domain, one has the following final synthetically refocused signal:

          2 ) ( exp ) , , ( 2 2       L j L H filt
                             2 2 2 1 2 exp ) , , ( ) , , ( ) , ( r
y x L H L H TF y x h filt R SA     (4)
In this expression TF -1 is the inverse Fourier transform. After numerical refocusing, the resolution is equal to r / 2 whatever is the initial defocus L. As a result the spatial/angular resolution trade-off (and thus the accessible depth of field/spatial resolution) no longer exists. These results are obtained at the price of a lower acquisition speed. Indeed instead of being limited by the CCD rate (around 1 ms for a full field image), one has a slow point by point scanning for LOFI with a 100 s integration time by pixel (~ 30 s for a 512*512 pixels image). Beside that there is one last major drawback with LOFI: the photometric balance which degrades with the defocus during raw recording [START_REF] Glastre | Synthetic aperture laser optical feedback imaging using a translational scanning with galvanometric mirrors[END_REF] (the coupling of the reinjected electric field with the laser cavity acts as a confocal pinhole). This last problem limits the accessible latitude of refocusing but is partially compensated by the ultimate shot noise sensitivity of LOFI [START_REF] Lacot | Comparative study of autodyne and heterodyne laser interferometry for imaging[END_REF][START_REF] Glastre | Deep and optically resolved imaging through scattering media by space-reversed propagation[END_REF] (this is an important advantage over the low sensitive CCD detector of conventional plenoptic setups).

Compensation of aberrations

Aberration compensation is an important concern for imaging through heterogeneous biological media [START_REF] Zawadzki | Ultrahigh-resolution optical coherence tomography with monochromatic and chromatic aberration correction[END_REF] or for obtaining highly resolved images with a cheap objective. To correct aberrations, the first solution consists in the introduction of adaptative optics resulting in an aberration-free laser spot in the target plane. Spatial light modulators or deformable mirrors [START_REF] Wright | Adaptive optics for enhanced signal in CARS microscopy[END_REF][START_REF] Facomprez | Accuracy of correction in modal sensorless adaptive optics[END_REF] can be used in this way. Another way to handle that problem is to use a plenoptic detector, the compensation is made after the recording of the raw image by a numerical postprocessing. This capability has already been demonstrated with conventional plenoptic setups with benefits both in the field of photography [START_REF] Ng | Digital light field photography[END_REF] as well as in microscopy [START_REF] Levoy | Recording and controlling the 4D light field in a microscope using microlens arrays[END_REF]. More precisely, defocus is just one particular second order (relatively to the aperture) aberration; if one totally controls light field rays, there is no reason that higher order aberrations could not been corrected too. In what follows we demonstrate aberration compensation in a LOFI plenoptic imaging setup. To simplify the situation, constant aberration in the image fieldi.e. the raw PSF does not depends on the field -is first considered, Eq. ( 2) is turned into:

) , ( ) , , ( ) , , ( '
      aber R R H L H L H  (5)
In this expression H aber (,) is the plane waves dephasing term responsible for the aberrations of the optics. In order to both refocus and correct aberrations, the filter function is turned into:

) , ( ) , , ( ) , , ( 1 '
        aber filt filt H L H L H (6)
It seems reasonable to only consider second and third order aberrations for correction; that is to say coma and astigmatism (in addition to defocus). The reason is that for higher order of aberrations the degradation of the PSF is mainly due to the extreme rays in the objective aperture (rays crossing objective lenses far from the center). In our case the laser beam has a Gaussian shape which implies that these aberrant extreme light rays have low amplitude and thus do not impact strongly the final PSF. Finally, with this approximation, the aberrant dephasing term of Eq. ( 5) and ( 6) is now given by the following simplified expression: [START_REF] Zawadzki | Ultrahigh-resolution optical coherence tomography with monochromatic and chromatic aberration correction[END_REF] In this expression,  and  are the polar coordinates in the spatial frequency space (the cartesian coordinates are (,)) and  0 is the phase origin related to the orientation of coma and astigmatism aberrations relatively to axes x and y. To finish, H astig , H coma , A and B are respectively the dephasing of plane waves due to astigmatism, to coma and the coefficients of these two aberrations.

        0 3 0 2 cos exp ) , , ( ) ( 2 cos exp ) , , ( ) , , ( ) , , ( ) , ( 
                    jB B H jA A H B H A H H coma astig coma astig aber
To illustrate the capabilities of LOFI technique and to compare it to the theory, spherical silica beads of 40 m diameter (fig. 3(a)) are imaged by the LOFI setup depicted on fig. 1. The advantage of this object is that a single bead acts as a punctual reflector (Dirac) located at its center. The raw PSF (Eq. ( 1)) is then directly accessible. Contrary to the paraxial lens used in [START_REF] Glastre | Synthetic aperture laser optical feedback imaging using a translational scanning with galvanometric mirrors[END_REF], an objective (lens L 4 of fig. 1) introducing large aberrations is now placed in the setup. This objective is a simple plano-convex spherical lens of focal length equal to 8 mm and is voluntarily tilted to induce important asymmetrical aberrations relatively to the optical axis (astigmatism and coma). Fig. 3(b) shows a raw image with a defocus L ~ 300 m resulting in a spot size of 100 m in the target plane. The numerical aperture of the beam is equal to 0.25 which implies an expected ultimate resolution of r / 2  1 m in the case of totally corrected aberrations (after numerical refocusing). From raw image of fig. 3(b), the defocus is first removed which leads to fig. 3(c) where only remain astigmatism, coma and higher order aberrations. Then astigmatism and coma are removed leading to fig. 3(d). This final image can be compared to the image (fig. 3(e)) obtained with an aberration-free objective of close numerical aperture and resolution (NA  0.15 and r / 2  1.5 m). To properly correct the defocus and compensate the aberrations it is necessary to evaluate precisely the defocus and the amount of astigmatism and coma (i.e. evaluating L, A and B, see Eq. ( 7)). Best coefficients are evaluated by minimizing the sum of the amplitude of each pixel. The pertinence of this criteria was firstly demonstrated by Dubois et al in [START_REF] Dubois | Focus plane detection criteria in digital holography microscopy by amplitude analysis[END_REF] for autofocus and here the possibility to extent it to aberration compensation is experimentally proved. This criteria is in fact related to the minimization of phase scrolling in the final PSF which explains that it works with all aberrations in general, including defocus. To finish, on fig. 3(f) |h SA (L,x,y)*h aber (L,x,y)| is represented with h aber the inverse Fourier transform of H aber (Eq. ( 6), A, B and L are obtained with [START_REF] Dubois | Focus plane detection criteria in digital holography microscopy by amplitude analysis[END_REF]). This illustrates the compensation of astigmatism and coma which were initially present in the objective. By comparing figs. 3(c) and 3(d), an important improvement in the image quality can be observed, which confirms the interest of our numerical aberration compensation technique. However by comparing figs. 3(d) and 3(e), it is also possible to see that these corrections are not totally perfect. Differences between these images can be essentially explained by higher order non-compensated aberrations.

Refocusing on a surface

Getting a sharp image of a 3D object everywhere in the field (i.e. increasing the depth of field while keeping both high resolution and large aperture) is not a new issue. The simplest solution is called focus stacking [11]; in this technique, several images are acquired, each corresponding to a different plane of focus. This stack of images is then computed in order to select in each zone of the field the sharpest image; this results in a globally sharp image everywhere in the field. The problem with that technique is the important amount of time and data resources required. Another solution is the use of non diffractive beam which automatically results in an extended depth of field. Such beams can be obtained with axicons [START_REF] Zhai | Extended depth of field through an axicon[END_REF], the accessible depth of focus is however limited by the lens size.

Here we propose to use plenoptic imaging, with the advantage that only one raw acquisition is needed. The refocus distance now simply depends on the image field. However, as it was stated before, the depth of field extension is limited up to N (the number of pixels behind each microlens) with a conventional plenoptic setup (fig. 2(b)). The use of the new LOFI plenoptic setup is thus interesting. In the previous case of a plane target perpendicular to the optical axis situated at a distance L of the laser waist (see fig. 1), the numerical refocusing is simply obtained by a convolution product of the raw signal s R (x,y) with the refocusing filter h filt (L,x,y) equal to the inverse Fourier transform of H filter (L,,) (see Eq. ( 2)-Eq.( 4)):

) , , ( * ) , ( ) , ( y x L h y x s y x s filt R SA  (8)
In this expression * is the convolution product. It is faster to work in Fourier space with Fast Fourier Transform (FFT, see Eq. ( 2) and ( 4)) and to use a fast multiplication (Eq. ( 3)) rather than making convolution product directly (Eq. ( 8)). Indeed in the first situation the computation time is only proportional to n log(n) (n the total number of pixel in the processed image) whereas it is proportional to n 2 for the second one.

When the ratio laser waist / target distance depends on the field, things are different and the convolution of Eq. ( 8) of raw defocused image s R (x,y) is now performed with a refocusing filter which changes with the field (L now depends on x and y). More quantitatively, Eq. ( 8

As stated before, the problem with a direct calculation of Eq. ( 9) is the excess computation time similar to the convolution product proportional to n 2 . For example for a 512*512 pixels image it is close to one hour with a 3 GHz dual core processor. In order to reduce computation time, another processing is used: from raw signal, images refocused in several planes spaced by the Rayleigh distance are calculated using quick calculation described in Eq. ( 2), ( 3) and [START_REF] Glastre | Synthetic aperture laser optical feedback imaging using a translational scanning with galvanometric mirrors[END_REF]. One just has to select in each region of the field the sharpest image exactly like in focus stacking, except that images of the stack are calculated from one raw image. Finally, the computation time is proportional to n log(n) and reduces from one hour down to 30 seconds.

In order to get an experimental demonstration, the object of fig. 3(a) is used again but now the film where beads are glued is curved in the y direction and simply tilted in the x direction. The objective lens is the same than in fig. 3 Once again, the optimal surface of refocusing shown in fig. 4(c) is calculated locally with the help of the Dubois's criteria [START_REF] Dubois | Focus plane detection criteria in digital holography microscopy by amplitude analysis[END_REF] (minimization of the sum of pixels amplitude). More precisely, the best refocusing distance is evaluated on four zones in the images and from these measurements, an optimal surface is fitted by a second order polynomial function L(x,y) of equation ax + by 2 + cy + d (in x direction, the film is simply tilted) leading to fig. 4(c). A restriction to the second order is in agreement with the smooth curvature of the target. Of course if the surface were unknown, the best fitting polynomial function could change but the approach would have been the same. However it is noteworthy that low quality bead images are obtained on the edge of fig. 4(d). This is due to both an incomplete scan of these beads (not all plane waves are present) and to a non perfectly planar laser scan which acts as a phase drift. These two effects results in a degradation [START_REF] Glastre | Limitations of synthetic aperture laser optical feedback imaging[END_REF] of the final image in the radial direction.

Conclusion

Plenoptic imaging setups are devices having the ability to record both position and angle of arrival of photons incoming on a detector. In this paper, a new kind of plenoptic microscope based on LOFI technique is demonstrated and its performances are compared with those of a previously existing device based on microlens array. It is shown that one benefits from a shot noise ultimate sensitivity limitation at low laser output power of 10 mW. This new plenoptic microscope is thus well suited for imaging through fragile scattering media (like biological media). Another advantage of LOFI over microlens array techniques is an extended depth of focus and an absence of trade-off spatial / angular resolution. However this is possible at the price of slow point by point galvanometric scanning and of degradation of photometric performances when the target is not in the objective's imaging focal plane. Unique possibilities related to plenoptic properties are also demonstrated in the paper like numerical compensation of constant aberration (astigmatism and coma of the objective) or refocusing on a curved surface after the acquisition. In the future, we plan to extend this work by implementing a correction of local aberrations beyond local defocus. Such an improvement will enable aberrations-free images through both scattering and heterogeneous media, which paves the way to biological sample imaging [START_REF] Zawadzki | Ultrahigh-resolution optical coherence tomography with monochromatic and chromatic aberration correction[END_REF].

Fig. 1 :

 1 Fig. 1: Experimental setup of the synthetic aperture LOFI-based imaging system. The laser is a 10 mW cw Nd:YVO4 microchip collimated by lens L1. A beam splitter sends 10% output power on a photodiode connected to a lock-in amplifier which gives access to the amplitude and phase of the signal. The frequency shifter is made of two acousto-optic modulators which diffract respectively in orders 1 and -1 and give a round trip total frequency shift of Fe = 3 MHz. x-y plane is scanned by galvanometric mirrors MX (scan in the x direction) and MY (scan in the y direction) conjugated by a telescope made by two identical lenses L3. f3 is the focal lengths of L3. αX and αY are the angular positions of MX and MY. L4 is the objective lens, f4 is its focal distance and r is the final waist of the laser after L4. When a punctual target located at a distance L from the image plane of the objective is scanned (fig. 1), one obtains the blurred Point Spread Function (PSF) [4]: 2

  )) where only the position of impacting photons is recorded (one pixel = one position) by a Charge Coupled Device (CCD), imaging setup showing such properties are called total, light field or plenoptic. A first configuration was presented by Ren Ng et al. in Stanford University[START_REF] Ng | Digital light field photography[END_REF]. In the case of Ng's setup (fig.2(b)) a microlens array is placed just before the CCD, each microlens covering N pixels. The impacting photons position are thus sampled by each microlens while the angle of arrival of the photon is sampled by the N pixels behind each microlens (see fig.2(b)).

Fig. 2 :

 2 Fig. 2: (a) Conventional and (b) Plenoptic imaging setup based on microlens array.

Fig. 3 :

 3 Fig. 3: Images of 40 m diameter silica beads. Observation (a) with classical bright field microscope; (b), (c), (e), (f) with LOFI-based microscope. Image (b) is a raw defocused image (L = 300 m resulting in a spot size of 100 m in the target plane) with highly aberrant lens (NA = 0.25). From (b), image (c) is obtained after numerical refocusing and (d) with both refocusing and numerical aberration compensation (expected resolution of r / 2  1 m). On (e) comparison with a numerically refocused image acquired with an aberration-free commercial objective (NA = 0.15, r / 2  1.4 m). (f) is the PSF used for filtering which reflects the aberrations which are compensated. Images (b)-(f) have a size of 256*256 pixels. Images from figs. 3(a), 3(b)-(d) and 3(e) have not been acquired on the same zone of the field: this explains why no correlations on the placement of beads can be observed. Color map of images (b) to (f) have no physical significance and thus are not displayed; false color are automatically scaled between the minimum and the maximum of each image's amplitude.

  (e) (no aberrations). Raw image is presented on fig. 4(a) with a defocus comprised between 250 and 1000 m depending on the field. After numerical refocusing on an intermediate plane corresponding to L = 750 m (see fig. 4(b)), one can observe that only a small number of beads are sharp. These sharp beads are distributed along an arc of a circle because of the curved-tilted support film. To conclude fig. 4(c) shows the best focusing surface and fig. 4(d) represents the final computed image. In this last image, all beads are sharp despite their initial different positions along optical axis. Because of a Rayleigh distance Z R ~ 10 m, the calculation of 50 stack images in intermediate plane is needed.

Fig. 4 :

 4 Fig. 4: Demonstration of numerical refocusing on a curved surface. Images size is 512*512 pixels. Target is composed of a curved flexible film with silica bead glued on it (identical to fig. 2(a). (a) Raw defocused image (250 m < L < 1000 m depending on the field). (b) Image after numerical refocusing over a distance 750 m (on a plane). (c) Calculated surface of refocusing and (d) image after refocusing on this curved plane.