Thomas Carle

Hugues Cassé

Static Extraction of Memory Access Profiles for Multi-core Interference Analysis of Real-time Tasks

Keywords: Multicore architectures, Worst-Case Execution Time, Static analysis

We present a static analysis framework for real-time task systems running on multi-core processors. Our method analyzes tasks in isolation at the binary level and generates worst-case timing and memory access profiles. These profiles can then be combined to perform an interference analysis at the task system level, as part of a multi-core Worst-Case Response Time (WCRT) analysis. In this paper we introduce a formal description of the models and algorithmic building blocks composing our framework. We also discuss how the memory access profiles generated by our method could be used to feed existing state-of-the-art WCRT frameworks. To the best of our knowledge, it is the first time that a method is documented on how to produce sound, safe and precise inputs for interference analysis methods.

Introduction

Worst-Case Execution Time (WCET) and Response Time (WCRT) analysis methods have existed for decades and are currently being used in the industry to provide static guarantees that tasks running in real-time systems will respect their deadlines. Such methods have been improved over the years in order to integrate the effects of complex hardware [START_REF]absInt aiT[END_REF][START_REF] Ballabriga | OTAWA: an Open Toolbox for Adaptive WCET Analysis[END_REF] (e.g. pipelines, caches, branch predictors) and software [START_REF] Altmeyer | Cache-related preemption delay via useful cache blocks: Survey and redefinition[END_REF] (e.g. preemption, mutual exclusion) mechanisms, but have mainly targeted single-core processors. The ongoing adoption of multicore architectures for the implementation of hard real-time systems raises new challenges for the research community. Indeed previously unseen phenomena appear in such architectures, which can have a significant impact on the execution time of the tasks that run in parallel. This so-called timing interference stems from the fact that while tasks run in parallel on separate cores, they share some hardware resources such as memories and interconnects. Classical timing analysis methods make the hypothesis that tasks run in isolation (either on completely isolated hardware or on the same core but at separate times). When tasks run in parallel, this hypothesis no longer holds, and some additional delay can be experienced if they try to access a shared resource simultaneously. As a result their actual execution time may exceed the WCET computed in isolation, thus voiding all timing guarantees. Different methods have been developed to handle this phenomenon, such as predictable hardware components [START_REF] Hahn | Design and analysis of SIC: a provably timing-predictable pipelined processor core[END_REF][START_REF] Liu | A pret architecture supporting concurrent programs with composable timing properties[END_REF][START_REF] Schoeberl | T-crest: Time-predictable multi-core architecture for embedded systems[END_REF], extensions of previously existing WCRT analysis [START_REF] Davis | An extensible framework for multicore response time analysis[END_REF][START_REF] Pellizzoni | Worst case delay analysis for memory interference in multicore systems[END_REF] and interference-free execution models enforced through careful static scheduling and synchronization [START_REF] Durrieu | Predictable flight management system implementation on a multicore processor[END_REF]. In this paper we focus on a mixed analysis/compilation framework based on the notion of Time Interest Points (TIPs), which was first introduced in [START_REF] Carle | Reducing Timing Interferences in Real-Time Applications Running on Multicore Architectures[END_REF]. In this framework, tasks are initially analyzed in isolation in order to pinpoint the instructions which may cause or suffer from interference at runtime (the TIPs). The result of this first phase is a representation of the worst case memory access profiles of the tasks in time, under the form of timed execution traces. This information is then abstracted as sequences of segments (one sequence per task) characterized by a worst-case duration and a worst-case number of memory accesses, to be combined in a static scheduling phase in which an interference analysis is performed. Finally, synchronizations are automatically injected in the code in order to enforce the schedule/response-time computed in the second phase. The main advantages of this method are:

-It is applicable to a wide range of Commercial Off The Shelf (COTS) processors. Some restrictions apply, but are not as strict as the ones imposed in predictable hardware components, -It is applicable to legacy code, with minimal automatic code modifications, where existing methods based on static scheduling require heavy transformations of the source code of the tasks to make it comply with the execution model, -The byproducts of the analyses of the first phase (traces and segments) can be used to feed state-of-the-art WCRT and Real-Time calculus [START_REF] Davis | An extensible framework for multicore response time analysis[END_REF][START_REF] Pellizzoni | Worst case delay analysis for memory interference in multicore systems[END_REF] analyses rather than a static scheduling back-end, in order to allow more dynamic implementations of the system, with minimal code adjustment (timetriggered or lock-based synchronizations, or thread yielding mechanisms and scheduler configuration).

In this paper we focus on architectures where all cores have a private scratchpad memory in which their code is loaded and are equipped with private L1 data caches and a shared memory bus implementing a greedy first-come first-served policy. In this context we provide a formal description of the models and algorithms which allow the abstraction of tasks binary code into time and memory access profiles, and discuss how these profiles can be fed to state-of-the-art analysis techniques for which, to the best of our knowledge, no method was yet provided to produce inputs.

Related Works

The real-time systems community has been working on the problem of multicore interference for nearly two decades now. A comprehensive survey on the topic has been published in [START_REF] Maiza | A survey of timing verification techniques for multi-core real-time systems[END_REF]. In this section we position our work within the state-of-the-art, and focus on two existing analysis frameworks for which our results can be particularly useful.

Reduction of Interference Through Predictable Execution:

The framework we present here can be seen as a generalization of the PRedictable Execution Model (PREM) [START_REF] Pellizzoni | A predictable execution model for cots-based embedded systems[END_REF] for multi-core architectures, or as a relaxation of the constraints of the Acquisition-Execution-Restitution (AER) [START_REF] Durrieu | Predictable flight management system implementation on a multicore processor[END_REF][START_REF] Rouxel | Hiding Communication Delays in Contention-Free Execution for SPM-Based Multi-Core Architectures[END_REF] execution model. The original idea of PREM was to avoid interference between memory accesses and asynchronous I/O traffic on a bus by carefully scheduling and enforcing the execution of tasks so that it does not occur in parallel with I/O interrupts or DMA transfers. The TIPs framework leverages this idea to the problem of multi-core interference analysis: the primary objective is to generate timing and memory access profiles of real-time tasks in order to statically schedule them on multi-core processors in a way that carefully accounts for, and possibly reduces the interference between them. The AER execution model aims at suppressing all interference by construction. The idea is to separate the execution of each task into three consecutive parts: the acquisition (A) of code and data for the task, the execution (E) of the task, and the restitution (R) of the outputs of the task to the shared memory. This separation is ensured either by the programmer or by the compiler [START_REF] Pagetti | Automated generation of time-predictable executables on multicore[END_REF]. Then the tasks are statically scheduled in a way that ensures that the A and R parts of the different tasks never occur in parallel. The TIPs framework implements the same idea, but the granularity at which it works (single memory accesses) is much finer, and it does not require to compile the task as three separate parts. This has multiple advantages such as the possibility to analyse and deploy legacy code with only small, automatic modifications (for synchronizations), and the limitation of the memory overhead due to static reservation in the AER model. Another difference is that TIPs allow the construction of programs in which some amount of interference can be tolerated (and statically quantified for compositionable processors [START_REF] Hahn | Enabling compositionality for multicore timing analysis[END_REF]).

