Mahsa Ghafarianzadeh

Matthew B Blaschko
email: matthew.blaschko@inria.fr

Gabe Sibley
email: gsibley@gwu.edu

Unsupervised Spatio-Temporal Segmentation with Sparse Spectral Clustering

Spatio-temporal cues are powerful sources of information for segmentation in videos. In this work we present an efficient and simple technique for spatio-temporal segmentation that is based on a low-rank spectral clustering algorithm. The complexity of graphbased spatio-temporal segmentation is dominated by the size of the graph, which is proportional to the number of pixels in a video sequence. In contrast to other works, we avoid oversegmenting the images into super-pixels and instead generalize a simple graph based image segmentation. Our graph construction encodes appearance and motion information with temporal links based on optical flow. For large scale data sets naïve graph construction is computationally and memory intensive, and has only been achieved previously using a high power compute cluster. We make feasible for the first time large scale graph-based spatio-temporal segmentation on a single core by exploiting the sparsity structure of the problem and a low rank factorization that has strong approximation guarantees. We empirically demonstrate that constructing the low rank approximation using a subset of pixels (30%-50%) achieves performance exceeding the state-of-the-art on the Hopkins 155 dataset, while enabling the graph to fit in core memory.

Introduction

The analysis of videos is a central task in computer vision. Videos encode useful information in appearance, motion cues, and temporal continuity, which can be leveraged to improve performance on a wide range of applications. Spatio-temporal segmentation is defined as the problem of partitioning a sequence of images into coherent regions using motion and appearance cues. This is an essential step for different tasks such as tracking, object detection and video analysis. Spatio-temporal segmentation is a challenging problem because there are different factors that may change throughout the sequence such as camera motion, scale, illumination, viewpoint, deformation of objects, occlusion, etc.

In this paper, we improve on a family of methods that generalize graph based image segmentation to the spatio-temporal case. Such methods build a graph based on an affinity connecting pixels in a neighborhood and also to their temporal correspondences. This is represented as a very large sparse matrix, from which we select a random subset of columns (which correspond to pixels) -only the randomly selected pixels are used for graph and matrix construction. d) we employ spectral clustering segmentation based on efficient and accurate low rank factorization based on the Nyström method to approximate the graph Laplacian.

matrix with connectivity between neighboring pixels, as well as temporal links computed using optical flow. While such techniques have been applied previously, they have either been based on a sparse set of pixels, super-pixel segmentation, or have made use of extensive computing power to exhaustively compute a dense graph based segmentation. The central contribution of this paper is to introduce a set of strategies that enable us to compute a dense graph based segmentation using a single processor and fitting in core memory. This is done primarily through two innovations: (i) we exploit the sparsity structure of the spatio-temporal graph, and (ii) we make use of an efficient and accurate low rank factorization based on the Nyström method to approximate the graph Laplacian in a spectral clustering approach. We randomly sample a subset of pixels in the sequence of images and build a sparse approximation of the full affinity matrix. The sampled pixels are the nodes of the similarity matrix and edges encode appearance and motion cues. Next we perform spectral clustering on this subset similarity matrix as it provides optimal global solutions and use the solution to segment all the pixels. Our results show that not all of the pixels contain meaningful information about images, and just a subset of pixels can be a good representation of the entire scene.

Related work

A number of video segmentation techniques has been proposed in recent years. Existing video segmentation approaches may differ in the assumptions they make and can be divided into different categories. One category include layered motion segmentation techniques [START_REF] Abhijit S Ogale | Motion segmentation using occlusions. Pattern Analysis and Machine Intelligence[END_REF][START_REF] Smith | Layered motion segmentation and depth ordering by tracking edges[END_REF][START_REF] Xiao | Accurate motion layer segmentation and matting[END_REF] which assume that a scene can be represented with a number of overlapping layers, and each layer has a motion model that can be obtained by some reasoning process and therefore is capable of capturing moving objects in the scene. One drawback is that they may not perform very well on non-textured regions. However, they are capable of segmenting objects when there is occlusion in the scene due to their global motion model.

Another category of video segmentation techniques which is based on sparse number of locations with point trajectories [START_REF] Brox | Object segmentation by long term analysis of point trajectories[END_REF][START_REF] Shankar R Rao | Motion segmentation via robust subspace separation in the presence of outlying, incomplete, or corrupted trajectories[END_REF]. For example in [START_REF] Shankar R Rao | Motion segmentation via robust subspace separation in the presence of outlying, incomplete, or corrupted trajectories[END_REF], a robust subspace separation method has been developed that uses Agglomerative Lossy Compression (ALC) for motion segmentation. The advantage of this method is that it can successfully deal with corrupted or incomplete trajectories where most of other methods in this category fail. In [START_REF] Brox | Object segmentation by long term analysis of point trajectories[END_REF], long term point trajectories based on dense optical flow are clustered into temporally coherent moving segments. However, a major drawback is that these techniques only consider up to 4% of the pixels in the video which is not enough for segmenting the entire video including non-moving objects. Hence, compared to dense video segmentation, the number of labeled points obtained by these sparse methods is considerably smaller.

There are also spatio-temporal segmentation methods which represent the video in a spatio-temporal domain and generalize from image segmentation. For example, in segmentation by hyper-graph cuts [START_REF] Huang | video object segmentation by hypergraph cut[END_REF], first the image is over-segmented into regions, then the regions are taken as the nodes of a hyper-graph that represents the relationship among regions and finally the hyper-graph is partitioned to get the video segmentation. Also in hierarchical graph-based segmentation [START_REF] Grundmann | Efficient hierarchical graph-based video segmentation[END_REF], a volumetric video graph is over-segmented into spatiotemporal regions, and then a region graph is constructed by merging the similar regions obtained from previous over-segmentation. This process is repeated iteratively to create a hierarchical segmentation of the video. Hierarchical mean shift analysis [START_REF] Dementhon | Spatio-temporal segmentation of video by hierarchical mean shift analysis[END_REF], is another method that maps pixels of a video stack into a 7D feature space (three color, two motion and two position component) and then uses mean shift repeatedly to achieve hierarchical clustering of the points in the video. In [START_REF] Levinshtein | Spatiotemporal closure[END_REF], the concept of spatio-temporal closure over a graph of temporal super-pixels has been introduced which is an extension of closure detection techniques in 2D images. Motion aware hierarchical image segmentation (MAHIS) [START_REF] Galasso | Video segmentation with superpixels[END_REF], is another technique that obtains super-pixels from different variety of pixel-based affinities such as brightness,color and texture and then in the next step obtains video segmentation from a variety of super-pixel-based affinities. Another recent work [START_REF] Chang | A video representation using temporal superpixels[END_REF], finds temporally consistent super-pixels in a video using a generative probabilistic model.

Spectral clustering has been used for image and video segmentation, and our work is closely related to [START_REF] Fowlkes | Spectral grouping using the nystrom method[END_REF][START_REF] Shi | Normalized cuts and image segmentation[END_REF]. We aim to partition a graph that encodes appearance and motion similarity in a sequence of images and obtain robust perceptual group which are coherent in time. However, spectral clustering is computationally expensive which makes it often impractical for large scale data set such as spatio-temporal data. In [START_REF] Fowlkes | Spectral grouping using the nystrom method[END_REF], the Nyström method has been applied as a numerical solution for speeding up the spectral clustering. They consider the video as a space time volume and partitioned the volume into coherent segments. In [START_REF] Sundaram | Long term video segmentation through pixel level spectral clustering on gpus[END_REF], they gained more computational power by parallelism in GPUs and perform video segmentation more efficiently. We instead use an algorithm proposed by [START_REF] Li | Time and space efficient spectral clustering via column sampling[END_REF], which is a time and space efficient spectral clustering algorithm for estimating the rank k approximation of normalized Laplacian. This forms the core of our method to perform spatio-temporal segmentation in memory on a single CPU.

Approach

Creating the Graph

Given an image I, we create a graph G = (V, E,W ), where the graph nodes V are the pixels in the image and are connected by edge E if they are within distance r from each other. As discussed in [START_REF] Cour | Spectral segmentation with multiscale graph decomposition[END_REF], in general choosing a smaller r makes the segmentation faster and choosing a larger r helps with propagating local grouping cues across larger neighborhood which is useful for detecting faint contours in the image. W measures the similarity of pixels connected by an edge. We define W as the following:

W i j = exp -d 2 (s i , s j ) σ i σ j
where W i j = 0 for i =j, and s i denotes pixel color and d(s i , s j ) is the Euclidean distance. σ is a local scaling parameter [START_REF] Zelnik | Self-tuning spectral clustering[END_REF] which takes into account the local statistics of the neighborhood around pixels i and j. Local scaling parameter is defined by:

σ i = d(s i , s k )
where s k is the K'th neighbor of pixel i. In order to extend this to video, we make use of optical flow and add temporal motion information to the graph. We use optical flow [START_REF] Liu | Beyond pixels: exploring new representations and applications for motion analysis[END_REF] to compute the motion vectors between frames. Then we connect pixel (x, y) in frame t to its 9 neighbors along the backward flow (u, v) in frame t -1, e.g. (x + u(x, y) + δ x , y + v(x, y) + δ y ) for δ x , δ y ∈ {-1, 0, 1}. The similarity matrix for the video is a sparse symmetric block diagonal matrix of the size n = number of frames × number of pixels in one frame.

Partitioning the Graph

Spectral Clustering

Spectral clustering has been widely used for partitioning the similarity graphs. Given the similarity matrix of n points,W , first the normalized similarity matrix L = D -1 2 W D -1 2 is computed, where D is the diagonal degree matrix defined as D ii = ∑ n j=1 W i j . Then it finds the k largest eigenvectors of L and form V = [v 1 , ..., v k ] ∈ R n×k . Next V will be re-normalized row-wise to have unit length, U, and then k-means is applied to cluster the rows of U. Finally, each pixel s i will be assigned to the corresponding cluster for row i of U. However, the O(n 3 ) complexity of eigen decomposition of L and also storing the similarity matrix often makes it expensive and impractical when n is large.

Nyström Method

The Nyström method used by [START_REF] Fowlkes | Spectral grouping using the nystrom method[END_REF] is a popular technique to tackle the mentioned problem by approximating a low-rank matrix of W , since we only need a number of eigenvectors. Given a symmetric matrix W ∈ R n×n , m ≪ n columns of the matrix will be sampled. For simplicity, we assume that these sampled columns are the first m columns of W . Matrix W can be represented as follows:

W = A B B T C
where A ∈ R m×m , B ∈ R n-m×m , and C ∈ R n-m×n-m . The Nyström approximation of W is defined by:

Ŵ = A B T A -1 A B
Then if A is positive definite, we can define S = A + A -1 2 BB T A -1 2 and perform singular value decomposition (SVD) on S to get S = U S Λ S U S T . The approximated L then can be defined as L = V ΛV T and V T V = I, where

