Pierre Caserta
email: pierre.caserta@loria.fr

Olivier Zendra
email: olivier.zendra@inria.fr

JBInsTrace: A Tracer of Java and JRE Classes at Basic-Block Granularity by Dynamically Instrumenting Bytecode

Keywords: Java, Tracing, Dynamic, Bytecode, Instrumentation, Profiling, Program Analysis

Understanding what happens during the runtime of a Java program is difficult. Tracking runtime flow can bring valuable information for program understanding and behavior analysis. Polymorphism, thread concurrency or even simple facts like the number of method invocations and the number of executed bytecodes are valuable information to track, but are difficult to compute outside the Java Virtual Machine (JVM) on running programs. In this paper, we present JBInsTrace, a new tool that instruments and traces Java bytecode. It produces static information about source code and a very fine grained trace of Java software execution, combining them to allow detailed analysis of the runtime. Our tool differs from others because it does not only trace program classes but also JRE classes, and does so at basic block level, without altering the JVM and without statically modifying class files. We explain JBInsTrace design, focused towards efficiency, which results in reasonable performance penalty.

Introduction

Dynamic analysis gives information on a particular execution of software and brings significant information to analyze and understand program behavior [START_REF] Ball | The concept of dynamic analysis[END_REF]; [START_REF] Cornelissen | A systematic survey of program comprehension through dynamic analysis[END_REF]), helping focus on modules which play a predominant role during the runtime of a program. Dynamic information such as reflection and dynamic class loading can be statically approximated [START_REF] Christensen | Precise analysis of string expressions[END_REF]; [START_REF] Sreedhar | A framework for interprocedural optimization in the presence of dynamic class loading[END_REF]), but dynamic analysis is needed for a real accuracy [START_REF] Hirzel | Pointer analysis in the presence of dynamic class loading[END_REF]). On the other hand, static analysis gives insight about the source code itself and is thus useful to manage its quality. Static analysis tools are more common than dynamic analysis tools; we think this is because static analyses are easier to implement.

One important goal of JBInsTrace1 tool and this paper is thus to provide a technique to perform fine-grained dynamic analysis of Java software with reasonable performance penalty. To do this, we divide profiling in two stages. The first stage consists only in tracing the exact runtime control flow and saving static information about the source code. The second stage performs an off-line analysis of the trace (a.k.a. post mortem trace analysis) and of the static information from the previous step, in order to create the call graph and to compute complex dynamic metrics [START_REF] Dufour | Dynamic metrics for Java[END_REF]; [START_REF] Arisholm | Dynamic coupling measurement for object-oriented software[END_REF]). One advantage of this two-stage technique is that metrics are not computed at runtime. In addition, when new metrics have to be computed, the instrumentation does not need to be modified.

Most of the work on dynamic analysis focuses on tracing runtime at the method level. Even though some works are very consistent in tracking method flow and do provide a lot of valuable information, we wanted the extra accuracy of a finer level. Our analysis of Java programs is thus made at the basic block2 level, which means that execution is tracked even within methods. We can also analyze call sites within basic blocks.

Our tracing tool JBInsTrace instruments program classes as well as JRE classes because we want a full coverage of the whole runtime of the program, to obtain the most complete trace of the control flow, which we think is necessary to fully understand the behavior of a program. Of course, JBInsTrace allows switching off the tracing of JRE classes when desired.

The emergence of new paradigms such as Aspect-oriented programming (AOP) makes it easier to do rapid prototyping of profilers, debuggers, tracers, and reverse engineering tools [START_REF] Villazon | Flexible Calling Context Reification for Aspect-Oriented Programming[END_REF]). Indeed, dynamically adding new functionalities on running programs by modifying classes at runtime through bytecode instrumentation is a convenient way to add behavior to programs. Our tracing tool and technique rely on this paradigm to instrument Java programs at strategic join points and extract information from their runtime. However, the main difficulty in instrumentation lies in finding an efficient design for the profiler [START_REF] Mock | Dynamic analysis from the bottom up[END_REF]).

In this paper, we present the design and some implementation details of our dynamic instrumentation and tracing technique. In section 2 we give an overview of the technique and tool we created to dynamically instrument Java classes. Section 3 details our technique. Then, in section 4 we discuss the static and dynamic outputs of the tracer and how they are analyzed to compute valuable information. Section 5 shows performance results that validate our approach to take user-perceived efficiency into account from the very beginning. Finally, section 6 concludes and presents future work.

Overview of the JBInsTrace tracing technique

Since version 1.5, Java provides services that allow Java agents to instrument programs running on the JVM. The agent is executed in the same JVM, loaded by the same class loader, and governed by the same security policy and context as the program. We thus intercept Java classes when they are loaded to instrument them "on-the-fly", even if the program has custom class loaders.

The design of the JRE is complex, with a lot of coupling between classes. Mechanisms to realize simple operations thus imply many underlying method calls and execute many bytecodes. This make understanding what happens when calling a JRE method difficult. As an example, the first program most people write is "hello world", coded in Java as: System.out.println("Hello world ...");. The runtime of this Java program executes around 42550 bytecodes, a major part of which are program initialization bytecodes. The println method call alone generates about 160 method calls and executes around 2660 bytecodes. I/O methods are known to be complex, but this example shows the importance of JRE methods execution.

A specificity of our profiling technique is that it performs a fully dynamic instrumentation on all Java classes, including core JRE classes (see section 3.1). Classes loaded because of reflection are also treated like the ones automatically loaded by the JVM: code Class.forName("X") causes class X to be loaded and initialized, hence instrumented.

Since JRE classes may be used by both the observed program and the tracer, instrumented JRE classes could lead to trace pollution by the tracer. We explain our mechanism to avoid this pollution in detail in section 3.3. However, like any Java instrumentation technique injecting bytecodes into classes, ours still impacts timing and thread scheduling because of the extra tracing bytecodes to execute, even though they are not put in the trace.

Runtime flow is traced at basic block level. This granularity gives insight about intraprocedural control flow. With our JBInsTrace trace, we know how many times every single bytecode is executed. Since call sites are traced, information about polymorphism can be extracted during the analysis step by comparing the call site static type and the method called at runtime. Native method calls and thread switches are detected as well and feed additional information to our tracer (see section 3). In fact, JBInsTrace provides enough information to precisely re-build the whole runtime flow of the profiled program and thus perform very precise dynamic analyses.

Figure 1 summarizes the global functioning of our profiling technique. When a class is loaded by the ClassLoader, a Java agent catches the bytecode of this class on the fly (see section 3). The Instrumentor parses the class, injects new bytecodes and saves static information about basic blocks, methods and classes (see section 3.2). The injected bytecode adds a new (tracing) behavior, without altering the functional semantics of the original bytecode but with an impact on timing and scheduling (see section 3.2). At runtime, the injected bytecode calls back our Tracer (see 3.2.3), which must take care not to pollute the trace with tracer events (see section 3.3). At the end of the runtime, all the static information and the dynamic execution trace are available in files, ready to be processed by any appropriate analyzer (see section 4).

Dynamic instrumentation of classes

Our dynamic instrumentation transforms only classes which are actually used during a specific execution, without altering the original .class bytecode files on disk. We made this choice because static instrumentation would have implied transforming all classes that might be used by the program, whereas only a small percentage of this huge set of classes would actually be used at runtime. Moreover, static instrumentation requires keeping both an instrumented and a non-instrumented version of class files on disk, more than doubling the space. Since classes appear in JAR files, it would have been constraining to first extract classes, then instrument them and finally re-JAR them, and do so at each new version of a particular class. Furthermore, static instrumentation cannot be done on dynamically created classes. For all these reasons we decided to instrument code dynamically.

The following three subsections detail our technique. Section 3.1 explains how the Java agent works. Section 3.2 details how we instrument classes. Section 3.3 shows how we avoid polluting the trace of the running program with our JBInsTrace tracer.

Instrumenting classes with a Java agent

A Java agent, deployed as a JAR file, is a service that allows instrumenting programs running on the JVM through the -javaagent command line option. Our JBInsTrace tracer has been tested with the Java HotSpot VM and the IBM z/VM, but it should work on any JVM that provides the agent service.

The premain method is a method called by the Java agent before the main method of the application but after the JVM has initialized. JBInsTrace uses premain to register its class transformer (Instrumentor) which will be applied on all future class definitions.

However some classes are loaded during the initialization phase of the JVM (we name them "core JRE classes"), hence before our class transformer is registered. Indeed, the JVM follows a lazy loading process, which means that classes are loaded only when needed. When the JVM starts without the Java agent, it loads the initial class, initializes it, and invokes its class method main(String[]) [START_REF] Lindholm | Java virtual machine specification[END_REF]). With the Java agent, the initialization process is changed and all the classes needed for the execution of the premain method are also loaded by the JVM during its initialization. In our case this adds around 350 core JRE classes, mainly from the java.lang package.

We instrumented such core JRE classes with the retransform method from Sun-Oracle implementation, but had to overcome some limitations: the instrumentation must not add, remove or rename fields or methods, change the signatures of methods, or change inheritance (Sun Microsystems (2008)).

This differs from other existing dynamic techniques which usually introduce new methods acting as wrappers and adding extra arguments to methods to track the calling context [START_REF] Villazon | Flexible Calling Context Reification for Aspect-Oriented Programming[END_REF]). With such techniques, core JRE classes have to be treated differently. The solution proposed by [START_REF] Binder | Advanced Java bytecode instrumentation[END_REF]) in their "JP" tool was to statically instrument the .class files. They have reimplemented their tool, now called "JP2" [START_REF] Sarimbekov | Complete and Platform-Independent Calling Context Profiling for the Java Virtual Machine[END_REF], which is much more advanced and does not suffer from many of the restrictions the older version imposed. The new tool does not perform structural change on classes anymore but still has to statically modify the Thread class. Instead of adding new parameters to methods, they add a new field to the Thread class and use it as a reference to save the calling context tree of the running program. To trace native calls, they still need to statically instrument all the JRE classes.

Obviously, with bytecode instrumentation, no tracing occurs during the JVM initialization phase. The only way to do such tracing would be by modifying the JVM.

Bytecode instrumentation

