Kevin Loquin
email: kevin.loquin@telecom-paristech.fr

Isabelle Bloch
email: isabelle.bloch@telecom-paristech.fr

Kiyoko Nakashima

Florence Rossant
email: florence.rossant@isep.fr

Michel Paques

Photoreceptor detection in in-vivo Adaptive Optics images of the retina: Towards a simple interactive tool for the physicians

Keywords: Adaptive Optics, Photoreceptor detection, in vivo diagnosis, retina imaging

INTRODUCTION

The light-sensitive tissue of the eye is the retina, which covers the eye fundus. At a microscopic level, the retina is a stack of several neuronal layers whose aim is to transform the incident light in an electrical signal (phototransduction) that will be transmitted to the brain (Figure 1) [START_REF] Polyak | The Retina[END_REF][START_REF] Curcio | Human photoreceptor topography[END_REF]. The photosensitive neurons are the photoreceptor cells (although some ganglion cells are responsive to light). These are of two types: rods and cones. Rods function mainly in dim light and provide black-and-white vision, while cones support daytime vision and the perception of colour. Many inherited and acquired diseases or disorders may affect the retina [3,[START_REF] De | Age-related macular degeneration[END_REF]. Some of these provoke progressive degeneration of the retina. While the loss of rods is usually well tolerated in everyday life, loss of cones leads to severe visual handicap. In some diseases such as retinitis pigmentosa, visual loss can occur very early in life. Yet, in most cases, visual acuity usually deteriorates gradually over decades and may end up to mere light perception.

In recent years, it became possible to image the living human retina at the microscopic scale by means of new imaging systems: Adaptive Optics (AO) [START_REF] Roddier | Adaptive Optics in Astronomy[END_REF][START_REF] Roorda | Adaptive optics scanning laser ophthalmoscopy[END_REF]. The AO images presented in this article are obtained with a prototypic AO system (rtx1 Adaptive Optics Retinal Camera manufactured by Imagine Eyes, Orsay, France), which is currently in operation in a clinical setting at the Clinical Investigation Center 503 of the XV -XX hospital. Each acquisition is a sequence of forty 1392 × 1040 en-face images. In each image, considering a mean axial length of 24.7mm, one pixel corresponds to a 1.6µm × 1.6µm spot on the retina. The images that we are dealing with are shown in Figure 2. They are the result of a registration/summation algorithm applied to the forty images acquired successively, in order to increase the signal-to-noise ratio.

The obtained AO images are very useful to a physician. This remark is illustrated in Figure 2, where a non degenerative (image a) and degenerative cases (image b, c, d) can easily be differentiated. In [START_REF] Wolfing | High-resolution retinal imaging of cone-rod dystrophy[END_REF], the authors illustrate the diagnosis power of retinal AO images on degenerative cones distrophy by comparing manual diagnosis procedures from these AO images with other usual diagnosis procedures.

Automated retinal AO image analysis methods for assisting in vivo diagnosis and follow up of early degenerative cones dystrophy are essential for physicians. In this paper, we propose a novel algorithm that aims at detecting the photoreceptors in the AO images (Figure 2). The cornerstone of all the AO images based diagnosis procedures presented in [START_REF] Wolfing | High-resolution retinal imaging of cone-rod dystrophy[END_REF] is the detection of the cones in these images. From these detected cones, we can derive photoreceptor density maps or extract other statistical features useful for the diagnosis.

Automatic photoreceptors detection algorithms have already been proposed in the literature. In [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF], the authors proposed a multistep algorithm. First, they propose to work with an above-threshold version of the image in order to eliminate the dim spots. Afterwards, the regional maxima are detected from a low-pass filtered version of the thresholded image in order to eliminate the noise. Finally, a morphological dilation is used to merge close enough detected spots. The structuring element is chosen according to the minimal cone spacing. In [START_REF] Xue | Photoreceptor counting and montaging of enface retinal images from an adaptive optics fundus camera[END_REF], the same kind of procedure is proposed but the regional maxima detection is performed by decreasing greylevel regions in order to favor the brightest detected spots and avoid the threshold step of the method in [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF].

We distance these articles by providing a new low-level image processing algorithm which is computationally very efficient. The geometrical arguments from which we construct our algorithm are the same as the ones used in [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF][START_REF] Xue | Photoreceptor counting and montaging of enface retinal images from an adaptive optics fundus camera[END_REF]: elimination of dim spots and merging of too close spots. Our algorithm is simpler and more computationally efficient because directly constructed from the image geometry and these geometric arguments. It is based on a recursive construction of thresholded connected components when the seeds of the recursions are the regional maxima of the image. The stopping rules of this recursive construction are deduced from the previous geometrical arguments and the finally detected photoreceptors are the highest greylevel spots in a merged maximum area (a thresholded connected component).

The efficiency of our algorithm and its simplicity allows a visual and adaptive control of the algorithm parameterization, which makes it easily tunable by non computer scientists like physicians. An important facet of our work stands in the usability of our tool, implemented as an ImageJ plugin in collaboration with physicians.

In Section 2, we present and geometrically justify our photoreceptor detection algorithm. Then, Section 3 details the main user interactions provided with our implemented tool. Finally, in Section 4, we illustrate the results of our method when applied to various images and propose a simple visual comparison with the method proposed in [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF], before concluding.

PHOTORECEPTOR DETECTION ALGORITHM

In order to illustrate the algorithm functioning, we accompany the following discussion about the algorithm with Figure 3 which shows the algorithm result for five different tolerance parameter (explained below) choices on a unique 1-dimensional grey level signal and with Figure 4 which shows the result of our labeling algorithm when applied to a medical AO image of the retina.

In the used AO images, the photoreceptors that we aim at detecting are bright spots, i.e. pixels with high grey levels. Hence a first natural step of our algorithm is a simple regional maxima detection. A regional maximum is a plateau of the grey level function, i.e. a connected component of a level α in which every pixel has a grey level equal to α.This terminology is reused from existing works on the component trees found in the literature [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF]. These regional maxima appear as RM on the top-left diagram of Figure 3. This raw regional maxima processing is not enough to qualify these regional maxima as the photoreceptors that we aim at detecting. Indeed, a regional maximum with a low greylevel (i.e. a dim spot) is not generally a photoreceptor. The thresholding of the AO image is the solution proposed in [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF] to eliminate these spots. Besides, too close regional maxima can be due to a unique photoreceptor. This idea led the authors of [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF][START_REF] Xue | Photoreceptor counting and montaging of enface retinal images from an adaptive optics fundus camera[END_REF] to morphologically filter the extracted regional maxima. In our approach, we define a geometrical concept of influence area of a maximum that we will call a maximum area. The construction of these maximum areas induces the elimination of the undesired spots. These maximum areas appear as colored rectangles in Figure 3 for different tolerances.

A maximum area of a regional maximum is defined up to a tolerance parameter, denoted by T . This tolerance parameter finds its roots in the component tree decomposition of an image [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF]. In this fructuous research field, an image is seen as a decomposition of αlevel connected components. The level range corresponding to the greyscale range of the image, an image is decomposed in a stack of connected components. The tolerance parameter T that we are using in our algorithm tunes the height of the top part of this stack that we consider as a maximum area. The top parts of the stacks appear as colored bars in Figure 3. When the color bar (or the tolerance top level) intersects a pixel, this pixel is incorporated in the maximum area of the associated regional maximum.

It is trivial to remark that the higher the tolerance value, the bigger the maximum areas of the regional maxima. Going further, for a given T , two or more regional maxima can lead to the same maximum area. In Figure 3, when the tolerance is T = 40 we can observe the fusion of Ph1 and Ph2 into one photoreceptor: Ph1 which is of higher grey level. This is the first criterion that we use to eliminate regional maxima of smallest grey levels, seeding the same fused maximum area.

More technically, our algorithm constructs the connected component rooted by all the regional maxima, taken in a decreasing greylevel order. For each regional maximum, this construction is recursive around the considered regional maximum and the recursivity is stopped if we reach a pixel • whose greylevel is out of the top stack of the image, i.e. out of the connected component defined from the considered regional maximum by means of the tolerance; or • already processed from another regional maximum: it means that we reach a connected component of higher greylevel regional maximum (the fusion of too close spots is performed here by keeping only the brightest spot in a connected component. Besides a filtering of the dimer spots is implicitely performed); or • whose greylevel is higher than the greylevel of the initial regional maximum: it occurs when the reached higher greylevel pixel was not itself a regional maximum because one of its direct neighbors is higher (Once again, an implicit filtering of the dim spots is performed here).

More than filtering the roughly detected regional maxima, the implemented algorithm labels the pixels since it is followed by the application of a marker-controlled watershed segmentation algorithm. Therefore, at each pixel is associated one of the following labels:

• M AX_P OIN T : a finally detected photoreceptor,

• EQU AL: inside a regional maximum,

• M AX_AREA: inside a maximal area.

Regional Maxima MAX_AREA EQUAL MAX_POINT Fig. 4. Pixel labeling illustration.

Note that the proposed algorithm is a truncated adaptation of the watershed by immersion algorithm [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[END_REF] when applied to the inverted image. In the inverted image, the origins of the basins used in the watershed by immersion are the regional maxima of the original image.

In our case, we are not interested in reaching a complete separation of the basins (the maxima) but only in elicitating these basins until a given height, i.e. the maximum areas up to a tolerance parameter).

ALGORITHM IMPLEMENTATION: AN INTERACTIVE TOOL FOR THE PHYSICIANS

This algorithm is very rapid: after a complete cover of the image for finding regional maxima, the recursive coverage is only applied to the extracted regional maxima which highly diminishes the computational cost. When java-implemented as an imageJ plugin (which is not the most rapid computational framework) the algorithm result is instantly displayable. This remark is very important because it allows us to provide an interactive tool to the final user, the physician. This non algorithmic part of our work is greatly appreciated by the physicians who can thus easily interact with the algorithm and adapt its parameters to his visual expertise in an intuitive way. This feature highly contributes to the fact that, even if still in development, this tool is already in use in a clinical context. Among the proposed interactions, the parameterization (more precisely the tolerance value choice) is the most important. The highest the tolerance value, the smallest the number of detected cones. As illustrated in Figure 2, the images are highly different according to the involved disease or the acquisition conditions. It does not seem possible to find an automatic robust optimal tolerance selection method, but it is still one of our research interests. Thus, in the tool that we develop and implement in collaboration with the physicians, we proposed the solution of incorporating a scrollbar to tune the tolerance parameter. Each time the user moves this scrollbar, the detected photoreceptors directly appear on a refreshed image. Thus, the physician is the judge of the best trade-off between wrong detected cones and undetected cones.

Another way to interact on the photoreceptors detection is to incorporate in our imageJ plugin the possibility to add non automatically detected cones and to remove wrongly detected cones by the hand. Note that in a following work, we propose to perform a statistical study of these two indices: the number of added cones and the number of removed cones in order to evaluate our detection algorithm. We assume that these indices quantitatively reflect the subjective judgement of the physician about our algorithm.

ILLUSTRATIVE EXPERIMENTS

In order to illustrate the results of our algorithm on AO images of the retina, we work with the images of Figure 2. All of them have the same size: 300 × 300 (i.e. 480µm × 480µm). The first column in Figure 5 presents the photoreceptors automatically detected with our algorithm for respective tolerance of T = 20, T = 30, T = 30 and T = 34. The tolerance has been interactively chosen by the authors with an aim of visualizing the results at best. No manual removal or addition of photoreceptors have been made here. The second column presents the results obtained with the algorithm of [START_REF] Li | Automated identification of cone photoreceptors in adaptive optics retinal images[END_REF]. The parameters are fc: the cutoff frequency of the used lowpass filter and thr is the chosen threshold. For proper parameterizations, the obtained results are very similar for these methods. Nevertheless, with our tool, the physicians have a less esoteric parameterization procedure that leads to more usability.

Figure 6 presents the results that we obtain with the markercontrolled watershed segmentation results using the same tolerance parameters as in the previous experiment. The segmentation result is very useful to physicians. First it exhibits parts of low cone density when the segmented parts are big (see Figure 6 b, c andd). In other words, the segmentation is a representation of information dual to a cone density map. Besides, the segmentation result provides information about the cone mosaic configuration (e.g. hexagonal) in the high cone density parts (see Figure 6 a). The extracted quantitative topographical information is of primary importance especially in the context of longitudinal follow-ups of cone degenerative distrophies. This is a real advantage for the physician to accompany his visual judgment with a quantization of the disease evolution.

CONCLUSION

In this paper, we proposed a new photoreceptor detection algorithm applied to in-vivo AO images of the retina. This algorithm is based on a recursive construction of thresholded connected components when the seeds of the recursions are the regional maxima of the image. Our method is implemented as an imageJ plugin already in use in a clinical setting at the XV -XX hospital. This use shows its early clinical value before any objective evaluation. In a raw comparative study, our algorithm seems to provide at least as good results (in terms of detection efficiency) as already existing approaches. The main advantage of our approach stands more in its computational efficiency, which is highly appreciated by the physicians especially because we provide an interactive parameterization, allowing an optimal (according to the physician expertise) detection. Note that a semi-subjective evaluation protocol is currently in progress with the XV -XX hospital: a statistical study of correction indices derived from the number of manually removed and added cones is established.

Fig. 1 .

 1 Fig. 1. Schematic eye and retina. This figure is based on materials obtained from http://www.mrawde.com/blog/2008/05/readingin-the-dark.html and http://webvision.med.utah.edu/sretina.html.

Fig. 2 .

 2 Fig. 2. AO images of a: a healthy eye. b: Retinitis Pigmentosa (RP) c: Age-related Macular Distrophy (AMD) due to Drusen. d: Occult Macular Distrophy (OMD).

Fig. 3 .

 3 Fig. 3. Illustration of the algorithm.

Fig. 5 .

 5 Fig. 5. Comparative experiment: our algorithm (left column) vs. Li and Roorda algorithm [8] (right column) results when applied to the AO images of Figure 2.

Fig. 6 .

 6 Fig. 6. Marker-controlled watershed segmentation results when applied to the AO images of Figure 2.

* This work has been supported by French National Research Agency (ANR) through Tecsan Program (project iPhot n o ANR-09-TECS-009).


Adrien Le Coënt
email: adrien@cs.aau.dk

˚˚D Fribourg ˚˚jonathan Vacher
email: fribourg@lsv.ens-cachan.fr

˚˚˚a Lbert

Control Synthesis for Stochastic Switched Systems using the Tamed Euler Method

Keywords: Stochastic systems, numerical simulation, control system synthesis, switched control systems, nonlinear control systems

In this paper, we explain how, under the one-sided Lipschitz (OSL) hypothesis, one can find an error bound for a variant of the Euler-Maruyama approximation method for stochastic switched systems. We then explain how this bound can be used to control stochastic switched switched system in order to stabilize them in a given region. The method is illustrated on several examples of the literature.

INTRODUCTION

Symbolic methods for the verification and control synthesis of hybrid systems (and, particularly, "switched systems") have received significant attention in the past few years.

One distinguishes two main classes of symbolic methods for hybrid systems: indirect methods and direct methods [START_REF] Asarin | Effective synthesis of switching controllers for linear systems[END_REF].

Indirect methods proceed by constructing a finite abstraction of the original system by discretization of the dense state space R n (where n is the dimension of the state space). Among the indirect methods, one of the most successful proceeds by approximate bisimulation [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF]. This method originally designed for deterministic switched systems has been recently extended for stochastic switched systems [START_REF] Zamani | Symbolic models for stochastic switched systems: A discretization and a discretization-free approach[END_REF][START_REF] Zamani | Symbolic control of stochastic systems via approximately bisimilar finite abstractions[END_REF][START_REF] Zamani | Towards scalable synthesis of stochastic control systems[END_REF]. This approach relies on the hypothesis of incremental stability of the stochastic switched system (or existence of a common/multiple Lyapunov function).

A direct method proceeds by working directly at the level of the dense state system R n ; it computes "trajectory tubes", which are over-approximations of the set of all the controlled trajectories starting at a given subregion of R n . In previous work, We have followed such a direct approach (eg, [START_REF] Fribourg | Finite controlled invariants for sampled switched systems[END_REF]). The idea is to start with two given hyperrectangles R and S of R n , (with R Ď S): one covers R with a finite number of subregions (of the form of balls or sub-rectangles), and finds by exhaustive search, for each subregion, a "control pattern" (i.e., a finite control sequence) such that the trajectories starting from the subregion and controlled by the pattern goes back to R while never leaving S. Such a direct method ensures the so-called property of "pR, Sq-stability". We have recently applied such a direct method in the deterministic framework, using the Euler approximation scheme for calculating over-approximations of tubes of trajectories Le [START_REF] Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF]. We show here how to extend this direct method in order to treat stochastic switched systems. The method is a simple extension of the deterministic method, but replaces the classical Euler approximation scheme, by a variant of the stochastic Euler-Maruyama scheme [START_REF] Hutzenthaler | Strong convergence of an explicit numerical method for SDEs with nonglobally lipschitz continuous coefficients[END_REF]. The correctness of these Euler-based methods does not rely on the hypothesis of incremental stability as in [START_REF] Zamani | Symbolic models for stochastic switched systems: A discretization and a discretization-free approach[END_REF][START_REF] Zamani | Towards scalable synthesis of stochastic control systems[END_REF], but on the hypothesis of 'onesided Lipschitz (OSL)' condition with constant λ P R d (also called 'monotonicity'/'dissipativity', see von Renesse and Scheutzow (2010)). It can be seen that if a stochastic switched system satisfies an OSL condition with λ ă 0, then the function V px, x 1 q " }x ´x1 } 2 is a common incremental Lyapunov function in the sense of [START_REF] Zamani | Symbolic control of stochastic systems via approximately bisimilar finite abstractions[END_REF], from which it follows that the switched system is incrementally stable, and can be treated by approximate bisimulation. However, Euler-based methods also apply when the system is not incrementally stable, in which case the constant λ is necessarily positive. We thus consider a class of systems different from that of [START_REF] Zamani | Symbolic control of stochastic systems via approximately bisimilar finite abstractions[END_REF].

