Jérôme Fran- Çois Jean

Keywords: lemand, roumain, espagnol, arabe

à nos esprits hermétiques d'informaticiens.

D�������� ��� ���������������� �� ��������� ������ �� ����� ��������� vii Résumé en français N ���� ������� de thèse porte sur la détection et la modélisation de configura- tions dynamiques dans des séquences d'images en recourant à des approches statistiques locales sans apprentissage supervisé. Deux cas peuvent se présenter :

1. les objets étudiés n'interagissent pas, et les dynamiques individuelles peuvent être analysées indépendamment ; 2. les objets étudiés interagissent, et la dynamique à analyser est celle du groupe d'objets entier.

En ce qui concerne les dynamiques individuelles, nous nous intéressons à certaines structures de la cellule, les vésicules, qui sont observées à la frontière de la cellule -où elles jouent un rôle essentiel pour la vie et les échanges cellulaires -en microscopie de fluorescence par réflexion totale interne (TIRF).

La dynamique de groupe est notamment rencontrée dans les mouvements de tissus cellulaires, le développement embryonnaire ou les mouvements de foules.

Ces derniers sont particulièrement intéressants en termes d'évaluation car de nombreuses méthodes d'analyse de mouvement de foule ont été proposées dans la littérature. Ce domaine d'application représente donc un excellent support expérimental.

Dans les deux cas d'étude, nous nous intéressons aux mouvements de collections d'objets similaires -molécules, cellules ou piétons -évoluant sur un fond immobile, et imagés par un dispositif statique. Du point de vue du traitement d'image, de nombreuses questions sont abordées, notamment la détection (spatio)temporelle d'objets ou d'évènements, la caractérisation de dynamiques ou encore l'estimation de grandeurs physiques. Dans l'un ou l'autre cas d'étude, nous abordons les problématiques selon une démarche commune, essentiellement dirigée par les données et mettant en oeuvre des tests statistiques. Par ailleurs, nous avons le souci de proposer des méthodes nécessitant le réglage d'un faible nombre de paramètres, soit peu sensibles, soit calibrés avec des règles statistiques.

Nous nous intéressons en premier lieu à l'analyse de dynamiques membranaires de la cellule en microscopie optique de fluorescence. Pour cause, la vie de la cellule est sujette à un équilibre très précis des différents éléments chimiques qui la peuplent.

La régulation de cet équilibre passe nécessairement par des échanges avec le milieu extérieur, donc par le franchissement de la membrane plasmique qui enceint la cel-viii lule. Notre objectif est de caractériser les dynamiques de vésicules évoluant dans la cellule. Les vésicules sont des intermédiaires de transport véhiculant divers composants dans travers la cellule. Dans le cas de l'exocytose, les vésicules déplacent ces composants à l'extérieur de la cellule. Pour ce faire, les composants sont transportés jusqu'à la membrane plasmique de la cellule, puis la vésicule fusionne avec celle-ci pour libérer son contenu.

Pour analyser ces dynamiques membranaires, un premier travail consiste à détecter les vésicules d'intérêt. Pour cela, nous proposons une méthode de détection de spots performante et simple à paramétrer. Une analyse, requérant des informations plus biologiques et biophysiques, permet ensuite de détecter et reconnaître certains évènements dynamiques dans la cellule, en particulier la fusion de la vésicule avec la membrane plasmique. Les paramètres biophysiques associés sont ensuite estimés.

Dans ce cadre, nous proposons des modèles de fusion, dont nous estimons les paramètres pour plusieurs conditions expérimentales. En particulier, nous mettons en évidence les comportements différents de deux protéines d'intérêt : le récepteur à la transferrine (TfR) et la Langérine.

Dans la troisième partie de la thèse, nous étudions la classification de mouvements de groupes, la détection des chemins les plus empruntés dans la vidéo et la détection de comportements anormaux rares, parfois suspects. Comme indiqué plus haut, nous avons principalement évalué nos méthodes d'analyse de mouvements de groupe sur des vidéos de foules. Aucune comparaison de méthode n'est actuellement possible sur des images de biologie. En revanche, nos approches trouvent un intérêt certain pour analyser des comportements collectifs souvent observés en imagerie biologique.

La thèse comprend trois parties : la détection de spots, l'analyse de dynamiques individuelles et l'analyse de dynamiques de groupes.

Partie I. Détection de spots

Dans le cadre de l'analyse d'images de microscopie, la détection automatisée des éléments à étudier constitue la plupart du temps une étape préliminaire essentielle qui conditionne l'ensemble des analyses ultérieures, que ce soit le suivi de vésicules [START_REF] Meijering | Methods for cell and particle tracking[END_REF], la classification de mouvements [START_REF] Sage | Automatic tracking of individual fluorescence particles: Application to the study of chromosome dynamics[END_REF] ou la caractérisation de dynamiques membranaires [START_REF] Mele | Automatic identification of fusion events in TIRF microscopy image sequences[END_REF]. Dans ce dernier cas, nous cherchons à segmenter les vésicules imagées en microscopie de fluorescence à réflexion totale interne (TIRF). Ces vésicules sont généralement de tailles semblables, et il convient donc d'estimer l'échelle caractéristique de ces objets.

Des études comparatives [START_REF] Smal | Quantitative comparison of spot detection methods in fluorescence microscopy[END_REF], Rezatofighi et al., 2012] ont permis d'évaluer les performances de nombreuses méthodes de détection non supervi-P����� I. D�������� �� ����� ix sées : produit multi-échelle d'ondelettes (WMP [Olivo-Marin, 2002]), détection multiéchelle après stabilisation de variance (MS-VST [START_REF] Zhang | Multiscale variance-stabilizing transform for mixed-Poisson-Gaussian processes and its applications in bioimaging[END_REF]), détection de « dômes » (HD [START_REF] Smal | A new detection scheme for multiple object tracking in fluorescence microscopy by joint probabilistic data association filtering[END_REF], MPHD [START_REF] Rezatofighi | A new approach for spot detection in total internal reflection fluorescence microscopy[END_REF]), filtre de rehaussement de contraste de spots (SEF [START_REF] Sage | Automatic tracking of individual fluorescence particles: Application to the study of chromosome dynamics[END_REF], TH [START_REF] Bright | Two-dimensional top hat filter for extracting spots and spheres from digital images[END_REF], MTH [Soille, 2003]), extraction de points caractéristiques (IFD [ter Haar Romeny, 2003])... La plupart de ces méthodes nécessitent un ajustement fin de plusieurs paramètres pour obtenir de bons résultats [START_REF] Rezatofighi | A new approach for spot detection in total internal reflection fluorescence microscopy[END_REF],Smal et al., 2010]. Plus précisément, ces détecteurs ont en commun deux types de paramètres : l'échelle des objets à détecter et le seuil de détection. Afin de réduire le nombre et la sensibilité des paramètres, nous développons une méthode de segmentation de vésicules avec sélection d'échelle automatique et seuillage adaptatif, dénommée ATLAS 1 . Elle s'appuie sur un filtre laplacien de gaussienne (LoG), dont la variance est automatiquement sélectionnée dans un ensemble fini de valeurs prédéfinies. Un seuillage adaptatif local permet alors de décider si un pixel appartient ou non à une vésicule, même lorsque l'arrière-plan est très complexe. Sous l'hypothèse d'une distribution localement gaussienne des intensités de l'image, le seuil est inféré, en chaque pixel, d'une probabilité de fausse alarme choisie par l'utilisateur pour l'ensemble des images traitées. deberg, 1998]. Cependant, pour des signaux discrets comme les images, il est recommandé d'utiliser des noyaux également discrets analogues au noyau gaussien, appelés plus simplement noyaux gaussiens discrets. On les obtient à l'aide des fonctions de Bessel B n [Lindeberg, 1990] :

Sélection automatique de l'échelle

∀ s ∈ S , ∀ n ∈ Z , G s (n) = e -s B n (s ) .
On applique ensuite un opérateur laplacien normalisé noté α(s ) ∆ qui permet d'augmenter le contraste des vésicules par rapport à l'arrière-plan [START_REF] Sage | Automatic tracking of individual fluorescence particles: Application to the study of chromosome dynamics[END_REF].

Le scalaire α(s ) désigne un facteur de normalisation de G s . Par associativité et commutativité de l'opération de convolution, le laplacien peut être appliqué une seule fois 1 Adaptive threshold of LoG with auto-selected scale x directement sur I , et non pas à chaque échelle. On obtient ainsi le LoG multiéchelle H :

∀ s ∈ S , H s = α(s ) G s * ∆I .
Nous cherchons ensuite, parmi les échelles de S, celle qui permet d'augmenter au mieux le contraste des vésicules. Pour la déterminer, on s'intéresse aux minima locaux de H dans l'espace S × Ω I , appelés blobs. Ces blobs sont principalement localisés en deux lieux précis :

• au centre des spots gaussiens correspondant aux vésicules ;

• au niveau des pixels brillants induis par le bruit.

Nous cherchons l'échelle pour laquelle la plus grande proportion de blobs est issue des vésicules, et non du bruit. La distribution du nombre de blobs par unité de surface en fonction de l'échelle nous informe sur l'échelle optimale. En effet, on peut comparer cette distribution calculée sur l'image I à la distribution obtenue sur une image g ne présentant que du bruit. Les courbes obtenues pour I et g se distinguent autour d'une certaine échelle s � , caractéristique de l'image I , qui offrira les meilleurs résultats de détection. Pour la détecter, nous proposons de maximiser le rapport entre les deux distributions. Ceci revient à faire un compromis entre taux de bonnes détections (spots) et taux de fausses détections (bruit).

Segmentation adaptative des spots

� p � = � W * H s � � � p � , σ 2 � p � = � W * H 2 s � � � p � -µ 2 � p � .
En supposant que H s � a une distribution locale gaussienne, on peut alors déduire le seuil local τ � p � :

τ � p � = Φ -1 � P FA � σ � p � + µ � p � ,
où Φ est la fonction de répartition de la loi normale centrée réduite et P FA est une probabilité de fausse alarme choisie par l'utilisateur.

P����� II. A������ ��� ���������� ������������� xi

Résultats expérimentaux

Nous utilisons différents bancs de test introduits dans [START_REF] Smal | Quantitative comparison of spot detection methods in fluorescence microscopy[END_REF], Pécot et al., 2015,Ruusuvuori et al., 2012], et nous en proposons quatrième2 , pour comparer notre méthode aux détecteurs de l'état de l'art : WMP [Olivo-Marin, 2002], MS-VST [START_REF] Zhang | Multiscale variance-stabilizing transform for mixed-Poisson-Gaussian processes and its applications in bioimaging[END_REF], TH [START_REF] Bright | Two-dimensional top hat filter for extracting spots and spheres from digital images[END_REF], MTH [Soille, 2003], HD [START_REF] Smal | A new detection scheme for multiple object tracking in fluorescence microscopy by joint probabilistic data association filtering[END_REF],

MPHD [START_REF] Rezatofighi | A new approach for spot detection in total internal reflection fluorescence microscopy[END_REF], SEF [START_REF] Sage | Automatic tracking of individual fluorescence particles: Application to the study of chromosome dynamics[END_REF], IDF [ter Haar Romeny, 2003], C-CRAFT [START_REF] Pécot | Background fluorescence estimation and vesicle segmentation in live cell imaging with conditional random fields[END_REF] et LR-MRF [START_REF] Ruusuvuori | Image segmentation using sparse logistic regression with spatial prior[END_REF]. Dans la quasi-totalité des cas, nous obtenons un meilleur score que la méthode la plus performante.

Qualitativement, nous mettons en évidence que la précision de la segmentation des spots sur images réelles, en particulier en comparant nos cartes de détection à celles de la méthode obtenant les meilleurs résultats sur simulations après ATLAS, à savoir MS-VST.

Le logiciel ATLAS peut être testé en ligne sur le portail web de l'équipe-projet Serpico3 .

Partie II. Analyse des dynamiques individuelles

Cette partie de la thèse est menée en collaboration avec l'équipe de Jean Salamero à l'Institut Curie, UMR 144 CNRS.

Nous étudions des complexes moléculaires impliqués dans les dernières étapes de l'exocytose, à savoir la fin de la phase de transport par les vésicules (le long des microtubules, puis dans le réseau de filaments d'actine), l'arrimage à la membrane et la fusion membranaire. Nous nous intéressons à trois protéines associées aux vésicules :

Rab11, associée aux vésicules durant la phase de transport, et Langérine et TfR, deux protéines transmembranaires transportées par des vésicules.

Pour imager la cellule vivante, on utilise généralement en microscopie optique des techniques de microscopie de fluorescence. Un faisceau laser est émis en direction de la cellule pour exciter des marqueurs fluorescents qui ont été liés aux protéines d'intérêt (pHluorine pour Langérine et TfR, mCherry pour Rab11). Dans notre étude, nous utilisons la microscopie TIRF, qui permet de n'exciter qu'une épaisseur très faible de la cellule au niveau de la membrane, et ainsi d'isoler la région où ont lieu les dernières étapes de l'exocytose [Axelrod, 2008]. En pratique, le faisceau laser est émis avec un angle d'incidence supérieur à l'angle limite de réfraction, de sorte qu'il est totalement réfléchi par l'interface entre la lamelle de verre et la cellule. Cependant, une faible quantité d'énergie (l'onde évanescente) pénètre dans la cellule, et son amplitude décroît exponentiellement avec la profondeur. Ainsi, seules les protéines marquées proches de l'interface sont principalement excitées.

xii

Classification des dynamiques membranaires de Rab11

Notre objectif est de caractériser le comportement dynamique des vésicules et des molécules cargo, c'est-à-dire transportées par la vésicule, au moment de l'exocytose.

Dans un premier temps, il faut détecter le lieu et le moment où se produit l'exocytose.

Nous proposons une méthode de classification instantanée simple s'appuyant sur un � la valeur au point p dans la t -ième image. Nous proposons deux modèles dynamiques pour décrire l'évolution spatio-temporelle des vésicules :

• Translation pour les vésicules en phase de transport ou d'amarrage, avec conservation de l'intensité :

∂I ∂t � p, t � = -∇I � p, t � • w , où w = � u, v
� T est la translation de la vésicule segmentée, valable pour tous les points de la vésicule, et ∇ désigne l'opérateur gradient ;

• Diffusion pour les vésicules fusionnant avec la membrane : Ces modèles dynamiques ne sont valides que pour des vésicules présentant un changement de fluorescence dans le temps. Une étape préalable à cette classification consiste à détecter, parmi les vésicules segmentées par ATLAS, celles qui se déplacent ou changent d'apparence, vésicules que nous qualifierons de dynamiques. Une vésicule segmentée à l'instant t est considérée comme dynamique si au moins un de ses pixels subit un changement en t [Boulanger et al., 2010a].

Une fois les vésicules dynamiques détectées, nous estimons les paramètres des modèles (respectivement, les deux composantes de translation u , v et le coefficient de diffusion D ) par une technique de moindres carrés. Pour distinguer l'état dynamique (translation ou diffusion) de chaque vésicule, nous procédons à un test de vraisemblance généralisé pour comparer les deux modèles ainsi estimés et retenons le modèle qui a la meilleure vraisemblance.

Modèle de fusion pour TfR et Langérine

Pour comparer les comportements de TfR et Langérine, deux protéines transmembranaires, nous proposons un modèle de diffusion augmenté d'un paramètre, dit temps P����� III. A������ ��� ���������� �� ������� xiii de résidence, τ, qui coïncide avec la durée nécessaire à une protéine pour quitter la vésicule :

I � p, t � = A 0 /τ σ 2 PSF exp � � � - t τ - � p � 2 2 2σ 2 PSF � � � + � t 0 A 0 /τ 2D (t -u ) + σ 2 PSF exp � � � - u -t 0 τ - � p � 2 2 4D (t -u ) + 2σ 2 PSF � � � du ,
où A 0 est l'amplitude initiale du spot et σ 2 PSF sa variance initiale. Les évènements de fusion membranaire sont plus aisément détectés pour TfR et Langérine que pour Rab11. Ces protéines transmembranaires sont en effet associées à un marqueur pH-sensible. Comme le pH augmente subitement au moment de la fusion, nous détectons fiablement les spots apparaissant au fil de la séquence, toujours à l'aide d'ATLAS. Parmi les évènements de fusion, on observe que certains ne donnent pas lieu à une diffusion. Pour éliminer ceux-ci, nous proposons un test statistique d'adéquation au modèle de fusion.

Une fois les évènements diffusifs détectés, les paramètres du modèle de diffusion sont également estimés par une technique de moindre carrés, mais certaines améliorations sont proposées pour obtenir des résultats performants. En particulier, nous mettons en concurrence plusieurs estimations obtenues avec plusieurs initialisations, et nous retenons celle qui minimise la somme des carrés des résidus.

Des expériences quantitatives sur simulation mettent en évidence les performances de notre méthode d'estimation. Nous traitons par la suite un grand jeu de séquences réelles, pour lesquelles est marquée la TfR ou la Langérine. Ceci nous permet finalement d'exhiber les différences de comportement de ces deux protéines au moment de la fusion, en comparant les distributions du temps de résidence τ dans ces deux cas.

Partie III. Analyse des dynamiques de groupes

Dans cette dernière partie, nous nous intéressons au mouvement collectif d'entités mobiles similaires. Comme nous l'avons expliqué en introduction, même si les méthodes sont applicables à des images biologiques, nous utilisons des vidéos de foule comme support de validation expérimentale.

La vidéosurveillance génère d'énormes quantités de données, qui ne peuvent pas être intégralement visionnées, et nécessitent donc des outils d'analyse automatisés [START_REF] Zhan | Crowd analysis: a survey[END_REF]. En particulier, l'analyse des mouvements de foule dans les vidéos est un domaine extrêmement vaste et en pleine expansion, qui englobe notamment la sécurisation des évènements publiques, la surveillance des lieux des transports en commun et l'analyse des comportements individuels au sein de groupes. xiv Dans ce contexte, nous abordons trois aspects importants :

• la classification de mouvements de groupes [START_REF] Garate | Group tracking and behavior recognition in long video surveillance sequences[END_REF], Solmaz et al., 2012, Wang et al., 2011, Wu et al., 2014, Zhou et al., 2013],

• la recherche des chemins les plus empruntés [START_REF] Wang | Trajectory analysis and semantic region modeling using nonparametric hierarchical Bayesian models[END_REF], Zhou et al., 2011],

• la détection de comportements anormaux [START_REF] Benezeth | Abnormality detection using low-level co-occuring events[END_REF],Kim and Grauman, 2009, Kratz and Nishino, 2009, Mehran et al., 2009].

Pour répondre à ces problématiques, nous proposons une nouvelle approche reposant sur l'analyse d'une paire ou un triplet d'images seulement, alors que les méthodes usuelles d'analyse de foules requièrent des intervalles de temps beaucoup plus longs (au moins une dizaine d'images, et parfois même des minutes de vidéo). Ces approches exploitent en général des cuboïdes spatiotemporels [START_REF] Feng | Online learning with selforganizing maps for anomaly detection in crowd scenes[END_REF], Kratz and Nishino, 2009, Rodriguez et al., 2011], des « tracklets » [START_REF] Garate | Group tracking and behavior recognition in long video surveillance sequences[END_REF], Zhou et al., 2011] ou plus généralement des trajectoires [START_REF] Cheriyadat | Detecting dominant motions in dense crowds[END_REF], Rodriguez et al., 2009, Solmaz et al., 2012, Wang et al., 2011, Zhou et al., 2012]. Contrairement à ces méthodes, les nôtres reposent sur la mesure de champs affines de mouvement 2D entre images successives, calculées sur une collection de fenêtres : nous parlerons d'approche instantanée. Nous évitons ainsi d'une étape coûteuse et souvent difficile de suivi individuel des piétons [START_REF] Baumgartner | Tracking people and their objects[END_REF], Cupillard et al., 2002,Hu et al., 2008,Idrees et al., 2014,Kratz and Nishino, 2012,Rodriguez et al., 2009] ou d'advection de particules à partir d'un flot optique [START_REF] Solmaz | Identifying behaviors in crowded scenes using stability analysis for dynamical systems[END_REF], Wu et al., 2010]. Dans nos méthodes, nous n'utilisons ni intégration temporelle, ni calcul de trajectoire, ni phase d'apprentissage.

Sélection ponctuelle de types de mouvement 2D

Nous abordons les trois problématiques sus-citées avec un descripteur commun de mouvement, basé sur une sélection ponctuelle de types de mouvements affines 2D.

Notre approche n'exploite pas de modèle de comportement humain, et peut s'étendre à l'analyse de groupes d'entités en mouvement : véhicules, animaux, entités biologiques...

Soit W = � W i � i ∈N
une collection de fenêtres de tailles variées et se chevauchant.

Nous considérons trois types de mouvements affines 2D, respectivement translation (T), divergence (D) et rotation (R). Nous notons θ k, i les paramètres du modèle k ∈ � T, D, R � pour la fenêtre W i . En chaque point p = �

x, y � de W i et pour chaque modèle k ∈ � T, D, R � , la vitesse w k, i est donnée par les formules suivantes :

• Pour la translation :

w T,i = � � � 0 0 0 0 � � � � � � x y � � � + � � � b 1 b 2 � � � avec θ T,i = (b 1 , b 2 ) T ;
P����� III. A������ ��� ���������� �� ������� xv • Pour la divergence :

w D,i = � � � a 1 0 0 a 1 � � � � � � x y � � � + � � � b 1 b 2 � � � avec θ D,i = (a 1 , b 1 , b 2 ) T ;
• Pour la rotation :

w R,i = � � � 0 -a 2 a 2 0 � � � � � � x y � � � + � � � b 1 b 2 � � � avec θ R,i = (a 2 , b 1 , b 2 ) T .
Cette façon de mesurer le mouvement permet d'éviter de calculer le flot optique, et facilitera la classification ultérieure du mouvement.

Les paramètres des modèles sont estimés par une méthode robuste dans un schéma d'estimation multirésolution [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF] avec le logiciel Motion2D4 . Nous obtenons ainsi de l'ordre d'une centaine de vecteurs de vitesse candidats par point (trois par fenêtre dans une trentaine de fenêtres). Nous proposons de sélectionner le meilleur candidat selon un critère de maximum de vraisemblance pénalisé basé sur l'hypothèse de conservation d'intensité. La conformité d'un point p au modèle de vitesse de paramètre θ k, i est représentée par la variable aléatoire ε k,i définie par :

ε k, i � p � = I � p + w k, i , t + 1 � -I � p, t � .
En supposant que chaque variable suit une loi gaussienne centrée, et qu'elles sont indépendantes, on peut calculer la vraisemblance jointe dans un patch centré en p. Le meilleur candidat en p est ensuite sélectionné selon le critère d'information d'Akaike corrigé pour tenir compte du petit nombre d'observations fournies par le voisinage local du point p [START_REF] Burnham | Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach[END_REF].

Classification de mouvements d'ensemble

Des trois types de mouvement affine considérés, nous déduisons huit classes de mouvement de groupe en fonction des valeurs des paramètres estimés. Les mouvements divergents correspondent au rapprochement du groupe (divergence positive) ou à son éloignement (divergence négative). Nous distinguons par ailleurs les mouvements de rotation directe et indirecte, et quatre directions principales de translation pertinentes dans le référentiel de l'image : nord, ouest, sud et est. Nous obtenons ainsi une carte de classification préliminaire � c , forcément bruitée. Pour obtenir la classification finale C des mouvements de foule, nous régularisons � c par une méthode de votes majoritaires. De plus, on voit souvent apparaître côte-àcôte dans � c des mouvements de divergences ou de rotations opposés. Pour contourner cette difficulté, nous introduisons un processus de vote en deux tours. Au premier xvi tour, nous introduisons des interactions entre classes, appelées inhibition et renforcement ; il permet de sélectionner le type de mouvement : translation, divergence ou rotation. Le second tour permet d'affiner la classification pour obtenir la carte finale régularisée avec huit classes de mouvement de foule.

Détection des chemins principaux

Connaissant le type de mouvement affine sélectionné en chaque point, nous déduisons le vecteur de vitesse en chaque point. Nous obtenons ainsi une approximation pertinente du flot optique. Nous estimons alors la vitesse moyenne dans chaque cellule d'une grille régulière et sur un court intervalle de temps. Cette mesure suppose un mouvement stationnaire dans la vidéo traitée, hypothèse cohérente avec la recherche des chemins les plus empruntés dans la scène. Pour retrouver le chemin le plus emprunté, nous discrétisons les directions de ces vitesses moyennes (en huit secteurs) de façon à pointer vers une cellule voisine.

Nous obtenons ainsi un graphe orienté où les noeuds correspondent aux cellules de la grille et les arêtes sont définies par les directions locales du flot optique moyen.

Ce graphe est en fait réduit à un arbre, car de chaque noeud ne part qu'une arête au plus. La recherche du plus long chemin dans ce graphe est donc très simple (il suffit de vérifier qu'il n'existe pas de cycle). Sous l'hypothèse d'un mouvement stationnaire, nous détectons ainsi les chemins les plus suivis dans la scène. extrêmes généralisées [START_REF] Embrechts | Modelling extremal events for insurance and finance[END_REF], dont nous estimons les paramètres P����� III. A������ ��� ���������� �� ������� xvii avec la méthode décrite dans [Boulanger et al., 2010a]. Nous déduisons ainsi un seuil sur les distances à partir d'une probabilité de fausse alarme choisie par l'utilisateur, qui contrôle précisément le taux de fausses détections.

Détection de comportements anormaux

Les expériences sur des séquences simulées et réelles mettent en évidence quantitativement l'efficacité de la méthode, tant en termes de détection que de localisation des anomalies.

Contents

Résumé en français vii

Partie I. I � ���� ������, we investigate statistical methods to detect, estimate and charac- terize dynamical events in image sequences. There are two cases:

• Studied objects do not interact, and individual dynamics can be independently analyzed;

• Studied objects interact, and group dynamics must be analyzed as a whole.

In the case of individual dynamics, our primary focus is on biological image sequences showing proteins evolving in a cell, and more precisely at the cell frontier named plasma membrane. The images are acquired by the means of total internal reflection fluorescence microscopy (TIRFM), an observation technique dedicated to plasma membrane dynamics analysis. We model the protein dynamics and estimate the biophysical parameters in TIRFM image sequences for further biological analysis.

We first propose a new spot detection method with automatic scale selection, aimed at localizing dynamical events of interest.

Eventually, we focus on image sequences at lower magnifications, that is, depicting groups of cells, for example cell tissues, instead of an isolated cell. Since they are free of appearance model, the developed methods are quite general and extended to other applications including crowd motion analysis in videos.

Whether it is for spot detection, protein dynamics estimation or group motion analysis, a common philosophy is ubiquitous, however. First, statistical arguments are used to automatically infer the method parameters. Therefore, the method does not need to be specified for every single image sequence, and the method can be used easily by neophytes in image processing.

Context and motivation

Human body is composed of 10 13 to 10 14 cells, millions of which are renewed every second. Correct cell functioning is an essential prerequisite to life, and a segment of biology in which much remains to be discovered. A key actor of life understanding, cell biology lies at the intersection of the very applied and fundamental life sciences, for example as a cornerstone of evolution theory. As such, it also interacts with many research fields, ranging from therapy-related fields (e.g., medicine, drug development, epidemiology) to technology-related ones (e.g., optics, electronics, physics),

and relies more and more on computer sciences and applied mathematics, specifically image processing.

In order to detect variations between normal and pathological situations, intraand inter-cellular dynamics must be understood. This is done through the observation and the characterization of molecular processes, among others, which can be tackled in various ways, from biochemistry to microscopy imaging. We are interested in the latter, which encompasses microscope technology, genetics engineering and processing tools. While optical and biological sciences are out of our reach, we are developing image processing methods, which can profit the biology research.

A need for automatic methods

Key role of proteins in the recycling process

Vital processes of the cell, from feeding to motility and reproduction, are subject to the precise equilibrium between the different chemical compounds present in the cell.

The equilibrium state is named homeostasis. While from Latin homeo, equilibrium, C������ 1. G������ ������������ and stasis, static, the term refers to a very dynamical process, or more precisely a set of interacting processes. Absolute and relative quantities of proteins, sugars, ions, among others, are therefore continuously regulated by the cell. The control of the total amount of a given compound is carried out by generating incoming and outgoing flows. Hence, the processes occurring at the cell frontier, the plasma membrane, are of utmost importance to the cell life.

Among various and complex dynamics taking place to ensure homeostasis, exocytosis is the process by which compounds are expelled out of the cell, by crossing the plasma membrane. While smallest structures like ions can autonomously diffuse through the membrane, larger ones rely on transport and exocytosis intermediates called secretory vesicles, or recycling vesicles in our specific case study. To fulfill their function, the vesicles are involved in a variety of processes, the latter being vesicle fusion to plasma membrane, when their contents is effectively expelled. Numerous actors are working to ensure the smooth functioning of exocytosis, most notably proteins, which act and interact at every level of the process in a very strict order.

Proteins are large molecules made of amino acids playing a key role in almost every cell functions. Modeling proteins dynamics and interactions is thus a major component of the cell biology understanding. In particular, proteins are at the heart of every stage of the process, encompassing the vesicle formation, its transport, tethering, docking and, ultimately, fusion.

In this thesis, we study the behavior of three of them: Rab11, Transferrin receptor (TfR) and Langerin. Rab11 is known to be involved in numerous aspects of exocytosis, but pieces of the puzzle are missing, especially regarding the vesicle fusion, and diffusion of the protein. TfR dynamics are well established. As a transmembrane protein, it is constrained to stay in the membrane, thus reducing its degrees of freedom, which facilitates its study. However, while Langerin is also a transmembrane protein, a different behavior is observed in the image sequences, which raises questions about the underlying mechanics. To bring to light quantitative differences between TfR and Langerin, dynamical models must be proposed, estimated and compared. Characterizing the dynamics and discovering the roles and interactions of Langerin or Rab11 is a priority of UMR 144.

Group motion is more than a collection of individual motions

In biological images with smaller magnification than in Part II, several cells move collectively as a single population. As an example, during embryo development, a group of cells divide, move together and interact. Also, in image sequences depicting cell tissues, cells move and deform, thus imposing constraints to their neighbors. In these cases and many others, the motion of single cells is embedded in a higher-level coherent group motion, that of the whole embryo or tissue.

O��������� ��� �������������

In biology, group motion analysis is generally tackled by estimating the dynamics of each cell individually. Actually, many methods were proposed to study the dynamics of embryos or tissues by detecting and tracking individually each cell. The tracking of individual cells is then performed in order to further derive the relevant information at the group level. However, it is not straightforward to take into account the interactions in individual dynamical models without globally estimating the whole group dynamics. Therefore, it is natural to directly study the dynamics of the whole group, that is, not to differentiate the individual cells.

While not used in biology so far, the latter approach is very common to analyze pedestrian group motions or vehicle traffic, among others. For this reason, we extend the scope of the last part of the thesis to crowd motion analysis, a very active field from where approaches can be adapted to analyze the motion of cell groups.

Objectives and contributions

Detection of individual spatiotemporal events

As stated earlier, we target two types of dynamics, namely that of independent individuals and that of groups. We first investigate the modeling of a cell regulation process named exocytosis and related control parameters estimation from microscopy image sequences. We specifically focus on the last stage of the process: diffusion. This investigation encompasses the detection of the exocytosis in the image sequence, which itself relies on spot detection because vesicles appear in the images as small spots. To this end, we design a spot detection method with very few parameters which takes advantage of a specificity of the image used for exocytosis characterization, that is, spots to detect have similar sizes in TIRFM image sequences. As opposed to classical spot detectors, we propose an approach for automatically selecting the best detection scale. Then, the detection amounts to a local thresholding which automatically adapts to the local image statistics. The detection sensitivity only depends on a userspecified probability of false alarm, which is set for the whole sequence (or set of sequences) to be processed. For the thresholding to be efficient, the spots in the input image have to be enhanced with a filter at a particular scale, whose choice is critical and is explicitly handled as aforementioned.

Modeling and estimation of fusion dynamics

As for the vesicle fusion characterization, we propose new diffusion models suited to Rab11, TfR and Langerin, which are more realistic than the classically used point source model. Indeed, in the latter, it is assumed that the vesicle size is smaller than the pixel pitch, and that its content is instantaneously released. We relax those hy-C������ 1. G������ ������������ potheses, enabling a spatial extent for the vesicle, and a continuous release, thus better fitting the data. Moreover, in the case of Rab11, it is not known if diffusion occurs in two and/or three dimensions. Therefore, we investigate both possibilities. Once detected, the fusion events with the aforementioned spot detector, the model parameters must be estimated. We rely on an estimation method able to cope with the complexity of the data. Actually, the number of data points is quite low with respect to the noise level, so classical diffusion coefficient estimation methods behave poorly for vesicle fusion dynamics. Furthermore, the biochemical parameters are then compared for different experimental conditions. Parameter statistics are extracted from the whole set of detected fusion events in sequence collections.

This leads us to the last part of the thesis. In most biological applications, the behavior of a single sample (e.g., the dynamics of one fusion event in the cell) is not the critical issue. Instead, it has to be compared to the mean, normal or dominant behavior (e.g., the average diffusion coefficient estimated on a set of sequences). Sometimes, the behavior is expected to be similar for all samples, so that coherent or consistent dynamics can be expected. This knowledge can be exploited to improve modeling and estimation, as done to select the optimal scale for spot detection. Therefore, the group properties (velocity statistics, common behavior, count...) are central in the study under concern.

Characterization of group dynamics

To capitalize on the strength of crowd motion analysis methods, we propose to adopt the philosophy in biological image sequence analysis. Specifically, we aim at characterizing the cell crowd dynamics without distinguishing individual cells. However, we do not simply transpose crowd motion analysis methods to biology: to study human crowds or road traffic, existing methods generally use hundreds of frames, or even the whole image sequence, which comes with a very high computation cost. In contrast, we develop instantaneous methods for group dynamics analysis in the remaining of the thesis, that is, methods using only a couple of consecutive frames to characterize the crowd dynamics.

Moreover, in order to gather more information on motion than classical methods based on optical flow, we design a method which provides quantitative and qualitative motion information at the same time. The motion is estimated as a labeled affine flow, which is a map conveying both a velocity vector and motion type information. Statistical arguments are used to select a relevant affine motion model at each point of the image, which can be classified owing to simple rules on the affine parameters.

Then, the estimates are used to characterize the group motion, or to detect and segment local anomalous dynamics inside the group. To validate the methods, we rely 1.3. O����������� �� ��� ������ on available benchmarks, which generally gather crowd and vehicle traffic videos, and on biological images showing groups of cells.

Organization of the thesis

After an introduction to basic knowledge on the cell and microscopy, the thesis is or- The thesis organization is synthetically presented below.

Chapter 2. On the cell and its observation 13 The cell organization and components are described to introduce the thesis.

The observation model is formulated while we describe the microscopy technique used throughout the thesis, namely TIRFM. Quantitative evaluations and comparisons demonstrate the efficiency and robustness of the method. Real TIRFM images are analyzed to exhibit the differences between TfR and Langerin dynamics.

Part I. Spot detection

Part III. Group dynamics

Chapter 9. Labeled affine flow for group motion characterization 127 Labeled affine flow is proposed as an augmented motion data for group motion analysis. At each point, the motion measurement conveys two types of information: a quantitative motion vector and a qualitative motion type. As a proof of concept, we rely on the proposed labeled affine flow to characterize the group motion, and to recover the dominant paths followed by the group.

Chapter 10. Anomaly detection and localization 147 The proposed labeled affine flow is also used to detect local abnormal events. A patch-based approach allows us to detect and locate anomalies, by comparing local motion class histograms in a center-surround setting.

Chapter 11. General conclusion 163 Main results are summarized. We discuss the contributions of the thesis with a critical analysis of all the developed methods. This leads us to envisage improvements and foresee future work.

D�������� ��� ���������������� �� ��������� ������ �� ����� ��������� 13

2

On the cell and its observation T � ������� image sequences showing cellular dynamics, a basic knowledge of the cell and imaging system is required. We shall not plunge into an exhaustive enumeration of the cell components, but some primary structures must be presented to understand the context of our work, as well as the microscopy technique, which actually relies on genetic processes. Indeed, numerous biological processes and associated cell compartments take care of parsing genetic information, contained in the cell, to infer which proteins, vital compounds of the cell, should be synthesized. Proteins are precisely at the heart of both the project and observation system.

Along with nucleic acids, lipids and sugars, proteins are the very basic organic components of the cell. They are involved in a wide variety of intra-and intercellular processes, where they carry out numerous functions: they generate mechanical forces, transport chemical materials, are key actors of homeostasis and cell defense, catalyze chemical reactions... Regarding the research project, three proteins are studied in particular, namely Transferrin receptor (TfR), Langerin and Rab11. More precisely, we will focus on their dynamics in the so-called exocytosis-recycling process, which is part of the homeostasis regulation process.

To be observed in total internal reflection fluorescence microscopy (TIRFM), proteins we are interested in are genetically associated to fluorescent tags. When excited, they emit light detected by the microscopy setup, thus allowing us to locate and track the proteins.

The chapter is composed of three sections. A quick overview of the cell biology chronology is given in Section 2.1, thus briefly introducing vital processes of the cell. Biological processes involved in protein synthesis and in exocytosis are then detailed in Section 2.2. Finally, the process of fluorescent tagging and TIRFM are presented in Section 2.3.

A bit of history

Today's cell biology knowledge is the result of a 350-year research. Back in the XVIIth century, Robert Hooke -an architect, philosopher, mathematician, physicist, and an optical engineer in his spare time -builds one of the first microscopes ever made to study the porosity of various materials [Magner, 2002]. Analyzing cork bark among other surfaces, he observes a tessellation made of tiles he terms cells [START_REF] Hooke | Micrographia: Or Some Physiological Descriptions of Minutes Bodies made by Magnifying Glasses with Observations and Inquiries thereupon[END_REF] by analogy to honeycomb cells (Figure 2.1). A few years later, specifically studying living beings, Antonie van Leeuwenhoek observes cells in various organisms: bacteria, algae, muscle tissues... [Hoole, 1800] But no generalization is made yet.

It will be two centuries before Matthias Jakob Schleiden draws the basement of modern cell theory, defining every plant as a set of cells [Schleiden, 1838]. The following year, Theodor Schwann shows that the theory also applies to animals [Schwann, 1839]. Cell theory basements are completed in 1852 by Robert Remak, who demonstrates that every cell originates from the division of a preexisting cell [Remak, 1852].

In turn, not only is cell a component of every living being, but also a living being itself: cell autonomously carries out the functions necessary to live and reproduce. Among others, it contains the genetic information, synthesizes complex vital molecules from basic compounds, and regulates its content to stay in homeostasis, that is the equilibrium state of its content: sugars, lipids, proteins, ions... In the second half of the century, genetic sciences emerges, pioneered by Gregor Mendel, who formalizes discrete traits inheritance in 1865 by identifying dominant and recessive traits [Mendel, 1866], and later Friedrich Miescher, who discovers deoxyribonucleic acid (DNA) in 1869 [Miescher, 1871]. The link between both, however, that the genetic information is contained in the chromosomes. The proof is given in 1915 by Thomas Hunt Morgan [START_REF] Morgan | The Mechanism of Mendelian Heredity[END_REF], a discovery for which he is awarded the Nobel Prize in Physiology or Medicine. In 1944, Oswald Avery finally locates the genetic information in DNA [START_REF] Avery | Studies on the chemical nature of the substance inducing transformation of pneumonococcal types. Induction of transformation by a deoxyribo-nucleic acid fraction isolated from pnuemococcus type III[END_REF].

Genetics understanding eventually leads to fluorescence microscopy, which will be described in more details hereafter. In 2008, Osamu Shimomura, Martin Chalfie and Roger Tsien receive the Nobel Prize in Chemistry for the discovery of the green fluorescent protein (GFP) and its application to cell biology [START_REF] Chalfie | Green fluorescent protein as a marker for gene expression[END_REF]. True to its name, the GFP is a protein which emits green light when excited, e.g., by an incident laser beam. By synthesizing a DNA fragment, they managed to bind GFP to another protein naturally present in the cell, and follow the latter in the living cell using an optical microscope.

The number of discoveries owing to fluorescence microscopy is phenomenal and continues growing fast. Ground-breaking technologies now allow combining the dynamical acquisition of light microscope with the resolution advantage of electron microscopy by freezing the cell during the experiment. In the meantime, the envelope of optical microscopy is expanded every year using multiple-exposure techniques like structured illumination microscopy (SIM) [Neil et al., 1997,Barlow andGuerin, 2007], which shall approach electron microscopy resolution in the near future, while still enabling dynamical recording. C������ 2. O� ��� ���� ��� ��� �����������

ADN A T G G C G T T C A G A A C T G A T A C G T A A | | | | | | | | | | | | | | | | | | | | | | | | Transcription ARN A U G G C G U U C A G A A C U G A U A C G U A A �������� ��������

Organization and processes of the cell

Let us now briefly describe the cell and its organization, by focusing on the relevant biological processes for the thesis. In Figure 2.2, the cell is represented along with some organelles which will be presented below. Cells can be either prokaryotic or eukaryotic. Eukaryotic cells (from Greek eu, true, and karyon, nucleus), which we study, contain a nucleus and cytoplasm enclosed within a plasma membrane. In contrast, prokaryotic cells (from pro, before) lack nucleous. The cytoplasm itself is mainly composed of an aqueous solution named cytosol, in which numerous organelles carry out specific functions.

Additionally, the cell is supported by the cytoskeleton, a set of molecular structures which give its shape and rigidity to the cell. Cytoskeleton is also an actor of numerous dynamical process, some of which will be presented in Section 2.2.3.

From nucleus to cytosol: From genetic information to proteins

Nucleus is the host of most of the genetic material of the cell. A nuclear envelope surrounds the nucleolus and nucleoplasm, the latter containing DNA. Essentially, DNA is a long string of nucleotides, the atomic bricks of genetic information, which can be of four types: Adenine (A), Thymine (T), Cytosine (C) or Guanine (G). Fragments of DNA, genes are sequences of nucleotides that define which proteins are synthesized by the cell, as explained hereafter.

The predominant component of the cell in volume, cytosol is the fluid in which nucleus, organelles and cytoskeleton lie, and the place where proteins are synthesized. The so-called protein biosynthesis process is illustrated in Figure 2.3. To produce a protein, DNA is first transcripted to ribonucleic acid (RNA), a molecular structure similar to DNA that can leave the nucleus to convey genetic information to the cytosol. In RNA, Thymine is replaced with Uracil (U), as shown in Figure 2

The plasma membrane, a place of exchange

Far from being a simple frontier, the plasma membrane is involved in different aspects of the cell life, from mechanical structure to motility and through chemical regulation, which we are mostly interested in.

As schematized in Figure 2.4, the membrane is composed of phospholipids, molecules with a hydrophilic-lipophobic head and hydrophobic-lipophilic tails.

Phospholipids are arranged in two layers where the hydrophilic head points toward the aqueous cytosol or in the aqueous medium outside the cell, which isolates the hydrophobic tails inside the membrane.

Every chemical flow entering or leaving the cell must pass through the plasma membrane. Thus, the membrane is of first importance in the protection from extracellular attack and regulation of the cell homeostasis. The plasma membrane is said to be selectively permeable, because only specific material can pass through, sometimes in limited quantity. While small molecules can move across the membrane, either by passing between phospholipids, or by specific biological gates, a more complex process is required to absorb or secrete bigger structures. The latter process, studied in Part II, relies on transport intermediates named vesicles. They will be presented in more details in the next subsection, and later in Part II.

The cytoskeleton and intracellular transport

Apart from the plasma membrane, the cell mechanical properties are due to the cy- Indeed, functions of the cytoskeleton are not restricted to structural properties.

For example, dynamics of the cytoskeleton are involved in the cell motility and cell division in the mitosis process. However, we are particularly interested in the conveying functions of the microtubules and Actin filaments, and more precisely in vesiclemediated exocytosis.

In the second part of the thesis, we study the dynamics of vesicles during exocytosis. A vesicle is a closed bag formed by a lipid bilayer membrane similar to the plasma membrane, as depicted in ported by the vesicle, from a donor compartment, or from the extracellular space in the case of endocytosis, to a target compartment, or to the extracellular space in the case of exocytosis. To this end, the vesicle forms from the membrane of donor organelle at the beginning of the transport, and fuses to the target membrane at the end of the transport. In between, the vesicles is moved by molecular motors along microtubules, which act as highways in the cell. Then, in the case of exocytosis transport, in order to reach the plasma membrane for expelling its content, the vesicle has to make its way through the Actin meshwork, before fusing to the membrane to release its content.

Furthermore, several types of structures travel with the vesicle, including:

• Transmembrane proteins, which span across the vesicle membrane during transport and are released in the target membrane after fusion;

• Vesicular motors, which pull the vesicle along the cytoskeleton. Vesicular transport will be described in more details in Chapter 6, as it is a fundamental component of exocytosis, that is the subject of Part II. Specifically, we will focus on two transmembrane proteins, TfR and Langerin, and on Rab11, a protein which is associated to the vesicular motor during transport.

TIRFM for imaging membrane dynamics

Fluorescence microscopy

As the name implies, total internal reflection fluorescence microscopy (TIRFM) is part of the fluorescence microscopy acquisition techniques. The latter is a type of microscopy where the structure of interest, some protein in the frame of the thesis, is tagged with a fluorescent compound called fluorophore. By tagged, we mean that the fluorophore, here fluorescent proteins, bounds to the protein of interest, so that locating the fluorophore roughly corresponds to locating the protein itself. As mentioned in Section 2.1, fluorophores emit photons when excited by incident light. Thus, in fluorescence microscopy, the tag itself emits the light which is collected by the microscope sensor. Structures of interest can hence be selectively observed, while all other structures remain invisible to the microscope sensor.

As for the imaging device, the principle of the fluorescent microscope is schematized in Figure 2.7, which represents an inverted microscope, meaning that the objective points upward. The objective top is immersed in some optical oil that has the same refractive index n as the front lens element and cover slip, about 1.5. This ensures that the contact surface between the cell and cover slip is the sole optical interface between the microscope objective and the cell itself. Nevertheless, fluorescence has drawbacks as well, principally phototoxicity and photobleaching, which are closely related to the laser beam power. Phototoxicity is the killing of certain cellular structures due to the relatively high-energy photons hitting them, which affects cell functioning and eventually leads to its death. In order to preserve the functions of the cell during the image sequence acquisition, a compromise must be made between the signal-to-noise ratio (SNR) and toxicity, that is, the laser power must be carefully set.

Photobleaching is the phenomenon by which excited fluorophores can either emit a photon or enter a chemical reaction which eventually kills its fluorescence ability.

Over time, more and more fluorophores are bleached and do not emit light anymore.

Therefore, bleaching is visible in the images as a slow fading of global intensity. Hopefully, we will see that it is mostly negligible in TIRFM.

Total internal reflection

In TIRFM, the fluorophores are excited by a so-called evanescent wave [Axelrod, 2008]. The latter forms inside the cell when the incident light beam is totally reflected at the interface between the cover slip and the cell itself. Total reflection occurs when the incident angle of the laser beam exceeds a critical angle. This is illustrated in Fig-

ure 2.8. Energy of the evanescent wave decreases exponentially with the distance to the cover slip, so that practically no energy is transmitted above a certain distance called penetration depth; its order of magnitude is 10 2 nm [START_REF] Steyer | A real-time view of life within 100 nm of the plasma membrane[END_REF].

Therefore, in TIRFM, only fluorophores close to the coverslip are illuminated and excited. This allows to make visible only the molecules that are inside or very close to the plasma membrane -approximately up to the penetration depth. As a consequence, TIRFM is particularly well suited to study the dynamics occuring at the plasma membrane, such as vesicle fusion studied in Part II [START_REF] Reichert | Total internal reflection fluorescence (TIRF) microscopy. I. Modelling cell contact region fluorescence[END_REF], Deng et al., 2009, Letinic et al., 2010, Burchfield et al., 2010].