V = A B T A -1 2 U S Λ S -1 2
But if A is indefinite, a more complicated process is required. We define Z = ŪS Λ 1 2 and ŪS T = U T S Λ -1 S U T S B so L = ZZ T . If singular value decomposition of ZZ T = FΣF T , then V = ZFΣ -1 2 has the orthogonalized eigenvectors of L. This algorithm is expensive for large scale data sets, since the time complexity of this algorithm is O(nm 2 + m 3 ) and it requires keeping A and A B T in the memory. Hence, we are using a time and space efficient spectral clustering via column sampling [START_REF] Li | Time and space efficient spectral clustering via column sampling[END_REF], that is similar to Nyström method, but with a further rank-k approximation of the normalized Laplacian using the sampled sub-matrix of the similarity matrix. This algorithm has shown promising results, since it reduces the time and space complexity of Nyström method and also it is able to recover all the degree information of the selected points. Here we briefly review the algorithm proposed in [START_REF] Li | Time and space efficient spectral clustering via column sampling[END_REF]. First we form the sampled sub-matrix A ∈ R m×m and compute the degree matrix D * for A. Next, we define M * ∈ R m×m as

M * = D -1 2 * AD -1 2 *
and S * ∈ R m×m based on the rank-k approximation M * ,k of M * as

S * = D -1 2 * M -1 * ,k D -1 2 *
M * ,k is computed by performing eigen decomposition of M * and find the k largest eigenvalues Λ and the corresponding eigenvectors V . Finally the approximated Ŵ is obtained by

Ŵ = A B T S * A B T
T By having M * ,k = V ΛV T and replacing it in above formulas, we get Ŵ = QΛQ T , where

Q = A B T D -1 2 * V Λ -1
Lastly, rank-k approximation of L is obtained by following

Lk = D-1 2 Â D-1 2 = UΛU T
where U = D-1 2 Q and D is the degree matrix of Ŵ . The time complexity of this algorithm is O(nmk) and there is no need to store large similarity matrix W or its sampled columns in the memory. Also we are using the proposed inexpensive algorithm by [START_REF] Li | Time and space efficient spectral clustering via column sampling[END_REF] to orthogonalized estimated eigenvectors U.

Post Processing

After performing spectral clustering on the similarity graph and obtaining the clusters, we first quantize each cluster into 256 bins (16 bins for each channel) and compute the RGB histogram. Then we find the Euclidean distance between adjacent clusters and merge them repeatedly if their similarity is more than a threshold τ to achieve the final segmentation.

Results

To evaluate the performance of our segmentation method, we run several experiments and provide quantitative and qualitative results. We compare our method against motion segmentation techniques such as [START_REF] Brox | Object segmentation by long term analysis of point trajectories[END_REF][START_REF] Shankar R Rao | Motion segmentation via robust subspace separation in the presence of outlying, incomplete, or corrupted trajectories[END_REF], which cluster point trajectories in the sequence of images and also other dense video segmentation techniques such as [START_REF] Grundmann | Efficient hierarchical graph-based video segmentation[END_REF][START_REF] Sundaram | Long term video segmentation through pixel level spectral clustering on gpus[END_REF]. In all of our experiments, we randomly select a subset sample of size 30%-50% for approximating the similarity matrix. The quality of our segmentation highly relies on the selected number of clusters for k-means and also distance r for computing pixel similarities. Choosing larger numbers for these parameters makes the overall segmentation better. However, the disadvantage is that it creates unsmooth boundaries and captures most of the undesired details in the scene and also increases the size of matrix rapidly. We found that by choosing a distance size of 11-15 and cluster size of 300-500, we achieve reasonable performance.

We use an annotated dataset by [START_REF] Brox | Object segmentation by long term analysis of point trajectories[END_REF] for motion segmentation. This dataset includes 26 sequences which is annotated sparse in time and dense in space. We ran our method on three scenarios and show the numerical result in Table 1. The first section shows the segmentation result for the first 10 frames of all 26 available sequences in the dataset, the second section shows the segmentation result for the first 50 frames of 15 selected sequences in the dataset and the third scenario shows the result for the first 200 frames of 7 selected sequences. The density shows the percentage of pixels for which a label is reported. The overall error is the number of mislabeled pixels over the total number of labeled pixels. the overall accuracy is the number of correctly labeled pixels over the total number of pixels. The average error is similar to the overall error but computed separately for each region first and then the average across all the regions is reported. The over-segmentation error indicates the number of clusters that needs to be merged to fit the ground truth regions. Also the number of regions with less than 10% error is reported as the extracted objects.

In this experiment, we aimed to compare our method against the method in [START_REF] Sundaram | Long term video segmentation through pixel level spectral clustering on gpus[END_REF]. The reason is that in contrast with their method which achieves good performance by using GPUs and heavily depending on parallelization, we exploit numerical algorithms to make our method fast and efficient. The density of both methods is 100% but we achieved comparable errors while using just a subset of pixels (30%-50%) to label all of the pixels as Figure 2: Comparison of our technique vs motion segmentation by [START_REF] Sundaram | Long term video segmentation through pixel level spectral clustering on gpus[END_REF] and [START_REF] Grundmann | Efficient hierarchical graph-based video segmentation[END_REF]. Left-toright shows video frames 1 and 50. Second row: results from Hierarchical graph segmentation [START_REF] Grundmann | Efficient hierarchical graph-based video segmentation[END_REF]; third segmentation with the full graph using a cluster of GPUs [START_REF] Sundaram | Long term video segmentation through pixel level spectral clustering on gpus[END_REF]; fourth row, our result runs on a single core, does not employ contour smoothing or texture similarity. compared with their method that uses all the pixels for clustering. Compared to the motion segmentation methods [START_REF] Brox | Object segmentation by long term analysis of point trajectories[END_REF][START_REF] Shankar R Rao | Motion segmentation via robust subspace separation in the presence of outlying, incomplete, or corrupted trajectories[END_REF], our over-segmentation error is higher but that is due to the fact that these methods only use sparse number of points to label about 4% of the data while we are labeling all the pixels and achieve better accuracy.

Discussion

There are several avenues for future research based on the methods presented in this work. Sparse spectral clustering algorithms are designed to be agnostic to the problem domain. However, we have a substantial amount of domain knowledge we can exploit to expedite the clustering.

The matrix factorization employed uses uniform sampling to select columns. However, selective sampling strategies generally decrease the approximation error for a given number of columns. As the error rate of our method is comparable to using 100% of the columns, this indicates that we may achieve an even higher savings in memory using selective sampling. An additional likely source of improvement is through the use of similarity measures based on feature functions that encode texture or other information sources in addition to color.

Conclusions

In this work, we have demonstrated a novel method for spatio-temporal segmentation of dense pixel trajectories based on spectral clustering. We found that fully connecting pixels to their spatial neighbors within a given radius is an effective strategy for improving segmentation accuracy. Additionally, we use optical flow to more accurately compute temporal connectivity than a simple method based on an interpretation of the video sequence as a 3D volume. In contrast to previous work, we do not resort to super-pixel segmentation to achieve computational tractability and memory efficiency. Instead, we exploit the natural sparsity structure of the graph, and employ a low rank approximation of the Laplacian closely related to the Nyström method. We have found that sampling 30-50% of the pixels to index columns of the low rank approximation leads to comparable performance with a method that uses 100% of the columns. This strategy results in a spectral clustering method that can run on a single processor with the graph representation fitting in core memory. For instance, a 100 frame 640x480 sequence naively requires 6.8GB of memory, while with our method we are able to reduce this to 2.1GB and still retain accuracy. We have demonstrated the effectiveness of the approach on the Hopkins 155 data set, where we have achieved the best reported results for dense segmentation using an order of magnitude less computation than [START_REF] Sundaram | Long term video segmentation through pixel level spectral clustering on gpus[END_REF].

Figure 1 :

 1 Figure 1: Approach overview: given a sequence of images a) we construct a graph b)connecting pixels in a neighborhood and also to their temporal correspondences. This is represented as a very large sparse matrix, from which we select a random subset of columns (which correspond to pixels) -only the randomly selected pixels are used for graph and matrix construction. d) we employ spectral clustering segmentation based on efficient and accurate low rank factorization based on the Nyström method to approximate the graph Laplacian.

Figure 3 :

 3 Figure 3: Results form the Army Men sequence. Notice the fine detail around the top-right gun figure.

Figure 4 :

 4 Figure 4: Results from Ice-skater Yu-Na Kim, 2009 World Championships, c 2009 NBC Olympics.

Table 1 :

 1 Comparison of dense and sparse segmentation methods

	Name	Density Overall Error Overall Accuracy Average Error Over-segmentation Extracted Objects
				First 10 frames (26 sequences)
	Ours (Sample:50%)	100%	5.06%	94.94%	22.31%	6.12	22
	Sundaram et al[16]	100%	6.97%	93.03%	19.72%	7.15	20
	Brox et al[1]	3.34%	7.75%	3.08%	25.01%	0.54	24
	Rao et al[13] corrupted	2.98%	7.88%	2.74%	24.05%	0.15	26
	Rao et al[13] incomplete 3.34%	11.20%	2.97%	26.73%	0.54	19
				First 50 frames (15 sequences)
	Ours(Sample:50%)	100%	7.63%	92.37%	32.24%	14.9	5
	Sundaram et al[16]	100%	8.42%	91.58%	28.76%	13.0	5
	Brox et al[1]	3.27%	7.13%	3.04%	34.76%	0.53	9
	Rao et al[13] corrupted	1.53%	7.91%	1.41%	42.13%	0.36	8
	Rao et al[13] incomplete 3.27%	16.42%	2.73%	49.05%	6.07	2
				First 200 frames (7 sequences)
	Ours(Sample:50%)	100%	13.24%	86.76%	47.83%	5.3	6
	Sundaram et al[16]	100%	15.46%	84.54%	45.05%	5.86	3
	Brox et al[1]	3.43%	7.64%	3.17%	31.14%	3.14	7
	Rao et al[13] corrupted	0.20%	0.00%	0.20%	74.52%	0.40	1
	Rao et al[13] incomplete 3.43%	19.33%	2.77%	50.98%	54.57	0
	Grundmann et al[7]	100%	20.77%	79.23%	40.37%	10.42	0


Arnaud Tonnelier

Hana Belmabrouk

Dominique Martinez

Event-driven simulations of nonlinear integrate-and-fire neurons

Event-driven strategies have been used to simulate spiking neural networks exactly. Previou work is limited to linear integrate-and-fire neurons. In this note we extend event driven schemes to a class of nonlinear integrate-and-fire models. Results are presented for the quadratic integrate-and-fire model with instantaneous or exponential synaptic currents. Extensions to conductance-based currents and exponential integrate-and-fire neurons are discussed.

Introduction