Bytecode instrumentation injects new bytecodes into class. Several tools and libraries exist to do this. We chose the ASM framework (Bruneton et al. (2002a,b); [START_REF] Bruneton | Asm 3.0 a java bytecode engineering library[END_REF]; Kuleshov (2007a)) which provides libraries to parse Java classes and add bytecode, allowing complex transformations, with very low memory requirements. According to [START_REF] Kuleshov | Using the asm framework to implement common java bytecode transformation patterns[END_REF]), it is much faster than other bytecode transformers such as BCEL [START_REF] Dahm | The bytecode engineering library[END_REF][START_REF] Dahm | Byte code engineering with the bcel api[END_REF]), SERP [START_REF] White | SERP, an Open Source framework for manipulating Java bytecode[END_REF]) and Javassist [START_REF] Chiba | Javassist: Java bytecode engineering made simple[END_REF]).

In subsection 3.2.1 we show how we extract static information from loaded classes. Then in subsection 3.2.2 we describe where we inject bytecodes in classes and what the purpose of these new bytecodes is.

Static information extraction

When classes are loaded, JBInsTrace parses them and performs three operations:

• unique identifiers are assigned to each class, method and basic block.

• static information is extracted from the class source code.

• the class is instrumented with new bytecodes.

These identifiers are used to identify parts of the source code, when the trace is saved on disk. All the important static information about classes, methods and basic blocks is also saved. This information is then used during the analysis phase, when identifiers are read from the trace and related the corresponding static information.

The static information is saved in XML (with the Xerces library) or CSV format. In this paper we use XML for the static information because it is more readable. Here is an example of the static information saved for a basic block: <basic_block id="5095440" metrics="4 0 0 0 1 1 0 1 0 0 0 0 0 0 0 1 0 1 0"> <call_site_list> <call_site call="java/nio/CharBuffer:hasRemaining()Z" type="INVOKEVIRTUAL"/> </call_site_list> </basic_block>

The metrics field is a list of basic block metrics. The latter contains the numbers of: bytecodes, array read operations, array write operations, floating point operations, reference loads, control bytecode instructions, changing control bytecode instructions, read operations on object fields defined in the class, write operations on object fields defined in the class, read operations on primitive field defined in the class, write primitives on object field defined in the class, read operations on object field defined in another class, write operations on object fields defined in another class, read operations on primitive fields defined in another class, write primitives on object fields defined in another class, invocations of another class, new bytecode operations, virtual calls and finally static calls. In future versions of JBInsTrace we plan to save the entire bytecode of each basic block instead of just a list of metrics. More information can be found on the JBInsTrace web page (http://www.loria.fr/˜casertap/jbinstrace.html).

Code segmentation

JBInsTrace traces 3 kinds of events: when a method starts and ends, like any other tracing tool, but also when a basic block is executed. Because a basic block contains only one jump instruction at its end, it is executed from beginning to end, unless an exception occurs. Section 4.2 explains how JBInsTrace follows the runtime flow with exceptions.

An event is recorded as a single integer to maximize performance. Figure 2 shows how we divide a 32-bit integer, named event number, to encode 3 pieces of information:

• The 2 most significant bits encode the event type identifier : method start, method end, basic block execution.

• The 19 bits in the middle encode the unique identifier of the method.

• If the event type is the execution of a basic block, the 11 least significant bits represent the unique identifier of the basic block within the method (unused otherwise).

Note that the event number can be negative, since the most significant bit is used for event type encoding. Together the method identifier and the basic block identifier form a unique identifier for the basic block. In fact, every event number is unique. This coding can identify 524288 different methods that can have at most 2048 basic blocks each. We think these limitations are reasonable. Indeed, ArgoUML, which is a large open source UML modeling tool loading around 5100 classes at startup, counts around 12000 instrumented methods, with at most 180 basic blocks in a method. Moreover, trying about 100 runs of ArgoUML with many different inputs led to a cumulated number of used methods of 50000. If needed, we could switch to using a 64 bit long integer to lift these limitations, but this would increase memory usage.

Instrumented code

Java programs are multi-threaded so runtime flow has to be put in the context of the actual thread where execution takes place. Our tool thus also records for each event (method start, method end, basic block start) the running thread number. To do it, every time an event occurs, the instrumented code calls the Tracer (via the notifyEvent static method) and passes it the associated event number and the current thread id. The Java code corresponding to our instrumented bytecode is:

Tracer.eventNotifier(eventNumber, (int)Thread.currentThread().getId()); Figure 3 shows an example of how a piece of source code is instrumented. The left part 3(a) is the Java code before instrumentation, while the right part 3(b) represents the source code after instrumentation.

As illustrated in figure 3, instrumentation is done:

• at the beginning of each method. There, eventNumber is composed of event type '01' binary meaning 'method start', followed by the method identifier.

• at the end of each method. There, eventNumber is composed of event type '10' binary meaning 'method end', followed by the method identifier.

• at the beginning of each basic block. There, eventNumber is composed of event type '11' binary meaning 'basic block start', followed by the method identifier, followed by the basic block identifier within the method.

With this technique, the numerous calls to the Tracer cost time, but we observed good performance results (see section 5). Since it adds no new method nor any extra parameter, it preserves the semantics of the original program and the injected code is the same for all classes, including JRE core classes. However, like any technique which injects new bytecodes into classes, ours impacts timing and thread scheduling.

Saving dynamic instance types and class loading information

Class identifiers have two purposes: depict the dynamic type of an instance on virtual and interface call sites, and identify the loaded classes.

To depict the dynamic type of instances on virtual and interface calls, the instrumented code is slightly different. When a method starts, the receiver dynamic type is obtained using the reflexive method getClass():java.lang.Class on this [START_REF] Sundaresan | Practical virtual method call resolution for java[END_REF]). Then the instrumented code calls method eventNotifierWithInstanceType of the Tracer, which saves the corresponding class identifier in the output trace.

The pseudo-Java code inserted at the beginning of methods is thus:

Tracer.eventNotifierWithInstanceType(eventNumber, (int)Thread.currentThread().getId(), this.getClass());

When a class is loaded, the Instrumentor notifies the Tracer which saves this information in the trace buffer. The event number is composed of event type '00' binary, followed by the identifier of the loaded class.

Trace neutral tracer

Because both program and tracer run on the same JVM, JBInsTrace uses many JRE classes which are themselves instrumented. We must separate the execution of the tracer and the execution of the program to avoid having traces of the tracer in the trace files. To solve this issue, we resorted to a system of boolean flags that switch tracing on or off. These flags avoid creating infinite recursion between JRE classes and the Tracer.

First, instrumented bytecode calls the Tracer using the notifyEvent method. The Tracer changes the authorization boolean to false to forbid tracing. Then the code of our Tracer is executed. This code relies on the JRE classes, which are themselves instrumented, causing these JRE classes to call the Tracer again. This time the authorization flag is false and forbids tracing, since we are in the Tracer. At the end of the notifyEvent method, the authorization flag is switched back to true to allow tracing again.

The same boolean mechanism is used by the Instrumentor to forbid its own tracing, as shown in figure 4.

The principle is trivial but care must be taken because Java programs are multithreaded, so we must constantly track the number of the current thread and manage one authorization flag per thread.

The Tracer

Our Tracer class is a singleton which can be called by any class of the program. This Tracer simply buffers the list of observed events (integers) before writing them to disk. The Tracer can manage the event buffer using two modes:

Mode 1one unique, shared trace to save the pairs (eventNumber, threadId) for all threads.

Mode 2one trace of events per thread (no sharing),

In mode 1 (one unique trace for all threads), synchronization is mandatory to prevent pollution and/or dead locks. However, this may drastically slow down the tracer because of waits to reach the critical zone. The tracer in mode 1 (one unique trace) is thus much slower than in mode 2 (one trace by thread), but it preserves the thread concurrency scheduling.

Note that the Tracer is called intensively by the instrumented program so methods notifyEvent and notifyEventWithInstanceType must be as fast as possible. In fact the only operations in these methods are:

• check the value of the authorization boolean that corresponds to the current thread (see figure 4)

• save the (event number, thread ID) pair in the shared buffer (in mode 1), or save only the event number in the buffer according to the thread ID (in mode 2)

• in notifyEventWithInstanceType an extra operation consists in saving the dynamic instance type in the appropriate buffer.

Exploiting the program trace and static information

Our JBInsTrace tracer provides files that contain the exact control flow of each thread at basic block level, plus static information about these basic blocks. We exploit this information with a trace analyzer (see the rightmost part of figure 1). In this section we explain how to develop such a trace analyzer for JBInsTrace traces. In subsection 4.1 we explain how to build a new JBInsTrace trace analyzer. Then in subsection 4.2 we describe how the tracing of exceptions is handled in JBInsTrace. Finally, subsection 4.3 details how we detect native calls.

Building a JBInsTrace Trace Analyzer

The first step is to read the static information and construct the data structure representing the observed program. An interesting meta-model of Java programs can be found in [START_REF] Hoffmann | A case study for program refactoring[END_REF]): since it takes basic blocks into account it is adapted to model the static information from JBInsTrace. We used a similar meta-model in our analyzer. Then, a map that links the integer identifiers of the JBInsTrace trace (see section 3.2.2) to their corresponding elements in the meta-model instance has to be built.

The functioning of our analyzer is as follows. We relate each event of the trace with the static information saved. For instance if we find that basic block b has been executed in the trace, we search for the static information of basic block b and see what kind of bytecode basic block b has. This technique allows to re-simulate the whole execution call stack and to compute very precise dynamic metrics. For example, for each method call event in the trace, we can compare the static type of a virtual call site (in the static information) and the dynamic type of the instance (in the trace), to compute metrics related to polymorphism such as the metric that counts how many times receiver types change on call sites.

The authors of [START_REF] Dufour | Dynamic metrics for Java[END_REF]) have defined several interesting dynamic metrics which are implemented in our analyzer. These metrics can be used to categorize programs according to their dynamic behaviour in five areas: size, data structure, memory use, concurrency, and polymorphism. Dynamic metrics related to program size and structure, use of data structures, use of polymorphism, memory footprint and concurrency are computable with our analyzer. We also use it to compute precise call graphs of programs and libraries, to feed in our advanced visualization tool [START_REF] Caserta | 3D Hierarchical Edge Bundles to Visualize Relations in a Software City Metaphor[END_REF]).

Tracking exceptions

The only way the execution flow of a basic block may be disrupted is by an exception, which causes program control to jump to an exception handler. In some cases, the call stack has to be popped until the corresponding catch block is found. The try/catch structure thus has to be recorded, to have each basic block aware of the basic block identifiers it can potentially jump to. This information allows simulating the call stack during the analysis phase.