The plan of the paper is as follows: In Section 2, we give an explicit upper bound on the mean square error of the tamed Euler method for SDEs under OSL condition. We apply the result in order to ensure properties of stochastic switched systems, such as "pR, Sq-stability" (Section 3). We conclude in Section 4.

BOUNDING THE ERROR OF THE TAMED EULER METHOD

Assumptions

The symbol } ¨} denotes the Euclidean norm on R d .The symbol x¨, ¨y denotes the scalar product of two vectors of R d . Given a point x P R d and a positive real r ą 0, the ball Bpx, rq of centre x and radius r is the set ty P R d | }x ý} ď ru.

Let τ P p0, 8q be a fixed real number, let pΩ, F, Pq be a probability space with normal filtration pF t q tPr0,τ s , let d, m P N :" t1, 2, . . . u let W " pW p1q , . . . , W pmq q : r0, RsˆΩ Ñ R m be an m-dimensional standard pW t q tPr0,τ s -Brownian motion and let x 0 : Ω Ñ R d be an F 0 {BpR d qmeasurable mapping with Er}x 0 } p s ă 8 for all p P r1, 8q. Moreover, let f : R d Ñ R d be a continuously differentiable and globally one-sided Lipschitz continuous function whose derivative grows at most polynomially and let g " pg i,j q iPt1,...,du,jPt1,...,mu : R d Ñ R dˆm be a globally Lipschitz continuous function.