WCRT Analysis Frameworks: In [START_REF] Davis | An extensible framework for multicore response time analysis[END_REF] the authors present a WCRT analysis framework for sporadic task systems scheduled on multi-core processors using a preemptive fixed priority algorithm (and static partitioning of tasks on the cores). The authors consider that each possible execution trace of each task in the system is available for analysis, and from this set provide precise formulas to quantify the effect of interference between tasks on the shared elements of the target processor (memories, busses, processor time). This work extends classical WCRT analyses [START_REF] Joseph | Finding Response Times in a Real-Time System[END_REF] by introducing new interference terms to cover the particularities of multi-core processors, and by making it possible to precisely account for the execution context of the tasks (i.e. which other tasks are running on the same core, or in parallel). These terms are computed by extracting worst-case information for any time interval of any given size on the execution traces of tasks. In [START_REF] Davis | An extensible framework for multicore response time analysis[END_REF] the authors discuss the empirical complexity of obtaining and manipulating the entirety of the execution traces for a task system corresponding to an industrial application. Their conclusion is that traces are a desirable abstraction of the tasks execution behavior since they can be easily manipulated and they express precisely the relation between the task and the shared resources. In particular they emphasize the fact that the worst case behavior of a task depends on its execution context, and that traces allow to exploit this. They conclude that although working on all execution traces is unfeasible for arbitrary applications, it is possible to feed the framework with a set of abstract traces which overestimate the worst case behaviors of the tasks. However nothing is said on how to obtain such an abstraction, nor on the potential costs of the various abstraction methods that could be used.

In [START_REF] Pellizzoni | Worst case delay analysis for memory interference in multicore systems[END_REF] the authors provide a method close to real-time calculus [START_REF] Thiele | Real-time calculus for scheduling hard real-time systems[END_REF] in order to compute the WCRT of a task system on a multi-core processor. Each task is represented as a sequence of time intervals, and for each time interval, a bound on the worst case number of memory accesses performed by the task is assumed to be known. Using this information, memory access arrival curves are derived and then combined to upper-bound the interference effect in time. A method is briefly sketched to derive the time intervals, which assumes precise knowledge on the tasks behavior (in particular local best and worst case execution times), but nothing is said on how this knowledge can be acquired in practice, nor on the abstraction cost of building the time intervals this way.

In Sec. 4 we discuss how the worst-case traces and the temporal segments that are generated by the TIPs framework could be good candidates to feed the analyses of [START_REF] Davis | An extensible framework for multicore response time analysis[END_REF] and [START_REF] Pellizzoni | Worst case delay analysis for memory interference in multicore systems[END_REF]. This discussion is preceded by a precise description of these models, how they can be generated, and on the various optimization objectives that can be used to tune the analysis and their potential impact on the precision of the abstracted representations of the tasks.

Static Analysis Framework

In this section we first provide an overview of the TIPs static analysis framework, and then focus on each of the separate transformations that compose it.

Overview of the Method and Models

The TIPs static analysis framework processes a real-time task system by a sequence of analyzes and transformations, which are detailed in the next sections:

-In a first step (Sec. 3.2), each task is analyzed in isolation. Starting from the disassembled binary of a task, a Control Flow Graph (CFG) is constructed.

The CFG is analyzed in order to extract TIPs, that is to say instructions which can produce or suffer from interference. In our current implementations, we focus on instructions which may generate traffic on the memory bus due to a data cache miss, but the method could be easily extended to misses from instruction caches. Other potential sources of interference such as shared L2 caches or effects from cache coherence protocols can also be modeled in the same framework, but are left for future work. -Once the TIPs have been obtained, the CFG is transformed into a TIPs-Graph (Sec. 3.2 as well): a simplified control flow graph where the nodes correspond to the TIPs of the task, and the edges represent the possible control flow between the TIPs, in an abstract version. Nodes are labelled with the number of memory accesses made by the corresponding TIP, and edges are labelled with the worst case execution time of any execution paths linking the source TIP and the destination TIP of the edge. This representation is TIP-centered, and simplifies the CFG while allowing the following analyses and transformations to remain conservative. -The TIPsGraph is then used to enumerate execution traces using a working list algorithm (Sec. 3.3). The enumerated traces exhibit the occurrence of the TIPs in all possible executions of the task. For each trace, the TIPs execution dates are a worst-case approximations. The enumerated traces can be used as timed memory access profiles for the tasks in WCRT analyses, but may remain too complex to be used in practice for other methods (such as static scheduling). -For uses for which the enumerated traces are too complex to be exploited, the traces for each task are then transformed into a sequence of so-called "time segments" (Sec. 3.4): each segment has a duration and a worst case number of memory accesses, and the sequence of segments represents an overapproximation of the number of memory accesses that can be performed by the task in the corresponding time windows. -In the TIPs framework, the tasks of the system are then subjected to static scheduling, using their representation as sequences of segments (Sec. 4). During this step, an interference analysis is performed, which assumes that the processor architecture is time-compositionable [START_REF] Hahn | Enabling compositionality for multicore timing analysis[END_REF], and its results are included in the schedule. Once an acceptable schedule (i.e. which respects all real-time constraints) has been found for the whole tasks system, synchronizations are automatically inserted in the binary code of the tasks to enforce the schedule.

In the remainder of this section we will provide more details and a formal representation for each of the aforementioned steps and models.

Extracting a TIPsGraph from a CFG

The analysis of each task τ in isolation starts working on the CFG CF G τ = {N , E} of τ , where N is the set of nodes called Basic Blocks (BBs) of the graph, and E is the set of edges e ∈ N × N which represent the control flow of the application. In this model BBs are sequences of instructions i 0 , i 1 , ..., i n ∈ I with a single entry point and a single exit point. Using MUST and MAY cache analyses [START_REF] Lv | A survey on static cache analysis for real-time systems[END_REF], TIPs are pinpointed from the rest of the instructions. As stated before, a TIP is an instruction which may create or suffer from interference.

Recall that in the scope of this paper we focus on multi-core architectures in which each core has a private L1 data cache, a private scratchpad holding the code to execute and all cores share a memory bus. In this context TIPs are the memory instructions which cannot be statically determined to always result in a hit (called in short Always Hit -AH) in the L1 data cache of the core which executes them. The objective of the first step of the analysis is to build for each task τ a TIPsGraph T G τ = {T , E T G } where T ⊆ I × N is the set of TIPs of where P AT HS(i j , i k ) is the set of possible execution paths between instructions i j .i and i k .i, and W CET (p) is a conservative approximation of the WCET of the code portion composed of the instructions of p, which can be computed using a static analysis tool.

To ensure that a TIPsGraph covers the possible executions of the whole task it represents, we add two fictive nodes i start and i end which represent the entry and exit points of the task. Both i start .µ and i end .µ are equal to 0. Fig. 1a shows a TIPsGraph along with the CFG from which it was extracted. The TIPsGraph starts with node i start and ends with node i end . The rest of the nodes composing the TIPsGraph is extracted from the CFG: in this example we assume that four memory instructions may access the bus (the cache analysis did not result in AH for these). Each of them is represented in the TIPsGraph, as well as the possible control flow between them. Each arc records such a possible transition, and is labelled with the WCET of the portion(s) of code that are executed between the TIP instructions.

In order to correctly handle loops, a TIP i head , with i head .µ = 0 is also created to represent the loop header BB, if and only if there exists at least a TIP i inside the loop with i.µ > 0. When there is no TIP inside the loop, the loop gets abstracted in the TIPsGraph, like illustrated in Fig. 1b : the control flow of the loop is no longer detailed in the TIPsGraph, but the edge representing the transition between the last TIP before the loop and the first TIP after the loop accounts for the worst case loop duration.

Enumeration of Timed Execution Traces

The next step of the analysis is to enumerate execution traces from the TIPs-Graph. The result of this enumeration is an abstract representation of the possible execution traces of the task, with two interesting properties for our analysis purposes:

-It is centered around memory accesses: only memory access instructions are represented (and loop headers, when the loop body contains memory access instructions) in the traces. In particular, control flow divergence which does not lead to memory accesses is abstracted away, and accounted for in the WCETs between TIPs. This reduces the empirical complexity of the subsequent analyses. -All transitions between TIPs are labelled with local WCETs. This guarantees that the abstraction used to represent the tasks execution is conservative: although not all actual execution traces are detailed in the analysis, the subset on which we work is a sound conservative approximation for WCET analysis. Moreover, in combination with the following steps (static scheduling, interference analysis, injection of synchronizations), this model is also sound for the analysis of interference.