Yet, TIRFM scores a better spatial and temporal resolution than classical optical microscopy techniques. Moreover, as the light energy transmitted to the cell is confined to a very thin region, the power of the incident laser beam can be quite low to get a sufficient intensity at the plasma membrane. Compared to classical fluorescence microscopy, the energy transmitted to the plasma membrane can be higher without being toxic, so TIRFM features a very high SNR among fluorescence microscopy techniques. For the same reason, photobleaching is much reduced, so that when using local approaches which rely on few consecutive frames, it can merely be neglected.

To acquire a color image (e.g., Figure 2.5), fluorophores with different emission wavelengths can be used to tag different proteins. Using a color splitting technique named Dual view [START_REF] Gidon | A Rab11A/myosin Vb/Rab11-FIP2 complex frames two late recycling steps of langerin from the ERC to the plasma membrane[END_REF], the wavelengths are then projected side-byside on the sensor using a set of mirrors or prisms. This technique can be employed to compare the localization of different structures, or to use one structure property to ease working on another structure, as done in Chapter 7.

Point spread function

Point spread function (PSF) is the optical transfer function. Even for an ideal microscope, the image of a point in focus by the optical system is not a point but rather a blurred spot. The more out-of-focus the point, the wider the spot. Theoretically, a point in focus is projected through a pinhole or diaphragm as an Airy disc [Airy, 1835, Marian et al., 2007]. Naturally, the model becomes more complex with more complex optical formulas.

Moreover, even with a theoretical, ideal lens, a PSF model should take into account the discrete nature of the microscope sensor. Indeed, the photons collected at a given pixel are spread over the photosite area. While small in absolute terms, the photosite area is not negligible compared to the PSF width. Therefore, instead of considering a sampled Airy function to model the amount of light gathered by a photosite, the function should be spatially integrated over the pixel surface [START_REF] Small | Fluorophore localization algorithms for super-resolution microscopy[END_REF].

Then come optical aberrations introduced by the lens, such as astigmatism and coma. Not only are these aberrations hard to infer from optical specifications, but they also vary spatially [Aguet, 2009]. The variation theoretically only depend on the distance from the optical center, but real optics are not precisely enough manufactured to stick to this assumption. Therefore, in order to precisely model the PSF over the whole image, classical methods consist in fitting parametric models over the image domain, and then interpolating the parameter field in order to get a continuous PSF map [Aguet, 2009].

In practice, apart from very specific applications as deconvolution [Sibarita, 2005, Sarder andNehorai, 2006] or super-resolution [START_REF] Carrington | Superresolution three-dimensional images of fluorescence in cells with minimal light exposure[END_REF], Zhang et al., 2006], a Gaussian approximation is often adequate to model the PSF in TIRFM images [START_REF] Small | Fluorophore localization algorithms for super-resolution microscopy[END_REF], so we rely on this approximation throughout the thesis.

Furthermore, in TIRFM, objects at positive heights are out of reach for the evanescent wave. Therefore, they do not appear in the image and all visible objects are in focus; if not, they are at about the same distance from the focus plane. Therefore, objects of similar size in the cell are projected as spots of similar size in the TIRFM image.

Digital image model

Notation and units

The image sequence, denoted I , is a succession of digital photographs termed frames.

Frame at time t will be denoted as I (t ). When not otherwise stated, temporal unit will [START_REF] Zhang | Multiscale variance-stabilizing transform for mixed-Poisson-Gaussian processes and its applications in bioimaging[END_REF]. Actually, this model is not restricted to fluorescence microscopy, and is generally assumed in digital images when the following two random components are taken into account:

• Photon noise, which results from the stochastic process of photon emission and conversion;

• Electronic noise, which is generated by the sensor and electronics.

In As for electronic noise, it results from different sources, yielding the following main components:

• Dark-current noise reflects the generation of electrons by the sensor even in the absence of photons, due to the non-null electron excitation in the sensor and electronics;

• Similar to photon noise, electronic shot noise accounts for the quantization of the electric potential levels, as the electric charge representing the pixel signal is quantified as the number of electrons;

• Flicker noise, which depends on the signal frequency, arises from direct current passing through the photosite.

Generally, as it is extremely difficult to tell apart the different intermingled components, the whole electronic noise is modeled as a Gaussian process [Boulanger, 2007].

Finally, the signal is amplified by a multiplicative gain g 0 . As a result, the intensity

I � p, t
� of a pixel p at time t is modeled as [Boulanger, 2007]:

I � p, t � = g 0 n � p, t � + ε � p, t

Variance stabilization

Most image processing methods assume a white Gaussian noise. Therefore, they are not well suited to the MPG model. To deal with fluorescence microscopy images, an intensity transform is often used, which gives a white Gaussian noise in the transformed images [START_REF] Zhang | Multiscale variance-stabilizing transform for mixed-Poisson-Gaussian processes and its applications in bioimaging[END_REF].

In the remaining of the thesis, we will rely on the generalized Anscombe transform (GAT) T , defined as [START_REF] Boulanger | Patch-based non-local functional for denoising fluorescence microscopy image sequences[END_REF]:

∀ p ∈ Ω I , ∀ t ∈ N , T � I � p, t � � = 2 g 0 � g 0 I � p, t � + 3 8 g 0 + σ 2 ε -g 0 µ ε . (2.2)
Parameters g 0 , µ ε and σ 2 ε will be estimated with a patch-based method proposed in [Boulanger et al., 2010a].

Other variance stabilization methods exist, as the multiscale variance stabilizing transform (MS-VST) proposed in [START_REF] Zhang | Multiscale variance-stabilizing transform for mixed-Poisson-Gaussian processes and its applications in bioimaging[END_REF], but the GAT is computationally lighter, and is sufficient in the situations we encounter in the thesis, as demonstrated by the experimental results reported in the remaining of the thesis.

I Spot detection

The landscape of spot detection S ���� ��� ����� ���� of protein tagging with green fluorescent protein (GFP) [START_REF] Chalfie | Green fluorescent protein as a marker for gene expression[END_REF], microscopy investigations at the single cell level have been faced with the problem of determining the location and behavior in space and time of spots, such as microtubule end tips, adhesion molecular complexes, or vesicles as illustrated in Figure 3.1. Detecting such subcellular particles in fluorescence microscopy is indeed of central interest for further quantitative analysis as particle counting [START_REF] Byun | Automated tool for the detection of cell nuclei in digital microscopic images: Application to retinal images[END_REF], particle pattern recognition [START_REF] Jackson | Model building and intelligent acquisition with application to protein subcellular location classification[END_REF],

particle tracking [START_REF] Sbalzarini | Feature point tracking and trajectory analysis for video imaging in cell biology[END_REF], Miura, 2005, Jaqaman et al., 2008, Chenouard et al., 2014] or dynamics classification [START_REF] Sage | Automatic tracking of individual fluorescence particles: Application to the study of chromosome dynamics[END_REF], Boulanger et al., 2010a, Sironi et al., 2011, Basset et al., 2014b]. All these subcellular analyses C������ 3. T�� ��������� �� ���� ��������� start with a reliable, accurate and efficient detection of particles in fluorescence microscopy images.

Our goal is to segment exocytotic vesicles in cell images acquired in total internal reflection fluorescence microscopy (TIRFM). Among fluorescence microscopy image modalities, TIRFM is the perfect tool to investigate processes occurring close to or at the cell surface such as endocytosis and exocytosis processes [START_REF] Reichert | Total internal reflection fluorescence (TIRF) microscopy. I. Modelling cell contact region fluorescence[END_REF], Deng et al., 2009, Letinic et al., 2010, Burchfield et al., 2010]. The physical size of exocytotic vesicles spans across a limited range. Given the limited depth of field (DOF) of TIRFM, the variation of the scale of these fluorescently labeled objects in the 2D images is also limited. In this part, we will focus on M10 cells showing the cargo proteins Langerin and Transferrin receptor (TfR) tagged with pHluorin, or the Rab11

GTPase tagged with mCherry. These proteins are associated to transport intermediates such as vesicles recycling to the cell surface and appearing as bright spots, which can be round or elongated, as depicted in Figure 3.1. Another application of the presented method could be the identification, detection and quantification of adhesion molecular complexes, in cells migrating or not. These biological architectures are relatively small and regular at the single cell, composed of multiple molecular partners.

As a consequence, it is worth developing a spot detection method able to automatically find the average object size or the most frequent one. We propose a segmentation framework with automatic scale selection and local adaptive thresholding. Our method exploits the Laplacian of Gaussian (LoG) of the intensity image and automatically detects the characteristic scale of the objects of interest. To cope with inhomogeneous background, thresholding is adapted to local statistics, while a single probability of false alarm (PFA) is set for the whole image or even the collection of images to be processed. In short, we will automatically infer from image data the optimal parameters usually left to the user guidance in other methods, that is, LoG scale and detection threshold. We name adaptive thresholding of LoG images with auto-selected scale (ATLAS) the method described in this part.

ATLAS comprises several significant improvements and extensions compared to the preliminary method SLT-LoG we introduced in [Basset et al., 2014a]:

• We now resort to a discrete filter for the scale-space representation and we can deal with any arbitrary scale, i.e., with scales of any precision;

• We have designed four original scale selection criteria;

• We have produced and made publicly available a new benchmark dataset for spot detection methods;

• We have conducted an extended comparative evaluation with existing methods on several datasets, and we have evaluated our method on a larger range of real images.

While our primary goal is to detect exocytic vesicles in two-dimensional TIRFM images, the ATLAS method can be applied to other types of images as well, provided objects to be detected are of similar size in the image or of a couple of sizes at most.

Comparisons of spot detection methods were reported in [START_REF] Smal | Quantitative comparison of spot detection methods in fluorescence microscopy[END_REF], Ruusuvuori et al., 2010], providing with a broad overview of state-of-the-art methods.

However, the dataset used in these previous experiments remains limited in terms of content and challenges. Indeed, real TIRFM images are far more complex than images of this dataset, specifically, the signal-to-noise ratio (SNR) is generally lower in real images and objects to be detected are smaller and often darker. We have then constructed a more realistic and more challenging dataset with ground truth to quantitatively evaluate and compare methods. In addition, we have used images supplied by the simulators designed in [START_REF] Rezatofighi | A framework for generating realistic synthetic sequences of total internal reflection fluorescence microscopy images[END_REF], Boulanger et al., 2009].

Part I is organized as follow. The present chapter gives an overview of existing spot detection methods in fluorescence microscopy, and compares the different approaches for setting the detection scale(s). Our automatic scale selection is presented in Chapter 4, and the actual adaptive detection is described in Chapter 5, where we also compare our detection results with those of state-of-the-art methods.

Existing spot detectors

In [START_REF] Smal | Quantitative comparison of spot detection methods in fluorescence microscopy[END_REF], the authors provide a broad panorama of spot detection methods, and thoroughly evaluate the performance of a dozen methods. As explained by [START_REF] Smal | Quantitative comparison of spot detection methods in fluorescence microscopy[END_REF], the common detection framework consists in first denoising the image and enhancing the spots to be detected. Then, highest or lowest values of the enhanced signal, corresponding to spots, are extracted.

The simplest way of detecting spots in a gray level image is to threshold the image intensities from the intensity histogram. The threshold value can be automatically selected by techniques such as Otsu's method [Otsu, 1979] or entropy minimization [START_REF] Kapur | A new method for gray-level picture thresholding using the entropy of the histogram[END_REF], Sahoo et al., 1997]. However, a single global threshold cannot tackle complex images where variation in background intensities may exceed spot intensity magnitude. Therefore, numerous space-varying thresholding methods were proposed [Sahoo et al., 1988, Sezgin and[START_REF] Sezgin | [END_REF]. In particular, local threshold values are deduced from local statistics to detect cell nuclei in [START_REF] Phansalkar | Adaptive local thresholding for detection of nuclei in diversity stained cytology images[END_REF].

More advanced methods, such as detectors based on top hat filter (TH) [START_REF] Bright | Two-dimensional top hat filter for extracting spots and spheres from digital images[END_REF]Steel, 1987, Breen et al., 1991] or LoG filter as in the spot enhancing filter (SEF) method [START_REF] Sage | Automatic tracking of individual fluorescence particles: Application to the study of chromosome dynamics[END_REF], not only smooth the image, but also enhance the underlying signal. More specifically, the LoG filter (which we will rely on) is a band-pass filter which enhances objects of a particular size, reduces noise and lowers low-frequency background structures. Sage et al. [START_REF] Sage | Automatic tracking of individual fluorescence particles: Application to the study of chromosome dynamics[END_REF] empirically observed that the C������ 3. T�� ��������� �� ���� ���������

LoG filter is close to the optimal whitened matched filter for Gaussian spots in fluorescence microscopy images, that is, the SNR of the filtered image is maximized at the spot center. Yet, the choice of the LoG variance is critical and highly dependent on the spot size. Similarly, the bandwidth of the TH filter is adjusted with two critical parameters, the top and brim radii. Ideally, they should correspond to the spot size and distance between neighboring spots, respectively. In the so-called morphological top hat filter (MTH) version of TH [Soille, 2003], the image background is estimated by an opening operation which removes objects smaller than the structuring element.

In order to reduce noise, a Gaussian blur is initially performed. The background estimate is then subtracted to the image to detect spots by thresholding.

In [Olivo-Marin, 2002, Zhang et al., 2007], an iteratively undecimated wavelet transform (IUWT) [Mallat, 1989] of the image is exploited to detect objects of various sizes. A wavelet multiscale product (WMP) operation is performed in [Olivo-Marin, 2002], which consists, for every point, in multiplying the wavelet coefficients of different scales to reveal correlations across the scales. Indeed, from a given wavelet scale, spots respond more strongly to IUWT than uncorrelated noise. However, for low SNRs, noise has a higher response than spots at smallest scales, inducing wrong detections. Hence, smallest scales -up to a characteristic scale -must be discarded to lower the false detection rate. Finally, the WMP map is thresholded to get the binary detection map. The multiscale variance stabilizing transform (MS-VST) method relies on variance stabilization to rule out insignificant coefficients of the IUWT [START_REF] Zhang | Multiscale variance-stabilizing transform for mixed-Poisson-Gaussian processes and its applications in bioimaging[END_REF]. Then, the image is reconstructed without taking into account the coarsest scale, corresponding to the background structures, nor the smallest ones corresponding to noise. The spots are finally detected by thresholding the reconstructed image.

Therefore, with both IUWT-based methods, the set of wavelet scales must be chosen accordingly to the spot size.

Finally, h-dome (HD) methods [START_REF] Smal | A new detection scheme for multiple object tracking in fluorescence microscopy by joint probabilistic data association filtering[END_REF], Rezatofighi et al., 2012, Vincent, 1993] detect local maxima, called domes, in a LoG-or Gaussian-filtered image. The kernel must be chosen smaller than the spots. Peaks of the filtered image with an amplitude greater than a given height h (hence, the name of the method) are extracted. The so-built "dome map" comprises small domes corresponding to noise, domes corresponding to spots, and large domes corresponding to background structures. To discard irrelevant large and small domes, samples are generated according to the domes map seen as an importance sampling function. Domes containing too few samples are removed since they probably correspond to noise. Domes where samples are too scattered are also removed, because they probably correspond to large background structures. Thus, the maximum dome size must be carefully set. However, the objects to detect do not often have the same magnitude h, so that the method sometimes merges very bright neighboring spots, and sometimes misses dark spots. To tackle this problem, Rezatofighi et al. [START_REF] Rezatofighi | A new approach for spot detection in total internal reflection fluorescence microscopy[END_REF] proposed a method called maximum possible height dome (MPHD) for locally detecting the best height threshold h. Then, the norm of the spatial image gradient is thresholded, which is more robust to strong background variations than directly thresholding intensity.

Table 3.1 collects the main parameters of the aforementioned methods to be set by the user. All these methods have in common one or two critical scale parameters C������ 3. T�� ��������� �� ���� ��������� whose optimal value is closely related to the size of the objects to detect. If the detection scale is chosen too small, overdetection occurs due to noise. If it is set too high, objects are smoothed out or merged when close to each other. Fortunately, the scale parameter, whatever it is, can often be inferred from image data. We will address this issue in Chapter 4. Moreover, let us point out that our characteristic scale detection method is not limited to our detection framework. Indeed, it could be applied as a preprocessing step to most of the aforementioned methods involving a scale-related parameter.

Most methods end up delivering a binary detection map after thresholding a filtered or reconstructed image. However, when processing fluorescence microscopy sequences, the statistics of the image may vary in time due to photo-bleaching, so that one threshold should be set for each frame according to the current image intensity range. Obviously, this approach is not applicable to sequences containing hundreds of frames, or to a large datasets containing images of various dynamic range. In contrast, as described in Chapter 5, we propose a locally adapted threshold automatically inferred from local intensity statistics. The user on his/her side only fixes once for all a PFA value which can be used for all the images of the conducted experiment.

4

Automatic selection of the detection scale T �� ������ ����� selection step consists in determining the most represented scale in the input image I , that is the scale that most spots share. To this end, we rely on the framework developed by Lindeberg for scale-space analysis [Lindeberg, 1993]. In the remaining of the thesis, I will denote a 2-dimensional image, that is I :

Ω I ⊂ Z 2 → R,
where Ω I is the image domain. In this thesis, we have focused on 2D TIRFM images, but the proposed spot detection method can be straightforwardly extended to three dimensions.

Scale-space representations

Definition 1. The scale-space representation L of an image I is a 3-dimension map defined in [Lindeberg, 1993] as:

L : Ω I × S → R � p, s � � → � G s * I �� p � , (4.1)
where

p is a point in Ω I , � G s �
s ∈S is a family of convolution kernels of scale s , and the set S of scales is a subset of R * + .

Gaussian kernels are often proposed to build the scale-space representation. As a matter of fact, they are the best kernels for representing continuous signals [Lindeberg, 1993]. However, in the case of discrete signals such as digital images, sampled Gaussian kernels are not optimal, because they fulfill the scale-space representation for restrictive conditions on S only.

Proposition 1 (from [Lindeberg, 1990]). If

� G s � s ∈S is a family of sampled Gaussian ker- nels of variance s , L is a scale-space representation of I if and only if S = � s 0 r n � n ∈ N �
where s 0 is a strictly positive real number and r an odd integer.

Thus, the smallest possible ratio r between two consecutive scales is 3 when using sampled Gaussian kernels. We first adopted this approach in a preliminary work described in [Basset et al., 2014a] with S = � 1, 3, 9, 27, 81 � . However, since vesicles are usually under-resolved in TIRFM images, small scales (say, between 1 and 3) are of primary interest. Hence, we must be able to deal with scale sampling finer than in [Basset et al., 2014a] to further improve detection results. As a consequence, we prefer to employ a non-Gaussian kernel, called discrete analogous of the Gaussian kernel or discrete Gaussian in short. It will allow us to introduce scales of arbitrary precision. It is a separable filter, based on the family of modified Bessel functions of integer order