Then consider the Stochastic Differential Equations (SDE): dX t " f pX t qdt `gpX t qdW t , X 0 " x 0 (1) for t P r0, τ s. The drift coefficient f is the infinitesimal mean of the process X and the diffusion coefficient g is the infinitesimal standard deviation of the process X. Under the above assumptions, the SDE (1) is known to have a unique strong solution. More formally, there exists an adapted stochastic process X : r0, τ s ˆΩ Ñ R d with continuous sample paths fulfilling X t,x0 " x 0 `ż t 0 f pX s qds `ż t 0 gpX s qdW s for all t P r0, τ s P-a.s. (see, e.g., [START_REF] Oksendal | Stochastic Differential Equations: An Introduction with Applications[END_REF]).

We denote by X t,x0 the solution of Equation (1) at time t from initial condition X 0,x0 " x 0 P-a.s., in which x 0 is a random variable that is measurable in F 0 .

We suppose that f behaves polynomially and g is Lipschitz, i.e.: there exist constants D P R ě0 , q P N and L g P R ě0 such that, for all x, y P R d }f pxq ´f pyq} 2 ď D}x ´y} 2 p1 `}x} q `}y} q q (H1)

}gpxq ´gpyq} ď L g }x ´y} (H2)
We also assume that the SDE (1) satisfies the following one-sided Lipschitz (OSL) condition with constant λ P R:

Dλ P R @x, y P R d : xf pyq´f pxq, y ´xy ď λ }y ´x} 2 (H3) Remark 1. Constants λ, L g and D can be computed using (constrained) optimization algorithms (see Le [START_REF] Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF]).

Tamed Euler approximation

The standard way to extend the classical Euler method for ordinary differential equations to the SDE (1) is the Euler-Maruyama scheme [START_REF] Maruyama | Continuous Markov processes and stochastic equations[END_REF]. More precisely, given z : Ω Ñ R d an F 0 {BpR d q-measurable mapping with Er}z} p s ă 8 for all p P r1, 8q, the explicit Euler-Maruyama (EM) method for the SDE (1) is given by the mappings

Y N n,z : Ω Ñ R d , n P t0, 1, . . . , N u, which satisfy Y N 0,z " z and Y N n`1,z " Y N n,z `τ N ¨f pY N n,z q `gpY N
n,z qpW pn`1qτ {N ´Wnτ{N q for all n P t0, 1, . . . , N ´1u and all N P N. See [START_REF] Maruyama | Continuous Markov processes and stochastic equations[END_REF]. Unfortunately, the convergence results for the EM scheme does not hold when the drift function f of the SDE (1) behaves polynomially (and not linearly). For the sake of generality, we will now adopt a refined scheme, which has been proposed recently in order to overcome this difficulty [START_REF] Hutzenthaler | Strong convergence of an explicit numerical method for SDEs with nonglobally lipschitz continuous coefficients[END_REF]

. Let X N n,z : Ω Ñ R d , X N n`1,z " X N n,z `τ N ¨f pX N n,z q 1 `τ N ¨}f pX N n,z q} `gpX N n,z qpW pn`1qτ N ´W nτ N q
(2) for all n P t0, 1, . . . , N ´1u and all N P N. We refer to the numerical method (2) as a tamed Euler scheme [START_REF] Hutzenthaler | Strong convergence of an explicit numerical method for SDEs with nonglobally lipschitz continuous coefficients[END_REF]. In this method the drift term τ n ¨f pX N n,z q is "tamed" by the factor 1{p1 `τ N ¨}f pX N n,z q}q for n P t0, 1, . . . , N ´1u and N P N in (2).

A time continuous interpolation of the time discrete numerical approximations ( 2) is also introduced in Hutzenthaler et al. ( 2012) as follows. Let XN z : r0, τ s ˆΩ Ñ R d , N P N, be a sequence of stochastic processes given by

XN t,z " XN n,z `pt ´nτ {N q ¨f p XN n,z q 1 `τ {N ¨}f p XN n,z q} `gp XN n,z qpW t ´W nτ N q
for all t P r nτ N , pn`1qτ N s, n P t0, 1 . . . , N ´1u and all N P N. Note that XN t,z : r0, τ s ˆΩ Ñ R d is an adapted stochastic process with continuous sample paths for every N P N.

Let us define X N t,z by

X N t,z :" X N n,z for t P r nτ N , pn `1qτ N q.
Note that XN t,z " X N t,z " X N n,z at time t " nτ N for n P t0, 1, . . . , N u.

The following theorem is proven in [START_REF] Hutzenthaler | Strong convergence of an explicit numerical method for SDEs with nonglobally lipschitz continuous coefficients[END_REF]: Theorem 1. [START_REF] Hutzenthaler | Strong convergence of an explicit numerical method for SDEs with nonglobally lipschitz continuous coefficients[END_REF]). Let us suppose (H1) (H2) and (H3). Let the setting in this section be fulfilled, and z : Ω Ñ R d be an F 0 {BpR d q-measurable mapping with Er}z} p s ă 8 for all p P r1, 8q. Then, for all p P r1, 8q

sup N PN sup nPt0,1,...,N u Er}X N n,z } p s ă 8
For the sake of simplicity, the number N of subsampling steps is now left implicit. From Theorem 1, it follows (cf. Lemma 4.3, [START_REF] Higham | Strong convergence of Euler-type methods for nonlinear stochastic differential equations[END_REF]): Lemma 1. Let us suppose (H1) (H2) and (H3). Let the setting in this section be fulfilled, and z : Ω Ñ R d be an F 0 {BpR d q-measurable mapping with Er}z} p s ă 8 for all p P r1, 8q. Then, for any even integer r ě 2, there exist two constants E r,z and F r,z such that sup 0ďtďτ E}X t,z ´X t,z } r ď p∆ t q r 2 pE r,z p∆ t q r 2 `Fr,z dq.

with ∆ t " τ {N and:

E r,z " 2 r p}f p0q} r `D2 r`1 2 p1 `E sup 0ďtďτ }X t,z } qr q 1 2 pE sup 0ďtďτ }X t,z } 2r q 1 2 q, F r,z " 2 r p}gp0q} 2r `Lr g E sup 0ďtďτ }X t,z } r 2 q.
Proof. see Appendix.

Remark 2. Constants E r,z and F r,z are computed using the constants λ and L g (see Remark 1), and the expected values of X t,z at each time t " 0, ∆t, 2∆t, . . . , N ∆t. These expected values are computed using a Monte Carlo method (by averaging here the value of 10 4 samplings).

Mean square error bounding

The following Theorem holds for SDE (1). This corresponds to a stochastic version of Theorem 1 of Le Er}x 0 ´z} 2 s ď δ 2 0 . Then, we have, for all τ ě 0:

Er sup 0ďtďτ }X t,x0 ´X t,z } 2 s ď δ 2 τ,δ0 ,
with δ 2 τ,δ0 :" βpτ qe γτ , where:

γ " 2p ? ∆ t `2λ `L2 g `128L 4 g q, and βpτ q " 2δ 2 0 `2τ

∆ t L 2 g p1 `128L 2 g qpF 2,z d `E2,z ∆ t q `4τ a ∆ t DpF 4,z d `E4,z ∆ 2 t q 1 2 p1 `4E sup 0ďtďτ }X t,z } 2q `4E sup 0ďtďτ } Xt,z } 2q q 1 2 .
(3) with ∆ t " τ {N .

Proof. The proof closely follows the proof of Theorem 4.4 in [START_REF] Higham | Strong convergence of Euler-type methods for nonlinear stochastic differential equations[END_REF]. Let e t " X t,x0 ´X t,z . We have, for all 0 ď t ď τ : de t " pf pX t,x0 q ´f pzqqdt `pgpX t,x0 q ´gpzqqdW t . (4) Then, by using Equation ( 4) and the integral version of Itô formula applied to function x Þ Ñ }x} 2 we obtain

}e t } 2 " }e 0 } 2 `ż t 0 2xe s , f pX s,x0 q ´f pX s,z qyds `ż t 0 }gpX s,x0 q ´gpX s,z q} 2 ds `M ptq, (5)
where e 0 " x 0 ´z, and M ptq " ż t 0 2xe s , gpX s,x0 q ´gpX s,z qydW s .

So we have using (H2):

}e t } 2 ď }e 0 } 2 `ż t 0 2xe s , f pX s,x0 q ´f p Xs,z qyds `L2 g ż t 0 }X s,x0 ´Xs,z } 2 ds `ż t 0
2xe s , f p Xs,z q ´f pX s,z qyds `M ptq.

(6) So we have using (H3) and Young's inequality:

}e t } 2 ď }e 0 } 2 `ż t 0 p2λ}e s } 2 `L2 g }e s } 2 qds `L2 g ż t 0 }X s,z ´X s,z } 2 ds `ż t 0 p 1 ? ∆ t }f p Xs,z q ´f pX s,z q} 2 `a∆ t }e s } 2 qds `M ptq. (7)
So we have using (H1), for all 0 ď t ď τ :

}e t } 2 ď }e 0 } 2 `pa ∆ t `2λ `L2 g q ż t 0 }e s } 2 ds `L2 g ż t 0 }X s,z ´X s,z } 2 ds `D ? ∆ t ż t 0 p1 `}X s,z } q `} Xs,z } q q}X s,z ´X s,z } 2 ds `M ptq. (8)
It follows using Lemma 1 for r " 2, and Cauchy-Schwarz inequality:

Er sup 0ďsďt }e s } 2 s ď E}e 0 } 2 `pa ∆ t `2λ `L2 g q ż t 0 E}e s } 2 ds `L2 g τ ∆ t pE 2,z ∆ t `F2,z dq `D ? ∆ t ż t 0 pEp1 `}X s,z } q `} Xs,z } q q 2 q 1 2 pE}X s,z ´X s,z } 4 q 1 2 ds `mptq, (9)
where mptq " Er sup 0ďsďt }M psq}s.

Hence, using using Lemma 1 for r " 4, and inequality pa `bq r ď 2 r pa r `br q:

Er sup

0ďsďt }e s } 2 s ď E}e 0 } 2 `pa ∆ t `2λ `L2 g qq ż t 0 E}e s } 2 ds `L2 g τ ∆ t pE 2,z ∆ t `F2,z dq `2Dτ a ∆ t pE 4,z ∆ 2 t `F4,z dq 1 2 p1 `4E sup 0ďtďτ }X t,z } 2q `4E sup 0ďtďτ } Xt,z } 2q q 1 2
`mptq.

(10)

On the other hand, from the Burkholder-Davis-Gundy inequality, we get:

mptq ď 16Er ż t 0 }e s } 2 }gpX s,x0 q ´gpX s,z q} 2 dss 1 2
Hence, using (H2):

mptq ď 16L 2 g Er sup 0ďsďt }e s } 2 ż t 0 }X s,x0 ´Xs,z } 2 dss 1 2
Then, using Young's inequality (for any α ą 0):

mptq ď 8L 2 g pαEr sup 0ďsďt }e s } 2 s `1 α Er ż t 0 }X s,x0 ´Xs,z } 2 dssq.
Hence, by using Lemma 1 for r " 2:

mptq ď 8αL 2 g Er sup 0ďsďt }e s } 2 s `8L 2 g α ż t 0 Er sup 0ďrďs }e r } 2 sds `8L 2 g α τ ∆ t pE 2,z ∆ t `F2,z dq. (11)
Hence, letting α "1 16L 2 g , we have by replacing in (10):

1 2 Er sup 0ďsďt }e s } 2 s ď δ 2 0 `pa ∆ t `2λ `L2 g `128L 4 g q ż t 0 Er sup 0ďrďs }e r } 2 sds `τ pL 2 g `128L 4 g q∆ t pE 2,z ∆ t `F2,z dq `τ 2D a ∆ t pE 4,z ∆ 2 t `F4,z dq 1 2 p1 `4E sup 0ďtďτ }X t,z } 2q `4E sup 0ďtďτ } Xt,z } 2q q 1 2 . ( 12
)
It results from Gronwall's inequality:

Er sup 0ďtďτ }e t } 2 s " βpτ qe γτ , with γ " 2p ? ∆ t `2λ `L2 g `128L 4 g q, and βpτ q " 2δ 2 0

`2τ p∆ t L 2 g p1 `128L 2 g qpF 2,z d `E2,z ∆ t q `4τ a ∆ t DpF 4,z d `E4,z ∆ 2 t q 1 2 p1 `4E sup 0ďtďτ }X t,z } 2q `4E sup 0ďtďτ } Xt,z } 2q q 1 2 . ( 13
)
It follows from Theorem 2 and Jensen's inequality: Proposition 1. Consider two points x 0 and z of R d ,and a positive real number δ 0 . Suppose that x 0 P Bpz, δ 0 q (i.e. }x 0 ´z} ď δ 0 ). Then EX t,x0 P Bp Xt,z , δ t,δ0 q for all t P r0, τ s.

It also follows from Theorem 2: Proposition 2. In the setting of Theorem 2, the expression δ τ,δ0 tends to δ 0 ? 2e 2λτ `L2

g `128L 4

g when ∆ t tends to 0 (i.e., when N tends to 8).

Implementation

This method has been implemented in the interpreted language Octave, and the experiments performed on a

-5.5 -5 -4.5 -4 -3.5 -3 -2.5 -2 -6 -4 -2 0 2 4 X1 X2
Fig. 1. Example 1 with mode u " 1, τ " 1, ∆ t " 10 ´4, initial ball Bpz, δ 0 q with z " p´4, 3.8q and δ 0 " 0.5, final ball Bpz 1 , δ τ,δ0 q with z 1 " p´3.6, 2.56q and δ τ,δ0 " 1.17

2.80 GHz Intel Core i7-4810MQ CPU with 8 GB of memory. The implementation is an adaptation of the program described in Le [START_REF] Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF] for controlling deterministic switched systems, but makes use of the tamed Euler scheme for SDEs (with the error function δ given in Theorem 2) instead of the classical Euler scheme.

Example 1. Consider the following system, corresponding to the example in Section 6.2 of [START_REF] Zamani | Symbolic models for stochastic switched systems: A discretization and a discretization-free approach[END_REF] (cf. [START_REF] Zamani | Symbolic control of stochastic systems via approximately bisimilar finite abstractions[END_REF]) for mode u " 1:

dx 1 " p´0.25x 1 `x2 `0.25qdt `0.05x 1 dW 1 t dx 2 " p´2x 1 ´0.25x 2 ´2qdt `0.05x 2 dW 2 t
The program gives (for τ " 1, ∆ t " τ {10 4 ): q " 0, D " 1.36, L g " 0.05, λ " 0.25; and for z " p´4, ´3.8q: E 2,z " 893.3, E 4,z " 2.14 ¨10 5 , F 2,z " 0.002, F 4,z " 4.9 10

´6.

Consider now the system corresponding to the example of [START_REF] Zamani | Symbolic models for stochastic switched systems: A discretization and a discretization-free approach[END_REF] for mode u " 2:

dx 1 " p´0.25x 1 `2x 2 ´0.25qdt `0.05x 1 dW 1 t dx 2 " p´x 1 ´0.25x 2 `1qdt `0.05x 2 dW 2 t
The program gives (for τ " 1, ∆ t " τ {10 4 ): q " 0, D " 1.36, L g " 0.05, λ " 0.25, and, for z " p0, 3q: E 2,z " 543.2, E 4,z " 7.94 ¨10 4 , F 2,z " 0.0442, F 4,z " 0.00178.

Both computations take less than 10 s. of CPU time. Simulations of the two systems are given in Figure 1 for mode u " 1 and starting point z " p´4, 3.8q, and in Figure 2 for mode u " 2 and starting point z " p0, 3q. On each figure, the initial ball (t " 0) is depicted in black, the final ball (t " τ ) in red, and 200 random sampling trajectories in blue for t P r0, τ s. 1

0 1 2 3 4 5 6 -1 0 1 2 3 4 X1 X2
Fig. 2. Example 1 with mode u " 2, τ " 1, ∆ t " 10 ´4, initial ball Bpz, δ 0 q with z " p0, 3q and δ 0 " 0.5, final ball Bpz 1 , δ τ,δ0 q with z 1 " p0.79, ´0.63q and δ τ,δ0 " 1.17

SAMPLED STOCHASTIC SWITCHED SYSTEMS

Stochastic switched system as a finite collection of SDEs

We now consider a finite number of SDEs. Each SDE is referred to as a mode j, and the set of modes is referred to as U " t1, . . . , M u. We will denote by X j t,x0 the solution at time t of the system: dxptq " f j pxptqq `gj pxptqqdW j t , xp0q " x 0 .

(

) 14
where x 0 is a random variable that is measurable in F 0 . Hypotheses (H1-H2-H3), as defined in Section 2, are naturally extended to every mode j of U . Accordingly, constants L g , λ, F associated to SDE (1) in Section 2, now become L gj , λ j , F j respectively, for each j P U .

Likewise, for each j P U , the nonnegative real pδ t,δ0 q2 becomes pδ j t,δ0 q 2 for each mode j; the approximate continuous-time solution of ( 14) starting from z, is denoted by Xj t,z , and the approximate staircase solution by X j t,z .

Control patterns

The control laws that we now consider are "piecewise constant of duration τ " in the sense that, every τ seconds, they select a given mode (see [START_REF] Zamani | Symbolic models for stochastic switched systems: A discretization and a discretization-free approach[END_REF]). We call "(control) pattern of length k" a sequence of k modes (i.e., an element of U k ). Each pattern π of the form j 1 j 2 ¨¨¨j k corresponds to the selection of mode j 1 for time t P r0, τ q, then mode j 2 for t P rτ, 2τ q, and so on, until t " kτ . We assume that the solution of the system is continuous at sampling instants t " τ, 2τ, . . . (which means that there is no "reset" of the system at sampling instants).

Given a stochastic switched system, a pattern π of length k and an initial random variable z, one constructs the "approximate solution controlled by π" by composing together the approximations obtained by successive application of the modes of π. Formally, the "continuous" approximate solution Xπ t,z is defined at time t P r0, kτ s as follows:

' Xπ t,z " Xj t,z if π " j P U , k " 1 and t P r0, τ s, and ' Xπ pk´1qτ `t1 ,z " Xj t,z 1 with z 1 " Xπ 1 pk´1qτ,z if k ě 2, t 1 P r0, τ s, π " π 1 ˚j for some j P U and π 1 P U k´1 .

The "staircase" approximate solution X π t,z is defined analogously. Likewise, given an initial error radius δ 0 ą 0 and a pattern π of length k ě 1, one defines the error radius δ π t,δ0 as follows: 2 ' δ π t,δ0 " δ j t,δ0 if π " j P U , k " 1 and t P r0, τ s, and ' δ π pk´1qτ `t1 ,δ0 " δ j t 1 ,δ 1 with δ 1 " δ π 1 pk´1qτ,δ0 , if k ě 2, t 1 P r0, τ s, π " π 1 ˚j for some j P U and π 1 P U k´1 .

Controlled pR, Sq-stability

Given a rectangle R Ă R d and a rectangle S Ă R d such that R Ď S, we now extend the problem of "controlled pR, Sq-stability", as defined in Le [START_REF] Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF] for deterministic switched systems, to SDEs, as follows:

For any starting point x 0 P R, find a pattern π of length k such that ' EX π t,x0 P R for t " kτ ' EX π t,x0 P S for all t " τ, 2τ, 3τ, . . . . It is easy to see that, in order to solve this problem, it suffices to exhibit a finite set of points z 1 , . . . , z p of S, and a positive real δ 0 ą 0 such that:

(1) all the balls Bpz i , δ 0 q, i " 1, . . . , p, cover R, and are included into S (i.e. R Ď Ť p i"1 Bpz i , δ 0 q Ď S); (2) for each i " 1, . . . , p, there is a pattern π of length k such that: ' B i,π,t Ď S for t " τ, 2τ, . . . , pk ´1qτ , and ' B i,π,t Ď R for t " kτ . where B i,π,t :" BpE Xπ t,zi , δ π t,δ0 q. The program mentioned in Section 2.4, has been extended in order to find, by exhaustive enumeration3 , patterns that make the balls covering R return to R, and such that the intermediate balls (at t " τ, 2τ, . . . ) belong to S. Please refer to Appendix B for illustrations and to Le [START_REF] Coënt | An improved algorithm for the control synthesis of nonlinear sampled switched systems[END_REF] for more information on the algorithm.

We now give an application of this program. Example 2. Consider the system (see [START_REF] Zamani | Symbolic control of stochastic systems via approximately bisimilar finite abstractions[END_REF][START_REF] Zamani | Symbolic models for stochastic switched systems: A discretization and a discretization-free approach[END_REF]): dx 1 " p´0.25x 1 `ux 2 `p´1q u 0.25qdt `0.01x 1 dW 1 t dx 2 " ppu ´3qx 1 ´0.25x 2 `p´1q u p3 ´uqqdt `0.01x 2 dW 2 t where u " 1, 2. For τ " 0.5, ∆ t " 10 ´4, one finds (for all modes u " 1, 2): q " 0, D " 1.36, L g " 0.01, λ " 0.25; for z " p´4, ´3.8q: E 2,z " 893.31, E 4,z " 2.14 ¨10 5 , F 2,z " 0.002, F 4,z " 4.9 10 ´6; and for z " p0, 3q: E 2,z " 543.22, E 4,z " 7.94 ¨10 4 , F 2,z " 0.0442, F 4,z " 0.00178. Our program shows pR, Sq-stability of the system for R " r´5, 5s ˆr´4.4s and S " r´8, 8s ˆr´7, 7s: given a covering of R with balls of radius δ 0 " 0.1, the program finds, by exhaustive search, patterns of length ď 5 that make the balls return to R. It takes 6 hours of CPU time. Figures 3,4, 5 and 6 depict in black the initial balls (at t " 0) centered at the corners of R; and for each initial ball, the pattern that sends the ball back to R (at time t " kτ ); the intermediate balls (at t " τ, 2τ, . . . , pk ´1qτ ) are depicted in red, and 200 sampling trajectories drawn in blue. Fig. 3. initial ball Bpz, δ 0 q with z " p´5, 4q and δ 0 " 0.1;

pattern " p1 ¨1 ¨1q; τ " 0.5 Fig. 4. initial ball Bpz, δ 0 q with z " p5, 4q and δ 0 " 0.1; pattern " p2 ¨2 ¨2q; τ " 0.5

Other applications

Our Euler-based method can also be used to control systems in order to achieve reachability properties. We sketched out this point in the following example. Example 3. (the slit problem) Fig. 5. initial ball Bpz, δ 0 q with z " p5, ´4q and δ 0 " 0.1; pattern " p2 ¨2q; τ " 0.5 Fig. 6. initial ball Bpz, δ 0 q with z " p´5, ´4q and δ 0 " 0.1; pattern " p1 ¨1 ¨1 ¨1 ¨1q; τ " 0.5

The problem is adapted from [START_REF] Morzfeld | Implicit sampling for path integral control, Monte Carlo localization, and SLAM[END_REF]. The controlled dynamics is:

dX " udt `dW, X 0 " 1 with mode u P t´6, ´5, ´4, ´3, ´2, 1, 0, 1, 2, 3, 4, 5, 6u. We have (at t " 0.5) a slit at x P r´1, ´4s. The objective is thus to control the system so that xptq P S " r´1, ´4s at t " 0.5.

One has, for all modes: q " 0, D " 0, L g " 0, λ " 0. For δ 0 " 0.5, an initial point z " 1 and a sampling time τ " 0.5 with subsampling ∆ t " 10 ´3, one has for mode u " ´6: E 2,z " 144, E 4,z " 20736, F 2,z " 4, F 4,z " 16; and for mode u " 0: E 2,z " 0, E 4,z " 0, F 2,z " 4, F 4,z " 16.

Suppose that all the trajectories start at x 0 with x 0 P Bpz, δ 0 q (i.e., |x 0 ´z| ď 0.5), with z " 1 and δ 0 " 0.5. When there is no control (u " 0), at time t " 0.5, the expected value of X t,x0 is in Bpz 1 , δ t,δ0 q with z 1 " 1 and δ t,δ0 " 2. From Markov's inequality, it follows that the trajectories pass by S " r´1, ´4s at t " 0.5 with low probability: see Figure 7. On the other hand, with control u " ´6, at time t " τ " 0.5, the expected value of X t,x0 is now in Bpz 1 , δ t,δ0 q with z 1 " ´2 and δ τ,δ0 " 2. This . Example 3 without control (u " 0) for t P r0, τ s; initial ball Bpz, δ 0 q with z " 1 and δ 0 " 0.5; final ball Bpz 1 , δ 1 q (at t " τ " 0.5) with z 1 " 1, δ 1 " 2 0 0.2 0.4 0.6 0.8 1 -6 -4 -2 0 2 Time X Fig. 8. Example 3 with control pattern p´6 ¨0q; initial ball Bpz, δ 0 q with z " 1 and δ 0 " 0.5; intermediate ball Bpz 1 , δ 1 q (at t " τ " 0.5 ) with z 1 " ´2, δ 1 " 2; final ball Bpz 2 , δ 2 q (at t " 2τ ) with z 2 " ´2, δ 2 " 3.6 explains why the trajectories now pass by S " r´1, ´4s at t " 0.5 with high probability: see Figure 8.

FINAL REMARKS AND FUTURE WORK

We have explained how to use an Euler-based method in order to control stochastic switched systems. We have focused our work on the property of pR, Sq-stability, but it can also be used for achieving reachability properties.

In the future, we plan to experiment the method with examples where the drift functions behave polynomially. We would like also to find bounds not only for the expected values of the solutions, but for their variance.

  Fig.7. Example 3 without control (u " 0) for t P r0, τ s; initial ball Bpz, δ 0 q with z " 1 and δ 0 " 0.5; final ball Bpz 1 , δ 1 q (at t " τ " 0.5) with z 1 " 1, δ 1 " 2

Note that, in the figures, all the end points (at t " τ ) of the sampling trajectories lie in the final ball, but this is not true in general; we only know by Proposition 1 that, for any starting point x 0 in the initial ball, the expected value of the end point lies in the final ball.

For the sake of notation conciseness, we suppose that the number of steps of subsampling N is the same whatever the mode j of the pattern π is, hence the stepsize of the subsampling is always equal to ∆t " τ {N ; in full generality, we should write N j instead of N to express the dependence.

The enumeration can be accelerated using different branch-andbound heuristics (see[START_REF] Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF]).

Proof. Let t P rk∆ t , pk `1q∆ t q. Then (using the inequality pa `bq r ď 2 r pa r `br q): }X t ´X t } r " }pt ´tk qf pX k q `gpX k qpW t ´Wt k q} r ď 2 r pp∆ t q r }f pX k q} r `}gpX k q} r }W t ´Wt k } r q }X t ´X t } r ď 2 r pp∆ t q r p}f pX k q ´f p0q} r `}f p0q} r q `p}gpX k q ´gp0q} r `}gp0q} r q}W t ´Wt k } r q ď 2 r pp∆ t q r pDpp1 `


N Perrissin

S Lequeux

N Strelkov

L Vila

L Buda-Prejbeanu

S Auffret

R C Sousa

I L Prejbeanu

B Dieny

Spin transfer torque magnetic random-access memory: Towards sub-10 nm devices

Magnetic Random-Access Memory (MRAM) is a non-volatile class of solid-state storage device where the information is stored in the magnetic state of a ferromagnetic layer. Microelectronic industry has recently shown a strong interest for MRAM as they are very promising for embedded RAM applications and particularly embedded FLASH replacement 1- 6 . The main building bloc of an MRAM is a trilayer FM1/I/FM2 (Fig. 1.a) named magnetic tunnel junction (MTJ) 7 , where FM1(2) refer to ferromagnetic layers and I to a thin insulating tunnel barrier (~ 1.2 nm). Nowadays researches are mainly focused on perpendicularly (p) magnetized tunnel junctions [8][9][10][11] written by spin transfer torque (STT) 12,13 , where the perpendicular axis is defined by the growth direction. FM1 is called the reference layer, its magnetization is pinned in one specific direction (for example, the up direction). FM2 is called the storage layer or free layer, its magnetization is free to be moved between its two stable states (up and down). The resistance of the stack depends on the relative orientation of the two magnetizations. The parallel (P) state (up/up) is a low resistance state while the antiparallel (AP) state (up/down) is a high resistance state, thus respectively coding the "0" and "1" of the binary logic. The tunneling magnetoresistance (TMR) is defined by the ratio = .

The P and AP states are separated by an energy barrier Eb, used to define the thermal stability of the memory ∆= E /k T with kBT the thermal energy (Fig. 1.b). For practical applications, one aims for ∆ values from 60 to 100 [START_REF] Dieny | Introduction to magnetic random-access memory[END_REF][START_REF] Apalkov | [END_REF] .

In STT-MRAM, the storage layer has commonly an almost cylindrical shape of diameter D and thickness L. When the memory cell is sufficiently small to behave in a macrospin manner (typically below 20nm in diameter), its stability is given by Eq. 1, where µ0 is the vacuum permeability, A the storage layer area and MS the saturation magnetization. The first term inside the bracket is the shape anisotropy, where Nxx and Nzz are the demagnetizing factors. The second and third terms respectively represent the interfacial anisotropy arising at the FeCoB/MgO interface and a possible uniaxial bulk anisotropy of magnetocrystalline or magneto-elastic origin. The critical voltage Vc is defined as the minimum voltage required to write the memory state at 0 K (Eq. 2), where RA is simply the resistance times area product of the memory point. The write voltage at non-zero temperature can be lower than Vc and depends on the pulse duration 14 .

(Eq. 1)

∆ = E k T = A k T 1 2 μ M L N -N ! + K $ + K % L& V ( = 4e ℏ α RA η k T ∆ A
In conventional p-MTJ, as the thickness of the free layer (≃ 1.4 nm) is much lower than its diameter (> 30 nm), so that 1 22 ~0, 1 44 ~1. The shape anisotropy (first term in Eq.1) thus favors in-plane anisotropy. The perpendicular magnetization is then obtained thanks to a strong interfacial anisotropy at the FM/I interface which overcomes the in-plane shape anisotropy. The MTJ stack comprises the central block: FeCoB(~1.3nm)/MgO(~1.2nm)/FeCoB(~1.4nm) which combines several advantages: After annealing, it provides a good crystallographic structure match between the FeCo(001) magnetic electrode and the MgO(001) tunnel barrier 16 , a quite low Gilbert damping constant 17,18 , a strong interfacial anisotropy 19,20 and a large spin polarization leading to a high TMR ratio 11,21 . However, at sub-20nm dimensions (macrospin behavior), ∆ scales with D², inevitably leading to excessively low values of ∆ at small diameters meaning too high bit error rate in standby or during read. In conventional STT-MRAM, one can insure ∆ > 60 down to diameters of around 20 nm 22 . It can be improved by capping the free layer with another MgO layer to increase the total interfacial anisotropy. This enables to achieve ∆ > 60 down to diameters of ~15 nm [START_REF] Hu | IEEE International Electron Devices Meeting[END_REF][START_REF] Cuchet | [END_REF] .

To go to even smaller diameters while maintaining sufficiently high values of ∆, we propose to take advantage of a perpendicular shape anisotropy (PSA) by patterning the storage layer with an aspect ratio L/D of the order or greater than one. By doing this, the shape anisotropy becomes positive and no longer fights against the interfacial anisotropy. To keep good interfacial properties (large TMR and large interfacial anisotropy), one can keep a thin FeCoB layer at the interface with the barrier, then deposit a thick FM layer with low Gilbert damping on top. This approach allows to have both good interfacial properties and to give to the storage layer the desired properties (low damping, weak thermal variation of anisotropy and TMR, etc). Let us for example consider the case of MgO/FeCoB(1.4nm)/Co(L-1.4nm) with a total thickness L (Fig. 2). The smallest thicknesses correspond to the conventional p-STT-MRAM regime (L < 2 nm) while the larger ones correspond to the novel PSA-STT-MRAM regime (t > 3 nm). The diagram clearly illustrates the inability of p-STT-MRAM to maintain ∆ > 60 at sub-20 nm nodes. On the other hand, once a thick Co layer is added on top of the storage layer, a strong perpendicular anisotropy can be recovered originating from the shape itself. By changing the thickness and/or the diameter of the storage layer, one can tune the stability over a very wide range, even for sub-10 nm nodes. For this example, a stability of ∆ = 60 can be achieved down to a diameter of 4 nm with a total storage layer thickness of 32 nm. We superimpose to the stability diagram a boundary evaluated by micromagnetic simulations, indicating the validity range of Eq. 1. Finally, we also superimpose some iso-V ( /α RA lines to emphasize the need of low RA tunnel junctions. Indeed, by considering a write voltage of V ( = 0.5 V and a damping constant α = 0.01, one needs a RA of 1 Ω. μm to switch a Co-based PSA-STT-MRAM with D = 7 nm, L = 16 nm (⇒ ∆ = 60). The smaller the diameter, the smaller α and RA have to be.

Co-based and FeNi-based PSA-STT-MRAM with thickness of 60 nm have been patterned by IBE in a two-step process. The storage layer was first etched vertically thanks to a hard mask of Ta with D = 30 nm. The pillar was then trimmed to sub-20 nm diameters by etching at a (Eq. 2) grazing angle. The lower diameter observed after trimming is 8.1 nm. Hysteresis resistance versus magnetic field loops (Fig. 3) show two well defined resistance states with TMR values around 50 % and ∆ values coherent with micromagnetic simulations.

Fig. 1 .Fig 2 .

 12 Fig. 1. (a) Sketch of a conventional p-MTJ. In red the storage layer FeCoB(1.4 nm) with two stable states highlighted by the white double-arrow, in yellow the MgO(1.2 nm) tunnel barrier and in green the reference layer FeCoB(1.3 nm) with the magnetization pinned in one direction highlighted by the white single-arrow. (b) Energy profile of the p-MTJ versus the orientation of the magnetization of the storage layer. The stability factor is defined by the ratio ∆= : ; < = > .

Fig. 3 .

 3 Fig. 3. (a), (b)Resistance R versus perpendicular field for MTJ of various diameters D in the cases of a 60 nm thick storage layer of NiFe and Co respectively. (c) Distribution of coercive fields Hc as a function of D for NiFe (black circles) and Co (red circles) storage layer. (d) Distribution of the TMR ratio as a function of D for NiFe (black circles) and Co (red circles) storage layer.

This work was funded by ERC Adv grant MAGICAL 669204. The authors would like to thank Jannier Roiz for fruitful discussions and valuable suggestions.


Gautier Berthou

Kevin Marquet
email: kevin.marquet@insa-lyon.fr

Tanguy Risset

Arnaud Carer

Henri-Pierre Charles

Steven Derrien

Ivan Miro-Panades

Davide Pale

Isabelle Puaut

Fabrice Rastello

Erven Rohou

Guillaume Salagnac

Olivier Sentieys

Bahram Yarahmadi

Inria Socrate (Leader, https://team.inria.fr/socrate/): embedded system design, communication systems, operating systems and memory management [START_REF] Berthou | Peripheral state persistence for transiently-powered systems[END_REF]. Inria Cairn (https://team.inria.fr/cairn/): energy-efficient system-on-chip and reconfigurable architecture, energy harvesting [START_REF] Berder | PowWow: Power Optimized Hardware/Software Framework for Wireless Motes[END_REF], compilers and design tools. Inria Pacap (https://team.inria.fr/pacap/): compiler optimization, compilation-time energy evaluation.

Inria Corse (https://team.inria.fr/corse/): link between compilation and runtime. CEA Leti (Lisan, http://www.leti-cea.fr/): ultra-low power IOT platform [START_REF] Beign | L-IOT: a Flexible Energy Efficient Platform Targeting Wide Range IoT Applications[END_REF].

CEA List (Lialp, http://www-list.cea.fr/): just-in-time compilation.

Expected Results

After 4 years, the project should have achieved the following: A NVRAM based platform demonstrator powered by harvesting. An innovating NVRAM based processor architecture designed for VLSI integration.

A solid system-level simulator for a generic NVRAM platform. ) by an FRAM-based MSP430FR5739 will provide a first and quickly available prototype. NVRAM-based processor architecture: we propose a collaboration between Cairn and CEA Leti to explore new techniques at architectural and design level to store the intermediate computation states before power failures occur, through architecture improvements. The core of the micro-architecture will be based on a RISC-V processor (https://riscv.org/), and we will use CEA Leti experience on the L-IOT architecture [START_REF] Beign | L-IOT: a Flexible Energy Efficient Platform Targeting Wide Range IoT Applications[END_REF].

Task 2: Compiler Optimizations and Worst-Case Energy Consumption

Energy Model: we will derive a precise energy model of a NVRAM-based processor. Code Generation: we will provide a Worst Case Energy consumption (WCEC) driven code generation . Checkpoints: The code generator could influence the checkpointing strategy. Quality of Service and Approximate Computing: We will also explore many ways to reduce energy consumption when the battery runs low.

Task 3: Runtime memory management

Provide support to the NVRAM architecture: the runtime system will have to share the non-volatile memory across the various applications, according to their needs and priorities. Choose how energy is used energy policy must define the way the system should spend the available energy depending on the length of each task to execute for instance. Ensure memory consistency: problems to be solved: processor state is volatile, devices' states are volatile and a priority between applications when saving data of must be established [START_REF] Berthou | Peripheral state persistence for transiently-powered systems[END_REF].

Task 4: Experimental platforms and use cases Use cases: Use cases will include a sensing application, a cryptography application, and a video decoding application. Experimental platform: the ZEP prototype will consist in the integration of recent FRAM MSP430 from Texas Instrument (MSP430FR5739 for instance) integrated in a new version of the PowWow sensor [START_REF] Berder | PowWow: Power Optimized Hardware/Software Framework for Wireless Motes[END_REF] node in which an the Actel Igloo FPGA also includes NVRAM. Simulation platform: for more advanced design issues, we will use simulation tools provided by CEA (for L-IOT) and other vendors for NVRAM to provide a precise system-level simulator

Figure :

 : Figure: ZEP Teams involved with their correspondants scientists.

Figure :

 : Figure: Overview of the ZEP resulting prototype

  Examples of batteryless RFID tag and typical On/Off cycles of platforms supplied by energy harvesting, with checkpointing and restoration Insa-Lyon or University of Rennes are involved in Inria project-teams.The ZEP project is an Inria Project Lab (IPL), funded by Inria with more than 15 person-years hired over 4 years, including 3 PhDs. The Insa-Spie IoT Chair is also participating to the research effort of Citi Lab on IoT.

				The INRIA ZEP project:
			NVRAM and Harvesting for Zero Power
					Computations
						Partners and Methodology
		Supply voltage			Inria and CEA
		Vboot				Inria, the French National Institute for computer science and applied mathematics, promotes scientific
						excellence for technology transfer and society. Inria is structured in research teams called project-teams.
						CEA, the French Alternative Energies and Atomic Energy Commission is a key player in technological
		Vsave				research for industry. CEA research institutes (List and Leti) are structured in research teams called labs.
	Figure: Context	Vdeath	Off time	Lifecycles Challenges	Time	Other institutes such as Planning
	The Internet of Things scenario predicts billions of communicating objects in a near future.	Provide a new computation model adapted to NVRAM based micro-systems.	T0.2: workshop Task O: IPL management T0.1: meeting	M0.1 W0.1	2017 M0.2	M0.3	2018 M0.4	M0.6 W0.2	2019 M0.7	M0.8	2020 M0.9 W0.2	M0.10
	Maintenance cost and pollution issues are	Master harvesting techniques and energy	Task 1: NVRAM architecture T1.1: ZEP prototype			D1.1		D1.2		PhD 1
	major bottleneck for IoT. Emerging non-volatile memory (NVRAM)	control in ultra-low power domain. Provide real computer science contributions	T1.2 PhD 1 T2.1: PhD 2 Task 2: Worst Case Energy Consumption							PhD 2
	technologies and harvesting techniques predict new kinds of batteryless embedded sensors. Checkpointing [1] of platform state will be a	(langage, compiler, operating systems, etc...) Deliver innovative prototypes and a useful system-level simulators.	Task 3: NVRAM OS T4.3: demonstrators T3.1: Phd 3 Task 4: Demonstrators T4.2: platform integration					D4.3	M4.1	D4.5	PhD 3 D4.4 D4.4 D4.6
	central technique is these systems.
		Teams and Skills
		Inria Pacap	Inria Socrate	Inria Corse
				CEA Leti
		Rennes
	Inria Cairn
			Lyon
			Grenoble
				CEA List

ZEP project: Context and Objectives


Sandeep Manandhar

Patrick Bouthemy

Erik Welf

Philippe Roudot

Charles Kervrann

A SPARSE-TO-DENSE METHOD FOR 3D OPTICAL FLOW ESTIMATION IN 3D LIGHT-MICROSCOPY IMAGE SEQUENCES

Keywords: 3D PatchMatch, Interpolation, Fluorescence Microscopy, 3D Optical Flow

We present a two-stage 3D optical flow estimation method for light microscopy image volumes. The method takes a pair of light microscopy image volumes as input, segments the 2D slices of the source volume in superpixels and sparsely estimates the 3D displacement vectors in the volume pair. A weighted interpolation is then introduced to get a dense 3D flow field. Edges and motion boundaries are considered during the interpolation. Our experimental results show good gain in execution speed, and accuracy evaluated in computer generated 3D data. Promising results on real 3D image sequences are reported.

INTRODUCTION AND RELATED WORK

Estimating motion of cells and of subcellular particles is crucial for deciphering cell mechanisms and understanding cell behaviours. Modern 3D light microscopy (LM) for cell biology enables to observe cell dynamics at a good resolution both in space and time, motivating the development of 3D optical flow methods. However, the acquired 3D image sequences exhibit several specificities making 3D motion computation a difficult problem. The 3D stacks are not isotropic, the axial resolution is coarser than the lateral ones, and the stack usually comprises a small number of slices compared to the 2D image size. As a consequence, particles undergoing 3D movement might disappear in between two slices. Images are noisy; in most cases, due to the phototoxicity issue, the photon budget per image is reduced to favour an as long as possible acquisition over time. Finally, images are poorly textured, especially fluorescence image sequences which usually contain bright spots over cluttered background.

2D optical flow has been an active research field in computer vision for several decades. Optical flow methods estimate velocity components by relying on the brightness constancy constraint and imposing spatial smoothness of the velocity field, either locally or globally, while attempting to preserve motion discontinuites and to detect occlusions [START_REF] Fortun | Optic flow modeling and computation: a survey[END_REF]. Recent performing methods proceed in two steps, computing a sparse field with an efficient matching method, especially to cope with large displacements and small moving objects, then interpolating the latter to recover a dense flow field as in [START_REF] Revaud | EpicFlow: Edge-preserving interpolation of correspondences for optical flow[END_REF]. Application of 2D optical flow methods to 2D fluorescence microscopy image sequences were reported in [START_REF] Delpiano | Performance of optical flow techniques for motion analysis of fluorescent point signals in confocal microscopy[END_REF][START_REF] Fortun | Aggregation of patch-based estimations for illumination-invariant optical flow in live cell imaging[END_REF].

3D optic flow computation has been far less investigated. A related problem is 3D image registration in medical imaging as surveyed in [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]. Let us also quote the local approach for 3D optical flow based on the all-pass filtering method defined in [START_REF] Gilliam | 3D motion flow estimation using local all-pass filters[END_REF] and applied to 3D MRI data to handle respiratory motion. However, 3D registration in medical imaging mainly involves smooth deformation estimation over a dense isotropic volume, which is not the case for biological microscopy image stacks. Only a few investigations have been undertaken in living cell microscopy. Variational optic flow methods were tested for the first time in [START_REF] Hubeny | Estimating large local motion in live-cell imaging using variational optical flow[END_REF]. In [START_REF] Amat | Fast and robust optical flow for time-lapse microscopy using super-voxels[END_REF], an algorithm is proposed where smoothness constraint is imposed by Markov Random Fields over the graph of presegmented super-voxels in the 3D volume. The method assumes that objects do not undergo complex motion, and computes a velocity vector per super-voxel. The method is used for the purpose of visualization by [START_REF] Kappe | Reconstruction and visualization of coordinated 3D cell migration based on optical flow[END_REF]. The authors of [START_REF] Boquet-Pujadas | BioFlow: a non-invasive, image-based method to measure speed, pressure and forces inside living cells[END_REF] integrate biophysical models in the estimation framework, and simultaneously compute force, pressure and velocity fields. The method, called BioFlow, originally applied to 2D microscopy image sequences, is extended to 3D data, and supplies convincing results on cell migration image sequences. Nevertheless, the method is effective if biophysical models are available, and takes about an hour for execution. In [START_REF] Tektonidis | Diffeomorphic multi-frame non-rigid registration of cell nuclei in 2D and 3D live cell images[END_REF], local optical flow computed in the log-domain is used within a non-rigid registration framework based on diffeomorphic transformations and applied to 2D and 3D live cell image sequences. Computation time is reported to be around 33 minutes per frame.

In the 2D optical-flow paradigm, patch-matching based methods have been gaining state-of-the-art performance [START_REF] Revaud | EpicFlow: Edge-preserving interpolation of correspondences for optical flow[END_REF][START_REF] Hu | Efficient coarse-to-fine PatchMatch for large displacement optical flow[END_REF][START_REF] Bailer | CNN-based patch matching for optical flow with thresholded hinge loss[END_REF]. Following this idea and motivated by the design of a general-purpose fast 3D optical flow estimation method able to deal with diverse 3D LM modalities and the variety of cell image contents, we propose a new 3D optical flow method for 3D LM image sequences. It is structured in two stages: sparse estimation and dense interpolation. For the first step, we develop a 3D version of PatchMatch [START_REF] Barnes | PatchMatch: A randomized correspondence algorithm for structural image editing[END_REF] on supervoxels. The second step is achieved with a weighted interpolation scheme.

The remainder of the paper is organized as follows. In Section 2, we explain how we handle the 3D stack for matching. Section 3 describes our 3D PatchMatch method to compute the sparse velocity field. Section 4 presents the weighted interpolation to get the dense flow. In Section 5, we report and comment experimental results on 3D synthetic and real images. Concluding remarks are given in Section 6.

HANDLING OF THE 3D STACK

A LM image volume consists of a stack of 2D slices where each slice depicts information corresponding to a given depth layer of the cell. Considering this volume structure, we compute the 3D displacement vectors only at the points of the stack grid, that is, at the pixels of every slice grid. In what follows, the term displacement vector is interchangeably used as flow vector.

Our framework starts with the computation of a set of correspondences in a pair of LM image volumes I 1 and I 2 . To be efficient, matching must be performed on a sparse set of points. Superpixels are good candidates for that purpose, since we can easily control their density, they are almost regularily distributed over the image while taking into account the intensity patterns. We adopt the SLIC method [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]. We obtain K non overlapping superpixels in the slices of the volume I 1 . Here, supervoxels restrict to superpixels taken as their intersection with slice of the stack. Thus, we will still use the term superpixel in the sequel. Let S be the set of superpixels and s k be the k th superpixel.

For each s k , we define a set of 3D neighbouring superpixels N k = {s ki , i = 1..n} in I 1 which can lie in different slices. This neighborhood will be used in the propagation step of the 3D PatchMatch method as explained in the next section. N k comprises the n superpixels with minimum distance to s k according to the similarity measure defined by:

d k,j = exp(- (c j -c k ) 2 ν d - (I 1 (c j ) -I 1 (c k )) 2 ν I ) (1)
where c k and c i are the centres of superpixels s k and s j , and ν d and ν I are respectively the spatial distance and intensity value normalizers.

SPARSE MATCHING

Correspondences between points of I 1 and I 2 are defined as labels {l} such that I 1 (x) and I 2 (x + l(x)) form a matched pair, where label l(x) is the 3D displacement vector of a voxel located at x in I 1 . To establish the correspondences, we opt for PatchMatch [START_REF] Barnes | PatchMatch: A randomized correspondence algorithm for structural image editing[END_REF] which is a fast method to compute correspondence between two images. We extend Patchmatch to 3D volumes, where correspondences act between voxels. The complexity of PatchMatch largely depends on the image grid size which is even bigger in case of volumetric data. We reduce the grid size by sparsely sampling the volumetric data into superpixels as explained in the preceding section. A similar strategy was used by [START_REF] Lu | PatchMatch filter: Edge-aware filtering meets randomized search for visual correspondence[END_REF] to compute approximate nearest neighbouring field (ANNF) between two 2D images. A smooth flow field suggests that similar voxels in a close vicinity move similarly. Then, we estimate a label for only the center c k of superpixel s k , as representative of its vicinity. A more sophisticated method is proposed in [START_REF] Giraud | SuperPatchMatch: An algorithm for robust correspondences using superpixel patches[END_REF], introducing a specific patch neighborhood.

As with the original PatchMatch [START_REF] Barnes | PatchMatch: A randomized correspondence algorithm for structural image editing[END_REF], the framework requires iterations of propagation and random search. For initialization, a label for c k is chosen as the 3D vector to the center c k of superpixel s k in I 2 selected as the most similar one according to the distance [START_REF] Fortun | Optic flow modeling and computation: a survey[END_REF], where I 2 (c k ) is substituted for I 1 (c j ). It proved to be better than the random initialization done in [START_REF] Barnes | PatchMatch: A randomized correspondence algorithm for structural image editing[END_REF]. Then, we iterate propagation and random search. Scanning of the volume is performed in scan-line order in each slice, visiting each c k , and slices are visited in a top-to-bottom way through the stack, as shown in Fig. 1.

We now briefly explain propagation. To update the label at c k , we take the current labels of the centers of the superpixels in N k . This provides n possible correspondences for c k . We select the label, among the n labels plus the current label of c k , with a patch-based similarity measure. The patch is a local 3D patch of size 9 × 9 × 3 around the point. For patch similarity, we adopt the L 2 norm of intensities between two patches. Other norms (L 1 , Lorentzian) did not improve results. After the best label lk has been inferred for c k in the propagation stage, the random search proceeds by first creating a search region of radius R centered at c k + lk in I 2 . Then, a set of additional label candidates {l r } are sampled with an exponentially decreasing R through iterations. The superpixel-based sparse sampling allows us to take floating point label during the random search, which would have been computationally prohibitive otherwise. A new label among the l r 's is assigned to c k if the similarity is improved, compared to label lk . At the very end of the alternate iterations, a sparse 3D flow field is obtained. A median filter is applied to the 3D vectors of the field to remove outliers that significantly disobey smoothness in the flow field. The filtered vectors are then used to interpolate displacements for the remaining voxels in I 1 as described in the next section.

SPARSE-TO-DENSE INTERPOLATION

To get a dense 3D flow field v between I 1 and I 2 , the 3D displacement vectors for all the points in each superpixel segment need to be computed. To do so, we locally interpolate the sparse set of displacement vectors supplied by the first stage. Of course, we do not want to just duplicate the flow vector of c k to all points of superpixel segment s k , the undergone motion being not necessarily locally translational.

To this end, we design an interpolation scheme inspired by EpicFlow interpolation introduced in [START_REF] Revaud | EpicFlow: Edge-preserving interpolation of correspondences for optical flow[END_REF]. EpicFlow interpolation is an edge-aware interpolation technique that has been known to perform well in 2D optical flow benchmarks. It assumes that motion discontinuities are subset of intensity edges, and derives weights using an approximation of geodesic distance between points in Voronoi cells with intensity edges as cost map. In our case, we exploit the superpixel segmentation instead of Voronoi tessellation. In addition, motion boundaries will be taken into account in the interpolation stage in a different manner.

Let us consider a point x in superpixel s k . The question is how to involve known vectors in the neighbourhood to interpolate the displacement or flow vector at x. For this, we introduce a flow-aware weight that encodes the respective orientation of known vectors. We allow interpolation with vectors of neighbours only when the flow vector at c k and those of the superpixel centers of its neighbours in N k are converging. Here, N k is the subset of neighbouring superpixels lying in the same slice as s k . Given the anisotropy of the volume, we do not involve off-the-plane superpixels.

We build an undirected graph G = (V, E) where V is the set of superpixel centers c k 's , and E is the set of edges between adjacent superpixels in S. Let a and b be any two vertices (centers of superpixel s a and s b respectively). We define edge weights w e using an approximation of geodesic distance between two vertices in V defined by:

D g (a, b) = pi∈L a,b q pi , (2)
where p i is a point in the line segment L a,b between a and b. q pi is the intensity edge response at p i which can be derived from an intensity edge detector or simply from local intensity contrast. Noisy responses are filtered out. We aim to nullify the interaction between two flow vectors if they are moving away from each other. Assuming no strong local curl in v, we encourage interpolation from v(a) and v(b) to get flow vector at x in superpixel s a , if the flow vectors v(a) and v(b) are converging as illustrated in Fig. 2. Now, we define edge weights on E as:

w e (a, b) = 1 Dg(a,b), if d f > 1 0, otherwise (3)
where d is the Euclidean distance between a and b, and f the one between points a + v(a) and b + v(b) as shown in Fig. 2b. Given the weights derived from the respective orientation of the flow vectors and the spatial location of the superpixels, we interpolate for the flow vector at x as follows:

v(x) = w a,x v(a) + i∈N a w e (a, i)v(i) w a,x + i∈N a w e (a, i) , ∀x ∈ s a , (4)
where w a,x = 1/D g (a, x). The interpolation of the 3D flow is performed slice by slice.

EXPERIMENTAL RESULTS

First, we objectively evaluated the method performance on computer-generated sequences with ground-truth, obtained by applying known transformations on real images issued from cell tracking challenge [START_REF] Masaka | A benchmark for comparison of cell tracking algorithms[END_REF]. The image volume is of size 275×231×35. 100 uniformly random transformations were applied for each class of transformation, from 3 to 10 units for xand y-translation, from 2 to 3 units for z-translation, rotation range around the z-axis was [-10 • , 10 • ], and scaling range [0.5, 3] in x and y, and [0.5, 1.5] in z. We compared the 3D flow field obtained by our method against [START_REF] Amat | Fast and robust optical flow for time-lapse microscopy using super-voxels[END_REF]. Table 1 contains average end-point error (AEE) and average angular error (AAE). Our method handles motion including rotation or scaling better than [START_REF] Amat | Fast and robust optical flow for time-lapse microscopy using super-voxels[END_REF]. The latter method, unlike ours, computes a vector for a supervoxel, and assigns the same vector for the involved voxels which is not appropriate to handle complex motion. Figure 3b contains a portion of flow fields in a slice in the case of rotation and translation. The flow field computed by our method closely follows the ground-truth flow field.

We have studied the tradeoff between runtime and accuracy with a varying number of superpixels. The method is implemented in C++ (no optimized) and run on a computer with 2.8 GHz Intel i7 processor and 16 GB of RAM. The parameters are set as follows: the size of initial random search region is [40 × 40 × 7], ν I = 0.3 and ν d = 100. Plots in Fig. 3a show a gain in accuracy with smaller size of superpixels (larger the K, smaller the superpixel size). The convergence is reached after 5-10 iterations (best viewed in electronic form for zooming or in color printed form).

We also applied our method in two real sequences [START_REF] Dean | Deconvolution-free subcellular imaging with axially swept light sheet microscopy[END_REF][START_REF] Dean | Diagonally scanned light-sheet microscopy for fast volumetric imaging of adherent cells[END_REF]. The first sequence depicts blebbing in a MV3 cell (Fig. 4a-b). Blebbing is a phenomenon where cell membrane protrudes increasing the surface area of the cell. Protrusions, referred to as blebs, appear and disappear in interval of minutes. The appearance of a bleb is accompanied by stretching of a local region of the cell membrane. The region gradually retracts to adhere to conformity with the cell surface. Figure 5 plots the displacement vectors in one of the protrusions in the cell. The vectors tend to stretch the cell membrane. The computation time to estimate the 3D flow field between two volumes was 163s, with resp. 19s, 120s and 24s for superpixel generation, 3D patch matching, and interpolation.

The second real sequence is a MV3 cell undergoing large scale protrusion depicted in Fig. 4c-d. The main body of the cell in middle section squeezes as protrusions abruptly appear. A cross-section view of its 3D flow field is plotted in Figure 6 which suggests shrinking in the middle region and stretching in protrusions. The total computation time for one 3D flow field was 101 seconds, with 49, 44 and 8 seconds for resp. superpixel generation, 3D patch matching, and interpolation.

CONCLUSION

We have presented a method for fast 3D optical flow estimation in 3D LM image volumes. It involves the interpolation of a sparse 3D flow field obtained by matching superpixels per slice. We have developed a 3D PatchMatch algorithm to find correspondences. For densification, we have designed a weighted interpolation taking intensity edges and local flow variations into consideration. Tests performed on synthetic cases show its effectiveness, and experiments with real cases look promising. The 3D computed flow field could be further refined by post-processing it with a variational algorithm.

Fig. 1 .

 1 Fig. 1. Scanning order in the stack for 3D PatchMatch. The slices are regularly tessellated with cell size adapted to superpixel mean size. Each cell encapsulates a superpixel center.At the very end of the alternate iterations, a sparse 3D flow field is obtained. A median filter is applied to the 3D vectors of the field to remove outliers that significantly disobey smoothness in the flow field. The filtered vectors are then used to interpolate displacements for the remaining voxels in I 1 as described in the next section.

Fig. 2 .

 2 Fig.2. Interpolation of flow vector at x (in red) according to the central flow vector (in purple) and the neighbouring flow vectors (in green). The purple and green points are the centers of their respective superpixels delineated by black lines.

  3D field in a slice.

Fig. 3 .Fig. 4 .

 34 Fig.3. Left: Accuracy converges to a value after few iterations (number of iterations is mentioned above each plotted points). Right: 3D flow vectors estimated in a slice with ground-truth in green, ours in blue,[START_REF] Amat | Fast and robust optical flow for time-lapse microscopy using super-voxels[END_REF] in red (ground-truth: rotation of 10 • around z-axis, 3D translation of[10,-20,3]).

Fig. 5 .

 5 Fig. 5. Cross section view of flow field for the 1st sequence.

Fig. 6 .

 6 Fig. 6. Cross section view of flow field for the 2nd sequence.

Table 1 .

 1 Error evaluation on 3D synthetic tests.

	Method	Ours	Amat et al.[8]
	Error type	AEE AAE AEE	AAE
	Translation	0.19	0.26	0.17	0.02
	Rotation+Translation	1.74	0.14	1.93	0.11
	Rotation+Scale	4.30	0.67	7.83	1.17