Our current theoretical understanding of the properties of neural systems is mainly based on numerical simulations from single cell models to neural networks. Recent experimental evidence has accumulated which suggests that precise spike-time coding is used in various neuronal systems [START_REF] Vanrullen | Spike times make sense[END_REF]. Individual spikes can be highly reliable and precisely timed with a submillisecond accuracy [START_REF] Berry | The structure and precision of retinal spike trains[END_REF][START_REF] Mainen | Reliability of spike timing in neocortical neurons[END_REF]. Synaptic plasticity depends critically on the relative timing of pre-and postsynaptic spikes; a synapse is strengthened if the presynaptic spike occurs shortly before the postsynaptic neuron fires, and the synapse is weakened if the sequence of spikes is reversed [START_REF] Bi | Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type[END_REF]. It is therefore important to have accurate numerical schemes to calculate spike times. Integrate-and-fire neurons reproduce many features of the neuronal dynamics [START_REF] Izhikevich | Simple model of spiking neurons[END_REF][START_REF] Tonnelier | Piecewise linear differential equations and integrate-and-fire neurons: insights from two-dimensional membrane models[END_REF][START_REF] Gerstner | Spiking neuron models[END_REF] and are widely used in the numerical simulations of spiking neural networks. Two strategies have been used for the simulation of integrate-and-fire neural networks: time-stepping methods that approximate the membrane voltage of neurons on a discretized time and event-driven schemes where the timings of spikes are calculated exactly. By definition time-stepping approximations are imprecise and it has been shown that time steps have to be chosen correctly to reproduce the synchronization properties of networks of spiking neurons [START_REF] Hansel | On the numerical simulations of integrate-and-fire networks[END_REF]. Standard time-stepping algorithms (Euler, Runge-Kutta) have to be modified to give an accurate approximation of firing times [START_REF] Shelley | Efficient and accurate time-stepping schemes for integrate-and-fire neuronal networks[END_REF]. A fundamental limitation on the accuracy of any time-stepping methods is imposed by the smoothness of the postsynaptic potentials [START_REF] Shelley | Efficient and accurate time-stepping schemes for integrate-and-fire neuronal networks[END_REF]. High-order time-stepping algorithms can be constructed only if the onset of postsynaptc conductance changes has smooth derivatives. Exact simulations avoid these problems. The term exact method or exact simulation means that spike timings are analytically given or are derived from the analytical formulation of the membrane potential. Thus it is possible to have an arbitrary precision (up to the machine precision) of spike timings. This method has become increasingly popular [START_REF] Mattia | Efficient event-driven simulation of large networks of spiking neurons and dynamical synapses[END_REF][START_REF] Makino | A discrete-event neural network simulator for general neuron models[END_REF][START_REF] Rochel | An event-driven framawork for the simulation of networks of spiking neurons[END_REF][START_REF] Brette | Exact simulation of integrate-and-fire models with synaptic conductances[END_REF][START_REF] Rudolph | Analytical integrate-and-fire neuron models with conductance-based dynamics for event-driven simulation strategies[END_REF] but it applies to a limited class of neuron models, mainly the linear integrate-and-fire models. It is known that the leaky integrate-and-fire model has some limitations: it has an unrealistic behavior close to the threshold and reproduces only some characteristics of conductance based neuron [START_REF] Fourcaud-Trocmé | How spike generation mechanisms determine the neuronal response to fluctuating inputs[END_REF]. More realistic models include nonlinear spike-generating currents that allow replacement of the strict voltage threshold by a smooth spike initiation zone. Here we simulate exactly the quadratic integrate-and-fire neuron [START_REF] Ermentrout | Parabolic bursting in an excitable system coupled with a slow oscillation[END_REF]. In the quadratic integrate-and-fire (QIF) model, the membrane potential follows

C dV dt = q(V -V th ) 2 -I th + I s (t), ( 1
)
where V is the membrane voltage, C is the membrane capacitance, q is a parameter characterizing the frequency-current response curve, I th is the threshold current and I s is the synaptic current and V th is the voltage threshold, that is, the largest steady voltage at which the neuron can be maintained by a constant input. Without synaptic currents, I s = 0, the QIF neuron presents two distinct regimes. When I th > 0 there are two fixed points. The stable one defines the resting state V rest of the neuron

V rest = V th - I th q .
(2)

The unstable one is the threshold below which trajectories tend towards infinity in finite time that defines the spike time. When I th < 0 the neuron fires regularly. This model represents the normal form of any type I neurons near the saddle-node bifurcation [START_REF] Ermentrout | Type i membranes, phase resetting curves, and synchrony[END_REF][START_REF] Ermentrout | Parabolic bursting in an excitable system coupled with a slow oscillation[END_REF] and is related to the so-called θ neuron [START_REF] Ermentrout | Type i membranes, phase resetting curves, and synchrony[END_REF][START_REF] Gutkin | Dynamics of membrane excitability determine interspike interval variability: a link between spike generation mechanisms and cortical spike train statistics[END_REF]. Since the quadratic integrate-and-fire neuron is expected to reproduce the characterisitcs of any type I neuron close to bifurcation, it has been widely used as a realistic neuron model [START_REF] Hansel | Existence and stability of persistent states in large neuronal networks[END_REF][START_REF] Latham | Intrinsic dynamics in neuronal networks, i. theory[END_REF][START_REF] Brunel | Firing rate of the noisy quadratic integrateand-fire neuron[END_REF][START_REF] Fourcaud-Trocmé | How spike generation mechanisms determine the neuronal response to fluctuating inputs[END_REF]. The action potential is defined as a divergence of the voltage. In numerical simulations, one has to introduce a cutoff at a finite voltage V peak . After a spike, the membrane potential is instantaneously reset to V reset . The synaptic current I s is induced by presynaptic spikes. An incoming spike at time t f triggers a postsynaptic current

I s (t) = w δ(t -t f ), ( 3
)
where w is the synaptic weight of the synapse and δ is the dirac delta function.

More realistic models use exponential currents

I s (t) = w exp(-(t -t f )/τ s ), (4)
where τ s is the synaptic time constant. An event-driven simulation requires an analytical expression for the membrane potential or, at least, a closed form expression for the firing times. In many cases a given neuron in the network will not fire (its firing time is +∞). In such situation the efficiency of the simulation is improved using a spike-test, i.e. an algorithm that checks quickly whether a neuron will fire. Note that the efficiency of the spike-test crucially depends on the overall activity of the network.

In this note, we describe a method to simulate exactly the quadratic integrateand-fire model (1) with the synaptic currents (3)-( 4). The membrane potential is solved analytically and a spike-test is derived.

Exact simulation

We consider the dimensionless QIF model

dV dt = V 2 -I th , (5)
obtained from (1) by the change of variables V ← q(V -V th )/C and I th ← q I th /C 2 (for convenience we do not change the notations). We further consider an instantaneous synaptic current so that V → V + w when a spike is received, V → V reset when a spike is emitted. We consider excitable neurons and we take a positive threshold current. In an interval with no spike (i.e., no incoming spike and no spike emitted by the neuron itself), equation ( 5) can be solved analytically with the separation-of-variables technique. The membrane potential evolves according to

V (t) = -I th tanh I th (t + c) , (6)
where c is a constant of integration given by the initial condition V (0) = -√ I th tanh(c √ I th ). For V (0) < √ I th the neuron goes back to its resting value and no spike is emitted. A spike is emitted for V (0) > √ I th , i.e. for an initial value greater than the unstable fixed point of the neuron. In this case, we calculate

1 c √ I th = -atanh( √ I th /V (0)) + iπ/2. Using tanh(x + iπ/2) = coth(x) we have V (t) = -I th coth( I th t -atanh( I th /V (0))). ( 7
)
The firing time is obtained when V (t) crosses V peak and is explicitly given using (7). In the limit V peak → ∞ the firing time is obtained equaling to zero the argument of the coth function and is simply given by

t f = 1 √ I th atanh √ I th V (0) . ( 8
)
An event-driven scheme for a network of QIF neurons with instantaneous coupling can be easily implemented using (8). The generalization of the event-based simulation to more realistic synaptic currents is not trivial.

A network of QIF neurons with exponential current follows the equations

dV dt = V 2 -I th + I s , (9)
τ s dI s dt = -I s . ( 10
)
When a spike is received by a synapse, the current I s is instantaneously modified I s → I s + w. Excitatory or inhibitory synapses are accounted through the sign of w provided that they have an identical time constant τ s . System (9-10) is analytically solvable in intervals with no spike and V (t) is expressed using Bessel functions (see Appendix 5.1). Thus we can calculate very accurately the value of the membrane potential at any time.

It is possible to speed up the simulation if neurons that will not spike are detected. We will derive a test that quickly check whether a neuron will spike.

Depending on the initial conditions V (0) ,I s (0) either a spike is emitted or the neuron goes back directly to its resting state. There is a curve V * (I s ) in the phase plane that defines the threshold curve, i.e. the neuron spikes if and only if its state is above this curve. The threshold curve is given by the stable variety of the unstable fixed point (I s , V ) = (0, √ I th ) (see Fig. A). When the spike test is positive we need to compute the firing time. This computation can be done very efficiently with standard methods of roots finding. These methods require an initial guess. We define V + and V -that evolve according to the following differential equations

dV + dt = V peak |V + | -I th + I s , (11)
dV - dt = 2 |V -| -2 -I th + I s , (12)
where I s follows (10). Here, we have used piecewise linear bounds of the quadratic nonlinearity (Fig. B, left). is an arbitrary parameter and we choose = 1. By construction we have, ∀t ≥ 0, V -(t) ≤ V (t) ≤ V + (t) (see Eqs. 9,11,12 and Fig. B,right). Since V + (t) and V -(t) are the membrane potentials of (piecewise) linear integrate-and-fire models we can calculate accurately the times of the threshold crossing [START_REF] Brette | Exact simulation of integrate-and-fire models with synaptic conductances[END_REF]. These times give an upper and a lower bound for the firing times of the QIF neuron that allow an efficient root finding. However when the network activity is high, the spike-test is most of the time positive and this method becomes time consuming. Moreover in the high activity regime, V (0) is near the threshold curve and it is possible to derive a more accurate bound of the firing time. Another approximation as initial guess is obtained using the bounding QIF

dV + dt = V + 2 -I th + max(0, I s (0)), (13)
that gives a lower bound for the firing time. Likewise, an upper bound can be obtained by using min(0, I s (0)) in the equation above. The firing time of the bounding QIF model ( 13) can be calculated quickly and accuratly using an equation similar to (8). The comparison with the lower bound given by the bounding LIF is shown figure C for different values of the synaptic current. We calculate the relative error E = |t f -t f b |/t f where t f is the firing time of the QIF and t f b is the approximation given by either the bounding LIF (11) or the bounding QIF (13), respectively. There exists a critical value Ṽ0 (I s (0)) of the membrane potential such that for V > Ṽ0 the bounding QIF model gives a better approximation of the firing time. Note that for positive values of V 0 the bounding QIF model is always better. This case is achieved when the average activity is high, i.e. neurons in the network are most of the time near the threshold curve.