� B n � n ∈N 1 .
When exploiting this filter family, the scale-space representation holds with scales arbitrarily close to each other [Lindeberg, 1990]. Precisely, we will consider the scale set S = �

s 0 r n � n ∈ N � ,
where the scale ratio r is a strictly positive real number and s 0 = 1. In our TIRFM images, the pixel size (about 160 nm) is smaller than the optical resolution (about 200 nm), so that no objects of scale smaller than 1 can be found, thus the choice of s 0 = 1. The one-dimensional discrete Gaussian kernel G s of scale s is given by:

∀ n ∈ Z , G s (n) = e -s B n (s ) . (4.2)
The two-dimensional kernel is obtained by convolution, resulting in the following scale-space formulation:

∀ � p, s � ∈ Ω I × S , L � p, s � = G s * G s T * I � p � , (4.3)
where T denotes the transposition operation. This formula straightforwardly extends to higher dimension.

The Bessel functions have no closed form and B n is defined as the canonical solution of the following differential equation [START_REF] Abramowitz | Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables[END_REF]:

∀ z ∈ C , z 2 ∂ 2 B n ∂z 2 (z ) + z ∂B n ∂z (z ) + � z 2 -n 2 � B n (z ) = 0 . (4.4)
In practice, recursive approximations of the Bessel functions (such as the one proposed in [START_REF] Abramowitz | Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables[END_REF]) are precise enough to cope with our problem.

The discrete Gaussian filter is a low-pass filter which smooths out a certain amount of noise depending on the filter scale. In order to enhance the objects to detect, we apply the Laplacian operator ∆. For two-dimensional images, we use the most isotropic

Laplacian kernel [Lindeberg, 1990] defined as:

∆ = � � � � � � � 1/6 2/3 1/6 2/3 -10/3 2/3 1/6 2/3 1/6 � � � � � � � . (4.5)
To save computation time and owing to convolution commutativity, we apply the Laplacian operator once for all to the input image, before applying the scale-space operator. Moreover, the contrast of the LoG map decreases when scale increases. A normalization operation is performed to countervail this effect, which gives the following normalized discrete LoG filter H :

∀ � p, s � ∈ Ω I × S , H � p, s � = α(s ) ∆L � p, s � = α(s ) � G s * G s T * ∆I � � p � , (4.6)
where α is the normalization function. Lindeberg proposed to normalize the discrete LoG kernel with its L 1 -norm [Lindeberg, 1990], which gives:

α -1 (s ) = +∞ � n=-∞ � ∆G s (n) � . (4.7)
Due to the presence of absolute values, the normalization function (4.7) is not derivable at some particular points of the Bessel functions. Instead, we normalize with the L 2 -norm of the discrete Gaussian kernel.

Proposition 2. The proposed normalization function is given by:

α -1 (s ) = +∞ � n=-∞ � G s (n) � 2 = e -2s B 0 (2s ) . (4.8)
Proof. The discrete Fourier transform (DFT) of G is defined as [Lindeberg, 1990]:

F � G � (θ) = e s (cos θ-1) .
Using Parseval's theorem, we get:

α -1 (s ) = 1 2π � π -π � e s (cos θ-1) � 2 dθ = 1 2π � π -π e 2s (cos θ-1) dθ = 1 e 2s 1 2π � π -π
e 2s cos θ dθ .

Then, we conclude using the integral representation of B 0 [START_REF] Abramowitz | Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables[END_REF]:

B 0 (2s ) = 1 2π � π -π
e 2s cos θ dθ .

�

We will use this expression in the remainder of the thesis. The scale-space representation is illustrated in Figure 4.1 for a TIRFM image.

(a) Input image I (b) LoG L � 1 � at scale 1 (c) LoG L � 2 � at scale 2 (d) LoG L � 4 � at scale 4 (e) LoG L � 8 � at scale 8 (f) LoG L � 16 � at scale 16

Blob distribution and scale selection

Since the LoG filter is negative and spots are bright in the input image, we focus on the negative extreme values of H , which leads us to the concept of (negative) blob as defined in [Lindeberg, 1998].

Definition 2. A blob b = � p b , s b � ∈ Ω I × S is a local minimum of H : ∀ � p, s � ∈ ν p , H � p b , s b � � H � p, s � , (4.9)
where ν p is a 3

N +1 neighborhood of b in Ω I × S for an image of N dimensions.
The set of blobs (resp. blobs at scale s ) of an image I is denoted as

B � I � (resp. B s � I � ).
Proposition 3. The set of blobs B � I � of an image I is unchanged when adding a constant a 0 and/or multiplying I by a positive constant a 1 :

∀ a 0 ∈ R , ∀ a 1 ∈ R * + , B � a 1 I + a 0 � = B � I � . (4.10) Proof. Let b = � p b , s b � ∈ B � I � be a blob of I and take � p, s � ∈ ν p . Definition 2 gives: H � p b , s b � � H � p, s � .
Since a 1 > 0, and by linearity of the convolution operation,

H a 1 I +a 0 � p b , s b � � H a 1 I +a 0 � p, s � . Therefore, b is also a blob of a 1 I + a 0 , so B � I � ⊂ B � a 1 I + a 0 � . Analogously, we get B � a 1 I + a 0 � ⊂ B � I � and we conclude that B � a 1 I + a 0 � = B � I � .

�

As we showed in [Basset et al., 2014a], the blob detection itself is not sufficient to satisfyingly extract spots from noisy images. Instead, blob detection will be used to select the LoG scale. In order to select the optimal scale s � , likely to correspond to the average size of the spots of interest, we have to find the scale at which the blob number is the highest, while discarding blobs due to noise. Precisely, we take into account the number of blobs normalized by the image area, that is, the empirical blob density.

The empirical blob density in I at scale s is then given by:

∀ s ∈ S , ρ I (s ) = � � � B s � I � � � � |Ω I | , (4.11)
where |Ω I | is the area of I evaluated in square pixels (px 2 ).

If the image I was noise-free, the maximum of ρ I (s ) would be attained close to the most frequent spot scale. However, due to noise, some detected blobs do not correspond to real spots. As explained in Section 2.3.4, we suppose that the image is corrupted by an additive Gaussian noise (possibly, after noise variance stabilization).

Then, I decomposes into a noise-free image I 0 containing only structures (including spots) and an additive Gaussian noise component g :

I = I 0 + g . (4.12)
Should I exhibit no structure, the density of blobs detected at each scale s would be ρ g (s ). Therefore, we aim at selecting the scale s � at which the dissimilarity between ρ I (s ) and ρ g (s ) is the largest. We first propose four possible scale selection criteria while motivating their formulation. Then, we will compare them to retain one criterion.

First, we can take the difference or the ratio of the densities, which leads to the two following selection criteria, respectively:

C D : s � = arg max s ∈S � ρ I (s ) -ρ g (s ) � ; (4.13) C R : s � = arg max s ∈S ρ I (s ) ρ g (s ) . (4.14)
Since we are dealing with a white Gaussian noise, pixel intensities in the noise image g are assumed to be independent, and all the pixels have the same probability of being blobs at scale s , which precisely means that the number of blobs follows a binomial distribution [Kingman, 1992]. Then, if the number of blobs is high enough, the C������ 4. A�������� ��������� �� ��� ��������� ����� binomial distribution leads to a Poisson distribution (which is the limiting case of the binomial distribution [Kingman, 1992]) for the counting variable

� � � B s � g � � � � . Therefore,
we can also resort to the symmetric Pearson distance [START_REF] Pearson | On the criterion that a given system of deviations from the probable in the case of a correlated system of variables is such that it can be reasonably supposed to have arisen from random sampling[END_REF], Belongie et al., 2002] and Kullback-Leibler divergence [START_REF] Hannig | On Poisson signal estimation under Kullback-Leibler discrepancy and squared risk[END_REF] to evaluate the dissimilarity between the two densities. Two more selection criteria can then be formulated:

C P : s � = arg max s ∈S � ρ I (s ) -ρ g (s ) � 2 2 � ρ I (s ) + ρ g (s )
� ;

(4.15)

C KL : s � = arg max s ∈S � � ρ I (s ) -ρ g (s ) + ρ I (s ) log ρ g (s ) ρ I (s ) � � . (4.16) While exact for � � � B s � g � � � �
, the independence assumption is only an approximation for

� � � B s � I � � � �
, which is precisely the reason why the proposed scale selection criteria work:

we are looking for the scale at which the independence hypothesis does not hold, that is, spots lie at this scale.

To evaluate ρ g (s ), we prefer to avoid a time-consuming denoising step, which could also introduce artifacts. Instead, we simulate g by generating an image � g containing only white Gaussian noise.

According to Proposition 3, it is unnecessary to estimate the mean and variance of g , so that we merely sample � g from the standard normal distribution. Therefore, we can estimate once for all ρ g (s

) ≈ ρ � g (s ) with ρ � g (s ) = � � � B s � � g � � � � /|Ω � g |, instead
of generating one noise image per input image. To limit quantification noise in the blob counting, it is preferable to generate a large image � g . In practice, we use an image of size 4096 x 4096. In our experiments, using an even larger noise image did not impact the scale selection results.

As shown in Figure 4.2b, blob densities in I and � g are close to each other except in the vicinity of the spot characteristic scale s � . Thus, the optimal scale is in general correctly determined by all the criteria. We have studied the sensitivity of the method to the ratio r between consecutive tested scales. Results on a real TIRFM image are presented in Figure 4.2c-f for the four criteria. The sensitivity of the scale ratio parameter r is very low. Nevertheless, using an extremely small ratio value such as 1.1 can give bad results due to the limited global amount of blobs. To achieve a precise enough while reliable scale selection, we fix the scale ratio r to 1.2 in all the experiments and we take S = � 1, 1.2, 1.44, ..., 18.49, 22.19 � .

To evaluate the scale selection precision of the different criteria, we have generated synthetic images with the ImageJ plugin presented in [START_REF] Smal | Quantitative comparison of spot detection methods in fluorescence microscopy[END_REF]. They contain Gaussian spots of various sizes and are corrupted by a Poisson noise for different SNRs, as illustrated in Figure 4.3. Table 4.1 demonstrates the strong correlation between the size of the spots (defined by the Gaussian variance) and the selected scale. When SNR is high enough (greater than 3), the same scale is selected by all the s ∈ � proposed criteria. The selected scale is different however, when objects are small and noise level high. In this case, C D , C P and C KL tend to underestimate the characteristic scale, while C R tends to overestimate it. If a scale lower than the spot variance is selected, the image is less filtered and less irrelevant structures are smoothed out, then, more false detections arise. On the contrary, a scale slightly higher than the spot variance is selected by C R , so that a larger kernel further reduces noise, resulting in fewer false detections. As a consequence, for detection purpose, we will prefer C R for higher precision. However, if a post-processing step is affordable to eliminate erroneous detections as in [Basset et al., 2014b], C D and C KL may be preferable for higher recall, since less spots will be smoothed out. In the approach we proposed in [Basset et al., 2014a], the scale ratio between successive scales was constrained to be an odd integer (3 in practice). In [Basset et al., 2014a], we also set s 0 = 1 as the finest scale. Choosing another value for s 0 could have helped obtaining a scale closer to the true spot variance for some variance, but the overall accuracy of the scale selection would not be increased. Really getting an arbitrary accuracy, as with the proposed method ATLAS, would have required to first automatically select the optimal s 0 value, which depends on the processed image, thus being an open issue by itself. For example it would be impossible to get an accurate scale both for an image containing spots of variance 9 and for an image containing spots of variance 16, with the same s 0 . In contrast, as shown in Table 4.1, with our new scale selection method and still with s 0 = 1, the first image would be correctly processed at scale 8.92, and the second image at scale 15.41.

r n � n ∈ N � r = 1.1 r = 1.2 r = 1.44 r = 2.1 (c) Criterion C D 1 2 4 8 16 1 2 3 4 5 s ∈ � r n � n ∈ N � r = 1.1 r = 1.2 r = 1.44 r = 2.1 (d) Criterion C R 1 2 4 8 16 0 0.4 0.8 1.2 •10 -4 s ∈ � r n � n ∈ N � r = 1.1 r = 1.2 r = 1.44 r = 2.1 (e) Criterion C P 1 2 4 8 16 0 1 2 3 •10 -4 s ∈ � r n � n ∈ N � r = 1.1 r = 1.2 r = 1.44 r = 2.1 (f) Criterion C KL

Spot segmentation

O ��� ��� ������ ����� �� ����������, we can proceed to actual spot segmen- tation. Since the scale selection step relies on LoG, it is natural to detect vesicles based on this particular filter. Furthermore, it has been shown in [START_REF] Sage | Automatic tracking of individual fluorescence particles: Application to the study of chromosome dynamics[END_REF] that LoG is close to the optimal filter in applications like ours, that is detecting subresolved objects in fluorescence microscopy images.

Adaptive thresholding

As explained in Section 4.2, our goal is to extract the lowest values of the selected LoG

ψ � p � = Φ � � H � p, s � � -µ H σ H � � , (5.1)
where Φ is the Gaussian probability cumulative density function.

Equation equation (5.1) can be inverted to get a threshold value below which a point is detected, according to a user-selected probability of false alarm P FA , or pvalue:

τ H : Ω I → R p � → σ H Φ -1 � P FA � + µ H . (5.2)
Let us point out that we need to compute Φ -1 only once.

The local thresholding can thus automatically adapt to the local image statistics, while the PFA setting does not depend on the image intensity range. As a consequence, the spot detection is not affected by photobleaching when processing fluorescence microscopy image sequences. Indeed, the PFA is a parameter which is not directly related to the image properties but to the desired performance of the algorithm. Thus, it can be set once for all for a whole set of images in a given experiment.

In contrast, most aforementioned detection methods have at least one threshold parameter, which directly depends on the image characteristics and has to be manually set.

Estimation window

Generally, a square window of given radius r W centered at point p ∈ Ω I is used [START_REF] Kapur | A new method for gray-level picture thresholding using the entropy of the histogram[END_REF], Sahoo et al., 1988, Sahoo et al., 1997, Phansalkar et al., 2011, Sezgin and Sankur, 2004, Wilkinson and Schut, 1998]. However, while simple and fast, this type of window presents several drawbacks. First, the choice of the window size is critical due to the discontinuity at the window border; increasing or reducing the window radius by one pixel only can affect the estimated mean and variance. For the same reason, the resulting threshold map τ can be spoiled, exhibiting strong blocky artifacts as illustrated in Figure 5.2b. Also, the square window is not invariant to image rotation. To cope with these undesirable effects, we prefer to utilize a Gaussian weighted window W p with a bandwidth σ W :

∀ q ∈ Ω I , W p � q � = 1 � 2πσ 2 W exp � � � - � q -p � 2 2 2σ 2 W � � � . (5.3)
To speed up the statistics estimation, we rely on the recursive Deriche filter [Deriche, 1990]. The good properties of the Gaussian window are kept even if the window profile is not strictly Gaussian. As reported in the next section, the smooth weight decrease allows for a very low sensitivity of the window size parameter. Furthermore, giving more importance to central points naturally yields a better estimation of local statistics. for the square window.

The computation load with Gaussian windows will obviously be higher than with square patches, but not dramatically. Tests were performed on a laptop with 2.3 GHz Intel i7 processor. For 512 × 512 images, the computation time is 50 ms per image with a quasi-Gaussian window, and 15 ms with a square one relying on a moving average algorithm [START_REF] Kenney | Mathematics of Statistics[END_REF], Crow, 1984, Shafait et al., 2008]. For both windows, the numerical complexity is linear with the size of the image.

Finally, in an image sequence, frames are segmented independently, so that the program could even be sped up by processing the frames in parallel. Let us also stress that, if the image background is stationary, the estimation window could be extended to a spatiotemporal window in order to estimate the background statistics still more robustly.

Experimental results

We have compared ATLAS to state-of-the-art spot detection methods in a wide variety of cases. Comparative quantitative evaluation was carried out on several datasets with ground-truth. The first dataset is generated with the Synthetic Data Generator ImageJ plugin introduced in [START_REF] Smal | Quantitative comparison of spot detection methods in fluorescence microscopy[END_REF] and later used in [Basset et al., 2014a,Rezatofighi et al., 2012]. Twelve methods were compared on this dataset, which is (to our knowledge) the most complete comparison of spot detection methods to date, but the images remain somewhat too artificial and too simple.

As mentioned in the Introduction section, we have conducted comparative experiments on three other datasets involving more complex contents with the most competitive detection methods, namely MS-VST, MPHD, HD and conditional random fields for protein transport carriers segmentation (C-CRAFT). First, Boulanger et al. [START_REF] Boulanger | A simulation and estimation framework for intracellular dynamics and trafficking in video-microscopy and fluorescence imagery[END_REF] and Rezatofighi et al. [START_REF] Rezatofighi | A framework for generating realistic synthetic sequences of total internal reflection fluorescence microscopy images[END_REF] proposed particle dynamics simulators, referred in the sequel as Traffic simulator and TIRFM simulator, respectively. The Traffic simulator was used in [START_REF] Pécot | Background fluorescence estimation and vesicle segmentation in live cell imaging with conditional random fields[END_REF] to evaluate the performance of several methods. Secondly, we have constructed another image dataset named Spot in M10 where image backgrounds are extracted from real TIRFM images.

As stated in Section 2.3.4, all the images processed in this section are first stabilized using the generalized Anscombe transform (GAT)-based variance stabilization method described in [START_REF] Boulanger | Patch-based non-local functional for denoising fluorescence microscopy image sequences[END_REF].

Performance measures

ATLAS delivers a binary detection map. In order to evaluate the performance of the method and compare it to other ones, we compute the centroid of every segmented connected component, resulting in a set of locations {δ}. Then, following [START_REF] Smal | Quantitative comparison of spot detection methods in fluorescence microscopy[END_REF], an object ω of the ground-truth is correctly detected if and only if: (1) its nearest neighbor δ in the set of detected centroids is closer than 4 pixels away, and (2) ω is also the nearest neighbor of δ in the ground-truth set of locations. Let us denote N TP the number of true positives, N FP the number of false positives and N FN the number of false negatives.

We can evaluate different scores for every image and parameter setting. As in [Basset et al., 2014a, Rezatofighi et al., 2012, Smal et al., 2010] to compare methods in [Basset et al., 2014a, Rezatofighi et al., 2012, Smal et al., 2010].

Moreover, to compare ATLAS with the detection methods tested in [START_REF] Pécot | Background fluorescence estimation and vesicle segmentation in live cell imaging with conditional random fields[END_REF],

namely HD, MS-VST and C-CRAFT, we compute the precision Prec = N TP /(N TP + N FP )

and recall Rec = N TP /(N TP + N FN ). Varying the threshold parameter for the existing methods or the PFA value for ATLAS, we can plot the free receiver-operator characteristic (FROC), that is the TPR-versus-FPR* curve, and the precision-versus-recall curve.

That way, the behaviors of the methods can be evaluated more thoroughly. Additionally, we compute the area under the FROC curve as a performance score over a wide range of thresholds or PFA values. We also resort to the F-measure defined by the harmonic mean of precision and recall F = 2 Prec . Rec/(Prec + Rec), and more precisely to the best reachable F-measure F*.

Synthetic Data Generator

In [Basset et al., 2014a, Rezatofighi et al., 2012, Smal et al., 2010], twelve methods are evaluated over six image sets of 16 images each. They are depicted in Figure 5.3. Two object shapes are considered: isotropic Gaussian spots of standard deviation 2 pixels, and elliptic Gaussian spots of standard deviations 5 and 2 pixels along the two principal axes, respectively. Three types of background are generated: uniform intensity (type A), horizontal intensity gradient (type B), and large random structures (type C). rounded values, while 1 means that every spot is correctly detected. In the later case, the 4096 spots of the sequence are all recovered even for a FPR* value lower than 0.01.

In [START_REF] Rezatofighi | A new approach for spot detection in total internal reflection fluorescence microscopy[END_REF], Smal et al., 2010], in order to perform a fair comparison, all the methods were run with various parameter settings and, for each method, the setting yielding the best performance was selected, that is the highest TPR for FPR* = 0.01. Therefore, reported results reflect the best possible performance of the compared methods. As for ATLAS, we have fixed the size of the background estimation window σ W once for all to 15 pixels, that is all the sequences were segmented using the same window parameter. Notwithstanding, as demonstrated in Table 5.1, ATLAS performs better than any other method described in [Basset et al., 2014a, Rezatofighi et al., 2012, Smal et al., 2010] on every image subset of the benchmark dataset.

Finally, we have conducted a complementary experiment to demonstrate the ability of the method to detect two different scales in an image. While not relevant in our study, the ability to deal with images containing objects of different sizes could be appealing in some other applications. Therefore, a series of 16 synthetic images containing 2048 spots of variance 4 px 2 and 2048 spots of variance 16 px 2 was generated.

A Poisson noise was added to get a SNR of 2 (which converts to 6 dB). The C R criterion As for the supervised method logistic regression with Markov random field (LR-MRF) described in [START_REF] Ruusuvuori | Image segmentation using sparse logistic regression with spatial prior[END_REF], the authors reported better results than those of 11 unsupervised methods on a synthetic dataset of 20 images named Subcell [START_REF] Ruusuvuori | Evaluation of methods for detection of fluorescence labeled subcellular objects in microscope images[END_REF]. In addition, we have evaluated our method on the same synthetic dataset. We give the results in Table 5.2, which demonstrates that our method ATLAS also outperforms the 12 supervised method of [START_REF] Ruusuvuori | Image segmentation using sparse logistic regression with spatial prior[END_REF].

Spot in M10 dataset

To further challenge spot detection methods, we have generated another synthetic image dataset to carry out comparative objective evaluation. The so-called Spot in M10 dataset was constructed by mixing real background images and synthetic spots, as illustrated in Figure 5.5. To obtain realistic backgrounds, the rolling-ball extraction method [Sternberg, 1983] was first applied to one hundred real TIRFM images. Other background extraction methods could have been used as well. Actually, we do not need a very precise estimation of the background to construct the dataset. Furthermore, in order to ensure a fair comparative evaluation, it is preferable to use a background extraction method which is not part of any spot detection methods. Then, for each background image, 45 to 80 vesicles, modeled as Gaussian spots, were randomly added inside the cell region. Two spot variance values were tested: 1 px 2 and 1.44 px 2 .

Finally, the resulting images are corrupted by various Poisson-Gaussian noises to form images very similar to real ones.

We prefer to rely on the peak signal-to-noise ratio (PSNR) instead of SNR to evaluate the noise level in fluorescence microscopy images. Indeed, it is more straightforward and easier to compute PSNR in real images. After stabilizing the variance, the PSNR can be estimated as the difference between the minimum and maximum intensities, since the noise variance after stabilization is supposed to equal 1. The estimated PSNR of the real TIRFM images we are dealing with ranges from 23 to 31 dB. In order to evaluate the robustness to noise of the methods, we have generated even noisier images, down to PSNR = 18 dB. Note that PSNRs are higher than aforementioned SNRs since we are considering the brightest spot in the sequence. As shown in Figure 5.5, we end up with a set of six image collections, each containing first, except in one case, where MPHD is slightly better around FPR* = 0.01. Whatever the evaluation measure, the best overall performance is achieved by ATLAS.

The sensitivity of the parameters involved in MS-VST and ATLAS was low enough to keep the same setting for the whole dataset. For MS-VST, scale levels 2 and 3 were used. For ATLAS, the window radius σ W was set to 60 pixels. Moreover, as illustrated in Figure 5.7, this window size is not critical. To achieve the best performance with MPHD, we varied the LoG standard deviation from 0.83 pixel for highest PSNR to 1 pixel for lowest PSNR. Results were unchanged when varying the maximum dome radius (Table 3.1) from 5 to 20 pixels.

To further investigate the behavior of ATLAS when varying the parameters, we have conducted additional experiments. First, Figure 5.7 demonstrates the extremely low sensitivity of the estimation window size. FROC curves are indistinguishable when doubling σ W from 40 to 80 pixels. Secondly, we evaluated the precision of the scale selection as reported in Figure 5.8. In the M10-1-30 image subset, spots have a variance of 1 px 2 , so that one might think that best results would be achieved with a matched LoG filter of variance 1 px 2 . Yet, the scale selected with criterion C R is higher, namely 1.44 px 2 . In order to check the relevance of the selected scale, we run ATLAS with several LoG variances s taken in � 1, 1.2, 1.44, 1.73, 2 � . For the sake of visibility, we only plotted results with s ∈ � 1, 1.44, 2 � in Figure 5.8. The FROC curve obtained with the scale 1.44 px 2 corresponding to s � dominates others, which demonstrates the actual relevance of the selected scale and the efficiency of ATLAS. Besides, this shows the benefit of applying the LoG at the appropriate scale, and justifies the proposed method.

TIRFM simulator

In [START_REF] Rezatofighi | A framework for generating realistic synthetic sequences of total internal reflection fluorescence microscopy images[END_REF], the authors propose a method for simulating realistic TIRFM images containing spots undergoing a Brownian motion. The background of a real TIRFM image is first extracted using the h-dome method MPHD [START_REF] Rezatofighi | A new approach for spot detection in total internal reflection fluorescence microscopy[END_REF]. Then, objects are randomly spread inside the cell and a Poisson-Gaussian noise is added. We have evaluated the performance of our method on two sequences provided by the authors (Figure 5.9a-d), and compared it to the two best ranked methods in [START_REF] Rezatofighi | A new approach for spot detection in total internal reflection fluorescence microscopy[END_REF], Smal et al., 2010], namely MS-VST and MPHD. FROC curves are plotted in Figure 5.9e,f for both sequences. While negligible in the first sequence, photobleaching is strong in the second one, so that all the methods performed worse in the latter case. ATLAS obtained better results than MS-VST on both sequences, and than MPHD on the second sequence. However, MPHD performed significantly better than ATLAS and MS-VST on the first sequence for low FPR*. Let us nevertheless remind that the backgrounds were extracted using MPHD.

Traffic simulator

We can also compare our method to a robust patch-based method named C-CRAFT, which jointly estimates background and segments vesicles in 2D or 3D fluorescence image sequences, in the framework of conditional random fields [START_REF] Pécot | Background fluorescence estimation and vesicle segmentation in live cell imaging with conditional random fields[END_REF].

Vesicle segmentation and background estimation are formulated as a global energy minimization problem.

Even if spot detection is only one of the goals of C-CRAFT, the authors also compared their method to several spot detectors, including HD and MS-VST, on image sequences generated with the vesicle traffic simulation method [START_REF] Boulanger | A simulation and estimation framework for intracellular dynamics and trafficking in video-microscopy and fluorescence imagery[END_REF].

Therefore, we can compare the spot detection results of ATLAS on the dataset used in [START_REF] Pécot | Background fluorescence estimation and vesicle segmentation in live cell imaging with conditional random fields[END_REF]. It is composed of four sequences of 120 frames. Two real backgrounds were extracted for Circle-and Crossbow-shaped micropatterned cells. Vesicle motions are then simulated, and both sequences are corrupted with two different Poisson-Gaussian noises. Since no FROC curves were given in [START_REF] Pécot | Background fluorescence estimation and vesicle segmentation in live cell imaging with conditional random fields[END_REF], precision-recall curves of ATLAS are plotted in Figure 5.10 and compared to the results provided in [START_REF] Pécot | Background fluorescence estimation and vesicle segmentation in live cell imaging with conditional random fields[END_REF] for C-CRAFT, HD and MS-VST.

On the two Crossbow sequences, ATLAS outperforms the three methods. On the Circle sequence, which exhibits extremely contrasted background, ATLAS and C-CRAFT are on par for the highest PSNR (33 dB), but C-CRAFT performs better when noise level is high and ATLAS ranks second. This can be explained by the temporal integration used by C-CRAFT to jointly estimate background and detect vesicles.

Real images

Finally, we have conducted experiments on real images. As in most experiments on microscope images, it is not possible to extract the full ground-truth, but only to get an expert annotation, which is not perfect by nature. As a consequence, a quantitative comparison would not be completely reliable. Therefore, qualitative comments are given to analyze results on real data.

We have processed two real sequences of M10 cells expressing various fluorescently tagged proteins, namely Langerin and TfR tagged with pHluorin, and Rab11 tagged with mCherry.

The first 300-frame sequence was acquired with the Dual View optical beamsplitting technique [START_REF] Gidon | A Rab11A/myosin Vb/Rab11-FIP2 complex frames two late recycling steps of langerin from the ERC to the plasma membrane[END_REF]. As explained in Chapter 2, with this technique, several fluorescence channels are simultaneously acquired side-by-side on the sensor, resulting in images including two parts with very different intensity ranges, as shown in Figure 5.11a. The upper half of the image shows Rab11 GTPase. Overall, it is much brighter than the lowest half showing Langerin. Even in that complicated case, the adaptive thresholding approach produces very satisfactory results in both parts of the image. In the lower part, spots even darker than the upper part cell background are detected, while no false detection appears even in the brightest region of the upper part.

As stated in the introduction, various images fulfill the characteristic scale assumption. To show the ability of the method to cope with a very different type of images, ATLAS was also successfully applied to tissue micro-array (TMA) images, where the cores must be first detected, before being individually analyzed. In these images, the difficulty mainly comes from the high variability of the cores intensity. Indeed, as illustrated in Figure 5.12, spurious bright dots have a much higher intensity than most cores. To smooth out these dots while still detecting low-intensity cores, the LoG scale must be precisely set. Cores are correctly detected by ATLAS, and only one false positive arises corresponding to a broken but intense core. However, this kind of error can be easily discarded by post-processing [START_REF] Nguyen | Automatic core segmentation and registration for fast tissue microarray de-arraying[END_REF].

We have also compared the binary segmentation maps supplied by ATLAS and MS-VST on a second TIRFM image sequence of 300 images. We only compared ATLAS

Discussion

We have proposed a novel and efficient vesicle segmentation method called ATLAS which involves an automatic scale selection and a local threshold setting. It is dedicated to situations where most of the visible structures share about the same size in the image. The selected scale can be of any arbitrary precision. After determining the optimal scale, a LoG operator is applied on the images. The segmentation threshold is automatically and locally set according to a given PFA value. Overall, ATLAS outperforms state-of-the-art methods on various datasets, including a new one we have constructed and made publicly available for further comparison. Satisfactory segmentation results on several challenging real TIRFM images have been reported. We have shown that ATLAS is not sensitive to the Gaussian window size in the segmentation step. Moreover, the PFA value is a user-friendly parameter which allows the user to adapt the method to the targeted detection sensitivity according to the application needs and the further exploitation of the detection results. Thus, no specific knowledge is required on the algorithm itself, that is, the method can be used as a black box by someone non-expert in image processing. We have shown that ATLAS can be successfully applied to different kinds of images. We have also demonstrated that ATLAS can deal with a couple of scales if needed. We will further investigate the detection and exploitation of a wider set of scales if one or two scales are not sufficient to accurately describe the structures of interest. We also plan to apply ATLAS to three-dimensional images.

II Individual dynamics

Modeling of the vesicle fusion dynamics R ���������� dynamic protein behaviors in live cell fluorescence microscopy is of paramount importance to understand cell mechanisms. In the case of membrane traffic, cargo molecules are transferred from a donor to an acceptor compartment [START_REF] Prydz | Arrivals and departures at the plasma membrane: Direct and indirect transport routes[END_REF]. During the exocytosis process, a vesicle conveys cargo molecules to the plasma membrane, and then opens to expel them from the cell. At each step, dedicated molecular platforms are acting to transport, steer and deliver selected proteins. In microscopy imaging, this sequence of processes leads to a series of different dynamics which need to be untangled in order to understand the spatiotemporal coordination of the molecular actors. Total internal reflection fluorescence microscopy (TIRFM) is particularly well suited for focusing on the late steps of exocytosis, which occur at the plasma membrane [Brown, 2006]. However, even with this modality, it is still a challenging task to classify dynamics of vesicular movements and protein diffusion.

In this part, we investigate the dynamics of three proteins associated with vesicles during the exocytosis process of vesicle recycling: Rab11, Transferrin receptor (TfR)

and Langerin.

Rab11 is a so-called GTPase protein. Among others, it is required during the vesicle transport. However, Rab11 is also known to be associated with many other actors (see, e.g., [Schwartz et al., 2007, Hutagalung and[START_REF] Hutagalung | [END_REF]), and one objective is to decipher its dynamics in order to better understand the recycling process.

TfR and Langerin are transmembrane proteins, that is, they are inserted in biological membranes, including the plasma membrane, where they are involved in several biological processes. After they have completed a process in the cell, TfR and Langerin are recycled by a recycling endosome [START_REF] Uzan-Gafsou | Rab11A controls the biogenesis of Birbeck granules by regulating Langerin recycling and stability[END_REF], which means that they are sent to another place, where they have to carry out some function. If sent to the plasma membrane, they are transported by a recycling vesicle, which eventu-ally fuses to the plasma membrane. At this point, the transmembrane proteins diffuse in the plasma membrane and it can be reused for the following endocytosis.

To the best of our knowledge, apart from this thesis, the issue of diffusion after vesicle fusion was only addressed in [START_REF] Mele | Automatic identification of fusion events in TIRF microscopy image sequences[END_REF], Burchfield et al., 2010] with a simple model and estimation method. In these papers, very restrictive hypotheses are assumed, yielding non-realistic results. Moreover, the estimation method, which we recall later, does not exploit all the available information and lacks accuracy. Therefore, we propose to go further both in terms of modeling and estimation of the vesicle fusion dynamics.