We define a trace tr as a sequence of couples (t, d) ∈ T × N, where t is a TIP and d is a conservative approximation of the worst case execution date of t.i in trace tr. For a trace tr = [p 0 , p 1 , ..., p n ] with ∀i ∈ [|0, n|], p i = (t i , d i ), we denote by last(tr) the element p n . We also use tr :: p k+1 to denote the trace obtained by concatenating trace tr with element p k+1 . A basic enumeration algorithm is described in Algo. 1. It is a working list algorithm which performs a depth-first traversal of the TIPsGraph of a task. The working list contains triplets composed of a trace currently under construction, a TIPsGraph edge and a stack containing information regarding the current iteration of loops that are being traversed. The algorithm iteratively builds the set T races of the enumerated traces. Initially, W L and T races are empty. At each step of the process, the algorithm gets a trace under construction from W L, along with an edge from the TIPsGraph whose source node is the current last node of the trace, and the corresponding loop iteration context. From this, the trace is extended with the destination instruction of the edge, and pushes this new state on W L, along with all possible successor edges of the new last node of the trace. One trace is completed and thus added to the T races set when the node i end has been reached.

(b) (a1) (a2) (a3) (a4) i start ,0 i start ,0 i start ,0 i start ,0 i start ,0 i 1 ,5 i 1 ,5 i 1 ,5 i 1 ,
The tricky cases concern loop headers (L.8 to L.26): in order for the algorithm to finish, it is mandatory for the number of iterations of each loop of the task to be bounded (which is a basic requirement for WCET computation). When the trace enumeration reaches an edge whose destination node corresponds to a loop header (L.10), the algorithm checks (L.11) whether the arc in question is a return arc from inside the loop (marking the end of an iteration of the loop), or not (meaning the enumeration is entering the loop for the first iteration). If the enumeration just enters the loop, a new loop iteration context is created by pushing 0 (corresponding to the first iteration of the loop) on the context stack (L.18). The algorithm uses a stack so it can handle nested loops. If on the other hand, the current arc is a return arc, the algorithm checks if the current iteration corresponds to a valid execution : it must not exceed the maximum iteration bound for the loop. If the execution is invalid, the current trace is simply discarded (L.13), and the algorithm pops a new element from W L. In order to work, the algorithm must also be able to pop an element from the context stack when exiting a loop. This is done by detecting that the current edge exits from the loop (L.21), and by checking that the minimum iteration bound has been reached in the current stack (L.22). This minimum iteration bound is set to 0 by default, but the more precise it is, the better the outcome of the analysis. 1b. The first element of the trace, i start , corresponds to the start of the execution of the task at date 0. The next elements are the execution of i 1 .i at date 5, the execution of i 2 .i at date 693 and finally the end of the task at date 707. Traces (a1) to (a4) are a subset of all possible enumerated traces from the TIPsGraph of Fig. 1a. In order to enumerate them, we assumed that the number of loop iterations varied at least between 0 iteration (trace (a1)) and 2 iterations (trace (a3)). Trace (a1) corresponds to the execution of the task when the loop is not executed. Traces (a2) and (a3) correspond to the execution of the task when the left branch of the loop is taken respectively once and twice before exiting the loop. Trace (a4) corresponds to the execution of the task when the right branch of the loop is taken once before exiting. Notice that the dates for each element of these traces are worst case dates, meaning that the corresponding instructions can in practice execute before that date, but are guaranteed to never execute after. This means that any such trace covers multiple execution patterns, which reduces the empirical complexity of the next steps of the analysis (regardless of the analysis framework). However, three important issues must be raised at this point:

-Knowing only a worst case date for the memory accesses may increase the imprecision of the interference analysis, since it must consider that an access can occur at any time before the worst case date. One way to mitigate this issue is to inject synchronizations inside the code of the tasks to reduce the size of the time intervals during which accesses may occur. -As illustrated by the 4 traces (a1) to (a4), the same instruction can have a different worst case date in different traces (e.g. i 4 .i), which can also lead to imprecision in the interference analysis. Methods must be found to mitigate this issue, either at the code generation level (once again, synchronizations) or at the analysis level (careful accounting of the worst case number of memory access of the task on a given time interval). -Enumerating traces when a loop has a different minimum and maximum iteration bounds dramatically increases the empirical complexity of the enumeration algorithm: the enumeration of all possible sub-traces after the exit of the loop must be performed entirely as many times as there are ways to exit the loop (i.e. for each iteration between the minimum and maximum loop bound), even though the enumeration of these sub-traces is exactly the same each time, since they correspond to exactly the same portion of the TIPsGraph (when infeasible paths are not considered). In the given example, the portion of the TIPsGraph located after the loop is very small, but in practice we have noticed that it is not the case for arbitrary applications, and that the enumeration may become infeasible in acceptable time when the min and max bounds for a loop differ too much.

The enumerated traces are a first, rather raw representation of the timing and memory access profile of the tasks of the analyzed system. They can be used to perform a WCRT analysis following the method described in [START_REF] Davis | An extensible framework for multicore response time analysis[END_REF], even though caution must be taken: these are not real execution traces, but worst-case approximations of execution traces. This means that the method of [START_REF] Davis | An extensible framework for multicore response time analysis[END_REF] will have to be adapted to take into account this specificity, or that synchronizations will have to be added to the task code in order to enforce some of the worst case dates for the memory accesses.

In the TIPs framework however, the objective is to perform static scheduling, in order to analyze and try to limit the interference between tasks. To do so, we need to transform the enumerated traces into entities that will be practical to schedule, such as temporal segments.

Temporal Segments

We now present the kind of temporal segments that are used in the TIPs framework in order to represent the time and memory access profile of tasks and to statically schedule them. A temporal segment s i is a triplet (s i .start, s i .dur, s i .µ) where s i .start is the start date of the segment, s i .dur is its duration, and s i .µ is a map which contains the number of memory accesses that can happen on the time interval [s i .start, s i .start + s i .dur] for each trace. In the following, we also denote by s i .end = s i .start + s i .dur the end date of segment s i .

Any task τ (resp. any enumerated trace tr ∈ T races(τ )) can be abstracted using a sequence of segments Segs τ = [s 0 (τ ), ..., s n (τ )] (resp. Segs tr = [s 0 (tr), ..., s k (tr)]), with the property that segments of a sequence do not overlap and the first segment starts at date 0 i.e. s 0 (τ ).start = 0 and ∀i ∈ [|1, n|], s i (τ ).start ≥ s i-1 (τ ).end.

The shape of the segments sequence of each task will have an impact on the scheduling and interference analysis phase. A trade-off must be found between:

-The number of segments for each task. Scheduling elements (tasks, or segments) on a multi-core target is a NP-hard problem, so increasing the number of segments to schedule can increase the time it takes to build a schedule, potentially to a point where it is no longer feasible in practice. -The length of the segments. During the interference analysis, any two segments from different tasks scheduled on overlapping time intervals on different cores will be considered as being in interference. By definition, smaller segments occupy a core for less time than larger segments, and are thus less exposed to interference from other cores. Moreover, smaller segments offer more flexibility to the scheduler to reduce the impact of interference. -The worst case number of memory accesses on each segment. The length of the segments and the number and position of the synchronizations used to enforce them have an impact on the number of memory accesses attributed to each segment. This number must be conservative for each segment, so a memory access from a single instruction can be counted in multiple segments if the execution date of the instruction cannot be proven to happen in the time interval of only one segment. This can increase the imprecision of the method if one is not careful when shaping the segments and selecting the synchronization points. -The number of synchronizations that will be required to guarantee that the code corresponding to the segments does not start before it is intended to. Each synchronization corresponds to additional code for the task, so their number must remain limited. Without optimization, code must be added (automatically) to the task code to ensure that in each execution trace a synchronization will be executed to enforce the start date of each segment.