Numerical results

To illustrate our method, we simulate a network of N identical QIF neurons receiving an external excitatory spike train modeled as a Poisson process with a constant rate of 10 kHz2 . The QIF neurons are coupled all-to-all through GABAergic inhibition. This scenario has been used to reproduce the oscillatory synchronization observed in early olfactory systems [START_REF] Martinez | Oscillatory synchronization requires precise and balanced feedback inhibition in a model of the insect antennal lobe[END_REF][START_REF] Ambard | Inhibitory control of spike timing precision[END_REF]. Numerical values of QIF neurons (Eq. 1) are taken from [START_REF] Ambard | Inhibitory control of spike timing precision[END_REF]: C = 0.2 nF, V rest = -65 mV, V th = -60.68 mV, q = 0.00643 mS.V -1 , I th = 0.12 nA, V peak = 30 mV and V reset = -70 mv. The synaptic time constant is τ s = 6ms and the synaptic strength is w = 0.05/N for the inhibitory synapses and w = 5.10 -5 for the excitatory poissonian input. We wrote our event-driven simulator in C++ based on the event-driven library MVASpike [START_REF] Rochel | An event-driven framawork for the simulation of networks of spiking neurons[END_REF] available at http://mvaspike.gforge.inria.fr. The Bessel functions needed for the exact computation of V (t) (see appendix 5.1) were implemented by using the GNU Scientic Library (http://www.gnu.org/software/gsl). The threshold curve defining the spike test was very accurately fitted by a polynomial of degree three. Each time the spike test is positive, the firing time is computed by finding the root of V (t) -V peak . This is accomplished by a few steps of a Newton-Raphson method starting from an initial guess obtained from the bounding QIF (Eq. 13). The precision achieved in the simulations is in the order of 10 -7 ms. We first simulate a network of N = 100 neurons. Initial conditions of the membrane potentials V (0) are taken randomly between V reset and V peak so that the firing times of uncoupled neurons are uniformly distributed. This way, the network starts in a completely desynchronized state. When inhibition is blocked, individual neurons fire at 380 Hz on average (see Figure D, left for the spike raster of the first 100 ms of the simulation). The simulation time was about 20 min for 1 sec of biological time on a portable PC running Linux at 1.86 GHz. This is a reasonnable time in regards to the large number of events encountered during the simulation (in the order of 10 6 ). Moreover, we found that the spike test is positive 100% of the time so that any event from the poissonian input requires spike timing computation of the N QIF neurons. When lateral inhibition is taken into account, the spike test is positive in 50% of the cases and neurons are synchronized with precise spiking activity at about 10Hz (see Figure D, right). In such a case, all neurons fire almost at the same time. Rather than searching for the next firing time within the entire population of neurons at each occurance of an event, it might be sufficient to perform the search among a random subset of k neurons. This leads to a considerable speedup as k << N , albeit at a cost of introducing errors in the simulation. A similar probabilistic speedup has been used in kernel methods [START_REF] Scholkopf | Learning with kernels[END_REF]. In Appendix 5.3, it is shown that the error probability3 is given by 0.95 k for any N . Therefore, a subset of 59 neurons is already sufficient to get a low error probability. For k = 59, the error probability equals 0.05 which is small enough to not perturb the synchronization as shown in Figure E, left for a network of N = 1000 neurons coupled all-to-all (10 6 synapses). For comparison, Figure E, right shows the dynamics obtained by simulating the same network with k = 30 (error probability of 0.21).

Discussion

Due to its simple dynamics the behavior exhibited by the leaky integrate-andfire model is limited. The QIF model has seen increasing interest in recent years, primarily because it reproduces the properties of detailed conductancebased neurons near the threshold. At low firing rate, the response of any type I neuron is described by the QIF neuron. In particular the QIF model reproduces the f-I curve of any type I neurons near the threshold. We have proposed a method to simulate the QIF neuron in an event driven way. We have presented simulations of the quadratic integrate-and-fire model with exponential synaptic currents.

In the numerical simulation of neural networks it is difficult to predict if errors on firing times will create numerical artefacts or will remain irrelevant. Event-driven methods are more precise than traditional time-stepping integration algorithms because the spike timings are calculated with an arbitrary desired precision. In general, implementation of an event-driven method require non-linear root-finding algorithms, for instance the Newton-Raphson method which converges exponentially. Thus calculations are exact in the sense that the limit of the computer accuracy is reached with few iterations of the method. The error done with Euler or Runge-Kutta methods is fundamentally different. It is inherent to these approximated methods and the precision is fixed with a polynomial decreasing with respect to the time-step. However when the number of events is large, event-driven schemes are time-consuming. This issue is implementation dependent and, in particular, depends on how the event queue is performed. Another important limitation of exact simulations is that it cannot be applied to any model. Implementation of exact simulations is possible if the membrane equations are analytically solvable. Extensions of our event driven scheme to other types of synaptic connections or spike-generating currents have to be studied case-by-case :

• A QIF neuron with exponential conductance-based synaptic currents follows

dV dt = V 2 -I th + g(V -E s ), τ s dg dt = -g,
where g and E s is the synaptic conductance and the synaptic reversal potential, respectively. The membrane voltage is analytically solvable and can be written using Whittaker functions (see Appendix 5.2). Thus, an event-driven scheme, similar to the one proposed in this note, can be used to simulate exactly QIF neurons with synaptic conductances.

• The recently proposed exponential integrate-and-fire (EIF) model [START_REF] Fourcaud-Trocmé | How spike generation mechanisms determine the neuronal response to fluctuating inputs[END_REF][START_REF] Brette | Adaptative exponential integrate-and-fire model as an effective description of neuronal activity[END_REF] is given by

dV dt = -(V + E) + e V , (14)
where V and t are normalized variables and E = (V T -E l )/∆ T with V T a threshold voltage, E l the leak potential, ∆ T the spike slope factor. Like the QIF model, the EIF model has a soft threshold. However the spike-generating current is no longer quadratic but exponential. It has been shown that the f-I curve of the EIF model matches the f-I curve of detailed conductance-based models for a range of input currents larger than for the QIF model. The EIF model with instantaneous current (3) has an implicit solution and the firing time is given by

t f = V peak V (0) du -(u + E) + e u ,
for V (0) greater than the unstable fixed point of ( 14). As for the QIF neuron, V peak controls the shape of the spikes. It is possible to use precalculated tables with an arbitrary precision for the integral for an eventdriven simulation of the EIF model. We do not know if the EIF model with exponential currents is analytically tractable or if an implicit solution can be found.

Appendices

The QIF with exponential synaptic current

The QIF neuron with an exponential synaptic current can be rewritten as a nonlinear and nonautonomous differential equation

dV dt = V 2 -I th + I s (0)e -t/τs . ( 15
)
The first step is to transform (15) into a linear ODE. We use the change of variables

V (t) = - 1 y(t) dy dt ,
to obtain d 2 y dt 2 = (-I th + I s (0)e -t/τs )y(t). ( 16)

We use the change of variables u(t) = 2τ s I s (0)e -t/(2τs) which transforms the equation into

u 2 d 2 y du 2 + u dy du + (u 2 -4τ 2 s I th )y = 0. ( 17
)
Let J α (x) and Y α (x) be the two independant solutions of the so-called Bessel's equation

u 2 d 2 y du 2 + u dy du + (u 2 -α 2 )y = 0,
where α = -2τ s √ I th . We have

y(u) = J α (u) + cY α (u),
where c is a constant defined by the initial condition. Using

dJ α (u) du = -J α+1 (u) + α u J α (u),
and a similar equation for Y α we find that the membrane potential of the QIF neuron with exponential synaptic current (9-10) is analytically solvable. The membrane potential is given by

V (t) = -I th - u(t) 2τ s J α+1 + cY α+1 J α + cY α .
where J and Y are Bessel functions evaluated at u(t) = 2τ s I s (0)e -t/(2τs) . The Bessel functions can be expressed as a series of Gamma functions. Note that calculations with infinite series are not fundamentally different from calculating non-polynomial functions (like the exponential or logarithmic functions). Numerically, the first terms of the serie are necessary to reach a good precision.

It is also possible to use precalculated tables to speed up the calculations.

The QIF with exponential synaptic conductance

The analytical integration of the QIF with exponential synaptic conductance proceeds along the same lines as the integration of the QIF with exponential currents. The analytical expression of the membrane potential uses Whittaker functions instead of Bessel functions. More precisely, let W µ,ν (x) and M µ,ν (x), the Whittaker functions, be the solutions of the differential equation

d 2 y dx 2 + - 1 4 + µ x + 1/4 -ν 2 x 2 y = 0.
The membrane potential of the QIF model with exponential synaptic conductance is given by

V (t) = 1/(cW α,β + M α,β ) * (E s M α,β + ( I th -E s )M α-1,β + . . . +cE s W α,β + cτ s (E 2 s -I th )W α-1,β
), where W and M are Whittaker functions evaluated at u(t) = τ s I s (0)e -t/τs , α = 1/2 + τ s E s , β = τ s √ I th and c is a constant defined by the initial condition. Whittaker functions have power series expansions that could be used to enhance the efficiency of the calculations.

The probabilistic speedup

Advancing the event-driven simulation to the next firing time requires the computation of the firing times of the N neurons. This can be very time-consuming as N is large. In case of synchronization, the neurons fire almost at the same time and it might be sufficient to compute the firing times t i , i = 1 • • • k, among a random subset of k neurons and advance the simulation to the next firing time min(t i ) obtained from this subset. This leads to a considerable speedup as k << N , albeit at a cost of introducing errors in the simulation. We define the error probability p e as the probability that min(t i ) is greater than the 5% first firing times of the entire population of neurons. We then have p e = P (min(t i ) > t ) where t is the maximum of the 5% first firing times of the N neurons. Considering the t i s as identically and independently distributed firing times, we have p e = (P (t i > t ))

k . Moreover, P (t i < t ) = 0.05 because we have considered the 5% first firing neurons in the definition of t , and thus p e = 0.95 k .

Figure

  Figure. A. The (I s , V )-phase plane for the QIF neuron with exponential synaptic currents. The V -nullcline is shown (bell-shape curve). The stable manifold of the saddle-node point (0, √ I th ) defines the threshold curve above which the neuron fires. Trajectories starting under the threshold curve tend toward the stable resting state (0, -√ I th ). B. (left) Piecewise linear sector bounds on the nonlinearity. The nonlinearity is shown in full line. The dashed lines show piecewise linear sector bounds. (right) Bounding of the membrane potential of the quadratic integrate-and-fire neuron using piecewise linear integrate-and-fire models. Parameters are V peak = 3, I th = 2. C. The relative error on the spike timing of the QIF as a function of the initial value of the membrane potential V 0 for different values of the initial synaptic current. Comparison is done between the approximation obtained by the bounding LIF (bLIF) and by the bounding QIF (bQIF). (left) I s (0) = 5. (middle) I s (0) = 7. (right) I s (0) = 10. Other parameters are those of B. D. Event-driven simulations of a network of QIF neurons (N = 100 neurons). (left) Spike trains produced without inhibition. (right) Spike trains produced with lateral inhibition (10 4 synapses). E. Event-driven simulations of a fully connected network of QIF neurons (N = 1000 neurons, 10 6 synapses). (left) Simulation with probabilistic speedup (k = 59). (right) Simulation with probabilistic speedup (k = 30).