Here is an example of the static information saved on a basic block with two possible exception handlers: <basic_block id="5095441" metrics="7 0 0 1 0 1 0 1 0 0 0 0 2 0 0 1 0 1 0"> <exception_handler_list> <exception_handler basic_block_id="3764234" type="java/io/InterruptedIOException" /> <exception_handler basic_block_id="3764235" type="java/io/IOException" /> </exception_handler_list> </basic_block>

Tracking native calls

With our tracing technique, it is the called method that informs the trace about its beginning and its end, not the caller. Native method calls are traced differently, because native methods are not coded in Java, hence are not instrumented. Moreover our technique does not add any new method to classes, because of the limitation on core JRE classes. Using wrapper methods to detect native calls is thus impossible for us. As a matter of fact, we use a very simple technique to detect native calls. When analyzing runtime, we still follow static information about each basic block. When a method call is found in the basic block static information, we check whether the corresponding call is present in the actual dynamic trace. If not, we know this call is a native one. Exceptions thrown in native methods are tracked like other exceptions (see section 4.2).

Performance

In this section, we show the time performance of our JBInsTrace v0.6 tracer on several Java benchmarks. We first present results obtained on well-known Open Source Java software, then results on the SPECjvm2008 benchmark suite. Our tests were made on a 8-core, 64-bit Intel(R) Xeon(R) CPU E5440 2.83GHz workstation, with 16GB RAM, running Windows Vista SP1. All tests were performed with our tracer in mode 2 (one trace by thread).

The benchmarks are:

• ArgoUML v0.28, the leading Open Source UML modeling tool including support for all standard UML 1.4 diagrams. 5173 classes are loaded during startup.

• JEdit v4.3pre17, a mature programmer text editor with hundreds of person-years of development behind it, including developing plug-ins development. 2822 classes are loaded during startup.

• Columba v1.0, an email client written in Java, featuring a user-friendly GUI with wizards and internationalization support. The runtime of the startup of this interactive benchmark loads 3646 classes.

• Apache Ant v1.8.2, a Java library and command-line tool to drive processes described in build files as targets and extension points dependent upon each other.

Ant is mainly used to build Java applications. The runtime of the startup of Ant (performing a compilation of itself) loads 1438 classes. We chose to build Ant itself because it does activate the most common features that are used to build a typical java project according to (Andy [START_REF] Andy Zaidman | Automatic identification of key classes in a software system using webmining techniques[END_REF]).

• SPECjvm2008, a benchmark suite that focuses on the performance of the JRE executing a single application; it reflects the performance of the hardware processor and memory subsystem [START_REF] Shiv | SPECjvm2008 performance characterization[END_REF]). The SPECjvm2008 workload mimics a variety of common general purpose application computations [START_REF] Benchmarks | Standard performance evaluation corporation[END_REF]).

An issue with instrumentation is the instrumented program slowdown. In JBInsTrace, three elements contribute to this: the instrumentation process itself, the instrumented bytecode which calls the tracer back, and the recording of the static information and trace to the disk.

Our JBInsTrace tool has two ways of writing the trace files to disk. When memory is constrained, the trace is incrementally written to the disk, at "stop-the-world" garbage collector times. This allows a larger trace to be built, without having to retain it all in memory. Without memory constraints, the entire trace can be kept in memory and dumped to disk at the very end of the execution.

Table 1 presents performance results on ArgoUML, JEdit, Columba and Ant, showing the impact of JBInsTrace on the benchmark runtime. Our goal is to demonstrate that the instrumented software is still fully usable even with our tracer switched on. The first three benchmarks are interactive programs, but to avoid disturbing performance measurements with human interaction and for the sake of repeatability we only traced their startup. The second column in the table shows the benchmark execution time without tracer. The 3rd and 4th columns show the runtime figures with the tracer on.

Column 3 is the execution time, while column 4 is the time to save data to disk (we had enough memory to keep the entire trace in memory and dump it at the end of the runtime). Regarding total execution times with and without the tracer, the slowdown factor is 21.8 for ArgoUML, 4.8 for JEdit, 5.1 for Columba and 42 for Ant. Obviously, there is still work to do to improve the performance of our tracer.

But speed is not our main goal: we first want to acquire a very detailed trace of the execution, including JRE classes, while keeping the tool easy to use and very flexible regarding the information that can be collected at runtime. This of course implies some performance penalty.

However, looking at how the different stages contribute to the slowdown in Table 1 shows that most of the time is spent saving data to disk: the larger the trace saved to disk, the higher the slowdown factor. This explains the difference of performances between these benchmarks, because slowdown factors correspond to the trace sizes magnitude. Note that the size of the trace depends on the size of the scenario and operation performed in the software. The time to write the data to disk depends strongly on hardware. For example, it should be significantly lower with a SSD hard drive.

Although the time spent to write to disk may not be neglected, is it important to remember that our main goal is to maintain a good perceived user experience, keeping the instrumented software fully usable at runtime, even when the tracer runs.

Moreover, in our tests the trace fitted in memory until the end of the tracing, so recording did not bother the user during the execution of the software per se.

For larger software systems, nonetheless, the trace does not fit in memory. In such cases, JBInsTrace has to periodically pause the execution to dump the trace to disk.

When the time spend to write data to disk is factored out, the -perceivedslowdown factors during the execution phase are 4.9 for ArgoUML, 2.3 for JEdit, 2.4 for Columba and 8 for Ant. This shows that the software remains usable with our tracer switched on, which is further confirmed by usability tests performed with actual users.

Figure 6 presents JBInsTrace performance on the famous SPECjvm2008 benchmark suite. These are much smaller, non-interactive and computationally intensive benchmarks. The overhead is computed as operations second and the slowdown factor is computed as operations/second without prof iling operations/second with prof iling . For this experiment the time to write the trace to the disk at the end of the profiling is not accounted. The tests were performed with both HotSpot Client and HotSpot Server JVMs, with little difference in the results.

The overall results with SPECjvm2008 are consistent with the ones obtained with our 5 larger Java applications, although the intensity of computations degrades performance a bit. The sunflow benchmark is a pathological case, not very representative of the actual performance of JBInsTrace, that we have to investigate.

Conclusion and future work

In this paper, we presented our tracing technique and tool, JBInsTrace. JBInsTrace makes it possible to trace what happens at runtime in Java programs at basic block level, on program classes as well as JRE classes using only dynamic bytecode instrumentation.

This complete tracing including JRE classes, at basic block level, generates a very detailed and complete trace of the runtime, representing a huge amount of information to process. Experimental results show that our tracer has a reasonable performance penalty. However pure speed is not our main goal: perceived user experience is more relevant and JBInsTrace provides a good one. Furthermore, JBInsTrace is mainly a research-oriented tool; nonetheless its scalability limits can be found in industrial software [START_REF] Zaidman | Scalability solutions for program comprehension through dynamic analysis[END_REF]).

Our implementation is based on the Java agent service bundled with Java 1.6. As a result, JBInsTrace is very easy to install and ready to use, with no impact on existing Java library files, which makes life easier for users.

Several ways exist to improve our technique and tool. The first direction of our future work will imply performing extensive benchmarking of Java programs with our tracer to gather broader results. The second direction will add new methods in classes when limitations to transform core JRE classes will be lifted. We will especially add a finalize() to every class of the program, to detect and trace object destructions when they occur. Another future direction concerns how to effectively exploit runtime information to have a better understanding of software runtime and to optimize programs.

Figure 1 :

 1 Figure 1: Functioning of our profiler

Figure 2 :

 2 Figure 2: Event number structure

Figure 3 :

 3 Figure 3: Java bytecode instrumentation principle (pseudo-Java code)

Figure 5

 5 presents the UML class diagram of our JBInsTrace tracer. The LockManager class plays a central role in the design of the tracer. It provides a way to differentiate, in the executed bytecodes, those from the program itself and those

Figure 4 :

 4 Figure 4: Sequence diagram illustrating the technique used to prevent trace pollution

Figure 6 :

 6 Figure 6: Performance results on the SpecJVM Benchmark.

Table 1 :

 1 Impact of our JBInsTrace tracer (execution times, in seconds).

	Software	Without	With the tracer
		the tracer	Execution Recording
	ArgoUML	9	44	152
	JEdit	6	14	15
	Columba	7	17	19
	Ant	5	40	170

Available at: http://www.loria.fr/˜casertap/jbinstrace.html

A basic block is a sequence of bytecode instructions which ends by a jump instruction.


Sébastien Beyou

Thomas Corpetti

Sai Gorthi

Etienne Mémin
email: etienne.memin@inria.fr

Fluid flow estimation with multiscale ensemble filters based on motion measurements under location uncertainty

Keywords: Data assimilation, stochastic filter, particle filters, fluid motion estimation

. The data assimilation proposed in this work incorporates measurement brought by an efficient multiscale stochastic formulation of the well-known Lucas-Kanade (LK) estimator. This estimator has the great advantage to provide uncertainties associated to the motion measurements at different scales. The proposed assimilation scheme benefits from this multiscale uncertainty information and enables to enforce a physically plausible dynamical consistency of the estimated motion fields along the image sequence. Experimental evaluations are presented on synthetic and real fluid flow sequences.

Introduction

The analysis of geophysical fluid flows is of the utmost importance in domains such as oceanography, hydrology or meteorology for applications of forecasting, studies on climate changes or for monitoring hazards or events. The forecasting of such flows requires the precise knowledge of an initial condition which may be only accessible through the measurements of the system's state variables such as pressure, temperature, or fluid flow velocity. These data may be provided through dedicated probes or Lagrangian drifters launch in the ocean or in the atmosphere. However, the coverage of such measurements is usually irregular and sometimes very sparse in underdeveloped regions or across oceans. At the opposite orbital or geostationary satellites provide a huge amount of image data, with a still increasing spatial and temporal resolution. Compared to in situ measurements (i.e. measure with local probes located inside the flow), satellite images provide a much more denser observation field. However they unfortunately offer only an indirect access to the physical quantities of interest, and give rise consequently to difficult inverse problems to estimate characteristic features of the flow such as velocity fields or vorticity maps.