In the remainder of this section, we provide algorithms that enable the extraction of valid segment representations for tasks. These are baseline algorithms which do not perform any optimization with regard to the aforementioned tradeoffs. In the description of the algorithms we use the empty sequence ([]) and concatenation of an element e at the right-end of a sequence seq (::).

These algorithms rely on the Intersect operator which is defined in Def. 1. This operator computes the intersection of two segments: if the segments correspond to non-overlapping time intervals, the return value is empty. Otherwise, the operator returns a segment whose time interval is the intersection of the time intervals of the two input segments, and its summary of worst case number of memory accesses is the union of the summaries of worst case memory accesses of the input segments. Our algorithms also use the Segments procedure described in Algo. 3. This procedure transforms a trace tr of a task τ into a sequence of segments in the following manner: for each node n in the trace, it creates two segments: s 1 which starts at the date of the node, spans the worst case duration of the accesses of this node and has s 1 .µ = {tr : n.µ} (marking that on this time interval trace tr makes at most n.µ accesses), and s 2 which starts just after and spans until the date of the next node and has s 2 .µ = {tr : 0}. The procedure also called with parameter d max which is the maximum of the dates of the last nodes of all traces of τ (i.e. the WCET of τ in the absence of interference). This is used to extend the last segment so that it spans until d max .

Definition 1. ∀s i , s j ∈ Segs, Intersect(s i , s j ) =                ∅ if s i .start ≥ s j .end ∨ s j .start ≥ s i .end (s i .start, s i .dur, s i .µ ∪ s j .µ) if s i .start ≥ s j .start ∧ s j .end ≥ s i .end (s i .start, s j .end -s i .start, s i .µ ∪ s j .µ) if s i .start ≥ s j .start ∧ s i .end ≥ s j .end (s j .start, s i .end -s j .start, s i .µ ∪ s j .µ) if s i .start < s j .start ∧ s i .end ≤ s j .end (s j .start, s j .dur, s i .µ ∪ s j .µ) if s i .start < s j .start ∧ s i .end > s j .end
The top-level algorithm is described in Algo. 2: starting with an arbitrary trace tr 1 from the set of traces of τ , it transforms tr 1 into a segments representation using procedure Segments (described in Algo. 3): Segs τ . Then each other trace tr i of τ is transformed into a sequence of segments, and Segs τ is updated with the intersection of the current segments of Segs τ and the segments that represent tr i . When this is done, a procedure tries to reduce the number of preserving all segments s with max access(s) = 0 and s.dur ≥ ∆, for all other segments, fusing consecutive segments until the result of the fusion has a length of at least ∆ or there is no more available segment to fuse. When fusing segments, the information about the worst case number of memory accesses is combined trace-wise instead of blindly summed in order to limit over-approximations.

We illustrate this algorithm in the examples of Fig. 3. Segment sequences (S a1 ) and (S a2 ) are extracted directly from traces (a1) and (a2) of Fig. 2 using Algo. 3. The result of their intersection is provided as (S a1+a2 ). In this sequence, the first access is displayed in gray to show that this segment corresponds to either one access from trace (a1) or one access from trace (a2). The sequence labelled (S a ) is obtained by iterating the intersection of the traces (a1), (a2), (a3) and (a4). Different colors mean that accesses from different traces may occur. Finally trace (S ′ a ) is obtained by fusing together the smaller segments and preserving large segments which are guaranteed to not perform any memory access.

Algorithm 3 Extract a segments representation for a single trace

Require: tr1 ∈ T races(τ ), tr1 = (i0, d0), (i1, d1), ..., (in, dn); dmax Ensure: Segstr

Exploitation of the Memory Access Profiles

In this section we conclude the description of the TIPs framework by a quick discussion about the uses that can be made of the enumerated traces and segments in order to perform an interference-aware analysis of the task system. We thus discuss how the enumerated traces and segments representations can be fed to various "back-ends" such as static schedulers or the WCRT analysis frameworks described in Sec. 2.

Interference Analysis for Static Scheduling

Regardless of the task model (single period, sporadic, multi-periodic, dependent or independent tasks) any existing static scheduling method can be adapted to the TIPs model: instead of scheduling one time interval for a task as in classical models, all segments for a task are scheduled in order. The potential interference is computed using the information contained in the segments, and must reflect the bus arbitration policy. Only segments scheduled on separate cores and whose time intervals overlap are considered to interfere. The potential interference is accounted for either by increasing the size of the segments on-thefly (e.g. [START_REF] Carle | Reducing Timing Interferences in Real-Time Applications Running on Multicore Architectures[END_REF][START_REF] Dupont De Dinechin | Scaling up the memory interference analysis for hard real-time many-core systems[END_REF]), or by consuming an interference budget which can be added to the tasks' WCET prior to scheduling [START_REF] Didier | Correct-by-construction parallelization of hard real-time avionics applications on off-the-shelf predictable hardware[END_REF] (e.g. a 10% overhead on the computed WCET for interference tolerance).

Multi-core WCRT Analysis Techniques

In order to use the interference formulas presented in [START_REF] Davis | An extensible framework for multicore response time analysis[END_REF] (at least for the part regarding interference on the buses and memories), an upper bound on the number of accesses that can be made in any time-interval of any size must be found for each task. In the TIPs framework, this information is available at 2 different abstraction levels: the enumerated worst-case traces (Sec. 3.3) which offer a finer level of granularity, and thus can lead to more precise bounds at the cost of a more complex computation, and the segments (Sec. 3.4) which represent the memory access profiles at a higher level of abstraction. In both cases, since the timings obtained in the TIPs framework are worst-case dates, synchronizations must be added to the code to be able to lower-bound with certainty the occurrence of accesses in time. The same is true for the real-time calculus method of [START_REF] Pellizzoni | Worst case delay analysis for memory interference in multicore systems[END_REF], for which the segment representation obtained through the TIPs framework is a natural input format.

Conclusion and Future Works

We presented the TIPs framework: a collection of models and algorithms for the extraction of precise timing and memory access profiles of real-time tasks. For each level of abstraction, we provided a formalization of the corresponding models and algorithms, as well as a discussion on the cost of the presented abstractions. We finally discussed how the obtained memory access profiles could be used as inputs for existing state-of-the-art analysis frameworks and tools. In the future, we will work on optimizations and the evaluation of their impact on the different "back-ends". In particular, we are currently working on a multicriterion optimization of the transformation of enumerated traces into segments, trying to minimize the overestimation of the number of accesses in each segment and the number of synchronizations to add to the tasks' code.

Fig. 1 :

 1 Fig. 1: Example of CFGs and their corresponding TIPsGraphs

Figure 2

 2 Figure 2 displays 5 traces enumerated from the TIPsGraphs of Fig. 1. The last trace (at the bottom), labelled (b) is the only trace that can be enumerated from the TIPsGraph of Fig.1b. The first element of the trace, i start , corresponds to the start of the execution of the task at date 0. The next elements are the execution of i 1 .i at date 5, the execution of i 2 .i at date 693 and finally the end of the task at date 707. Traces (a1) to (a4) are a subset of all possible enumerated traces from the TIPsGraph of Fig.1a. In order to enumerate them, we assumed that the number of loop iterations varied at least between 0 iteration (trace (a1)) and 2 iterations (trace (a3)). Trace (a1) corresponds to the execution of the task when the loop is not executed. Traces (a2) and (a3) correspond to the execution of the task when the left branch of the loop is taken respectively once and twice before exiting the loop. Trace (a4) corresponds to the execution of the task when the right branch of the loop is taken once before exiting. Notice that

