Guillaume Doux
email: guillaume.doux@inria.fr

Patrick Albert
email: albertpa@fr.ibm.com

Gabriel Barbier
email: gbarbier@mia-software.com

Jordi Cabot
email: jordi.cabot@inria.fr

Marcos Didonet
email: marcos.ddf@inf.ufpr.br

Del Fabro

Scott Uk-

Jin Lee

An MDE-based approach for solving configuration problems: An application to the Eclipse platform

Keywords: Configuration, MDE, Eclipse, Plug-in, Cartography

Most of us have experienced configuration issues when installing new software applications. Finding the right configuration is often a challenging task since we need to deal with many dependencies between plugins, components, libraries, packages, etc; sometimes even regarding specific versions of the involved artefacts. Right now, most configuration engines are adhoc tools designed for specific configuration scenarios. This makes their reuse in different contexts very difficult. In this paper we report on our experience in following a MDE-based approach to solve configuration problems. In our approach, the configuration problem is represented as a model that abstracts all irrelevant technological details and facilitates the use of generic (constraint) solvers to find optimal solutions. This approach has been applied by an industrial partner to the management of plug-ins in the Eclipse framework, a big issue for all the technology providers that distribute Eclipsebased tools.

Introduction

Complex software systems are built by assembling components (components in a broad sense, i.e. COTS, libraries, plug-ins,…) coming from different repositories. This simplifies the development of the system but inevitably introduces an additional complexity dimension due to the need of managing these components. Each component can evolve independently and new releases can introduce/break dependencies with other components.

In particular, this is becoming a huge problem in the Eclipse community where new tools are built on top of several other plug-ins already available in the platform, many times requiring a specific version of the plug-ins. Therefore, releasing a new Eclipse tool implies a precise build definition for the tool that must be continuously evolved.

Therefore, technology providers commercializing Eclipse tools are looking for solutions that help them to automate and optimize the build definitions for their tools so that end-users do not need to suffer all these configuration problems. Right now, this very costly and time-consuming task requires a dedicated engineer in the provider company. This engineer needs to manually provide all the information regarding the tool dependencies, the plug-ins that can satisfy those dependencies and also the repositories where the plug-ins are available. Moreover, once everything is defined, the generated build definition needs to be empirically tested. Clearly, for non-trivial projects, this process does not scale.

In this paper, we propose to overcome this situation by means of using Model Driven Engineering and Constraint Programming techniques to automate the generation of build definitions. This work has been done in collaboration with two industrial partners: Mia-Software 1This paper is structured as follows. Section , a well-known technology provider in the Eclipse community that leads several Eclipse projects and IBM that has contributed its expertise in commercial constraint programming tools.

2 discusses the motivations of our solution in an industrial environment. Section 3 presents the overall approach used to manage our Eclipse plug-ins configuration problem. Section 4 is focused on the management of the configuration as a Constraint Satisfaction Problem (CSP) whereas Section 5 describes the decision tree approach for finding configurations, and Section 6 illustrates the tool used to visualize the configurations. Section 7 presents the implementation, a comparison between the resolution approach described and the lessons learnt. Section 8 focuses on the related works and Section 9 concludes this study.

Motivation: Industrial Challenge

This work has been motivated by the need of Mia-Software to configure and control build definitions for its tools and to be able to update such definitions during the tool lifecycle.

A second (and more complex) requirement is to be able to tailor build definitions to different scenarios, such as targeting the minimal subset of elements to run the application in a headless mode (using scripts on a server), or selecting only non-GPL component to allow integration in proprietary applications.

The results of this work are being integrated in the MoDisco2 and EMF Facet3

As an example, the MoDisco Eclipse project alone contains 94 plug-ins (without the 30 test plug-ins) and depends directly or indirectly on around 920 additional plugins. For the time being, the MoDisco project has a dedicated plug-in to configure the build definition. This plug-in contains more than ten types of files to do so (ant files, xml files, properties files, cspec files, cquery files, mspec files, rmap files, xsl files, sh Eclipse projects and in the custom developments the company builds internally for its clients.

files and txt files). To be able to maintain all of these artifacts and to reproduce in a server environment the behavior of a development environment, a dedicated engineer is now assigned to the task. To initially configure the build definition one full month of the engineer was required. Unfortunately, due the continuous evolution of the Eclipse platform this is not just a one-time effort. Every time there are new relevant plug-in versions or a new release of the Eclipse platform available, two or three additional days are spent in adapting the configuration (e.g. to update the locations of update sites for dependencies). Furthermore, the correctness of the process cannot be detected until the application is rebuilt again.

Therefore, MoDisco is clearly a tool that could benefit from the results of our work. It is worth to note that this kind of complex dependencies scenario is not the exception but the norm and thus, any tool that improves the current state of the art could have a real impact on current industrial practices of Mia-Software and similar technology providers.

Overall approach

This section gives an overview of our Model Driven approach for solving Eclipse plug-ins configuration problems. Adopting an MDE approach has several advantages. First, it provides a homogeneous representation for all the technologies involved in the solution. Secondly, it allows designers to deal with the problem at a higherabstraction level where some irrelevant details are omitted. And finally, the own Eclipse platform is moving more and more towards the adoption of model-based solutions (as the b3 Eclipse models used to define build generations) so an MDE approach perfectly fits in this scenario. Our approach follows a three step process (see Fig. 1).

Figure 1 -Overview of the Eclipse plug-ins build generation process

In the first step, a discovery phase allows the designer to express the requirements for the tool/component/plug-in she wants to build and the possible locations (i.e. repositories) where to find plug-ins that satisfy those requirements. The information concerning these plug-ins such as the dependencies they need, the name, the version or every other useful information is stored in a plug-in model conforming to the metamodel presented in Figure 2. This metamodel allows the representation of the different elements needed for the plug-ins configuration representation. The main entity of this metamodel is the Plugin class whereas the main relationship between plug-ins is represented by the PluginDependency class. This class allows linking the plug-ins according to their dependencies. A second set of entities and relations is expressed with the JavaPackage and the JavaPackageDependencies classes. A JavaPackageDependency element allows the representation of the relation between a Plugin element and the imported JavaPackage elements. At this stage, we just have the "raw data", i.e. we have the candidate plug-ins but not yet the selected configuration.

This is done in the second step: a possible combination of the candidate plug-ins (i.e. respecting all their dependencies) is created either manually (visualization option), interactively (decision tree option) or automatically (constraint programming option). When several configurations are possible, the final selection can be driven by additional search criteria like newest versions of the plug-ins (default option), license or cost. We propose these three different ways to obtain a configuration since each one offers a different trade-off as explained later on. As a result of this step we get a refined plug-in model with information from the selected configuration for the build generation.

Figure 2 -Metamodel to represent the plug-ins

The last step is the generation of the final configuration file from the refined plugin model. In our case the configuration will be expressed as a b3 model but it could easily be expressed as a Maven4 file or Ant 5 the development of a new generation of Eclipse technology to help building and assembling software. It proposes an approach using model driven engineering to represent the different artifacts that are relevant for building applications. More specifically, b3 proposes a metamodel to represent all of artifacts needed for the build, and execution support for these build models. Therefore in our approach the generation of the final configuration file is realized using model to model (or model to text in the case of Maven and ANT) transformations.

Configuration as a CSP

Constraint Programming [START_REF] Jackson | Alloy: a lightweight object modelling notation[END_REF] is a declarative problem solving paradigm where the programming process is limited to the definition of the set of requirements (constraints). A constraint solver is in charge of finding a solution that satisfies the requirements. Problems addressed by Constraint Programming are called constraint satisfaction problems (CSPs). A CSP is represented by the tuple CSP = <V,D,C> where V denotes the finite set of variables of the CSP, D the set of domains, one for each variable, and C the set of constraints over the variables. A solution to a CSP is an assignment of values to variables that satisfies all constraints, with each value within the domain of the corresponding variable. We can represent the problem of configuration of Eclipse plug-ins as a CSP. This solution is a practical instantiation of the approach called Model Search [START_REF] Kleiner | Model Search: Formalizing and Automating Constraint Solving in MDE Platforms[END_REF].

The problem can be stated as follows: given a set of partially-connected Eclipseplug-ins and a set of constraints that must be satisfied, find one (or the optimal, according to a given property) valid and executable build distribution. The constraints may be of different nature. For instance, defining version dependencies between the plug-ins, or specifying one desired vendor.

The constraints are written using OCL++. OCL++ is an adaptation of OCL (Object Constraint Language) [6]. OCL++ simplifies OCL for writing CSP problems. For instance, it enables writing multi-class invariants, which is a common construct in CSP problems. It also enables writing optimization functions.

More specifically, the re-expression of a configuration Eclipse problem in terms of a CSP is implemented as a chain of transformation operations over the initial plug-in model. Since Existing CSP solvers cannot directly read EMF models and OCL++ constraints as input, we need to translate the input artifacts into the CP language supported by the constraint solver of the ILOG OPL-CPLEX development bundle [START_REF]ILOG OPL-CPLEX development bundle[END_REF] engine, which is OPL (Optimization Programming Language). The OPL engine enables adding optimization functions, i.e., to find the best solution given an optimization criterion. The main steps of the process are:

1. Transformation of the Eclipse plug-in metamodel and the constraints into the OPL language.

2. Transformation of the input model into the OPL data format. The separation of the input model and metamodels into two transformations enables having independence between the problem specification (metamodels + constraints) and the input models with the initial data to start the CP process.

3. Execution of the CSP engine. This operation is called model search. During this phase the input model is extended with the solutions found by the engine. However, the result produced by the CP engine is expressed as sets of integers, String and floats (the OPL output format).

4. Transformation of the output into a model conforming to the Eclipse plug-ins metamodel.

To facilitate the execution of these different steps a predefined script in charge of chaining the transformations is provided.

This approach combines the benefits of CSP with the advantages of expressing the problem at the model level (e.g. writing the constraints in OCL++). Moreover, the transformation chain makes the CP solver transparent to the user who provides and receives models as input/output of the process. Clearly, it would be even better that the transformation from the models to the CSP included an additional intermediate step where the CSP is expressed as instance of a solver-independent CSP metamodel. This would facilitate the utilization of different CSP solvers.

Decision Tree

Decision Trees is a strategy used in the field of Software Product Line (SPL) to illustrate all possible product configurations in terms of decisions on variations. It enables interactive configurations where user selects an appropriate decision for each variation to configure a particular product. The main benefits of utilizing decision tree in configuration are the clear presentation of all possible configurations and the ability to customize the configuration by allowing each decision made on variation to be based on different criteria. As only valid configurations are proposed in the decision tree, the build engineer work becomes simpler and safer using this approach. On the other hand, the strategy main limitation, compared to the CSP one, appears when there are a large number of choice criteria involving an important number of choices for the engineer. In that case, it can become difficult to manage efficiently a big configuration.

In order to take advantage of these benefits, we adapt the concept of decision tree and Sequoia, a UML based SPL tool embedded in Papyrus, to our configuration challenge. The main difference with the CSP implementation is in the generation of several build configurations instead of only one. This characteristic involves some user interactions for the final configuration choice. The process of Eclipse plug-ins configuration with Sequoia consists of different steps as described below:

1. Construction of the initial bundle model in UML -The initial bundle model obtained from the discovery phase must be transformed into a UML model since Sequoia is a UML specific tool. As in the Eclipse plug-ins problem, several versions of the same plug-in can exist, and thus, a way to identify plug-ins conforming to the same unique plug-in definition is needed. The plug-in definition can be seen as a "formal" plug-in having several "instance" plug-ins, corresponding to the different available version of this plug-in.

In the transformation, formal plug-ins are defined as classes to type all the instance plugins of the model. Then, plug-in locations are defined as packages to group all plugins in the same location. Plug-in classes include useful metadata such as version, price and license. Each possible instantiation of the plug-in (i.e. different versions or vendors for the plug-in) instantiate these classes with the appropriate information to be considered during the configuration. Dependencies for a plug-in are represented as a dependency relationship from the depender plug-in to the dependee plug-in class. 2. Extraction of dependency constraints -Once the initial bundle model is constructed in UML, the dependency constraints are extracted following the specific profile defined in Sequoia. In addition, the extraction process can also accommodate the dependencies with constraints on criteria by allocating all the instances of the class that meet the constraints. For example, a dependency from the plug-in instance 'a' to a plug-in class 'B' with the constraint limiting the version of 'B' to be less than 3.0 will be converted into the set of dependencies from 'a' to 'b' with version 1.0 and 'b' with version 2.0. 3. Computation of dependency constraints -Once all dependencies are expressed as constraints, Sequoia uses the formal verification tool Alloy Analyzer [START_REF] Jackson | Alloy: a lightweight object modelling notation[END_REF] to produce all feasible configurations. The result of the calculation represents all possible configuration of the plug-ins computed based on their dependencies and represented in a textual format. 4. Decision tree creation -After the computation, the extracted dependency constraints are analyzed against the textual result of the calculation to construct a decision tree with decision nodes representing dependency constraints and its resolution edges representing all the configuration decisions that satisfy that dependency. In addition, values of various plug-in criteria are calculated and indicated for each resolution edge. Users can use these values to make more informed choices when interacting with the tree.

Transformation of decision tree into final bundle model -Finally, an Eclipse

plug-in configuration interactively generated from the decision tree is transformed into a final bundle model.

Visualization

The visualization mechanism allows quickly checking if the obtained configuration fits the user needs and, if several possible configurations have been produced, the user can choose the one he prefers from the visualizations. For simple configuration problems, the visualization of the plug-ins suffices to manually define the optimal configuration. Nevertheless, this kind of approach cannot be used to manage configurations involving an important number of elements, as the generated graph becomes too complex to be understandable. As an example, a visualization of an initial bundle model is shown Figure 3.

This visualization component relies on the cartography plug-in Portolan 6

To visualize plug-in data we just need to link the plug-in metamodel with the generic cartography metamodel provided with Portolan and, optionally, configure the view definitions that filter the input data and specify how this data will be visualized. The relationship between the plug-in and the cartography metamodels is done by defining the plug-in metaclasses as subclasses of the two main cartography metaclasses (entity and relationship). Once this is done, transforming data conforming to the plug-in metamodel to data conforming to the cartography metamodel is trivial (it is mainly a simple copy transformation).

. Integration with Portolan is easy since Portolan uses a model driven approach for the cartography analysis and visualization.

Implementation and preliminary results

The MDE-based approach presented here for solving Eclipse configuration problems has been implemented as an Eclipse set of plug-ins that provide the discovery of available plug-ins, the computation of possible dependencies, their visualization and the final build generation services. These functionalities are briefly presented in this section.

The discovery functionality is implemented as an Eclipse file creation wizard. This wizard proposes the creation of a model which is conforming to our metamodel dedicated to the Eclipse bundles representation (presented in Figure 2). To this aim, it allows selecting several plug-ins present in the workspace and then choosing the update sites that have to be considered when discovering candidate plug-ins for the bundle model.

The connection with the configuration engines (both the CSP and the decision tree versions) has been implemented as described in their respective sections. Also, as indicated, the visualization service is implemented using a model driven cartography tool called Portolan. A specific extension of Portolan has been designed for this study to be able to visualize configuration models.

The build generation functionality takes the feedback from the previous configuration plug-ins and generates a b3 model representing the selected configuration. This is mainly done by executing an ATL transformation between the internal model conforming to the plug-in metamodel (presented in Figure 2) and the final b3 model conforming to the b3 metamodel. This b3 model will be processed by the b3 engine to drive the build generation (by retrieving the needed plug-ins and launching the different steps of the application build). An excerpt of the b3 metamodel (Figure 4) presents the main elements used for the build generation. The BeeModel class represents the build model root; this class contains references to the BuildUnit and Repository elements used. In the model, a BuildUnit represents something to build with b3, in our case it will be a bundle (in the general case, it can also be a library or any other kind of component). The repositories reference of BuildUnit allows knowing which repositories can be used for the build unit's resolution. The BuildUnitRepository class allows the declaration of a build unit repository location in b3. A specific type of it is represented by the BeeModelRepository class; this repository declaration refers to the BeeModel to use for building the components contained in the repository.

Figure 4 -Representative Excerpt of the b3 Metamodel

Of course, other alternative implementations (e.g. Maven, ANT) of this service could easily be provided using the own Eclipse extension mechanisms.

After the initial set of experiments we have been able to validate that all three strategies (CSP, decision tress and purely visualization) can be used to solve the configuration problem. Each one has its own trade-offs and is best suited to address a specific kind of configuration problem. This is exactly the reason why we decided to keep the three of them in the framework without clearly favoring any of them.

The CSP-based approach is the best option when looking for a completely automatic solution. It is also recommended when looking for a single solution according to a specific criteria and when dealing with very complex problems (on which human interaction is not feasible).

Decision trees is an intermediate solution. It does most part of the job automatically (calculating all possible solutions) but still gives some flexibility to the designer to influence the final choice.

Visualization per se is only useful for simple solutions as an aid for the designer but it is a good complement to the other two as a visualization tool for the computergenerated solutions.

Besides this, the realization of this project has also shown the benefits of MDE when used as a tool for the unification of heterogeneous domain, such as the Eclipse plug-ins and the constraint programming domains. By expressing both domains (i.e. technical spaces) as models, we obtain a homogeneous representation that facilitates the transformation/communication between them.

Nevertheless, these first experiments have also pointed out some challenges that need to be addressed in the future. Reasoning tools usually suffer from scalability problems and our scenario is not an exception. Sometimes user interaction is required just because the search space is too big to get an answer from the solver in a reasonable time and the designer must help to reduce it by providing additional constraints to limit the search. Also, our approach suffers from the lack of standards in the constraint programming domain. Even if part of the transformation chain is generic, the last steps are solver-dependent and need to be reimplemented if the development team decides to use a different solver in the future.