In this chapter, we present the biological processes taking place during exocytosis in Section 6.1, before and after the vesicle fusion to plasma membrane. After presenting existing diffusion models in biology in Section 6.2, dynamical models for vesicle fusion dynamics are proposed in Section 6.4 for Rab11, TfR and Langerin. Strengths and limitations of the proposed models are discussed in Section 6.5. Fusion dynamics will then be detected and estimated in Chapter 7 and Chapter 8, respectively.

6.1 Exocytosis process

Sequence until vesicle fusion

Let us first present very briefly the successive steps constituting the process, which are represented in Figure 6.1. The selected material being secreted is initially contained in a donor organelle. In order to leave the latter, a vesicle forms by invagination of the organelle membrane (Figure 6.1b). The material is enclosed when the nascent vesicle closes and detaches from the donor membrane.

As illustrated in Figure 6.1c, a vesicular motor then pulls the vesicle along the cytoskeleton, which allows the vesicle and its contents to move toward the plasma membrane. Let us mention that the Rab11 protein is bound to the vesicle and most probably to a molecular motor during the transport step, and stays with the vesicle for a while after transport, as explained in Section 6.1.2.

The overall process follows a docking-tethering-fusion sequence. While the vesicle approaches the plasma membrane, specific receptors in the latter, target SNAREs or t-SNAREs in short, bind with the vesicle own receptors, v-SNAREs, constituting a SNARE complex represented in Figure 6.1d. Then, the complex pulls the vesicle toward the plasma membrane until fusion.

Fusion consists in melting the vesicle and plasma membranes. Specifically, after docking at the plasma membrane, the vesicle opens while its membrane becomes part of the plasma membrane. As depicted in Figure 6.1e, the interior of the vesicle therefore becomes the exterior of the plasma membrane, and the material which was ). The behavior of Rab11 is not easy to establish: it could be that (i)

Rab11 diffuses along the plasma membrane and later dissociates from the membrane to the cytosol; or (ii) it dissociates and diffuses in the cytosol immediately after fusion.

We will focus on this final step in the remaining of Part II.

As mentioned above, the behavior of TfR and Langerin on the one hand, and of Rab11 on the other hand, are very different. Specifically, TfR and Langerin are two transmembrane proteins spanning across the vesicle membrane, while Rab11 binds at the outer side the vesicle membrane, as depicted in Figure 6.1. We shall now describe the behavior and specificity of those proteins.

Diffusion of Rab11 after vesicle fusion

The exact sequence of events followed by Rab11 and the associated dynamics are not well established. However, the interaction of Rab11 with many proteins during exocytosis has been demonstrated [START_REF] Lindsay | Rab11-FIP2 functions in Transferrin recycling and associates with endosomal membranes via its COOH-terminal domain[END_REF],Hales et al., 2001,Lapierre and Goldenring, 2005, Gidon et al., 2012, Boulanger et al., 2014, Novick et al., 2006],

suggesting that Rab11 is required for the successive steps of exocytosis to be completed from transport to fusion, and shows that Rab11 lies in the vicinity of the vesicle at least until fusion.

As opposed to transmembrane proteins, Rab11 can leave the membrane for the cytosol, an event named dissociation. While it is known that Rab11 dissociates from the plasma membrane, at some time point of the exocytosis or later, it is unclear if

Rab11 dissociates before, during, or after vesicle fusion. For example, Rab11 may first move along the plasma membrane before leaving to the cytosol, that is, diffusing in two dimensions before dissociating and diffusing in three dimensions, or on the contrary, it may only diffuse in three dimensions. In both cases, deciphering the sequence of consecutive steps during and after vesicle fusion remains an open question.

Owing to thermal agitation in the membrane or cytosol, Rab11 proteins are assumed to undergo a Brownian motion [Einstein, 1905]. This model holds if the plasma membrane and cytosol are isotropic homogeneous, and free of obstacles. Moreover the membrane is supposed to be flat in the thesis. However the hypotheses must only be valid in a small region around the vesicle fusion location, as we will rely on very local estimation methods.

Furthermore, we focus on image sequences depicting micropatterned cells, which means that they are "glued" on the coverslip. This ensures that the cell does not drift during the acquisition and the membrane itself is supposed to be immobile.

Proteins undergoing a Brownian motion may go away from the location where vesicle fusion occurred, and we will see that the group motion of the proteins is a lateral diffusion process, that is, a process by which the protein concentration tends to homogenize in the membrane or cytosol [Einstein, 1956].

To select a dynamical model for Rab11, different candidates must be proposed.

We are principally interested in the dissociation process, i.e., we aim at distinguishing membrane and cytosol dynamics. For example, the following models are plausible, but others can be considered:

• When the vesicle fuses to the plasma membrane, Rab11 dissociates to diffuse in the cytosol;

© Adapted from an illustration by Mariana Ruiz Villarreal, Wikimedia Commons In the first case, we will see that a simple image model can be inferred, with only one biophysical parameter: the diffusion coefficient. In the other case, however, two concentration components superimpose in the image: (i) a component diffusing in two dimensions in the membrane; and (ii) a component diffusing in three dimensions in the cytosol. In addition, a flow model should be proposed to account for dissociation. We end up with two diffusion coefficients plus one or several parameters for the dissociation model.

Diffusion of TfR and Langerin after vesicle fusion

Transmembrane proteins span across the membrane, which they never leave. Therefore, their motion is constrained by the two-dimensional surface formed by the membrane they belong to. More precisely, transmembrane proteins move in the membrane by sliding between phospholipids, as illustrated in Figure 6.2. In the frame of exocytosis, TfR and Langerin either belong to the vesicle membrane or to the plasma membrane, but never dissociate to the cytosol.

Like for Rab11, the type of dynamics undergone by transmembrane proteins in the plasma membrane is assumed to be a Brownian motion [START_REF] Clegg | Translational diffusion of proteins and lipids in artificial lipid bilayer membranes: A comparison of experiment with theory[END_REF], which holds under the aforementioned hypotheses (isotropic, homogeneous, obstacle-free membrane). More complex models have been proposed to take into account interactions between proteins and membrane structures [START_REF] Sako | Compartmentalized structure of the plasma membrane for receptor movements as revealed by a nanometerlevel motion analysis[END_REF], Schuss et al., 2007,Rehfeldt and Stichlmair, 2007] or to allow for rotation [Cherry, 1979,Swaminathan et al., 1997]. However, the rotational component has been shown to be negligible with respect to the lateral component [START_REF] Almeida | Lateral diffusion in membranes[END_REF]. For an extensive survey on the diffusion in plasma membrane, we refer the reader to [START_REF] Clegg | [END_REF]Vaz, 1985, Almeida andVaz, 1995].

Existing diffusion models

For thirty years, numerous biological studies focused on membrane diffusion [START_REF] Vaz | Lateral diffusion of lipids and proteins in bilayer membranes[END_REF], Clegg and Vaz, 1985, Jacobson et al., 1987, Almeida and Vaz, 1995, Chen et al., 2006, Guo et al., 2008, Ramadurai et al., 2009] Generally, the models were proposed to analyze 2D diffusion in the plasma membrane, as considered in what follows. The extension to three dimensions is generally straightforward. If this is not the case, the three-dimensional model will be explicitly given.

Intensity model break-down

As represented in Figure 6.3, the mathematical model of I � p, t � for each point p ∈ Ω I and time t ∈ R + is fully determined by three components:

• The source particle distribution;

• The evolution model;

• The observation model.

The source distribution characterizes the way particles are spread over space before they start diffusing. The source distribution defines both the spatial distribution of the particles before they start diffusing, and the law governing their release time to the plasma membrane or cytosol. This will be thoroughly explored in the remaining of the chapter.

The particle evolution model is the mathematical description of the motion of the proteins after fusion. As explained above, it is assumed to be Brownian. Then, resulting from Brownian motion of each individual particle, lateral diffusion is the dynamical model governing the evolution of the whole particle population.

As for the observation model, we have explained in Chapter 1 that it is itself subdivided into components, including different noises and the optical transfer function or point spread function (PSF). In this chapter, we will first consider a noise-free observation model to derive the intensity model.

Since the evolution and observation models are fixed once for all, we investigate in this chapter the source distribution modeling, that is the distribution of the lightemitting particles before they start diffusing.

From Brownian motion to lateral diffusion 1

As a fundamental stochastic process, involved for instance in fluid dynamics [Fick, 1855, Maxwell, 1867, Philibert, 2005], heat transfer [Fourier, 1822, Carslaw andEger, 1959] or stellar dynamics [Reid andBrunthaler, 2004, Merritt, 2013], Brownian motion was extensively studied by mathematicians and physicists. While named after botanist Robert Brown, who observed such motion without a mathematical description [Brown, 1828], most mathematics of Brownian motion were formulated by Albert

Einstein in the second of the Annus mirabilis papers [Einstein, 1905]. More recently, biophysicists and biologists introduced Brownian motion to model various dynamic structures inside the cell [Hellriegel andGratton, 2009,Schuss, 2012], most extensively focusing on lipid and protein motion in the plasma membrane [START_REF] Vaz | Lateral diffusion of lipids and proteins in bilayer membranes[END_REF], Jacobson et al., 1987, Almeida and Vaz, 1995, Saxton and Jacobson, 1997, Chen et al., 2006, Guo et al., 2008, Ramadurai et al., 2009].

Mathematically, a given particle i undergoes a Brownian motion if and only if its position over time k i (t ) follows a Wiener stochastic process [Wiener, 1966]. This characterization shows that the particle displacement between arbitrary times t 1 and t 2 follows a normal distribution of mean 0 and variance 2D �

t 2 -t 1 � : ∀ t 1 , t 2 ∈ R + , k i (t 2 ) -k i (t 1 ) ∼ N � 0, 2D � t 2 -t 1 � I � , (6.1)
where I denoted the 2 × 2 or 3 × 3 identity matrix depending on the diffusion space dimension.

For the problem to be well-posed, some constraints must be added. To this end, we assume that particle i is static and located at point k 0i from t = 0 to a given release C������ 6. M������� �� ��� ������� ������ �������� time t 0i . At t 0i , particle i starts moving following a Wiener process:

∀ t � t 0i , k i (t ) = k 0i , (6.2) ∀ t > t 0i , k i (t ) ∼ N � k 0i , 2D � t -t 0i � I � . (6.3)

Local concentration

While Brownian motion describes the dynamics of individual particles, the dynamical model of a population of such particles is named lateral diffusion [Einstein, 1956].

To cross the bridge from Brownian motion to lateral diffusion, let us now introduce the concept of local concentration. In the vesicle, and later in the cytosol or plasma membrane, proteins of interest are numerous, so that, in TIRFM image sequences, we do not observe a single particle, but a population of N particles. Concentration is generally defined as the number of particles in a given region. Let ν p be a region around any point p. We define the local concentration C � p, t

� by reducing the region area � � � ν p � � � to zero, that is, reducing the region ν p to point p:

C � p, t � = lim � � � ν p � � � →0 � 1�i �N 1 ν p � k i (t ) � � � � ν p � � � (6.4)
where 1 is the indicator function.

If the particles are numerous enough, independent and identically distributed, the law of large numbers yields:

C � p, t � ≈ N f K (t ) � p � , (6.5)
where K (t ) is the random variable associated to the k i (t )'s, and f K (t ) denotes the probability density function (PDF) of K (t ). f K (t ) can be expressed as a function of the conditional PDF f K (t )|K 0 =p 0 , termed transition PDF, and the source PDF f K 0 :

f K (t ) � p � = � R + � R 2 f K (t )|K 0 =p 0 � p � f K 0 � p 0 � f T 0 � t 0 � dp 0 dt 0 , (6.6)
where the Wiener characterization equation ( 6.3) shows that f K (t )|K 0 =p 0 is the Gaussian PDF of mean p 0 and variance 2D � tt 0 � I. This transition PDF is specific to Brownian motion. Therefore, the whole dynamical model is only dependent on the source distribution f K 0 .

Source and diffusive concentrations

The total concentration C results from the sum of the source concentration C s and the diffusive concentration C d :

C = C s + C d , (6.7)
where, by definition, the particles considered for C s do not move, while those considered for C d undergo a Brownian motion. Specifically, C s is defined as:

C s (p, t ) = lim � � � ν p � � � →0 � 1�i �N 1 ν p � k 0i (t ) � 1 [0,T0i ] (t ) � � � ν p � � � (6.8)
and the law of large numbers gives:

C s (p, t ) ≈ N f K (t ) � p � � 1 -F T 0 (t ) � , (6.9)
where F T 0 is the cumulative density function (CDF) of the release time.

On the other hand, the diffusive concentration C d is defined as:

C d (p, t ) ≈ N f K (t ) � p � F T 0 (t ) .
(6.10)

Diffusion equation

Now that local concentration has been introduced, let us present the lateral diffusion model which governs evolution over time of the diffusive concentration.

If we assume that t 0i is identical for all the particles, we merely get, with t 0 = t 0i :

C � p, t � =            C s (p, t ) for 0 � t � t 0 , C d (p, t ) for t > t 0 .
(6.11)

In this case, the Fick's second law gives the evolution over time and space of the local concentration as a function of the diffusion coefficient D [Fick, 1855]: (6.12) where ∆ = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 (∂ 2 /∂x 2 + ∂ 2 /∂y 2 + ∂ 2 /∂z 2 in three dimensions) denotes the Laplace operator. Based on the law of large numbers, equation (6.12) is the macroscopic counterpart of the single particle description of Brownian motion (6.1) in the sense that it characterizes the evolution of a population of such particles observed at a higher scale.

∂C d ∂t = D ∆C d ,
On the other hand, in the general case where t 0i is different for all particles, as in the proposed "small-extent source with exponential decay release" (SSED) model, a

flow is introduced from C s to C d : ∂C d ∂t = D ∆C d - ∂C s ∂t . (6.13)
To our knowledge, we are introducing this flow in the study of membrane dynamics.

Fick's second law (6.12) can be solved by Fourier analysis, which yields the following closed form Green's function Φ [Fourier, 1822]:

∀ t > t 0 , ∀ p ∈ R 2 , Φ � p, t � = 1 4πD � t -t 0 � exp � � � - � p -p 0 � 2 2 4D � t -t 0 � � � � . (6.14)
Equation (6.14) can also be interpreted from the stochastic point of view as reflecting the probability of finding particles at position p and time t , if they undergo a Brownian motion of diffusion coefficient D and are initially concentrated at p 0 . Actually, instead of relying on the Fick's second law equation (6.12), the Green's function can also be viewed as a transition PDF:

∀ t > t 0 , ∀ p ∈ R 2 , f K (t )|K 0 =p 0 � p � = 1 4πD � t -t 0 � exp � � � - � p -p 0 � 2 2 4D � t -t 0 � � � � . (6.15)
This shows the equivalence of the microscopic, stochastic process and the macroscopic, deterministic partial differential equation (PDE).

Stationary model

In the stationary model [Sergeev, 2004], particles are initially spread with a uniform distribution, and start moving at t 0 :

∀ t ∈ R + , ∀ p ∈ R 2 , C s (p, t ) =            C 0 for 0 � t � t 0 0 for t > t 0 (6.16) and ∀ p ∈ R 2 , C d (p, t 0 ) = C 0 .
(6.17)

The Fick's second law then gives the temporal evolution of C d : (6.18) which yields:

∀ t > t 0 , ∀ p ∈ R 2 , C d (p, t ) = C 0 ,
∀ t ∈ R + , ∀ p ∈ R 2 , C � p, t � = C 0 . (6.19)
In turn, the stationary model is the multi-particle extension of the Brownian motion model, where particles themselves undergo a Brownian motion, but the local concentration C is expected not to vary in space and time, thus the stationary moniker.

When the particle number N is high, and the law of large numbers holds, the intensity I is thus expected to be constant, as illustrated in Figure 6.4a.

In contrast, when the particle density is low, intensity is not constant. Indeed, in this case, the microscope PSF gives the particles a Gaussian-like profile in the image, resulting in a sum of Gaussian spots, whose center positions are sampled from the uniform distribution (Figure 6.4b). Parameters of the model are C 0 , the diffusion coefficient D and the radius σ PSF of the PSF.

As presented in Chapter 8, many papers dealt with the problem of estimating the diffusion coefficient D under the stationary model, but it can only be used for low enough particle densities; otherwise, the intensity model does not even depend on D .

-2 -1 0 1 2 0 0.5

1 x C (x, t ) -2 -1 0 1 2 0 0.5 1 x I (x, t )
(a) Uniform concentration (left) and resulting intensity (right)

-2 -1 0 1 2 0 0.5 1 x C (x, t ) t 0 t 0 + 1 t 0 + 2 -2 -1 0 1 2 0 0.5 1 x I (x, t ) t 0 t 0 + 1 t 0 + 2 (b)

Point source model

In the point source model, we assume that all the particles are initially concentrated at p 0 and that they all start diffusing at time t 0 , that is, the source distribution is proportional to a spatiotemporal Dirac distribution:

∀ t ∈ R + , ∀ p ∈ R 2 , C s (p, t ) = C 0 δ � p -p 0 � δ � t -t 0 � . (6.20)
The Green's function (6.14) is the solution of the Fick's second law for a space-time Dirac source distribution. Then, by linearity of the Fick's second law, the concentration dynamics is merely obtained by multiplying the Green's function by C 0 :

∀ t > t 0 , ∀ p ∈ R 2 , C � p, t � = C 0 4πD � t -t 0 � exp � � � - � p -p 0 � 2 2 4D � t -t 0 � � � � . (6.21)
Then, to derive the intensity model, we need to incorporate the observation model, reduced to the PSF and gain of the microscope. In two dimensions, that is, for membrane diffusion, the PSF is modeled as a two-dimensional Gaussian function of variance σ 2 PSF [START_REF] Small | Fluorophore localization algorithms for super-resolution microscopy[END_REF], so I is obtained by convolving the con-C������ 6. M������� �� ��� ������� ������ �������� centration (6.14) with a Gaussian kernel of variance σ 2 PSF , which gives:

-6 -4 -2 0 2 4 6 0 0.5 1 x C (x, t ) t 0 t 0 + 1 t 0 + 2 -6 -4 -2 0 2 4 6 0 0.5 1 x I (x, t ) t 0 t 0 + 1 t 0 + 2
∀ t � t 0 , ∀ p ∈ R 2 , I � p, t � ∝ C 0 4πD � t -t 0 � + 2πσ 2 PSF exp � � � - � p -p 0 � 2 2 4D � t -t 0 � + 2σ 2 PSF � � � . (6.22)
For the sake of simplicity, we introduce constant A 0 such that:

∀ t � t 0 , ∀ p ∈ R 2 , I � p, t � = A 0 2D � t -t 0 � + σ 2 PSF exp � � � - � p -p 0 � 2 2 4D � t -t 0 � + 2σ 2 PSF � � � . (6.23)
The source distribution and expected concentration of the point source model are plotted in Figure 6.5.

Since it is very practical for computation and estimation purposes, the point source model was used in [START_REF] Mele | Automatic identification of fusion events in TIRF microscopy image sequences[END_REF],Burchfield et al., 2010] to estimate the dynamics of vesicle fusion. To this end, two supplementary hypotheses were assumed:

H1 The vesicle is pointwise;

H2 The whole transmembrane protein material is released instantaneously at t 0 . However, we will see that these hypotheses are not always justified, which will lead us to relax them and propose a more elaborated model.

To our knowledge, the extension of [START_REF] Mele | Automatic identification of fusion events in TIRF microscopy image sequences[END_REF], Burchfield et al., 2010] to the 3D diffusion has not been done yet, but it will be proposed in Section 6.3.3 as an ingredient of the Rab11 dynamics modeling.

FRAP model 2

In FRAP experiments, transmembrane proteins are supposed to be uniformly distributed with an average concentration denoted C 0 . At time t 0 , a leaser beam locally kills fluorescence by photobleaching [START_REF] Axelrod | Mobility measurement by analysis of fluorescence photobleaching recovery kinetics[END_REF], Cherry, 1979, Kapitza and Jacobson, 1986, Seiffert and Oppermann, 2005]. Bleached transmembrane proteins then stop emitting photons, so that they do not appear anymore in C or I .

2 FRAP theory and tools were presented to us by François Waharte at UMR 144.

-6 -4 -2 0 2 4 6 0 0.5

1 x C (x, t ) t 0 t 0 + 1 t 0 + 2 -6 -4 -2 0 2 4 6 0 0.5 1 x I (x, t ) t 0 t 0 + 1 t 0 + 2 Figure 6.6 -FRAP model at three time steps for C 0 = 1, D = 1, p beam = 0, σ beam = 2
To ease computation, the beam profile is often approximated by a Gaussian function of center p beam and radius σ beam [START_REF] Seiffert | Systematic evaluation of FRAP experiments performed in a confocal laser scanning microscope[END_REF]. After t 0 , fluorescent (resp. bleached) particles, which undergo a Brownian motion, repopulate the bleached (resp. fluorescent) region, until both bleached and fluorescent particles are uniformly spread. Then, the intensity becomes uniform again. Some extended models were proposed, e.g., to account for uniformly distributed static particles. We refer the reader to [START_REF] Carrero | Using FRAP and mathematical modeling to determine the in vivo kinetics of nuclear proteins[END_REF] for more details on extended FRAP models.

Owing to the superposition principle, the solution of the Fick's second law (6.12) for an arbitrary source profile s : R 2 → R + is obtained by convolving the Green's function with the profile. In FRAP experiments, the laser beam is generally assumed to have a Gaussian profile, resulting in a Gaussian "hole" fluorescence profile presented in Figure 6.6:

s � p � = 1 2πσ beam - 1 2πσ beam exp � � � - � p -p beam � 2 2 2σ 2 beam � � � . (6.24)
The intensity model is straightforwardly deduced by convolution with the PSF, as the convolution of two Gaussian functions is the Gaussian function whose variance is the sum of the original variances:

I � p, t � = A 0 σ 2 PSF + σ 2 beam - A 0 2D � t -t 0 � + σ 2 PSF + σ 2 beam exp � � � - � p -p 0 � 2 2 4D � t -t 0 � + 2σ 2 PSF + 2σ 2 beam � � � .
I � p, t � = exp � - t τ � ×         A 0 σ 2 PSF + σ 2 beam - A 0 2D � t -t 0 � + σ 2 PSF + σ 2 beam exp � � � - � p -p 0 � 2 2 4D � t -t 0 � + 2σ 2 PSF + 2σ 2 beam � � �         .
(6.26)

Proposed fusion models for Rab11

While the sequence of events of Rab11 is not totally established yet, we propose a meta-model to encompass different hypothetical models. It is illustrated in Figure 6.7.

Numerous models compatible with current biological knowledge can be instantiated from it by defining flows between a stationary state (which accounts for the proteins remaining in the vesicle after fusion) and three dynamical states, namely translation (from transport to docking), membrane diffusion and cytosol diffusion. We consider the following flows:

• Fusion specifies the beginning of the diffusion(s);

• Release characterizes the flow between static concentration and concentration diffusing in the plasma membrane;

• Dissociation is the process by which Rab11 leaves the membrane, regardless of its previous state; as explained in [START_REF] Carrero | Using FRAP and mathematical modeling to determine the in vivo kinetics of nuclear proteins[END_REF], Michelman-Ribeiro et al., 2009,Im et al., 2013], it is modeled as an exponential decay of the observed concentration.

Let us detail the dynamical models for translation and the two diffusions.

Vesicle motion before fusion

We assume that the frame rate (10 frames/s in the processed sequence) is sufficient to neglect rotation and deformation of the vesicle. Hence, the vesicle is supposed to undergo a translational motion between two successive time points t and t + 1. LetV � i (t ) denote the spatial support at time t of such a vesicle. By definition of the translation model, all the pixels of the connected component V � i (t ) formed by the vesicle at time t share the same displacement w i :

∀ p ∈ V � i (t ) , w i � p, t � = w i (t ) , (6.27)
where

w i � p, t
� denotes the displacement of p ∈ V � i (t ) between t and t + 1.

Two-dimensional diffusion along the plasma membrane

In [START_REF] Mele | Automatic identification of fusion events in TIRF microscopy image sequences[END_REF], the point source model was used to model the diffusion in the plasma membrane after vesicle fusion. In contrast, in order to take into account the size of the vesicle, we consider a finite spatial support (e.g., a spherical vesicle) of diameter 2r , that is, the spatial source profile s : R 2 → R is null outside a disc B r ⊂ R 2 of radius r . Without loss of generality, let us assume that � R 2 s = 1, and that B r is centered in (0, 0). Let I r (resp. I δ ) be the solution of the Fick's second law for the source of profile s (resp. for the point source case).

Proposition 4. The difference between I r and I δ is small when r is small. More precisely:

� � � I r -I δ � � � = O � r σ PSF � , (6.28)
by introducing the Landau O notation.

Proof. Let us up-bound the absolute difference between I r and I δ :

I r � p, t � -I δ � p, t � = � R 2 A 0 2πσ 2 (t ) exp � � � - � p -q � 2 2 2σ 2 (t ) � � � s � q � dq - A 0 2πσ 2 (t ) exp � � � - � p � 2 2 2σ 2 (t ) � � � , with σ 2 (t ) = 2D � t -t 0 + σ 2 PSF � .
Since � R 2 s = 1, we get:

I r � p, t � -I δ � p, t � = A 0 2πσ 2 (t ) � R 2         exp � � � - � p -q � 2 2 2σ 2 (t ) � � � -exp � � � - � p � 2 2 2σ 2 (t ) � � �         s � q � dq .
And since s � q � = 0 when q � B r , the integration is made over B r :

I r � p, t � -I δ � p, t � = A 0 2πσ 2 (t ) � B r         exp � � � - � p -q � 2 2 2σ 2 (t ) � � � -exp � � � - � p � 2 2 2σ 2 (t ) � � �         s � q � dq .
We will now exhibit an upper bound of the bracket. To this end, let us introduce the following auxiliary function:

h : R + → R + x � → exp � � - x 2 2σ 2 (t ) � � .
From the mean value theorem, we get:

� � � � h � � p -q � 2 � -h � � p � 2 � � � � � � � � � � p -q � 2 - � p � 2 � � � max � h � � � � q � 2 max � h � � , with � � � � h � � � p � 2 � � � � � = � p � 2 σ 2 (t ) exp � � � - � p � 2 2 2σ 2 (t ) � � � < 1 σ(t ) .
Thus:

� � � I r � p, t � -I δ � p, t � � � � = A 0 2πσ 2 (t ) � B r � � � � h � � p -q � 2 � -h � � q � 2 � � � � � s � q � dq < A 0 2πσ 2 (t ) � B r � q � 2 σ(t ) s � q � dq < A 0 2πσ 2 (t ) r σ(t ) � B r s � q � dq < A 0 2πσ 2 (t ) r σ(t ) < A 0 2πσ 2 (t ) r σ PSF .

�

This shows that, whatever the vesicle profile, if its diameter is small enough with respect to the PSF width, the point source model remains accurate. Otherwise, the model becomes accurate only after a given time interval, when σ(t ) is high enough with respect to the vesicle diameter.

Three dimensional diffusion in the cytosol for Rab11

Proteins diffusing in the cytosol also follow the Fick's second law (6.12), but the Laplacian operator now acts in a 3D space, that is:

∂C d ∂t = D ∆C d = D � � � ∂ 2 C d ∂x 2 + ∂ 2 C d ∂y 2 + ∂ 2 C d ∂z 2 � � � . (6.29)
This yields a Green's function analogous to (6.14):

Φ � p, t � = 1 � 4πD � t -t 0 � � 3/2 exp � � � - � p -p 0 � 2 2 4D � t -t 0 � � � � , (6.30)
where p = �

x, y, z � T .

However, as fluorophores are present at z > 0, the PSF cannot be reduced to a simple two-dimensional Gaussian profile, as done for diffusion along the plasma membrane. Let us derive the intensity model for the point source in three dimensions observed with TIRFM.

The exponential decay of the excitation field of TIRFM must be integrated, which gives, with d the penetration depth:

C d (x, y, t ) = � ∞ 0 C d (x, y, z, t ) exp � - z d � dz = C 0 � 4πD � t -t 0 � � 3/2 exp � � � - � x -x 0 � 2 + � y -y 0 � 2 4D � t -t 0 � � � � � ∞ 0 exp � � - z 2 4D � t -t 0 � - z d � � dz = C 0 � 4πD � t -t 0 � � 3/2 exp � � � - � x -x 0 � 2 + � y -y 0 � 2 4D � t -t 0 � � � � × � πD � t -t 0 � exp � � D � t -t 0 � d 2 � �           1 -erf � � � � � � D � t -t 0 � d � � � � �           = C 0 8πD � t -t 0 � exp � � � - � x -x 0 � 2 + � y -y 0 � 2 4D � t -t 0 � � � � exp � � D � t -t 0 � d 2 � �           1 -erf � � � � � � D � t -t 0 � d � � � � �           . (6.31)
This equation could be integrated numerically, so that parameters could be estimated, but an approximation will exhibit the issue of studying the 3D+t model in 2D+t image sequences. Let us rely on the following exponential approximation of the error function [START_REF] Chiani | New exponential bounds and approximations for the computation of error probability in fading channels[END_REF]: As plotted in Figure 6.8, the approximation is precise when x � 1. For x = D � tt 0 � /d 2 , this condition corresponds to tt 0 � d 2 /D, which rapidly holds after t 0 ; orders of magnitude of d and D are 1 px and 1 px 2 /frame so that tt 0 � d 2 /D at t = t 0 + 1 already. In this case, we get:

erf � √ x � ≈ 1 - 1 2 exp(-x) . (6.32) 0 1 2 3 4 5 0 0.5 1 x erf � √ x � 1 -1 2 exp(x) 0
erf � � � � � � D � t -t 0 � d � � � � � ≈ 1 - 1 2 exp � � - D � t -t 0 � d 2 � � . (6.33)
This yields:

C d (p, t ) ≈ C 0 16πD � t -t 0 � exp � � � - � x -x 0 � 2 + � y -y 0 � 2 4D � t -t 0 � � � � , (6.34)
that is, we get the same solution as for the point source problem in two dimensions, up to a factor. While this means that we can apply the same estimation procedure to get the diffusion coefficient of Rab11 for this model, this raises an issue. Actually, we are not precisely interested in the diffusion coefficient estimation, but in the selection of a dynamical model. Therefore, as two-and three-dimension diffusion models have very close forms, 2D TIRFM it not sufficient to reliably discriminate the models, even for the simplest two-and three-dimensional point source models. Naturally, the situation would be even more complex and intractable if we would have also considered a continuous release and dissociation models.

For now, we have no 3D+time TIRFM images with a frame rate sufficient to study sub-second dynamics. However, we are confident that the technology will be ready in a few months [START_REF] Boulanger | Fast high-resolution 3D total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging[END_REF], thus allowing to conduct the study of Rab11 dynamical behavior in the near future.

Proposed SSED model for TfR and Langerin

In order to take into account the non-instantaneous release of the proteins at fusion, we propose a new fusion model for TfR and Langerin. The proposed "small-extent source with exponential decay release" (SSED) model is briefly compared to the aforementioned existing models in Table 6.1. The SSED model consists in relaxing both the pointwise source (H1) and instantaneous release (H2) hypotheses of the point source model.

Continuous release of the concentration

In Figure 6.9, we compare the evolution of concentration over time in true and simulated diffusion image sequences. Precisely, fusion events in sequences showing TfR and Langerin are represented by the means of kymographs. A kymograph gives the evolution in time of a given image line, by appending its successive profiles. Namely, in Figure 6.9b, the line x = 136 is used to display two TfR fusion events. The evolution of line x = 161 in the second sequence shows another fusion event, for Langerin this time. Then, simulated kymographs are presented, both for the point source and SSED models. Apart from quantitative aspects, we observe that the respective intensity profiles are very different, especially long after vesicle fusion. In particular, the central intensity peak observed after t 0 in real image sequences (mainly for Langerin) cannot be obtained with the point source model.

As we have shown that the pointwise source hypothesis is indeed valid, the other hypothesis, the instantaneous release, must be wrong. Therefore, instead of assuming that all the particles start diffusing at t 0 , we now consider a continuous release of the particles, where each particle is expected to stay at p 0 during a certain amount of time τ called residence time: (6.35) This yields an exponential decay of the source concentration: This model is illustrated in Figure 6.10. Different underlying dynamics can yield such an exponential decay release, e.g., a narrow escape [START_REF] Schuss | The narrow escape problem for diffusion in cellular microdomains[END_REF], Singer et al., 2008,Schuss, 2012], corresponding to the release of proteins through a small aperture, or dissociation-like process [START_REF] Carrero | Using FRAP and mathematical modeling to determine the in vivo kinetics of nuclear proteins[END_REF], Michelman-Ribeiro et al., 2009, Im et al., 2013]. Other hypotheses could yield other continuous release models, but the exponential decay has the advantage of being easy to handle from an estimation point of view. Actually, another (maybe even simpler) model is the constant flow model, where a constant amount of concentration is released at each time instant. However, this linear decrease has to stop when all the material is released, resulting in non differentiable functions. We are aware that the exponential decay is a simple model, but it is both a plausible and practical choice.