Algorithm 2

 2 Segments creation for a taskRequire: T races(τ ), tr1 ∈ T races(τ ), ∆ ∈ N, dmax ∈ N Ensure: Segsτ 1: Segsτ = Segments(tr1, dmax) 2: for all tri ∈ T races(τ ), tri ̸ = tr1 do 3: Segsτ ← Intersect(Segsτ , Segments(tri, dmax)) 4: end for 5: Segsτ ← F usion(Segsτ , ∆) 6: return Segsτ segments using a minimum size ∆, by :

1 1 : 4 : 1 ← 8 : 1 (Fig. 3 :

 141813 Fig.3: Examples of memory access profiles obtained from the traces of Fig.2


Noémie Debroux
email: noemie.debroux@uca.fr

Carole Le Guyader
email: carole.le-guyader@insa-rouen.fr

Luminita A Vese
email: lvese@math.ucla.edu

Multiscale Registration

Keywords: Multiscale analysis, registration, nonlinear elasticity, Ogden materials, bi-Lipschitz homeomorphisms, Γ -convergence

In the seminal paper E. Tadmor, S. Nezzar and L. Vese, A multiscale image representation using hierarchical (BV, L 2 ) decompositions, Multiscale Model. Simul., 2(4), 554-579, (2004), the authors introduce a multiscale image decomposition model providing a hierarchical decomposition of a given image into the sum of scale-varying components. In line with this framework, we extend the approach to the case of registration, task which consists of mapping salient features of an image onto the corresponding ones in another, the underlying goal being to obtain such a kind of hierarchical decomposition of the deformation relating the two considered images (-from the coarser one that encodes the main structural/geometrical deformation, to the more refined one -). To achieve this goal, we introduce a functional minimisation problem in a hyperelasticity setting by viewing the shapes to be matched as Ogden materials. This approach is complemented by hard constraints on the L ∞ -norm of both the Jacobian and its inverse, ensuring that the deformation is a bi-Lipschitz homeomorphism. Theoretical results emphasising the mathematical soundness of the model are provided, among which the existence of minimisers, a Γ -convergence result and an analysis of a suitable numerical algorithm, along with numerical simulations demonstrating the ability of the model to produce accurate hierarchical representations of deformations.

Introduction

Assuming that the considered grey-level images are L 2 -observations that may include scale-varying objects, ranging from large homogeneous regions that are L. Vese acknowledges support from the National Science Foundation under Grant # 2012868. This project was co-financed by the European Union with the European regional development fund (ERDF, 18P03390/18E01750/18P02733) and by the Haute-Normandie Régional Council via the M2SINUM project.

faithfully modelled by the smaller functional space BV , to oscillatory patterns/ texture that require more involved functional spaces, multiscale representation can be viewed as the task which aims to accurately describe these different levels of details. This means representing the noticeable characteristics of a given image into suitable intermediate spaces/subclasses lying in between the rougher space L 2 (or L p ) and the smaller space BV of functions of bounded variation. While, as previously mentioned, this latter space is the proper one to describe homogeneous regions with sharp edges, finer details (whether it be oscillatory patterns or noise) are generally captured by a function or distribution belonging to a larger space such as L p or a dual space of distributions (such as {div g | g ∈ L ∞ }, {div g | g ∈ BMO (Bounded Mean Oscillation)} or {∆g | g Zygmund function}, as suggested by Meyer ([11])). To identify such intermediate functional spaces, a standard strategy is interpolation (see [START_REF] Bennet | Interpolation of Operators[END_REF] or [START_REF] Bergh | Interpolation Spaces. An Introduction[END_REF]). The theory of interpolation studies the family of spaces Y that are intermediate spaces between Y 0 and Y 1 , Banach spaces continuously embedded in the same Hausdorff topological vector space Z, in the sense that

Y 0 ∩ Y 1 ⊂ Y ⊂ Y 0 + Y 1 .
The model investigated in [START_REF] Tadmor | A Multiscale Image Representation Using Hierarchical (BV , L 2 ) Decompositions[END_REF] falls within this framework and proves to be a special instance of the canonical form

J p (f, η; X, Y ) = inf u+v=f { v p X + η u Y }
(closely related to the so-called K-functional (cf. [START_REF] Vese | Variational Methods in Image Processing[END_REF]Chapter 5])) for which the intermediate spaces (X, Y ) θ , θ ∈ [0, 1] with Y ⊂ X, ranging from (X, Y ) 0 = X -the larger space -to (X, Y ) 1 = Y -the smaller one -are quantified by the behaviour of J p (•, η; •, •) as η tends to 0. More precisely, in [START_REF] Tadmor | A Multiscale Image Representation Using Hierarchical (BV , L 2 ) Decompositions[END_REF], the authors consider the larger functional space X to be L 2 , while Y is chosen to be BV , and intend to measure how accurately an L 2 -object can be approximated by its BV characteristics by considering the family of J-functionals

J(f, λ) = J 2 (f, λ; BV, L 2 ) = inf u+v=f λ v 2 L 2 + u BV ,
with increasing λ's. Note that the focus is shifted in comparison to the canonical form, from smaller η's to larger λ's.

The weighting parameter λ serves as a scale level to discriminate properly the two components: u, that captures the structural features of the observation f , while v encodes texture and oscillatory patterns, the discrimination between the two constituents being dictated by the scale parameter λ. This principle is embedded in a dyadic refinement process with λ 0 a given initial scale, that reads as:

f = u 0 + v 0 , [u 0 , v 0 ] = arg min u+v=f J(f, λ 0 ), v j = u j+1 + v j+1 , [u j+1 , v j+1 ] = arg min u+v=vj J(v j , λ 0 2 j+1 ), j = 0, 1, • • • ,
producing at the end of the k th step, the following hierarchical decomposition:

f = u 0 + v 0 = u 0 + u 1 + v 1 = • • • = u 0 + u 1 + • • • + u k + v k ,
the u k 's resolving finer edges. A remarkable result ([24, Theorem 2.2]) states that under the assumption f ∈ BV , the (BV, L 2 ) hierarchical decomposition of f , k j=0 u j , strongly converges to f in L 2 . In accordance with this framework, we propose transposing it to the case of registration, the task which aims to determine an optimal (in a sense to be specified) diffeomorphic deformation ϕ that aligns the structures visible in an image called Reference into their counterparts in another one -called Template -(see [START_REF] Modersitzki | Numerical Methods for Image Registration[END_REF][START_REF] Modersitzki | FAIR: Flexible Algorithms for Image Registration[END_REF][START_REF] Sotiras | Deformable medical image registration: A survey[END_REF] for a relevant analysis of the registration problem), the underlying goal being twofold: (i) obtaining a hierarchical decomposition of the sought deformation as

ϕ 0 • • • • ϕ k • • • • • ϕ n ,
ϕ 0 encoding the main geometry-driven/structural deformation, while the ϕ 0 • • • • • ϕ k 's capture more refined deformations; (ii) assessing to what extent this hierarchical decomposition conveys the hidden structure of the deformation to be recovered, and how it could be used to derive image statistics, to measure variability inside a population, to retrieve the inherent dynamics of a phenomenon, and so on. The introduced algorithm is patterned after the one in [START_REF] Tadmor | A Multiscale Image Representation Using Hierarchical (BV , L 2 ) Decompositions[END_REF] in the sense that it embeds successive applications of a refinement step -the composition of deformations is now a substitute for the sum of the scale-varying constituents of the seminal model [START_REF] Tadmor | A Multiscale Image Representation Using Hierarchical (BV , L 2 ) Decompositions[END_REF] -, and can be schematised as follows. Given R (resp. T ) the Reference (resp. Template) image and its related multilayered/hierarchical decomposition k j=0 R j (resp. k j=0 T j ), the algorithm reads as, F being a general functional that will be specified later on,