We are now working in both aspects. For instance, regarding the second one we are adopting the idea of a CSP solver-independent metamodel that abstracts until the last step the specificities of the solver to use. The translation of the configuration information present in a solver independent model into solver specific models becomes easy to specify using model transformations and should permit to choose the most appropriated solver for the configuration resolution.

Related Work

An alternative solution for plug-in management in Eclipse is called p2 [START_REF] Berre | Dependency management for the eclipse ecosystem[END_REF]. This solution proposes to use the metadata of the plug-ins to create a set of constraints that are solved with the SAT solver SAT4j 7 7 Sat4j website:

. In this approach there is no explicit modeling http://www.sat4j.org/ of the problem so designers cannot define additional constraints about the desired characteristics of the solution (e.g. to get an optimal configuration). Moreover, since the translation is adhoc cannot be reused. Besides this, this approach only focuses on one of the alternative strategies we have explored.

Another interesting proposal is [START_REF] Cicchetti | Towards a model driven approach to upgrade complex software systems[END_REF]. It proposes to use a model driven approach to represent the configuration and available packages for FOSS distributions. These models are then used to predict the effect of changes on the installed package base (e.g. upgrades) on the system configuration. As our own work, this approach brings the advantages of working at a higher abstraction level when compared with approaches that rely on a direct manipulation of the available metadata. Nevertheless, our approach is able to deal with a more general problem since we are able to create the entire configuration and not only simulate/predict what would happen if something changes.

The topic of translating models into other formalisms for an automatic analysis has been explored in several previous approaches (e.g. [START_REF] Brucker | The HOL-OCL book[END_REF][START_REF] Cabot | UMLtoCSP: a tool for the formal verification of UML/OCL models using constraint programming[END_REF][START_REF] Straeten | Using description logic to maintain consistency between. UML models[END_REF][START_REF] Anastasakis | Uml2alloy: A challenging model transformation[END_REF]) but they mostly focus on specific kinds of UML diagrams. Some of these techniques could be adapted to our configuration models and integrated in our framework to provide additional analysis capabilities.

Conclusions

This paper reports a collaboration between industrial and research partners to solve configuration problems faced by technology providers using a combination of modeldriven engineering and constraint programming techniques.

We have focused on the specific configuration problems for tools developed on top of the Eclipse platform that need to manage and solve a lot of plug-in dependencies. This use case has been provided by Mia-software, a software editor specialized in the application Model Driven approaches for the software lifecycle industrialization with plenty of experience in the development of Eclipse projects.

As further work, we plan to generalize our framework to deal with configuration problems in other domains. The core of the approach can be easily reused but specific metamodels (e.g. Linux packages metamodel for Linux distributions configuration) need to be developed for each specific application domain.

Figure 3 -

 3 Figure 3 -Screenshot of the visualization tool

Mia-Software: http://www.mia-software.com/

The MoDisco project: http://www.eclipse.org/MoDisco/

The EMF Facet project: http://www.eclipse.org/modeling/emft/facet/

The Maven project, http://maven.apache.org/ script. The b3 Eclipse project focuses on

The ANT project, http://ant.apache.org/

Portolan website: http://code.google.com/a/eclipselabs.org/p/portolan/

Acknowledgement: This work has been partially funded by ANR IdM++ project.


A New Statistical Similarity Measure for Change Detection in Multitemporal SAR Images and its Extension to Multiscale Change Analysis Jordi INGLADA, Grégoire MERCIER Abstract-In this paper, we present a new similarity measure for automatic change detection in multitemporal SAR images. This measure is based on the evolution of the local statistics of the image between two dates. The local statistics are estimated by using a cumulant-based series expansion which approximates the probability density functions in the neighborhood of each pixel in the image. The degree of evolution of the local statistics is measured using the Kullback-Leibler divergence. An analytical expression for this detector is given, allowing a simple computation which depends on the 4 first statistical moments of the pixels inside the analysis window only.

The proposed change indicator is compared to the classical mean ratio detector and also to other model-based approaches. Tests on simulated and real data show that our detector outperforms all the others.

The fast computation of the proposed detector allows a multiscale approach in change detection for operational use. The so-called multiscale change profile (MCP) is introduced to yield change information on a wide range of scales and to better characterize the appropriate scale. Two simple yet useful examples of applications show that the MCP allows the design of change indicators which provide better results than a monoscale analysis.

Index Terms-Change detection; multitemporal SAR images; Kullback-Leibler divergence; Edgeworth series expansion; multiscale change profile (MCP).

I. INTRODUCTION

R EMOTE sensing imagery is a precious tool for rapid mapping applications. In this context, one of the main uses of remote sensing is the detection of changes occurring after a natural or anthropic disaster. Since they are abrupt and seldom predictable, these events cannot be well temporaly sampled -in the Shannon sense -by the polar orbit satellites which provide the medium, high and very high resolution imagery needed for an accurate analysis of the land cover. Therefore, rapid mapping is often produced by detecting the changes between an acquisition after the event and available archive data.

This change detection procedure is made difficult due to the time constraints imposed by the emergency context. Indeed, the first available acquisition after the event has to be used whatever its modality which is more likely to be a radar image due to weather and daylight constraints.

The kind of changes produced by the event of interest are often difficult to model. The same kind of event -a floodcan have different signatures depending on where it happenshigh density built-up areas, agricultural areas, etc. -and on the characteristics of the sensor. Also, the changes of interest are all mixed up with normal changes, which can be the majority if the time gap between the two acquisitions is too long.

All these issues present us with a very difficult problem: detecting abrupt unmodeled transitions in a temporal series with only two dates 1 .

From this position of the problem, one can make the straightforward deduction that pixel-wise comparison between the two images will not be robust enough.

In the case of radar acquisitions, the standard detector is based on the ratio of local means [START_REF] Rignot | Change Detection Techniques for ERS-1 SAR Data[END_REF]. This detector is robust to speckle noise, but it is limited to the comparison of first order statistics. The classical model for SAR intensity introduced by Ulaby et al. [START_REF] Ulaby | Textural information in SAR images[END_REF] assumes that the texture is a zero-mean multiplicative contribution. Therefore, changes taking place at the texture level which preserve the mean value will not be detected by the mean ratio detector. One can thus assume a miss-detection behavior of the detectors using only the mean pixel values. This remark invites a more accurate analysis of the local statistics of the images to be compared. Bujor et al. [START_REF] Bujor | Application of log-cumulants to the detection of spatiotemporal discontinuities in multitemporal SAR images[END_REF] did very interesting work by analyzing the interest of higher order statistics for change detection in SAR images. They concluded that the ratio of means was useful for step changes and that the second and third order log-cumulants were useful for progressive changes appearing in consecutive images in multi-temporal series. Since higher order statistics seem to be helpful, one may want to compare the local probability density functions (pdfs) of the neighborhood of the homologous pixels of the pair of images used for the change detection.

Of course, this assumes that the pdfs are known, and that there exists a robust way to compare them. The estimation of pdfs can be made with different approaches, but the straightforward histogram method should be avoided due to the need a high number of samples for the estimation. Indeed, small analysis window sizes are required to yield high resolution change maps. In this paper, we will present several approaches for this estimation by using only a small number of samples for the local statistics estimation, up to order 4.

Once the pdfs are estimated, their comparison can also be performed using different criteria. Information theory shows that a good measure is the Kullback-Leibler divergence, also called information gain. We will use a symmetrical version of this measure and show that it is superior to the classical detector when the pdfs are correctly estimated.

Therefore, these measures will be based on the comparison of local neighborhoods where an analysis window for the computation of the local estimation of probabilities is used. The problem which arises here is the one of choice of window size. Since we are facing unmodeled changes, we cannot choose the window size to fit the size of the expected changes. An inappropriate window size can produce miss-and overdetections:

• when using a small window for a correlation analysis, no detection will be performed in a homogeneous area, which was globally changed to another homogeneous area; • on the contrary, when using a larger window size, change areas have to be of larger size or strong in intensity (relatively to the measure) to be detected. In these cases, it will produce a coarse resolution change map. One way to overcome this problem is by applying a multiscale change detection analysis.

Scale is to be understood in its geographic meaning, which is the spatial extent of the study area. It does not refer to the cartographic meaning of scale (the larger the scale the more detailed the information) [START_REF] Lam | On the issues of scale, resolution, and fractal analysis in the mapping sciences[END_REF]. For an interesting discussion on scale issues in remote sensing see [START_REF] Marceau | Remote sensing contributions to the scale issue[END_REF].

Image processing techniques for multiscale analysis often use the cartographic meaning and apply low-pass filtering and possibly sub-sampling. For change detection analysis, this filtering and sub-sampling can be justified in the case where the images are not perfectly registered [START_REF] Carvalho | Digital change detection with the aid of multiresolution wavelet analysis[END_REF]. In other cases we think that it is better to use all the available information, that is, maximizing the number of available samples by using increasing window sizes. Nevertheless, pyramidal multiscale decompositions can also be useful in the case of phenomena characterizations (see [START_REF] Mercier | Multiscale oil slick segmentation with Markov chain model[END_REF] for example).