∂C s ∂t = - 1 τ C s .
∀ t � t 0 , ∀ p ∈ R 2 , C s (p, t ) = C 0 δ � p -p 0 � exp � - t -t 0 τ � . ( 6
t 0 t 0 + 1 t 0 + 2 -6 -4 -2 0 2 4 6 0 0.5 1 x I (x, t ) t 0 t 0 + 1 t 0 + 2
The proteins which leave the source compartment diffuse in the plasma membrane. Thus, the flow between both compartments acts as a source in the diffusion equation (6.13), which yields:

∀ t > t 0 , ∀ p ∈ R 2 , ∂C d ∂t � p, t � -D ∆C d � p, t � = - ∂C s ∂t � p, t � = C 0 τ δ � p -p 0 � exp � - t -t 0 τ � . (6.37)
Still owing to the superposition principle, the solution of this equation is obtained by convolving the Green's function (6.14) with the source function (6.36):

∀ t > t 0 , ∀ p ∈ R 2 , C d (p, t ) = � t t 0 C 0 /τ 4πD (t -u ) exp � � � - u -t 0 τ - � p -p 0 � 2 2 4D (t -u ) � � � du . (6.38)
Finally, the observed concentration is the sum of both components: (6.39) and the intensity model, illustrated in Figure 6.10, is obtained after convolution by the PSF:

C = C s + C d ,
∀ t � t 0 , ∀ p ∈ R 2 , I � p, t � = A 0 /τ σ 2 PSF exp � � � - t τ - � p -p 0 � 2 2 2σ 2 PSF � � � + � t t 0 A 0 2D (t -u ) + σ 2 PSF exp � � � - u -t 0 τ - � p -p 0 � 2 2 4D (t -u ) + 2σ 2 PSF � � � du .
(6.40)

Arbitrary spatial source distribution

Interestingly, the proof of Section 6.4.2 holds for the SSED model, so that a pointwise source is not required either for equation (6.40) to be valid. The vesicle only needs to be small enough with no assumption on its shape.

This way, we have relaxed both the pointwise and instantaneous release hypotheses of the point source model used in [START_REF] Mele | Automatic identification of fusion events in TIRF microscopy image sequences[END_REF].

Discussion

In this chapter, we have constructed new diffusion models for vesicle fusion to plasma membrane. Both for Rab11 and transmembrane proteins, we have shown that a small vesicle model is equivalent to a pointwise model, which relaxes this generally assumed hypothesis, and enables the use of the point source model with fewer constraints than in the literature.

Unfortunately, we have also shown that 2D and 3D diffusion models observed by the means of 2D TIRFM are barely distinguishable. Therefore, future experiments should rely on the nascent 3D TIRFM [START_REF] Boulanger | Fast high-resolution 3D total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging[END_REF], which should provide superior data quality to allow for better dynamics classification. When such image sequences will show up, estimation methods will be extended to the depth dimension.

In the meanwhile, the diffusion coefficient can still be estimated in 2D TIRFM image sequences (see Chapter 8).

As for the continuous release introduced for transmembrane proteins, we have proposed an exponential decay model. While more complex by one parameter than the point source model, the so-built SSED better corresponds to observed dynamics, and can still be estimated, as proposed in Chapter 8.

Further extensions to the SSED model could target the so-called kiss-and-run dynamics [Rizzoli and Jahn, 2007,Miklavc et al., ], assuming that the vesicle only partially empties before reforming and leaving the membrane.

Detection of diffusion dynamics

B ����� ���������� the parameters of the models representing the membrane fusion, the events of interest, which we name fusion events, must be detected.

Fusion events are defined by the space-time point at which a vesicle fusion starts, that is, the location and time at which the protein of interest is released to the plasma membrane or cytosol.

Depending on the protein observed in the image sequence, there are three cases:

• Rab11 is tagged with mCherry;

• TfR or Langerin is tagged with pHluorin;

• Dual View is used to observe both Rab11-mCherry and TfR-pHluorin or both Rab11-mCherry and Langerin-pHluorin.

In the first case, the vesicle is observed while moving before diffusing, while in the second case, owing to the pH-sensitivity of pHluorin, the fluorescence intensity is very low before fusion. This imposes different strategies to detect fusion events.

In the case of Rab11-mCherry, we propose a likelihood test to distinguish translation from diffusion dynamics. For transmembrane proteins tagged with pHlurin, fusion events start with a spot appearance, so a simple spot detection can be applied to the temporal frame difference sequence. This method can also be employed to study events in Rab11 images when Dual View with TfR or Langerin is available.

The chapter is composed of two sections. The detection of fusion events in image sequences depicting only Rab11 is addressed in Section 7.1. In Section 7.2, the properties of pHluorin are exploited to detect the events in image sequences depicting transmembrane proteins.

Detection based on Rab11 dynamics classification

In order to detect Rab11 fusion events, we propose a method based on a likelihood test to select the most appropriate dynamical model representing the vesicle dynamics in each frame. Specifically, we test the translation model, valid for the transport to docking steps, against the diffusion model, valid from fusion onward. For each detected vesicle, the transition frame, that is the first frame when the diffusion model is selected, merely corresponds to t 0 . The fusion location p 0 corresponds to the vesicle location at t 0 .

Let us assume for now that the parameters for the translation and diffusion models have been estimated -this will be presented in the next chapter. Then, we will merely select the most likely model. We have thus developed a three-step detection method. First, spots are detected relying on the spot detector proposed in Part I. Secondly, static spots are discarded owing to a statistical argument; only moving and diffusing spots then remain. Dynamics are finally classified as either translation or diffusion. The fusion event begins when the spot dynamics changes from the former to the latter.

The whole classification process only relies on three consecutive frames used to detect time-varying vesicles and evaluate temporal gradients.

Segmentation of time-varying vesicles

We first need to extract the vesicles that apparently change over time, or time-varying vesicles, in the TIRFM image sequence, that is, the vesicles which are moving or diffusing. The segmentation of the time-varying vesicles combines the detection of spatial spots and of significant temporal changes.

Spatial spot segmentation

We naturally use our adaptive thresholding of LoG images with auto-selected scale (ATLAS) detector (see Part I) which provides a binary vesicle presence map. In each frame I (t ) of the sequence, we end up with a set of connected components

V(t ) = � V i (t ) � .

Detection of temporal changes

In order to handle the temporal dimension of the segmentation of the time-varying vesicles, we examine the evolution of the segmentation map over time by using the change detection algorithm introduced in [Boulanger et al., 2010a � centered at that location and we compute the minimum of the intensity similarity dis- The minimum distance d min for all the pixels in the image sequence is assumed to follow a generalized extreme value (GEV) distribution [Boulanger et al., 2010a], defined by its CDF as:

η i � p, t -1 � η � p, t � η i � p, t + 1 � x y t
F � d min � p, t � , α, β, κ � = exp        -� � 1 -κ d min � p, t � -β α � � 1/κ       , (7.1)
where α, β and κ are respectively the width, location and shape parameters of the GEV distribution. The parameters are estimated by using a mixed L-moments/maximum likelihood method [Boulanger et al., 2010a].

We set a probability of false alarm (PFA) in order to derive a threshold ( quantile of the GEV distribution) able to detect the significant temporal changes in the image between two time points. The resulting significant change map is denoted E (t ) (see [Boulanger et al., 2010a] for more details).

Map of time-varying vesicles

Finally, to obtain the set

V � (t ) = � V � i (t )
� of time-varying vesicles, we simply discard every vesicle of V(t ) whose intersection with E (t ) is empty:

V � (t ) = � V i (t ) ∈ V(t ) � � � ∃ q ∈ V i (t ) , E � p, t � = 1 � . (7.2)
The spatiotemporal segmentation workflow is illustrated in Figure 7.2. To save computation time, we evaluate the complete similarity distance map, estimate the associated GEV parameters, and deduce the change detection threshold from the first three frames of the sequence only. Afterwards, we evaluate the similarity distance only for points belonging to the vesicles V i (t ) in the subsequent images. Moreover, since the threshold has already been evaluated, not all distances must be evaluated for each vesicle.

Finally, only translating and diffusing vesicles belong to V � (t ).

Classification of dynamical events

As explained above, we have to decide for a given vesicleV � i (t ) at time point t , whether it undergoes a translational motion or a diffusion. The two competing hypotheses are the following ones:

• H 0 : the vesicle V � i (t ) is undergoing a translation with displacement w i , • H 1 : the vesicle V � i (t )
is undergoing a diffusion with coefficient D i . We consider statistical image models to derive the corresponding likelihood functions.

To assess local motion, we use the displaced frame difference (DFD) given by:

ε 0 � p, w i , t � = I � p + w i � p, t � , t + 1 � -I � p, t � . (7.3)
For diffusion modeling, we can straightforwardly infer the random variable ε 1 from the Fick's second law, which yields:

ε 1 � p, t � = I t � p, t � -D i ∆I � p, t � , (7.4)
where I t denotes the temporal gradient of I .

We assume that both variables ε 0 and ε 1 are independent and identically distributed, and follow a zero-mean Gaussian distribution. As a consequence, the likelihood functions are defined as the product of Gaussian densities over the vesicle con-

nected component V � i (t ).
Then, the decision relies on the following log-likelihood ratio test:

log � � � � ψ 1 � V � i (t ), D i � ψ 0 � V � i (t ), w i � � � � � > τ. (7.5)
If the ratio is greater than τ, H 1 is selected and H 0 otherwise.

Experimental results

Synthetic image sequences 1

In order to validate our approach, simple but realistic image sequences were simulated. The microtubules and cortical Actin are modeled with a set of 3D random splines. While one extremity of each microtubule is connected to a centrosome, the other one is located near the plasma membrane. The Actin lies mostly in the first 100 nm slice of the cell. Finally, the vesicles are modeled by a set of points whose evolution is governed by interactions with the cytoskeleton. The action of molecular motors such as myosin and kinesin is represented by an elastic force between the particle and the second nearest point to the cytoskeleton. The module of the force depends on the type of interaction between the particle and the cytoskeleton to account for the diversity of the molecular motors involved during vesicle traffic. The spatial confinement is handled with a dedicated potential and we consider an additional viscosity term in the evolution model. Let us note that only the vesicles evolve while the cytoskeleton remains still, which is true for this time scale. Finally, when a vesicle reaches the plasma membrane, a finite difference scheme is used to model the diffusion process by the Fick's second law.

To get a realistic amount of noise, we have estimated the peak signal-to-noise ratio Globally, the results on these challenging sequences are promising, but we will see that a more efficient and much simpler approach allows to reliably detect fusion events relying on Dual View system and pH-sensitivity of pHluorin.

Event detection using fluorophore pH-sensitivity

When TfR and Langerin are tagged with pHluorin, fluorescence increases with pH.

This tag is actually employed because the pH of the vesicle is about 5.5, while that of extracellular medium is expected to be 7. Acidic pH inside the vesicle before t 0 leads to very low pHluorin photon emission. When the vesicle fuses to the plasma

Fusion event detection

Hence, we aim at detecting localized rapid fluorescence increases of intensity in I . To this end, we rely on the temporal backward difference I δt defined as:

∀ p ∈ Ω I , ∀ t ∈ N * , I δt � p, t � = I � p, t � -I � p, t -1 � . (7.6) A fusion event e i = � p 0i , t 0i
� corresponds to a bright spot centered at p 0i in the map

I δt � t 0i
� . To detect events, we rely on an adapted version of the spot detection method ATLAS (see Part I).

First, the scale of the vesicles s � is automatically selected in a multiscale representation of the images I (t ). We use the first ten frames of the input sequence I , as it contains more spots than ten frames of I δt , and ten frames is more than enough to observe hundreds of spots.

In fact, scale selection in I (instead of I δt ) is appropriate, because numerous static spots are present, aside the few fusing vesicles. Therefore, the scale of the static spots is selected, which also correspond to the scale of diffusing spots, at the time they appear and start diffusing. The found scale corresponds to the scale of the spots we aim at detecting in I δt . Gaussian window, whose size is not critical. Its radius is set to 60 px, which is a tradeoff on background structure sizes in the processed images. The detection threshold is inferred pointwise from a PFA fixed to 1 × 10 -6 . We end up with a set of N fusion events � �

e i � i =1.
.N detected over the image sequence. Let ν i be a spatiotemporal neighborhood of e i and let I i denote the restriction of I to ν i . From observation of real TIRFM images and quantitative evaluation on simulations, we have chosen a spatial neighborhood of 21 × 21 pixels centered at � p 0i , and a temporal range of 5 frames from � t 0i to � t 0i + 4 for the point source model, and 20 frames from � t 0i to � t 0i + 19 for the SSED model. Background structures and static vesicles may exist close to the fusing spot, and then constitute outliers of the Gaussian spot model.

To get rid of them, we need to subtract the image background as follows. The background intensity b i � q � at image point q ∈ ν i is estimated as the median of the intensity values at q over the twenty frames acquired before � t 0i :

� b i � q � = median t 0i -20�t <t 0i I i � q, t � . (7.7) Then, subtracting � b i � q � for q ∈ ν i merely gives us the estimated foreground subse- quence � z i = I i -� b i .
As shown in Figure 7.8, most background structures are removed.

Prior selection of diffusion processes under the point source model

Nevertheless, we observe that detected fusion events are not always followed by a point source diffusion process. Actually, some detections correspond to spots which apparently only fade out without spreading. As shown in the previous chapter (Figure 6.9f), this case possibly corresponds to events with long residence time under the SSED model. However, to estimate the point source model parameters, we need to retain only spots which really satisfy it.

While classical goodness-of-fit tests can be employed after the model parameters have been estimated, we propose a method to select diffusing events before estimating all the parameters, which reduces computation time, as the diffusion coefficient will be estimated only for actual diffusion processes.

Let us focus on a given event e i detected at � p 0i in frame � t 0i . If e i corresponds to a point source diffusion process, equation (6.23) holds in z i . Let us subscript with i the parameters associated to the diffusion process in ν i .

• t 0i is the fusion event beginning;

• p 0i is the fusion event location;

• A i is the initial spot amplitude;

• D i is the diffusion coefficient.
According to equation (6.14), the protein concentration is Gaussian-distributed in every frame after t 0 . Let A i (t ) and σ 2 i (t ) respectively denote the amplitude and variance of the observed Gaussian spot at time instant t . They are supposed to vary during the diffusion process. We will use � t 0i as an approximation of t 0i -as shown below, � t 0i could even be chosen arbitrarily higher than t 0i to estimate D i . Then,

A i � � t 0i � and σ 2 i � � t 0i
� are the "initial" spot amplitude and variance, and equation (6.23) yields:

z i � p, t � = A i (t ) exp � � � - � p -p 0i � 2 2 2σ 2 i (t ) � � � , (7.8)
where

A i (t ) = A 0i 2D i � t -t 0i � + σ 2 0i (7.9) and σ 2 i (t ) = 2D i � t -t 0i � + σ 2 0i .
(7.10) Thus, we get: .11) In this case, the series �

σ 2 i (t ) = A 0i A i (t ) . ( 7
σ 2 i (t ) � t �t 0i and � A -1 i (t ) � t �t 0i
are proportional and, a fortiori, correlated. Therefore, we discard every detected fusion event for which the (empirical) correlation is insufficient, by hypothesis testing. To estimate σ 2 i (t ) and A i (t ), we resort to a Gauss-Newton algorithm to fit a Gaussian spot in each frame of � z i . Then, to recognize diffusion processes, we rely on the Spearman's rank correlation coefficient ρ, which is more robust to outliers than Pearson's test [McDonald, 2014]

�

The test statistics for the correlation is t s = � 3ρ 2 / � 1ρ 2 , which is t-distributed with three degrees of freedom under H 0 [McDonald, 2014]. Hence, we can infer a threshold on ρ from a given rate of type I error α. As mentioned above, we use 5 consecutive frames to estimate the diffusion coefficient. This results in a threshold on ρ of -0.9

for a rate of type I error of α = 5% [McDonald, 2014], which we have chosen. We obtained better results with this correlation test than with the variance increase test proposed in [START_REF] Cortes | Efficient annotation of vesicle dynamics video microscopy[END_REF].

The whole proposed detection procedure is summarized in Algorithm 1.

Quantitative evaluation of the diffusive event detection

To evaluate the ability of the method to only retain diffusion events, experiments were conducted on the TIRFM sequence introduced in Figure 7.7. 90 fusion events were detected and then manually labeled as diffusion, non-diffusion, or debatable if the class was unclear. α: rate of type I error; N : detected fusion events; #deb: debatable events; #diff: diffusion events (found/ground-truth); #n-diff: non-diffusion events (found/ground-truth); #miss: missed diffusion events; #err: non-diffusion events erroneously labeled as diffusion.

As summarized in Table 7.1, very satisfactory precision results are obtained (precision ranges from 87% to 92% depending on α). To estimate the diffusion coefficient, it is important to obtain as few erroneous detections as possible, in order to limit the number of irrelevant measurements. Recall is not crucial here. The few classification errors are due to a second spot moving inside ν i , thus affecting the Gaussian fit. As a tradeoff between precision and recall, we will fix the error rate to α = 5%. However, as shown in the next chapter, the choice of α is not critical.

Detection from pH-sensitivity for Rab11 with Dual View

When both Rab11 and a pHluorin-tagged protein are imaged using Dual View, another detection method than the translation/diffusion classification can be employed. Indeed, when detecting an event e i in the pHluorin channel, we directly obtain the space-time coordinates e i = � p 0i , t 0i � , which are indeed valid in the Rab11-mCherry channel as well.

Also, this method does not rely on any particular dynamical model, so it is compatible with any diffusion model for both imaged proteins, while the classification proposed in Section 7.1 does not handle Rab11 dissociation, for example.

Discussion

We have proposed two main ways of detecting diffusion dynamics. The first approach consists in classifying vesicle dynamics into translation or diffusion, so the fusion event starts at the first frame where the vesicle dynamics are classified as diffusion.

While well suited to Rab11 under the point source model, this solution is not adequate for proteins observed with pH-sensitive tags, since the translation is not visible in that case. Also, as it is specific to the point source model, another approach must be used for the SSED model.

Fortunately, fusion event detection is much easier for TfR and Langerin thanks to the properties of pHluorin, since a fusion event corresponds to a spot appearance in the image sequence. Hence, our spot detector ATLAS is applied to the sequence of temporal differences to locate fusing vesicles.

When working with the point source model instead of the SSED model, a selection is performed among the detected spot appearance events, relying on a statistical test.

This allows to discard non-point source events to later obtain reliable estimates of the diffusion coefficient.

Finally, using Dual View, this simple method enables the study of Rab11 diffusion after detecting the fusion events in the TfR or Langerin channel.

Estimation of the biological parameters I � ���� �������, we propose methods to estimate the biological parameters of the previously introduced translation, point source and SSED models.

First, in order to classify Rab11 dynamics with the method presented in Section 7.1, we propose an estimation method for the translation and point source models which can be performed over a short temporal interval. As explained in the previous chapter, using only three consecutive frames, the method enables to decide the type of dynamics undergone by the vesicles on a frame-by-frame basis.

Then, methods requiring a small temporal integration (5 frames) are presented to estimate the parameters of the point source model. In addition to existing estimation methods, we propose new ones and compare them, in order to select the best method. Quantitative results for the point source model show the superiority of the intensity fitting method, which will then be extended to the SSED model. Since the model is more complex, several improvements are introduced to better fit the data.

We will employ this method to compare the dynamical behaviors of TfR and Langerin in real TIRFM image sequences. Specifically, distributions of the residence time τ and diffusion coefficient D will be compared for the two transmembrane proteins.

The chapter organization is the following. Existing methods for estimating the diffusion coefficient are presented in Section 8.1. We propose a method to estimate vesicle drift (translation) in Section 8.2, and to estimate the point source model in Section 8.3. Experimental results show that our method performs best. As for transmembrane dynamics tackled in Section 8.4, the best estimation method for the point source model is improved to cope with the increased complexity. Then, we carry out several experiments which allow us to exhibit differences in the dynamical behaviors of TfR and Langerin. Finally, estimation methods are discussed in Section 8.5, where we also propose improvements for future experiments.

Existing diffusion estimation methods

Numerous methods were proposed to estimate the diffusion coefficient for all the existing models presented in Chapter 6. Nevertheless, three classes of methods can be distinguished:

• Methods based on single particle tracking (SPT);

• Fluorescence correlation spectroscopy (FCS), which relies on the spatial and/or temporal intensity correlation between spatially and/or temporally neighboring pixels;

• Intensity fitting methods in which an intensity model is formulated and estimated in a space-time volume of the microscopy image sequence.

Hypotheses on the diffusion models required by those methods are summarized in Table 8.1. Let us briefly describe each type of method to enlighten the limitations when applied to the modeling of vesicle fusion.

Single particle tracking methods

Here, we present the basic mean square displacement (MSD) approach to show SPT methods are inappropriate to study vesicle fusion. Nevertheless, more complex SPT models were proposed in the literature and were successfully applied to real images under more realistic hypotheses, although not at vesicle fusion. In particular, authors of [START_REF] Schuss | The narrow escape problem for diffusion in cellular microdomains[END_REF], Singer et al., 2008] propose to relax the planar membrane assumption, and to address the narrow escape issue introduced in the previous chapter. To our knowledge, however, it was never applied to vesicle fusion, because the followed protein cannot be selectively located in a fusing vesicle. As opposed to SPT methods, we will see that our SSED model relies on higher-level hypotheses to cope with this issue in the context of vesicle fusion, by incorporating a continuous release rate.

SPT-based methods require images where particles can be individually detected.

By detecting and tracking moving particles, the diffusion coefficient is generally recovered from the so-called MSD [START_REF] Sako | Compartmentalized structure of the plasma membrane for receptor movements as revealed by a nanometerlevel motion analysis[END_REF],Saxton and Jacobson, 1997,Kim et al., 2011]. MSD is a function of a time step ∆t , defined as the expected displacement of a particle i in an interval step ∆t :

MSD � ∆t � = � � k i (t ) -k i � t + ∆t � � 2 � t , (8.1)
where �•� t denotes the temporal averaging.

From the Wiener characterization of Brownian motion (6.1), it is established that the MSD linearly increases with time step:

MSD � ∆t � = 4D∆t , (8.2)
where D denotes the diffusion coefficient. Hence, if individual proteins are tracked over time, the diffusion coefficient can be straightforwardly estimated from empirical MSD measurements.

More recently, Hozé et al. [START_REF] Hozé | Heterogeneity of AMPA receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF] proposed a method for estimating the diffusion coefficient with a Eulerian approach. After tracking particles, statistical arguments are proposed to characterize the dynamical properties owing to a probability of passing at each position of a grid. This way, local properties are extracted instead of particle motion parameters.

However, in the images we deal with, single proteins cannot be resolved since they are too close from each other compared to the microscope resolution. MSD cannot be computed as a consequence. Other SPT-based methods are ruled out for the same reason.

Correlation fitting methods 1

Temporal image correlation spectroscopy (TICS) can be employed to estimate the diffusion coefficient of the lateral diffusion model even when single particles are not separated. They were successfully applied to images depicting to the stationary model, where a set of Gaussian spots undergo Brownian motion and/or drift [START_REF] Ohsugi | Lateral mobility of membrane-binding proteins in living cells measured by total internal reflection fluorescence correlation spectroscopy[END_REF], García-Sáez and Schwille, 2008, Macháň and Hof, 2010, García-Sáez et al., 2010, Kolin and Wiseman, 2007, Di Rienzo et al., 2014].

In TICS, the following correlation function G is introduced [Sergeev, 2004]:

G : R 2 → R � t, ∆t � � → � I � p, t � I � p, t + ∆t � � p � I � p, t � � 2 p , (8.3)
where �•� p denotes the spatial averaging.

C������ 8. E��������� �� ��� ���������� ����������

In the stationary model, following the Wiener characterization of the Brownian motion (6.1), G does not depend on t , id est:

G : N → R ∆t � → � I � p, t � I � p, t + ∆t � � � I � p, t � � 2 , (8.4)
where �•� is now the spatiotemporal averaging. However, while valid in classical membrane diffusion experiments [Sergeev, 2004, Kolin and[START_REF] Kolin | [END_REF], we will now show that the time independence assumption does not hold for the point source model and, a fortiori, for our vesicle fusion SSED model.

To our knowledge, the correlation function (8.3) has never been derived for the point source model. To simplify equations, denote σ2 (t ) = 2Dt + σ 2 0 . The numerator of (8.3) writes:

� I � p, t � I � p, t + ∆t � � p = 1 |Ω I | � Ω I A 0 σ 2 (t ) exp � � � - � p -p 0 � 2 2 2σ 2 (t ) � � � × A 0 σ 2 � t + ∆t � exp � � � - � p -p 0 � 2 2 2σ 2 � t + ∆t � � � � dp = A 2 0 |Ω I | σ 2 (t ) σ 2 � t + ∆t � � Ω I exp � � � - � p -p 0 � 2 2 2σ 2 (t ) - � p -p 0 � 2 2 2σ 2 � t + ∆t � � � � dp = A 2 0 |Ω I | σ 2 (t ) σ 2 � t + ∆t � � Ω I exp � � - σ 2 (t ) + σ 2 � t + ∆t � 2σ 2 (t ) σ 2 � t + ∆t � � p -p 0 � 2 2 � � dp
The exponential is negligible when p is not close to p 0 . Thus, the integral over Ω I is approximately the integral over R 2 , which is a Gauss integral, and we get:

� I � p, t � I � p, ∆t � � p ≈ πA 2 0 |Ω I | � D ∆t + σ 2 (t ) � .
Analogously, the denominator of G � t, ∆t � is a squared Gauss integral when integrating over R 2 instead of Ω I . Then, the following result is straightforward:

� I � p, t � � p ≈ 2πA 0 |Ω I | .
We finally obtain:

G � t, ∆t � ≈ |Ω I | 4π � D ∆t + 2Dt + σ 2 0 � , (8.5)
which depends on t . In the point source model, non stationarity comes from the pointwise initial particle distribution, which implies that a particle is expected to reach a point p � p 0 only after a given time which increases with �

pp 0 � Therefore, classical FCS methods cannot be applied to estimate the diffusion coefficient for vesicle fusion. However, a new correlation-based estimation method will be proposed in Section 8.3.2 and evaluated in Section 8.3.4.

Intensity fitting method

Intensity fitting was investigated in a few papers to estimate the diffusion coefficient of the stationary [START_REF] Fortun | Correlation and variational approaches for motion and diffusion estimation in fluorescence imaging[END_REF], point source [START_REF] Mele | Automatic identification of fusion events in TIRF microscopy image sequences[END_REF] and FRAP models [START_REF] Seiffert | Systematic evaluation of FRAP experiments performed in a confocal laser scanning microscope[END_REF].

In [START_REF] Fortun | Correlation and variational approaches for motion and diffusion estimation in fluorescence imaging[END_REF], an optical flow estimation method was proposed to estimate the diffusion coefficient in the stationary model. A global variational approach is used, where a regularization term penalizes high gradients, and a data fidelity term is derived from the Fick's second law. As opposed to correlation fitting which assume a spatially and temporally constant diffusion coefficient, the method proposed in [START_REF] Fortun | Correlation and variational approaches for motion and diffusion estimation in fluorescence imaging[END_REF] is able to estimate pointwise a varying diffusion coefficient.

Seiffert and Oppermann [START_REF] Seiffert | Systematic evaluation of FRAP experiments performed in a confocal laser scanning microscope[END_REF] proposed a diffusion coefficient estimation method for the FRAP model, based on frame-by-frame Gaussian fitting. As presented in Section 6.2.5, in the FRAP model, the fluorescence hole has a Gaussian profile of linearly increasing variance σ 2 (t ) and decreasing amplitude A (t ):

I � p, t � = A (0) -A (t ) exp � � � - � p -p 0 � 2 2 2σ 2 (t ) � � � (8.6) with A (t ) = A 0 2D � t -t 0 � + σ 2 PSF + σ 2 beam (8.7) and σ 2 (t ) = 2D � t -t 0 � + σ 2 PSF + σ 2 beam . (8.8)
After fitting a Gaussian spot model in each frame I (t ), the amplitude model (8.7) and variance model (8.8) are fitted to recover the diffusion coefficient D , that is merely the slope of A -1 and σ 2 . We will see later that best estimation results are obtained with our estimation method.

Finally, Mele et al. [START_REF] Mele | Automatic identification of fusion events in TIRF microscopy image sequences[END_REF] proposed an approach for estimating the diffusion coefficient for vesicle fusion. Intensity is spatially averaged around fusion events to get a 1D+time signal whose decay can be fitted. Specifically, the averaging is performed in a region whose size depends on σ 2 PSF , which must be estimated beforehand. The final estimate of D hence depends on two complex fitting steps. It turns out that much information is lost and the two-step estimator is lacking accuracy.

In the remaining of this chapter, we propose and evaluate a method which can be used in any situations where an intensity model can be derived. In particular, the proposed SSED model will be estimated to better analyze membrane diffusion after vesicle fusion in real TIRFM sequences.

Proposed estimation method for the translation

Vesicle motion is quite slow -a few pixels per frame at most. Moreover, since we observe only the membrane region, the spatial density of vesicles is low enough to neglect about potential crossings of vesicles. Therefore, instead of using complex methods such as particle tracking (for a survey, see [START_REF] Chenouard | Objective comparison of particle tracking methods[END_REF]), it is sufficient to rely on local techniques of optical flow estimation to recover their frameby-frame displacement. Namely, after detecting the vesicles, we merely estimate the displacement with a least squares method, exploiting the brightness constancy constraint equation (BCCE), which is valid since photobleaching is negligible in the short time interval considered (the frame period). BCCE leads to the well-known motion equation [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF]:

I t � p, t � + I x � p, t � u i + I y � p, t � v i = 0 , (8.9)
where I t , I x and I y are the temporal and spatial derivatives, and u i and v i represent the

x and y components of the displacement w i of every point of the vesicle, respectively. The temporal and spatial derivatives of the image intensities must be estimated at each point of the vesicle support. While simple, the finite difference approach is very noise-sensitive. Therefore, we first apply a low-pass filtering with a Gaussian kernel, whose variance is the previously selected scale s � (see Part I).

The velocity vector w

i = � u i , v i � T of the vesicle of spatial support V � i (t )
is easily obtained [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF]:

w i = -� � � � q ∈V � i (t ) I xx � q, t � � q ∈V � i (t ) I x � q, t � I y � q, t � � q ∈V � i (t ) I x � q, t � I y � q, t � � q ∈V � i (t ) I y y � q, t � � � � -1 � � � � q ∈V � i (t ) I x � q, t � I t � q, t � � q ∈V � i (t ) I y � q, t � I t � q, t � � � � , (8.10)
where I xx and I y y are the second spatial derivatives of the image intensity.

8.3 Proposed estimation method for the point source model

Frame-by-frame estimation by exploiting the Fick's second law

The diffusion coefficient D i of the vesicle of support V � i (t ) can straightforwardly be estimated from the Fick's second law (6.12), with a linear least square fitting. In order to estimate D i more robustly, we do not only consider the points ofV � i (t ) at time t , but also their corresponding points in the previous and next frames, respectively at times t -1 and t + 1: (8.11) We also apply beforehand the aforementioned low-pass Gaussian filter of variance s � to smooth the Laplacian field, which is very sensitive to noise. This filtering does not impact the estimation, as the Fick's second law is invariant to Gaussian convolution.

D i = � q ∈V � i (t ) � t +1 τ=t -1 ∆I � q, τ � I t � q, τ � � q ∈V � i (t ) � t +1 τ=t -1 � ∆I � q, τ � � 2 .

Correlation fitting method

We have shown in Section 8.1.2 that classical STICS methods cannot be applied to estimate the diffusion coefficient of vesicle fusion models. However, we propose a simple correlation-based method to estimate the diffusion coefficient under the point source model. From (8.5), we can infer that G -1 is a linear function of D and σ 2 0 as:

G -1 � t, ∆t � ≈ 8π |Ω I | Dt + 4π
|Ω I | D ∆t + 4π |Ω I | σ 2 0 . (8.12)
Interestingly, the result is independent of A 0 and p 0 , so that D and σ 2 0 can be estimated with a standard two-dimensional linear regression. This method will be evaluated in Section 8.4.1.

Unfortunately, we were not able to extend this method to the SSED model, whose correlation function is much more complex. Indeed, non-stationarity is even reinforced in the SSED model by the time-varying source, which introduces another dependence on t . Moreover, those correlation-based methods make the assumption that every particle undergoes a Brownian motion. This hypothesis cannot be assumed in the SSED model because a stationary component is introduced, which prevents FCS from correctly estimating the related parameters. We must resort to intensity fitting methods.

Intensity fitting method

Since we have already estimated A i (t ) and σ 2 i (t ), as explained in Chapter 7, to detect fusion events, D i could be straightforwardly deduced from (7.9) or (7.10), as proposed in [START_REF] Seiffert | Systematic evaluation of FRAP experiments performed in a confocal laser scanning microscope[END_REF] for FRAP experiments. However, as demonstrated in the next section, better results are obtained by directly fitting the diffusion model equation (6.23) to the estimated foreground � z i in the space-time cuboid ν i (see Chapter 7). Specifically, the parameter vector

θ i = � p 0i , A 0i , σ 0i , D i � T is estimated by using a Gauss-Newton algorithm initialized with θ (init) i = � � p 0i , s � max ν i � z i , √ s � , 0 � T : � θ i = arg min θ i � (q,t )∈νi r 2 i � q, t � (8.13) with r i � q, t � = � z i � q, t � -A i (t ) exp � � � - � q -p 0i � 2 2 2σ 2 i (t ) � � � .

Quantitative comparison of estimation methods 2

To carry out an objective and comparative evaluation, we have simulated a collection of 300 image sequences containing a single diffusing spot. Each sequence is made of five images of size 21×21 pixels to mimic real spatiotemporal patches ν i . As explained in Chapter 2, the noise in real TIRFM images is supposed to be Poisson-Gaussian, so the noise variance can be stabilized to end up with a Gaussian noise. We can work with simulated sequences corrupted by Gaussian noise. We create sequences with a random signal-to-noise ratio (SNR) between 1 and 10. The diffusion coefficient and the initial spot variance are randomly varied from 0.1 to 10 px2 /frame, and from 0.5 to 1.5 px 2 , respectively. As shown in Section 8.4.2, diffusion coefficients are of the order of magnitude of 1 px 2 /frame in TIRFM image sequences we are dealing with.

We have compared our results with the aforementioned diffusion coefficient estimation methods, namely:

• Direct use of Fick's second law presented in Section 8.3;

• Amplitude decay fitting based on (7.9);

• Variance increase fitting based on (7.10);

• Correlation fitting presented in Section 8.3.2.

Comparative results are reported in Figure 8.1. This clearly demonstrate that our new method outperforms the others. Overall, the estimation of the diffusion coefficient was improved, especially for fast diffusion processes (D � 1 px 2 /frame). This order of magnitude corresponds to real situations, as shown in Section 8.3.4. The mean absolute logarithmic errors (MALEs) of the five methods are respectively:

• 0.24 for the direct least-square estimation,

• 0.22 for the amplitude fitting,

• 0.06 for the variance fitting,

• 0.07 for the correlation fitting,

• 0.03 for our intensity fitting method.

Quantitative evaluation on real TIRFM image sequences

We report now results obtained in real sequences. Six TIRFM image sequences were acquired, depicting micropatterned M10 cells, as in Figure 7.6. TfR was tagged with pHluorin. For each sequence, 600 images of size 256×256 were acquired at 10 f/s. The overall processing time is less than 0.1 s/f on a 2.3 GHz 4-core laptop.

As reported in Figure 8.2a, varying the rate of type I error does not impact significantly the estimated diffusion coefficient statistics. For α = 5%, 124 diffusion events where detected through the six real sequences. The average diffusion coeffi- cient equals 1.1 px 2 /f. It converts to 0.28 µm 2 s -1 , which is coherent with related studies [Ohsugi et al., 2006, Sako andKusumi, 1994].

Finally, in order to estimate the robustness of the method, we have artificially reduced the frame rate by using only one frame over two. We use the same parameters (detection parameters, size of ν i , and α) as before. Estimation results are presented in Figure 8.2b. Again, the average diffusion coefficient is hardly modified, demonstrating the robustness of the proposed method.

However, the method is only valid for the point source model. Indeed, events with long residence time τ are discarded from the estimation, as explained in the previous chapter. In order to obtain statistics of D including fusion events with large τ, we will propose an extension of this estimation method to the SSED model.

Proposed estimation method for the SSED model

Let us now focus on the estimation of the SSED model parameters. The model is defined by (6.40), which has no closed form, but can be numerically solved, so that standard optimization algorithms could be used. However, the SSED model has one more parameter than the point source model, and we were not able to satisfyingly estimate the SSED model parameters in simulated sequences using the estimation procedure presented in Section 8.3.3. We need to design a more elaborate algorithm, described below.. First, a good approach is to estimate κ = τ -1 instead of τ, as it "reduces nonlinearity" in (6.40). Since the Gauss-Newton algorithm does not always converge to the global minimum, we adopt the Levenberg-Marquardt algorithm and the update scheme of [Nielsen, 1999]. � . This way, κ and D can be estimated with a regression operating in two dimensions only.

Finally, the initialization of the parameters, and in particular the initialization of κ, is crucial. Quantitative results show that the estimation of all parameters is accurate when κ is correctly initialized. Therefore, we propose another extension to the estimation procedure. Instead of estimating the parameters only once for each detected event, we start with different initialization vectors, and after running the optimization algorithm, we compare the associated residuals. To obtain the best estimation of κ and D , we select the run which minimizes the sum of squared residuals. In practice, as a tradeoff between accuracy and computation time, we have chosen the set � 0.1, 0.31, 1, 3.1, 10 � of initial values for κ (init) and � 0.1, 10 � for D (init) .

The overall estimation algorithm is given in Algorithm 2.

Quantitative evaluation on simulated sequences

The same procedure as for the point source model evaluation is used to evaluate the proposed estimation method for the SSED model. 300 synthetic image sequences of size 21 × 21 pixels and length 20 frames were generated with different parameters to mimic real � z i 's. Like for the point source model, we have randomly set the diffusion coefficient between 0.1 and 10 px 2 /f, and the PSF variance from 0.5 to 1.5 px 2 . As for the residence time τ, it varies between 0.1 and 10 frames. The SNR still ranges from 1 to 10. Logarithmic errors on the estimation of both κ and D are reported in Figure 8.3 for each sequence. As plotted in Figure 8.3a, the estimation of κ is less accurate than that of D , but we will see in the next subsection that the accuracy is largely sufficient to extract relevant information from real TIRFM images. Moreover, large errors are very rare. Over the 300 generated sequences, only 5 have an absolute logarithmic error higher than 0.5, and the MALE is quite low at 0.12. As for the estimation of D , results reported in Figure 8.3b are very good when κ is high enough. Indeed, this behavior is not a surprise, since, for low κ, the flow between C s and C d is very small. Consequently, few particles are available to estimate D (precisely the one undergoing a Brownian motion). On the contrary, when increasing κ, the SSED model tends to the point source model, and the estimation becomes more and more accurate as the amount of signal available to estimate D increases. When κ > 0.25, estimation of D is as precise as for the point source model with a MALE of 0.03. Including the worst estimates, the overall MALE for the diffusion coefficient is still very low at 0.06.

Comparison of TfR and Langerin behaviors

The behaviors of TfR and Langerin after fusion have been compared in [Cinquin, 2011], but no biophysical model was proposed. Instead, a simple 1D+time intensity signal was used to classify fusion events as slow or fast in collections of sequences depicting either TfR or Langerin. It was reported that the number of slow events in Langerin image sequences is much higher than in TfR images [Cinquin, 2011].

Similar results can be expected when estimating the SSED model with our method in similar image sequences. In the SSED model, slow events correspond to long residence times (i.e., small κ), while fast events correspond to short residence times (high κ).

We have applied the proposed detection and estimation procedure to sixteen real TIRFM image sequences, half of which depict TfR, and half of which depict Langerin.

As explained in Chapter 7, in order to only take into account reliable estimates, a chisquare goodness-of-fit test is performed [McDonald, 2014]. Like for the point source model in previous chapter, we fix the rate of type I error α to 5%. In the set of sequences depicting TfR, 3,147 diffusive events are detected, versus 4,223 for Langerin.

The results are gathered in Figure 8.4 in the form of four histograms of � κ and � D estimated in the sequences depicting TfR or Langerin. There are two main conclusions to be drawn. First, the histogram of � κ exhibit a strong peak around 10 2 for TfR which does not exist for Langerin. In contrast, much more slow events are found in Langerin sequences, around � κ ≈ 10 -1 . Regarding the release rate, our results are therefore consistent with those reported in [Cinquin, 2011], but supply a more precise description of the biophysical model with a quantified biophysical parameter.

The second conclusion concerns the diffusion coefficient statistics. Indeed, Langerin shows a much higher dispersion of the estimates than TfR. Various hypotheses could explain such different behaviors, related, for example, to structures present around the vesicle fusion location, interactions between proteins... These preliminary results demonstrate the relevance of our approach. Statements on κ and D raise biological questions, answers of which could lead to a better understanding of the structures and interactions involved in the recycling process.

Discussion and perspectives

We have designed methods to estimate three types of membrane dynamics using TIRFM, namely translation, 2D and 3D diffusions. First, we developed methods using only three frames to estimate the translation of the vesicle and its diffusion dynamics.

While simple and fast, the diffusion coefficient estimation is not sufficiently accurate to draw reliable conclusions on real TIRFM image sequences. As a consequence, we investigated a more efficient approach for the point source and SSED models, relying on a longer temporal integration (5 and 20 frames, respectively). After demonstrating the efficiency of the method on simulated sequences, we successfully applied it to real TIRFM images depicting TfR and Langerin.

The experiments demonstrated that the residence time and diffusion coefficient distributions of the two transmembrane proteins exhibit clear differences. A biological model to explain it remains to be elaborated. In collaboration with UMR 144, further work could focus on the design and validation of biological models which could explain the behavior differences. The estimation of the SSED model should also be carried out in sequences showing mutant and wild-type cells.

III Group dynamics

Labeled affine flow for group motion characterization G ���� �������� can be encountered in biological image sequences, for exam- ple images of cell tissues as depicted in Figure 9.1, embryogenesis, cell migration, and in many other types of image sequences, ranging from urban traffic surveillance to crowd motion monitoring. The latter is particularly interesting in terms of evaluation because many crowd motion analysis methods have been proposed in the literature. Therefore, this scope represents an excellent experimental support.

Analyzing videos of crowded scenes is of interest in many applications [START_REF] Zhan | Crowd analysis: a survey[END_REF], Chandola et al., 2009, Thida et al., 2013]. Needs can encompass crowd safety in big social events, exhibitions, sports events or musical shows, surveillance and monitoring in public transportation areas like subways, airports or railway stations, people behavior understanding in commercial venues. The huge amount of daily acquired videos urges to define automatic tools for processing the available data and/or assisting human operators depending upon the targeted applications. Goals can comprise classifying dynamic behaviors [START_REF] Hu | Learning motion patterns in crowded scenes using motion flow field[END_REF], Rodriguez et al., 2011, Solmaz et al., 2012, Zhou et al., 2013], identifying main followed paths [START_REF] Jodoin | Meta-tracking for video scene understanding[END_REF], Wang et al., 2011, Zhou et al., 2011], preventing critical situations [START_REF] Solmaz | Identifying behaviors in crowded scenes using stability analysis for dynamical systems[END_REF], and detecting abnormal behaviors or events [START_REF] Cong | Abnormal event detection in crowded scenes using sparse representation[END_REF], Kim and Grauman, 2009, Kratz and Nishino, 2009, Wu et al., 2014].

A similar image analysis process can be envisaged as well for other sets of moving elements such as urban traffic, animal flocks or groups of cells like those depicted in Figure 9.1. However, in contrast to most existing methods, our methods are said instantaneous, that is, they rely on short-term analysis and do not require any prior learning stage. More precisely, the proposed methods are based on simple affine motion models estimated from only two consecutive images. They do not involve long temporal integration or trajectory computation. Furthermore, in order to create methods as generic and reusable as possible, we do not introduce any individual ap- North, West, South, East. A finer subdivision could be handled as well if required.

These eight group motion classes can be related to the behaviors introduced in [START_REF] Solmaz | Identifying behaviors in crowded scenes using stability analysis for dynamical systems[END_REF], as summarized in Table 9.1. However, let us mention that our scheme is applicable to any point in the image, not only around few critical points as in [START_REF] Solmaz | Identifying behaviors in crowded scenes using stability analysis for dynamical systems[END_REF].

Our motion estimation method is divided into three main steps:

1. Detection of moving areas in the image;

2. Estimation of three affine motion models (translation, rotation, scaling) in each window of a collection of windows of different sizes; The first step is achieved with the motion detection algorithm [START_REF] Crivelli | Simultaneous motion detection and background reconstruction with a mixed-state conditional Markov random field[END_REF], which follows a background subtraction approach and involves a mixed-state conditional random field (see [START_REF] Crivelli | Simultaneous motion detection and background reconstruction with a mixed-state conditional Markov random field[END_REF], for details). The set of detected moving areas is denoted Ω � I with Ω � I ⊂ Ω I , where Ω I is the image domain. The motion detection algorithm is tuned by two parameters that we have kept fixed for all processed sequences, and their setting was not critical for our classification task. With step 2, we end up with a set of motion model candidates at every point p ∈ Ω � I .

Step 3 allows us to select at every point p ∈ Ω � I the most relevant motion model among these candidates with an information criterion.

The outline of the chapter is the following. After giving a brief overview of crowd motion analysis methods in Section 9.1, the proposed LAF map is introduced in Section 9.2. Two applications are then proposed to demonstrate the relevance of the LAF motion map: motion patterns recognition in Section 9.3 and recovery of dominant paths in Section 9.4. We finally discuss the methods in Section 9.5.

Related work

Important research efforts have been devoted to crowd analysis for several years [START_REF] Thida | A literature review on video analytics of crowded scenes[END_REF], Li et al., 2015]. Specialized descriptors have been designed to capture the dynamics of crowds motion from videos and have been used for a number of inference tasks in crowd analysis, such as pedestrian tracking [START_REF] Hu | Learning motion patterns in crowded scenes using motion flow field[END_REF], Rodriguez et al., 2009], group motion pattern or path classification [START_REF] Wang | Trajectory analysis and semantic region modeling using nonparametric hierarchical Bayesian models[END_REF], Zhou et al., 2011, Solmaz et al., 2012], or anomaly detection [START_REF] Basharat | Learning object motion patterns for anomaly detection and improved object detection[END_REF], Kratz and Nishino, 2009, Feng et al., 2010, Ryan et al., 2011, Chockalingam et al., 2013].

Existing crowd analysis methods usually exploit motion-based features computed on extended time intervals: spatio-temporal cuboids [START_REF] Feng | Online learning with selforganizing maps for anomaly detection in crowd scenes[END_REF], Kratz and Nishino, 2009,Rodriguez et al., 2011], tracklets [START_REF] Zhou | Random field topic model for semantic region analysis in crowded scenes from tracklets[END_REF] and mostly trajectories [START_REF] Cheriyadat | Detecting dominant motions in dense crowds[END_REF],Rodriguez et al., 2009,Solmaz et al., 2012,Wang et al., 2011,Zhou et al., 2012]. Moreover, most crowd analysis methods rely on dense optical flows [START_REF] Chockalingam | Localized anomaly detection via hierarchical integrated activity discovery[END_REF], Cong et al., 2013, Kim and Grauman, 2009, Mehran et al., 2009, Ryan et al., 2011, Solmaz et al., 2012, Zhou et al., 2013], temporal or spatiotemporal gradients [START_REF] Benezeth | Abnormality detection using low-level co-occuring events[END_REF], Boiman and Irani, 2007, Kratz and Nishino, 2009, Roshtkhari and Levine, 2013].

The recognition of prominent paths in the scene is investigated in [START_REF] Wang | Trajectory analysis and semantic region modeling using nonparametric hierarchical Bayesian models[END_REF] within the framework of hierarchical Dirichlet processes and latent topics, while in [START_REF] Zhou | Random field topic model for semantic region analysis in crowded scenes from tracklets[END_REF] random field topic models are introduced which account for spatial and temporal coherence between tracklets. The authors of [START_REF] Zhou | Measuring crowd collectiveness[END_REF] have proposed a collectiveness measure based on trajectories and local velocities.

Tracking pedestrians in dense crowds may require specific approaches. [START_REF] Idrees | Tracking in dense crowds using prominence and neighborhood motion concurrence[END_REF] made the assumption that neighboring pedestrians have comparable trajectories to track individuals in dense crowds. In [START_REF] Rodriguez | Datadriven crowd analysis in videos[END_REF], a dense optical flow is first computed to deduce words, composed of quantized positions and directions. Correlated topic models (CTM) are exploited to learn word priors on a large video database. These models are then exploited for tracking individuals in a dense crowd.

As for crowd behavior classification, [START_REF] Zhou | Coherent filtering: Detecting coherent motions from crowd clutters[END_REF] and [START_REF] Cheriyadat | Detecting dominant motions in dense crowds[END_REF] studied coherent and dominant crowd motions. [START_REF] Zhou | Coherent filtering: Detecting coherent motions from crowd clutters[END_REF] proposed to group moving points according to the so-called coherent neighbor invariance. The latter provides information both on the spatial proximity of data points and on the correlation over time of their velocity vectors. In [START_REF] Cheriyadat | Detecting dominant motions in dense crowds[END_REF], the trajectories are organized into clusters according to a longest common subsequence (LCSS) criterion. In [START_REF] Wang | Trajectory analysis and semantic region modeling using nonparametric hierarchical Bayesian models[END_REF], a method is developed for the recognition of semantic regions (prominent paths in the scene) within the framework of hierarchical Dirichlet processes and latent topics, while [START_REF] Zhou | Random field topic model for semantic region analysis in crowded scenes from tracklets[END_REF] introduced the so-called random field topic model for semantic region analysis to account for spatial and temporal coherence between tracklets.

To our knowledge, only [START_REF] Hu | Learning motion patterns in crowded scenes using motion flow field[END_REF] and [START_REF] Solmaz | Identifying behaviors in crowded scenes using stability analysis for dynamical systems[END_REF]

w θ � p � = � � � a 1 a 2 a 3 a 4 � � � �������������� A � � � x y � � � + � � � b 1 b 2 � � � ���� B , (9.1) with θ = � a 1 , a 2 , a 3 , a 4 , b 1 , b 2 � the model parameter vector.
In order to characterize the eight previously introduced group motion classes, only three specific affine motion models are necessary: translation (denoted T), scaling (S), and rotation (R) motions. They respectively correspond to the following 2 × 2 matrices A , as explained in [START_REF] François | Derivation of qualitative information in motion analysis[END_REF]]:

A T = � � � 0 0 0 0 � � � A S = � � � a 1 0 0 a 1 � � � A R = � � � 0 a 2 -a 2 0 � � � . (9.2)
The vector B is considered in any case, since it corresponds to the displacement of the origin of the coordinate system. Hence, for each motion model, only two (translation case) or three coefficients (scaling and rotation cases) have to be estimated, respectively:

θ T = (b 1 , b 2 ) θ S = (b 1 , b 2 , a 1 ) θ R = (b 1 , b 2 , a 2 ) . (9.3)
Since we do not know in advance the appropriate spatial support to estimate the motion models, we consider a collection W = �

W i � � � 1 � i � � W � �
of overlapping windows of various sizes -typically, 25%, 50%, and 100% of the image dimensions. For a given size, the overlap rate is 50%, so that a given point p belongs to four windows of that size (apart from border effects). An example is given in Figure 9.2.

We estimate the three motion models defined in (9.2) in every window, using the robust method [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF] based on a multiresolution and incremental scheme, which is available as an open-source software1 . The robust estimation allows us to capture the dominant motion if several motions are present inside the window, and to tolerate errors of the motion detection stage. Since the minimization of the robust penalty function amounts to an iteratively reweighted least squares (IRLS) procedure [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF], each point p is assigned at the end of the estimation process a weight representing its influence in the robust estimation. A point whose weight is close to 1 (namely greater than 0.5) is called an inlier. Let θ i, k be the parameters of the motion model k ∈ � T, S, R � , estimated in the window W i ∈ W .

The set of inliers for the model of parameter vector θ i, k is denoted by X i, k .

The conformity evaluation of a point p to a given motion model of parameters θ i, k is based on the displaced frame difference (DFD) and is defined by:

ε i, k � p, t � = I � p + w i,k � p, t � , t + 1 � -I � p, t � , (9.4)
where

w i, k � p, t
� is the displacement of p between frames t and t + 1 deduced from θ i,k according to (9.1). Conformity corresponds to ε i,k close to 0.

For every motion model k , in every window W i , we estimate both the motion parameters θ i, k and the empirical variance σ 2 i, k computed over the inliers, given by:

σ 2 i, k = 1 � X i, k � � q ∈X i,k ε 2 i, k � q � , ( 9

Candidate selection

The optimal motion model at p should best fit the real (unknown) local motion at p while being of the lowest possible complexity. We consider a local patch η � p � centered in p and we exploit the fitting variable (9.4) which is likely to be close to 0 for the correct velocity vector. Let us assume that the ε i, k � q, t � 's are i.i.d. variables over points q ∈ η � p � and follow a zero-mean Gaussian law of variance σ 2 i,k . Then, we can write the joint likelihood in the patch η � p � :

ψ � p, θ i, k � = 1 � 2πσ 2 i,k � � � � η (p)∩Ω � I � � � � � q ∈η (p)∩Ω � I exp � � � - ε 2 i, k � q � 2σ 2 i,k � � � . (9.6)
To penalize the complexity of the motion model, we resort to the Akaike information criterion corrected for small sample size (AICc) [START_REF] Burnham | Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach[END_REF].

The correction is especially useful when the sample size is small, which is our case here. The criterion is given by:

AICc � p, θ i, k � = -2 ln � ψ � p, θ i, k � � + 2n k + 2n k (n k + 1) � � � � η � p � ∩ Ω � I � � � � -n k -1 , (9.7)
where n k is the dimension of the motion model k , that is n k = 2 for T-motion model and n k = 3 for S-and R-motion models. Finally, we select the optimal motion model � k at p by minimizing the criterion:

θ i, � k = arg min θ i,k ∈Θ(p) AICc � p, θ i, k � . (9.8)

From motion candidates to motion classes

At every point p ∈ Ω � I , we have selected the motion model candidate, which is represented by:

• a motion type � k ∈ � T, S, R � ,
• a parameter vector θ i, � k and, consequently, an affine flow vector w

TN North b 1 + b 2 > 0 & b 1 -b 2 < 0 TW West b 1 + b 2 < 0 & b 1 -b 2 < 0 TS South b 1 + b 2 < 0 & b 1 -b 2 > 0 TE East b 1 + b 2 > 0 & b 1 -b 2 > 0 S SN Convergence a 1 < 0 SP Divergence a 1 > 0 R RN Clockwise a 2 < 0 RP Counterclockwise a 2 > 0
Southward translation (TS);

Eastward translation (TE).

The scaling motion type corresponds to two classes:

Convergence, or negative scaling (SN); Divergence (SP).

And the rotation motion type naturally yields two classes:

Clockwise rotation, or negative rotation (RN);

Counterclockwise rotation (RP).

The different classes are characterized by criteria summarized in Table 9.

Motion patterns recognition

The initial classification map � c is noisy, due to the pixelwise selection process, as illustrated in Figure 9.2c. In order to characterize group motions, we propose to regularize � c . As we will see in the subsequent subsections, by regularizing the map with simple weighted votes, we also directly obtain the motion patterns at the desired scale. While the voting procedure is lighter than classical regularization approaches, based for example on Markov fields or morphological operations, very good results are reported below. Moreover, the regularization algorithm complexity is linear with the image size and does not depend on the regularization scale.

Local class histograms

For each point p ∈ Ω � I , we first compute a local class histogram H � p � = � H c � p � � c ∈Γ
counting the number of points of each class in a neighborhood of p. To avoid block artifacts induced by square neighborhoods, and to put more weight on the pixels q closer to p, we rely on a window function g r of radius r . Hence, the bin H c of class c ∈ Γ is evaluated as:

H c � p � = � q ∈Ω � I g r � � q -p � 2 � δ c � � c � q � � , (9.9)
where

δ c � � c � q � � =            1 if � c � q � = c 0 otherwise.
Typically, we use a Gaussian window of standard deviation r for g r , and the computation time can be reduced by relying on recursive approximations of the Gaussian filter, such as the Deriche separable filter [Deriche, 1990]. Using this approximation, the computation of the histograms is proportional to the image size, or more precisely the area of Ω � I . It is about three times longer than using an unweighted square window.

Class map regularization with inhibition

We propose a two-step method to regularize the local class map � c . We first build a regularized motion type map K from the local type map e st i mk , before inferring the regularized class map C. This is accomplished with a two-round voting procedure.

The notation of the types and classes, and of the associated maps, is summarized in To obtain K , accumulating votes of the same motion type regardless of the underlying motion class is unavailing and even counterproductive in some situations.

As a consequence, points in V(p) locally assigned with classes clockwise rotation and counterclockwise rotation, are not considered as bearing corroborating information on rotation but rather contradictory information. The same holds for class pairs divergence and convergence. Therefore, we introduce the notion of inhibition and reinforcement in the regularization procedure.

The algorithm proceeds as follows. In a given V(p) at time instant t , we compute the occurrences of the eight group motion classes of Γ from the values the inhibition factor to them would be possible but more complex than for scaling and rotation types. Moreover, it would not be very useful as this case is very rare in practice. Conversely, classes of the T-motion type benefit from the reinforcement factor.

� c � q, t � , q ∈ V(p).
Nevertheless, non-coherent translations will be tackled in the next chapter.

The first round of the voting procedure consists in selecting the motion type K � p � , as illustrated in Figure 9.3. Three motion type scores H T , H S and H R are computed as follows:

H T � p � = H TN � p � + H TW � p � + H TS � p � + H TE � p � , (9.10) H S � p � = � � � H SN � p � -H SP � p � � � � , (9.11) and H R � p � = � � � H RN � p � -H RP � p � � � � . (9.
12)

The regularized motion type K � In turn, we are able to recover the motion patterns of the group behavior on a frame-by-frame basis without any prior learning stage and without any critical parameter setting. Indeed, for an image sequence of τ images, we end up with τ -1 instantaneous maps C(t ) of group motion classes.

K � p � max H T � p � + H TN � p � H TW � p � H TS � p � H TE � p � H S � p � | -| H SN � p � H SP � p � H R � p � | -| H RN � p � H RP � p � | -|

Experimental results

We have carried out experiments on real image sequences to validate our new group motion pattern recognition method. Classification errors are avoided by using both the AICc (instead of the maximum likelihood (ML) criterion used in our preliminary method [START_REF] Basset | Frame-byframe crowd motion classification from affine motion models[END_REF]) and the inhibition procedure. Irrelevant local motions generally comprise opposite motions in that case. As an example, if a runner goes away from someone, he/she comes closer to another neighbor. The tendency would be to classify those neighboring points of � c as Convergence and Divergence. In such a case, the inhibition factor helps reducing the score of the scaling motion type, and recovering the correct translation type.

In Figure 9.5, we investigate the impact of the size of the regularization window. Its radius r is varied from 40 to 80 pixels for the Marathon bend sequence, where different types of motion can be distinguished. To sum up, the window size affects the classification in the sense that motions smaller than the window size are smoothed out, so that this parameter must be set according to the desired classification "granularity".

Finally, let us present the classification results of the proposed method for the Tissue sequence introduced in Figure 9.1. It shows a group of cells, whose membrane is fluorescently labeled, evolving as a group. The membranes appear in the image as moving edges with similar apparent widths. Since the Laplacian of Gaussian (LoG)

filter is also an edge enhancer, we have used our ATLAS method to detect the membranes.

The regularization radius for the classification was set to 50 pixels, which is the In addition, smaller events are found over the sequence, such as cell divisions, as shown in Figure 9.6g,h. While not directly detected, due to the very low intensity of the nascent membranes, cell divisions deform the neighboring cells. This results in localized detection of divergence ( ) and rotation classes ( ).

Recovery of principal paths

Local path map

In this section, we investigate the recovery of the principal paths followed by the group, with a Eulerian approach. It involves the introduction of local paths determined, in each cell of a grid superimposed on the image, from the space-time average (c,f,i) Class color code is given in Table 9.2.

of the affine flow over five frames. We make the assumption that the motion is stationary, meaning that the paths followed by the group do not vary over the sequence. This is a natural assumption in the frame of the recovery of the principal paths defined as the most followed paths.

Principal paths are simply recovered by accumulating the local paths, and the longest principal path is the dominant one. Thus, as opposed to classical methods, the whole path recovery procedure only relies on local estimates obtained without tracking or advection, and using only five consecutive frames.

Let us first introduce the local path map P . At each position β of P will be associated a local path π i . The set of local paths Π = � π i � i =1 .. 8 corresponds to eight oriented translations: North, Northeast, East, Southeast, South, Southwest, West, Northwest.

Indeed, the definition of the π i 's precisely allows for a 8-connectivity walk in the grid.

The resolution of P is coarser than the pixelwise resolution of � c . Since we seek the principal paths in the observed scene, such a spatial integration is required, and it also permits to get more reliable local paths. We generally take blocks β of size 31×31, but the sensitivity of this parameter is very low, as shown below.

We consider the motion subfields formed at each time instant t by the ve-

locity vectors w i, � k � p, t
� given by the motion models selected at time t at each point p ∈ β according to equation (9.8). Then, we compute in each block β the mean velocity vector w � β, t � in space ant time from the motion subfields

� w i, � k � p, u � � � � � � p ∈ β, t -2 � u � t + 2 � over the video sequence: w � β, t � = t +2 � u =t -2 � p ∈ β w i, � k � p, u � . (9.15)
The short temporal integration is necessary to ensure that the dominant motion is recovered, as some frames may locally exhibit an "outlier" motion which we do not want to take into account in this section.

Principal and dominant paths

The principal paths followed by the group in the observed scene are then retrieved from the local path map P . Starting from one given block in the image, we straightforwardly reconstruct a global path by concatenating the local paths from block to block following at each step the direction given by each local path among the eight possible ones (N, NW, W, SW, S, SE, E, NE). Then, we cluster all these trajectories and the main clusters supply the principal paths in the observed scene. In particular, the dominant path followed by the group is simply the longest global path.

Let us stress that our approach does not require any tracking of moving points which is a difficult issue in videos of densely crowded scenes. It does not rely on an advection stage from optical flow vectors estimated on a regular grid and on a subsequent trajectory clustering. We can recover the whole set of global paths by starting in turn from every block in the image. Then, the longest path is automatically selected as the dominant one.

The method is illustrated in Figure 9.7.

Experimental results

We have tested our dominant path recovery approach on a wide range of video sequences. Results on a synthetic sequence (Obstacles) and several real sequences are presented in Figures 9.7 and 9.8. The accuracy of the detected global path is constrained by the connectivity of P (8 neighbors) and its resolution, which can be increased by simply reducing the size of the blocks β at the cost of higher computation time.

Figure 9.8 highlights the performance of the method in very different situations.

The effect of the block resolution can be observed in the Marathon lane sequence, where the true motion direction (East-Northeast) lies between two elements of Π.

The processing of the sequence results in the piecewise linear curve for the recovered dominant path (Figure 9.8e), but the dominant path is very close to a straight line for 9 × 9 blocks (Figure 9.8f ). In the other sequences (Marathon bend, Obstacles, Tissue and Shoal), the paths are very well recovered, regardless of the block size. Computing � c (t ) and global paths takes about 10 seconds per frame (on the other hand, computing C(t ) from � c (t ) takes less than 1 second) on a 4-core 2.3 GHz laptop, and could be easily parallelized both in space and time.

Discussion and perspectives

We have proposed a new data-driven overall approach for group motion analysis, which in that context is original both in terms of motion measurements used (LAF, local paths) and in terms of motion classification and path recovery criteria. In contrast to most existing methods which are based on trajectories or tracklets, our approach provides a group motion classification on a frame-by-frame and pixelwise basis. Candidate motion models are estimated in a collection of overlapping windows for each pair of images. An information criterion then allows us to select the best motion candidate in each point of the image domain, allowing to jointly extract the velocity vector and classify the local motion type.

We have then proposed an original method to classify coherent group motions in videos on a frame basis. The group motion classification is achieved with a decision tree regularized with majority votes and involving inhibition between opposed classes. Moreover, since the whole classification process only requires two consecutive frames, even short events can be captured. The algorithm is fast and does not require any learning stage, no fine parameter tuning, and no trajectories computation. The experiments we have carried out demonstrate the accuracy and efficiency of our approach in various real situations.

Finally, we have proposed a method to extract principal paths in a Eulerian way.

Local paths are easily obtained by averaging the velocity information of the LAF in small blocks. The principal paths followed by the group can be straightforwardly inferred from the space-time stationary local paths, without any tracking, clustering or learning stage.

Only three parameters are involved in the whole method: the sizes of the windows Future work should include addressing the issue of perspective effects, which can be undesirable in some applications, and focusing on a higher level analysis of the dominant paths.

D�������� ��� ���������������� �� ��������� ������ �� ����� ��������� 147 10

Anomaly detection and localization I � ���� �������, we aim at demonstrating the potential of LAF by relying on it for a much more demanding application than the two group motion analyses presented in the previous chapter. Indeed, we target the detection and precise localization of so-called abnormal events. The method will be evaluated on videos of crowded scenes which represent very challenging application data, but for which ground-truth can be built by hand. Moreover, apart from the interest we have for videos of crowded scenes for the purpose of method validation, automatic anomaly detection is a much demanded tool in the world of videosurveillance.

Crowd analysis usually requires intense human supervision. When the objective is to analyze anomalies in the scene, sparse events must be identified. This demands an all-time extreme attention from the watcher. Actually, to decide the necessary actions to counteract those anomalies, one has to focus on special events for further analysis, ignoring a vast majority of normal occurrences. This task becomes even more difficult in crowded scenes, where the behavioral complexity in different parts of the video can cause confusion and distraction. Thus, the need for automatic systems that are able to assist the monitoring process of crowded scenes has been growing steadily.

There is no unique definition of an abnormal event. It may depend on the context and the application. As in [START_REF] Chandola | Anomaly detection: A survey[END_REF], we consider in this work that anomalies are events that cannot be fit in a proper model obtained from a set of surrounding data. Some samples are given in different situations in Figure 10.1. This formulation is general enough to be of large practical interest. Anomaly is taken here in a broad acceptance of a different behavior with respect to context. It does not mean that the so-called abnormal behavior is necessarily malicious, dangerous, or forbidden. Specifically, we deal with pixel-based action detection and localization, which is finer than frame-based anomaly detection. We want to determine where and when anomalies occur in the successive images of the video sequence, and to be able to do it at the pixel level.

The desired solution, however, has to comply with a number of requirements.

First, the devised model has to be simple and generic enough so that it can be used in a wide range of applications. Secondly, the algorithm has to be fast. Computational performance is an important criterion looking towards real-time implementation [START_REF] Lu | Abnormal event detection at 150 fps in matlab[END_REF]. Finally, an anomalous event detection at the frame level does not provide enough information to a human supervisor in order to quickly take actions on the possible anomaly, as pointed out in [START_REF] Li | Visual abnormal behavior detection based on trajectory sparse reconstruction analysis[END_REF]. This is even more glaring in densely crowded scenes where it may be utterly difficult to quickly find the localization of subtle anomalies even when it is known beforehand that a given frame contains one. To tackle this problem, the method has to be able to localize anomalies both temporally and spatially. Specifically, we target the pixelwise segmentation of the anomalies, which has never been done as far as we know.

The common infrastructure for videosurveillance is a network of fixed cameras.

Thus, we aim at designing a general efficient solution for anomaly detection and localization in crowded scenes applicable to static cameras. We want it to be instantaneous (or on-line) by delivering frame-by-frame output with a one-frame lag only. Indeed, we use three successive frames to come to a decision at every pixel of every image. We do not formulate any explicit modeling of what a normal situation is supposed to be nor an abnormal one. We are only searching for local discrepancy configuration. As a consequence, we do not need any prior supervised learning. Thus, our method is selfadaptive by being data-driven. It exploits a dense map of local crowd motion classes obtained by maximizing a penalized likelihood criterion and using affine velocity vectors computed in a collection of overlapping windows, as explained in Chapter 9. This chapter presents a novel method for detecting and localizing anomalies in videos of crowded scenes. It is based on the extraction of local histograms of crowd motion classes over a dense set of patches. These motion class histograms, extracted from LAF, can be seen as relevant motion descriptors. A dedicated distance histogram will be specified and modeled, so that a simple probability of false alarm (PFA) is set by the user to fix the detection sensitivity.

The rest of the chapter is organized as follows. In Section 10.1, we review the related literature and previous work on crowd anomaly detection. The ad-hoc histogram distance is proposed in Section 10.2. Then, in Section 10.3, we fully describe our motion-based anomaly detection method and give insights about its main properties. In Section 10.4, we report a comparative and objective evaluation on numerous generated and real images sequences. Finally, we offer concluding comments in Section 10.5.

Related work

Several approaches have been investigated for anomaly detection in crowd videos.

Some methods target specific scenarios, or are specialized for certain types of video data. For instance, escape behaviors can be considered as a specific case of anomaly in surveillance videos [START_REF] Wu | A Bayesian model for crowd escape behavior detection[END_REF]. However, this is a global type of anomaly since the escape behavior is shared by all the people in the scene from a given starting time.

Here, we are interested in local anomaly detection, which means that only a (small) part of the crowd behaves in a different way. Other works are able to detect anomalies locally in videos and without an explicit definition of what the abnormality is. Among these, two main classes are found: trajectory-based methods [START_REF] Li | Visual abnormal behavior detection based on trajectory sparse reconstruction analysis[END_REF], Piciarelli et al., 2008,Stauffer and Grimson, 2000,Wu et al., 2010] and feature-based ones [START_REF] Adam | Robust real-time unusual event detection using multiple fixed-location monitors[END_REF], Antić and Ommer, 2011, Cong et al., 2013, Kim and Grauman, 2009, Kratz and Nishino, 2009, Li et al., 2014].