These kinematical information can be estimated from image sequences through motion estimation techniques. Motion estimation is an old problem in computer vision and a huge number of techniques has been proposed to estimate the motion on the image plane of a 3D scene. The best state-of-the-art approaches perform efficiently for the recovery of rigid scene motions [START_REF] Papenberg | Highly accurate optic flow computation with theoretically justified warping[END_REF][START_REF] Sun | Layered segmentation and optical flow estimation over time[END_REF]. They are generally built on strong photometric and motion hypothesis which prevent them to be sufficiently accurate for deformation metrology. Several motion estimators dedicated to the measurement of specific phenomenon such as fluid flows have been proposed in the literature (see [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF] for a detailed overview). These estimators differ mainly on the smoothness prior they are handling: first order penalization [START_REF] Ruhnau | Variational optical flow estimation for particle image velocimetry[END_REF], second order div-curl regularization [START_REF] Corpetti | Fluid experimental flow estimation based on an optical-flow scheme[END_REF][START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF], data-dependent [START_REF] Becciu | Cardiac motion estimation using covariant derivatives and Helmholtz decomposition[END_REF][START_REF] Duits | A variational approach to cardiac motion estimation based on covariant derivatives and multi-scale helmholtz decomposition[END_REF][START_REF] Nir | Over-parameterized variational optical flow[END_REF] or power law auto-similarity principles [START_REF] Heas | Bayesian selection of scaling laws for motion modeling in images[END_REF][START_REF] Heas | Power laws and inverse motion modeling: application to turbulence measurements from satellite images[END_REF]. These methods provide accurate instantaneous displacements as they generally implement additional constraints imposed by the physics, as in [START_REF] Corpetti | Fluid experimental flow estimation based on an optical-flow scheme[END_REF][START_REF] Florack | The intrinsic structure of optic flow incorporating measurement duality[END_REF] and most of them are embedded into a multiscale formalism that enables capturing efficiently the large scales deformations [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF]. However, those techniques generally still exhibit difficulties for mid to small scales measurements. Such artifacts may reveal particularly pronounced in regions with poor photometric contrasts where the smoothing prior takes the lead in the solution elaboration. In those regions the data come into play only at the boundary. Bad estimations or instabilities on the boundary vicinity are immediately echoed on inside such regions. For large regions this may reveal problematic and constitutes a potential factor of instabilities along time of the estimates.

Dynamical consistency of the velocity measurements can be enforced by embedding the estimation problem within an image based assimilation process. Variational assimilations of image information have been recently considered for the estimation of fluid motion fields [START_REF] Corpetti | Pressure image asimilation for atmospheric motion estimation[END_REF][START_REF] Papadakis | An optimal control technique for fluid motion estimation[END_REF]. Those optimal control methods, though efficient, constitute batch methods, which requires forward and backward integrations of the dynamical system and the adjoint of the tangent linear dynamics respectively. The latter relies implicitly on a linearization of the dynamics and is adapted in practice for short time horizon. The constitution of this adjoint dynamics may turn out quite tedious in practice for complex dynamical models.

Stochastic filters are also well known alternative techniques for data assimilation. Opposite to the variational data assimilation framework, stochastic filtering has the great advantage to couple noisy data and a stochastic dynamics incorporating the unavoidable uncertainties we have on the system evolution. Such filters are also generally set up in a recursive way and through a Markovian property of the dynamics, they are at least theoretically less dependent on the initial condition. This relative independency with respect to the initial condition and their recursive structure along time are two important advantages. Obviously those filters require a precise modeling of the conditional distribution of the data given the true state, which is not always an easy task. Besides, the constitution of sound stochastic dynamics is in itself a difficult research problem for which no general theory exists. This filtering problem is in practice implemented either through particle filters or ensemble Kalman filters. The latter performs efficiently for linear systems but does not converge toward the true filtering distribution for nonlinear likelihood or dynamics [START_REF] Le Gland | Large sample asymptotics for the ensemble Kalman filter[END_REF]. As for the former, it has the advantage to deal with complex dynamics but is in its basic version limited to state spaces of moderate dimension [START_REF] Snyder | Obstacles to high-dimensional particle filtering[END_REF], typically inferior to 100 and in any case three or four orders of magnitude lower than the dimension we are facing (corresponding to the image dimension e.g. ∼ [10 5 , 10 6 ]). Recently, a data assimilation procedure with the attempt of coupling the advantages of both filters has been proposed in embedding an Ensemble Kalman filter (EnKF) [START_REF] Evensen | The ensemble Kalman filter, theoretical formulation and practical implementation[END_REF] into the particle filter (PF) framework. The approach is referred as Weighted Ensemble Kalman filter (WEnKF) [START_REF] Papadakis | Data assimilation with the weighted ensemble Kalman filter[END_REF]. The objective of this work consists to specify such a procedure from local noisy velocity measurements and their uncertainties and to define a multiscale extension of the filter.

To that end, we organize the paper as follows: Section 2 introduces the stochastic formulation of a local motion estimator based on the Lucas and Kanade technique [START_REF] Lucas | An iterative image registration technique with an application to stereovision[END_REF] and referred as SLK in the rest of the paper. This modified version supplies the motion measurement at specific scale/resolution, along with their corresponding uncertainties. In Section 3, we present the Weighted Ensemble Kalman Filter (WEnKF) based fluid flow assimilation, relying on Navier-Stokes state dynamics and a linear observation model that uses SLK observations and uncertainties. Our novel multi-level WEnKF scheme, which is based on the same dynamical model but assimilates the observations at different scales, is proposed in section 4. Some experimental validation and comparisons with the state-of-art fluid flow estimators are brought out in section 5 and finally, section 6 concludes the paper with a discussion to future directions.

Stochastic Lucas-Kanade Estimator

Many possibilities are available to derive motion estimates from image data. This task, commonly named "optical flow", has been widely studied over the past three decades and huge amount of techniques are available. The reader can refer to [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF] for a recent review on estimators dedicated to fluid flows and more generally to [START_REF] Baker | A Database and Evaluation Methodology for Optical Flow[END_REF][START_REF] Barron | Performance of optical flow techniques[END_REF][START_REF] Galvin | Recovering motion fields: an analysis of eight optical flow algorithms[END_REF] for presentation and comparative performance evaluations of some state-of-the-art techniques in computer vision. Interested readers may also refer to [START_REF] Becciu | Cardiac motion estimation using covariant derivatives and Helmholtz decomposition[END_REF][START_REF] Duits | A variational approach to cardiac motion estimation based on covariant derivatives and multi-scale helmholtz decomposition[END_REF][START_REF] Nir | Over-parameterized variational optical flow[END_REF], which provide recent trends for data based adaptive deformation estimation schemes.

Apart from some very specific applications where the observation process can be modeled in a precise way, almost all existing approaches are based on the well-known "opticalflow constraint equation" (OFCE), assuming a global conservation of the image luminance:

d f (x, t) d t = ∂ f (x, t) ∂ t + v(x, t) • ∇ f (x, t) = 0, (2.1)
where x = (x, y) corresponds to the spatial coordinates, ∇ is the gradient operator in the x and y directions and f denotes the luminance function.

Lucas-Kanade estimator

From relation (2.1), the seminal works of Horn and Schunck [START_REF] Horn | Determining optical flow[END_REF] and Lucas and Kanade [START_REF] Lucas | An iterative image registration technique with an application to stereovision[END_REF] have been proposed. In the former method, the solution is obtained by minimizing a global energy functional. This functional combines the OFCE with a first-order spatial regularizer that promotes smooth motion fields over the whole image domain. The Lucas and Kanade technique is at the opposite a local approach that assumes a locally constant velocity on neighborhoods centered at each location x of the image domain. At each point, the velocity is estimated through a local least squares formulation:

(x) = ∂ f (x, t) ∂ t + v(x, t) • ∇ f (x, t) 2 and v = min v=(u,v) T Ω [g λ (•) * (•)](x)dx, (2.2)
where Ω is the image domain, g λ is an isotropic Gaussian windowing function of standard deviation λ in which the velocity v is assumed to be homogeneous. Solving the Euler-Lagrange equations associated to cost function (2.2), one gets (see appendix A for details):

v = -g λ * f 2 x f x f y f x f y f 2 y -1 g λ * f x f t f y f t , (2.3)
where

f • = ∂ f /∂ •.
To guarantee a good conditioning of the previous matrix to invert, the spatial gradients must not vanish or be identical. The gaussian smoothing of the OFCE aims in fact at alleviating such a bad conditionning by capturing the spatial information at a scale related to λ. It is important to outline that the estimated velocity is hence intrinsically related to this scale λ.

From Horn & Schunck and Lucas & Kanade estimators, a huge number of methods based on diverse variations of these approaches has been proposed in the literature. Among them, the techniques focused first on the design of new regularization terms (able for instance to deal with occlusions, discontinuities or relying on physical grounds [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF]) and second on the application of advanced minimization strategies. Surprisingly, apart for some specific applications devoted to some particular types of imagery (fluid, biology, infrared imagery, tomography, IRM, ...), only very few authors have worked on a generic alternative to the classical brightness consistency assumption, and this despite to the fact it plays a crucial role in the motion estimation process.

In this paper, instead of relying on the OFCE to derive local motion measurements, we propose to formulate the motion estimation upon a stochastic version of the image luminance variation and to embed it in a local least squares formulation similar to the Lucas & Kanade technique. This leads naturally to a continuous multiresolution formulation and enables not only to extract the motion fields at different resolutions but supplies uncertainties of those estimates as well. This is presented in the next sections.

Luminance variation with uncertainties

We first write the image luminance as the function of a stochastic process related to the position of image points. If one assumes, as illustrated in figure 1, that the velocity v to be estimated transports from time t -1 to time t the grid of 2D points X t-1 ∈ 2m to X t ∈ 2m up to an isotropic Brownian motion, we can write:

X t = X t-1 + t t-1 v(s)ds + t t-1 Σ 1/2 (s)dB s , (2.4)
where B s = (B 1 s , ..., B m s ) T is a multidimensional standard Brownian motion of 2m , Σ a (2m × 2m) covariance matrix-valued function and X t , X t-1 represents the grid positions The initial grid at time t -1 (white) is transported by the velocity field v to reach the configuration at time t represented (black), up to some uncertainties (dashed lines).

at time t -1 and t. Let us note the second right-hand integral has to be interpreted in a stochastic sense as Brownian motion is almost nowhere differentiable. In this work we shall rely on the Itô interpretation of such integral. The luminance function f usually defined on spatial points (x, y) at time t is now defined on the grid as a map from 2m × + into m and is assumed to be C 1,2 ( 2m , + ). Its differential is obtained following the differentiation rules of stochastic calculus (the so called Itô formulae) that gives the expression of the differential of any continuous deterministic function of an Itô diffusion of the form 2.4 (see [START_REF]Stochastic differential equations[END_REF] for an introduction to stochastic calculus):