Therefore, the main point of the problem is how to choose the largest window size which robustly detects the changes but which is small enough to preserve the resolution of the final map without miss-detections.

We propose to use multiscale change profiles, which are defined as the change indicator for each pixel in the image as a function of the analyzing window size. The computation of the change detection for each window size can be very time consuming. We present here a method for the computation of these profiles which allows the change indicator at scale n to be computed from the value obtained at scale n -1 plus a correction term which takes into account the addition of new samples only. Analytical expressions are given for three different change indicators. This paper proposes three main contributions:

1) an Information Theory-based similarity measure which uses full local statistics; 2) the use of cumulant-based series expansions of similarity measures, which allow a robust and fast computation by using a small number of samples; 3) the concept of multiscale change profile and its fast implementation using recurrence evaluations. The paper is organized as follows: section II presents the problem formulation; section III introduces the measures used for the production of a change image; in section IV we introduce the concept of multiscale change profile and present the mathematical formulation allowing its optimized computation; sections V and VI present the results obtained on simulated and real data respectively, and section VII concludes the paper and proposes some directions for future work.

II. PROBLEM FORMULATION

Let us consider two co-registered SAR intensity images I X and I Y acquired at two different dates t X and t Y respectively. Our objective is to produce a map representing the changes occurring in the scene between t X and t Y . The final goal of a change detection analysis is to produce a binary map corresponding to the two classes: change and no change. The problem can be decomposed into two steps: the generation of a change image and the thresholding of the change image in order to produce the binary change map. Figure 1 The overall performance of the detection system will depend on both the quality of the change image and the quality of the thresholding. In this work, we choose to focus on the generation of an indicator of change for each pixel in the image. For interesting approaches in the field of unsupervised change image thresholding, the reader can refer to the works of Bruzzone and Fernández Prieto [START_REF] Bruzzone | Automatic Analysis of the Difference Image for Unsupervised Change Detection[END_REF], [START_REF]An adaptive parcel-based technique for unsupervised change detection[END_REF], Bruzzone and Serpico [START_REF] Bruzzone | An Iterative Technique for the Detection of Land-Cover Transitions in Multitemporal Remote-Sensing Images[END_REF] and Bazi et al. [START_REF] Bazi | An Unsupervised Approach Based on the Generalized Gaussian Model to Automatic Change Detection in Multitemporal SAR Images[END_REF]. The reader may note that some of these approaches need a statistical modelling of the detectors' response, which is not presented here.

The change indicator can also be useful by itself. Indeed, often the end user of a change map wants not only the binary information, given after thresholding, but also an indicator of the intensity of the change and eventually a confidence level. In order to evaluate the quality of a change image independently of the choice of the thresholding algorithm, the evolution of the detection probability, P det as a function of the false alarm probability, P fa , may be evaluated in the case where a set of constant thresholds are applied to the whole image. These are the so-called Receiver Operating Characteristics (ROC) and the plots of P det (P fa ) are called the ROC plots.

III. DISTANCE BETWEEN PROBABILITY DENSITIES

The main difficulty in the multitemporal analysis of SAR images is the presence of speckle noise. When moving away from interferometric configurations, the speckle is different from one image to the other and it can induce a high number of false alarms in the change detection procedure. Because of the multiplicative nature of speckle, the classical approach in SAR remote sensing involves using the ratio of the local means in the neighborhood of each pair of co-located pixels. The Mean Ratio Detector, MRD, is usually implemented as the following normalized quantity:

r MRD = 1 -min µ X µ Y , µ Y µ X , (1)
where µ X and µ Y stand for the local mean values of the images before and after the event of interest. The logarithm of eq. ( 1) may also be used. Nevertheless, this operation does not modify the performance of the detector in terms of ROC even if the contrast of the image of change indicator is modified. However, the logarithm is used since it modifies the initial pdf of the image of change indicator and then facilitates the development of Bayesian thresholding approaches [START_REF] Bazi | An Unsupervised Approach Based on the Generalized Gaussian Model to Automatic Change Detection in Multitemporal SAR Images[END_REF]. This detector assumes that a change in the scene will appear as a modification of the local mean value of the image. If the change preserves the mean value but modifies the local texture, it will not be detected.

The change detection algorithm proposed in this paper extends the MRD by analyzing the modification of the statistics of each pixel's neighborhood between the two acquisition dates. A pixel will be considered as having changed if its statistical distribution changes from one image to the other. In order to quantify this change, a measure, which maps the two estimated statistical distributions (one for each date at a co-located area) into a scalar change index is required. Several approaches could be taken into consideration: the mean square error between the two distributions, the norm of a vector of moments, etc. We have chosen to use a measure derived from Information Theory called the Kullback-Leibler divergence [START_REF] Kullback | On information and sufficiency[END_REF].

A. Kullback-Leibler divergence

Let P X and P Y be two probability laws of the random variables X and Y . The Kullback-Leibler divergence from Y to X, in the case where these two laws have the densities f X and f Y , is given by :

K(Y |X) = Z log f X (x) f Y (x) f X (x)dx. (2)
The measure log f X (x) f Y (x) can be thought of as the information on x for discrimination between hypothesis H X and H Y if hypothesis H X is associated with pdf f X (x), and H Y with

f Y (x). Therefore, the Kullback-Leibler divergence K(Y |X) can
be understood as the mean information for discrimination between H X and H Y per observation. This divergence appears to be an appropriate tool to detect changes when we consider that changes on the ground induce different shapes of the local pdf.

Since the Kullback-Leibler divergence can be understood as the entropy of P X relative to P Y , it is also called information gain. It can easily be proved that K(Y |X) 0; K(Y |X) vanishes only when the two laws are identical. K(Y |X) can be used as a measure of the divergence from P Y to P X . This measure is not symmetric as it stands: K(Y |X) = K(X|Y ), but a symmetric version may be defined by writing:

D(X,Y ) = D(Y, X) = K(Y |X) + K(X|Y ) (3)
that will be called the Kullback-Leibler distance (KLD).

In order to estimate the KLD, the pdfs of the two variables to be compared have to be known. As stated in the introduction, the processing of high resolution change maps requires analysis windows of small size, which makes impossible the use of local histogram estimations. In the following sections, we will introduce several approaches which allow the estimation of the pdfs by using a limited number of samples only. This requires some a priori information on the data which can be introduced by using models of local statistics.

B. Gaussian KLD

As seen in section III, the classical detector of eq. ( 1) uses first order statistics only. Yet, second order statistics are often used for SAR image processing. For instance, many speckle reduction filters [START_REF] Lee | Digital Image Enhancement and Noise Filtering by Use of Local Statistics[END_REF], [START_REF] Frost | A model for radar images and its application to adaptive digital filtering of multiplicative noise[END_REF], [START_REF] Kuan | Adaptive noise smoothing filter for images with signal-dependent noise[END_REF] are based on the contrast coefficient σ 2 X /µ 2 X , that is, the ratio between the variance and the square of the mean value. If the local statistics have to be compared up to the second order, the local random variables, X and Y may be assumed to be normally distributed (i.e. of Gaussian law). Then, the pdf of P X can be written as:

f X (x) = G(x;µ X , σ X ) = 1 2πσ 2 X e - (x -µ X ) 2 2σ 2 X . (4)
An analogous expression holds for f Y (x). Fig. 2(b) shows the Gaussian approximation of the probability distribution of a small region of interest (Fig. 2(a)) extracted from a SAR image.

If this Gaussian model is used in eq. ( 3), it yields the Gaussian Kullback-Leibler detector (GKLD):

r GKLD = σ 4 X + σ 4 Y + (µ X -µ Y ) 2 (σ 2 X + σ 2 Y ) 2σ 2 X σ 2 Y -1. (5)
It can be seen that even in the case of identical mean values, this detector is able to underline the shading of textures which is linked to the local variance evolution. Nevertheless, the reader should note that the Gaussian model should not be used since SAR intensity values are always positive. However, this example has been given as a simple case of a parametric model which takes into account second order statistics. Since some Gaussianity may be introduced into the data when resampling and filtering the images during the pre-processing step, the Gaussian model may nevertheless be justified.

C. KLD using the Pearson system

The drawback of the GKLD is that SAR intensity statistics are not normally distributed, and the use of a bad model can induce bad performance of the detector, whatever the accuracy of the parameter estimation. In the absence of texture, the radar intensity follows a Gamma distribution:

f X (x) = 1 Γ(L) L σ X L e - Lx σ X x L-1 . (6)
The Gamma distribution is characterized by the following parameters: L is the number of looks and σ X , the square-root of the intensity SAR image. Γ(•) is the Gamma function.