Trajectory-based methods make use of the relevant information embedded in object tracks [Porikli andHaga, 2004, Stauffer andGrimson, 2000]. Nevertheless, these methods are usually constrained to scenes where it is possible to perform foreground tracking, otherwise they are subject to a large amount of false positives, as pointed out by [START_REF] Adam | Robust real-time unusual event detection using multiple fixed-location monitors[END_REF]. In [START_REF] Wu | Chaotic invariants of Lagrangian particle trajectories for anomaly detection in crowded scenes[END_REF], representative trajectories are first extracted after particle advection and chaotic features are exploited. The normality is modeled by a Gaussian mixture model. A ML estimation with comparison to a predefined threshold enables to determine normal and abnormal frames. Then, anomalies are located within frames identified as abnormal. A different approach was investigated in [START_REF] Mehran | Abnormal crowd behavior detection using social force model[END_REF], still based on particle trajectories. Interaction forces between particles are introduced, which yield a force flow in every frame. Recognizing normal frames and abnormal ones in the video sequence is achieved using a bag-of-C������ 10. A������ ��������� ��� ������������ words approach involving a latent Dirichlet allocation (LDA) model. Anomalies are delineated in abnormal frames as regions with high force flow.

The method described in [START_REF] Cui | Abnormal detection using interaction energy potentials[END_REF] relied on tracked keypoints to calculate interaction energy potentials, and to separate normal and abnormal crowd behaviors with a support vector machine (SVM) classifier. A non parametric Bayesian framework is designed in [START_REF] Wang | Trajectory analysis and semantic region modeling using nonparametric hierarchical Bayesian models[END_REF], which can be used to detect anomalous trajectories. Trajectories are described as bags of words, composed of quantized positions and directions. A dual hierarchical Dirichlet process (Dual-HDP) is defined to cluster both words and trajectories. Unlikely trajectories are considered as anomalous ones.

On the other hand, feature-based approaches are less prone to depend on specific scenarios. In [START_REF] Kratz | Anomaly detection in extremely crowded scenes using spatio-temporal motion pattern models[END_REF], spatiotemporal intensity gradients are used, whose distribution over patches in normal situations is supposed to be Gaussian. The Gaussian parameters are learned on the training set. In [START_REF] Kim | Observe locally, infer globally: A space-time MRF for detecting abnormal activities with incremental updates[END_REF], a mixture of probabilistic principal component analysis (MPPCA) aims at modeling normal flow patterns, estimated over patches of the training video set.

The method [START_REF] Chockalingam | Localized anomaly detection via hierarchical integrated activity discovery[END_REF] relies on probabilistic latent sequential models (PLSM) learned from 15 to 45 minutes of the sequence to analyze, where the user knows that few anomalous events happen. The spatiotemporal compositions (STC) method [START_REF] Roshtkhari | An online, real-time learning method for detecting anomalities in videos using spatiotemporal compositions[END_REF] requires about a hundred initialization frames to start learning weights of so-called codewords representing normal behaviors. Afterwards, weights are updated on-line so that no other training sequences are required.

In [START_REF] Benezeth | Abnormality detection using low-level co-occuring events[END_REF], co-occurrence matrices for key pixels are embedded in a Markov random field formulation to describe the probability of abnormalities.

Mixture of dynamic texture (MDT) are introduced in [START_REF] Li | Visual abnormal behavior detection based on trajectory sparse reconstruction analysis[END_REF] with conditional random fields (CRF) to represent crowd behaviors and reported successful results on several datasets, but at the cost of sophisticated models that require intensive learning and high computation time. Other authors focused on giving explicit inclusion of spatial awareness, by subdividing the image in local regions or blocks, in order to obtain a good detection performance with less learning requirements [Adam et al., 2008, Boiman and[START_REF] Boiman | [END_REF].

Another approach was explored in [START_REF] Antić | Video parsing for abnormality detection[END_REF]. Vectors of spatiotemporal derivatives were utilized as input of a SVM classifier with linear kernel to support the foreground separation process. The latter feeds a graphical probabilistic model.

It is interesting to note that such method depends heavily on how well the foreground elements of a video dataset are separated, undermining a possible application for very crowded scenes. Social force models based on optical flow of particles, as introduced in [START_REF] Mehran | Abnormal crowd behavior detection using social force model[END_REF] is another example of descriptor used to detect anomalies.

H�������� ��������

Sparse representations have been increasingly adopted for anomaly detection, as the problem can be elegantly modeled with sparse linear combinations of representations in a training dataset [START_REF] Cong | Abnormal event detection in crowded scenes using sparse representation[END_REF], Zhao et al., 2011, Zhu et al., 2014]. Explicit image space subdivision can also benefit anomaly localization performance in sparse representation-based methods [START_REF] Biswas | Sparse representation based anomaly detection with enhanced local dictionaries[END_REF]. It is shown in [START_REF] Mo | Adaptive sparse representations for video anomaly detection[END_REF] that, by introducing nonlinearity into the sparse model, better data separation can be achieved. Also, some modifications can be made to the usual construction of the sparsity models by introducing small-scale least-square optimization steps [START_REF] Lu | Abnormal event detection at 150 fps in matlab[END_REF], sacrificing accuracy for the benefit of a fast implementation.

However, although elegant and sound, sparse representation methods for anomaly detection have not shown high performance in demanding videos.

To summarize, most methods need a large number of frames to perform anomaly detection, which is often stated at the frame level first. They usually require manually labeled sequences or frames to learn models of normal situations, anomalies being defined as outliers of the normal behavior model. In contrast, we aim at detecting and at the same time localizing anomalies directly on a pixelwise basis in every incoming frame following a data-driven "agnostic" approach. We do not build any model of the normal situation and rely on local statistical tests inspired from the center-surround biological vision process [START_REF] Cavanaugh | Nature and intearction of signals from the receptive field center and surround in macaque V1 neurons[END_REF], Tadin et al., 2006].

Histogram distance

From the local crowd motion classification maps presented in the previous chapter, we will locate the abnormal behaviors at every pixel p and at every time instant t .

As we seek for local anomalies, the classification is not regularized, that is, we rely on � c instead of C. The detection is based on simple statistics on the

� c � p, t
� values, inspired from the center-surround interaction operating in biological visual motion mechanisms [START_REF] Cavanaugh | Nature and intearction of signals from the receptive field center and surround in macaque V1 neurons[END_REF], Gao et al., 2008, Tadin et al., 2006]. The centersurround interaction has been recently exploited for action recognition [START_REF] Escobar | Action recognition via bio-inspired features: The richness of center-surround interaction[END_REF]. It was also influential in the anomaly detection method defined in [START_REF] Li | Anomaly detection and localization in crowded scenes[END_REF], but in a different way than ours. We consider spatial surrounding only, while it may be spatiotemporal in [START_REF] Li | Anomaly detection and localization in crowded scenes[END_REF] with an extended time interaction. We only evaluate distances between histograms whereas elaborated models are designed (and previously learned) and associated with a Markovian framework in [START_REF] Li | Anomaly detection and localization in crowded scenes[END_REF].

For every point p, we compute the occurrences of the eight motion classes of Γ in a patch ν �

p � centered in p, that is, the histogram H � p � of the � c � q, t � values for q ∈ ν � p � .
To increase the temporal stability of the method, we compute a spatiotemporal histogram over three consecutive frames so that ν � p � is a m × m × 3 space-time patch. For the translation class sub-histograms, we resort to the modulo distance introduced in [START_REF] Cha | On measuring the distance between histograms[END_REF] for sets of modulo measurements, that is, measurement values forming a ring, which is typically the case for the translation classes and their associated compass directions. It is expressed by: [START_REF] Cha | On measuring the distance between histograms[END_REF]. The distance d mod between elements is given by:

D mod � H T i � p � , H T � p � � = min U i ,U � � � � n � l,l � =1 d mod � u i (l ), u � l � � � � � � � , ( 10
d mod � u i (l ), u � l � � � =            � � � u i (l ) -u � l � � � � � if � � � u i (l ) -u � l � � � � � � |ΓT| 2 , � Γ T � -� � � u i (l ) -u � l � � � � � otherwise.
(10.3) where Γ T is the subset of the four translation classes, and � Γ T � = 4. Since the distance between two histograms can be expressed in terms of the distances of element measurement values, criterion (10.2) allows us to find the minimum difference of pair assignments between sets U and U i .

The interest of the modulo distance is that the distance between opposite translation directions (e.g., North and South directions) will be higher than the one between adjacent translation directions (e.g., North and West directions). Such a behavior is recommendable in our case, all the more so as we deal with a coarse direction quantization. Indeed, points undergoing a given slanted translation (e.g., North-West) could be shared out among two adjacent translation classes (North and West translation classes to continue our example).

After testing several usual histogram distances, we have retained the L 1 distance for comparing histograms H S, R i and H S, R : with equally weighted distances, since the ranges of the modulo and L 1 distances are similar as explained in [START_REF] Cha | On measuring the distance between histograms[END_REF].

D L 1 � H S, R i � p � , H S, R � p � � = � c ∈Γ S, R � � � � H S, R i � p, c

Pixel-level anomaly detection and localization

Since we have defined an abnormal behavior as revealed by an irregular motion pattern with respect to the surrounding motion patterns, an anomaly will be detected at point p if and only if: Ω � I follows a generalized extreme value (GEV) distribution [START_REF] Embrechts | Modelling extremal events for insurance and finance[END_REF], which is defined by:

F � ζ � p � , α, β, κ � = exp        -� � 1 -κ ζ � p � -β α � � 1/κ       , (10.7)
where α, β and κ are respectively the width, location and shape parameters of the GEV distribution. First, we have to estimate its parameters. We have adopted the mixed L-moments/maximum likelihood method described in [Boulanger et al., 2010a]. An example is given in Figure 10.2f. Then, we set a PFA to automatically set the threshold value λ as a quantile of the GEV distribution, in order to detect the local abnormal behaviors in the observed crowd. In practice, PFA will be set between 0.5 and 2%. As explained later, this value can be seen as the expected false positive rate (FPR), so the detection sensitivity can easily be set according to the targeted application.

Experimental results

To demonstrate the performance of our method, we need to show that it detects anomalies in abnormal situations, the anomalies are well segmented, and it does not detect anomalies in normal situations.

Simulated cases

We first deal with simulated cases to make an easy objective evaluation of different aspects of our method. We have taken a sequence depicting a normal behavior of We have computed receiver-operator characteristic (ROC) curves, that is, the evolution of true positive rate (TPR) versus FPR, for different sizes of local neighborhood ν p used to compute the AICc, and different sizes of patches ν � p � used to compute local histograms, as reported in Figure 10.3d,e. Since our method supplies pixelwise anomaly localization, TPR is the proportion of anomaly pixels that are correctly detected, and FPR is the proportion of normal points that are detected as anomalies.

The optimal size for the affine model selection neighborhood ν p is 3 × 3. The same size will be adopted for all the processed sequences. Concerning square patches ν � p � , best results are obtained with patch width between 31 and 51 pixels, which is close to the cropped region size. For this patch width range, the area under ROC curve is equal to 0.99.

Real scenes

To the best of our knowledge, there is no benchmark available for pixelwise anomaly localization in dense crowds, and such an evaluation has not been attempted so far because most existing methods are not designed for both detecting and localizing anomalies on a pixel basis. The UMN dataset and the web dataset of [START_REF] Mehran | Abnormal crowd behavior detection using social force model[END_REF] are concerned with globally abnormal crowd behaviors (as escape panics). The UCSD dataset [START_REF] Li | Anomaly detection and localization in crowded scenes[END_REF] mainly involves very sparsely crowded scenes. Besides, the pixel-level evaluation used in [START_REF] Li | Anomaly detection and localization in crowded scenes[END_REF] remains a frame-based evaluation taking into account a sufficient recall (40%) of abnormal regions defined by bounding boxes in the ground truth. We propose a more demanding evaluation criterion: we want to carry out a pixel-based objective evaluation. To that end, we have manually determined the pixelwise anomaly ground truth for two real sequences: the Wrong way sequence and one sequence from the PETS dataset. Thus, we are able to supply a truly pixel-based objective evaluation for real anomaly cases.

In Figure 10.5, we report results of our anomaly localization method on the Wrong way sequence (Figure 10.1b). In this crowd scene, one pedestrian is walking against the crowd. As shown at four different time instants, respectively in The bounding box of the man walking downward has a width of about 60 pixels all over the sequence, while its height in the image varies between 50 and 100 pixels due to partial occlusions. The best detection rate is obtained with square patches of 61pixel width. It is not possible to compare with [START_REF] Mehran | Abnormal crowd behavior detection using social force model[END_REF] since the authors have only provided ROC curves for abnormal frame detection and color-coded force flow maps, they did not supply any binary detection maps.

In Figure 10.6, we demonstrate the robustness of our method on a sequence of the PETS'09 dataset1 . In this sequence, we aim at detecting the only pedestrian walking toward the left of the frame and passing behind a group of people going in the opposite direction. The task is particularly challenging here, because the group often occludes the lonely pedestrian and the pedestrian is very small in the image. Yet, the lonely pedestrian is accurately segmented and no false alarm occurs as displayed in As expected, no anomalies were detected as shown in Figure 10.1a (apart from a few spurious points isolated in space and time). Finally, the Wrong way sequence was presented in Figure 10.5. In every case, the obtained FPR is close to the PFA parameter value since each plot in Figure 10.8 is close to the straight line of slope 1. The same holds for all the other sequences we have analyzed. Thus, a clear advantage of the proposed method is that the FPR can be easily controlled by the user through the setting of PFA.

Parameter setting and computation time

Overall, the anomaly localization is accurate and stable over time. The method involves two parameters only, easy to set: PFA value which can be seen as an expected FPR, and dimensions of ν which must be close to those of the expected anomaly bounding box -typically the size of the image projection of a pedestrian. It can be deduced from the scene characteristics and the camera setting. Furthermore, as reported in Figure 10.2e, Figure 10.5c and Figure 10.6c, the sensitivity of the patch size is low, so that it does not need to be finely tuned.

C������ 10. A������ ��������� ��� ������������

The computation time for the largest sequences we have processed (of image size 720 × 576) is of 6 to 10 seconds per frame on a laptop with a 4-core 2.3 GHz processor. Let us notice that using running sums, the computation load of the criterion AICc remains low and does not depend on the patch size. Even if our method is not real-time so far, we have described it as on-line, since by construction it delivers the anomaly map at each time instant with a one-frame lag only. It can be highly parallelized to eventually attain real-time processing. On-line detection is unreachable with approaches somehow relying on temporal integration, regardless of the processing power.

Discussion

We have developed an original data-driven method for anomaly detection in videos of dense crowds. It yields a pixelwise anomaly localization map at every time instant using only three successive frames. The method relies on the computation of parametric motion models, on the local selection of the crowd motion classes with the AICc penalized likelihood criterion, and on the bio-inspired evaluation, at every pixel, of a distance between center-surround histograms of the local crowd motion classes.

The method is simple to perform, there is no normal or abnormal behavior models to learn and it is almost parameter-free. It can handle a large category of situations. We have reported experimental results on both simulated and real cases. The objective evaluation has demonstrated the accuracy and reliability of the method.

Further work should investigate the impact of different geometric centersurround configuration of patches ν � p � and the weighting of the class histograms by the velocity vector magnitude, to fully take advantage of the LAF data. The latter could enable to discriminate not only with respect to the motion class, but also between different speeds in a same class.

various applications, including the characterization of the motion of the group as a whole, the recovery of main paths in the scene, and the detection and localization of abnormal behaviors in the group.

11.1 Main contributions and results

Spot detection

First, we have investigated methods for the detection of spots in TIRFM images. In such images sequences, vesicles appear as small bright spots of similar size. Therefore, our spot detection method relies on the similar size hypothesis as a prior knowledge to improve the detection performance. Considering the Laplacian of Gaussian (LoG) to enhance the spots, the best detection scores are obtained if the LoG scale is accurately selected beforehand. The optimal filter scale is automatically selected in a scale-space representation of the image. The spot detection then amounts to a pointwise thresholding of the LoG image computed at this scale, and the threshold adapts to the local image statistics. While designed in view of the exocytosis detection in TIRFM, the proposed adaptive thresholding of LoG images with auto-selected scale (ATLAS) method can be applied to a wide variety of images and delivers a leading performance.

Individual dynamics analysis

In the frame of the exocytosis study, there are two cases depending on the tagged protein. As for Rab11, it undergoes a translational motion before the vesicle fusion to plasma membrane, and then diffuses. Therefore, we have proposed a classification method to distinguish between the two dynamic situations. On the other hand, when observing a transmembrane protein tagged with pHluorin, vesicle fusion to plasma membrane is visible in the image sequences as spot appearances. In this case, we simply applied ATLAS to the sequence of temporal difference images to detect these spots.

Given the detected space-time fusion events, the biological models for the proteins diffusion can be explored. In particular, we have proposed a new fusion model, termed "small-extent source with exponential decay release" (SSED). As opposed to the existing point source model for which all the contents of the vesicle start diffusing at a given time t 0 , the vesicle is supposed to slowly release the proteins after fusion in the SSED model. This motivates the introduction of another biophysical parameter, the residence time, in addition to the diffusion coefficient.

Furthermore, the new SSED vesicle fusion model is then estimated with a method which proved to perform better than existing estimation methods on the simpler point source model. Dynamical behaviors have been characterized in a set of image sequences, with the use of the SSED model, to exhibit the different behaviors of the two transmembrane proteins TfR and Langerin.

Group dynamics analysis

In the last part of the thesis, we have focused on the development of LAF, a new motion measurement, which is more appropriate to group analysis than classical detectand-track methods, in the sense that it captures both the type and the numerical value of the local velocities without requiring any segmentation. Thus, LAF combines the advantage of quantitative and qualitative motion descriptors.

As an application, we have shown the relevance of this motion descriptor to characterize group motion by classifying motion patterns and by recovering the principal paths followed by the group. We have finally proposed a new method to detect and segment local anomalous behaviors, still building on LAF. All these methods only rely on very short temporal integration, and are therefore suitable for on-line processing, as opposed to most existing group motion analysis methods.

Ongoing and future work

We are currently working on a multiscale extension of our spot detector, relying on the same scale-space representation, but with a more elaborated selection criterion. Preliminary results were presented in Part I, but a more advanced method is envisaged.

It relies on a probabilistic formulation of the selection criterion, which allows us to set a meaningfulness threshold. This way, the user does not need to specify the number of scales. Preliminary results on the scale set selection itself are already promising, but the adaptive thresholding must be extended to cope with multiple scales. Then, detection performance must be evaluated to verify that the automatic scale set selection behaves better than multiscale detectors like the multiscale variance stabilizing transform (MS-VST).

The proposed SSED model will be used in other experimental conditions to investigate possible differences between mutant and wild cells. In particular, the team of Jean Salamero at Institut Curie aims at investigating the role and interactions of Actin in the recycling process. Ultimately, the dynamics of Rab11 are targeted, since open questions remain for this protein.

This leads us to the future work to be accomplished regarding Rab11. The behavior of Rab11 should be modeled and estimated thoroughly in 3D+time image sequences acquired with multi-angle TIRFM [START_REF] Boulanger | Fast high-resolution 3D total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging[END_REF]. At this time, preliminary results and discussions exposed in this thesis will constitute some basis for the development of 3D+time estimation methods.

It is already known that Rab11 finally returns to the cytosol at a given, but unknown, time. Understanding the timeline of the membrane fusion process, and in particular the position in time of Rab11's dissociation, is of primary interest for cell biology. Dissociation should therefore be modeled and estimated in 3D+time image sequences, constituting the short-term target of the study. In particular, the different Rab11 dynamical models presented in Part II should be challenged. Then, efforts should be put in the identification of interactions between proteins, in particular

Rab11 with TfR and Rab11 with Langerin, in order to point out potential behavioral differences.

At a low level modeling, other assumptions could be relaxed, for example the membrane isotropy, homogeneity and planeity. These improvements would naturally impact both the models and estimation procedures, and open up a large area of possible studies to end up with a fine model, and reliable estimation and classification methods.

As for group motion analysis, we have presented interesting results for crowd motion analysis throughout Part III, and shown some qualitative examples of cell group dynamics.

In video analysis, working on anomaly detection should be further investigated, for example with the use of the magnitude of the velocity vectors in the histogram descriptors or with another decision paradigm. This is currently investigated with Juan-Manuel Pérez-Rúa, a PhD student in the Serpico team.

In addition, the group motion study is still a proof of concept which we think could find a vast set of biological application domains. The agnostic approach we have proposed is interesting from an image processing point of view, but better results would probably be obtained by adding specific priors to cope with a particular issue. Anomaly detection is a classical topic in biology, found for example in tissue genesis.

As pointed out by a reviewer of this thesis, an interesting application could be active polar gels, where biological knowledge is available to enrich the affine models and therefore improve the biological relevance of the method. Given the immensity of the biological research landscape, we are confident that other applications could find their way in PhD theses to come. 28,32,44,47,62,89,95,137,141,142 PSF point spread function. 19,20,71,74,75,77,80,81,83,85,54,58,91,92 RNA ribonucleic acid. 14,138,140,141 SEF spot enhancing filter. 29 SIM structured illumination microscopy. 13

MALE

SNR signal-to-noise ratio. 18,19,29,30,38,41,48,52,106,107 SPT single particle tracking. 19, 99, 100

SSED "small-extent source with exponential decay release". 73,82,83,86,97,101,102,104,107,108 STICS spatiotemporal image correlation spectroscopy. 103 Transferrin receptor. 4,6,8,11,17,28,58,[65][66][67]69,82,83,87,94,99,108,111,143 TH top hat filter. 29, 30

TfR

TICS temporal image correlation spectroscopy. 100 TIRFM total internal reflection fluorescence microscopy. 1, 5, 8, 11, 17-20, 28, 29, 33-35, 38, 46, 50, 52, 55, 59, 62, 65, 72, 77, 81, 82, 86, 88, 91, 92, 95, 102, 104, 106, 108, 111, 143 TMA tissue micro-array. 59

TPR true positive rate. 92, 137, 138

WMP wavelet multiscale product. 30

  La sélection d'échelle s'appuie sur une représentation multi-échelle de l'image I à segmenter. La représentation multi-échelle est une famille d'images � L s � s ∈S obtenues par convolution avec des noyaux G s d'échelles s croissantes prises dans un ensemble de valeurs prédéfinies S ⊂ R * + : ∀ s ∈ S , L s = G s * I , où * désigne l'opération de convolution. Pour construire la représentation multiéchelle d'un signal continu, on peut utiliser des noyaux gaussiens de variance s [Lin-

  test statistique d'hypothèses. Ces hypothèses sont élaborées à partir d'une représentation mathématique du mécanisme d'exocytose. Avant d'atteindre la membrane puis de libérer leur contenu, les vésicules transitent sur le cytosquelette. Certaines vésicules s'amarrent alors à la membrane et peuvent ensuite fusionner avec elle. Notons I la séquence d'images à analyser et

  le coefficient de diffusion pour cette vésicule et ∆ désigne l'opérateur laplacien 2D.

A

  partir de la carte de classification préliminaire � c , nous pouvons également détecter et localiser les comportements localement anormaux. Un mouvement est considéré comme anormal s'il diffère significativement du mouvement environnant. Nous introduisons une comparaison d'histogrammes des classes locales de mouvement de foule. En chaque point p, on mesure la distribution des huit classes de mouvement dans un voisinage η � p � de ce point. Nous comparons cette distribution empirique à celles de huit patchs voisins η i � p � répartis autour de p. Une mesure de distance entre histogrammes est calculée, prenant en compte la spécificité des classes. Les translations ont un caractère circulaire, tandis que les divergences, comme les rotations, s'opposent entre elles. Nous considérons qu'une anomalie, associée à un mouvement local, peut être caractérisée en chaque point p à travers la distance minimale entre l'histogramme des classes dans le patch central η � p � et ceux dans les patchs voisins η i � p � . Autrement dit, un mouvement local est anormal s'il ne ressemble à aucun mouvement situé à proximité. Le minimum de ces distances est supposé suivre une distribution des valeurs

  Joint research efforts in optics, electronics, chemistry and biology, among others, came to tremendous advances in observation technologies, cornerstone of modern biology research. Progress in microscopy imaging make possible the study of cell mechanisms in vivo and at sub-micron dimensions. Complex observation techniques can even go further in terms of resolution (e.g., super-resolution methods based on multiple exposures or electron microscopy), but, for now, they do not provide a fine temporal resolution, or do not even allow for acquiring temporal sequences or observing living cells. In contrast, conventional fluorescence microscopy features a very 1.1. C������ ��� ���������� fast acquisition rate, thanks to a continuous improvement of the sensitivity of digital camera sensors. This enables the use of fluorescent tags which themselves emit the light acquired by the microscope when excited. Researchers in genetics engineering developed protocols to tag the structures of interest, in order to locate and follow the latter in living cells, a breakthrough in cell biology. Modern cameras are so sensitive that sub-second exposure times are sufficient to detect and track fluorescent tags, and enable the fine recording of (intra)cellular dynamics. As an order of magnitude, in this thesis, we typically consider 100-nanometer structures traveling 10 microns per second, which represents a few pixels per frame at most. However, advances in microscopy bring new challenges, both in terms of biological queries and massive data handling. With the improvement of image resolution and with the advances in cell biology knowledge come more specific and complex questions. Expert visual analysis is no more sufficient to reliably answer new queries in biology, while human annotation is not fine enough to deal with the required accuracy. As an example, a classical way of characterizing the actors involved in a biological process is to inhibit the said actors one after the other and see if the process has changed. Changes can be either qualitative -something new or different has happened, or something does not occur anymore -or quantitative -some biological parameter value has changed. While visual inspection is sometimes sufficient to exhibit qualitative changes, parameter variations are generally impossible to assess quantitatively without the help of image processing and statistical tools. Moreover, dramatic improvement of spatial and temporal resolutions yields a so-fast increase of the amount of data to process, and computing methods must now cope with emerging big big data at the petabyte scale. The challenge of confronting biological models to such rich observations has yet to be faced. Hence, the need of automatic image processing tools in support of cell biology.For the three years of the thesis, we have collaborated with Jean Salamero's team at UMR 144 (CNRS, Institut Curie), a team interested in deciphering the dynamics of intracellular transports by the means of spatiotemporal imaging. In response to the need of digital tools, the purpose of the collaboration is to develop methods for the fully automatic analysis of the dynamics of biological structures observed by the means of optical microscopy imaging. To this end, dynamical models representing biological processes are proposed and evaluated on numerous microscope image sequences.

  ganized in three parts. Spot detection is handled in Part I, where we focus on fluorescence microscopy images. Then, proteins behaviors after vesicle fusion are modeled and estimated in Part II. Group motions are analyzed in Part III, spanning from microscopy to crowd analysis in videos. In Part I, an overview of existing spot detection methods is given in Chapter 3. The proposed two-step method is presented in the two following chapters: the scale selection is described in Chapter 4, and the spot segmentation in Chapter 5. Part II splits into three chapters. A presentation of the exocytosis process is given in Chapter 6, where we provide an overview of existing diffusion models in biology, and present the proposed fusion models. The exocytosis detection is presented in Chapter 7, and the issue of estimating biological parameters is tackled in Chapter 8. As for Part III, after introducing group motion analysis and related methods, the proposed labeled affine flow is presented in Chapter 9. and motion descriptors for the group motion characterization are introduced. Anomaly detection is addressed in Chapter 10. Thesis contributions are summarized and possible future works are proposed in Chapter 11. In particular, limits of the acquisition setups and methods designed in this thesis are pointed out, and improvements are proposed.

Chapter 3 .

 3 The landscape of spot detection 29 Before analyzing dynamics, objects of interest and space-time events must be detected. The (very active) domain of spot detection is introduced and existing spot detection methods are presented. The scale parameters of the existing detectors are outlined. C������ 1. G������ ������������ Chapter 4. Automatic selection of the detection scale 37 We propose a method for automatically selecting the optimal detection scale, which is related to the spot size. Efficiency of the scale selection is shown, and an extension to multiscale detection is proposed. Chapter 5. Spot segmentation 47 The spots are segmented by thresholding an enhanced image. A local adaptive method is proposed to cope with complex backgrounds. It is thoroughly compared to state-of-the-art spot detectors. Part II. Individual dynamics Chapter 6. Modeling of the vesicle fusion dynamics 71 The biological context of Part II is introduced. We describe the exocytosis process with a focus on vesicle fusion and protein diffusion. Classical diffusion models in biology are presented and discussed. The fusion dynamics of different proteins, namely Rab11, TfR and Langerin, are mathematically modeled. Chapter 7. Detection of diffusion dynamics 95 In order to further investigate biological parameters of the vesicle fusion, the latter must be detected. We propose detection methods for the three proteins of interest, which exploit the methodology presented in Part I. Chapter 8. Estimation of the biological parameters 109 After presenting existing methods for estimating diffusion in microscopy images, we propose an estimation method for the proposed vesicle fusion models.

  C������ 2. O� ��� ���� ��� ��� ����������� © Engraving by Robert Hooke[START_REF] Hooke | Micrographia: Or Some Physiological Descriptions of Minutes Bodies made by Magnifying Glasses with Observations and Inquiries thereupon[END_REF]

Figure 2

 2 Figure 2.1 -Sections of a cork bark showing cells

  Figure 2.2 -Some components of the cell

Figure 2

 2 Figure 2.5 -The cell cytoskeleton

  Figure 2.6 -Vesicle section showing the lipid bilayer

Figure 2

 2 Figure 2.7 -Principle of fluorescence microscopy

Figure 2

 2 Figure 2.8 -Principle of total internal reflection generating the evanescent wave

  Figure 3.1 -Cell images depicting particles of similar scale. (a,b) Tagged vesicles (bright spots) are of almost constant size over the image. Rab11 is tagged with mCherry in (a), (b) TfR is tagged with pHluorin in (b).

Figure 4

 4 Figure 4.1 -Scale-space LoG-transform L of a real TIRFM image I of M10 cell (Rab11-mCherry)

  ρ I (s ) (b) Blob density in the Gaussian noise image � g and in the input image I of (a)

Figure 4 . 2 -

 42 Figure 4.2 -Statistics on the number of blobs and behavior of the scale selection criteria. (a) Input TIRFM image. (b) Distribution of the number of blobs with respect to scale for the TIRFM image and a noise image. (c,f) Output of the four scale selection criteria. In all cases, the maximum is reached for s � ≈ 2 in this example.

  Figure 4.3 -Synthetic images with Gaussian spots corrupted by a Poisson noise

  Figure 5.1 -Segmentation maps obtained with global and local thresholding. (a) Gaussian spots are added to a varying background so that contrast increases from left to right. (b) With a global threshold, segmentation maps contain both false positives (red) and false negatives (yellow). (c) With a locally adapted threshold, far better performance is achieved.

  Figure 5.2 -Threshold maps with different windows. (a) Input image. (b) Use of a square window. (c) Use of a Gaussian window. Bright areas of the input image should correspond to low thresholds as the LoG is a negative filter. Strong blocky artifacts appear with the square window and background statistics are better estimated with the Gaussian window.

  Figure 5.2 demonstrates the advantage of such a smooth windowing compared to the square one. The bright background structures of the input image are well reflected in the threshold map of Figure 5.2c, while large errors are encountered in Figure 5.2b

  Figure 5.3 -Sample images from the Synthetic Data Generator benchmark for SNR = 2. Types are defined in the main text.

  Figure 5.4 -Case of images containing spots of two different sizes. (a) Input image with spots of two different sizes (resp. 4 and 16 px 2 ). (b) Two local maxima are located by C R at about 4 and 16 px 2 , corresponding to the spots variances. (c) Segmentation map when the LoG filter scale is set to the first local maximum in C R , that is ŝ1 = 4.29, then to the second local maximum ŝ2 = 15.41.

  Figure 5.5 -(a) Constructed image containing spots of variance 1.44 px 2 ; estimated PSNR: 30 dB. (b-f) Constructed images containing spots of variance 1 px 2 ; PSNR decreases from 30 dB (b) to 18 dB (f ).

  Figure 5.6 -Comparison of FROC curves obtained with ATLAS, MPHD and MS-VST on the proposed Spot in M10 dataset

( a )

 a Figure 5.9 -Comparison of FROC curves obtained with ATLAS, MPHD and MS-VST on the TIRFM simulator dataset

  Figure 5.11 -Segmentation results supplied by ATLAS on a real two-channel TIRFM image sequence. (a) First frame of the sequence. The two channels are displayed one above the other, Rab11-mCherry channel on top, Langerin-pHluorin below. (b) Segmentation of vesicles by ATLAS.

Figure 5

 5 Figure 5.12 -Segmentation of cores with ATLAS. (a) Input TMA 16-bit image (courtesy of Innopsys Company). The brightest core is 76 times brighter than the darkest one. Some spurious spots are 36 times brighter than the darkest spot. (b) Segmentation results supplied by ATLAS. Only one error is made (in red), by detecting a broken core of small spatial extent.

  Figure 5.13 -Comparison of segmentation results on a real 300-frame TIRFM sequence presenting elongated objects. (a) Two input TIRFM images (Rab11-mCherry) out of the 300frame sequence at time instants 93 and 300. (b) Segmentation results with MS-VST for the two frames of (a). A total of 24,481 objects are detected over the sequence. Some elongated objects of (a) are split (red) by MS-VST. (c) Segmentation results with ATLAS for the two frames of (a). A total of 24,195 objects are detected over the sequence. Few very close objects of (a) are merged (marked in yellow) by ATLAS, and the elongated objects are well recovered.

Figure 6

 6 Figure 6.1 -Main steps of exocytosis

Figure 6 . 2 -

 62 Figure 6.2 -Diffusion of a transmembrane protein (in green) in the plasma membrane

Figure 6

 6 Figure 6.3 -Break-down of the intensity model

  . Different models correspond to different applications, observed objects or microscopy techniques. Before proposing a new diffusion model for vesicle fusion, let us comment three widely-used diffusion models, namely the stationary model, the point source model and the fluorescence recovery after photobleaching (FRAP) model.

  Figure 6.4 -Stationary model at three time steps for C 0 = 1, D = 1, σ PSF = 1. Arrows represent the Dirac distribution.

Figure 6 . 5 -

 65 Figure 6.5 -Point source model at three time steps for C 0 = 1, D = 1, p 0 = 0, σ PSF = 1. Arrow represents the Dirac distribution.

Figure 6

 6 Figure 6.7 -Particle evolution meta-model for the Rab11 dynamics

  Figure 6.8 -Approximation of the error function

  Figure 6.9 -Comparison of real image sequences with simulations of the point source and SSED models

Figure 6 .

 6 Figure 6.10 -SSED model at three time steps for C 0 = 1; τ = 2. D = 1; p 0 = 0; σ PSF = 1. Arrows represent the Dirac distribution.

Figure 7 . 1 -

 71 Figure 7.1 -Illustration of the temporal block neighborhood

( a )

 a Figure 7.2 -Segmentation of the time varying vesicles

  Figure 7.3 -Classification results for several synthetic sequences. Top row: First frame of the sequence for different PSNRs. Wide spots are vesicles that have already diffused. Bottom row: Vesicles classified as diffusing are in red, vesicles classified as translating are in green. In all these examples, classification is performed with no error, except for the lowest noise level (a), where two vesicles are misclassified.

(Figure 7

 7 Figure 7.4 -FROC curves of the translation class for different noise levels. The curves are obtained by varying τ. Since only two classes are considered, the FROC curve for the diffusion class can be deduced from this one.

  Figure 7.5 -Examples of translating and diffusing vesicles in a real TIRFM sequence depicting a M10 cell. (a) First frame of the sequence. (b) Patch in successive frames depicting an elongated translating vesicle cluster. (c) Patches in successive frames depicting two diffusing vesicles.

  Figure 7.6 -Classification results for a real TIRFM sequence, whose estimated PSNR is 28.6. Results are displayed for a representative frame. The only classification error -framed in cyan -is a diffusion classified as translation. However, this vesicle has a very low intensity and changes its shape while diffusing. Two vesicles framed in green are detected as a single connected component. The diffusing vesicle framed in red corresponds to the upper one of Figure 7.5c.

( a )

 a Figure 7.7 -Temporal evolution of the region around a fusion event

  Figure 7.8 -Background subtraction in the neighborhood of membrane events detected in the TIRFM sequence of Figure 7.7

Figure 8

 8 Figure 8.1 -Comparison of the performance of several estimation methods for the point source model. The residuals log 10 � D -log 10 D are plotted, where D and � D are the true and estimated diffusion coefficients, respectively.

  Figure 8.3 -Results obtained for the SSED model on simulated sequences

  Figure 8.4 -Comparison of the histograms of the biophysical parameters � κ and � D estimated in 16 TIRFM image sequences depicting TfR or Langerin

  Figure 9.1 -Tissue sequence depicting cells whose plasma membranes are fluorescently labeled

  have focused on classifying structured group motions. The former determined motion patterns by clustering 4D flow vectors (2D position and velocity of points) in each frame according to proximity and similarity rules. The latter proposed to extract trajectories and accumulation points from the advection of flow fields over video sequences.

9. 2

 2 Labeled affine flow 9.2.1 Motion model candidates As aforementioned, to propose motion candidates, we only consider 2D parametric motion models. Specifically, at any point p = � x, y � ∈ Ω � I , the optical flow vector w � p � is approximated by an affine flow vector w θ � p � defined by:

( a )

 a Figure 9.2 -Overview of the motion estimation method applied to the Marathon bend sequence, where runners follow a 'U' from the upper left corner to the upper right corner. (a) Sample windows from W are plotted over the first frame of the sequence (red: large windows; green: medium ones; blue: small ones). (b) Map of the moving points detected in the first frame. (c) Motion models are estimated in each window, and candidates are pointwise selected (class color code in Table 9.2). (d) Optical flow deduced from the affine parameters of selected candidates. (e) Color code of (d).

  k � denotes the number of inliers. Let W � p � ⊂ W be the subset of windows containing a given point p, M � p � the set of motion model candidates for p, and Θ � p � the set of estimated parameters of the candidates. In our experiments, using the previously mentioned windows collection, 33 motion model candidates are available for each pixel (only 30 for pixels lying on the image borders).

  are further classified into eight motion classes, represented in Table 9.2 by colors, depending on the sign of the parameters and some combinations of them. Namely, the translation motion type is subdivided into four classes defined from view-based directions: Translation towards the top of the frame, or Northward translation (denoted TN); Westward translation (TW);

�

  2. The set of classes is denoted Γ = � TN, TW, TS, TE, SN, SP, RN, RP � , and we finally obtain the local classification map � c , as depicted in Figure 9.2c.In the meantime, we also obtain an estimate w of the velocity vector at every point p. The flow map will be denoted by F (t ). We term LAF the combination of Algorithm 3 Label affine flow estimation � Motion detection Determine the moving regions Ω � I � Motion estimationfor each W i ∈ W do Estimate the motion parameters θ i, k , k ∈ � T, S, R � � Model selection for each p ∈ Ω � I do for each θ i, k do Evaluate the variable ε i, k � p � .Evaluate the likelihood ψ � p, θ i,k � .Select the best model according to the AICc.the affine flow and motion class maps. The steps for obtaining it are summed up in Algorithm 3.

  From the computed occurrences, we build the signed histogram H � p, t � where bins corresponding to opposite classes (the pairs listed above) have opposite signs, that is, we implement the inhibition factor. As for translations, the inhibition could be added to the pair of North translation and South translation classes, and West translation and East translation classes, respectively. However, in contrast to scaling-and rotation-motion type classes, we rarely observe neighboring opposite directions of the translation-type classes. Adding

  second round of the regularization procedure consists in assigning the group motion class C � p � which is the most represented among the classes associated to the

  Figure 9.3 -Classification regularization with reinforcement and inhibition factors selected motion type K � p � : C � p � = arg max c ∈Γ s.t. k (c )=K (p)

  Figure 9.4 demonstrates the advantage of penalizing the complexity of the motion models and inhibiting local opposite motion classes. In that example, people do not run perfectly straight ahead from the left to the right of the image, resulting in local irrelevant rotation or scaling decisions.

  Figure 9.4 -Classification results on the Marathon lane sequence where runners move roughly toward the right ( ). Class color code is given in Table9.2.

  scale at which interesting behaviors (e.g., cell divisions) are found. As shown in Figure 9.6, the cells are evolving over time. Globally, the whole cell group is moving toward the lower left corner of the frame. This motion is well recovered as a predominant translation class, divided in both Southward ( ) and Westward ( ) translations.From the 60th frame to the end of the sequence, there are large deformations of the cell tissue on top of the tissue translation. Particularly visible is a large stretch of the tissue in the middle of the image (Figure9.6d,e), which is well recovered as a group divergence ( ). Conversely, in the upper part, contracting cells are recovered ( ). The large deformation also induces other types of motion around the central stretching.

Figure 9

 9 Figure 9.6 -Classification results on the Tissue sequence of Figure 9.1. (d) The arrow represents the approximate stretch direction. (h) Regions where cell divisions occur are framed in red.

Finally, we assign

  to β the local path label π i ∈ Π which is the closest to the orientation of vector w � β, t � . Under our stationarity assumption, the local path labels do not depend on time, which will be verified below. Indeed, the eight-quantization of local paths smooths out the slight variations in time of w � β, t � . The local path then expresses the space-time stationary information on the dominant orientation of the local displacement observed in the block β.

( a )

 a Figure 9.7 -Dominant path recovery on the Marathon bend sequence (color code is given in Figure 9.2e). (a) The velocity flow is deduced from the selected motion models. (b) The image is divided into small blocks, where the mean velocity vector in space and time is computed. (c) The global paths are recovered from P l and the dominant path is outlined in red.

  Figure 9.8 -Dominant path recovery (red) in several examples. (a-c) The Obstacles sequence, where pedestrians slalom between walls. (d-f) The Marathon lane sequence. (g-i) Cells of the Tissue sequence. (j-l) Fishes swirling in the Shoal sequence.

  of the blocks β. The parameter sensitivity is low in practice, except for the regularization size r , which can be set to attain a user-selected classification granularity.

Figure 10

 10 Figure 10.1 -Normal scene and variety of local anomaly cases (framed in red) in group motion scenes. (a) A normal situation where people run together. (b) A man is walking against the crowd; abnormal behavior is localized on him. (c) The cyclist crossing a group of pedestrians is the anomaly. (d) People starting a "circle pit" during a music festival form the abnormal behavior. (e) The car driving into oncoming traffic is the anomaly, along with the one maneuvering to avoid it.

Figure

  Figure 10.2 -Anomaly detection on the Wrong way sequence, where a man is walking to the bottom of the image (south translation ) against a crowd moving in the opposite direction (north translation ). (a) The central patch ν � p � is plotted in red at p = � 190, 140 � , along with

�

  defined in (10.1). The minimum distance ζ � p � for all the points p ∈

a

  dense moving crowd (Figure10.3a), and we have added a simulated anomaly. Specifically, we have cropped a rectangular small region from the synthetic sequence Corridor where pedestrians are walking from left to right. After mirroring it, we have inserted it again in the original sequence but with a translation from right to left. Thus, one pedestrian is walking in the opposite direction of the whole group and constitutes an anomaly in the resulting Modified corridor sequence. The anomaly is correctly detected (Figure10.3c). Conversely, very few false detections occur in the original Corridor sequence without anomaly. More precisely, as explained later, the FPR is indeed really close to the user-selected PFA.

Figure 10

 10 Figure 10.3 -Anomaly detection results and impact of the size of η � p � and the size (m ×m ×3) of

Figure 10

 10 Figure 10.4 -Temporal analysis of anomaly localization over the synthetic Escape sequence. (ab) People are leaving the room toward the right through a single door. (c) Map of the number of times an anomaly is detected at each point over 400 frames with PFA = 2% (color scale is given on the right).

  Figure 10.5 -Results on the Wrong way sequence. (a,b) Detection results on the first and last labeled frames, with m = 61 and PFA = 1% (green: TP; blue: FN; yellow: FP). (c) ROC curves for m = 41 and m = 61. Area under curve is equal to 0.97 for m = 61.

  Figure 10.1b, Figure 10.2d, and Figure 10.5a,b, the man heading to the bottom of the frame is well segmented. ROC curves established for the Wrong way video are very good. Area under the curve for m = 61 is equal to 0.97 and equal to 0.96 for m = 41. The camera being closer to the crowd than in previously processed sequences, pedestrians look bigger.

Figure 10 .

 10 Figure 10.6a,b. The area under ROC curve is 0.93 with a patch size of 21 pixels, similar results were obtained for a size of 31 pixels (Figure 10.6c).Supplementary results obtained on other real sequences involving various behaviors (presented in Figure10.1) are reported in Figure10.7. They confirm the ability of

  Figure 10.6 -Results on a sequence of the PETS'09 dataset. (a,b) Top row: Input frames (the red arrow points to the pedestrian to be detected). Bottom row: Anomaly maps with m = 31 and PFA = 1% (green: TP; blue: FN; yellow: FP). (c) ROC curves for m = 21 and m = 31 (area under curve is equal to 0.93 for both).

Figure 10

 10 Figure 10.7 -Anomaly detection for the samples presented in Figure 10.1. Top row: input images. Bottom row: Anomaly maps computed with our method with the same parameter setting (PFA = 1% and m = 61) apart from (e) where m = 81. (a) No anomaly is detected in this normal situation. (b) The man who is walking against the crowd, (c) the cyclist crossing a group of pedestrians, (d) people running in the crowd, and (e) the car driving into oncoming traffic and the one maneuvering to avoid it are correctly detected.

Figure 10 .

 10 Figure 10.8 contains plots expressing the relationship between PFA and FPR for different sequences (synthetic and real ones, without and with anomalies) and for two different sizes of patches ν � p � . Corridor and Modified corridor were introduced in Figure 10.3a,b. In the Marathon lane sequence, people run from left to right, and nothing anomalous happens (one image of the sequence is given in Figure 10.1a).

  Détection de spots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Exocytosis process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Existing diffusion models . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Proposed fusion models for Rab11 . . . . . . . . . . . . . . . . . . . . . 6.4 Proposed SSED model for TfR and Langerin . . . . . . . . . . . . . . . . 6.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Detection based on Rab11 dynamics classification . . . . . . . . . . . . 7.2 Event detection using fluorophore pH-sensitivity . . . . . . . . . . . . .

	xx II Individual dynamics 6 Modeling of the vesicle fusion dynamics 6.1 7 Detection of diffusion dynamics Bibliography List of publications Acronyms D�������� ��� ���������������� �� ��������� ������ �� ����� ��������� 3 C������� xxi 7.1 C������� 1 General introduction

viii Partie II. Analyse des dynamiques individuelles . . . . . . . . . . . . . . . . . xi Partie III. Analyse des dynamiques de groupes . . . . . . . . . . . . . . . . . xiii 1 General introduction 3 1.1 Context and motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Objectives and contributions . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3 Organization of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2 On the cell and its observation 13 2.1 A bit of history . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Organization and processes of the cell . . . . . . . . . . . . . . . . . . . 16 2.3 TIRFM for imaging membrane dynamics . . . . . . . . . . . . . . . . . 19 I Spot detection 27 3 The landscape of spot detection 29 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4 Automatic selection of the detection scale 37 4.1 Scale-space representations . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.2 Blob distribution and scale selection . . . . . . . . . . . . . . . . . . . . 40 5 Spot segmentation 47 5.1 Adaptive thresholding . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.2 Estimation window . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.3 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 7.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Estimation of the biological parameters 109 8.1 Existing diffusion estimation methods . . . . . . . . . . . . . . . . . . . 8.2 Proposed estimation method for the translation . . . . . . . . . . . . . 8.3 Proposed estimation method for the point source model . . . . . . . . 8.4 Proposed estimation method for the SSED model . . . . . . . . . . . . . 8.5 Discussion and perspectives . . . . . . . . . . . . . . . . . . . . . . . . . III Group dynamics 9 Labeled affine flow for group motion characterization 127 9.1 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.2 Labeled affine flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.3 Motion patterns recognition . . . . . . . . . . . . . . . . . . . . . . . . . 9.4 Recovery of principal paths . . . . . . . . . . . . . . . . . . . . . . . . . 9.5 Discussion and perspectives . . . . . . . . . . . . . . . . . . . . . . . . . 10 Anomaly detection and localization 147 10.1 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.2 Histogram distance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.3 Pixel-level anomaly detection and localization . . . . . . . . . . . . . . 10.4 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 General conclusion 163 11.1 Main contributions and results . . . . . . . . . . . . . . . . . . . . . . . 11.2 Ongoing and future work . . . . . . . . . . . . . . . . . . . . . . . . . . .

  more easily and more accurately than global modeling, for which the complexity rapidly explodes if one wants to finely model the entire image. Indeed inevitable stationarity assumptions of global approaches are much reduced in local approaches. Moreover, the local characteristics of the approach is also true in the temporal dimension when we process image sequences. Such instantaneous methods,

	using only a couple of consecutive frames to estimate the dynamics, offer a lot of ad-
	vantages, including the possibility of ensuring real-time computation if needed. Local
	approaches also allow neglecting low frequency variations such as spatially varying
	background contrast or, in fluorescence microscopy, temporal fading known as pho-
	tobleaching.

Secondly, we rely on local approaches which have the advantage of being computationally efficient. Their complexity being typically proportional to the image sequence data volume, they are particularly recommendable in a context of increasingly big data, for which global approaches cannot provide reasonable computation time. Local modeling handles spatially varying im-age statistics much

Table 2

 2

	.1 -Biological and microscope orders of magnitude
	Description		Dimension Weight*	Duration
	Cell Nucleus Vesicle Actin filament (diameter) (diameter) (diameter) (diameter)	≈ 10 2 µm ≈ 10 1 µm ≈ 10 2 nm 7 to 8 nm	≈ 1 ng
	Microtubule (diameter) 15 to 25 nm
	Membrane (thickness)	4 to 5 nm
	Membrane phospholipid			< 1 kDa
	Rab11			24 kDa
	Langerin			37 kDa
	TfR			85 kDa
	mCherry			29 kDa
	pHluorin			27 kDa
	Spatial optical resolution		200 nm
	Physical pixel side		16 µm
	Pixel side in the image		160 nm
	Penetration depth Exposure time		≈ 10 2 nm	100 ms
	Frame rate			10 f/s
	Sequence length			30 to 120 s
	1 kDa ≈ 1.66 × 10 -21 g *Source: PhosphoSitePlus (www.phosphosite.org)
	be the frame period (f ), or frame in short, in the remaining of the thesis, such that I (0)
	is the first frame of the sequence and I (t + 1) is the frame following that of time t .
	The photosites collect photons and convert them into electrons. Photosites are ar-
	ranged in a grid, and are represented by their Cartesian coordinates x and y . The dis-
	tance between horizontally and vertically adjacent pixels is the pixel pitch, denoted
	px, which is the natural distance unit in the image. This way, pixels have integer co-
	ordinates. The set of pixel locations, called domain, is denoted Ω I and is therefore a
	subset of Z 2 .
	The intensity of a pixel p = tion of the number of photons n � x, y � p, t �	∈ Ω I in frame I (t ) is denoted I � received by the associated photosite from time � p, t � . It is a func-
	t to time t + ∆t , where ∆t is the exposure time. In the images we deal with, there is
	a negligible time interval between consecutive exposures, so ∆t is approximately the
	frame period: ∆t ≈ 1 frame.

Table 3 .

 3 1 -Main parameters of state-of-the-art spot detection methods (see main text for acronym meaning)

	Method	Main parameters
	TH
	[Bright and Steel, 1987]

Table 4

 4

		Spot		Selection criterion
	SNR variance SLT-LoG	C D	C R	C P	C KL
	1	4.00	3.00	2.99	6.19	1.00	2.99
	2	4.00	3.00	3.58	5.16	1.00	3.58
	3	4.00	3.00	4.30	4.30	4.30	4.30
	4	4.00	3.00	4.30	4.30	4.30	4.30
	5	4.00	3.00	4.30	4.30	4.30	4.30
	5	1.00	3.00	1.44	1.44	1.73	1.44
	5	2.25	3.00	2.49	2.99	2.49	2.49
	5	4.00	3.00	4.30	4.30	4.30	4.30
	5	9.00	9.00	8.92	8.92	8.92	8.92
	5	16.00	9.00 15.41 15.41 15.41 15.41
	5	25.00	27.00 22.19 22.19 22.19 22.19
	5	36.00	27.00 31.95 31.95 31.95 31.95

.1 -Scale selected by five criteria on synthetic images

  from local image statistics. To this end, we assume that the distribution of the image background is smooth and corrupted by white Gaussian noise. It holds because low frequency background structures are locally constant if the neighborhood is small enough, while noise is supposed to be normally distributed. Then, H is obtained by finite convolution of I , so that this assumption also holds for H . For every point p ∈ Ω I , the local mean µ H and variance σ 2 H are estimated over a window W p centered in p.

	Then, we can infer the likelihood ψ	�	p	� of the background model N	�	µ H , σ H	�	given

H � p, s � � :

Table 5 .

 5 1 -Comparison of TPR* values with state-of-the-art methods on the Synthetic Data Generator dataset for SNR = 2

	Object shape		Round			Elongated
	Background type	A	B	C	A	B	C
			TPR*s of supervised detection methods
	AB	[Jiang et al., 2007] 0.99 0.94 0.94 0.99 0.99 0.99
	FDA	[McLachlan, 1992] 0.99 0.99 0.96 0.99 0.99 0.99
		TPR*s of unsupervised detection methods
	TH	[Bright and Steel, 1987] 0.99 0.88 0.48 0.99 0.96 0.56
	HD	[Smal et al., 2008] 0.99 0.97 0.90 0.99 0.99 0.97
	WMP	[Olivo-Marin, 2002] 0.81 0.37 0.30 0.31 0.17 0.18
	IFD1	[ter Haar Romeny, 2003] 0.98 0.67 0.89 0.53 0.31 0.31
	IFD2	[ter Haar Romeny, 2003] 0.99 0.46 0.71 0.59 0.23 0.19
	MTH		[Soille, 2003] 0.99 0.87 0.88 0.99 0.98 0.91
	SEF	[Sage et al., 2005] 0.99 0.91 0.95 0.99 0.99 0.95
	MS-VST	[Zhang et al., 2007] 0.99 0.99 0.93 0.99 0.99 0.96
	MPHD [Rezatofighi et al., 2012] 0.99 0.97 0.94 0.99 0.99 0.95
	SLT-LoG	[Basset et al., 2014a] 1.00 0.99 0.98 1.00 1.00 1.00
	ATLAS		1.00 1.00 0.99	1 †	1 †	1.00

† 1 means that every spot is correctly detected, while 1.00 corresponds to a rounded value A Poisson noise is added to obtain a SNR of 2. Table

5

.1 summarizes the TPR* values obtained by the detection methods on this benchmark. Number 1.00 corresponds to

Table 5 .

 5 2 -Comparison with LR-MRF on the Subcell dataset

		Prec	Rec	F-score
	LR-MRF 0.876 0.858	0.867
	ATLAS	0.994 0.957	0.969

Table 5 .

 5 3 -Comparison with MS-VST and MPHD on the Spot in M10 dataset

	Spot variance 1.44	1	1	1	1	1
	PSNR	30	30	25	23	21	18
				TPR at FPR* = 0.01
	MS-VST	0.69 0.71 0.46 0.40 0.29 0.13
	MPHD	0.64 0.66 0.50 0.44 0.35 0.18
	ATLAS	0.70 0.71 0.50 0.45 0.34 0.18
			Area under FROC curve
	MS-VST	0.79 0.80 0.65 0.60 0.52 0.36
	MPHD	0.77 0.78 0.66 0.61 0.53 0.39
	ATLAS	0.80 0.80 0.69 0.64 0.55 0.44
			Maximum F-score F*
	MS-VST	0.82 0.83 0.68 0.64 0.56 0.41
	MPHD	0.80 0.81 0.70 0.66 0.58 0.44
	ATLAS	0.83 0.83 0.71 0.67 0.58 0.48

one hundred images and more than six thousand spots. The dataset is available at https://serpico.rennes.inria.fr.

FROC curves are plotted in Figure

5

.6 for each image collection. MS-VST performed better than MPHD at high PSNR, while the opposite holds for noisier images. Our method ATLAS is on par or slightly better than MS-VST at high PSNR, and demonstrates better robustness to noise. Therefore, FROC curves obtained with AT-LAS nearly always dominate others. Table

5

.3 summarizes TPR* scores, areas under FROC curves and best F-scores F* of the three methods. Our method always ranks

Table 6

 6 C 0 : initial local concentration; D : diffusion coefficient; σ PSF : radius of the PSF; p beam : center of the laser beam; σ beam : radius of the laser beam; p 0 : position of the vesicle; τ: residence time.

	.1 -Characterization of the diffusion models
	Source distribution
	Model Spatial distribution	Temporal distribution Parameters*
	Stationary Uniform	Dirac	C 0 , D , σ PSF
	Point source Dirac	Dirac	C 0 , D , p 0 , σ PSF
	FRAP Gaussian hole	Dirac	C 0 , D , p beam , σ beam
	SSED Arbitrary	Exponential decay	C 0 , τ, D , p 0 , σ PSF
	(with small support)

*

  Detection and estimation procedure for the point source model � Fusion event detection Select detection scale s � in I Detect spots e i in I δt at scale s � Estimate model parameters p 0i , A i , σ i and D i in �

	hypotheses are: Estimate foreground � z i � Non-diffusing event discarding for each t ∈ � t 0i , t 0i + 4 � do Fit a Gaussian spot model to � z i (t ) if � σ 2 i (t ) � t �t 0i and A i (t ) are not correlated then Discard e i • H for each e i do � Fusion model estimation	z i	� t 0i	. The competing

0 : σ 2 i and A i are not correlated; • H 1 : σ 2 i and A i are correlated.

Algorithm 1

Table 7

 7

		.1 -Diffusive event detection performance
	α	N	#deb	#diff	#n-diff #miss #err
	1% 90 6/90 13/34 49/50	21/34 1/50
	5% 90 6/90 18/34 48/50	16/34 2/50
	10% 90 6/90 20/34 47/50	14/34 3/50

Table 8 .

 8 1 -Requirements of the diffusion estimation methods

	Method	Observed quantity	Source distribution
	Particle tracking	Particle trajectory	Sparse
	Correlation fitting Local concentration Uniform
	Intensity fitting	Local concentration Known

  Algorithm 2 Detection and estimation procedure for the SSED model for each e i do Estimate foreground z Estimate Gaussian spot parameters p 0i , A i and σ i in � Moreover, as the intensity model at t = t 0 is a Gaussian spot, we can reliably estimate p 0 , A 0 and σ 0 by fitting a Gaussian spot model in frame I �

	Estimate κ i and D i with various θ (init) i Retain best fit parameters	z i in �	z i	� t 0i	�
	if goodness-of-fit is too low then
	Discard e i
					t 0

Table 9 .

 9 1 -Group motion classes compared to behaviors from[START_REF] Solmaz | Identifying behaviors in crowded scenes using stability analysis for dynamical systems[END_REF]

	Motion	Motion class	Behavior from
	type	direction	[Solmaz et al., 2012]
		North
	T	West South	Lane
		East
	S	Convergence Divergence	Bottleneck Fountainhead
	R	Clockwise Counterclockwise	Ring
	3. Pointwise selection of the optimal motion models.

Table 9

 9

			.2 -Motion classes definition
	Motion Motion Class	Class
	type	class	color	direction	Criterion
	T

Table 9 .

 9 3.

Table 9

 9

			.3 -Classification notation
			Regularized
	Motion type Motion class	Local map � k � c	map K C	Possible values S TN TW TS TE SN SP RN RP T R

  .2) where U , resp. U i , designates the set of translation class values in the patch ν � � � , resp. u i (l ), represents any value of the set U , resp. U i , once recorded as � 0, 1, 2, 3 � for � TN, TW, TS, TE � . n is the number of elements of both U and U i . In case that the numbers of elements of sets U and U i are different, a normalization step is added

	in the patch ν i	�	p	� , resp.

p � . u � l

  mean absolute logarithmic error. 106 ML maximum likelihood. 133 MPG mixed Poisson-Gaussian. 22, 23 MPHD maximum possible height dome. 31, 46, 52, 54, 55, 59 MSD mean square displacement. 100 MS-VST multiscale variance stabilizing transform. 23, 30, 46, 47, 52, 55, 58-60 MTH morphological top hat filter. 30 PDE partial differential equation. 74 PDF probability density function. 72, 74 PFA probability of false alarm.

serpico.rennes.inria.fr

mobyle-serpico.rennes.inria.fr

www.irisa.fr/vista/Motion2D

The Bessel functions are historically denoted I n but we want to avoid confusions with the image sequence I

Probabilistic formulation of lateral diffusion was derived with Vincent Briane during his PhD with the Serpico team at Inria.

Synthetic sequences were generated by Jérôme Boulanger at UMR 144.

Correlation methods were investigated with Anca Caranfil during her Master internship in the Serpico team at Inria Rennes.

, thus breaking the stationary assumption for the correlation. The same argument holds for spatiotemporal image correlation spectroscopy (STICS) methods, which include a spatial component in the correlation function (8.3)[START_REF] Hebert | Spatiotemporal image correlation spectroscopy (STICS) theory, verification, and application to protein velocity mapping in living CHO cells[END_REF], Kolin and Wiseman, 2007, Di Rienzo et al., 2014].

Correlation fitting was evaluated with Anca Caranfil during her Master internship in the Serpico team at Inria Rennes.

Motion2D: http://www.irisa.fr/vista/Motion2D/

http://www.cvg.rdg.ac.uk/PETS2009

Remerciements

General conclusion

T �� ������ aimed at detecting, estimating and characterizing dynamical events in image sequences. We have focused on two types of dynamics, namely individual and group dynamics. In each case, we have proposed automatic methods to calibrate algorithm parameters by using statistical arguments. Local approaches allowed us to design computationally light and efficient methods able to cope with complex image sequences. We have mainly been concerned with fluorescence microscopy image sequences and videos of crowded scenes.

Within the scope of fluorescence microscopy, we have principally focused on the exocytosis, a dynamical process by which some molecules are transported from the cell interior to the extracellular medium. Specifically, we have studied the behavior of three proteins, Transferrin receptor (TfR), Langerin and Rab11, during the late steps of exocytosis. This includes transporting the molecules and passing through the cell frontier (plasma membrane). Molecule transport is mediated by the so-called vesicles, which enclose the molecules and then literally fuse to the plasma membrane to free molecules outside the cell.

By using total internal reflection fluorescence microscopy (TIRFM), we have studied in particular the vesicle fusion dynamics. In this analysis workflow, vesicles and fusion events, appearing as bright spots in the images, must first be detected. Therefore, we have developed a new spot detection method ATLAS based on adaptive thresholding of LoG images with auto-selected scale. After detecting the fusion events, we have proposed models and methods to estimate the motion of the proteins at the vesicle fusion.

We have finally investigated methods for the study of dynamics involving interactions between individuals. We have introduced a motion measurement data, labeled affine flow (LAF), which involves both the pointwise velocity vectors and the local motion type, using only two consecutive frames. This hybrid data allowed us to target