df(X t , t) = ∂ f ∂ t d t + i=(1,2) ∂ f(X t , t) ∂ x i d X i t + 1 2 (i, j)∈(1,2)×(1,2) ∂ 2 f(X t , t) ∂ x i ∂ x j d < X i t ,X j t > . (2.5)
The term < X i t ,X j t > denotes the joint quadratic variations of X i and X j defined as the limit in probability over a partition {t 1 , . . . , t n } of [0, t] with t 1 < t 2 < • • • < t n , and where denoting δt i = t it i-1 , the partition spacing defined as

|δ n | = max i δt i is such that |δ n | → 0 when n → ∞: < X i t ,X j t >= lim |δ n |→0 t k ≤t (X i t k-1 -X i t k )(X j t k-1 -X j t k ), for n → ∞.
For Brownian motion, this quantity can be evaluated through the expectation of the righthand sum, an almost sure convergence is then obtained from the strong law of large numbers (which implies convergence in probability). We have then the following identities:

< B i ,B j >= δ i j t a.s. < h(t), h(t) >=< h(t), d B i >=< B j , h(t) >= 0 a.s. , (2.6)
where δ i j is the Kronecker symbol (δ i j = 1 if i = j, δ i j = 0 otherwise), and h(t) is a deterministic function. Compared to classical differential calculus, new terms related to the Brownian random terms have been introduced in this stochastic formulation. In this work the stochastic part of (2.4) is defined as an isotropic uncertainty variance map σ(X t , t)

: 2m × + → m Σ 1/2 (X t , t) dB t = diag( σ(X t , t)) ⊗ 2 dB t , (2.7)
where 2 is the (2 × 2) identity matrix, and ⊗ denotes the Kronecker product. Applying Itô formula (2.5) to the isotropic uncertainty model yields a luminance variation defined as:

df(X t , t) = ∂ f(X t , t) ∂ t + ∇f(X t , t) • v(X t , t) + 1 2 σ(X t , t)∆f(X t , t) d t + σ(X t , t)∇f(X t , t) • dB t . (2.8)
It is straightforward to remark that the standard brightness consistency assumption is obtained from (2.8) using zero uncertainties (σ = 0). This stochastic formulation enables thus to introduce a softer constraint that takes into account an inherent motion measurement due to the discrete grid representation for instance. In this simple model the set of grid points locations at time t are thus conditionally independent with respect to their position at time t -1.

From relation (2.8), the conservation of the image luminance can be quite naturally expressed from the conditional expectation df(X t , t)|X t-1 between t -1 and t. To compute this term, we exploit the fact (as shown in appendix B) that the expectation of any function Ψ(X t , t) of a stochastic process dX t (as in (2.4)) knowing the grid X t-1 reads:

(Ψ(X t , t)|X t-1 ) = Ψ(•, t) * g σ (•) (X t-1 ), (2.9)
where g σ (X t ) ∼ (X t , Σ) is a multidimensional Gaussian centered in X t with covariance Σ(X t , t) † and Ψ(X t , t) = Ψ(X t + v, t). Relation (2.9) indicates that the expectation of a function Ψ(X t , t) knowing the location X t-1 up to a Brownian uncertainty of variance σ is obtained by a convolution of Ψ(X t-1 + v, t) with a Gaussian kernel of mean X t and covariance Σ.

Assuming σ known, our new luminance variation model (df(X t , t)|X t-1 ) is hence defined as (removing the indexes (X t , t) for sake of clarity):

df(X t , t)|X t-1 = g σ * ∇f • v + ∂ f ∂ t + 1 2 σ∆f . (2.10)
Instead of using the OFCE as the basis constraint in the motion estimation process, we prefer to build our least squares estimator upon relation (2.10), which constitutes an OFCE model under location uncertainty. In the next section, we present a way to estimate such uncertainties σ(X t , t).

Uncertainty estimation

Assuming an observed motion field v obs that transports the luminance is accessible (we will describe in the next section an incremental framework involving such a motion), it is possible to estimate the uncertainty σ(X t , t) by computing the quadratic variation of the luminance function df between t -1 and t. Using the identities in (2.6), this yields:

d f(X t , t), f(X t , t) = σ(X t , t) ∇f(X t , t)) 2 .
(2.11)

This quadratic variation can also be approximated from the luminance f by:

d f(X t , t),f(X t , t) ≈ f(X t , t) -f(X t-1 , t -1) 2 .
(2.12)

As convergence in probability implies convergence in distribution, the conditional expectation of both previous terms should be identical, and one can estimate the variance by:

σ(X t , t) = f(X t , t) -f(X t-1 , t -1) 2 ∇f(X t , t)) 2 . (2.13)
The expectation in the numerator and denominator are then computed at the displaced point X t-1 + v obs (X t-1 ) through the convolution of variance Σ(X t-1 , t -1). A recursive estimation process is thus emerging from equation (2.13). Before entering into the details of the complete framework in section 2.6, let us now rewrite the Lucas-Kanade estimation associated to our OFCE model under location uncertainty (2.10).

Stochastic Lucas-Kanade estimator

A local motion estimator relying on the new observation model defined in (2.10) can be now expressed as the following minimization problem (dropping the coordinates (X t , t)):

min v g λ * g σ * ∇f • v + ∂ f ∂ t + 1 2 σ∆f 2 .
(2.14)

In a similar way to the standard Lucas and Kanade estimator (see appendix A) this yields the following relation, for the velocity estimate of a given grid point:

g λ * g σ * f 2 x f x f y f x f y f 2 y v = -g λ * g σ * ( 1
2 σ∆f+f t ) f x f y . (2.15)
In this model the Gaussian windowing function g λ of mean X t and covariance λ 2m can be interpreted as the distribution of a new isotropic constant uncertainty term related to the grid resolution and independent of the motion uncertainties depending on the image data.

Let us note otherwise that the data model (2.14) introduces a Laplacian diffusion with a minus sign (antidiffusion effect). This numerically unstable effect is however counterbalanced by the Gaussian smoothing g σ with a variance that is twice the antidiffusion coefficient. The supplementary filtering due to the Gaussian windowing function g λ introduced to counter the illposed nature of the brightness consistancy assumption, further stabilizes the process. In practice, the choice of the convolution kernel standard deviation λ is crucial: a large value of λ will remove all details whereas a small value is likely to be unstable and may lead in an extreme case to an illposed problem. In addition, such motion estimation procedure, based on a linearized version of the displaced frame difference, leads to inaccurate measurements of large displacements when the linear assumption of the brightness consistency breaks (high photometric gradients and/or large displacements).

To prevent such limitations many authors have proposed to embed such estimation procedures within a pyramidal setup. However, a pyramidal representation requires Gaussian filtering, sampling and interpolation of the input data which introduce potential artifacts that may spoil the estimation. Interpreting the windowing function convolution as associated to the computation of the expectation of a grid uncertainty random term will allow us defining an original continuous multiresolution framework. This is presented in the next paragraph.

Multiresolution version

A multiresolution scheme consists in redefining the problem on a grid X ℓ which can be viewed as a coarse representation of the initial grid X 0 = X with a Brownian isotropic uncertainty of constant variance λ ℓ : from the initial pixel grid X 0 = X, which corresponds to the whole image domain, the problem is redefined on a grid X ℓ associated to a coarse representation of X. This reads

X ℓ = X 0 + λ ℓ 2 dB t .
(2.16)

The motion v(X ℓ t-1 , t -1) on this grid should minimize the expectation E(df 2 (X ℓ t , t)|X 0 ) which is equivalent (see appendix B) to a convolution of df 2 (•, t) with the isotropic gaussian (0, λ ℓ ). Therefore, one gets exactly the system in (2.14) which is solved locally by inverting the system in (2.15). The main advantage of such a formulation of the multiresolution setup is to naturally get rid of the use of a pyramidal image representation. Instead of dealing with successive decimations of factor 2 of the initial image to fix the different multiresolution levels, the evolutions of the levels ℓ are much flexible here. This framework may be interpreted as a scale space representation of the extended brightness consistency assumption [START_REF] Lindeberg | Scale-space theory: A basic tool for analysing structures at different scales[END_REF] stemming in our case rigorously from an uncertainty analysis of the point location. Let us note the heat equation, which correspond to the backward Kolmogorov equations associated to (2.16), constitutes formally the connection point between the traditional scale space representation and this multiresolution setup.

From the Stochastic Lucas-Kanade normal equations (2.15), the estimation of uncertainties from σ of relation (2.13) and the above multiresolution framework, one can define a complete incremental framework presented in the next section.

Incremental framework

The local estimator proposed here comprises a specific multiresolution scheme where at each level, an incremental framework is defined to estimate as accurately as possible the uncertainties. The overall scheme is figured in algorithm 2.1. Algorithm 2.1 Stochastic Lucas-Kanade: incremental framework 1. Initializations :

• Fix an initial resolution level ℓ = L • Define f(X t-1 , t) := f(X t-1 , t) ; v = 0;
2. Estimation for the level ℓ (a) Initializations : This framework is a natural and simple implementation of a local motion estimation technique using the proposed isotropic uncertainty model for the evolution of the luminance. The values σ provide us a spatial distribution of the motion estimate uncertainties. These uncertainties constitute a very usefull piece of information when coping with the filtering of noisy motion measurements. Let us now turn to the presentation of the stochastic filter considered in this work.

• n = 1; v 0 = 0; • Fix an initial σ 0 (b) Estimate σ n by relation (2.13) (c) Find v n by local inversions of the system (2.15) (d) Update motion field : v := v + v n (e) Warp the image f(X t , t) : f(X t-1 ) = f(X t-1 + v, t) (f) n := n + 1

Monte Carlo implementation of stochastic filtering with the Weighted-Ensemble Kalman Filter

In this section we briefly review the main principles driving the construction of the Weighted Ensemble Kalman filter, proposed in [START_REF] Papadakis | Data assimilation with the weighted ensemble Kalman filter[END_REF], and discuss its advantages and limitations in the context of fluid flow analysis. This technique is a particle implementation of a nonlinear stochastic filtering problem build upon an ensemble Kalman update stage. Let us first introduce basic elements of stochastic filtering.

Stochastic filtering, filtering distribution