In the presence of texture, the local statistics can deviate from the Gamma distribution. For instance, if the texture is modeled by a Gamma distribution with a shape parameter ν, the resulting intensity distribution follows a K-law [START_REF] Ward | Compound representation of high resolution sea clutter[END_REF]:

f X (x) = 2 x Lνx µ X L+ν 1 Γ(L)Γ(ν) K ν-L 2 Lνx µ X 1/2 , (7)
where K(•) is the modified Bessel function of the second kind and µ X is the mean of X.

More generally, it is now accepted that the statistics of SAR images can be well modeled by the family of probability distributions known as the Pearson system [START_REF] Delignon | Statistical modelling of ocean SAR images[END_REF]. It is composed of eight types of distributions among which the Gaussian and the Gamma distributions may be found. The Pearson system is very easy to use since the type of distribution can be inferred from the following parameters:

β X;1 = µ 2 X;3 µ 3 X;2 and β X;2 = µ X;4 µ 2 X;2 ,
where µ X;i is the centered moment of order i of variable X.

That means that any distribution from the Pearson system can be assessed from a given set of samples by computing the first 4 statistical moments. Any distribution can therefore be represented by a point on the (β X;1 , β X;2 ) plane. For instance, the Gaussian distribution is located at (β X;1 , β X;2 ) = (0, 3), and the Gamma distributions lie on the β X;2 = 3 2 β X;1 + 3 line. Details about the theory of the Pearson system can be found in [START_REF] Johnson | Distributions in statistics: continuous univariate distributions[END_REF].

Figure 2(c) shows an example of distribution estimation. The Pearson approximation fits the data better than the Gaussian one (Fig. 2(b)). The example shown corresponds to a Beta distribution of the first type with parameters β 1 = 2.51 × 10 -6 , β 2 = 1.87.

The Pearson-based Kullback-Leibler Detector (PKLD) was originally introduced in [START_REF] Inglada | Change detection on SAR images by using a parametric estimation of the Kullback-Leibler divergence[END_REF]. It does not have a unique analytic expression, since 8 different types of distribution may be hold. Therefore, 64 different possibilities for the couples of pdf exist. Once the couple of pdfs is identified, the detection can be performed by numerical integration:

r PKLD (X,Y ) = Z log f X (x; β X;1 , β X;2 ) f Y (x; β Y ;1 , β Y ;2 ) f X (x; β X;1 , β X;2 ) + log f Y (x; β Y ;1 , β Y ;2 ) f X (x; β X;1 , β X;2 ) f Y (x; β Y ;1 , β Y ;2 ) dx. (8)
The correct way to proceed to use the Pearson system is to choose a pdf using the estimated moments and then estimate the parameters of the distribution by maximum likelihood. While this can improve the results of the pdf estimation, the effect is not noticeable in terms of the estimation of the change indicator. This approach was not used in the experiments in order to reduce the computation cost.

The reader should note that in the case of single-look high resolution SAR data (better than 10 m) other statistical models may be more appropriate, mainly on urban areas. Nicolas et al. have proposed a new model based on the log-statistics and a set of pdfs coming from the Fisher system of distributions [START_REF] Nicolas | Introduction aux statistiques de deuxième espèce: applications des log-moments et des log-cumulants à l'analyse des lois d'images radar[END_REF], [START_REF] Nicolas | Comparison of Beta1-Beta2 diagram and Logcumulant2-Logcumulant3 diagram for SAR image texture classification[END_REF]. It has been applied to high resolution SAR images on dense urban areas with promizing results [START_REF] Tison | Accuracy of Fisher distributions and log-cumulant estimation in a markovian segmentation of high resolution SAR images over dense urban areas[END_REF], [START_REF] Tison | A new statistical model of urban areas in high resolution SAR images for Markovian segmentation[END_REF].

D. Cumulant-based KL approximation

Instead of considering a parameterization of a given density, or set of densities, it may be of interest to describe the shape of the distribution. Such a description is based on quantitative terms that may approximate the pdf itself. The cumulants themselves do not provide such a pdf estimation directly but are necessary to describe its shape: for example, third order (κ 3 ) is linked to the symmetry (i.e. skewness), while the fourth (κ 4 ) is linked to the flatness (i.e. kurtosis). The density is then estimated through a series expansion. In fact, the cumulant generating function is used for such an estimation. By definition, the cumulant generating function K X (•) of a random variable X is defined by:

K X (ω) = ln M X (ω) = ∑ n κ X;n ω n n!
with M X (•) being the moment generating function defined by:

M X (ω) = Z e ωx f X (x)dx = Z 1 + ωx + ω 2 2 x 2 + • • • f X (x)dx.
For the case of the four first order cumulants, the following expressions hold [26, p.8]:

κ X;1 = µ X;1 κ X;2 = µ X;2 -µ 2 X;1 κ X;3 = µ X;3 -3µ X;2 µ X;1 + 2µ 3 X;1 κ X;4 = µ X;4 -4µ X;3 µ X;1 -3µ 2 X;2 + 12µ X;2 µ 2 X;1 -6µ 4 X;1 . (9)
Let us assume that the density to be approximated is not too far [START_REF] Stuart | Kendall's advanced theory of Statistics[END_REF] from a Gaussian pdf (denoted as G X to underline the fact that it has the same mean and variance as X), that is, with a shape similar to the Gaussian distribution. The difference between K X (•) and K G X (•), can be written in terms of the difference of the cumulants κ X;n -κ G X ;n . By inversion, the density may be expressed by a formal Taylor-like series:

f X (x) = G X (x) + c 1 dG X dx + c 2 d 2 G X dx 2 + • • • Since a Gaussian density is used, it yields f X (x) = ∞ ∑ r=0 c r H r (x)G X (x),
with H r (x) known as the Chebyshev-Hermite polynomial of order r [START_REF] Stuart | Kendall's advanced theory of Statistics[END_REF]. When choosing a Gaussian law so that its first and second cumulants agree with those of X, the number of terms of the series expansion is greatly reduced. This is the so-called Edgeworth series expansion. Its expression, when truncated to order 6, is the following:

f X (x) = 1 + κ X ′ ;3 6 H 3 (x) + κ X ′ ;4 24 H 4 (x) + κ X ′ ;5 120 H 5 (x) + κ X ′ ;6 + 10κ 2 X ′ ;3 720
H 6 (x) G X (x). (10)
It can be thought of as a model of the form X = X G + X ′ where X G is a random variable with Gaussian density with same mean and variance as X, and X ′ , a standardized version of X [START_REF] Mccullagh | Tensor Methods in Statistics[END_REF] with:

X ′ = (X -κ X;1 )κ -1/2 X;2 .
Fig. 2(d) shows an example of such an approximation of a histogram.

The Edgeworth series expansion of the two pdfs f X and f Y may be introduced into the Kullback-Leibler divergence (eq. ( 2)). It yields an approximation of the Kullback-Leibler divergence by Edgeworth series, truncated at a given order. In [START_REF] Lin | Edgeworth approximation of the Kullback-Leibler distance towards problems in image analysis[END_REF], such an approximation has been truncated to order 4 by using the equality

f X f Y = f X G X G X G Y G Y f Y , where G X (resp. G Y ) is
a Gaussian density of same mean and variance as f X (resp. f Y ). Then,

KL Edgeworth (X,Y ) = 1 12 κ 2 X ′ ;3 κ 2 X;2 + 1 2 log κ Y ;2 κ X;2 -1 + 1 κ Y ;2 κ X;1 -κ Y ;1 + κ 1/2 X;2 2 -κ Y ′ ;3 a 1 6 + κ Y ′ ;4 a 2 24 + κ 2 Y ′ ;3 a 3 72 - 1 2
κ 2 Y ′ ;3 36 c 6 -6 c 4 κ X;2 + 9 c 2 κ 2 Y ;2 -10 κ X ′ ;3 κ Y ′ ;3 (κ X;1 -κ Y ;1 ) (κ X;2 -κ Y ;2 ) κ 6 Y ;2 (11)
where

a 1 = c 3 -3 α κ Y ;2 a 2 = c 4 -6 c 2 κ Y ;2 + 3 κ 2 Y ;2 a 3 = c 6 -15 c 4 κ Y ;2 + 45 c 2 κ 2 Y ;2 - 15 κ 3 Y ;2 c 2 = α 2 + β 2 c 3 = α 3 + 3αβ 2 c 4 = α 4 + 6α 2 β 2 + 3β 4 c 6 = α 6 + 15α 4 β 2 + 45α 2 β 4 + 15β 6 α = κ X;1 -κ Y ;1 κ Y ;2 β = κ 1/2 X;2 κ Y ;2 .
Finally, the cumulant-based Kullback-Leibler detector (CKLD) between two observations X and Y is written as:

r CKLD = KL Edgeworth (X,Y ) + KL Edgeworth (Y, X). ( 12
)
The reader should note the fact that, like for the Pearson-based detector, despite the apparent complexity of the formulas, and thanks to eq. ( 9), only the moments up to order 4 have to be computed.