Depending on the initial condition V (0), the constant c can be a complex number

This mimics the interaction with 1000 excitatory presynaptic neurons firing at 10 Hz on average

Error probability is defined as the probability that the firing time computed with the probabilistic speedup is greater than the 5% first firing times of the entire population of neurons.

Acknowledgments

Research supported by the INRIA cooperative research initiative RDNR.


Olivier Chesneau

Florentin Millour

Orsola De Marco

S N Bright

Alain Spang

D P K Banerjee

N M Ashok

T Kaminski

John P Wisniewski

Anthony Meilland

T Kamiński

E Lagadec

V838 Monocerotis: the central star and its environment a decade after outburst

Keywords: high angular resolution, individual: V838 Mon, Stars: circumstellar matter, Stars: mass-loss

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

V838 Mon is a nova-like object, which erupted in 2002 in a series of outbursts (reaching V∼6.8). It subsequently developed a light echo which was extensively studied by Bond et al. (2003); Banerjee et al. (2006); Tylenda & Kamiński (2012). The eruption was unlike classical novae as the effective temperature of the object dropped and the spectral type evolved into a very late Ltype supergiant (Evans et al. 2003;Loebman et al. in press). The event has been interpreted as the merger of a ∼8 M ⊙ star with a sub-solar mass star (Soker & Tylenda 2003;Tylenda et al. 2005;Tylenda & Soker 2006). The mid infrared flux of V838 Mon increased by a factor of two between 2004 and 2007, suggesting that new dust was forming in the expanding ejecta of the outbursts (Wisniewski et al. 2003a;Tylenda 2005;Wisniewski et al. 2008). The expanding ejecta engulfed a companion close to the central source (Bond 2006;Munari et al. 2007;[START_REF] Kolka | The Astronomer's Telegram[END_REF].

Send offprint requests to: Florentin.Millour@oca.eu ⋆ Based on observations at Paranal Observatory under programs 088.D-0005, 090.D-0011, 091.D-0030 and 093.D-0056.

O. Chesneau passed away shortly before submitting this letter. We express our profound sadness on this untimely demise and convey our deepest condolences to his family.

The progenitor of V838 Mon was found in a variety of surveys, including 2MASS, and appears consistent with a somewhat reddened early main sequence star. Post-outburst observations have found a faint, blue component in the spectrum and Tylenda et al. (2005) argued that the pre-outburst spectral energy distribution (SED) is well-matched by a pair of early main sequence stars (B3V + B1.5V or B4V + A0.5V), making V838 Mon a triple system de facto (main star + sub-solar merging star + B3V companion). Lynch et al. (2004) obtained numerous IR spectra of V838 Mon during and after the eruption (2002)(2003). They fit their data to a model consisting of a cool (T eff = 2100K, R * = 8.8 AU) stellar photosphere surrounded by a large, absorbing molecular cloud. A more recent work (Loebman et al. in press) provides precise estimates of the stellar and enshrouding cloud parameters, with an effective temperature of the star of 2000-2200 K, and a radius of the shell of R = 263 ± 10 AU. Lane et al. (2005) first used long-baseline near-IR (2.2 µm baselines smaller than 85 m) interferometry in November-December 2004 to provide the first direct measurement of the angular size of V838 Mon using the two-telescope recombiner Palomar Testbed Instrument (PTI). At this epoch, they measured an angular diameter for the central source of Θ = 1.83 ± 0.06 mas. Reasonably accurate distances were determined from the light echo and using other methods (Munari et al. 2005;Sparks et al. 2008;Kamiński et al. 2011). We adopt the distance D = 6.1 ± 0.6 kpc from Sparks et al. (2008). Using this distance, the interferometric diameter translates to a linear radius for the supergiant of 1200 ± 150 R ⊙ (5.6 ± 0.7 AU). Lane et al. (2005) could also fit the data assuming an elliptical structure with a major axis oriented at P.A. 15 • whose extensions are 3.57×0.07 mas (i.e. the minor axis is unresolved).

V838 Mon is embedded in a dense, large scale environment (Tylenda & Kamiński 2012;Kamiński et al. 2011;[START_REF] Exter | Asymmetrical Planetary Nebulae VI conference[END_REF], potentially contaminating observations of the central star. High angular resolution studies with a field-of-view (FOV) of a few arcsec are an asset in that domain to isolate the central regions from the extended dusty cloud. A great advantage of optical interferometry studies is to fully isolate the measurements from the extended environment.

This letter presents optical interferometry measurements obtained with the Very Large Telescope Interferometer (VLTI). The observations are presented in Sect. 2. In Sect. 3 we analyze the MIDI mid-IR measurements followed by the AMBER K band measurements by means of simple geometrical models. The results are then discussed in Sect. 4.

Observations

Mid-IR interferometric data were obtained with the twotelescope recombiner MIDI (Leinert et al. 2004) between October 2011 and February 2012 using the 8m Unit Telescopes (UTs) of the Very Large Telescope. These observations provided dispersed fringe visibilities between 8 and 13 µm in addition to classical spectrophotometric capabilities.

Near-infrared observations were obtained with AMBER (a three-telescopes combiner located at the VLTI: Petrov et al. 2007) in November 2012 and February 2013, using the low spectral resolution mode (R=35) with the UTs. Unfortunately, these data are unusable due to a bad calibrator (see details in appendix A) and therefore not presented here. In April 2013, observations were repeated with the 1.8m Auxiliary Telescopes (ATs) under photometric conditions, using the compact configuration (B ≤ 35m) and more suitable calibrators. The extended configuration (B ≤ 140m) was subsequently used in May 2014 to observe V838 Mon. The great advantage of both instruments is that the FOV of the interferometric measurements are within the point spread function of a single telescope, i.e. the same value of ∼0.3" at 2.2 µm for the 1.8m ATs with AMBER, and at 10 µm for the 8m UTs with MIDI. This isolates the central source from the extended emission of the cloud in which V838 Mon is embedded. The observation log is presented in Table 1 and the (u, v) plane coverage is plotted in Fig. 1.

We reduced the AMBER data using the standard data reduction software amdlib v3.0.8 (Tatulli et al. 2007;Chelli et al. 2009). Near-infrared JHK photometry was obtained on a regular basis from the 1.2m telescope at the Mt. Abu Observatory, India (Banerjee & Ashok 2012). These measurements helped to prepare the interferometric observations. On March 1, 2013, the J, H, K band magnitudes were 7.12, 5.91, 5.43, respectively, almost unchanged since November 2011. Visibilities are a proxy to the object's size and shape. The AMBER visibilities are presented in Fig. 2 at selected wavelengths for all baselines, and for the longest baseline in Fig. 3 together with the JHK Mt. Abu spectrum. Strong variations as a function of wavelength are seen, especially at the edge of the H and K bands where they decrease relative to the bands centers (1.7 µm and 2.2 µm). All baseline lengths show this curved shape. We are confident this is not a flux bias or instrumental effect, but comes from the source itself.

Closure phases can be qualitatively interpreted as a proxy to asymmetries: a non-zero (modulo π) closure phase is a sign of asymmetries in the object's shape. The 2014 closure phases are of good quality and compatible with zero with a mean value of -1. • 04 ± 1. • 12. The object is therefore likely centro-symmetric as seen by AMBER. The spectral-range of the AMBER data is 1.5-2.5 µm.

The MIDI data were processed using the MIA+EWS software1 (Chesneau 2007). The data were secured in High Sens mode implying that the photometry is obtained subsequently to fringes. The quality of the data ranges from high quality (error level ∼10%) to poor quality for some baselines (error level reaching ∼25%). The spectrophotometric calibration of the mid-IR flux was performed with the calibrator HD 52666 (M2III, IRAS F12=18.3 Jy with a 5% error). The ten MIDI spectra from two telescopes were merged leading to a spectrophotometry with an 8% error level, representing the mid-IR flux of the source between October and February 2011. The spectrum was compared with measurements at the same epoch from the WISE2 and AKARI3 satellites obtained with much larger apertures. The

Analysis

The MIDI dispersed visibilities were translated to the simplest possible geometrical ad-hoc model -a Gaussian brightness distribution, as described in Leinert et al. (2004). The resulting fits shown in Fig. 4 depict the general appearance of the mid-infrared structure. The semi-major and semi-minor axes of the drawn ellipses represent the half-width at half maximum (HWHM) of the Gaussian. The extent of that structure increases from a FWHM of 25 mas at 8 µm to 70 mas at 13 µm, with a high flattening ratio. Despite the heterogeneity of the quality of the MIDI visibilities, the flattening of that structure is certain, but we can only provide a loose range for flattening ratio values between 1.1 and 2.5 near 9 µm, and between 1.3 and 60 near 12 µm. The major axis is oriented at a P.A. close to -10 • with an uncertainty of more than ±30 • . The flattening increases with the wavelength, from a roughly round structure at 8 µm to a very flattened structure at 13 µm. This may come from the contribution of the star, stronger at 8 µm. The flux from WISE agrees fairly well with the MIDI flux. AKARI data is offset in time from the MIDI data by 5 to 6 years. Given that the source's SED has been evolving in the IR (Wisniewski et al. 2008) it is not unexpected or surprising that there is some difference in the value of the 9 micron flux measured by AKARI and MIDI. This implies that most of the mid-IR flux in 2011 originates from the vicinity of the central star (≤ 0.4").

The 2013 short-baseline AMBER visibilities are all very close to unity in the band centers at 1.7 µm and 2.2 µm and are consistent with a completely unresolved object. The longbaseline, 2014 visibilities are close to V 2 = 0.8 at the same wavelengths.

The visibilities as a function of wavelength of the H and K bands decrease relative to the bands centers (1.7 µm and 2.2 µm), as shown in Fig. 3. This is indicative of an object whose shape changes as a function of wavelength. At the same time, the visibilities at the edge of the bands (1.5 µm, 2.0 µm and 2.4 µm) show a "plateau-like" shape as a function of spatial frequencies (see Fig. 2), indicative of an additional resolved component. We therefore analysed first the K band center using a single-component model, and then added an additional component to analyze the full AMBER dataset. We used the fitOmatic software (Millour et al. 2009) to derive an angular diameter assuming a uniform disk model around 2.2µm. If we include all the 2013 (short baselines) and 2014 (long baselines) observations, we find a diameter of 1.17 ± 0.37 mas. If we keep only the 2014 data, we obtain a diameter of 1.15 ± 0.20 mas. We adopt the latter diameter of 1.15 ± 0.20 mas as the 2013 data (short baselines) do not resolve the central star. This measured diameter is significantly smaller than the one published by Lane et al. (2005).

The AMBER closure phases are all equal to zero within the error bars. We can therefore constrain the flux of an hypothetical companion star in addition to the central merger. To achieve that, we added a point-source to our model and used a simulated annealing algorithm to constrain its parameters using the observed visibilities and closure phases and starting with random initial parameters. Our conclusion is that its flux was always smaller than 3% of the total flux, i.e. there would be at minimum 4 magnitudes difference between primary and secondary companion. The AMBER observations do not exclude the presence of the companion, but put important constraints on its NIR contribution to the SED.

The decreasing visibilities at the edges of the H and K bands indicate the presence of extended emission with a varying spectrum (as shown in lower panel in Fig. 3, blue dotted line). We cannot constrain the exact size of this contribution due to the plateau-like shape of visibilities as a function of spatial frequen-cies, but we can provide a lower limit of its size of θ ≥ 20 mas. We tentatively associate this emission to the MIDI elongated dusty structure. The extended component would correspond to the H and K bands contribution of the MIDI elongated structure. It is noteworthy that adding this additional source of emission does not significantly change the estimated diameter for the central star.

To complement these measurements, we used a JHK spectrum of the star from Mt Abu (India), obtained on 11 December 2013, to derive a low-resolution spectrum of the central star and the extended emission (Fig. 3, lower panel). The stellar spectrum is typical of a cool luminous star (see Lanc ¸on & Wood 2000; Banerjee & Ashok 2002, for a comparison) with prominent first overtone CO bands at 2.29 µm and beyond. Similar, albeit weak, CO second overtone features are also seen in the H band. The deep down-turns at the band edges are due to the presence of water in its atmosphere (Banerjee et al. 2005). Other molecular features of VO and TiO are also seen, including some rather uncommon bands of AlO. The spectrum of the extended structure shows emission bands at the H and K band edges, where the water absorption occurs in the total spectrum.

Discussion

The P.A. of the major axis, as inferred from the dusty flattened structure discovered by MIDI is ≈ -10 • . The MIDI size increase is expected for a non-truncated dusty disk. At the adopted distance, the dust is distributed between 150 and 400 AU from V838 Mon. Such findings are also in agreement with the earliest accurate spectropolarimetric measurements reported by Wisniewski et al. (2003a,b). They reported an intrinsic polarisation with a P.A. of 127. • 0 ± 0. • 5 interpreted as scattering by a disk with a major axis at a P.A. of 37 • .

The 2004 PTI interferometric data had been interpreted as either a uniform disk or an elliptical model, both consistently explaining the data well. We note that the P.A.=15 +3 -27

• derived from PTI is roughly aligned with our longest AMBER baselines, i.e. we would be sensitive to the major axis size in case the correct Lane et al. (2005) model had been that of an elongated structure. Our AMBER data suggest a different picture. The supergiant has decreased in angular size by ≈ 40% in ≈ 10 years since the measurements of Lane et al. (2005). The linear radius is now estimated to 750 ± 200 R ⊙ (3.5 ± 1.0 AU). Such a radius is comparable to the famous M2 Iab supergiant Betelgeuse (Haubois et al. 2009, 885 ± 90 R ⊙ ). Such an evolution has already been suggested by [START_REF] Morisset | [END_REF]). This diameter decrease means the star's photosphere has shrinked during the intervening period at an approximate rate of 1km s -1 .

The good quality of the closure phases provides crucial constraints and excludes any secondary star brighter than a few percent of the total flux in the 6-600 AU separation range. Yet, we caution that a B3V star whose K band magnitude is fainter than 15 would be totally undetectable by the interferometer.

The AMBER interferometric observations obtained 10 yr after PTI suggest a contraction of the central star by about 40%, which is in line with the spectral change of V838 Mon toward a 'normal' M-type supergiant. No clear sign of deformation of the photosphere is detected.

The extended emission found in the AMBER data is enhanced relative to the central star around the CO overtone absorption features and the H 2 O absorption at the band edges. This means there is water emission around the star at an angular scale larger than 20 mas, and smaller than 250 mas (field of view of one telescope at 2.2 µm). We tentatively associate this emission to the MIDI elongated dusty structure.

Several hypotheses could be proposed to explain the detection of the flattened structure at large distance, as observed by MIDI and supported by AMBER.

The most likely hypothesis is that the flattened structure is simply the relic of the large dust formation event as a consequence of the merging of the two stars (Wisniewski et al. 2008(Wisniewski et al. , 2003a;;Nicholls et al. 2013). The ejecta velocity was low at outburst time -less than 200km s -1 as derived from P-Cygni profiles by Lynch et al. ( 2004) from several lines. In the 10 years since outburst, ejected material would not have travelled beyond 400 AU. This scenario is compatible with our finding.

The interferometric observations may also be discussed in the context of the B3V companion discovered around V838 Mon (Desidera & Munari 2002;Wagner & Starrfield 2002;Munari et al. 2002Munari et al. , 2005) ) and observed in the postoutburst spectra. The presence of the companion would explain the X-ray activity discovered in the vicinity of V838 Mon (Antonini et al. 2010). It is important to note that the distance of this companion from the central star was investigated by Munari et al. (2007), and they estimated the minimum separation to about 28 AU. On the other hand, Tylenda et al. (2009) derived a separation of ≈ 250 AU. This is well within the separation range detectable by AMBER. The MIDI measurements suggest a large separation of about 100 AU. In 2009, the B3V hot signature was no longer observed [START_REF] Kolka | The Astronomer's Telegram[END_REF]. We would propose then that the companion is building a circumbinary disk via the Lagrangian point L2, as proposed in some observational and theoretical studies on other sources (Plets et al. 1995;Meilland et al. 2010;Millour et al. 2011).

To conclude, we find a flattened dusty structure around V838 Mon. We also find that the central star decreased its diameter by nearly 40% in 10 years.

The interferometric picture of V838 Mon as of today is the following: it is now slowly becoming an anonymous red supergiant, surrounded by a flattened, probably transitory, dusty environment extending up to several hundreds of AU.

Fig. 1 .

 1 Fig. 1. Left: (u, v) coverage of the AMBER observations. Small dots are observations with ATs and larger circles with UTs. Right: (u, v) coverage of the MIDI observations.

  Fig.2. AMBER data shown as a function of spatial frequencies for selected wavelengths. The black line shows our best-fit model of a uniform disk plus an extended component (50 mas FWHM here).

Fig

  Fig. 3. Top: The 0.85-2.5 micron spectrum of V838 Mon on 11 December 2013 corrected for extinction using E(B-V) = 0.5. For comparison, we show the AMBER longest baseline visibility (thick black lines, right axis), and indicate with arrows the wavelengths of Fig. 2. Bottom: Spectra extracted from the AMBER data according to a downgraded resolution version of the above spectrum (black line). Red dashed line: the star spectrum. Blue dotted line: extended emission spectrum.

Fig. 4 .

 4 Fig. 4. Top: MIDI spectrum obtained in 2011. We added for comparison AKARI (triangle) and WISE (square) fluxes. Bottom: 2D Gaussian fit to the whole MIDI dataset. The dashed lines represent the direction of the observation baselines.

  The selected calibrator for the 2012 observations, HD 45299, is a previously-unseen visual binary star (see Fig.A.1) that we detect at 16.0 mas separation, position angle (P.A.) -67 • and 0.6

  Fig. A.1. AMBER closure phases and differential visibilities (no calibrated V 2 available) for HD45299 (color points with error bars) together with a best-fit model of a binary star (black lines). The spatial frequencies are projected along the binary direction (-67 • ) to show the binary modulation.

Table 1 .

 1 Log of AMBER and MIDI observations of V838 Mon.

	Date	Stations	Calibrators	Wavelength	Nb. obs.
	MIDI
	10/10/2011	U3-U4	1 •	8-13 µm	2
	13/12/2011	U2-U3	1 •	8-13 µm	1
	14/12/2011	U2-U3	1 •	8-13 µm	1
	11/02/2012	U3-U4	1 •	8-13 µm	1
	AMBER
	14/04/2013 A1-B2-D0	6 ⊳ , 7 ⊲	1.54-2.50 µm	3
	15/04/2013 A1-C1-D0	6 ⊳ , 8 △	1.56-2.50 µm	2
	06/05/2014	A1-G1-J3	9

▽ 1.53-2.46 µm 1 Notes. Calibrator angular diameters from SearchCal/JMMC (Bonneau et al. 2006) and getCal/NExScI. • : HD 52666 2.53 ± 0.21 mas, ⊳ : HD 52265 0.44 ± 0.03 mas, ⊲ : HD 64616 0.51 ± 0.04 mas, △ : HD 63660 0.96 ± 0.07 mas, ▽ : HD 54990 0.57 ± 0.59 mas.

http://home.strw.leidenuniv.nl/∼nevec/MIDI/

taken in April 2 and October 10, 2010

obtained between May 2006 and August 2007

Acknowledgements. O. Chesneau had wished to express his deepest gratitude to the hospital staff for their professionalism and devotion during the writing of this letter. This letter makes use of the CDS, the JMMC and NExScI services.


Emmanuel Soubies

Laure Blanc-Féraud

Sébastien Schaub

Gilles Aubert

SPARSE RECONSTRUCTION FROM MULTIPLE-ANGLE TOTAL INTERNAL REFLECTION FLUORESCENCE MICROSCOPY

Keywords: Inverse problems, Image reconstruction, Highresolution imaging, Evanescent wave microscopy

Super-resolution microscopy techniques allow to overstep the diffraction limit of conventional optics. Theses techniques are very promising since they give access to the visualisation of finer structures which is of fundamental importance in biology. In this paper we deal with Multiple-Angle Total Internal Reflection Microscopy (MA-TIRFM) which allows to reconstruct 3D sub-cellular structures of a single layer of ∼ 300 nm behind the glass coverslip with a high axial resolution. The 3D volume reconstruction from a set of 2D measurements is an ill-posed inverse problem and a regularization is essential. Our aim in this work is to propose a new reconstruction method for sparse structures robust to Poisson noise and background fluorescence. The sparse property of the solution can be seen as a regularization using the 'ℓ 0 norm'. In order to solve this combinatorial problem, we propose a new algorithm based on smoothed 'ℓ 0 norm' allowing minimizing a non convex energy, composed of the Kullback-Leibler divergence data term and the ℓ 0 regularization term, in a Graduated Non Convexity framework.

INTRODUCTION

Understanding tumor organisation and development is a major challenge in current biological researches. Super-resolution microscopy techniques are widely used in this context since they allow the visualisation of small sub-cellular structures near the cell membrane. Indeed, these new and powerful tools allow to overstep the limitations of conventional optics for which the axial resolution is restricted to a few hundred of nanometers (e.g. ∼ 500 nm for confocal microscope).

In this study we deal with Total Internal Reflection Fluorescence Microscopy (TIRFM) [START_REF] Axelrod | Total internal reflection fluorescence[END_REF] which gives access (in the better configuration) to a single layer of ∼ 100 nm behind the glass coverslip (5× better, in terms of axial resolution, than confocal microscopy) and provides low background fluorescence and high signal-to-noise ratio images. While single TIRFM acquisition is actually a unique projection (2D image) of the observed volume and can not provide 3D information, Multiple-Angle TIRFM (MA-TIRFM) provides a set of TIRFM acquisitions at different incident angles. Through the development of dedicated algorithms, these acquisitions allow to reconstruct the 3D volume of the sample and overcome the axial resolution of single TIRFM. Indeed the set of 2D images obtained from Contact: soubies@i3s.unice.fr, blancf@i3s.unice.fr, sebastien.schaub@unice.fr, gilles.aubert@unice.fr MA-TIRFM corresponds to different projections (with different integration depth) of the 3D volume on the (x, y) plane. In the following we study the ill-posed inverse problem of finding axial positions of the observed structures from this set of images. This is a reconstruction problem.

Some authors have proposed different methods to estimate the 3D sub-cellular volume from MA-TIRFM acquisitions. All of these methods use a strong shape prior on the unknown structures. In [START_REF] Ga Truskey | Total internal reflection fluorescence microscopy (tirfm). II. topographical mapping of relative cell/substratum separation distances[END_REF][START_REF] Bence P Ölveczky | Mapping fluorophore distributions in three dimensions by quantitative multiple angle-total internal reflection fluorescence microscopy[END_REF] axial cell membrane position is estimated for each pixel of the (x, y) plane using a simple one-dimensional geometry (top hat). Estimating position and diameter of isolated secretory granules have been studied in [START_REF] Rohrbach | Observing secretory granules with a multiangle evanescent wave microscope[END_REF][START_REF] Loerke | Quantifying axial secretory-granule motion with variable-angle evanescent-field excitation[END_REF] where granules are modelled by spherical or cubic volumes. A statistical framework for curvilinear structures reconstruction (microtubules) is proposed in [START_REF] Yang | 3D reconstruction of microtubules from multipleangle total internal reflection fluorescence microscopy using Bayesian framework[END_REF]. More recently, 3D particle estimation is achieved in a MAP-Bayesian framework [START_REF] Liang | A Bayesian method for 3D estimation of subcellular particle features in multi-angle tirf microscopy[END_REF] and a Marked Point Process approach is used in [START_REF] Soubies | A 3D model with shape prior information for biological structures reconstruction using Multiple-Angle Total Internal Reflection Fluorescence Microscopy[END_REF] to deal with objects features estimation.

In this paper, we propose a new model without strong shape constraints on the observed structures. We formulate the reconstruction problem as a minimization of the sum of two terms ; a data fidelity term defined from the likelihood function of the observations given the solution, and, a regularization term. We focus particularly on finding sparse solutions to the problem since we are interested in reconstructing cell membranes and fine particles which are sparse structures in the medium. Moreover, sparse property could also be useful in the case where the solution is searched in a particular basis (e.g. in a wavelet basis). Sparsity of the solution can be modelled using the 'ℓ 0 norm'1 as regularization term. Following recent work [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed l0 norm[END_REF], this difficult combinatorial optimization problem is tackled by using the SL0 algorithm, based on smoothed 'ℓ 0 norm' and using a Graduated Non Convexity (GNC) approach [START_REF] Blake | Visual reconstruction[END_REF]. In [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed l0 norm[END_REF], the SL0 algorithm is developed for finding sparse solutions of underdetermined system of linear equations. In the following we propose a new algorithm based on this principle allowing the minimization of a non convex energy defined by the sum of the Kullback-Leibler divergence data term adapted to Poisson noise and the ℓ 0 regularization term enforcing sparsity of the solution.

The paper is organized as follows. Section 2 formalizes the TIRF image formation. We begin Section 3 with the formulation of the inverse problem as a minimization problem and we propose to solve this problem by minimizing a sequence of functional where the 'ℓ 0 norm' is approximated using smooth functions. We show that we can define a sequence of functional which starts with a convex one (on a large convex set) and introduce progressively the non convexity of the 'ℓ 0 norm'. Then we give a theoretical result on the limit of the minimizers with respect to the 'ℓ 0 norm'. Finally, some numerical results are presented in Section 4.

IMAGE FORMATION AND PROBLEM STATEMENT

Total internal reflection of a light beam, at the interface between two mediums of refractive indices ni (incident) and nt (transmitted), occurs when the incident angle α of the light beam becomes greater than a critical angle αc. This phenomenon produces an evanescent wave capable of exciting fluorophores near the dielectric surface. Theoretically [START_REF] Axelrod | Total internal reflection fluorescence[END_REF], the axial profile of the evanescent field can be modelled using a decreasing exponential depending on the incident angle. In the following, we denote by:

-

O d = {1, . . . , n} × {1, . . . , m} × {1, . . . , p} ⊂ N 3 the image space (3D) -M = {1, . . . , n} × {1, . . . , m} × {1, . . . , L} ⊂ N 3 the measures space (2D + angle)
where n and m are the images' size (in the (x, y) plane) indexed in the following by i ∈ {1, . . . , n} × {1, . . . , m}, p is the number of voxels of the 3D unknown image in the axial direction (indexed by j) and L the number of incident angles used (i.e number of acquired 2D images), indexed by l.

Let f : O d -→ R denote the 3D unknown fluorophore density and s : M -→ R the set of recorded images. According to the model given in [START_REF] Axelrod | Total internal reflection fluorescence[END_REF] and considering shot noise and background fluorescence degradations, TIRF image formation can be formulated as follows:

∀(i, l) ∈ M , s(i, l) = P (H(f + bg)(i, l)) (P)
where P denotes the Poisson distribution, which mimics the photon collection process, bg is a constant that models the background fluorescence produced by parasites light sources in the medium (e.g. autofluorescence, light reflection...) and H is the TIRF operator defined by: ∀f

∈ ℓ 2 (O d ), ∀(i, l) ∈ M , H(f)(i, l) = p j=1 f(i, j) exp - zj d l δz(j)µi . (1)
Here, δz(j) is the j th discretization step along the (Oz) axis, µi is the i th pixel area (2D) and d l models the penetration depth of the evanescent wave corresponding to the incident angle α l , l ∈ {1 • • • L}.

SPARSE RECONSTRUCTION METHOD

Optimization Problem Formulation

Finding a solution of the ill-posed inverse problem (P) (i.e determine f in (P) from acquisition s) can be formulated as a minimization problem. Assuming that observations s given the fluorophore density f are pairwise independent (which can be justified by the fact that for a given f, data depend only to the noise which is pairwise independent), then maximizing the likelihood function of the observations s given f leads to the following problem:

f = arg min f∈ℓ 2 (O d ) J d (f) (P1)
where J d is defined by2 :

J d (f) = (i,l)∈M -s(i, l) log [H(f + bg)(i, l)] + H(f + bg)(i, l) . (2)
Let's give some properties of J d which will be used, in the following, to prove theoretical results:

-It is easy to show that J d is continuous and convex -It has been proved in [START_REF] Saima | Space Variant Blind Image Restoration[END_REF], in the case where H is a convolution operator, that:

∃a > 0, b ∈ R s.t ∀f ∈ ℓ 2 (O d ), J d (f) ≥ a f 1 + b (3)
where . 1 is the ℓ 1 norm. TIRF operator (1) satisfies the conditions used in [START_REF] Saima | Space Variant Blind Image Restoration[END_REF] and then (3) holds which implies that J d is coercive.

Since we are interesting in finding a sparse solution of the problem, we will consider, in the following, the regularized problem:

f = arg min f∈ℓ 2 (O d ) G(f) := J d (f) + λ f 0 (P2)
where f 0 = # (i, j) ∈ O d ; f(i, j) = 0 and λ ∈ R + is a weight parameter.

SL0 Algorithm

Solving problem (P2) with a combinatorial search is NP-hard. However, it has been shown [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed l0 norm[END_REF] that a smoothed 'ℓ 0 norm' can be used for finding sparse solution of an undetermined system of linear equations. The main idea in [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed l0 norm[END_REF] is to approximate the 'ℓ 0 norm' by a suitable continuous function depending on a parameter σ > 0. For example, they propose the following smooth approximation:

f 0 ≃ Jσ(f) := N - (i,j)∈O d hσ(f(i, j)) (4)
where N = |O d | is the number of voxels in f and hσ a continuous function defined by3 :

hσ(s) = exp - s 2 2σ 2 . ( 5
)
It is easy to show that Jσ tends to the 'ℓ 0 norm' when σ tends to zero and consequently the smaller the values of σ the better the approximation. On the other side the larger the value of σ, the smoother the approximation. Following the idea in [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed l0 norm[END_REF], we propose to solve a series of optimization problems:

f σ = arg min f∈ℓ 2 (O d ) Gσ(f) := J d (f) + λJσ(f) (P σ 2 )
for a decreasing sequence of σ. At each step (i.e for each σ of the sequence), problem (P σ 2 ) is tackled using a descent algorithm 4 initialized with the solution obtained for the previous (larger) value of σ. These iterates allow to introduce gradually the non-convexity of the problem and escape from local minima (GNC approach [START_REF] Blake | Visual reconstruction[END_REF]). Let remark that our work differs from [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed l0 norm[END_REF] in the sense that [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed l0 norm[END_REF] focus on the resolution of linear systems under 'ℓ 0 norm' sparsity condition.

We will show that a large enough value of σ ensures that problem (P σ

2 ) has no local minima. It can be easily shown that Jσ is strictly convex on [-σ, σ] N . Then the following property states: Property 1. There exists a large enough value of σ denoted σ0 such that:

∀σ ≥ σ0 ∀(i, j) ∈ O d , |f σ (i, j)| ≤ σ (6)
where f σ is a solution of (P σ 2 ).

Proof. Let f be a solution of the problem (P1) (i.e a solution of (P σ 2 ) with λ = 0). Such a solution exists since J d is continuous, convex and coercive. Because f σ is a solution of (P σ 2 ) and that ∀f ∈ ℓ 2 (O d ), Jσ(f) ≤ N we have:

J d (f σ ) + λJσ(f σ ) ≤ J d ( f) + λJσ( f) ≤ J d ( f) + λN . (7)
Using result (3) and the fact that ∀f ∈ ℓ 2 (O d ), Jσ(f) ≥ 0 we have:

∃a > 0, b ∈ R such that, a f σ 1 + b ≤ J d (f σ ) ≤ J d (f σ ) + λJσ(f σ ) ≤ J d ( f) + λN . (8)
Consequently, ∀(i, j) ∈ O d , |f σ (i, j)| is bounded by a value independent of σ and then by taking σ0 equal to this bound, equation ( 6) is satisfied.

From Property 1, we ensure that by taking a large enough value of σ1 to initialize the algorithm, the solution f σ 1 of the approximated problem will be in the convex part of Jσ (i.e f σ 1 ∈ [-σ1, σ1] N ) which ensures f σ 1 to be unique since Jσ is strictly convex on [-σ, σ] N . Finally the proposed algorithm is described in Algo 1.

Input: f0, σ = [σ1, • • • , σK ], λ for k ← 1 to K do σ ←-σ k ; f σ 0 ←-fk-1 ;
Compute a finite number n of iterates of a descent algorithm on problem

(P σ 2 ) -→ f σ n ; fk ←-f σ n return fK
Algorithm 1: Iterative algorithm for solving (P2)

The decrease of the sequence of σ is chosen such as σ k = cσ k-1 , k ≥ 2 where c is a coefficient between 0.5 and 1. The first value, σ1 is chosen according to Property 1 in order, as explained before, to ensure the uniqueness of the minimizer of (P σ 2 ). We stop the algorithm when σ becomes smaller than the estimated bg (the bg 4 For example the regularized Richardson-Lucy algorithm [START_REF] Dey | Richardson-Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution[END_REF]. constant is estimated from an acquired image region, selected manually, without biological structures). Note that for each σ in Algorithm 1, only a finite number n of a descent algorithm iterates are performed which is sufficient. In practice it is not necessary to converge to the minimiser of problem (P σ 2 ) for each σ. In this paper we use the regularized Richardson-Lucy algorithm [START_REF] Dey | Richardson-Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution[END_REF] to perform the descent step. Finally, we initialize f0 by a uniform 3D image equal to the estimated value of bg.

A Theoretical Result on the Limit of (P σ

2 ) The following theorem ensures that the global minimizers of (P σ 2 ) are also global minimizers of (P2) when σ tends to zero.

Theorem 1. Let S0 be the set of global solutions of problem (P2) and Sσ the one of problem (P σ 2 ), then:

lim σ→0 Sσ ⊆ S0 (9)
in the sense of Painleve-Kuratowski convergence [START_REF] Kuratowski | Introduction to set theory and topology[END_REF].

Proof. It is easy to prove the existence of a solution of (P2) and (P σ 2 ) which implies that S0 and Sσ are a non-empty sets. We define by Eγ(J d ) := {f s.t J d (f) ≤ γ} the γ sublevel set of J d . Let f ∈ S0 and f σ ∈ Sσ , then we have the following inequality:

∀σ, Gσ(f σ ) ≤ Gσ( f) ≤ G( f) . (10)
The first inequality comes from the fact that f σ is a minimizer of (P σ 2 ) and the second inequality is true since ∀σ, ∀f ∈ ℓ 2 (O d ), Jσ(f) ≤ f 0. From (10) we get:

∀σ, f σ ∈ Sσ ⊆ E G( f) (Gσ) . (11)
Now, we show that:

E G( f) (Gσ) -→ σ→0 E G( f) (G) = S0 (12)
which is equivalent (in the sense of Painleve-Kuratowski convergence [START_REF] Kuratowski | Introduction to set theory and topology[END_REF]) to satisfy the two following conditions:

∀f ∈E G( f) (G), ∃σ0 > 0, ∀σ > 0, σ ≤ σ0 ⇒ f ∈E G( f) (Gσ) , (13) ∀f
/ ∈E G( f) (G), ∃σ0 > 0, ∀σ > 0, σ ≤ σ0 ⇒ f / ∈E G( f) (Gσ) . (14)
Using the fact that ∀f ∈ ℓ 13) is directly deduced.

2 (O d ), Jσ(f) ≤ f 0 =⇒ Gσ(f) ≤ G(f) it comes that ∀σ, E G( f) (G) ⊆ E G( f) (Gσ) from which (
In order to show condition [START_REF] Hadley | Bayesian-based iterative method of image restoration[END_REF]

, let f / ∈ E G( f) (G) which is equiv- alent to: G(f) > G( f) . (15)
On which condition on σ do we have

f / ∈ E G( f) (Gσ) ? f / ∈ E G( f) (Gσ) ⇔ Gσ(f) > G( f) ⇔ J d (f) + λJσ(f) > G( f) ⇔ J d (f) + λN -λ (i,j)∈O d hσ(f(i, j)) > G( f) ⇔ λ (i,j)∈O d exp - f 2 ij 2σ 2 < J d (f) + λN -G( f) . (16)
Let N 0 f be the number of zeros of f.

If N 0 f = N , (16) is equivalent to 0 < J d (f) -G( f) = G(f) -G( f)
which is satisfied ∀σ since we have inequality (15). If N 0 f < N , we define fm = min{|fij| = 0} (i,j)∈O d , then we can write:

λ (i,j)∈O d exp - f 2 ij 2σ 2 ≤ λN 0 f + λ(N -N 0 f ) exp - f 2 m 2σ 2 . ( 17
)
Thus, ∀σ for which:

λN 0 f + λ(N -N 0 f ) exp - f 2 m 2σ 2 < J d (f) + λN -G( f) (18)
is satisfied, (16) will be also satisfied thanks to (17). Let's now working on inequality (18). We can rewrite (18) as follows:

λ(N -N 0 f ) exp - f 2 m 2σ 2 < J d (f) + λ(N -N 0 f ) -G( f) . ( 19
)
Using inequality (15) we have J d (f) + λ(N -N 0 f ) = G(f) > G( f) which implies that the right hand side of ( 19) is strictly positive. Thus we can make the following development:

λ(N -N 0 f ) exp - f 2 m 2σ 2 < G(f) -G( f) ⇔ - f 2 m 2σ 2 < log 1 λ(N -N 0 f ) G(f) -G( f) ⇔ 1 σ 2 > - 2 f 2 m log (A) (20)
where A 2 = σ0 Consequently, there exists a σ0 > 0 such that inequality (18) holds ∀σ ≤ σ0 and, combined with (17), imply (16) which is equivalent to f / ∈ E G( f) (Gσ) and prove [START_REF] Hadley | Bayesian-based iterative method of image restoration[END_REF]. Both conditions ( 13) and ( 14) are then satisfied that is equivalent to (12) which, with [START_REF] Saima | Space Variant Blind Image Restoration[END_REF] finish the proof.

A = 1 λ(N -N 0 f ) G(f) -G( f) .

NUMERICAL RESULTS

In this section we present some numerical results on synthetic data. We use a simulated membrane (3D) shown on Figure 2. From this synthetic sample, a numerical simulator of the microscope is used in order to get the MA-TIRFM acquisitions (set of 2D images) according to the model in (P). Two images of this set are also presented on Figure 2. These acquisitions are performed using a critical angle equal to α = 61 • and 30 different incident angles from 62 • to 80 • . The noise level bg (on the acquired images) is set to 10% of the maximum intensity of the simulated sample.

Since we are interested in the estimation of the axial membrane position, we quantify the robustness of the method using the position errors along the axial direction as shown on Figure 3. For each pixel in the (x, y) plane we compute the absolute error between the axial membrane positions of the simulated sample and the reconstructed one. This error is quantified using the colormap of the Figure 3 (dark blue = 0 nm and dark red corresponds to an error > 50 nm). As we can see on Figure 3 right, we obtain a high accuracy reconstruction with the Algorithm 1 compared to the solution obtained without regularization (Figure 3 left). Indeed the axial position errors are between 5 and 15 nm for the membrane reconstructed using the ℓ 0 regularization whereas it varies from 30 to more than 50 nm for the reconstruction without regularization. Following the theoretical TIRFM model (P), it is easy to see that the exponential decay of the evanescent field causes that the deeper the structures, the lower the signal contribution on the acquired images (ill-posedness of the problem) which explain the large errors obtained without regularization (Figure 3 left). We have also tested to reconstruct the membrane using Total Variation (TV) regularization. The reconstruction obtained is similar, in term of axial position errors, to the one presented on the Figure 3 left (without regularization). Finally, using an ℓ 1 regularization, the axial position errors are of the same order of magnitude than the errors showed on Figure 3 right (ℓ 0 regularization) for structures below 200 nm in the axial direction but increase drastically above this limit while the errors with the ℓ 0 regularization remain between 5 and 15 nm until structures greater than 300 nm. However, due to page limit comparisons with TV and ℓ 1 regularizations are not presented here as well as results on the reconstruction of fine particles.

CONCLUSION

In this paper we proposed a new model for the MA-TIRFM reconstruction problem based on a ℓ 0 regularization to enforce the sparsity of the solution. Following the idea proposed in [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed l0 norm[END_REF] for finding sparse solutions of underdetermined system of linear equations, we propose to minimize a sequence of functional Gσ (Problem (P σ 2 )) which are smooth approximations of the objective function G (Problem (P2)) and for which, global minimizers when σ tends to zero are also global minimizers of G (Theorem 1). Thus if for each σ of the sequence, Algorithm 1 does not get trapped into a local minimizer of Gσ, it will converge to a global solution of P2. Finally, simulations results show the accuracy of this approach. Future work will focus on the possibility to construct an optimal sequence σ which ensures the convergence to a global solution of P2.

Figure 1

 1 Figure 1 contains an example of MA-TIRFM acquisitions from a sample of cortactin. It is readily seen from these images that when the penetration depth decreases (from left to right), some fluorophores are not excited anymore. These structures are the deeper ones and so we can get 3D information from this set of images.

Fig. 1 .

 1 Fig. 1. Multiple-angle TIRFM (MA-TIRFM) acquisitions (2D images s(i, l), (i, l) ∈ M ) from a cortactin sample (protein involved in cell adhesion). Left : deeper penetration depth (α ≃ αc). Right : thinner penetration depth (α ≃ αmax).

  We can discern two cases: • If A ≥ 1 then the right hand side of (20) is negative or zero and so, (20) holds ∀σ. • If A < 1 then the right hand side of (20) is strictly positive and (20) holds only if σ < fm log -1 2 1

Fig. 2 .Fig. 3 .

 23 Fig. 2. Left: Synthetic membrane. Right: Microscope acquisition (numerical simulations) with the angles α = 62 • and α = 65.1 • .

We write 'ℓ 0 norm' because the absolute scalability is not satisfied.

Note that minimizing J d in (2) is equivalent to minimize the well known Kullback-Leibler divergence.

In[START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed l0 norm[END_REF] they propose a general class of function for hσ. In this paper we will only deal with the function hσ defined in[START_REF] Loerke | Quantifying axial secretory-granule motion with variable-angle evanescent-field excitation[END_REF].