Stochastic filters aim at estimating the posterior probability distribution p(x 0:k |y 1:k ) of a state variable trajectory x 0:k starting from an initial state x 0 ∈ n up to the current time k given a complete measurements trajectory y 1:k . The state variable trajectory is obtained through the integration of a dynamical system:

x t = M(x t-δt ) + η t , (3.1)
where M denotes a deterministic linear/nonlinear dynamical operator, corresponding to a discrete representation (through numerical integration with time step δt) of a physical conservation law describing the state evolution. Parameter η t is usually a white Gaussian noise of covariance Q δt , that accounts for the uncertainties in the deterministic state model. However, as the true initial state is unknown, observation y k ∈ m of the state occurring at discrete instants are assumed to be available. These observations and the state variable are linked through a measurement equation:

y k = H(x k ) + γ k , (3.2)
where γ k , the observation noise, is a white Gaussian noise with covariance matrix R, and operator H stands for the linear/nonlinear mapping from the state variable space to the observation space. We note that the (integration) time step used for the state variable dynamics δt is usually much smaller (about 10-100 times), than the latency δk between two subsequent measurements. A sequence of measurements or observations from time 1 to k will be denoted by a set of vectors of dimension m as: y 1:k = {y i , i = 1, . . . , k} where the latency between two successive measurements is arbitrarily set to δk = 1.

A recursive expression of the filtering distribution p(x 0:k |y 1:k ), describing the distribution of the unknown hidden Markov process conditioned upon the whole set of past observations y 1:k , can be obtained from Bayes' law and the assumption that the measurements depend only on the current state (e.g. p(y k |x 0:k , y 1:k-1 ) = p(y k |x k )):

p(x 0:k |y 1:k )p(y k |y 1:k-1 ) = p(x 0:k , y k |y 1:k-1 ) = p(y k |x 0:k , y 1:k-1 )p(x 0:k |y 1:k-1 ) = p(y k |x k )p(x k |x 0:k-1 , y 1:k-1 )p(x 0:k-1 |y 1:k-1 )
from Markoviannity of the dynamic system, we finally get the sought recursive expression:

p(x 0:k |y 1:k ) = p(x 0:k-1 |y 1:k-1 ) p(y k |x k )p(x k |x k-1 ) p(y k |y 1:k-1 ) . ( 3.3)
For linear models, this recursive equation can be solved through the Kalman filter. This is described in the next section.

Linear Gaussian models and the Kalman Filter

For a Gaussian initial distribution, additive Gaussian noises, linear dynamics and linear measurement operator, denoted by M and H respectively, the distribution p(x k |y 1:k ) remains a Gaussian distribution whose first and second moment, x a k = (x k |y 1:k ) and P a k = ((xx a k )(xx a k ) T |y 1:k ), can be explicitly computed from the well known recursive Kalman equations [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]:

x f k = Mx a k-1 , P f k = MP a k-1 M T + Q k , (3.4)
and

K k = P f k H T (HP f k H T + R) -1 , (3.5)
x a k = x f k + K k (y k -Hx f k ), (3.6)
P a k = ( -K k H)P f k , (3.7)
here superscripts f and a on state variable expectation and covariance denote the respective quantities before and after analysis (update) at time k, respectively. The prediction or forecast step (3.4) brings forward the first two moments of the state vector, from its previous time step k -1, through the dynamical model parameters, while the analysis or the correction step (3.7) provides the first two moments of the state characterizing the Gaussian filtering distribution at time k. The matrix K k is referred to as the Kalman gain matrix and defines the weighting between the forecast state and the innovation brought by the observation : y k -Hx f k . In the case of nonlinear models, it is possible to approximate the posterior distribution with a set of particles. This is presented in the next section.

Particle implementation of the nonlinear filtering

For nonlinear dynamics or nonlinear measurement equation, a direct sampling from the filtering distribution is impossible since it would require the complete knowledge of the filtering distribution -which is in the general case a non Gaussian multimodal distribution -at a previous time.

Particle filtering techniques introduce a discrete approximation of the sought density as a sum of N weighted Diracs:

p(x 0:k |y 1:k ) ≈ N i=1 w (i) k δ x 0:k x 0:k , (3.8)
centered on hypothesized locations of the state space sampled from a proposal distribution, π, (also called the importance distribution) approximating the true filtering distribution p(x 0:k |y 1:k ). Each sample is then weighted by w

(i)
k , accounting for the ratio between the two distributions. A great variety of importance functions can be chosen (with the only restriction that its support contains the filtering distribution one). Obviously the closer it is to the targeted filtering distribution the better the results. Under a weak hypotheses assuming a recursive factorization form of the importance function, the importance ratio can be recursively defined as:

w (i) k ∝ w (i) k-1 p(y k |x (i) k )p(x (i) k |x (i) k-1 ) π(x (i) k |x (i) 0:k-1 , y 1:k ) .
(3.9) By propagating the particles from time k-1 through the proposal density π(x

(i) k |x (i) 0:k-1
, y 1:k ), and by weighting the sampled states with w

(i)
k , a sampling of the filtering law is obtained. When the proposal distribution π is set to the prior (i.e. π(x

(i) k |x (i) 0:k-1 , y 1:k ) = p(x (i) k |x (i) k-1
), the weights updating rule (3.9) simplifies to the data likelihood p(y k |x (i) k ). This particular instance of the particle filter is called the Bootstrap filter or sequential importance resampling (SIR) filter [START_REF] Gordon | Novel approach to non-linear/non-Gaussian Bayesian state estimation[END_REF]. Due to its simplicity it is the most commonly used particle filter. Nevertheless, such an importance function does not take into account the current observation and depends only weakly on the past data through the filtering distribution estimated at the previous instant. High dimensional probability distribution spaces being excruciatingly difficult to sample, it is very important to devise an importance function that enables focusing on the most meaningful areas of the state space. To that end it is essential to consider proposal distributions that take into account more significantly the past and current measurements. Along this idea, the weighted ensemble Kalman filter defines the proposal distribution from the sampling mechanisms of ensemble Kalman filtering techniques.

Ensemble Kalman Filtering

The Ensemble Kalman filter [START_REF] Evensen | Sequential data assimilation with a non linear quasi-geostrophic model using Monte Carlo methods to forecast error statistics[END_REF] can be interpreted as a Monte Carlo implementation of the Kalman filter recursion for the propagation of the two first moments. The Ensemble filter relies hence intrinsically on a Gaussian approximation of the filtering distribution.

More precisely, let us assume that we have sampled N members from the initial filtering distribution p(x 0 |y 0 ), denoted by x (i) 0 , i = 1, ..., N . Propagating these samples, iteratively, through the Kalman prediction and correction steps, provides us the Gaussian approximations of the prediction and filtering distributions.

The prediction step consists in propagating the ensemble members x a,(i)

k-1 and their associated uncertainties (noise) through the state dynamics in order to obtain a predicted particles or forecast ensemble as:

x f ,(i) k = k-δt t=k-1 M(x f ,(i) t ) + η (i) t+δt , x f ,(i) k-1 = x a,(i) k-1 .
(3.10)

From this, the empirical mean, x f k , of the forecast ensemble and the corresponding empirical forecast covariance matrix P , i = 1, ..., N } defined as:

x a,(i) k = x f ,(i) k + K e k y k + γ (i) k -Hx f ,(i) k . (3.11)
Here, we note that, in the Kalman gain or in the update stage, computation of the high dimensional covariance matrix or pseudoinverse of the n × n covariance term, (HP -1 , are never explicitly computed nor stored. Rather, Kalman gain and update are efficiently implemented by defining and employing matrices with ensemble of perturbations. In most of the geophysical applications, the size n of the state vector related usually to temperature, pressure or velocity fields is of much higher dimension than the number of samples N used in EnKF. i.e., n >> N . Thus, handling the perturbation matrices (instead of the actual corresponding covariance matrices) approximately brings down the number of operations from O(n 2 ) to O(nN ). The inverse needed in the Kalman gain can be efficiently computed through the singular value decomposition of a n × N matrix [START_REF] Evensen | The ensemble Kalman filter, theoretical formulation and practical implementation[END_REF].

f e k H T + R)

Weighted EnKF

Starting from the descriptions of the previous section, a hybrid filtering procedure that takes advantage of both the particle filter and the EnKF can be devised. We briefly describe the approach proposed in [START_REF] Papadakis | Data assimilation with the weighted ensemble Kalman filter[END_REF].

The importance sampling principle indicates that a wide range of proposal distributions can be considered. We will experimentally show that a proposal distribution defined by the EnKF procedure constitutes an efficient proposal mechanism for particle filter techniques in high dimensional spaces.

Relying on the usual assumption of the EnKF (i.e. considering the dynamics as a discrete Gaussian system), the conditional distribution p(x k |x (i) k , y k ) can be approached by a Gaussian distribution of respective mean and covariance [START_REF] Papadakis | Data assimilation with the weighted ensemble Kalman filter[END_REF]:

µ (i) k = ( -K e k H) k-∆t t=k-1 M (x f ,(i) t ) + K e k y k , Σ e k = ( -K e k H)P f e k . ( 3

.12)

This distribution provides us a natural expression for the proposal distribution. In order to make the estimation of the filtering distribution exact (up to the sampling), each member of the ensemble must be weighted at each instant, k, with appropriate weights, w

(i)
k , defined from (3.9). For high dimensional systems represented on the basis of a very small number of particles, the weights simplify as [START_REF] Papadakis | Data assimilation with the weighted ensemble Kalman filter[END_REF]:

w (i) k ∝ p(y o k |x (i) k ),

Multiscale WEnKF assimilation of SLK observations

In what follows, we first detail the dynamical model, the observation model and the strategy we adopt to incorporate the uncertainties supplied by the SLK estimator. In a second step, the section 4.2 presents its multiscale extension to extract consistent fluid motion estimates.

Models and uncertainties

The objective being to estimate velocities related to 2D incompressible fluid flows, we prefer to depict the velocity v = [u, v] T by the scalar ξ = ∂ v/∂ x -∂ u/∂ y = v xu y that represents the vorticity. For such incompressible flows, the velocity field can easily be recovered from its vorticity using the Biot-Savart kernel: v = ∇ ⊥ G * ξ where G denotes the Green kernel (G = 1 2π ln(|x|)) associated to the Laplacian operator.

Dynamical model:

As in this work we considered only 2D incompressible fluid flows, we will rely for the dynamics on the vorticity-velocity formulation of the Navier-Stokes equation with a stochastic forcing function:

dξ = -∇ξ • vd t + ν∆ξd t + ηd B t , (4.1)
where ν is the kinematic viscosity and ηd B is a random forcing term (introduced later in this section). The numerical simulation of this dynamical model is detailed in [START_REF] Papadakis | Data assimilation with the weighted ensemble Kalman filter[END_REF]. We will note ξ k ∈ |Ω| the finite dimentional state vector describing at instant k the vorticity map over the spatial domain Ω and v(ξ k ) will denote the corresponding velocity field (where an explicit dependance on the vorticity is emphasized to indicate the nonlinear nature of dynamics (4.1) ).

Observation model:

The measurements on which we will rely on are set directly as the curl map (i.e. vorticity) of the SLK velocity estimates (cf section 2). Assuming the observation is a corrupted version of the true vorticity map (state), we define the observation model as:

ξ k = ξ k + γ k , (4.2)
where γ k is a Gaussian random field whose variance is fixed to the spatially varying uncertainties associated to the measurements. The constitution of the random fields associated to the dynamics and measurement noises are detailed in the next sections.

Noise and uncertainties:

To simulate the random forcing term dB in the dynamics (4.1) and the random field of the observation model (4.2), homogeneous Gaussian fields, correlated in space but uncorrelated in time, are used. Their covariance have a general form given by:

Q(r, τ) = [dB(x, t)dB T (x + r, t + τ)] = g λ (r)d tδ(τ), (4.3)
where g λ (r) describes the spatial correlation structure with cutoff parameter λ beyond which two points are uncorrelated and δ denotes the Dirac distribution (indicating the random fields are temporally uncorrelated). These random fields are in practice sampled in the Fourier domain or on wavelet basis [START_REF] Elliot | A Fourier-wavelet Monte Carlo method for fractal random fields[END_REF][START_REF] Majda | Simplified models for turbulent diffusion: Theory, numerical modelling, and physical phenomena[END_REF].

Dynamics noise More precisely, for the noise associated to the dynamics, we consider a divergence free random field that mimics real turbulent flows. This is done by exhibiting within an inertial range of scales an energy spectrum defined as a power law model

x ζ 2
with an exponent 0 < ζ < 2. The covariance tensor associated to this random field, usually refered in the litterature as Kraichnan model [START_REF] Kraichnan | Small-scale structure of a scalar field convected by turbulence[END_REF], can be written more conveniently in terms of its spectral representation:

(ηd B t (x)ηdB s (y) T ) = δ(t -s) d t (2π) d d e ik(x-y) Q(k)dk, (4.4)
where the spectral correlation tensor in dimension d is defined as:

Q(k) i j = C 2 ζ |k| -ζ-d (δ i j - k i k j |k| 2 )( ψ κ γ ) 2 .
In the previous relation, C ζ corresponds to a dissipative constant fixed on dimensional considerations and function ψ γ κ defines a band-pass cut-off function, allowing to keep only a given range of scales, corresponding to the so-called inertial scale range defined between the short dissipative scale ℓ D = 1/κ and the large integral scale L = 1/γ at which the forcing takes place. The variance of this homogeneous random field is

C 2 ζ d t (2π) d d |k| -ζ-d (δ i j - k i k j |k| 2 )( ψ γ κ ) 2 dk, (4.5)
which is constant. For a pass-band filter ψ κ γ defined as a box filter: 1I

[γκ] (k) = k if k ∈ [γ, κ], 1I [γκ] (k) = 0 otherwise, this is given by d t C ζ (2π) d d -1 d κ γ k r = k |k| -ζ-d dk d k r δ i j = d t C ζ (2π) d d -1 d 2π d/2 Γ( d 2 ) ζ -1 (L -ℓ D )δ i j .
The constant C ζ has a dimension of ℓ 1-ξ 2 /t which follows from dimensional analysis (the energy transfer rate ε = (

u 2 t ) ∼ ℓ 2 t 3 and ε ∼ C 2 ζ ℓ ξ t ).
The value of the exponent can be fixed either from theoretical considerations or from experimental values.

Measurement noise

As for the observation noise involved in relation(4.2), it is defined from the the final uncertainty supplied by the local motion estimator of section 2 (eq. 2.13). The random field is specified as

γ k (x) = σ k (x)η k (x) where η k is a Gaussian random field with an exponential covariance g λ (x -x ′ ) = exp -(x-x ′ ) 2 λ 2
, and σ k (x) is defined by the uncertainties associated to the motion measurements at each point x. To mitigate the effect of outliers, a Gaussian smoothing with a low variance value is considered. This smoothing avoids the apparition of sharp discontinuities on the uncertainty map. The additive uncertainty γ k is hence a Gaussian random field with non-stationary variance

R(x, x) = σ 2 k (x), and covariance R(x, x ′ ) = σ k (x)σ k (x ′ )g λ (x -x ′ ).
Let us now turn to the description of the WEnKF implementation corresponding to our application.

WEnKF implementation:

The knowledge of the dynamics and observation models together with the corresponding uncertainties we have previously set up constitute all the ingredients needed for the implementation of the WEnKF. The following steps summarize the successive stages of the method:

• At k = 0, the ensemble of states {ξ a,(i) 0 , i = 1, ..., N } are initialized with noisy versions of the SLK vorticity map obtained from the two first images of the sequence using the technique of section 2;

• At the current time k = 0, the ensemble members ξ a,(i) k-1 obtained at the previous measurement instant are propagated through the stochastic state dynamics (4.1) to generate the forecast ensemble members ξ f ,(i) k .

• The EnKF update is then performed with the new observation in order to sample the proposal distribution. The importance sampling weighting (relation (3.13)) based on the likelihood and a resampling process of the particles with respect to those weights are performed. The empirical mean of the analysis ensemble provides the vorticity estimate at time k.

• The corresponding velocity field is finally obtained from the Biot-Savart law.

Although this direct WEnKF filtering of the SLK vorticity maps does provide good results as we shall see it, the estimation may fail for long range velocities. To overcome this limitation and to further improve the performance of the WEnKF, we propose in the next section a multiscale extension of WEnKF.

Multiscale SLK-WEnKF filtering

The idea of multiscale WEnKF consists to provide an improved proposal distribution from velocity measurements at different scales. The update step operates iteratively in an incremental coarse-to-fine way by introducing motion measurements obtained at different scales through the Gaussian smoothing parameter λ ℓ in (2.15). More precisely, at a given scale ℓ ∈ [ℓ 0 , L] the proposal ensemble is build for each ensemble members from successive analysis steps as follow (in the following for sake of notational convenience we drop the member superscript (i), nevertheless apart from mean vectors all the vector are meant to depend on a given member) :

ξ a,ℓ k = ξ f ,ℓ k + K ℓ k ξ ℓ k + γ ℓ k -Hξ f ,ℓ k , (4.6) ξ f ,ℓ k = ξ f ,ℓ+1 k -ξ a,ℓ+1 k . ( 4.7)
Here, ξ

f ℓ k (X k-1 ) = f(X k-1 + L j=ℓ v(ξ a, j+1 k ), k), ( 4.8)
and image f(X k-1 , k -1). The SLK estimator is run for Gaussian windowing standard deviation within the range [λ ℓ 0 , λ L ]. As for the values of parameter λ, a standard pyramidal multi-scale scheme would correspond to λ ℓ-1 = 0.5λ ℓ . In practice we have founded, as in [START_REF] Wu | Image registration using wavelet-based motion model[END_REF], that smaller values were more efficient and we use λ ℓ-1 = 0.3λ ℓ .

The quantity ξ a,ℓ k denotes the empirical mean of the analysis ensemble. The initial analysis ensemble is fixed to a null value (ξ a,L+1 k = 0) and the initial forecast is set to the forecast ensemble computed from the dynamics (x

f ,(i),L+1 k = ξ f ,(i) k
). At the coarsest scale,L, the motion estimation is hence carried out on the original images f k-1 and f k . These estimates are computed at a rough scale through a large value of the gaussian window function g λ L (2.14). At finer scales the motion estimation is performed for a decreasing value of this parameter on a pair of images constituted by the initial image f k-1 and the image f ℓ k obtained from backward registration of image f k with respect to the analyzed motion field at the previous scale, v(ξ a,ℓ+1 k

). At each scale, the Gaussian random fields, γ ℓ , attached to the measurements perturbations involved in the analysis step (4.6), are drawn with the uncertainties provided by the stochastic Lucas and Kanade formulation computed from the couple of images (f k-1 , f ℓ k ) and the current analysis ensemble (2.13). Let us note that compared to the previous single scale filtering where the proposal was based on a single ensemble Kalman update, here several updates associated to different Kalman gains computed from observations at different scales are considered. In the experimental section, three successive scales will be considered in such a filtering. The final proposal corresponds to the sum of the analysis ensemble obtained at the different scales:

ξ a k = L ℓ=ℓ 0 ξ a,ℓ
k . In the same way as the previous filter, these ensemble members are then resampled according to the importance weights computed from the likelihood associated to the original couple of images

(f k-1 , f k ).
The multiscale version of the Weighted ensemble Kalman filter (WEnKF) procedure can be summarized by the algorithm (4.1).

= 0; ξ f ,L+1 k = ξ f k ; ξ a k = 0; for ℓ = L → ℓ 0 do Compute the backward registered image f ℓ k (X k-1 ) = f(X k-1 + v(ξ a k ), k);
Compute the SLK estimate between images f k-1 and f ℓ k ; Get the observation covariance R ℓ ; Draw N samples of Gaussian random fields γ ℓ k ∼ (0, R ℓ ); Compute the analysis ensemble at scale ℓ and the forecast residual from

ξ f ,ℓ k = ξ f ,ℓ+1 k -ξ a,ℓ+1 k ; ξ a,ℓ k = ξ a,ℓ+1 k + ξ f ,(i),ℓ k + K ℓ k ξ ℓ k + γ ℓ k -Hξ f ,ℓ k ; ξ a,ℓ k = 1 N N i=1 ξ a,ℓ k ; ξ a k = ξ a,ℓ k ;
end for Compute the weights w

Experimental Results and Comparisons

In this section, we first validate and compare our technique with some state-of-theart approaches on a synthetic sequence built from a Direct Numerical Simulation (DNS) corresponding to a 2D sustained turbulence. In a second step, we test our approach on a real fluid sequence.

Synthetic DNS sequence