IV. MULTISCALE CHANGE PROFILE

Scale plays a strategic role in image analysis and more especially in change detection applications. In section I, it has been shown how an inappropriate scale of analysis can produce miss-or over-detections. In [START_REF] Bovolo | A Wavelet-Based Change-Detection Technique for Multitemporal SAR Images[END_REF], Bovolo and Bruzzone, stress the fact that the scale of analysis is a key parameter for better discrimination between change and no change areas. Such a point of view is implemented by a wavelet transform of the log-ratio estimated with a window of a user-defined size.

Instead of applying a multiscale analysis of the change image, the purpose here is to produce a set of change indicators estimated at various scales. We will call it multiscale change profile (MCP).

As stated in the introduction, the multiscale term refers here to the size of the analyzing window. The MCP will therefore involve computing the change indicator for a pixel by using neighborhoods of increasing sizes. The so-called profile corresponds to the sequence of change measures as a function of scale. We will restrict our formulation to the case of the CKLD. Given the fact that this detector needs the estimation of the statistical moments of the samples inside the analyzing window, we are interested in finding an approach which avoids the computation from scratch of the moments at every scale.

A. Optimized computation of the MCP

Let us consider the following problem: how to update the moments when an N + 1 th observation x N+1 is added to a set of N observations {x 1 , x 2 , . . . , x N } already processed. When considering raw moments of order r, the formulation comes easily as:

μr,[N+1] = N N + 1 μr,[N] + 1 N + 1 x r N+1 .
μr,[N] (resp. μr,[N+1] ) stands for the raw moment of order r estimated with N samples (resp. N + 1 samples). Since the analyzing window may contain textured areas, the mean value itself may be modified by the increase in the number of samples. Therefore, by using simple binomial properties, it can be shown that central moments may be characterized by:

µ 1,[N] = 1 N s 1,[N] ( 13
) µ r,[N] = 1 N r ∑ ℓ=0 r ℓ -µ 1,[N] r-ℓ s ℓ,[N] ,
where the notation s r,[N] = ∑ N i=1 x r i has been used. Hence, when considering a new sample x N+1 , each moment may be updated directly by using updates of s 1,[N+1] and then s r,[N+1] for increasing values of order r. The Edgeworth series is also updated by transforming moments to cumulants (by using eq. ( 9)) to be introduced in eq. ( 10) and then in eq. [START_REF]An adaptive parcel-based technique for unsupervised change detection[END_REF].

Fig. 3 shows an example of a pdf estimation on a homogeneous area (shown in fig. 2(a)) when the window increases from 9 × 9 to 17 × 17. In fact, the availability of updating the estimation of the distance between distributions from windows of any size without re-processing the overall data is the most interesting point for multiscale change detection purposes. This on-line multiscale moment estimation is the key for the operational use of the MCP concept.

For example, the computation of r CKLD with windows of size ranging from 5 × 5 pixels to 51 × 51 pixels (22 different window sizes) takes only 42% additional time with respect to the computation of a single detection with a window of median size of 29 × 29 pixels (300 s. versus 210 s. for a 800 × 400 pixel image).

B. MCP exploitation

The MCP computation produces a multichannel image (one scale per channel) whose pixels have to be transformed into scalar values in order to provide a change indicator. In order to exploit the information available at all scales, two approaches may be investigated. The first one consists in choosing the best scale for each image pixel. The second one consists in fusing the information available at all scales in order to provide a single change value.

The development of an optimal approach for the exploitation of the MCP may be application dependent. Indeed, multiscale fusion approaches could be tuned to a particular type of change -shape, nature, etc. In this section, two simple yet useful choices will be proposed which yield an improvement in comparison to the performance of a single scale detection.

1) In order to choose the best scale, we will choose the one which produces the highest KLD value. This assumes that this scale is the one that is associated with the largest window inside a homogeneous area with respect to the classes change and no change.

2) The fusion of the multiscale information will be performed by using Principal Component Analysis (PCA).

The first principal component of the MCP multichannel image will be considered as the change indicator. This corresponds to a linear combination of all scales which maximizes the contrast of the final image.

V. EXPERIMENTS WITH SIMULATED DATA A. Data set description Simualtions have been performed to better understand the behavior of the detectors relatively to a given kind of change and a given size of the change area. Since this study focuses on change detection on radar images, a speckle simulation is performed from a map of ground reflectivity. The simulated changes are applied on a small area, drawn as a circle, located in the center of the initial image.

The simulation procedure is based on the radar image formation mechanism. Each pixel is simulated with a given amplitude (coming from a SPOT, NIR band image, normalized to [0, 1]) and thousands phases coming from independent uniform generations in [0, 2π[ to characterize elementary wave scatterers. Taking the square of the modulus of each pixel yields a 1-look intensity image. A 4-look instensity image is obtained by averaging and subsampling two adjacent pixels along lines and rows.

Each simulation of change is applied to the initial image by using a change circle of given size taken from {5, 10, 15, 20}. Once the speckle simulation is performed (independently from one image to another), the speckled-changed images are mosaiked on a 2 × 2 grid as shown on fig. 4(b).

B. Simulation of changes

Three kinds of change were considered: 1) Offset change: fig. 4(c). The initial image is modified by applying an offset value (i.e. a shift) to the inital data. This is a very simple type of change which seldom occurs in reality, but is useful to characterize the behavior of the detectors. 2) Gaussian change: fig. 4(d). The initial image is modified by applying a zero mean gaussian additive noise to the initial data. This corresponds to a change in the state of the surface -field, vegetation. This is the main type of change that one can encounter in medium resolution SAR images.

3) Deteministic change: fig. 4(e). The initial image is modified by pasting values copied from another area of the image itself. This type of change can occur when there is a land-use change, anthropic activities, etc.

C. Results

1) Mono-scale detection: the results of the different detectors for a fixed analysis window size are analyzed.

Figure 5 shows the ROC plots for the case where the change consists in a shift of the reflectivity value (fig. 4(c)). In this case, all 4 detectors are able to detect the changes with high accuracy. There is a slight difference in performance between the pair CKLD -GKLD and the pair PKLD -MRD, but it is difficult to infer general behavior from this result. To draw a preliminary conclusion, for a simple change such as a reflectivity shift, the mean value criterion is efficient enough for good discrimination in the changes, even on speckled images.

Figure 6 shows the ROC plots in the case of a Gaussian change. The change is simulated by the addition of a Gaussian noise to the reflectivity (before speckle simulation). In this case, the mean value of the observed pixels remains approximately the same. It is difficult for this kind of change to be observed by a human operator. However, it is more likely to occur when modifications affect the surface without changing its nature. In this case, even if all the detectors show bad performance, in comparison to the offset case, the MRD and the PKLD are far below the GKLD and CKLD. The bad performance of the MRD is easy to understand, since the zeromean Gaussian noise added to the reflectivity slightly changes the observed mean value. For the PKLD, it can be argued that the type of law in the Pearson system is not very different from the initial case and the main difference is seen through the mean value, thus obtaining the same performance as the MRD. On the contrary, the GKLD assumes a simpler model than the PKLD and is able to take into account the mean and the variance modifications together. Finally, the ability of the CKLD to fit many different types of densities, allows better detection for this difficult type of change. The third type of change is that of a texture change which can occur when there is a land use change, anthropic activities, etc. In this case, as can be infered from figure 7, the mean value of the regions may or may not change and it is therefore interesting to analyze the shape of the density. The Pearson detector can be even worse than the MRD when the model does not fit the data, which is the case in presence of mixtures. 2) Analysis of the MCPs: some collected multiscale change profiles, obtained by applying r CKLD of eq. ( 12) to our data set are analyzed. Four different profiles are presented. They are extracted from a change area of the simulated data set for the case of a deterministic texture change and a radius of 10 pixels. These profiles are labeled as follows: Far for the case where the analysis window is located 30 pixels from the center of the change area; Outside border for a distance of 15 pixels; Inside border for a distance of 7 pixels and Inside centered for a distance of 0 pixel. Fig. 8(a) presents a diagram explaining how the profiles are extracted with respect to the change area and fig. 8(b) presents the profiles themselves.

The Far profile shows low values for small window sizes, and these values increase as the window size increases and it begins to include pixels from the change area. The values decrease for large window sizes, since the window stops including new change pixels while including no change pixels present in all directions. The Outside border profile has a similar behaviour, but the CKLD values are high for small scales since the pixel is nearer to the change area. The Inside border profile shows higher values for the change indicator for small window sizes. Finally, the Inside centered profile shows very high values of the detector for a large interval of window sizes. It is worth noting that the CKLD values are nearly the same for all detectors for the largest window sizes, since, at this scale, all detectors include the same proportion of change and no change pixels.

3) MCP exploitation: in this section, the interest of the use of the MCP is illustrated with respect to the selection of a fixed scale of analysis, (i.e. a fixed window size). The MCP allows the best scale to be selected for each pixel location in the images. Here the maximum of the profile is used as a means to select the appropriate scale.