ϕ k = arg min ϕ F (ϕ0 • • • • • ϕ k-1 • ϕ; k j=0 Rj , k j=0 Tj ).
Before depicting in depth our model, we would like to point out that prior related works ( [START_REF] Paquin | Multiscale Image Registration[END_REF][START_REF] Paquin | Hybrid multiscale landmark and deformable image registration[END_REF][START_REF] Paquin | Multiscale deformable registration of noisy medical images[END_REF], [START_REF] Modin | A multiscale theory for image registration and nonlinear inverse problems[END_REF]) suggest to foster the use of this image multiscale representation in a registration setting. The work [START_REF] Paquin | Multiscale Image Registration[END_REF] is focused on (rigid) multiscale image registration for the registration of medical images degraded by significant levels of noise. Where conventional methods fail to correctly align the image pair, i.e., to remove the artificial differences while highlighting the real differences due to intrinsic variations of the objects -whether it be ordinary registration techniques alone or sequential treatments including a prior step to denoise the image pair and then the application of a classical registration algorithm -, their proposed joint model achieves accurate registration. Given f an image with hierarchical expansion k j=0 u j -recall that for small k, k j=0 u j is a coarse representation of the image f that actually encodes the structural/geometrical information of f whilst removing small details -to be registered with g approximated by k j=0 v j , the authors suggest registering the truncated hierarchical representations rather than f and g themselves, to get an accurate estimation of the deformation aligning them. Two main differences can be observed in comparison to our model: (i) first, each hierarchical step k is treated independently of the previous ones, that is, without connecting the optimal deformations Φ 0 , • • • , Φ k-1 obtained at previous steps and the optimal deformation Φ k produced at step k, whereas in our approach, motivated by physical arguments stated below, deformation Φ k is built up by taking advantage of the Φ 0 , • • • , Φ k-1 's; (ii) second, and this is in some way a corollary of the first point, the optimal deformation Φ meant to bring f and g into spatial alignment is computed as a weighted average of the form 1 m m-1 l=0 b l Φ l , the b l 's being weights appropriately chosen, while we promote composition of deformations. This work is then extended to the case of landmark-driven registration ( [START_REF] Paquin | Hybrid multiscale landmark and deformable image registration[END_REF]) in a B-spline setting and then, to the case of non-rigid deformations ( [START_REF] Paquin | Multiscale deformable registration of noisy medical images[END_REF]). The main motivation in our work to promote composition of deformations stems from the following fact stated in [START_REF] Vercauteren | Diffeomorphic Demons: Efficient Non-parametric Image Registration[END_REF]. As mapping a point through a first transformation and then through a second one is equivalent to mapping the point through the composition of these two spatial transformations, the most natural and geometrically meaningful operation the space of non-parametric spatial transformations can be endowed with, is the composition. Note, on the contrary, that addition of spatial transformations has no geometric meaning. More recently in [START_REF] Modin | A multiscale theory for image registration and nonlinear inverse problems[END_REF], Modin et al. propose constructing analogous hierarchical expansions for diffeomorphims, the sum being replaced by composition of maps, in the context of image registration. Their method can be seen as a series of Large Deformation Diffeomorphic Metric Mapping (LDDMM) steps with varying weighting parameters (note that the adjective multiscale applies more to the way these weighting parameters are set at each iteration, i.e., on how strong penalisations on the data fidelity term and on the deviation of the deformation composition from the identity mapping are, since they do not exploit the multiscale expansions of the images contrary to us). Although our method and theirs have in common the composition of deformations to refine the registration process, there are several key differentiating points in addition to the above: among them, the deformation model we embrace to describe the setting in which the objects to be matched are interpreted and viewed. In our case, it originates from physical considerations (i.e. elastic models in which the shapes to be matched are considered as the observations of a same body before and after being subjected to constraints) contrary to the model developed by Modin et al. which, to our point of view, is rather rooted in pure mathematical considerations and disconnected in some way from the physics of the problem. This choice of devising a model connected to physical considerations and especially to the hyperelasticity framework was motivated by the fact that it proves to be suitable when dealing with large and nonlinear deformations and that many applied problems, e.g. biological tissue behaviour, are modelled within this setting. In [START_REF] Lam | Multiscale Representation of Deformation via Beltrami Coefficients[END_REF], the authors propose decomposing the transformation f of the domain D ⊂ R 2 in which the object is embedded, using quasi-conformal theories. As a desirable one-to-one representation of an orientation-preserving mapping f : D → D can be achieved via its Beltrami Coefficient (BC) µ(f ) : D → C (-measure of non conformality, i.e., to what extent a deformation deviates from a conformal map, and uniquely related to f by µ(f

)(z) = ∂f ∂ z / ∂f ∂z -)
, the authors suggest applying a wavelet transform to the Beltrami coefficient, yielding a decomposition of the BC into different components of different frequencies compactly supported in different sub-domains. Quasi-conformal mappings related to different components of the BC (different scales) can thus be reconstructed, those yielding a multiscale decomposition of the deformation. Finally and for the sake of completeness, we refer the reader to [START_REF] Athavale | Multiscale properties of weighted total variation flow with applications to denoising and registration[END_REF], [START_REF] Gris | A Sub-Riemannian Modular Framework for Diffeomorphism-Based Analysis of Shape Ensembles[END_REF], [START_REF] Risser | Simultaneous Multi-scale Registration Using Large Deformation Diffeomorphic Metric Mapping[END_REF] and [START_REF] Sommer | Sparse Multi-Scale Diffeomorphic Registration: The Kernel Bundle Framework[END_REF] for alternative approaches. We now turn to the mathematical formulation and analysis of the proposed physics-based multiscale modelling. We would like to point out that the core of the paper is on the fine theoretical properties exhibited by our modelling. The proofs being long, we have deliberately chosen to focus primarily on the asymptotic result, which, in our view, is the most significant one. A section is dedicated to the design of a suitable algorithm, but this does not constitute the crux of the manuscript. Extensive numerical simulations are currently being carried out to put this mathematical model to the test, and to assess its ability to unveil the hidden structure of the deformation to be retrieved.

Mathematical Modelling

Depiction of the Model

Let Ω be a connected bounded open subset of R 2 of class C 1 , thus satisfying the cone property. Let us denote by R : Ω → R the reference image and by T : Ω → R the Template image, assumed to belong to the functional space BV (Ω). For theoretical and numerical purposes, we assume that T is compactly supported on Ω to ensure that T • ϕ is always well-defined. Of course, in practice, the sought deformation should be with values in Ω, but from a mathematical point of view, if we work with such spaces, we lose the structure of vector spaces. Nonetheless, we can show that our model retrieves deformations with values in Ω -based on Ball's results [START_REF] Ball | Global invertibility of Sobolev functions and the interpenetration of matter[END_REF]. A deformation is a smooth mapping that is orientationpreserving and injective, except possibly on ∂Ω. The deformation gradient is ∇ϕ : Ω → M 2 (R), the set M 2 (R) being the set of real square matrices of order 2. The sought deformation is seen as an argument of minimum of a specifically designed cost function comprising a regularisation on ϕ prescribing the nature of the deformation, and a term measuring alignment, or how the available data are used to drive the registration process. To allow large deformations, the shapes to be matched are viewed as isotropic (exhibiting the same mechanical properties in each direction), homogeneous (showing the same behaviour everywhere inside the material), hyperelastic materials (they exhibit both nonlinear behaviour and large changes in shape), and more precisely, as Ogden ones. Note that rubber, filled elastomers and biological tissues are often modelled within the hyperelastic framework, which motivates our approach. This perspective dictates the way the regularisation on the deformation ϕ is devised, based on the stored energy function of an Ogden material. Recall that the general expression of the stored energy function of such a material (see [START_REF] Ciarlet | Elasticité Tridimensionnelle[END_REF]) is given by W O (F ) = K i=1 a i F γi + Γ (det F ), with ∀i ∈ {1, • • • , K}, a i > 0, γ i ≥ 0, and Γ :]0, +∞[→ R, being a convex function satisfying lim δ→0 + Γ (δ) = lim δ→+∞ Γ (δ) = +∞. The first term penalises changes in length, while the second one restricts the changes in area.