We present the results obtained by the application of the single scale and the multiscale WEnKF denoted as 1L-WEnKF and 3L-WEnKF respectively as the latter has been applied on a set of 3 three successive scale ranges. Those filters have been compared with stateof-the-art fluid motion estimators [START_REF] Heas | Bayesian selection of scaling laws for motion modeling in images[END_REF][START_REF] Papadakis | An optimal control technique for fluid motion estimation[END_REF][START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF] on a sequence of 100 images depicting a 2D turbulent flow seeded with particles. The flow has been computed with a direct numerical similation of the Navier-Stokes equations for a Reynolds number of 3000. This sequence can be downloaded at http://www.fluid.irisa.fr.

Quantitative comparisons of the different estimation techniques with the ground truth in terms of the Root-Mean-Square-Error (RMSE) of vorticity and velocity are both presented in figure 3.

For such particle images, it is noticeable to observe from fig. 3(a), that the RMSE in vorticity of the SLK approach is close to state of art approaches [START_REF] Heas | Bayesian selection of scaling laws for motion modeling in images[END_REF][START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF], though the RMSE values in velocity are higher (fig. 3(b)). The RMSE in vorticity by assimilating the SLK observation through 1L-WEnKF is much lower, while the error in terms of velocity estimates is close to the approach of Yuan et al. [START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF]. However, the 3L-WEnKF assimilation shows better results both in terms of vorticity or velocity. These errors are lower than all the fluid motion estimators that have been tested and are at the same level as the errors provided by the batch variationnal assimilation techniques ‡ [START_REF] Papadakis | An optimal control technique for fluid motion estimation[END_REF] (which corresponds thus to a smoothing filter as opposed to a recursive filter as in our case). Let us note that the RMSE criterion corresponds only to a large scale indicator as the errors on small scales motion structures contribute to a small amount to the whole error energy. Besides, ‡ Due to technical limits associated to batch techniques used in variational data assimilation, the results were unfortunately available only for 50 images the vorticity map is very sensitive to the presence of noise in the velocity fields. The higher errors in vorticity of the motion estimators may be observed from the different vorticity maps plotted in fig. 4. Only data assimilation techniques (either variationnal or based on stochastic filtering as in this work) enable to recover vorticity map that are closer to the ground truth maps even in the finer small scales details. The presence of noisy structures in the vorticity maps provided by all the motion measurement techniques reflects the emergence of small scale dynamical inconsistancies along time. For a visual comparison we show in fig. 4 the vorticity maps obtained by the different methods for the 50th images of the sequence. We in addition depicted the estimation at 3 various scales issued from SLK in order to observe the evolution of the structures along scales (fig.

(g-h-i))

. From these images, it appears that the vorticity estimated by the 3L-WEnKF assimilation (fig. 4 (f)) is less noisy and exhibits much finer scales vortical structures. The quantitative gain obtained by the proposed 3L-WEnKF method indicates that such fine structures are physically consistent. This constitutes a key contribution of the paper.

Let us now analyze the behavior of those techniques for real data.

Real fluid sequence

Our next set of results corresponds to a real world image sequence of a 2D turbulence generated from the wake of a soap film behind a comb. The flow is visualized through a Schlieren technique at a rate of 2500 frames per second. This experiments and the image acquisition have been performed at the Fluid dynamics laboratory of the University of Buenos Aires. A typical image of the sequence is shown in figure 5 in false color. In this experiment, we have compared a pure SLK estimation with its filtering using 1L WEnKF and 3L WEnKF techniques. The estimated vorticity maps and velocity fields corresponding to SLK, 1L WEnKF and 3L WEnKF are shown in figs. 6(a), (b) and (c), respectively. We note that though the 1L-WEnKF assimilation of SLK brings out some details at a smaller scale

Conclusion

In this paper, we have proposed an efficient multiscale extension of the Weighted Ensemble Kalman filter for fluid flow motion estimation problem. This filter is a particle filter relying on a proposal distribution built from the ensemble Kalman filtering mechanism. The particular instance we considered here incorporates measurements issued from a stochastic extension of the Lucas and Kanade estimator. It benefits from a natural multiscale formulation of the motion estimator. The overall multiscale data assimilation filter is defined in an incremental way.

We have tested our approach on synthetic and real fluid flow sequences. As for the the former situation, we have also compared our method with several state-of-the-art fluid motion estimators. From these experiences it can be pointed out that on the basis of SLK motion estimates and their associated uncertainties, the overall process presented in this article enables to recover finer scale structures and therefore, to dramatically improve the quality of the results. Considering only few levels of this hierarchical filter enables indeed to improve the errors of the corresponding single level filter. This promising behavior makes the algorithm an appealing technique for fluid motion estimation in image sequences. As for the computational load of the technique proposed, it relies on a local incremental motion estimator that is only slighly more expansive than the original Lucas and Kanade estimator due to the uncertainty computation. This estimator could be coded in real time on a Graphics Processing Unit (GPU) [START_REF] Champagnat | How to calculate dense piv vector fields at video rate[END_REF]. The Weighted ensemble particle filter is intrinsically parrallel and could be very efficiently implemented on a grid computer, in order to get small computational time. In this study, all the experiments have been performed through a very inefficient Matlab code. In the future we plan to investigate an efficient implementation of this code together with a 3D Navier Stokes.

By a change of variable Y t = X t-1 + v -X t , this expectation can be written as: where g Σ (X t ) ∼ (X t , Σ 1/2 ) is a multidimensional Gaussian of mean X t with covariance Σ and Ψ(X t , t) = Ψ(X t-1 + v, t).

E(Ψ(X t , t)|X t-1 ) = 1 2πdet(Σ)

Figure 1 :

 1 Figure 1: Displacement of the grid of points.The initial grid at time t -1 (white) is transported by the velocity field v to reach the configuration at time t represented (black), up to some uncertainties (dashed lines).

  (g) Loop to step (b) until convergence (|v n | < ε); 3. Decrease the multiresolution level : λ ℓ = αλ ℓ+1 where α < 1 is fixed by the user; 4. Loop to step 2 until convergence (λ ℓ < λ ℓ min ).

k

  are computed. Using this ensemble based forecast covariance, an ensemble based Kalman gain matrix K e k can be computed. With this Kalman gain and the observation model the forecast ensemble members are then corrected towards the current observation. This correction consists to update the forecast ensemble members x f ,(i) k , through the Kalman update equations (3.4 to 3.7), with a set of perturbed observation y k + γ (i) k obtained from samples of the observation noise {γ

k

  , i = 1, ..., N }. This supplies an analysis ensemble members {x a,(i) k

  and

Algorithm 3 . 1

 31 Weighted ensemble Kalman filter (WEnKF) procedure can be simply summarized by the algorithm 3.1. The next section presents the way we combine SLK (Stochastic The WEnKF algorithm, one iteration.Require: Ensemble at instant k -1: {x (i) k-1 , i = 1, . . . , N } observations y k Ensure: Ensemble at time k: {x (i) k , i = 1, . . . , N } EnKF step: Get x (i)k from the assimilation of y o k with an EnKF procedure; Compute the weights w (i) k according to (3.13); Resample: For j = 1 . . . N , sample with replacement index I( j) from discrete probability {w (i) k , i = 1, . . . , N } over {1, . . . , N } and set x in section 2) observations with WEnKF to perform multiscale fluid motion estimation.

k

  ) corresponds to the analysis step (respectively forecast) at scale ℓ of the particles set and the measurements ξ ℓ k are supplied by the stochastic Lucas and Kanade motion estimates between the backwarped image

Algorithm 4 . 1

 41 Multiscale WEnKF algorithm, one iteration.Require: Ensemble at instant k -1: {ξ (i) k-1 , i = 1, . . . , N } Ensure: Ensemble at time k: {ξ (i) k , i = 1, . . . ,N } EnKF step: Compute the forecast ensemble ξ f k according to (3.10); ξ a,L+1 k

k

  according to (3.13); Resample: For j = 1, . . . , N , sample with replacement index I( j) from discrete probability {w (i) k , i = 1, . . . , N } over {1, . . . , N } and set ξ

Figure 2 :

 2 Figure 2: (a): Particle image (50th of the sequence) and (b): the associated true vorticity map

Figure 3 :

 3 Figure 3: Comparison with State of the Art: RMSE in (a) Vorticity (b) Motion field .

Figure 4 :Figure 5 :

 45 Figure 4: Estimated vorticity maps for (a) Yuan et al. ; (b) Heas et. al. [14]; (c) Papadakis et al. [26]; (d) SLK; (e) 1-Level WEnKF assimilation; (f) 3-Levels WEnKF assimilation with the representation color scale. Last line : intermediate scales estimations during the SLK estimation process for three various scales (from coarser scale in (g) to finer in (i))

Figure 6 :

 6 Figure 6: Vorticity maps with their corresponding velocity fields (a) SLK (b) 1L-WEnKF and (c) 3L-WEnKF for images 76,86 and 96 (first, second and third row respectively).

  1/2 Ψ(X t-1 +v-Y t , t)exp -Y t Σ -1 Y t dY t = Ψ(•, t) * g Σ (•) (X t ).

	(B.3)

† Let us note that as the covariance Σ depends on the scalar σ(X t , t), we note the Gaussian function as g σ to lighten the notation.

Acknowledgements

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant PREVASSEMBLE (ANR-08-COSI-012).

A. Appendix: Lucas-Kanade technique

Starting from the from the optical-flow constraint equation:

the Lucas-Kanade estimator assumes the velocity v(x, t) at a given position x to be homogeneous inside a gaussian window g λ of variance λ centered in x. Therefore, the estimation of v(x, t) reads:

The Euler-Lagrange equation associated to the cost-function (v ⋆ ) reads (omitting indexes for clarity reasons):

The value of v(x, t) being constant inside the convolution window g λ , one gets:

B. Appendix: Expectation of a function of a stochastic process

The conditional expectation given X t-1 of any function Ψ(X t , t) of a stochastic process defined through Itô diffusion (2.5) and discretized through an Euler scheme

) may be written as:

where Σ is a fixed matrix-valued function. As the process X t is known up to the Brownian motion Σ 1/2 dB t , the probability p(X t |X t-1 ) is a multidimensional Gaussian of variance Σδt (δt = 1 here) and we get:

Ψ(X t , t) exp -(X t-1 + v -X t )Σ -1 (X t-1 + v -X t ) dX t .

(B.2)