The maximum of the MCP and 2 different scales, 5 × 5 and 17 × 17 are compared. The small window size is used in order to detect small changes, but its main drawback is that the false alarms may increase in the presence of noise. The larger window size gives a lower false alarm rate since the noise is averaged, and therefore its effect is reduced. But small changes can also be averaged and therefore the detection probability The results of the comparison are presented in figures 9, 10 and 11. As expected, small windows are able to give high detection rates. In the case of radiometric shift, the false alarm rates are low, for a given detection probability, since the type of change is easily detected by computing the mean value over a few pixels only. However, when more complex changes occur (figures 10 and 11), the false alarm rate is very high at a given detection probability. Another interesting effect can be observed in figures 9 and 11, where for the large window sizes, the false alarm rate increases without an increase in the detection probability. This is due to the fact that when the window is too large for the small changes and not as large as the larger changes (see the mask in fig. 4(b)), the new detections induce false alarms only in the neighborhood of the small changes.

In addition, the MCP gives results which do not suffer from these drawbacks without the constraint of choosing a window size without prior information on the size of the changes in the images. before and after the eruption of the Nyiragongo volcano (D.R. of Congo) which occurred in January 2002, were used. Fig. 12 shows the two images to be compared and a change map produced using ground measures. The images have a ground resolution of 10 m and cover an area of 4 km by 8 km. The images were ortho-rectified by IGN-F, the French National Geographic Institute, to a UTM35S projection, which was the same as the one used for the reference map. No filtering or calibartion was applied to the data. The 16 to 8 bit conversion was performed using a 3σ thresholding followed by a linear intensity rescaling. It is worth noting that the image resampling applied in the ortho-registration step modifies the local statistics of the image. Indeed, the image resampling implies local image interpolation which is based on approximate interpolators. A bicubic interpolation was used in this case. This type of filter has a smoothing effect which depends on the local shift [START_REF] Inglada | Analysis of artifacts in sub-pixel remote sensing image registration[END_REF]. Because of these radiometric artifacts introduced during the geometric pre-processing, the theoretical models for SAR statistics may not hold locally. The area at the bottom right hand corner of the ground truth mask corresponds to an area where a severe mis-registration exists due to the lack of a proper digital terrain model. Finally, one has to take into account the fact that the ground truth is not perfect and that all results presented in this section should be analyzed rather in a relative manner -one detector with respect to another -rather than in an absolute one -absolute value of detection probabilities.

A. Change indicator

The comparisons between the result coming from the classical image intensity ratio and the method proposed in this paper are shown in fig. 13. Fig. 14 gives the ROC plots using the ground truth of fig. 12(c). It shows that the use of KL approximation by the Edgeworth series outperforms any other methods such as model-based (Gaussian-based or Pearsonbased) KL distance, or the ratio measure. As stated in the introduction, a miss-detection behaviour of this detector can be observed because it uses the mean pixel values only. It is interesting to underline the fact that the ratio criterion is not always worse than pdf-based criteria. In fact, a density model has to fit the data in order to yield pertinent results.

For a detection probability below 0.3, it is more interesting to use the ratio criterion instead of a model-based one (by using Gaussian or Pearson assumption) in this example, even if a better change detection could have been expected by using Gaussian or Gamma laws coming from local analysis of the two Radarsat images.

This point confirms that it is more interesting, for operational use, to consider a more flexible pdf approximation by using the Edgeworth series instead of a pdf parameterization. The cumulant-based approximation may give equivalent results to the Pearson-based approximation if the estimated cumulants correspond to a pdf belonging to the Pearson system of distributions, even though it may be less appropriate in the case of heavy tailed distributions (single look data). If cumulants of orders 3 and 4 vanish, the Edgeworth series is equivalent to a Gaussian model. If the variance of X and Y are equivalent, the Edgeworth series yields the same behavior as the ratio measure. However, when the local observations X and Y to be compared do not fit an a priori model, the Edgeworth series becomes a more suitable tool. Fig. 15 draws the minimum distance of ROC curves to the point (P d = 1, P f a = 0). It is an interesting point of view to evaluate the threshold to be applied to obtain the best tradeoff between detection and false alarms. The best value of the threshold is to be found at the minimum of the curves. Fig. 15 shows that this minimum is lower -and therefore, more interesting -for the Edgeworth series than for the Pearson measure or the ratio detector.

When no ground truth is available, the end-user has no a priori knowledge to set the value of the threshold. In this case the Pearson measure seems to be better since a trivial value of zero could be used (i.e. pixels with values greater to 0 may be considered as a change). Unfortunately, simulations and comparisons with other sets of images have shown that this trivial threshold is very sensitive to noise and fluctuates. The same observations about sensitivity hold for the ratio measure. On the contrary, the cumulant-based measure takes its minimum for a wider range of values. Therefore, a threshold chosen a priori from the interval [40, 50] gives an almost optimal change map for all cases.

B. Multiscale change indicators

As stated in section IV-B, our goal here is not find the optimal way of exploiting the MCP, but only show the interest of the concept with simple examples. The results presented here use an MCP with window sizes ranging from 29 × 29 to 51 × 51.

In order to select the appropriate analysis window for each pixel in the image, we will choose the maximum of the MCP. The resulting change image is shown in fig. 16(a). Fig. 17(a) presents the histogram of the sizes of the selected analysis windows when using the maximum of the MCP. It is interesting to observe that there is a high variability of window sizes, meaning that no trivial choice exists, like for instance choosing the largest window in order to increase the number The ROC plots of fig. 18 show that this simple strategy improves the results with respect to the case where the 35 × 35 window was used.

As an approach to multiscale fusion, we propose here to use the first principal component of the stack of multiscale detection images. The obtained change image is presented in fig. 16(b). The ROC plot of fig. 18 shows that this approach also provides better performance than the monoscale detector.

VII. DISCUSSION AND CONCLUSION

In this paper a new similarity measure between images has been introduced in the context of multitemporal SAR The proposed original cumulant-based detector has been shown to have a more robust behaviour than other detectors in terms of receiving operator characteristics. The two simple yet useful schemes for the exploitation of the multiscale change profile provide better performance than the monoscale detector.

The main advantages of the proposed approach are the following: our detector needs only the computation of the first 4 statistical moments and can deal with a great variety of pdfs; the multiscale change profile provides change information over a wide range of scales at very low computation cost. Some improvements could be done in order to use this approach with single-look images, where the heavy tailed distributions may need other statistical models. The use of Gamma distributions instead of Gaussian for the series expansion seems to be a good starting point. Some questions still remain open about the use of multiscale change profiles. Indeed, it would be interesting to analyze if we could establish a classification of the profiles and thereby derive useful information, not only about the scale of the change, but also about its type. This task could be carried out by visual inspection, but automatic clustering techniques, like for instance the Self Organising Map [START_REF] Kohonen | The Self-Organizing Map[END_REF] could be used. The parametric modeling of the profiles by projection on an orthogonal basis could be envisaged.

Another issue remaining is the automatic thresholding of the change images. Whether it be for the case of a single scale, or for the case of a multiscale analysis, the statistics of the change indicators could be used in order to propose adaptive Bayesian thresholding techniques, as done in [START_REF] Bazi | An Unsupervised Approach Based on the Generalized Gaussian Model to Automatic Change Detection in Multitemporal SAR Images[END_REF].

Finally, direct classification of multiscale profiles by using Support Vector Machines seems an appropriate choice for the production of binary change maps in the case of supervised analyses. This approach has successfully been applied to the classification of hyper-spectral images [START_REF] Melgani | Classification of hyperspectral remotesensing images with support vector machines[END_REF].

All these aspects will be studied in future work.

Fig. 1 .

 1 Fig. 1. Block diagram for a classical change detection processing chain.

Fig. 2 .

 2 Fig. 2. Approximation of a histogram, coming from a 50 × 50 Region of Interest (ROI), using three different strategies. The Pearson fitting yields a Beta distribution of the first type.

Fig. 3 .

 3 Fig. 3. Example of pdf estimation update by increasing sample set from a window of size 9 × 9 to 17 × 17. The histogram has been estimated with a 17 × 17 window.

Fig. 4 .

 4 Fig. 4. Simulated data set.

Fig. 5 .

 5 Fig. 5. ROC plot comparison of the 4 detectors for a simulated change consisting in an offset on reflectivity.

Fig. 6 .

 6 Fig. 6. ROC plot comparison of the 4 detectors for a simulated change consisting in a Gaussian random modification of the reflectivity.

Fig. 7 .

 7 Fig. 7. ROC plot comparison of the 4 detectors for a simulated change consisting in a deterministic modification of the reflectivity.

Fig. 8 .

 8 Fig. 8. Typical examples of multiscale change profiles obtained from the Edgeworth approximation of the KL distance

Fig. 9 .Fig. 10 .

 910 Fig. 9. ROC plot comparison between MRD -2 scales -CKLD -2 scales -and MCP -maximum of the profile -for a simulated change consisting in a reflectivity offset.