The latter one also ensures preservation of topology by enforcing positivity of the Jacobian almost everywhere. In this work, we focus more specifically on the following energy, • denoting the Frobenius norm,

W Op (F ) = a 1 F 4 + a 2 (detF -1) 2 + a3 (detF ) 10 -4a 1 -a 3 if detF >0
+∞ otherwise, which fulfils the previous assumptions. The two latter terms govern the distribution of the Jacobian determinant. While the middle term promotes Jacobian determinants close to 1, the rightmost one prevents singularities and large contractions by penalising small values of the determinant. The constants are added to comply with the energy property W Op (I) = 0, where I stands for the identity matrix. In the following, we will set W Op (F ) = W Op (F, det F ). This stored energy function is complemented by the term

1 { . L ∞ (Ω,M 2 (R)) ≤α} (F ) +1 { . L ∞ (Ω,M 2 (R)) ≤β} (F -1
), with α ≥ 1, and β ≥ 1, where 1 A denotes the convex characteristic function of a convex set A, the underlying idea being to recover deformations that are bi-Lipschitz homeomorphisms.

Remark 1. In terms of functional spaces, if ϕ ∈ W 1,∞ (Ω, R 2 ) (suitable space owing to the L ∞ hard constraints), det∇ϕ is automatically an element of L ∞ (Ω).

Penalising the L ∞ norm of ∇ϕ thus entails control over the Jacobian determinant. This additional term implicitly gives an upper and lower bound on the Jabobian determinant ensuring thus topology preservation.

The aforementioned regulariser is then applied along with a classical L 2discrepancy measure, yielding the following minimisation problem

inf ϕ∈W F (ϕ) = λ 2 T • ϕ -R 2 L 2 (Ω) + Ω W Op (∇ϕ, det(∇ϕ)) dx (P) + 1 { • L ∞ (Ω,M 2 (R)) ≤α} (∇ϕ) + 1 { • L ∞ (Ω,M 2 (R)) ≤β} ((∇ϕ) -1 ) , with W = {ψ ∈ Id+W 1,∞ 0 (Ω, R 2 ) | ∇ϕ L ∞ (Ω,M2(R)) ≤ α, (∇ϕ) -1 L ∞ (Ω,M2(R))
≤ β, det(∇ϕ) > 0 a.e. in Ω}, and λ > 0.

Theoretical Results

The first theoretical result claims that problem (P) admits at least one minimiser. Due to the limited number of pages, we only sketch the proof here. Theorem 1. Problem (P) admits at least one minimiser in W.

Proof. The proof mainly relies on the following elements. Ball's results ( [START_REF] Ball | Global invertibility of Sobolev functions and the interpenetration of matter[END_REF]) enable one to conclude that any minimising sequence is such that ∀k ∈ N, ϕ k is a homeomorphism from Ω to Ω -and even from Ω to Ω owing to [3, Theorem 2]-, and more particularly, a bi-Lipschitz homeomorphism. Also, ϕ k * k→+∞ φ in W 1,∞ (Ω, R 2 ) and similar arguments to those previously used enable one to get that φ is a bi-Lipschitz homeomorphism from Ω to Ω and even from Ω to Ω. By continuity of the trace operator, one gets that φ

∈ Id + W 1,∞ 0 (Ω, R 2 ). Now, ϕ -1
k is uniformly bounded according to k in W 1,∞ (Ω, R 2 ) using Poincaré-Wirtinger inequality, thus there exist a subsequence still denoted by ϕ -1

k and ū ∈ W 1,∞ (Ω, R 2 ) such that ϕ -1 k * k→+∞ ū in W 1,∞ (Ω, R 2
). On the one hand,

ϕ -1 k • φ → k→+∞ Id in L ∞ (Ω, R 2
) due to the β-Lipschitz property of ϕ -1 k , thus almost everywhere in Ω up to a subsequence. On the other hand, ϕ -1

k uniformly converges to ū in C 0 ( Ω, R 2 ), yielding ϕ -1 k • φ → k→+∞ ū • φ pointwise in Ω. By
uniqueness of the limit, ū • φ = Id a.e., leading to ū = φ-1 a.e. in Ω. It remains to address the question of weak lower semicontinuity. Since T ∈ BV (Ω) and all ϕ k and φ are bi-Lipschitz homeomorphisms, we get that ∀k ∈ N, [START_REF] Ambrosio | DC calculus[END_REF]Theorem 2.6]). We then prove that