Fig. 11 .

 11 Fig. 11. ROC plot comparison between MRD -2 scales -CKLD -2 scales -and MCP -maximum of the profile -for a simulated change consisting in a deterministic modification of the reflectivity.

Fig. 12 .

 12 Fig. 12. Data and ground truth for the Nyiragongo volcanic eruption of January 2002.

Fig. 13 .

 13 Fig. 13. Change detection: comparison between the different change indicators using the same window size (35 × 35 pixels).

Fig. 14 .Fig. 15 .

 1415 Fig. 14. ROC plots for the different detectors. The cumulant-based Kullback-Leibler detector outperforms all other detectors. The Pearson-based detector gives results identical to the classical mean ratio. The Gaussian-based detector shows the worse behavior.

Fig. 17 .Fig. 18 .

 1718 Fig. 16. Change detection results obtained with the MCP.

Grégoire

  Mercier (M'01) was born in France in 1971. He received the Engineer Degree from the Institut National des Télécommunications (INT, France) in 1993 and his Ph.D. degree from the University of Rennes (France) in 1999. Since 1999, he has been with the Ecole Nationale Supérieure des Télécommunications de Bretagne (ENST Bretagne, Brest, France), where he is currently an Associate Professor in the Image and Information Processing department. His research interests are in remote sensing image compression and segmentation, especially in hyperspectral and Synthetic Aperture Radar. Actually, his research is dedicated to change detection and combating pollution.

In the case where a sequence of several images is to be processed, the approaches presented in[START_REF] Bruzzone | An advanced system for the automatic classification of multitemporal SAR images[END_REF],[START_REF] Lombardo | Maximum likelihood signal processing techniques to detect a step pattern of change in multitemporal SAR images[END_REF] may be applied.


T Lelong

P Thomas

I Magne

Riccardo Scorretti

Noël Burais

F Piriou

A critical

teaching and research institutions in France or abroad, or from public or private research centers.

I. INTRODUCTION

The finite element method (FEM) is used for evaluating the induced currents into the human body by extremely low frequencies (ELF) electric and magnetic fields. Many formulations, developed under different assumptions, are available. The most common approximations which are done are: i) the reaction field due to the presence of eddy currents into the human body may be neglected, ii) displacement currents inside the human body are negligible compared with conduction currents, iii) both exposure cases can be considered separately (electric and magnetic fields are decoupled) [START_REF] Hoang | Numerical Dosimetry of Induced Phenomena in the Human Body by a Three-Phase Power Line[END_REF]. It is observed that, with the same computational phantoms, these formulations produce somehow different results [START_REF] Ducreux | Benchmarks for computing induced currents in the human body by ELF electric and magnetic fields[END_REF]. The purpose of this work is to examine in details the validity of these commonly used approximations. To this aim, we have considered a general formulation, where none of these approximation is used, so that all the other studied formulations can be deduced from this one by negliging some terms.

II. FORMULATIONS

A. Complete formulation

We start from a general -→ Aφ formulation, where displacement currents are taken into account. All the -→ E -conform formulations can be deduced by this one, namely the classical -→ Aφ formulation [START_REF] Dular | Magnetostatic and magnetodynamic mixed formulations compared with conventional formulations[END_REF], the φ --→ A formulation [START_REF] Hoang | Numerical Dosimetry of Induced Phenomena in the Human Body by a Three-Phase Power Line[END_REF], and the dielectric formulation.

The -→ A -φ formulation in the frequency domain and without any approximation inside a bounded domain Ω ⊂ R 3 reads:

--→ curl (ν --→ curl - → A ) -(σ + iωǫ)(-iω - → A - --→ grad φ) = 0 div ( σ iω + ǫ)(-iω - → A - --→ grad φ) = 0 (1)
where

- → B = --→ curl - → A , - → E = -iω - → A - --→ grad φ.
Equation 1 can be written in a weak form as:

                   Ω (ν --→ curl - → A ). --→ curl -→ A ′ dΩ+ Ω (σ + iωǫ)(iω - → A + --→ grad φ). -→ A ′ dΩ = Γ (ν --→ curl - → A ) ∧ - → n . -→ A ′ dΓ Ω ( σ iω + ǫ)(iω - → A + --→ grad φ). --→ grad φ ′ dΩ = Γ ( σ iω + ǫ)(iω - → A + --→ grad φ). - → n .φ ′ dΓ (2)
where Γ = ∂Ω and -→ n the exterior unit normal to Γ. Boundary conditions have to be adapted to the studied exposure in order to take into account a source term:

-for a magnetic exposure:

( - → A ∧ - → n ) |Γ fixed , (-iω - → A - --→ grad φ) |Γ . - → n = 0 (3)
-for an electric exposure:

( - → A ∧ - → n ) |Γ = 0 , φ |∂Ω fixed on Γ d ⊂ Γ (4)

B. Magnetic formulations

For magnetic exposure, two main approximations are considered: displacement currents and reaction field from the human body are neglected. The classical -→ Aφ formulation only neglects displacement currents (ωǫ ≪ 1):

                 Ω (ν --→ curl - → A ). --→ curl -→ A ′ dΩ+ Ω σ(iω - → A + --→ grad φ). -→ A ′ dΩ = 0 Ω ( σ iω + ǫ)(iω - → A + --→ grad φ). --→ grad φ ′ dΩ = 0 ( - → A ∧ - → n ) |Γ fixed , (-iω - → A - --→ grad φ) |Γ . - → n = 0 (5)
The formulation φ --→ A neglects displacement currents and reaction field. Therefore computational domain Ω can be limited the human body, and

- → A = -→
A s (which is assumed to be known) is used as source term in the second equation:

Ω σ iω (iω -→ A s + --→ grad φ). --→ grad φ ′ dΩ = 0 (6)

C. Electric formulation

The "dielectric" formulation [START_REF] Hoang | Numerical Dosimetry of Induced Phenomena in the Human Body by a Three-Phase Power Line[END_REF] only considers the electric field and totally disregards the magnetic field by imposing -→ A = 0, so that φ must verify:

     Ω ( σ iω + ǫ) --→ grad φ. --→ grad φ ′ dΩ = 0 φ |∂Ω fixed on Γ d ( σ iω + ǫ) --→ grad φ |∂Ω . - → n = 0 on Γ\Γ d (7)
Starting from the same "complete" formulation discretized by using mixed finite elements, comparisons with dedicated formulations are possible in order to assess the impact of approximations on the solution.

III. MAGNETIC EXPOSURE

For each exposure case, the results are obtained with the ANSOFT human body model (Fig. 1), which is composed of 29 different tissues with different conductivities [START_REF] Gabriel | The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues[END_REF] and is discretized by using N E = 7 • 10 5 tetrahedra. The in-Fig. 1. Human body model duced currents are calculated with the "complete" (2) and the classical -→ Aφ formulation (5). In the latter case, a much faster convergency is observed. In order to compare these two formulations, we define the normalized discrepancy:

∆J = || - → J 1 ||-|| - → J 2 || max|| - → J || × 100%, where - → J 1 et - → J 2 are computed
with these two formulations, and max|| -→ J || is the maximum current density computed with (2). For each tetrahedron, the quality q defined as in [START_REF] Ducreux | Benchmarks for computing induced currents in the human body by ELF electric and magnetic fields[END_REF] and the discrepancy in the barycenter are computed, so we obtain N E couples {∆J e , q e } e=1...N E . These couples are ordered in ascending order with respect of q e . Tetrahedra are grouped in classes, the quality of which spans within [10k%, 10(k + 1)%] (with k = 0 . . . 9) and the average value of the discrepancy ∆J is computed for each class. If Fig. 2 the average ∆J is plotted for each class together with the quality q (the continuous line). One observes that ∆J is correlated with the quality of the elements (the bigger is the quality q, the smaller is the average discrepancy ∆J ). However, within a same class of elements very large variations of ∆J are observed. The same approach has been used for comparing the classical -→ Aφ formulation (5) with the φ--→ A formulation (6). Both formulations neglige displacement currents, but φ--→ A formulation negliges also the reaction field due to the human body. The average discrepancy is depicted in Fig. 3). In both cases the computed discrepancies are quite negligible: this confirms the validity of the approximations employed in (5) and (6).

IV. CONCLUSION

The calculation of induced currents in human body is traditionally performed by using different "dedicated" formu- lations, which introduce some approximations with respect of the original problem. These numerical experiments show that the error introduced by these approximations is negligible. Furthermore, high discrepancies between formulations are correlated with low quality elements -even if this correlation is rather weak, in that it holds only for the average discrepancy, but not for each element. In the full paper the case of exposure to electric field will also be adressed.

Fig. 2 .

 2 Fig. 2. Average discrepancy ∆J between the complete and classical-→ A -φ formulations, and quality of the elements.

Fig. 3 .

 3 Fig. 3. Average discrepancy ∆J between the classical and dedicated -→ A -φ formulations, and quality of the elements.