T • ϕ k ∈ BV (Ω) ⊂ L 2 (Ω) and T • φ ∈ BV (Ω) (see
ϕ k • φ-1 strongly converges to Id in C 0,α ( Ω, R 2 ), which is a cornerstone step to prove that T • ϕ k → k→+∞ T • φ in L 2 (Ω).
We now include a multiscale representation of the deformation in our model, relying on the hierarchical decomposition into the sum of scale-varying components introduced in [START_REF] Tadmor | A Multiscale Image Representation Using Hierarchical (BV , L 2 ) Decompositions[END_REF]. Let (T j ) j ∈ BV (Ω) and (R j ) j ∈ BV (Ω) be the sequence of varying scale structural features of respectively T and R coming from the following problems -S standing for either R or T below -:

     (S 0 , v 0 ) = arg min (u,v)∈BV (Ω)×L 2 (Ω) | S=u+v {λ 0 v 2 2 + T V (u)}, (S j+1 , v j+1 ) = arg min (u,v)∈BV (Ω)×L 2 (Ω) | vj =u+v {2 j+1 λ 0 v 2 2 + T V (u)}, j = 1, . . . ,
λ 0 being an initial scale parameter provided by the user. We refer the reader to [START_REF] Tadmor | A Multiscale Image Representation Using Hierarchical (BV , L 2 ) Decompositions[END_REF] for a thorough description and analysis of this multiscale representation of an image.

In this work, we assume that T and R have similar scale structures and that each level of the following hierarchical decomposition of T , k j=0

T j , can be matched to the corresponding level of hierarchical decomposition of R, k j=0 R j . We then derive a corresponding hierarchical decomposition of the deformation, going from global structural deformations to more locally refined ones, and based on the composition operator -a more natural and physically meaningful operator than the addition for transformations -. This multiscale decomposition is described as follows :

     ϕ 0 = arg min ϕ∈W {F(ϕ, T 0 , R 0 )}, ϕ k = arg min ϕ∈X k {F(ϕ 0 • ϕ 1 • . . . • ϕ k-1 • ϕ, k j=0 T j , k j=0 R j )}, (P k ) with F(ϕ, T, R) = λ 2 T • ϕ -R 2 L 2 (Ω) + Ω W Op (∇ϕ, det(∇ϕ)) dx + 1 { • L ∞ (Ω,M 2 (R)) ≤α} (∇ϕ) + 1 { • L ∞ (Ω,M 2 (R)) ≤β} ((∇ϕ) -1 ),
and

X k = {ψ ∈ Id + W 1,∞ 0 (Ω, R 2 ) | det∇ψ > 0 a.e., (∇(ϕ0•ϕ1•...•ϕ k-1 •ψ)) -1 ∈ L ∞ (Ω, M 2 (R))}. The next theoretical result is a Γ -convergence result.
Theorem 2. Problem (P k ) admits at least one minimiser and the associated functional Γ -converges to the functional related to problem (P) according to De Giorgi's definition.

Proof. For the existence of minimisers, the same arguments as those previously used can be applied since orientation-preserving bi-Lipschitz homeomorphisms form a group stable for the composition, and BV (Ω) is stable for the sum. The Γ -convergence result mainly relies on [24, Theorem 2.2] stating the strong convergence of hierarchical decompositions to the initial image in L 2 (Ω), as well as Ball's results [START_REF] Ball | Global invertibility of Sobolev functions and the interpenetration of matter[END_REF]. These enable one to prove both inequalities required for the Γ -convergence property. This result ensures that our multiscale decomposition gives a sequence of bi-Lipschitz homeomorphic deformations which model more and more locally refined distortions until it resembles the original deformation that maps T to R. We now turn to the numerical resolution of our problem (P k ).

Numerical Resolution

Due to its non-linearity in both ϕ and ∇ϕ and its non-differentiability, problem (P k ) is hard to solve numerically. Inspired by a prior work by Negrón Marrero [START_REF] Marrero | A numerical method for detecting singular minimizers of multidimensional problems in nonlinear elasticity[END_REF] followed by more related works [START_REF] Ozeré | Joint segmentation/registration model by shape alignment via weighted total variation minimization and nonlinear elasticity[END_REF][START_REF] Debroux | A joint segmentation/registration model based on a nonlocal characterization of weighted total variation and nonlocal shape descriptors[END_REF][START_REF] Debroux | A variational model dedicated to joint segmentation, registration, and atlas generation for shape analysis[END_REF] in which the proposed numerical method is adapted to the image registration problem in a nonlinear elasticity setting, we introduce auxiliary variables with quadratic penalty terms. Especially, this enables us to decouple the deformation Jacobian from the deformation itself and from its inverse. For the sake of conciseness, we denote by Tk =

F k = λ 2 Ω ( Tk • φ -Rk ) 2 dx + Ω W Op (V, detV ) dx +1 { . L ∞ (Ω,M 2 (R)) ≤α} (V ) + 1 { . L ∞ (Ω,M 2 (R)) ≤β} (W ) + γ1 2 V -∇φ 2 L 2 (Ω,M2(R)) + γ2 2 W -∇ψ 2 L 2 (Ω,M2(R)) + γ3 2 ζ -1 k-1 • φ -ϕ 2 L 2 (Ω,R 2 ) + γ4 2 ψ • φ -Id 2 L 2 (Ω,R 2 ) , (P k,d ) in which φ simulates ϕ 0 • ϕ 1 • . . . • ϕ k-1 • ϕ = ζ k-1 • ϕ, V mimics ∇φ while W simulates ∇(φ -1
), and ψ approximates φ -1 . To design the L ∞ -bound on W , we used the fact that if u is a homeomorphism of Ω into Ω, and the inverse function u -1 belongs to W 1,q (Ω), the matrix of weak derivatives reads ∇(u -1 ) = (∇u) -1 (u -1 ).

We then use a classical alternative scheme which consists of splitting the original problem into sub-problems that are more computationally tractable. The sketch of our numerical method of resolution is summarised in Algorithm 1.

We now present numerical simulations validating the accuracy of the multiscale representation of the obtained deformation.

Numerical Experiments

We tested our model on 4DMRI4 acquired during free-breathing of right lobe liver [START_REF] Von Siebenthal | 4D MR imaging of respiratory organ motion and its variability[END_REF]. We chose the images (195 × 166) corresponding to the liver in full exhalation and the liver in full inhalation to illustrate the capability of our model to handle large deformations. The parameters were chosen as follows : N = 10 levels, λ = 1, a 1 = 5.0, a 2 = 1000.0, a 3 = 4.0, γ 1 = 80000, γ 2 = 1.0, γ 3 = 1.0, γ 4 = 1.0, α = 100, and β = 100. The coefficients a 1 , a 2 and a 3 involved in the Ogden stored energy function affect respectively the averaged local change of length, and the averaged local change of area. This leads to the following phenomenon: the higher the a i 's are, the more rigid the deformation is. The γ i are chosen rather big to ensure the closeness between the auxiliary variables and the original ones. The λ weighs the fidelity term, and a trade-off has to be met between the accuracy of the registration with high values of λ and physically meaningful smooth deformations with small λ. The obtained results are illustrated in Fig. 1. We remark that for each level, the deformed Template is well-aligned with the Reference, demonstrating the ability of our model to deal with large deformations and localised small feature movements as the scale grows. The deformation grids do not exhibit overlaps and hence show that the obtained transformations are physically meaningful with a determinant remaining positive everywhere at all scales. By looking at the last column, one can see that the hierarchical decomposition of the deformations obtained with our model behaves as expected, that is to say the deformations on the first scales are global and represent the movements of the main organs and as k grows, the deformation becomes more localised and more refined to model the movements of small features, i.e. blood vessels imaged as small white dots.

Conclusion

To conclude, we have introduced a mathematically sound and physically relevant multiscale registration model inspired by the hierarchical decomposition presented in [START_REF] Tadmor | A Multiscale Image Representation Using Hierarchical (BV , L 2 ) Decompositions[END_REF]. It shows promising results on medical data, producing deformations from global coarse deformations to localised finer transformations

V 21,k ← proj { . L ∞ (Ω) ≤α} (V 21,k -(4 V k 2 V 21,k -2a 2 (detV k -1)
V at different levels. A theoretical analysis of our numerical scheme (asymptotic result) complemented by additional tests with adapting parameters constitute the next steps of our work, along with an extension to more involved fidelity terms and a discussion on the possible statistics unveiled by this hierarchical decomposition.

R

  j . The decoupled problem is thus defined by : inf ϕ,ψ,φ,V,W

Fig. 1 .

 1 Fig. 1. Multiscale registration results : each row represents a scale of the deformation; the first column displays the Template image at scale k, i.e. Tk , the second column shows the Reference image at scale k, i.e. Rk , the third one illustrates the deformed Template obtained at scale k, i.e. Tk •φ k , the fourth one exhibits the absolute difference | Tk • ϕ k -Rk | at scale k, the fifth column presents the inverse deformation at scale k, i.e. ψ k ≈ φ -1 k ≈ (ϕ0 • ϕ1 • . . . • ϕ k ) -1 , the sixth column represents the composition of deformations at scale k, i.e. φ k ≈ ϕ0 • ϕ1 • . . . • ϕ k , and finally the last column displays the deformation obtained at scale k, i.e. ϕ k .

  12,k + 10a 3 V 12,k detV 11 Lagrange equation with respect to ψ k using an L 2 gradient flow with implicit Euler time stepping ; , ψ k, V 11,k , V 12,k , V 21,k , V 22,k , W 11,k , W 12,k , W 21,k , W 22,k , ϕ k , Tk • φ k ;

		10a 3 V 11,k k detV 22	+ γ1(V 22,k -	∂φ 2,k
			∂ψ 2,k ∂x );
			∂ψ 2,k ∂y );
	17	for each pixel

k + γ1(V 21,k -∂φ 2,k ∂x ))); 12 V 22,k ← proj { . L ∞ (Ω) ≤α} (V 22,k -(4 V k 2 V 22,k +2a 2 (detV k -1) V 11,k -∂y ))); 13 W 11,k ← proj { . L ∞ (Ω) ≤β} ( ∂ψ 1,k ∂x ); 14 W 12,k ← proj { . L ∞ (Ω) ≤β} (

∂ψ

1,k ∂y ); 15 W 21,k ← proj { . L ∞ (Ω) ≤β} ( 16 W 22,k ← proj { . L ∞ (Ω) ≤β} ( 18 Solve the Euler -Lagrange equation with respect to φ k using an L 2 gradient flow with implicit Euler time stepping ; 19 Solve the Euler -20 ϕ k ← ζ -1 k-1 • φ k ; 21 return φ k

http://www.vision.ee.ethz.ch/ organmot/chapter download.shtml

Algorithm 1 Our Proposed Method (L ∞ constraints applied componentwise) Rj ;

7 for l = 1, . . . , nbIter:

for each pixel :

V 22,k -

V 21,k +