Tomas & Pati

Juliette Baptiste

Keywords: Géométrie Stochastique, Modèles à Base d'Énergies, Détection, Objets Multiples, Images Satellitaires, Très Haute Résolution. viii Learning Stochastic geometry models and convolutional neural networks. Application Stochastic Geometry, Energy Based Models, Detection, Multiple objects, Satellite Images, Very High Resolution. viii Stochastic Geometry, Energy Based Models, Detection, Multiple objects, Satellite Images, Very High Resolution

Apprentissage de modèles de géométrie stochastique et réseaux de neurones convolutifs. Application à la détection d'objets multiples dans des jeux de données aérospatiales. Résumé Les drones et les satellites en orbite basse, dont les CubeSats, sont de plus en plus utilisés pour la surveillance, générant d'importantes masses de données à traiter. L'acquisition d'images satellitaires est sujette aux perturbations atmosphériques, aux occlusions et à une résolution limitée. Pour détecter de petits objets, l'information visuelle est limitée. Cependant, les objets d'intérêt (comme les petits véhicules) ne sont pas uniformément répartis dans l'image, présentant des configurations spécifiques. Ces dernières années, les Réseaux de Neurones Convolutifs (CNN) ont montré des compétences remarquables pour extraire des informations, en particulier les textures. Cependant, modéliser les interactions entre objets nécessite une complexité accrue. Les CNN considèrent généralement les interactions lors d'une étape de post-traitement. En revanche, les Processus Ponctuels permettent de modéliser la vraisemblance des points par rapport à l'image et leurs interactions simultanément. La plupart des modèles stochastiques utilisent des mesures de contraste pour la correspondance à l'image ; elles sont adaptées aux objets à contraste fort et faible complexité du fond. Cependant, les petits véhicules sur les images satellitaires présentent divers niveaux de contraste et une grande variété d'objets de fond et de fausses alarmes. Cette thèse de doctorat propose d'utiliser les CNN pour l'extraction d'informations, combinées aux Processus Ponctuels pour modéliser les interactions, en utilisant les sorties CNN comme données. De plus, nous introduisons une méthode unifiée pour estimer les paramètres du modèle de Processus Ponctuel. Nos résultats montrent l'efficacité de ce modèle sur plusieurs jeux de données de télédétection, avec régularisation géométrique et robustesse accrue pour un nombre limité de paramètres.

Generic

A[i, j, k], A[ρ]

For a tensor A, the value at coordinates (i, j, k) (interpolated if not integer values). By extension indexing with pixel ρ returns the value at location (ρ i , ρ j ).

Introduction

La télédétection, dont les racines remontent jusqu'à la photographie en ballon à la fin du XIXe siècle, a considérablement évolué depuis, pour devenir essentielle pour la recherche, les applications commerciales ou la défense. Au fil des années, les capteurs sont devenus plus sophistiqués, cette tendance a été en partie motivée par l'escalade du renseignement militaire dans le contexte de la Seconde Guerre mondiale, puis de la guerre froide. Cela s'ajoute à un besoin accru de surveiller notre environnement à l'échelle mondiale, à mesure que les effets de l'humanité sur celui-ci deviennent plus préoccupants. Ce faisant, la quantité croissante de données produites par ces capteurs exige des solutions informatiques avancées pour le traitement. Les progrès de la télédétection correspondent à l'essor de la vision par ordinateur : depuis ses débuts dans les années 1950, en passant par le Perceptron en 1958 [START_REF] Rosenblatt | The perceptron: A probabilistic model for information storage and organization in the brain[END_REF] et la transformée de Hough en 1972 [START_REF] Duda | Use of the Hough transformation to detect lines and curves in pictures[END_REF] jusqu'à la résurgence des réseaux de neurones à la fin des années 1980 et l'essor des réseaux de neurones à convolution (CNN), enfin plus récemment l'adoption du Deep Learning avec AlexNet [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] ou RCNN [START_REF] Girshick | Rich Feature Hierarchies for Accurate Object Detection and Semantic Segmentation[END_REF] par exemple. À première vue, la télédétection peut être comprise comme "la collecte d'informations à distance" [START_REF] Campbell | Introduction to Remote Sensing, Fifth Edition[END_REF]. Cependant, cette définition trop large engloberait la microscopie jusqu'à l'enregistrement audio. Le domaine de la télédétection se concentre plutôt sur "l'observation des surfaces terrestres et aquatiques de la Terre au moyen de l'énergie électromagnétique réfléchie ou émise" [START_REF] Campbell | Introduction to Remote Sensing, Fifth Edition[END_REF]. En pratique, les données considérées sont constituées d'images optiques (ou proches-optiques comme l'infrarouge), ainsi que d'images SAR (Synthetic Aperture Radar, ou Radar à Synthèse d'Ouverture, RSO, en français) ou Lidar (Laser Imaging, Detection, and Ranging) prises depuis des drones, des avions, des ballons ou des satellites.

Les images satellitaires contiennent de nombreux artefacts et occultations inhérentes à leur acquisition (voir [START_REF] Tupin | Remote sensing imagery[END_REF]), ce qui les distingue des images usuelles ; c'est-à-dire des photographies que nous rencontrons quotidiennement, prises à la hauteur d'oeil, généralement avec un sujet proche du centre, où les objets d'intérêt représentent une proportion importante de l'image. Compte tenu de la faible résolution spatiale des images satellitaires (autour de 0, 5 m par pixel), les objets d'intérêt ne mesurent parfois que quelques pixels. Les informations visuelles étant limitées, les opérateurs humains ont souvent recours à leur connaissance a priori de l'objet d'intérêt lors de l'analyse d'une image. De la même façon, nous visons à construire des modèles qui prennent en compte cet a priori.

Le but de cette thèse est d'étudier la combinaison d'approches modernes de réseaux de neurones convolutifs (CNN) avec des modèles géométriques utilisant des a priori ou des contraintes. Bien que les modèles CNN se soient révélés très efficaces pour extraire des informations sur la texture, ils ne parviennent pas à prendre en compte les interactions à plus longue portée entre les objets sans augmenter considérablement la complexité et le nombre de paramètres. D'autre part, les modèles de Processus Ponctuels (PP) proposent de modéliser des configurations d'objets géo-6 métriques comme un processus stochastique : premièrement, la construction du Processus Ponctuel elle-même définit le type de géométries extraites. De plus, le modèle stochastique contient des a priori qui peuvent guider le modèle vers des configurations plus probables a priori.

Modèle d'observation de l'image. L'approche par Processus Ponctuel dérive de modèles stochastiques d'analyse d'images. En considérant l'image comme un ensemble de pixels dans S d ⊂ Z 2 , nous pouvons la modéliser comme une variable aléatoire X dans l'espace de probabilité (Ω, A, P), où les valeurs en niveaux de gris correspondent aux réalisations de la fonction X :

X : Ω → R S d .
(0.0.1)

Si l'on considère l'image comme une observation bruitée d'un phénomène sous-jacent Φ, on peut utiliser le modèle bayésien suivant : Processus Ponctuels. Les Processus Ponctuels permettent de modéliser de tels phénomènes Φ, comme un ensemble d'objets paramétriques représentant des objets géométriques tels que des segments, des cercles, des rectangles, des polygones, etc. Nous revenons plus en détail sur les Processus Ponctuels dans le chapitre 3. Ici, nous considérons les configurations y comme un ensemble {y 1 , . . . , y n(y) } de n(y)éléments y qui représentent nos objets d'intérêt. Dans le cadre bayésien, un modèle d'observation peut être construit comme :

P(X = X) =
P(X ρ = X ρ |Y = y) = 1(ρ ∈ y)P(X ρ = X ρ |ρ ∈ y) + 1(ρ / ∈ y)P(X ρ = X ρ |ρ / ∈ y). (0.0.5)
Dans cette approche, on applique un modèle d'objet (ρ ∈ y) ou d'arrière-plan (ρ / ∈ y) si un pixel ρ appartient à la silhouette de y ou pas. Un modèle courant si la silhouette correspond à des pixels clairs sur un fond plus sombre consiste à modéliser le premier plan P(X ρ = X ρ |ρ ∈ y) comme une loi gaussienne de grande valeur. L'arrière-plan P(X ρ = X ρ |ρ / ∈ y) se voit ensuite assigné une valeur inférieure [START_REF] Baddeley | Stochastic geometry models in high-level vision[END_REF][START_REF] Perrin | Tree crown extraction using marked point processes[END_REF].

Cependant, cette approche bayésienne présente certaines limites (Ben Hadj, Chatelain, Descombes, & Zerubia, 2010) :

• Les arrière-plans non homogènes sont plus difficiles, voire impossibles à modéliser. Le modèle gaussien simple pour P(X ρ = X ρ |ρ / ∈ y) n'est plus utilisable.

• La vraisemblance des pixels de la classe de premier plan ne dépend pas des propriétés morphologiques de la classe à extraire. Par exemple, le modèle peut essayer d'insérer autant de petits objets que possible dans une grande silhouette [START_REF] Craciun | Joint Detection and Tracking of Moving Objects Using Spatio-temporal Marked Point Processes[END_REF].

Pour ces raisons, [START_REF] Ben Hadj | Estimation des paramètres de modèles de processus ponctuels marqués pour l'extraction d'objets en imagerie spatiale et aérienne haute résolution[END_REF] propose un modèle détecteur fondé sur une mesure de contraste entre l'intérieur et l'extérieur de l'objet. Néanmoins, ces mesures de contraste reposent sur le contraste élevé de l'objet par rapport à son arrière-plan. Cette approche atteint ses limites lorsque le fond présente des éléments contrastés qui ne devraient pas être détectés. Nous montrons un exemple pratique de ces limites dans la partie 4.2.1.1.

Réseaux de Neurones Convolutifs. Alors que le modèle d'image ci-dessus était motivé en partie par les limitations de calcul, la disponibilité croissante de la puissance de calcul a facilité l'essor des réseaux de neurones convolutifs (CNN, Convolutional Neural Network). Bien qu'introduit pour la classification d'images en 1989 [START_REF] Lecun | Backpropagation Applied to Handwritten Zip Code Recognition[END_REF], ce n'est que dans les années 2010 que CNN a été largement adopté avec AlexNet [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] pour la classification ou RCNN [START_REF] Girshick | Rich Feature Hierarchies for Accurate Object Detection and Semantic Segmentation[END_REF] pour la détection d'objets.

Les réseaux de neurones convolutifs traitent les données d'une topologie de type grille via une séquence d'opérations de convolution et de regroupement (agrégation spatiale). Les modèles convolutifs profonds apprennent les filtres de convolution pour une tâche spécifique à partir de grandes bases de données d'images telles que ImageNet [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] par exemple. Un modèle de segmentation tel qu'Unet [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] effectue une classification par pixel. Ceux-ci produisent des cartes à valeurs réelles dans la plage [0, 1], qui sont interprétées comme des probabilités de classe. Par exemple, la probabilité que le pixel ρ appartienne à la classe C 1 (parmi les classes N C ) est donnée par :

 P (ρ ∈ C 1 |X) = exp (F (X)[ρ, C 1 ]) ∑ N C k=1 exp (F (X)[ρ, C k ])
, (0.0.6) avec F (X) est le tenseur de sortie du CNN appartenant à R S d × 1,N C , et N C le nombre de classes (2 dans le cas de la classification objet ou arrière-plan). Ce score ˆ P (ρ ∈ C 1 |X) mesure la présence ou non d'un objet, avec un modèle convolutif "léger" (faible nombre de paramètres) dont les noyaux sont appris à partir des données.

Contributions

Le point de départ de notre travail -inspiré de (T. [START_REF] Li | Feature Extraction and Tracking of CNN Segmentations for Improved Road Detection from Satellite Imagery[END_REF]) -est d'exploiter la partition ˆ P (ρ ∈ C 1 |X) et de l'utiliser pour construire l'énergie externe du Processus Ponctuel (attache aux données), remplaçant la mesure de contraste. Dans cette thèse de doctorat, nous proposons l'incorporation de modèles d'interaction dans les méthodes de détection d'objets, tout en tirant parti des capacités des réseaux de neurones convolutifs profonds. Les images satellitaires sont intrinsèquement imparfaites en raison des perturbations atmosphériques, des occultations partielles et de la résolution spatiale limitée. Pour compenser ce manque d'informations visuelles, il devient essentiel d'incorporer des connaissances préalables sur la disposition des objets d'intérêt.

D'une part, les méthodes basées sur des CNN sont excellentes pour extraire des motifs dans les images, mais ont du mal à apprendre des modèles d'interaction objet à objet sans avoir à y introduire des mécanismes d'attention tels que les Transformers qui augmentent considérablement la complexité du modèle : par exemple, certaines approches proposent d'inclure les a priori sous forme de texte descriptif des objets et leurs relations [START_REF] Lu | Few-Shot Object Detection in Aerial Imagery Guided by Text-Modal Knowledge[END_REF], tandis que d'autres utilisent des cascades de modules d'attention [START_REF] Zeng | Dynamic Cascade Query Selection for Oriented Object Detection[END_REF]. D'autre part, les Processus Ponctuels proposent de résoudre conjointement les vraisemblances relatives à l'image (terme d'attache aux données), et la cohérence de la configuration d'objets elle-même (terme a priori). Premièrement, les approches Processus Ponctuel modélisent les configurations sous forme de géométrie vectorielle, contraignant l'espace d'état par construction. De plus, ces modèles permettent de formuler des a priori explicites sur les configurations en tant que fonctions énergétiques. Cependant, dans la plupart de la littérature [START_REF] Verdié | Detecting parametric objects in large scenes by Monte Carlo sampling[END_REF][START_REF] Schmidt | Forest point processes for the automatic extraction of networks in raster data[END_REF], les termes d'attache aux données reposent sur des mesures de contraste entre les objets d'intérêt et l'arrière-plan. Nous montrons dans cette thèse les limites de ces mesures sur les données satellitaires.

Au lieu d'augmenter la complexité du modèle en y ajoutant par exemple des Transformers, nous proposons dans cette thèse de combiner l'extraction de motifs CNN avec l'approche Processus Ponctuel. Le fondement de cette approche est d'utiliser la sortie d'un CNN comme terme de données pour un modèle de Processus Ponctuel. Dans le Chapitre 4, nous montrons d'abord les limites des mesures basées sur le contraste utilisées comme attache aux données pour les images satellitaires, puis procédons à la construction d'un terme d'attache aux données à partir de la sortie d'un CNN.

Nous proposons plusieurs manières de construire une carte de potentiel :

• Dans un premier temps, nous proposons d'utiliser une méthode de type contraste (correspondance de motifs utilisant une mesure de corrélation) sur la sortie d'un CNN entraîné pour transformer les centres des objets en blobs dans sa sortie (Mabon, Ortner, & Zerubia, 2021).

• Nous proposons ensuite une nouvelle méthode utilisant un CNN entraîné pour inférer un champ de vecteurs pointant vers le centre de l'objet le plus proche. Pour obtenir une carte de probabilité des centres, il suffit de calculer la divergence de ce champ. Cela permet de séparer facilement des objets proches pour un faible coût de calcul (Mabon, Ortner, & Zerubia, 2022b).

• Enfin, nous montrons comment réinterpréter toute sortie CNN de type heatmap (voir partie 2.1.1 ou [START_REF] Law | CornerNet: Detecting Objects as Paired Keypoints[END_REF]) comme une énergie pouvant être utilisée pour notre modèle, à la fois pour la position et les marques de chaque objet (Mabon, Ortner, & Zerubia, 2023a).

Avec notre formulation de l'énergie d'attache aux données sous forme d'interpolation d'une carte de potentiel, nous transformons le calcul de la mesure de contraste en une simple extraction et interpolation de valeurs dans un tenseur 2D ou 3D.

Échantillonnage. Avec ces cartes de potentiel pré-calculées, nous procédons dans la partie 5.1.1 à l'adaptation des méthodes d'échantillonnage des Processus Ponctuels pour utiliser ces informations facilement disponibles. En effet, ces tenseurs définissent des densités sur un espace discret facile à échantillonner et qui se rapproche d'une version tronquée du modèle énergétique (qui est défini dans un espace continu). Tout d'abord, nous appliquons cet échantillonnage guidé par les données au noyau de perturbation de transformation locale (Mabon et al., 2021), puis construisons une carte de naissance pour le noyau de naissance et mort (Mabon et al., 2023a). Ces densités spatiales nous permettent également de concentrer l'échantillonnage parallèle du Processus Ponctuel sur les zones où le nombre de points est susceptible d'être plus élevé et nécessite donc plus de propositions. Enfin, nous exploitons les moteurs de différenciation automatique facilement disponibles pour calculer le gradient d'énergie et effectuer une diffusion pour une exploration plus efficace de l'espace d'état pour un nombre de points fixé.

Estimation des paramètres. Les méthodes précédentes utilisent la Programmation Linéaire pour estimer les poids relatifs des énergies du modèle (Q. [START_REF] Yu | Multiple-Target Tracking by Spatiotemporal Monte Carlo Markov Chain Data Association[END_REF][START_REF] Craciun | Joint Detection and Tracking of Moving Objects Using Spatio-temporal Marked Point Processes[END_REF] à partir d'un ensemble de contraintes. Cependant, cela est sujet à des sur-contraintes et ne prend en compte qu'un nombre limité de points de données. Notre première approche d'estimation des paramètres consiste à proposer une méthode utilisant un SVM, qui recherche la meilleure séparation possible entre les échantillons positifs et négatifs. Grâce à sa frontière perméable, cette méthode permet de considérer beaucoup plus de points de données, même ceux bruités qui contrediraient des contraintes strictes (Mabon et al., 2021). Ensuite, nous proposons d'adapter une méthode utilisée pour les modèles basés sur l'énergie (EBM) (Mabon, Ortner, & Zerubia, 2022a) : la divergence contrastive [START_REF] Hinton | Training Products of Experts by Minimizing Contrastive Divergence[END_REF][START_REF] Du | Implicit generation and modeling with energy based models[END_REF]. Dans un schéma de descente de gradient, nous alternons entre la génération d'échantillons contrastifs/négatifs à partir de l'énergie actuelle et la maximisation de la différence d'énergie entre les échantillons positifs (similaires à la vérité terrain) et les échantillons contrastifs. Cette procédure permet d'estimer non seulement les pondérations énergétiques, mais également les différents paramètres internes des termes énergétiques, qui autrement auraient été définis manuellement (ou avec une procédure distincte pour chacun). Ces deux approches sont détaillées dans la partie 5.2.

Intensité de Papangelou comme score de détection. Pour pouvoir comparer notre modèle avec d'autres modèles de détection d'objets utilisant des CNN, nous introduisons une fonction de notation pour les objets déduits par le modèle de Processus Ponctuel dans la partie 5.3. Ce score est calculé à partir de l'intensité conditionnelle de Papangelou, qui mesure la probabilité d'obtenir un point compte tenu du reste de la configuration. Il s'agit d'une différence essentielle par rapport aux fonctions de score classiques qui ne prennent en compte les points qu'individuellement, quel que soit leur voisinage.

Application. Dans le Chapitre 6, nous appliquons nos méthodes sur des ensembles de données de télédétection ; à la fois sur les benchmarks COWC [START_REF] Mundhenk | A Large Contextual Dataset for Classification, Detection and Counting of Cars with Deep Learning[END_REF] et DOTA [START_REF] Xia | DOTA: A Large-Scale Dataset for Object Detection in Aerial Images[END_REF] qui sont accessibles au public, ainsi que sur certaines données fournies par Airbus Defence and Space (ADS). Même avec la complexité limitée de nos modèles, nous obtenons de bons résultats qualitativement en assurant la régularité de la géométrie déduite. Quantitativement, nos modèles sont meilleurs sur ces données que les autres méthodes CNN auxquelles nous nous comparons. Cette différence de performance est plus prononcé avec des données d'entrée bruitées : nos modèles fournissent de bons résultats qualitatifs et quantitatifs même dans des conditions de bruit défavorables sur l'image d'entrée. Nous démontrons également que la méthode d'estimation des paramètres proposée correspond ou surpasse légèrement celle du réglage manuel des paramètres du modèle, qualitativement et quantitativement, avec l'avantage intrinsèque de limiter les fastidieux essais et erreurs manuels. Enfin, nous montrons que notre modèle peut se généraliser à de nouvelles données en examinant les résultats qualitatifs sur les données ADS après un entrainement du modèle sur DOTA.

Limitations. Cependant, comme notre implémentation n'est que légèrement optimisée (principalement pour permettre la modularité lors de l'expérimentation) et en raison de la nature itérative des méthodes d'échantillonnage par chaines de Markov, le temps d'inférence est significativement plus long que celui des méthodes de détection d'objets purement CNN. Une accélération pourrait être obtenue en approchant la méthode d'échantillonnage : par exemple en utilisant moins d'itérations, avec une inférence approximative comme état initial (par exemple la sortie CNN-LocalMax., voir la partie 6.2.2) et un recuit simulé plus rapide. Cependant, cela s'éloignerait davantage de la configuration d'échantillonnage optimale, mais donnerait des résultats plus rapidement. D'autres solutions existent telles que la discrétisation de l'espace d'état, puis la résolution d'un problème d'optimisation de variables binaires (T. T. [START_REF] Pham | Efficient Point Process Inference for Large-Scale Object Detection[END_REF].

Perspectives. Alors que, dans cette thèse, nous nous sommes concentrés sur la détection de petits véhicules, le cadre que nous proposons peut être facilement adapté à toute autre modalité d'objets où les interactions inter-objet sont clés. Par exemple, les réseaux routiers peuvent être modélisés avec des Processus Ponctuels de segments comme dans [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF], et avoir des a priori importants sur leurs interactions (par exemple un rayon de courbure limité, physiquement et légalement). Les Processus Ponctuels aussi peuvent être appliqués aux données temporelles ; [START_REF] Craciun | Joint Detection and Tracking of Moving Objects Using Spatio-temporal Marked Point Processes[END_REF] Remote sensing, a technology with roots dating back to photography from balloons in the late 19th century (Figure 1.1), has evolved significantly since, to become key for research, commercial use or defense. Over the years, sensors got more sophisticated, in part sparked by the escalation of military intelligence in the context of World War II and subsequently the Cold War. This coupled with an increased need to monitor our environment globally as humanity's effects on it became more concerning. In turns, the growing amount of data produced by those sensors demanded advanced computing solutions for processing, analysis and interpretation. The advancements in remote sensing matches the rise of computer vision: from the early foundation of computer vision in the 1950s, through the Perceptron in 1958 [START_REF] Rosenblatt | The perceptron: A probabilistic model for information storage and organization in the brain[END_REF] and the Hough transform in 1972 [START_REF] Duda | Use of the Hough transformation to detect lines and curves in pictures[END_REF] to the resurgence of neural networks in the late 1980s and the rise of Convolution Neural Networks (CNN), and most recently the adoption of Deep Learning with AlexNet [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] or RCNN [START_REF] Girshick | Rich Feature Hierarchies for Accurate Object Detection and Semantic Segmentation[END_REF] for instance.

At first glance, remote sensing can be understood as "the gathering of information at a distance" [START_REF] Campbell | Introduction to Remote Sensing, Fifth Edition[END_REF]. However, this too broad definition would encompass from microscopy to audio recording. The field of remote sensing, rather focuses on the "observation of the Earth's land and water surfaces by means of reflected or emitted electromagnetic energy" [START_REF] Campbell | Introduction to Remote Sensing, Fifth Edition[END_REF]. In practice, the considered data consists of optical (or near-optical such as infrared) imagery, as well as SAR (Synthetic Aperture Radar) or Lidar (Laser Imaging, Detection, and Ranging) taken from drones, planes, balloons or satellites; in short "from an overhead perspective" [START_REF] Campbell | Introduction to Remote Sensing, Fifth Edition[END_REF].

In this thesis our main goal is to detect small objects, such as vehicles, in satellite images. With resolutions around 0.5m per pixel and the perturbations inherent to the image acquisition process, we try to complement the lack of information in the image with priors on the extracted geometries. We aim at complementing modern CNN approaches with interaction models from Point Processes to regularize the extracted configurations of objects.

This introductory chapter, first presents some practical applications of remote sensing and applications of object detection. We then present the key elements of satellite image acquisition 16 CHAPTER 1 -Introduction with its challenges and constraints. Finally, we broadly review approaches for object detection to finally conclude on the thesis structure.

Computer vision in remote sensing

Thanks to constant technological advances and the rise of open data, the uses for computer vision in remote sensing are growing. As a motivation, we present a few examples of practical applications.

Land use. Satellite images provide a live view of a territory and helps keeping cartographic resources up to date (IGN, 2020). This allows for instance to monitor urban sprawl at a large scale (Une carte d'occupation du sol pour faciliter le suivi de l'artificialisation -Institut -IGN, 2023), which is of great importance as it is a key cause for loss of biodiversity [START_REF] Mckinney | Urbanization, Biodiversity, and Conservation: The impacts of urbanization on native species are poorly studied, but educating a highly urbanized human population about these impacts can greatly improve species conservation in all ecosystems[END_REF], pollution and heath issues [START_REF] Jones | Spatial Distribution of U.S. Household Carbon Footprints Reveals Suburbanization Undermines Greenhouse Gas Benefits of Urban Population Density[END_REF]. Combining object detection in satellite images with up-to-date land registers can have uses for public finances; recently, the French tax office has been experimenting with detection of undeclared swimming pools from satellite imagery ("La détection par intelligence artificielle de piscines non déclarées va être généralisée en France", 2022).

Disaster monitoring. A live view of large swathes of land can prove crucial in monitoring disasters, both for responding to humanitarian crises or in the long run. Satellites can help monitor the ongoing climate crisis, from measuring the rising sea surface temperature (European Sea Surface Temperature, 2023) to producing a live-view of forest fires across the world ("Cartes des incendies dans le monde : suivez en temps réel les feux de forêt observés sur Terre", 2023), (see this monitoring tool from NASA with only 3 hours delay (Earth Science Data Systems, 2015)). Optical satellite imagery is also in use for disaster response: for instance, in September 2023, in response to the earthquake in Morocco, the available satellites of 8 space agencies where programmed to produce pictures of the affected area [START_REF] Bronner | Séisme au Maroc : les satellites peuvent aider les secours à réagir au plus vite. The Conversation[END_REF]. This data allows identifying and measuring the extent of the damage, and helps to coordinate the relief efforts.

Military intelligence. As mentioned previously, remote sensing technologies are deeply rooted in military intelligence. From observation balloons in 1794 during French Revolutionary wars, to World War I with the first plane photographic reconnaissance missions, and the development of spy satellite during the cold-war: monitoring -preferably from an overhead perspective -the presence, amount and nature of military assets is key to any conflict or defense strategy. In recent years, the wide-spread availability of satellite image resources has lead to the rise of Open Source Intelligence (OSINT) (OSINT : aux sources d'un nouveau journalisme ?, 2023), which became a significant source of public information during the recent conflict in Ukraine. Domestic intelligence. However, it is important to note the global increase in use of automated acquisition and detection systems aimed towards domestic surveillance; from the NSA global communication eavesdropping programs unveiled in 2013 [START_REF] Gellman | [END_REF], to the ongoing growth of the Chinese surveillance state (facial recognition, phone trackers, voice recognition etc.) [START_REF] Xiao | Video: China's Surveillance State Is Growing. These Documents Reveal How[END_REF], or France's recent legislation on the usage of police drone during protests, or the experimentation on automatic facial recognition opening the way for the trivialization of such practices [START_REF] Verdon | Proposition de loi sur la reconnaissance faciale : Un pas de plus vers la surveillance généralisée ? The Conversation[END_REF]. Although non-military satellite technology is not capable enough yet for these practices, the use of drone systems coupled with automated recognition algorithms can raise ethical concerns.

Satellite image acquisition

To better understand the specificities and constraints of detection in satellite images, we present in this section a few key elements of the image acquisition process. Most of the information in this section have been sourced from [START_REF] Tupin | Remote sensing imagery[END_REF], to which we refer the reader for a more thorough explanation.

Orbits and cycles. Satellites for civilian imagery, usually orbit from 450 to 910 km around the Earth. The orbital plane is often inclined with respect to the equatorial plane in order to cover as much ground as possible. For instance a polar orbit passes above the North and South poles, describing an orbit orthogonal to the equatorial plane. Lower orbits require less fuel to launch and make for higher resolution images, however as satellites get lower in orbit the effects of atmospheric drag get higher. Even if satellites are loaded with some fuel for trajectory corrections, low orbits demand more corrections because of the drag. Thus, the lifespan of low orbit satellites is often limited. Due to the combination of the satellite's orbit and the Earth's rotation (eastward with North up), a position on the orbit does correspond to a single position on the ground; at every orbit the position of the satellite gets shifted to the west (see Figure 1.2). The time a satellite takes to come back to the same orbit position and location relative to the ground is called the cycle of the orbit (e.g. 1 day for Formosat-2, 26 days for SPOT or Pleiades NEO). This cycle time will determine the minimal time difference between two acquisitions of the same area; which plays a key role in CHAPTER 1 -Introduction change detection. This can be mitigated by having constellations of satellites, better covering the globe at all times.

Passive sensors: Optical. There exists two types of sensors: passive and active. The passive sensors measure the back scattering of sunlight (or artificial lights) on the surface. For thermal sensors it can measure the ground radiation. The first sensors where film stock. In the 1960s US satellite platforms such as GAMBIT-1 or CORONA would drop a reentry canister back into the atmosphere to be recovered midair by plane or on the ground (60th Anniversary of the First GAMBIT-1 Photoreconnaissance Satellite Flight, 2023). Modern satellites use a linear array of Charged Coupled Device (CDD) that scan an area as the satellites moves over it (as a document scanner does), this is referred to as pushbroom (see Figure 1.3). 2D CCD arrays (that we find in everyday digital cameras) are marginal in use up to now, although the upcoming CO3D constellation will embark two 2D CCD sensors, one for RGB bands, the other for NIR (CO3D Constellation, 2023). The interval between two rows of pixels ∆y is given by the time interval ∆t multiplied the satellite ground speed V S . Adapted from [START_REF] Tupin | Remote sensing imagery[END_REF].

The footprint of a satellite and sensor describes the intersection of the perspective cone and ground surface, in which two objects cannot be resolved. The footprint is dependent on the optical system and sensor, as well as incidence (angle between ground normal and sensor), altitude and ground inclination. For a horizontal ground, the lowest footprint is at low altitude and incidence; i.e. the optical satellite best acquires images at Nadir (satellite pointing towards the center of the Earth). Off-Nadir images have less spatial resolution, and are more prone to occlusion due to mountains, trees or tall buildings. Note that for a satellite at Nadir, as we move further from the optical axis (vector -→ SO in Figure 1.3), the light paths get progressively off Nadir.

Active sensors: Synthetic Aperture Radar (SAR). Synthetic Aperture Radar (SAR) satellites are active sensors; they emit electromagnetic wave with an antenna, that gets reflected on the surface and measured back by the antenna. The system measures the signal return time, thus the SAR footprint depends on the temporal resolution and the incidence. Here lower incidence increases footprint; i.e. SAR satellites are preferably pointing at the ground with an angle. The movement of the satellite relative to the ground introduces some more subtleties to this system; we will not go into further details about of SAR acquisition in this thesis, we rather refer the reader to [START_REF] Tupin | Remote sensing imagery[END_REF] or [START_REF] Campbell | Introduction to Remote Sensing, Fifth Edition[END_REF] for more detailed information.

Spectral resolution. Optical sensor measures the quantity of light within a given range of wavelengths. For instance panchromatic images are produced by capturing the visible spectrum from 400 to 700 nm, producing a grayscale image corresponding to the integrated intensity over the bandwidth, exposure time and sensor area. The bandwidth can be split into bands, to produce multispectral images. The most common we encounter every day is splitting into Red, Green and Blue bands (RGB) to be able to reproduce color images for the human eye. Satellites can have from ten to hundreds of bands. The fine spectral resolution allows discerning different materials and surfaces given their reflectance properties. Multispectral sensors often include bands beyond the visible spectrum such as near infrared (NIR), which is strongly reflected by vegetation. For instance Pleiades NEO satellites (PNEO-3 and PNEO-4) have 7 spectral bands: panchromatic (450 -800 nm), RGB (450 -520, 530 -590, and 620 -690 nm), NIR (770 -880 nm), deep blue (400 -450 nm, used for instance for oceanography) and red edge (700 -750 nm, used for vegetation) (Pléiades Neo, n.d.). As bands get finer, the flux of light gets lower, thus the signal-to-noise ratio worsens. As a compromise narrowband systems often increase sensor surface, at the cost of increased footprint (i.e. lower spatial resolution). For instance the panchromatic band in Pleiades NEO has a resolution of 30 cm while the other narrower bands have a resolution of 1.2 m.

Atmospheric Perturbations. "It is necessary to talk about the atmosphere, since it is situated between the satellite and the surface" [START_REF] Tupin | Remote sensing imagery[END_REF]. The most glaring hindrance to optical imaging is clouds; those hide the surface and project shadows. Their non-solid edges make them difficult to isolate and correct for. However, clear skies also interact with the light the satellite is capturing. The atmosphere both absorbs and scatters the light depending on its chemical composition and density. Both do not have uniform effect across the light spectrum. We shall also mention, surface effects, such as specular reflections on mirror like surfaces, which can cause lens flare artifacts.

Image transmission. Nowadays, satellites beam back the acquisitions by radio; both the bit rate and transmission window are limited as the satellite is not always in range of the ground antenna. The size of the transmitted signal depends on the number of bands, pixels and radiometric resolution (the number of quantized values for the intensity of one pixel in one band). Usually, the ground to satellite link has a low transmission rate as it only transmits commands to the satellite. Meanwhile, satellite to ground links have higher throughput to accommodate for the large amounts of data to send back. Some Earth observation satellites (such as Sentinel) can beam back data trough higher orbit relay satellites, such as the European Data Relay Satellites, EDRS-A and C (European Data Relay Satellite System (EDRS) Overview, n.d.). The transmission rate is still limited due to the distances (loss of signal energy) involved with geostationary satellites.

Image processing. Once on ground, the image can be processed. Modeling the atmospheric interactions allows correcting for part of the above-mentioned absorptions. The images are bespangled to conform to an orthoimage; i.e. transformed image that corresponds to a cartographic geometry. To compensate for the limited resolution of fine spectral bands, one can combine the high spatial resolution of panchromatic images with the fine spectral resolution of multispectral. This process is known as pansharpening.

From the multispectral output, spectral indices are built. Those combine the information of multiple bands into a single index that corresponds to a specific element of interest. For instance NDVI combines red (R R ) and near infrared (R N IR ) bands. It is a great tool to quickly identify vegetated areas without any complicated image processing:

NDVI = (R R -R N IR ) (R R + R N IR ) (1.1)

Object detection

As shown previously, satellite images contain many artifacts and occlusions compared to usual images; i.e. photographs we encounter every day, taken at human eye-level, usually with a close centered subject, where objects of interest represent a significant proportion of the image. Given the low spatial resolution of satellite images, the objects of interest may be only a few pixels large.

As the visual information is limited, human operators often resort to their prior knowledge of the object of interest when analyzing an image. Similarly, we aim to build models that factor in this prior.

The goal of this thesis is to study the combination of modern Convolutional Neural Network (CNN) approaches with geometrical models using priors or constraints. While CNN models have shown to be really efficient at extracting texture information, they fail to encompass longer range interactions between objects without drastically increasing complexity and parameters count. On the other hand, Point Process (PP) models propose to model configurations of geometrical objects as a stochastic process. First, the Point Process construction itself defines the type of geometries that are extracted. Also, the stochastic model contains priors that can guide the model towards configurations that are more likely a priori.

Image probability model

The Point Process approach derives from stochastic models of image analysis. Considering the image as a set of pixels in S d ⊂ Z 2 , we can model it as a random variable X in probability space (Ω, A, P), where grayscale values correspond to the realizations of function X:

X : Ω → R S d . (1.2)
If we consider the image as a noisy observation of underlying phenomenon Φ, a Bayesian model for the image is then:

P(X = X) = ∑ y P(X = X|Φ = y)P(Φ = y), (1.3)
where image X and configuration y are respective realizations of X and Φ. Here P(X = X) corresponds to the marginal law of observations, built from the law of observations P(X = X|Φ = (1.4)

The Bayes model provides a class of natural estimators: Bayes estimators. An estimator of Φ is a function Φ * (X) of observation X, which we want to best approach the underlying values of Φ that generated image X. This relies on Bayes cost:

R B (Φ * ) = E Φ [ E X|Φ [L(Φ * , Φ)] ]
.

(1.5)

For cost L(Φ * , Φ) = 1(Φ ̸ = Φ * )
, the optimal estimator is the Maximum A Posteriori (MAP)

estimator:

y * = arg max y P(Φ = y|X = X). (1.6)
If we consider all pixels to be independent of each other we get:

y * = arg max y ∏ ρ∈S d P(Φ = y|X ρ = X ρ ).
(1.7)

Point Process

Point Processes allow to model such phenomenons Φ, as set of parametric objects representing geometrical objects such as segments, circles, rectangles, polygons etc... We will go further into details of Point Processes later on. Here we consider configurations y as a set {y 1 , . . . , y n(y) } of n(y) elements y that represent our objects of interest. Within the Bayesian framework, an observation model can be built as: (1.8) In this formula, we apply a foreground (ρ ∈ y) or background (ρ / ∈ y) model whether a pixel ρ belongs to the silhouette of y or not. A common model if the silhouette corresponds to bright pixels on a darker background is to model the foreground P(X ρ = X ρ |ρ ∈ y) as a high value Gaussian law. The background P(X ρ = X ρ |ρ / ∈ y) is then mapped to a lower value [START_REF] Baddeley | Stochastic geometry models in high-level vision[END_REF][START_REF] Perrin | Tree crown extraction using marked point processes[END_REF].

P(X ρ = X ρ |Y = y) = 1(ρ ∈ y)P(X ρ = X ρ |ρ ∈ y) + 1(ρ / ∈ y)P(X ρ = X ρ |ρ / ∈ y).
However, this Bayesian approach has some limitations [START_REF] Ben Hadj | Estimation des paramètres de modèles de processus ponctuels marqués pour l'extraction d'objets en imagerie spatiale et aérienne haute résolution[END_REF]:

• Non-homogeneous backgrounds are harder or impossible to model. The simple Gaussian model for P(X ρ = X ρ |ρ / ∈ y) is no longer usable.

• The likelihood of pixels of the foreground class does not rely on the morphological properties of the class to extract. For instance the model may try to fit as many small objects it can in a large silhouette [START_REF] Craciun | Joint Detection and Tracking of Moving Objects Using Spatio-temporal Marked Point Processes[END_REF] as illustrated in Figure 1.4.

For those reasons, [START_REF] Ben Hadj | Estimation des paramètres de modèles de processus ponctuels marqués pour l'extraction d'objets en imagerie spatiale et aérienne haute résolution[END_REF] propose a detector model, formulating the external energy U d as:

U d (y, X) = w d ∑ y∈y V d (y)
, with (1.9) where w d is a scalar weight, and U p is the prior model formulated as an energy. Meanwhile, V d (y) measures the statistical difference between the inside and outside of object y, with a low value for high contrast (likely an object) and high values for low contrast (likely not an object). Still these contrast measures rely on the high and consistent contrast of the object against its background. This approach reaches its limits when the background start to exhibit contrasted elements that should not be detected for instance. We show a practical example of these limits in Section 4.2.1.1.

P(Φ = y|X = X) ∝ exp (-(U d (y, X) + U p (y))) , (1.10) (a) (b) (c)
Notably, (T. [START_REF] Li | Feature Extraction and Tracking of CNN Segmentations for Improved Road Detection from Satellite Imagery[END_REF] propose to perform such contrast based approach on the output of a CNN trained for segmentation; that way the authors fall back into a simple foreground and background where the contrast measure is efficient.

Convolutional Neural Networks

While the above image model with independent pixel was motivated in part by the computational limitations, the increasing availability of computational power eased the rise of Convolutional Neural Networks (CNN). Although introduced for image classification in 1989 [START_REF] Lecun | Backpropagation Applied to Handwritten Zip Code Recognition[END_REF], it was only in the 2010s that CNN gained widespread adoption with AlexNet [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] for classification or RCNN [START_REF] Girshick | Rich Feature Hierarchies for Accurate Object Detection and Semantic Segmentation[END_REF] for object detection.

Convolutional Neural Networks process data from a grid-like topology through a sequence of convolution and pooling (spatial aggregation) operations. We will go further in depth about CNN later in the thesis. Deep convolutional models learn the convolution filters for a specific task from large image databases such as ImageNet [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] for instance. Initially those were built towards classification of images (i.e. assigning a category to the whole image, usually describing its main content). Then [START_REF] Girshick | Rich Feature Hierarchies for Accurate Object Detection and Semantic Segmentation[END_REF] propose RCNN that classifies a series of bounding boxes to build the first two stages object detector. Segmentation model such as Unet [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] perform pixel-wise classification. Those produce real-valued maps in the range [0, 1], that are interpreted as class probabilities. For instance, the probability of pixel ρ belonging to class C 1 (amongst N C classes) is given as: .11) with F (X) is the CNN output tensor in R S d × 1,N C , and N C the number of classes. This approach is leveraged in heatmap based object detection such as (X. [START_REF] Zhou | Objects as Points[END_REF].

 P (ρ ∈ C 1 |X) = exp (F (X)[ρ, C 1 ]) ∑ N C k=1 exp (F (X)[ρ, C k ]) , ( 1

-Proposed approach 23

The heatmap corresponds to a centerness probability map, measuring the probability of each pixel to be an object center. One has to be careful about the interpretation of these as probabilities: as (C. [START_REF] Guo | On Calibration of Modern Neural Networks[END_REF] show, these measures are often poorly calibrated and should rather be considered as scores.

Nonetheless, this score measures the presence or not of an object, with a lightweight convolutional model which kernels are learned from data.

Proposed approach

The starting point of our work -inspired from (T. [START_REF] Li | Feature Extraction and Tracking of CNN Segmentations for Improved Road Detection from Satellite Imagery[END_REF] -is to leverage the score  P (ρ ∈ C 1 |X) and use it to build the external energy of the Point Process V d , replacing the contrast measure.

In this PhD thesis, we propose the incorporation of interaction models into object detection methods, while taking advantage of the capabilities of deep convolutional neural networks. Satellite images are inherently imperfect due to atmospheric disturbances, partial occlusions and limited spatial resolution. To compensate for this lack of visual information, it becomes essential to incorporate prior knowledge about the layout of objects of interest.

On the one hand, methods based on CNN are excellent for extracting patterns in images, but struggle to learn object-to-object interaction models without having to introduce attention mechanisms such as Transformers that considerably increase the complexity of the model. For example, some approaches propose to include a priori information in the form of descriptive text about objects and their relationships [START_REF] Lu | Few-Shot Object Detection in Aerial Imagery Guided by Text-Modal Knowledge[END_REF], while others use cascades of attention modules [START_REF] Zeng | Dynamic Cascade Query Selection for Oriented Object Detection[END_REF].

On the other hand, Point Processes propose to jointly solve the likelihood relative to the image (data term), and consistency of the object configuration itself (prior term). Firstly, Point Process approaches model configurations as vector geometry, constraining the state space by construction. In addition, these models allow explicit priors relative to configurations to be specified as energy functions. However, in most of the literature [START_REF] Verdié | Detecting parametric objects in large scenes by Monte Carlo sampling[END_REF][START_REF] Schmidt | Forest point processes for the automatic extraction of networks in raster data[END_REF], data terms rely on contrast measures between objects of interest and background. We illustrate the limitations of these measures on satellite data later in this thesis.

Instead of increasing the complexity of the model by adding, for example, Transformers, we propose in this thesis to combine CNN pattern extraction with the Point Process approach. The starting point of this approach is to use the output of a CNN as the data term for a Point Process model. From the latter we derive more efficient sampling methods for the Point Process that do not rely on application specific heuristics. Finally, we propose to bridge the gap in terms of parameters estimation using modern learning techniques inspired from Energy Based Models.

Thesis structure

The thesis is organized as follows. We choose to present our contributions structured into the three main elements of our models (energy model in Chapter 4, sampling and parameters estimation in Chapter 5) rather than going through the incremental improvements of each model we published (which would induce many repetitions) :

• Chapter 2 provides a review of the State Of The Art (SOTA) approaches for object detection both with CNN and PP based methods, with a focus on remote sensing applications.

• Chapter 3 goes through the fundamental theory on Point Process and the basics of Convolutional Neural Networks necessary to build our models.

• Chapter 4 introduces several ways to build our Point Process model, in which we contribute with several novel approaches:

-Building data terms for PP from CNN.

-Interpreting classical CNN output as energy to incorporate in the PP model.

• Chapter 5 develops the methods to sample and estimate the parameters of the model. Our main contributions are:

-Thanks to CNN output potential maps we adapt sampling methods for more efficiency (parallel sampling with cell picking, birth with density from truncated model).

-We propose two parameter estimation methods: the first based on Support Vector Machines ; the second based on gradient descent and contrastive divergence inspired from the literature related to Energy Based Models.

-We propose using the Papangelou intensity as a per object score value in order to compute precision recall curves and compare our method to the classical ones.

• Chapter 6 shows application of our model and methods on synthetic and real data on benchmarks and data provided by Airbus Defense and Space (ADS).

• Chapter 7 summarizes contributions, limitations and perspectives for future works.

CHAPTER 2

State of the art In this chapter we review the literature around object detection and its related methods. As our aim is to build a model combining Point Processes (PP) and Convolutional Neural Networks (CNN), we first focus on the methods using CNN models for object detection, and more specifically applied to remote sensing data. We then look into the approaches based on PP, applied to object detection. Finally, we delve into the realm of Energy Based Models as their approach will be of use when training our model later on.

-27

Looking back at the references used for this thesis, we identify three main topics, which we will explore in this chapter. Figure 2.1 shows relations between publications and topics, as well as these thematic clusters. As our aim is to build a model combining Point Processes (PP) and Convolutional Neural Networks (CNN), we first focus on the methods using a CNN for object detection, and more specifically applied to remote sensing data. We then look into the approaches based on Point Processes, applied to object detection. Finally, we delve into the realm of Energy Based Models (EBM) as their approach will be of use when training our model later on.

Convolutional Neural Networks for object detection

In this section we review the object detection first through the emergence of Convolutional Neural Networks (CNN) and then their use for remote sensing. We distinguish usual images -photographs usually taken at human eye-level with large mostly-centered objects of interest -from the images seen in remote sensing (or in microscopy for instance), as those have different characteristics that will be discussed further in this chapter.

The rise of CNN models for object detection

This first part presenting a short history of object detection is mainly sourced from [START_REF] Zou | Object Detection in 20 Years: A Survey[END_REF], in which the reader will find an even more thorough survey of object detection over the last 20 years. First it is important to mention that object detection is adjacent to many other computer vision tasks such as segmentation, classification, change detection etc... Initially, classical detection methods were based on handcrafted features such as for [START_REF] Viola | Rapid object detection using a boosted cascade of simple features[END_REF]. This real-time face detection algorithm uses a sliding window with some handcrafted rectangular features (see Figure 2.2) for which a clever use of integral images allows to simplify the algorithm to a few lookup operations.

As handcrafted features showed their limit, object detection went through a shift with the arrival of Convolutional Neural Networks. Thanks to previous key works on backpropagation [START_REF] Lecun | Backpropagation Applied to Handwritten Zip Code Recognition[END_REF], CNNs took off in computer vision in 2012 with [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] introducing a large scale training data (ImageNet) and proposing a CNN model (AlexNet) trained to perform image classification in usual images.

Two stage detectors. In 2014 Girshick et al. propose Region with CNN features (RCNN) in [START_REF] Girshick | Rich Feature Hierarchies for Accurate Object Detection and Semantic Segmentation[END_REF] to bridge the gap from classification to object detection. RCNN works by proposing a set of objects as boxes. Each box is then translated into a vector of features by a CNN model trained on ImageNet [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF], a Support Vector Machine (SVM) then predicts the presence and class of the potential object in the box (see Figure 2.3). The large number 2.1 -2.1.1 The rise of CNN models for object detection 29 of proposals that overlap is a limitation of that method (over 2000). In later works [START_REF] Girshick | Fast R-CNN[END_REF] propose FastRCNN, simultaneously training the detector and a bounding box regressor of the model, making for a great improvement in speed and performance. Even though [START_REF] Ren | Faster R-CNN: Towards real-time object detection with region proposal networks[END_REF] introduce the region proposal network, computational redundancy remains due to the high number of proposals. Such approach of box proposal then classification is referred to as two stage approach. [START_REF] Girshick | Rich Feature Hierarchies for Accurate Object Detection and Semantic Segmentation[END_REF].

From these works, improvements are proposed to alleviate the limitations of the initial two stage models: for instance, (K. [START_REF] He | Spatial Pyramid Pooling in Deep Convolutional Networks for Visual Recognition[END_REF] propose Spatial Pyramid Pooling (SPP) to bypass the fixed input size restriction of previous CNN models with a new pooling method. Later on Feature Pyramid network (FPN) proposed in (Lin, Dollar, et al., 2017), uses features of deeper layers -previous approaches would use only the last layer -introducing top down and lateral connections. This architecture integrated in FastRCNN gives improved results.

Single stage detectors. The first single stage object detector is introduced in [START_REF] Redmon | You Only Look Once: Unified, Real-Time Object Detection[END_REF] as YOLO (You Only Look Once). It shows greater performances in terms of object detection metrics and speed than the two-stage approaches at the time. While YOLO initially struggles with small objects, the following version, YOLOv4 [START_REF] Bochkovskiy | YOLOv4: Optimal Speed and Accuracy of Object Detection[END_REF], YOLOv7 [START_REF] Wang | YOLOv7: Trainable Bagof-Freebies Sets New State-of-the-Art for Real-Time Object Detectors[END_REF] and so on, incrementally alleviate these issues.

In their work [START_REF] Lin | Focal loss for dense object detection[END_REF] propose that the reason one stage detectors kept being outperformed by modern two stage detectors was the imbalance in foreground and background classes. They propose a new focal loss to put more focus on hard, misclassified elements when training the single stage CNN model.

For most previous methods, anchor boxes are used as reference box proposals for classification and regression. As the objects are varied in size and shape, the user has to set up these reference boxes tailored to the application. CornerNet [START_REF] Law | CornerNet: Detecting Objects as Paired Keypoints[END_REF] propose to instead view the bounding box estimation as a keypoint detection problem: for each object they infer the two points that would define a bounding box, the points are matched together by some extra embedding information. The estimation of the points' location is done though some heatmap: raster maps corresponding to a pseudo probability map of points presence. Those are used not only for object localization: for instance (Y. [START_REF] Guo | Density Map-based vehicle counting in remote sensing images with limited resolution[END_REF] use heatmaps to count vehicles in satellite images. CenterNet simplifies this approach further (X. [START_REF] Zhou | Objects as Points[END_REF], considering objects to be a single point (their center) extracted through a heatmap. From the center location is regressed the other attributes such as size, orientation, location, and pose (see Figure 2.4). Here the heatmap resolution is coarser than the image resolution, the local offset allows compensating for the discretization error. From (X. [START_REF] Zhou | Objects as Points[END_REF].

Attention and transformers. Attention mechanisms were introduced in Natural Language Processing (NLP) [START_REF] Bahdanau | Neural Machine Translation by Jointly Learning to Align and Translate[END_REF] to help learn relations between distant words in a sentence, and quickly gained popularity in their field. Works in [START_REF] Vaswani | Attention is all you need[END_REF] introduced transformers, a deep neural network block relying on attention mechanisms, that aims to use global dependencies between input elements, while avoiding recurrence based models1 . It did not take long for the object detection community to adapt this method to object detection. In 2020 [START_REF] Carion | Endto-End Object Detection with Transformers[END_REF] propose DEtection TRansformer (DETR) -which as the name suggests uses transformers for object detection -viewing detection as a set prediction problem and using transformers to learn relations between the elements of these sets; i.e. relations between objects. We discuss of the limitations of attention mechanisms in next Section 2.1.2.

CNN for detection in remote sensing

We now focus on the application of the object detection methods above to remote sensing data. We articulate this section through the multiple challenges that remote sensing data exhibits and the proposed methods to alleviate those.

Small objects. In remote sensing, and more specifically with satellite imagery, objects of interest can be of very small size (only a few pixels). A first approach is to use super-resolution methods to compensate for the object size limitations of the detection model as in (J. [START_REF] Zhang | SuperYOLO: Super Resolution Assisted Object Detection in Multimodal Remote Sensing Imagery[END_REF]. We see it used also for microscopy images in [START_REF] Mayo | Detection of Parasitic Eggs from Microscopy Images and the emergence of a new dataset[END_REF] using Pix2Pix [START_REF] Isola | Image-to-Image Translation with Conditional Adversarial Networks[END_REF] for rescaling images. One limitation of super-resolution approaches is the very large memory footprint of resulting images (images being large in the first place, even at low resolution). Moreover, deep-learning super-resolution models have a risk of hallucinating patterns; i.e. introducing nonexisting and non-relevant information into the scaled-up image. The Feature Pyramid Network (Lin, Dollar, et al., 2017) introduced previously, is great at improving detection at multiple scale, as the lateral connections in the architecture allow to maintain a good localization while detecting 2.1 -2.1.2 CNN for detection in remote sensing 31 object at higher scales in the deeper part of the model. This type of architecture is similar to the Unet introduced in [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF].

It is to note that small objects induce some issues on the metrics used to evaluate or train the models. Indeed, the Intersection Over Union (IOU) becomes instable for small objects; [START_REF] Jeune | Rethinking Intersection Over Union for Small Object Detection in Few-Shot Regime[END_REF] propose a more perceptually accurate measure for small objects.

Large images. Satellite images span huge areas of land, resulting in oftentimes very large images. It requires CNN approaches to be able to handle varied image size, and be able to process patched images; as memory is limited, one needs to be able to split the image, process each patch independently, and stitch back the results without inconsistencies in the output. Fully Convolutional Networks (FCN) formalized in [START_REF] Long | Fully Convolutional Networks for Semantic Segmentation[END_REF] for segmentation, then in [START_REF] Tian | FCOS: Fully Convolutional One-Stage Object Detection[END_REF] for object detection (see [START_REF] Sun | An Anchor-Free Detection Method for Ship Targets in High-Resolution SAR Images[END_REF] for an application in remote sensing), use CNN with only convolution and pooling operations, allowing for translation invariance. FCN are introduced in biomedical image segmentation with the Unet [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF].

Limited spatial resolution. In remote sensing, the sensor is intrinsically quite distant from the observed objects, implying a limited spatial resolution (usually measured in meters per pixel). This limited resolution in turn induces a limited amount of visual information to perform object detection from. Combines with sensor noise or atmospheric perturbation, remote sensing detection methods need to innovate in compensating the limited signal. Approaches such as (LaLonde, [START_REF] Lalonde | ClusterNet: Detecting Small Objects in Large Scenes by Exploiting Spatio-Temporal Information[END_REF] [START_REF] Corsel | Exploiting Temporal Context for Tiny Object Detection[END_REF] use the temporal information from image time series to improve detection of small objects. However, static small objects remain difficult to detect with that method. Others propose using the multiple modalities of data the satellite sensor can produce; using the multiple spectral bands instead of optical grayscale or RGB. For instance [START_REF] Belmouhcine | Multimodal Object Detection in Remote Sensing[END_REF] propose to fuse sensors with attention mechanisms, and (J. [START_REF] Zhang | SuperYOLO: Super Resolution Assisted Object Detection in Multimodal Remote Sensing Imagery[END_REF] combined this with super resolution. The method proposed in [START_REF] Lu | Few-Shot Object Detection in Aerial Imagery Guided by Text-Modal Knowledge[END_REF]) uses text-modal descriptors to introduce prior knowledge on the objects into the model and their relation. For instance "An airport consists of [...] runways for planes to take off and land" introduces a relation between the airport, runway and plane objects and their co-occurence. Similarly, the model proposed in (K. [START_REF] Zheng | Auto Learner of Objects Co-Occurrence Knowledge for Object Detection in Remote Sensing Images[END_REF] learns the co-occurence of objects in images in the training data. Another approach is to model the interaction between objects. Markov Random Field approaches [START_REF] Moser | Land-Cover Mapping by Markov Modeling of Spatial-Contextual Information in Very-High-Resolution Remote Sensing Images[END_REF]) model pixel spatial context, which is combined with CNN for image segmentation in [START_REF] Pastorino | Semantic Segmentation of Remote-Sensing Images Through Fully Convolutional Neural Networks and Hierarchical Probabilistic Graphical Models[END_REF]. Works in (Z. [START_REF] Zheng | FarSeg++: Foreground-Aware Relation Network for Geospatial Object Segmentation in High Spatial Resolution Remote Sensing Imagery[END_REF] propose to relate foreground and background elements while [START_REF] Cao | Vehicle Detection Algorithm Based on Background Features Assistance in Remote Sensing Image[END_REF] use transformer to learn object-level relations. Finally, [START_REF] Zeng | Dynamic Cascade Query Selection for Oriented Object Detection[END_REF] use prior on objects alignment (thus considering object level spatial interaction) within a cascade of attention modules.

Attention models applied to images remain limited in use due to their limitations in terms of memory footprint: often time the solution would be to sub-sample (loosing spatial precision) or consider attention locally (loosing longer range interactions) [START_REF] Cherel | Attention stochastique basée patchs pour l'édition d'images[END_REF]. Alternatively [START_REF] Cherel | Attention stochastique basée patchs pour l'édition d'images[END_REF] propose an efficient attention layer based on a stochastic algorithm for patch matching which is used for determining approximate nearest neighbors.

CHAPTER 2 -State of the art

The limited signal-to-noise ratio raises concerns about the robustness of the developed methods. Works in [START_REF] Mei | A Comprehensive Study on the Robustness of Image Classification and Object Detection in Remote Sensing: Surveying and Benchmarking[END_REF] and (H. [START_REF] He | On the Robustness of Object Detection Models in Aerial Images[END_REF] study the robustness of remote sensing computer vision methods. [START_REF] Mei | A Comprehensive Study on the Robustness of Image Classification and Object Detection in Remote Sensing: Surveying and Benchmarking[END_REF] show the effect of natural perturbations (Gaussian blur, fog, etc...) as well as adversarial attacks including background obfuscation patterns on which to place objects to make those virtually invisible to the detector.

Oriented objects. The detection of oriented objects raises some challenges compared to the detection of horizontal (non-oriented) bounding boxes. Importantly the angle parameter defining the oriented rectangle introduces a discontinuity (or cyclic aspect) to the parameter space which poses problem when computing the loss over the inferred configurations of rotated objects. To bypass this limitation [START_REF] Xu | Rotated Object Detection with Circular Gaussian Distribution[END_REF] propose reformulating angles as a circular Gaussian density, while [START_REF] Yao | On Improving Bounding Box Representations for Oriented Object Detection[END_REF] or (D. [START_REF] Yu | An Anchor-Free and Angle-Free Detector for Oriented Object Detection Using Bounding Box Projection[END_REF] propose redefining the oriented bounding box representation into a continuous one. Finally, [START_REF] Llerena | Gaussian Bounding Boxes and Probabilistic Intersection-over-Union for Object Detection[END_REF]) and (Z. [START_REF] Li | Gaussian synthesis for highprecision location in oriented object detection[END_REF] propose a fuzzy formulation of bounding boxes using Gaussian distributions and a probabilistic IOU. The above-mentioned IOU issues caused by small objects is even more present with oriented objects, where a small angle error can drastically change the IOU as [START_REF] Xu | Rotated Object Detection with Circular Gaussian Distribution[END_REF] note.

Annotated data is costly. As for any deep learning approach, annotated data is key to training models. As data annotation is a time-consuming and arduous task, some methods try to learn efficient model with limited amounts of data. Weakly supervised methods propose to learn models on sparsely annotated data to extract more fine-grained information at inference. For instance (M.-T. [START_REF] Pham | Weakly supervised marine animal detection from remote sensing images using vector-quantized variational autoencoder[END_REF]) train a Variational Auto-Encoder (VAE) to perform image reconstruction, and use it to perform detection of marine animals. Similarly, (J. [START_REF] Bai | Localizing From Classification: Self-Directed Weakly Supervised Object Localization for Remote Sensing Images[END_REF] train a model to classify images (requiring simpler annotations than object location Ground Truth (GT)). The authors then use Class Activation Maps (CAM, introduced in (B. [START_REF] Zhou | Learning Deep Features for Discriminative Localization[END_REF]) to extract the object(s) location in the image looking at which part of the image activates the classification output. While weakly supervised models are trained on sparser data, few shot models are built to learn with limited data samples, for instance being able to categorize never seen classes (as humans can do) [START_REF] Antonelli | Few-Shot Object Detection: A Survey[END_REF]. Previously mentioned works in [START_REF] Lu | Few-Shot Object Detection in Aerial Imagery Guided by Text-Modal Knowledge[END_REF] combine few shot model with the text-modal priors.

Lightweight models for onboard processing. A challenge for some applications in remote sensing is the ability to apply the algorithm in flight, aboard the satellite for instance. Onboard computing faces multiples challenges: the monetary and energetic cost of launching mass (thus processing power) in orbit, the necessity of systems to be stable and reliable (no on-site maintenance possible), limitations with power and cooling, etc... Thus, lightweight models for remote sensing computer vision are key to some applications. For that matter [START_REF] Hu | Efficient-Lightweight YOLO: Improving Small Object Detection in YOLO for Aerial Images[END_REF] propose to train a model based on YOLO with limited number of parameters. Other methods such as neural network reduction [START_REF] Vandame | Réduction de modèles neuronaux profonds par adaptation de l'architecture[END_REF] utilize distillation [START_REF] Hinton | Distilling the Knowledge in a Neural Network[END_REF] to transfer knowledge from a larger model to a sparser one. Lastly, quantization (limiting the number of bits the model parameters are stored on) is proposed in (J. Zhang, Lei, Xie, Li, et al., 2023) using distillation. From [START_REF] Schmitt | There is No Data Like More Data -Current Status of Machine Learning Datasets in Remote Sensing[END_REF].

Remote sensing datasets

As data is the key to deep learning approaches for remote sensing, [START_REF] Schmitt | There is No Data Like More Data -Current Status of Machine Learning Datasets in Remote Sensing[END_REF] provide a review of available datasets in remote sensing (see Figure 2.5). Amongst the object detection datasets we identify a few where our approach could bring improvements on:

• COWC [START_REF] Mundhenk | A Large Contextual Dataset for Classification, Detection and Counting of Cars with Deep Learning[END_REF]: satellite optical images with vehicle position annotations.

As this dataset only provides positions and no geometries, we only use this dataset for the first non-marked point process models.

• DOTA v2.0 [START_REF] Xia | DOTA: A Large-Scale Dataset for Object Detection in Aerial Images[END_REF]: aerial and satellite optical images with oriented rectangles annotations for multiples classes including ground vehicles, planes, boats, peers etc... It provides a large diversity of contexts and objects density, we use it extensively for training and estimation of our methods.

• spacenet 4 [START_REF] Weir | SpaceNet MVOI: A Multi-View Overhead Imagery Dataset[END_REF]: optical satellite images with annotated building footprints.

While not investigated in this thesis, applying point process models to building extraction is quite relevant as the interaction priors are strong [START_REF] Ortner | A Marked Point Process of Rectangles and Segments for Automatic Analysis of Digital Elevation Models[END_REF] • BIRDSAI [START_REF] Bondi | BIRDSAI: A Dataset for Detection and Tracking in Aerial Thermal Infrared Videos[END_REF]: aerial thermal images from drones, with person tracking annotations. Point processes are also useful to take into account priors on dynamics [START_REF] Craciun | Joint Detection and Tracking of Moving Objects Using Spatio-temporal Marked Point Processes[END_REF]. We thus mention this dataset for future works.

• AU-AIR [START_REF] Bozcan | AU-AIR: A Multi-modal Unmanned Aerial Vehicle Dataset for Low Altitude Traffic Surveillance[END_REF]: drone optical images with horizontal bounding boxes for multiple classes of objects. While this dataset would be interesting to model priors over dynamics, the viewing angle and horizontal bounding boxes make the elaboration of interaction priors less pertinent.

• RarePlanes [START_REF] Shermeyer | RarePlanes: Synthetic Data Takes Flight[END_REF]: synthetic optical images (reproducing satellite data) with fine plane location and type annotations. The non-oriented bounding boxes this dataset provides, and the limited interaction between neighboring planes make this dataset of limited interest to our approach.

• HRSID [START_REF] Wei | HRSID: A High-Resolution SAR Images Dataset for Ship Detection and Instance Segmentation[END_REF]: synthetic Aperture Radar (SAR) satellite images, with annotated oriented ships. This dataset provides a good challenge to test a model robustness to noise. However, the objects of interest are often isolated, thus the modeling of interaction priors is of limited interest. Still the dataset exhibits some high density samples that are worth looking into.

• VisDrone [START_REF] Zhu | Detection and Tracking Meet Drones Challenge[END_REF]: optical drone images with annotated tracks of bounding boxes. Same as for AU-AIR, the viewing angle and horizontal bounding boxes make the elaboration of interaction priors less pertinent.

Point Processes for object detection

Point Processes model the distribution of point in space [START_REF] Cox | Point Processes[END_REF][START_REF] Daley | An introduction to the theory of point processes[END_REF][START_REF] Lieshout | Markov Point Processes and Their Applications[END_REF]. These are used for a variety of tasks such as modeling the species distribution for ecological studies [START_REF] Renner | Point process models for presence-only analysis[END_REF], random tessellation for image modeling 2.2 -2.2.1 Point Process models 35 [START_REF] Bordenave | The dead leaves model: A general tessellation modeling occlusion[END_REF], building fractal-like Poisson cascades for super resolution [START_REF] Chainais | Virtual Super Resolution of Scale Invariant Textured Images Using Multifractal Stochastic Processes[END_REF], or Determinantal Point Processes (DPP), a repulsive Point Processes used for efficient sampling of a volume [START_REF] Macchi | The Coincidence Approach to Stochastic Point Processes[END_REF] used for column subset selection [START_REF] Belhadji | A determinantal point process for column subset selection[END_REF] or Kernel quadrature [START_REF] Belhadji | Kernel quadrature with DPPs[END_REF]). Here we focus on Point Process for image analysis/object detection [START_REF] Descombes | Marked point process in image analysis[END_REF][START_REF] Descombes | Stochastic geometry for image analysis[END_REF].

Point Process models

Multiple flavors of points. While simple Point Processes model the spatial distribution of points, Marked Point Processes add a random vector to each point modeling some attributes of the objects; size, angle, color, type of object etc . . . This allows modeling a variety of objects:

• Circles to model flamingos in [START_REF] Descamps | Automatic Flamingo detection using a multiple birth and death process[END_REF].

• Segments for road detection [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF], to model blood vessels (T. [START_REF] Li | An Unsupervised Retinal Vessel Extraction and Segmentation Method Based On a Tube Marked Point Process Model[END_REF], or networks in general [START_REF] Schmidt | Forest point processes for the automatic extraction of networks in raster data[END_REF].

• Rectangles in time, adding a persistent label to each object in space and time allows for tracking objects [START_REF] Crȃciun | Stochastic geometry for automatic multiple object detection and tracking in remotely sensed high resolution image sequences[END_REF].

• Contours as polygons for object detection in microscopy images [START_REF] Kulikova | A Marked Point Process Model Including Strong Prior Shape Information Applied to Multiple Object Extraction From Images[END_REF] or as a dictionary of precomputed shapes [START_REF] Descombes | Multiple objects detection in biological images using a marked point process framework[END_REF].

• Mix of processes; for instance [START_REF] Ortner | A Marked Point Process of Rectangles and Segments for Automatic Analysis of Digital Elevation Models[END_REF] use a mix of rectangles and segments to model buildings and the networks separating those.

Energy terms for the Point Process. Point Process models for object detection are defined from a density that on one side models prior on objects such as their shapes or interactions, and on the other the fitness of the point configuration against the image. Most of the approaches presented above use contrast measures. These measures are based on statistical tests between the inner and outer pixels of an object such as T-test (Student, 1908) or the Bhattacharyya distance [START_REF] Goudail | Bhattacharyya distance as a contrast parameter for statistical processing of noisy optical images[END_REF]. These measures perform well on high foreground-to-background contrast images with limited clutter in the background. Notably, (T. [START_REF] Li | Feature Extraction and Tracking of CNN Segmentations for Improved Road Detection from Satellite Imagery[END_REF] build a Point Process with contrast measures from the output of a CNN model trained for segmentation. It allows transforming the original image into a highly contrasted map to apply the contrast measure on.

Sampling the Point Process

Sampling a Point Process poses a challenge due to the non-fixed dimension of the space to sample in; as the number of points (i.e. number of objects) is a random variable too. Using the Reversible Jump Monte Carlo Markov Chain developed by [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF] and [START_REF] Geyer | Simulation Procedures and Likelihood Inference for Spatial Point Processes[END_REF] allow improving from the classical Metropolis Hasting algorithm [START_REF] Metropolis | Equation of State Calculations by Fast Computing Machines[END_REF] to jump properly in between dimensions. It involves adding a Birth and Death perturbation kernel that adds or removes one point at a time. To speed up sampling [START_REF] Descombes | Object Extraction Using a Stochastic Birth-and-Death Dynamics in Continuum[END_REF] introduce multiple births and deaths to propose more than a single point per step, later on improved with graph cuts in [START_REF] Eldin | Multiple Birth and Cut Algorithm for Multiple Object Detection[END_REF]. For more efficient fixed-dimension exploration [START_REF] Lafarge | Geometric Feature Extraction by a Multimarked Point Process[END_REF] and [START_REF] Tu | Image segmentation by data-driven Markov chain Monte Carlo[END_REF] (doing image segmentation with data driven MCMC) use Jump Diffusion from [START_REF] Grenander | Representations of Knowledge in Complex Systems[END_REF]. It uses the energy gradient to modify the current state instead of proposing perturbations independent of the energy (with a simple Gaussian perturbation for instance). This diffusion process stems from Langevin dynamics [START_REF] Welling | Bayesian learning via stochastic gradient Langevin dynamics[END_REF]. As Point Process exhibit spatial Markovianity, [START_REF] Verdié | Detecting parametric objects in large scenes by Monte Carlo sampling[END_REF] propose sampling the Point Process in parallel throughout the image, allowing to process more than one perturbation per step and making use of multiple core computing. Finally, (T. T. [START_REF] Pham | Efficient Point Process Inference for Large-Scale Object Detection[END_REF] propose to approximate the state space into a finite set, the sampling of the most likely configuration of the Point Process becomes a binary variable optimization problem

Point Process parameter estimation

As any model, Point Process models come with their set of parameters to estimate. One of the most important is the relative importance/weight of the several priors and data terms of the energy model that define the Point Process density. In their work [START_REF] Craciun | Joint Detection and Tracking of Moving Objects Using Spatio-temporal Marked Point Processes[END_REF] and (Q. [START_REF] Yu | Multiple-Target Tracking by Spatiotemporal Monte Carlo Markov Chain Data Association[END_REF] generate a set of linear constraints from generated bad configurations, and solve for the parameters with linear programming. For their model [START_REF] Chatelain | Parameter Estimation for Marked Point Processes. Application to Object Extraction from Remote Sensing Images[END_REF] use an approximation of the Expectation Maximisation algorithm [START_REF] Dempster | Maximum Likelihood from Incomplete Data Via the EM Algorithm[END_REF] to learn the Point Process parameters; using pseudo-likelihood estimation and stochastic EM (SEM) [START_REF] Celeux | On Stochastic Versions of the EM Algorithm[END_REF]. Lastly, [START_REF] Hurtut | Appearance-guided synthesis of element arrangements by example[END_REF] propose to synthesize textures (here arrangements of elements in space), from an example patch. To do that, they extract a set of perceptually meaningful statistics, and build a Point Process such that its statistics matches the observed ones by maximizing the pseudo log likelihood.

Energy Based Models

"Energy-Based Models (EBMs) capture dependencies by associating a scalar energy (a measure of compatibility) to each configuration of the variables" [START_REF] Lecun | A Tutorial on Energy-Based Learning[END_REF]. These models define a distribution over a state space from a Gibbs distribution (also known as Boltzmann distribution) as proportional to the exponential of an energy h ∝ exp(-U ).

As most of the Point Processes are defined through such a distribution, and aim to capture the dependency between a configuration of points and the image, it makes these Energy Based Models as defined above; we shall then consider methods developed in the EBMs literature for our Point Process model.

Energy Based Models as generative models

Recent interest in Energy Based Models sparks from their generative capabilities; as those define a density, one can sample new elements x ∼ exp(-U (x)). Initially generative models (e.g. image generation) were based on Generative Adversarial Networks (GAN) [START_REF] Goodfellow | Generative adversarial nets[END_REF] built from Variational Auto Encoders (VAE) [START_REF] Kingma | Auto-Encoding Variational Bayes[END_REF]. However, in recent year, popular image generation models such as [START_REF] Ramesh | Hierarchical Text-Conditional Image Generation with CLIP Latents[END_REF]) leverage diffusion models [START_REF] Ho | Denoising Diffusion Probabilistic Models[END_REF]) that progressively denoise an image into a new sample (see Figure 2.6). This denoising mimics Langevin dynamics, but instead of computing the gradient of a learned energy function, [START_REF] Ho | Denoising Diffusion Probabilistic Models[END_REF] estimate the gradient value directly.

Learning EBMs

In their work [START_REF] Lecun | A Tutorial on Energy-Based Learning[END_REF] propose to learn EBMs by maximizing their likelihood within a gradient descent scheme. However, this formulation produces some intractable integrals over the state space. While those can be approximated with Monte Carlo sampling, [START_REF] Hinton | Training Products of Experts by Minimizing Contrastive Divergence[END_REF] proposed the Contrastive Divergence method, in short taking a single sample of the Monte Carlo sampling, initialized at the ground truth to which only a few steps of the Markov Chain is applied.

Later on [START_REF] Tieleman | Training restricted Boltzmann machines using approximations to the likelihood gradient[END_REF] propose persistent divergence, initializing this contrastive sample with the sample obtained at the previous iteration of the gradient descent. We find this method used for image generation in [START_REF] Du | Improved Contrastive Divergence Training of Energy-Based Models[END_REF] and [START_REF] Song | Generative Modeling by Estimating Gradients of the Data Distribution[END_REF]. We will go further in details about these methods in Section 5.2.2.

EBM applications to computer vision

Generative models can seem quite remote from our object detection tasks. However, we see some interesting works proposing to use the composable nature of EBMs2 such as (R. [START_REF] Zhang | Robust and Controllable Object-Centric Learning through Energy-based Models[END_REF] performing unsupervised object discovery and controllable scene manipulation. EBMs are intrinsically bound to classification task, as [START_REF] Grathwohl | Your classifier is secretly an energy based model and you should treat it like one[END_REF]

PART

Model foundations CHAPTER 3

Foundations for Point Processes and Convolutional Neural Networks

Ce chapitre introduit les notions fondamentales sur les Processus Ponctuels et Réseaux de Neurones à Convolution. Celles-ci sont nécessaires à la construction de nos modèles de détection dans les chapitres suivants. Pour plus de détails sur les Processus Ponctuels, le lecteur peut se référer à [START_REF] Stoyan | Stochastic geometry and its applications[END_REF][START_REF] Lieshout | Markov Point Processes and Their Applications[END_REF], dont est issu en grande partie ce chapitre. Par la suite, on décrit comment modéliser la détection d'objets avec des Processus Ponctuels et les méthodes nécessaires pour échantillonner ces derniers. Enfin, sont présentés les outils relatifs aux réseaux de neurones à convolution qui seront utilisés dans nos contributions.

In this chapter we introduce the fundamental notions for Point Processes and Convolutional Neural Networks. These are needed to build our models in the later chapters. For more details on Point Processes please refer to [START_REF] Stoyan | Stochastic geometry and its applications[END_REF][START_REF] Lieshout | Markov Point Processes and Their Applications[END_REF] from which most of this chapter is sourced. We describe how to model object detection with Point Processes, and the necessary methods for sampling the latter. Finally, we define the necessary tools of Convolutional Neural Network needed for our contributions.

CHAPTER 3 -Foundations for Point Processes and Convolutional Neural Networks

Poisson Point Processes

We first introduce the fundamentals of Point Process. For more details please refer to [START_REF] Stoyan | Stochastic geometry and its applications[END_REF][START_REF] Lieshout | Markov Point Processes and Their Applications[END_REF].

We consider configurations of points in the space S ⊆ R d (e.g. a square in R 2 ). In the rest of this chapter, we associate a metric d to the space S; e.g. the Euclidean distance in R d . Definition 3.1.1. We call configuration, denoted y, a finite unordered set of points y k , k = 1, . . . , n(y) in S, with n(y) denoting the number of points in y: y = {y 1 , . . . , y n(y) }.

(3.1)

The set of all possible configurations of points in S with any number of points is denoted:

Y = ∞ ⋃ n=0 Y n , ( 3.2)
where

Y 0 = ∅ and Y n = {{y 1 , . . . , y n }, y i ∈ S, ∀i}.
A configuration y is said to be locally finite if for any bounded Borel set A ⊆ S, the number of points of y in A is finite. The set of locally finite configurations is denoted N lf . Definition 3.1.2. A Point Process (PP) on S is a mapping Φ from a probability space (Ω, A, P) into N lf such that for any bounded Borel set A ⊆ S. The number of points that fall in A, denoted N Φ (A), is a finite random variable.

Practically, a Point Process is a random variable the realization of which is a configuration of points. If the space S is bounded or N Φ (A) is almost surely finite, we call the Point Process finite.

For most Point Processes, a point pattern will not contain multiple points at exactly the same location, either because it is impossible, or the multiplicity might be encoded within a mark (see Definition 3.1.8); i.e. N Φ ({y}) ∈ {0, 1} for all y ∈ S. The set of locally finite configurations composed of distinct points is denoted N lf s .

Definition

1. N Φ (A) is a Poisson random variable of mean ν(A) 2. For k disjoint Borel sets A 1 , . . . , A k , random variables N Φ (A 1 ), . . . , N Φ (A k ) are in- dependent.
The Poisson Point Process is called homogeneous (see Figure 3.1 (a)) when ν is of the form λ|A|, with |A| the Lebesgue measure of A ⊆ S and λ a strictly positive scalar. The intensity measure ν can be used to control the point density over the space S as shown in Figure 3.1 (b).

CHAPTER 3 -Foundations for Point Processes and Convolutional Neural Networks

The law of a Poisson Point Process on intensity measure ν in a window S ⊂ R d is defined as

µ(A) = ∞ ∑ n=0 ν(S) n exp(-ν(S)) n! ⏞ ⏟⏟ ⏞ pn ∫ y 1 ∈S . . . ∫ y n ∈S ⏞ ⏟⏟ ⏞ n times 1 A ({y 1 , . . . , y n }) 1 ν(F ) n dy 1 . . . dy n (3.3) µ(A) = ∞ ∑ n=0 exp(-ν(F )) n! ∫ S n 1 A (y)ν(dy 1 ) . . . ν(dy n ), (3.4)
where A is part of the sigma-algebra B from Y.

Markov Point Processes

While this modeling of points in space accounts for a wider variety of distributions, there is still no consideration for point interactions. We can build a wider range of Point Process models by means of their probability density with respect to the Poisson process. Here we consider a metric space (S, d) and the Poisson Point Process on S with density µ, and intensity measure ν. The probability density f of a Point Process w.r.t. the law µ of a Poisson Point Process is a mapping from the configuration space Y to [0, ∞[ such that:

f : Y ↦ → [0, ∞[, ∫ Y f (y)µ(dy) = 1. (3.5)
With density f w.r.t. the measure of a Poisson process, we can now model interactions. While the density function can take into account interactions between any points in the configuration, real life spatial processes often only interact locally. This characteristic simplifies the modeling of the Point Process, and yields some properties that are useful when simulating it. The locality of the Point Process is formalized through Markov Point Processes: Definition 3.1.5. We define the neighborhood ∂ S A of a set A ⊆ S in S for reflexive and symmetric relation ∼ as:

∂ S A = {x ∈ S : x ∼ a for some a ∈ A}. (3.6) 3.1 -3.1.2 Markov Point Processes
By extension the set of neighbors of a single point y is ∂ S {y} , and the set of points of configuration y neighboring y is ∂ y {y} .

From the definition of neighborhood, we can define a Markov Point Process: Definition 3.1.6. Let Φ be a Point Process of density f . Point Process Φ is a Markov Point Process under the symmetric and reflexive relation ∼ if and only if, for every configuration y ∈ Y such that f (y) > 0:

1. ∀x ⊂ y, f (x) > 0 2.
For every point u ∈ S, f ({u} ∪ y)/f (y) only depends on u and its neighborhood in y, ∂ y {u} .

The Hammersley-Clifford theorem allows the density of the Markov process to be decomposed as the product of local functions defined on cliques; a clique being a subset of y where all points are neighbors of each other.

Theorem 3.1.1. A Point Process density f : Y ↦ → [0, ∞[ is Markov w.r.t. the relation ∼ if and only if there is a measurable function Ψ : Y ↦ → [0, ∞[ such that: f (y) = ∏ y∈Cy Ψ(y), (3.7)
where C y is the set of cliques

C y = {y ⊆ y : ∀{y, y ′ } ⊆ y, y ∼ y ′ }
In practice, we often write f as the normalized version of an unnormalized function h:

f (y) = h(y) Z , ( 3.8)
where Z = ∫ Y f (y)dy. Markov Point Processes are also known as Gibbs Point Process, since (3.7) can be written in an energy form, where V (y) represents the potential of a clique :

f (y) ∝ h(y) = exp ⎛ ⎝ - ∑ x∈Cy V (x) ⎞ ⎠ , ( 3.9)
with ∝ standing for "proportional to" (skipping the normalizing factor ∫ Y h(y)dy). The Markov property of a Point Process comes in useful when the above-mentioned normalizing factor becomes intractable. Indeed, the second property in definition 3.1.6 allows the normalizing constant to be canceled out; as long as the function h in integrable over Y, we can skip computing its integral. The density ratio f ({u} ∪ y)/f (y) is in fact the Papangelou intensity, which we now define: Definition 3.1.7. The Papangelou intensity λ(•; •) associated to a simple Point Process Φ, can be interpreted as:

λ(y; y)dy = p (N Φ (dy) = 1|Φ ∩ (dy) c = y ∩ (dy) c ) , ( 3.10)
i.e. the infinitesimal probability to find a point in region dy around y ∈ S, given the configuration y outside dy ((dy) c being the complement of dy in S).

CHAPTER 3 -Foundations for Point Processes and Convolutional Neural Networks

We then have the following:

Theorem 3.1.2. For a finite Point Process Φ specified by a density f w.r.t. a Poisson Point Process with intensity measure ν, the Papangelou conditional intensity is given as: .11) for u / ∈ y.

λ(u; y) = f ({u} ∪ y) f (y) , ( 3
Example: pairwise interaction process. A typical Markov Point Process is the pairwise interaction process, for which the density function is

h(y) = 1 z ∏ y∈y β(y) ∏ u,v∈y, u∼v γ(u, v), (3.12)
where z is the normalizing constant, β : S ↦ → [0, ∞[ the intensity function and γ : S ×S ↦ → [0, ∞[ the pair interaction function. The Strauss Point Process is one instance of pairwise interaction process. Its density is given as:

h(y) ∝ β n(y) γ s(y) , ( 3.13)
where β > 0 is an intensity parameter, and s(y) the number of pairs of points in y that are at a distance r or less apart. Depending on the value of parameter γ, the process exhibits different behaviors illustrated in Figure 3.2 :

• γ = 0: points cannot be closer than distance r. This is called a hardcore process.

• 0 < γ < 1: the points are repulsive to each other.

• γ = 1: this is the Poisson process of intensity βλ.

• γ > 1: the density is not integrable, the Point Process is not defined. Introducing a limit on the number of points n max such that h(y) ∝ 1{n(y) ≤ n max }β n(y) γ s(y) , allows defining this process for which the points are attractive.

Markov marked Point Processes

Marked Point Processes are Point Processes for which a mark is attached to each point. As such, they are useful for applications where each point or event has some properties attached to it. For instance marks can encode the radius of a point, or combining radius with height: a cylinder. The mark can also carry semantic information: the type of object (e.g. species of trees), or multiplicity (if one point can represent several events). All these descriptors can be combined to enrich the Point Process. Let Φ be a marked Point Process with positions in S and marks in M specified by a density h w.r.t. f ν×η . Then Φ is a Markov marked Point Process w.r.t. the symmetric reflexive relation ∼ on S × M if for all y such that f (y) > 0,

1. f (x) > 0 for all x ⊆ y; 2. for all (u, l) ∈ S × M, f (y ∪ {(u, l)})/f (y) depends only on (u, l) and its neighbor- hood ∂ y {(u,l)} .
The Hammersley-Clifford theorem remains valid for Markov marked Point Processes.

Example: Marked pairwise interaction process. We consider a marked pairwise interaction process on S × M, with S a compact window in R 2 and M = [0, r max ] a mark representing a radius, bounded by r max :

h(y) = ∏ (y,k)∈y β ∏ (u,k),(v,l)∈y γ(||u -v|| -k -l), (3.14)
for some intensity β k > 0 and measurable interaction function

γ : [0, ∞[↦ → [0, ∞[. If γ(r) = 1
for r > r 0 for some permeability distance r 0 else γ(r) = 0, the process is Markov with respect to the relation

(u, k) ∼ (v, l) ⇔ ||u -v|| -k -l ≤ r 0 . (3.15)
We show the effect of the permeability parameters r 0 on some realizations in Figure 3.3.

Stability conditions

As we have seen previously we can build the density of the Point Process f from an unnormalized density h. This density h needs to be normalizable w.r.t. the reference Point Process. While this can prove tedious, the Ruelle condition [START_REF] Ruelle | Superstable interactions in classical statistical mechanics[END_REF] allows us to simply bound the density to make sure it is normalizable. This condition is enough to ensure that h can be normalized:

∫ y∈Y h(y)dµ(y) ≤ ∞ ∑ n=0 M n ν(S) n! = exp(M ν(S)). ( 3

Point Process for object detection

Up to now, we settled the fundamentals of Point Processes, however our goal remains to perform object detection in an image. How can we use the Point Process to model our objects both with their geometrical properties, and their relation to the image ? First, we define configurations such that they represent correctly the objects of interest, where the point position informs on the location of the object (in space, but also in time if necessary), and marks add geometrical or semantic information. Tree crowns can be represented as Point Processes of circles [START_REF] Perrin | Tree crown extraction using marked point processes[END_REF], road networks as segments [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF] or buildings as rectangles [START_REF] Ortner | A Marked Point Process of Rectangles and Segments for Automatic Analysis of Digital Elevation Models[END_REF].

Priors on the interactions (and on the geometry of object themselves) are added; penalizing overlaps, forbidding some configurations and favoring others. We denote h priors (y) the priors' density, also referred as internal energy when using potentials.

Then we take into account the image data. The Bayesian approach is to compute the likelihood of the observation (or image) X given the configuration y [START_REF] Rue | Bayesian object identification[END_REF][START_REF] Perrin | Tree crown extraction using marked point processes[END_REF]: given the silhouette of the Point Process s ρ (y) at ρ. However, such approach requires both a foreground and background model to compute g. Thus, in cases where the background proves to be more complex (presence of other objects, clutter etc...), the approach is to use data terms for each object of configuration y as an external field. For instance one can evaluate the homogeneity of the image within the object boundary and its contrast with nearby background [START_REF] Descombes | Multiple objects detection in biological images using a marked point process framework[END_REF][START_REF] Schmidt | Forest point processes for the automatic extraction of networks in raster data[END_REF][START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF][START_REF] Ortner | Processus ponctuels marqués pour l'extraction automatique de caricatures de bâtiments à partir de modèles numériques d'élévation (These de doctorat[END_REF]:

h data (y) = h data (X|y) = ∏ pP X g(X[ρ]|s ρ (y)), ( 3
h data (y) = exp ⎛ ⎝ - ∑ y∈y V data (y, X) ⎞ ⎠ . (3.20)
Each potential gets lower as the object better fits the image. With this non-Bayesian approach we risk overlapping objects by accumulating low potentials, thus high penalization on overlaps is needed.

Finally, we get:

f (y) ∝ h data (y)h priors (y). (3.21)
In the Bayesian case (image model using the likelihood (3.19)) we get the posterior density f (y|X) instead of f (y).

At inference, one needs to sample the configuration y * that maximizes the density h(y) (or minimizes the energy U for a Gibbs model):

y * = arg max y∈Y h(y).
(3.22)

Energy model

In the following, we define our density h through an energy U (y, X) function of the configuration y and image X such that h(y) = exp(-U (y, X)). One way to build U is to define potentials for each object interactions [START_REF] Ortner | Processus ponctuels marqués pour l'extraction automatique de caricatures de bâtiments à partir de modèles numériques d'élévation (These de doctorat[END_REF][START_REF] Crȃciun | Stochastic geometry for automatic multiple object detection and tracking in remotely sensed high resolution image sequences[END_REF]:

U (y, X) = ∑ y∈y V ext. (y, X) + ∑ y,y ′ ∈y y∼y ′ v int. (y, y ′ ). (3.23)
In our work we rather use the following model, the choice of which we explain in Section 4.1:

U (y, X) = ∑ y∈y V ( y, X, N y {y} ) (3.24) = ∑ y∈y ∑ e∈ξ V e ( y, X, N y {y} ) , (3.25)
where V is the energy of a point composed of energy terms V e , ξ being the set of energy terms. We purposefully introduce a new neighborhood notation N y {y} as opposed to ∂ y {y} , as with this model, the object interaction relation and Markovianity relation may differ.

Markovianity relation

We denote ∼ (3.26)

h(y ∪ {u}) h(y) = exp (U (y, X) -U (y ∪ {u})) = exp ⎛ ⎝ ∑ y∈y [ V ( y, X, N y {y} ) -V ( y, X, N y∪{u} {y} )] -V ( u, X, N y∪{u} {u} ) ⎞ ⎠ = exp ⎛ ⎜ ⎝ ∑ y∈N y {u} [ V ( y, X, N y {y} ) -V ( y, X, N y∪{u} {y} )] -V ( u, X, N y {u} ) ⎞ ⎟ ⎠ using the fact that y / ∈ N y {u} =⇒ N y {y} = N

Local stability

To build converging simulations of our Point Process, we need to check for Condition 3.1.2; i.e. we want to find M ∈ R such that h(y ∪ {u})/h(y) ≤ M . We suppose all V e , e ∈ ξ are built such that:

∃A > 0 |V e (y, X, N y {y} )| < A, ∀y ∈ y, ∀y ∈ Y, ∀e ∈ ξ. (3.27)
Thus, |V (y, X, N y {y} )| < A|ξ|. From the previous results on the Markovianity relation we have that: .28) We are left to find an upper bound for the number of neighbors. [START_REF] Ortner | Processus ponctuels marqués pour l'extraction automatique de caricatures de bâtiments à partir de modèles numériques d'élévation (These de doctorat[END_REF]) use an exclusion energy that effectively forbids configurations with a number of neighbors above a certain threshold. In our case we use a limit on objects per cell (see Section 3.3.2.4 and 5.1.2 on parallelization) n c,max . In turn this sets the maximum number of neighbors a point can get up to 4n c,max -1 (the worst case being a point in the corner of a cell with neighbors in the 3 neighboring cells of that corner, as illustrated in Figure 3.4). Finally, we show we meet the Condition 3.1.2 as such:

h(y ∪ {u}) h(y) ≤ exp ( 2A|ξ| ⃓ ⃓ ⃓N y {y} ⃓ ⃓ ⃓ + A|ξ| ) . ( 3
h(y ∪ {u}) h(y) ≤ M, M = exp (A|ξ|(8n c,max -1)) . (3.29)
3.3 Sampling the Point Process

Markov chains

To sample configurations from the Point Process, we need to build a Markov chain that converges towards the Point Process distribution.

(Y n ) is a Markov chain if p(Y t+1 ∈ A|Y t = y t , . . . , Y 0 = y 0 ) = p(Y t+1 ∈ A|Y t = y t ), ∀A ∈ B.
(3.30)

In short, the evolution of the Markov chain only depends on the current state. We will consider homogeneous Markov chains, defined as such: Definition 3.3.2. A Markov chain is homogeneous if the evolution does not depend on the time parameter t. Definition 3.3.3. A transition kernel is a function K defined on Y × B such that:

• ∀x ∈ Ω, K(x, •) is a probability measure. • ∀A ∈ B, K(•, A) is measurable. The transition kernel K for Markov chain (Y t ) t∈N is such that K(y, A) = p(Y t+1 ∈ A|Y t = y).
(3.31)

The Markov chain needs to meet some other properties in order to use it for sampling:

• Stationary: a measure π is stationary for Markov chain (Y t ) with transition kernel K if:

π(A) = ∫ K(y, A)π(dy), ∀A ∈ B. (3.32)
• Reversible: a Markov chain is reversible if its transition kernel K is such as

∫ A K(y, B)π(dy) = ∫ B K(y ′ , A)π(dy ′ ), ∀A, B ∈ B. (3.33)
• Irreducible: a Markov chain is irreducible if there is a positive probability to reach any y from any other y ′ in a finite number of steps:

∀y, y ′ ∈ Y, ∃k < ∞, p(Y t+k = y|Y t = y ′ ) > 0. (3.34)
• Aperiodic: a Markov chain is aperiodic if it exhibits no cycles, i.e.:

GCD {d : p(Y d = y|Y 0 = y)} = 1, (3.35)
where GCD is the Greatest Common Divisor.

• Harris recurrent: a Markov chain is Harris recurrent, if for any A ∈ B such that π(A) > 0, p(∃t :

y t ∈ A|Y 0 = y) = 1, ∀y ∈ Y. (3.36)
• Ergodic: an aperiodic Harris recurrent Markov chain is ergodic and converges towards π:

⃦ ⃦ ⃦K t (x, •) -π(•) ⃦ ⃦ ⃦ ---→ t→∞ 0, ∀x ∈ Y, (3.37) with ∥•∥ such as ∥µ 1 -µ 2 ∥ = sup A |µ 1 (A) -µ 2 (A)|.

Reversible

Jump Monte Carlo Markov chain [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF] build this ergodic Markov chain with a mixture of perturbation kernels

Q m (y, A), m ∈ Q such that : Q(y, A) = ∑ m∈Q Q m (y, A). ( 3
y ′ ∼ Q m (y, •) Q m (y, Y) . (3.39)
3. Compute the Green ratio r(y, y ′ ), see (3.47).

4. Accept perturbation from y to y ′ with probability α(y, y ′ ) = min(1, r(y, y ′ ))

As such the transition kernel K for the Markov chain is a function of the perturbation kernel Q and the acceptance probability α. Thus, to be able to perform such a procedure, each kernel Q m needs to have the following properties:

1. We must know the probability of picking

Q m ; Q m (y, Y).
2. The probabilities of picking each kernel sum to at most 1;

∑ m∈Q Q m (y, Y) ≤ 1, ∀y ∈ Y 3.
We can sample state perturbations from the normalized kernel density Q m (y, .)/Q m (y, Y)

for any configuration y ∈ Y. The last condition above is here to ensure the reversibility condition of the Markov Chain is met (implying that π is the stationary measure of the Markov Chain Y t ). For the chain to be reversible, its probability to go from a set A ∈ B towards a set B ∈ B needs to be the same at going from B to A (see (3.33)). This is known as detailed balance (DB) condition. As the kernel Q is a sum of kernels Q m , we only need to show DB is maintained for each kernel Q m [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF][START_REF] Ortner | Processus ponctuels marqués pour l'extraction automatique de caricatures de bâtiments à partir de modèles numériques d'élévation (These de doctorat[END_REF]. The transition sub-kernel for

Q m is K m (y, A) = ∫ A Q m (y, dy ′ )α(y, y ′ ). (3.41)
To verify DB we then need :

∫ A ∫ B Q m (y, dy ′ )α(y, y ′ )π(dy) = ∫ A ∫ B Q m (y ′ , dy)α(y ′ , y)π(dy ′ ). (3.42)
To find the function α that verifies this equation, [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF] proposes to find a symmetric measure ψ m on Y × Y, such that πQ m is absolutely continuous w.r.t. ψ m ; i.e.:

ψ m (A, B) = 0 =⇒ πQ m (A, B). (3.43)
The Radon-Nikodym theorem then tells us a unique function D m exists called the Radon-Nikodym derivative (or density) of πQ m w.r.t. ψ m , such that :

∫ A Q m (y, B)π(dy) = ∫ A ∫ B D m (y, y ′ )ψ m (dy, dy ′ ). (3.44)
Then we can replace πQ m in (3.42) to get:

∫ A ∫ B α(y, y ′ )D m (y, y ′ )ψ m (dy, dy ′ ) = ∫ A ∫ B α(y ′ , y)D m (y ′ , y)ψ m (dy ′ , dy). (3.45)
Since we built ψ m to be symmetric, we only need:

α(y, y ′ )D m (y, y ′ ) = α(y ′ , y)D m (y ′ , y).
(3.46) [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF] proposes to take:

α(y, y ′ ) = min { 1, D(y, y ′ ) D(y ′ , y) } .
(3.47)

Balancing in terms of random numbers

In more recent works [START_REF] Green | Reversible jump MCMC[END_REF] propose an alternative derivation of α for the Markov Chain in order to maintain detailed balance. The perturbation is defined as follows: From current configuration y, choose deterministic function τ m with probability j m (y). It generates the new proposed state x as (x, δ ′ ) = τ m (y, δ), where δ is an r-dimentional random number to generate x from y, and δ ′ the r ′ -dimentional random number used to generate y from x via τ ′ , the reverse transform from τ m (i.e. τ ′ m (x, δ ′ ) = (y, δ)). Random numbers δ and δ ′ are respectively generated with known densities g and g ′ . Move from y to x accepted with probability α(y, x). The Detailed Balance is them given as:

∫ y,x∈A×B h(y)j m (y)g(δ)α(y, x)dydδ = ∫ y,x∈A×B h(x)j m (x)g ′ (δ ′ )α(x, y)dxdδ ′ . (3.48) If transform τ m from (y, δ) to (x, δ ′ ) and τ ′
m are differentiable, we can apply the change-ofvariable formula to the right-hand side of the equation. The equation now holds if

h(y)j m (y)g(δ)α(y, x) = h(x)j m (x)g ′ (δ ′ )α(x, y) ⃓ ⃓ ⃓ ⃓ ∂(x, δ ′ ) ∂(y, δ) ⃓ ⃓ ⃓ ⃓ .
(3.49)

A valid choice for α is then:

α(y, x) = min { 1, h(x)j m (x)g ′ (δ ′ ) h(y)j m (y)g(δ) ⃓ ⃓ ⃓ ⃓ ∂(x, δ ′ ) ∂(y, δ) ⃓ ⃓ ⃓ ⃓ } , (3.50) with ⃓ ⃓ ⃓ ∂(x,δ ′ ) ∂(y,δ) ⃓ ⃓
⃓ the determinant of the Jacobian matrix of τ m . This approach works across dimensions, provided that the transformation from (y, δ) to (x, δ ′ ) remains a diffeomorphism [START_REF] Green | Reversible jump MCMC[END_REF]. Denoting with n and n ′ the dimensions y and x respectively, the diffeomorphism requires that n + r = n ′ + r ′ . Note that either one or both of r or r ′ might be 0.

Perturbation kernels

Here we (re)build the basic perturbation kernels for the RJMCMC. Note that the birth and death kernel is the minimal requirement for the Markov chain to converge to the stationary distribution. The demonstration for the birth and death kernel is sourced from [START_REF] Geyer | Simulation Procedures and Likelihood Inference for Spatial Point Processes[END_REF][START_REF] Lacoste | Extraction de réseaux linéiques à partir d'images satellitaires et aériennes par processus ponctuels marqués (These de Doctorat[END_REF][START_REF] Ortner | Processus ponctuels marqués pour l'extraction automatique de caricatures de bâtiments à partir de modèles numériques d'élévation (These de doctorat[END_REF][START_REF] Crȃciun | Stochastic geometry for automatic multiple object detection and tracking in remotely sensed high resolution image sequences[END_REF].

Birth and Death. The birth and death kernel adds and removes point in the configuration. It is essential as it allows moving across dimensions in Y. This kernel itself is a mix of two sub kernels: a birth kernel (that adds new points) and a death kernel (that removes points). The most basic birth and death kernel proposes new points y with density ν(•)/ν(S) and removes points by picking one uniformly in y. The resulting kernel is as follows, with p B and p D the probabilities to pick birth or death respectively:

Q(y, y ′ ) = p B (y)Q B (y, y ′ ) + p D Q D (y, y ′ ) (3.51)
Both kernels are defined as:

Q B (y, A) = ∫ y∈S 1 A (y ∪ {y}) ν(dy) |S| , (3.52) Q D (y, A) = ∑ y∈y 1 A (y \ {y}) 1 n(y) . ( 3
Q(y, B n+1 )π(dy) = p B ∫ An (∫ y∈S 1 B n+1 (y ∪ {y}) ν(dy) |S| ) f (y)µ(dy) (3.54) ∫ B n+1 Q(y, A n )π(dy) = p D ∫ B n+1 ⎛ ⎝ ∑ y∈y 1 An (y \ {y}) 1 n(y) ⎞ ⎠ f (y)µ(dy) (3.55) Splitting ψ n into ψ + n and ψ - n respectively defined on Y n × Y n+1 and Y n+1 × Y n , we define ψ + n (A n , B n+1 ) = ∫ An ∫ y∈S 1 B n+1 (y ∪ {y})ν(dy)µ(dy) (3.56) ψ - n (B n+1 , A n ) = ∫ B n+1 ∑ y∈y 1 An (y \ {y})µ(dy) (3.57)
We show the symmetry of ψ + n and ψ - n by decomposing the Poisson process measure µ using (3.4):

ψ - n (B n+1 , A n ) = exp(-ν(S)) (n + 1)! ∫ S n+1 ∑ y∈y 1 An (y \ {y})1 B n+1 (y)ν n+1 (y) = exp(-ν(S)) (n + 1)! ∫ S n+1 (n + 1)1 An ({y 1 , . . . , y n })1 B n+1 ({y 1 , . . . , y n+1 })ν n+1 (y) = exp(-ν(S)) n! ∫ S n ∫ S 1 B n+1 (x)1 An (x ∪ {u})ν n (dx)ν(du) = ψ + n (A n , B n+1 ).
Then we can define ψ on Y × Y; more specifically on we define ψ on

⋃ ∞ n=0 {{Y n × Y n+1 } ∪ {Y n+1 × Y n }}: On each {Y n × Y n+1 } ∪ {Y n+1 × Y n }, ψ is equal to ψ n . We have ψ(A, B) = 0 =⇒ πQ(A, B) = 0.
And the Radon-Nikodym derivatives given for birth and death (deduced respectively from (3.54), (3.56) and (3.55), (3.57)), are respectively:

D(y, y ∪ {u}) = p B f (y) ν(S) , (3.58) D(y, y \ {u}) = p D f (y) n(y) . (3.59)
The associated Green ratios for birth and death are then

r(y, y ∪ {u}) = p D p B ν(S) n(y) + 1 h(y ∪ {u}) h(y) , (3.60) r(y, y \ {u}) = p B p D n(y) ν(S) h(y \ {u}) h(y) , (3.61)
with h the unnormalized density of the process. Local perturbation kernel. This perturbation does not change the dimension of the current state, the number of points is fixed. It modifies the parameters of an object, for instance translating, rotating or scaling as illustrated in Figure 3.5. While [START_REF] Ortner | Processus ponctuels marqués pour l'extraction automatique de caricatures de bâtiments à partir de modèles numériques d'élévation (These de doctorat[END_REF][START_REF] Crȃciun | Stochastic geometry for automatic multiple object detection and tracking in remotely sensed high resolution image sequences[END_REF] use the symmetric measure approach to derive acceptance probability α, we here use the formalism from [START_REF] Green | Reversible jump MCMC[END_REF] described in Section 3.3.2.2 to show the same result.

We consider a Marked Point Process with dim(M) marks (e.g. 3 for a rectangle or ellipse) and denote d = dim(M) + dim(S) the dimension of each point (e.g. d = 5 for rectangles in 2D space). We use the random number formalism as described in Section 3.3.2.2. The kernel to transition from y to x works as follows:

1. The transform kernel is picked with probability p T .

2. The k th object y is picked uniformly in y 3. Generate a perturbation δ ∈ R d+2 4. Proposed configuration is defined as:

(x, δ ′ ) = τ k (y, δ) = ({y 1 , . . . , y k + δ, . . .}, δ), (3.62) thus τ k (x, δ ′ ) = ({x 1 , . . . , x k -δ ′ , . . .}, δ ′ ): then δ ′ = δ.
Here we consider a set of n(y) sub kernels, each one modifying the k th object. Each sub-kernel is picked with probability j k (y) = j k (x) = p T n(y) . The density of random variable δ is g(δ) and the one for δ). Then Equation (3.50) gives us:

δ ′ is g ′ (δ ′ ) = g ′ (δ) = g(-
α(y, x) = min { 1, h(x)g(-δ) h(y)g(δ) ⃓ ⃓ ⃓ ⃓ ∂(x, δ ′ ) ∂(y, δ) ⃓ ⃓ ⃓ ⃓ } (3.63)
The Jacobian is equal to 1, as shown in Appendix A. Thus:

α(y, x) = min { 1, h(x)g(-δ) h(y)g(δ) } . (3.64)
For α to be non-zero for a given δ, we need g(-δ) > 0; i.e. there exists a way to reverse the move defined by δ. The basic approach is to pick a density symmetric around zero, such that g(δ) = g(-δ) (e.g. Gaussian distribution centered in 0); it also simplifies the acceptance probability. The above-described approach is still valid if we pick the object y to transform in y with a non-uniform distribution, requiring to change j k (y) and j k (x) accordingly. From the definition of a Markov Point Process (Definition 3.1.6) and the Green ratio (e.g. birth (3.60)), we see that the ratio h(y ′ )/h(y) only depends on added point y and its neighborhood according to ∼. [START_REF] Verdié | Detecting parametric objects in large scenes by Monte Carlo sampling[END_REF] show that two sequential perturbations for which the neighborhoods do not intersect can be processed in any order (the resulting Green ratio from each move not depending on the other ones), thus can be done in parallel (see Figure 3.6).

For a Point Process derived from an energy U (y) =

∑ y∈y V 1 (y) + ∑ y ′ ∼y V 2 (y, y ′ )
, two cells c and c ′ are independent if and only if:

min u∈c,y∈c ′ ∥u -y∥ 2 ≥ d max + 2δ max , ( 3.65)
with d max the interaction radius between points (for ∼) and δ max the maximum translation distance performed by any kernel. Each set of mutually independent cells is called a mic-set. While it works with a regular grid partitioning of space, [START_REF] Verdié | Detecting parametric objects in large scenes by Monte Carlo sampling[END_REF]) use a data driven quad tree structure to efficiently propose points for non-uniform distributions.

The general transition kernel Q is now a mixture of the usual sub-kernels, restricted to a cell c.

Kernel Q m restricted to c is denoted Q c,m : ∀y ∈ Y, Q(y, •) = ∑ c ∑ m∈Q p c,m Q c,m (y, •), (3.66)
with p c,m the probability to pick cell c and kernel m.

Jump diffusion

Intra-dimensional kernels (i.e. kernels that operate at a fixed dimension) such as the local perturbation kernel (see Section 3.3.2.3) propose transformations randomly to explore nearby configuration of points. However, the transformations are done with no consideration of the energy landscape; in (3.64) the perturbation proposal density g is not tied to the current configuration y or density h. In Section 5.1.1 we propose a data driven local transform kernel inspired from [START_REF] Ortner | Processus ponctuels marqués pour l'extraction automatique de caricatures de bâtiments à partir de modèles numériques d'élévation (These de doctorat[END_REF][START_REF] Descombes | Stochastic geometry for image analysis[END_REF]. It leverages a pre-computed position energy map to sample perturbation δ. However, this proves costly as it requires producing and sampling from a distribution derived from the position energy map around the current point y (which needs to be extracted from the whole map, and normalized locally). Moreover, it only takes into account a part of the energy model; interactions are not taken into account, thus it might propose perturbations δ that result in higher energies due to overlaps. Ideally we would use the local gradient of the energy to update y: on one hand, it avoids building a density from a large energy map around the current point at each iteration; the computation is local. On the other we consider the whole energy model, not just a subset of it.

Jump-diffusion is introduced by [START_REF] Grenander | Representations of Knowledge in Complex Systems[END_REF] and was successfully applied to object detection and tracking [START_REF] Srivastava | Jump-diffusion Markov processes on orthogonal groups for object pose estimation[END_REF], segmentation [START_REF] Han | Range image segmentation by an effective jump-diffusion method[END_REF], or geometric object extraction [START_REF] Lafarge | Geometric Feature Extraction by a Multimarked Point Process[END_REF] since then.

Jump-diffusion mixes Langevin dynamics [START_REF] Geman | Diffusions for Global Optimization[END_REF]) (also mentioned as diffusion) and the MCMC [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF]. The diffusion explores at fixed number of points, i.e. Y n , while jumps, passes between subspaces of Y.

The jump process consists of the birth and death kernel of the RJMCMC. Kernels such as the local transform kernel (Section 3.3.2.3) are no longer needed as the diffusion fulfills the same purpose.

The diffusion process differs from the usual MCMC procedure as it does not involve an accept/reject mechanism. For a Point Process such that h(y) ∝ exp(-U (y)), the following diffusion equation converges towards π(dy):

dy t = -∇U (y t )dt + W t , (3.67)
where W t is a Brownian motion. This equation is defined over continuous time, the discrete diffusion, with t ∈ N is given by:

y t+1 = y t + -γ∇U (y t ) + √ 2T (t)w t , w t ∼ N (0, γ), (3.68)
where γ is the step size and T (t) serves as a relaxation temperature; at higher T the Brownian motion allows exploring more of Y n , while at low T the Brownian motion becomes negligible and the diffusion dynamics acts as gradient descent. While the jump-diffusion is usually more effective than the standard MCMC sampler [START_REF] Descombes | Stochastic geometry for image analysis[END_REF], it requires being able to compute the gradient of the energy U at any configuration y, with the density function being Lipschitzcontinuous [START_REF] Geman | Diffusions for Global Optimization[END_REF]).

Definition 3.3.4. A real valued function f is Lipschitz-continuous if there exists a positive real constant k such that for all real x 1 and x 2 :

|f (x 1 ) -f (x 2 )| ≤ k|x 1 -x 2 |.
(3.69)

Simulated annealing

The RJMCMC procedure described here will -given the proper conditions -run a Markov chain for which the stationary measure is π. In order to get the best fitting configuration y * (see (3.22)) -the configuration that maximizes the posterior density in the case of a Bayesian model -we use simulated annealing to gradually collapse the Markov chain stationary distribution into the single global maximum of h. Simulated annealing consists in running a Markov chain (y t ) t∈N with density h t (y) = h 1/Tt (y) as the temperature T t decreases to 0. Simulated annealing is loosely based on annealing in metallurgy, where the heating and cooling of a metal can alter its physical properties.

-3.3.5 Stopping conditions 59

At high temperature the chain can explore more of the configurations space Y, the stationary density gets closer to the uniform density over. At temperatures near 0, only perturbations that increase the density h are accepted. The annealing scheme starts at high temperature to be able to reach all states of Y, then the temperature slowly decreases (as to maintain stability and explore different modes), as the chain reaches the global maximum only small steps are accepted until the stopping condition is reached.

The chain is proved to converge to the global maximum for Reversible Jump MCMC [START_REF] Stoica | Simulated annealing and object point processes: Tools for analysis of spatial patterns[END_REF] and birth and death process [START_REF] Robert | Diagnosing Convergence[END_REF], but for carefully chosen decrease function.

Logarithmic decrease.

T t = C log(t + 1)
. (3.70) This decrease scheme ensures the convergence if the constant C is larger than the deepest local maximum. In practice this decrease scheme is very slow, often time logarithmic decrease is replaced by geometric decrease.

Geometric decrease. Given as:

T t+1 = αT t , (3.71)
with α ∈ [0, 1] driving the speed of the decrease; usually α is picked very close to 1. Contrary to the logarithmic scheme, this one is not proven to converge. Lower α yields faster temperature decrease, at the cost of breaking away from the convergence guarantees of low speed temperature decrease.

Others. [START_REF] Ortner | An adaptive simulated annealing cooling schedule for object detection in images[END_REF] propose to use adaptive annealing, that adapts the schedule depending on some statistics of the chain at step t. [START_REF] Guilmeau | Simulated Annealing: A Review and a New Scheme[END_REF] propose to mix Fast Simulated Annealing (that modifies the acceptance probability function) and Sequential Monte Carlo (that simulates a weighted population of particles) into a new scheme.

Stopping conditions

In order to get the exact optimal configuration y * , the RJMCMC needs to run for an infinite amount of steps. In practice, we do not have that much time. For some cases exact simulations [START_REF] Lieshout | Markov Point Processes and Their Applications[END_REF] allow to ensure the simulated Markov chain has reached its theoretical optimum. While this might be necessary for some statistical inference tasks, the application to object detection requires some compromises towards lower inference times. [START_REF] Robert | Diagnosing Convergence[END_REF] propose using statistics on the configurations such as

s t ¯= 1 t t ∑ k=0 s(y t ), (3.72)
where s(y t ) is some measurable characteristic of the configuration at step t; e.g. the number of points n(y t ) or the energy U (y t |X). [START_REF] Lacoste | Extraction de réseaux linéiques à partir d'images satellitaires et aériennes par processus ponctuels marqués (These de Doctorat[END_REF]) compute this empirical mean every N (i.e. average for values t ≡ 0 mod N ) as to limit the autocorrelation effect of the Markov chain. Then, if the empirical mean s t ¯does not change more than a set ϵ s for more than n test sequential computations of s t ¯, the chain is stopped. For more stability one can monitor more than one statistic.

One might also monitor the normalized energy change (U (y ′ ) -U (y))/U (y), the running average of the acceptance ratio (or number of acceptances), and process those as the configurations' statistics above to produce a stopping condition.

For our application we rather set a fixed amount of iterations, although scaling it with the size of the image space S to explore and the amount of parallelization. While this means we sometimes might run more steps than necessary, it alleviates the running of the numerous parameters involved in using multiple statistics to define a good stopping condition. [START_REF] Robert | Diagnosing Convergence[END_REF]) also proposes running parallel chains and computing both interchains and intra-chains statistics to monitor convergence.

CNN fundamentals

Definition 3.4.1. A Convolutional Neural Network (CNN) is a type of artificial neural network that uses convolutions to extract information from data with grid like topology [START_REF] Goodfellow | Deep learning[END_REF]. These models use a combination of convolutions and pooling operations to extract patterns from the data.

For this approach the data need only be in a grid like structure; while 2D image are the most usual application, it can be generalized to 3D for instance with electron microscopy scans [START_REF] Hallou | Deep learning for bioimage analysis in developmental biology[END_REF]. The grid might also not be isotropic, for instance in image sequence segmentation [START_REF] Lalonde | ClusterNet: Detecting Small Objects in Large Scenes by Exploiting Spatio-Temporal Information[END_REF] the third dimension represents time, and in audio classification [START_REF] Hershey | CNN architectures for large-scale audio classification[END_REF] one axis of the spectrogram represents time while the other is frequency. Closer to our subject, multispectral satellite images (see Chapter 1) can be treated as 3D data, with two spatial dimensions on the last corresponding to electromagnetic frequency.

The convolution operator applies a cross-correlation of the input image (or tensor) and a kernel, the output is often referred to as a feature map [START_REF] Goodfellow | Deep learning[END_REF]. In short the crosscorrelation of the kernel on the input, returns a signal that can be interpreted as proportional to the presence of the pattern represented by the kernel at every location of the image. The input and output can be composed of one or several channels (also referred to as features): e.g. the color channels for RGB images. Output and input number of features may be different; an input tensor of size (H, W, C int ) (height, width, number of features/channels) with a kernel of size (K, K, C int , C out ) will produce a tensor of size (H, W, C out ), modulo some padding applied to the input.

The pooling operator reduces the size of the input image or tensor by applying an aggregation function on chunks of the image. For instance for every two-by-two block of pixels, the Max pooling operation keeps only the maximum value of the block to produce the output map. This reduces the spatial size of the input, allowing to aggregate spatial information.

The sequential combination of convolutions and pooling allows alternating between pattern extraction and spatial aggregation, so that steps after steps the extracted patterns grow from low level textures to larger features in the image. The surge in popularity of CNN models stems from the ability to learn the convolution kernels for a specific objective (defined by a loss function to minimize). This is enabled by the ability to back-propagate the gradient of from the loss to the kernel parameters that can them be updated within a gradient descent learning scheme.

One such model is the Unet [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF], it combines convolutions and pooling layers with upsampling (or up-convolutions) that allows to regain spatial resolution, in a sense reversing the pooling layers. The downward path, alternating convolution blocks with pooling, 3.4 -CNN fundamentals 61 progressively extracts bigger size patterns in the image. The downward path, using convolutions and upsampling, rebuilds the spatial resolution of the output progressively, with the help of feature maps transferred from the same depth level of the downward path. This simple architecture allows for fine detail outputs (the output is of same spatial resolution of the input), and higher scale pattern matching (thanks to the depth of the network). This makes it a good fit for pixel wise classification or dense feature extraction. For instance a classical approach in the first image classification networks such as [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF] is to flatten the image after the last convolution layer; a feature map of shape (H ′ , W ′ , F ) gets flattened to a vector of features of size (H ′ W ′ F ), with F the number of features. In turn this vector of features is fed through a sequence of dense layers (in our example C in = H ′ W ′ F ) to finally output a vector of class probabilities through a Softmax. First, this architecture sets the network input to be of fixed size; since H ′ and W ′ are a direct function of the CNN model and the input image size (H, W ). Also, as the processed image gets bigger, the number of elements in the weight W of the dense layer scales quadratically to the image size; for a bigger image of size (sH, sW ), we get weight tensor W of shape (s 2 H ′ W ′ F, C out ).

As satellite images are big and also of varying size, the use of dense layers is prohibitively expensive in terms of number of parameters and computational cost (if applied over all pixels). This is why Fully Convolutional Networks are proposed as such: Definition 3.4.3. A Fully Convolutional Network (FCN) is a type of CNN architecture making use of locally connected layers such as convolution pooling or upsampling [START_REF] Shelhamer | Fully Convolutional Networks for Semantic Segmentation[END_REF].

Avoiding densely connected layers allows for lower parameter counts, and ensures any input size can be used. Lastly the properties of convolutions and pooling blocks make this type of architectures inherently translation invariant. This last property is very useful in the context of biological imagery or remote sensing for instance, where images are large and the processing demands to be translation invariant, as the elements of interest are not necessarily centered as they most often are in usual images (e.g. photographs, taken at eye-level, with a centered subject).

Remark 3.4.1 -The translation invariance is valid if we disregard edge effects. These can be dealt with proper padding or cropping since the size of the edge effect can be computed. To compute a result on big images that do not fit in memory, one has to prepare patches of the source image with padding for each patch. The inference results on each patch can then be cropped to remove edge effects and stitched back together without any difference to inferring on the whole image at once (up to float precision errors). See Appendix B for more details.

By construction these FCN architectures are good at extracting patterns and texture. However, to build longer range interactions between objects within the image requires further increasing the depth at an increased cost in number of parameters and computation. Attention mechanismoriginally introduced for natural language processing [START_REF] Vaswani | Attention is all you need[END_REF] and then applied to object detection [START_REF] Carion | Endto-End Object Detection with Transformers[END_REF] -propose a mechanism to relate each pixel with every other pixel in a feature map. This kind of mechanism is done at a great parameter and computational cost too.

While a lot of work is being done on tuning CNN architecture to any task, in our work we focus on complementing a basic FCN model with the Point Process framework described in the first sections of this chapter; we will not go much further in depth in the design of neural networks but rather consider what we can extract from those.

Conclusion

In this chapter we first presented the building blocks of Point Process models. Those allow modeling configurations of points as a random variable that depend on the underlying image. This model is defined by a Gibbs density, i.e. through an energy function that measures the compatibility of a configuration y with the image X. At inference, for a given image X, the estimated configuration y ˆis obtained by approximating the minimum of that energy function, through a Monte Carlo Markov Chain sampling method combined with simulated annealing.

Meanwhile, Convolutional Neural Networks introduce a broad category of models to transform images into new representations. The key breakthrough of these approaches is to learn the convolution filter parameters on a training dataset through gradient descent.

While Point Processes model the high level objects and their interactions (here our objects of interest: vehicles), the CNN methods shine in performing pixel level transformations. Contrast measures try to bridge this pixel-to-object gap inherent to the Point Process approach, but remain quite tedious to tune and often fail on more complex applications (see Section 4.2.1.1). This mo- In this chapter we build the Point Process through its energy function. First we show the limitations of contrast based measures for satellite images, and proceed to build a data term from the output of a CNN. We propose several ways to build a potential map; we find that a center-pointing-vector to divergence approach gives the best results. We also show how to re-interpret a classical pseudo-probability output of a CNN into an energy for the Point Process. Furthermore, we propose a set of priors to regularize configurations of vehicles, while being careful about the choice of the aggregation operator over the multiple interactions of an object. Lastly, we discuss how the CNN inference is factorized within the energy computation pipeline for a limited computational cost.

PART

-Point

Process as an energy model 69

Point Process as an energy model

In this chapter we build a Point Process for detecting small objects in remotely sensed images. This Point Process is derived from an energy function via a Gibbs density (see (3.9)). For a given image X and configuration y, the unnormalized density of the Point Process is defined as:

h(y) = exp (-U (y, X)) . (4.1)
We now consider a Point Process of rectangles parametrized as follows: A point y in y, has two spatial coordinates in S, y i and y j . The marks (y a , y b , y α ) ∈ M encode respectively the width, length and angle relative to the -→ i vector as shown in Classical formulation of a Gibbs Point Process performs the sum of an energy function over the cliques of configuration y (3.9). But formulating the energy as a sum of energies per pointthat may depend on their neighborhood -is equivalent up to the size of the cliques in question (see Section 3.2.1.1):

U (y, X) = ∑ y∈y V ( y, X, N y {y} ) . (4.2)
For object detection, the energy of a point V (y, X, N y {y} ) should encompass both the fitness of point y against the image X, and the prior knowledge on objects and their interaction. For that matter we write the energy of a point, as a sum of several energy terms:

V ( y, X, N y {y} ) = w 0 + ∑ e∈ξ w e V e ( y, X, N y {y} ) , ( 4.3)
with w 0 , w e some weight parameters, and ξ the set of energy terms, split into the following categories:

• Internal / prior energies encode the internal coherence of the point configuration. Either at a single point level (e ∈ ξ prior-point ) or considering the interactions of a point with its neighborhood (e ∈ ξ prior-interact. ).

• External / data energies (e ∈ ξ data ) reflects how good the point fits the image data X. With this formulation, we measure the correspondence of a single point y against the image X.

The Bayesian approach mentioned in Section 3.2, measuring the likelihood of the image X given the configuration y, would require a term over the V e (y, X).

Formulating the energy of a point V as a sum of energy terms V e -albeit classical, see [START_REF] Ortner | An adaptive simulated annealing cooling schedule for object detection in images[END_REF][START_REF] Lafarge | Geometric Feature Extraction by a Multimarked Point Process[END_REF][START_REF] Descombes | Multiple objects detection in biological images using a marked point process framework[END_REF] -remains a modeling choice; we discuss other possibilities in Section 4.4.2.2.

In the following example we show how the above-described framework allows for easy combinations of simple Point Processes to create more complex patterns.

Example 4.1.1 -Let 3 points processes Φ 1 , Φ 2 , Φ 3 on S × M = [0, 1] 2 × [0.01, 0.1],
each defined to exhibit one basic behavior, and derived from a single energy term V 1 , V 2 , V 3 . We denote y i , y j , y r respectively the position in the two spatial axes, and the radius of each point. The processes are defined as such:

• For Φ 1 , points in the white pixels of image X are more likely to be rejected. For a point y the energy corresponds to the value of image X (values in [0, 1]) at location (y i , y j ):

V 1 (y, X) = X[y i , y j ].
• Process Φ 2 avoids overlaps; with d(y, u) the distance between y and u, energy term V 2 checks if the gap between the two circles is greater than zero with all neighbors, if not it assigns an energy of 1:

V 2 (y, N y {y} ) = max u∈N y {y} {1 d(y,u)-yr-y ′ r <0 (y)}.
• Process Φ 3 favors circles of growing size as points get higher in the image, within a margin 0.01:

V 3 (y) = 1 |yr-(0.99y j +0.01)|>0.01 (y).
Finally, we define a composite process Φ 4 , from the sum of the energies V 1 , V 2 , V 3 (or product of densities). As illustrated in Figure 4.2, Point Process Φ 4 exhibits a combination of all the constraints set on each Φ 1 , Φ 2 , Φ 3 .

Data terms from a CNN

The data terms of the Point Process aim to measure the fitness of each point y ∈ y against the image data X. First we look at contrast measures and their limitations, then we describe our contribution to use Convolution Neural Networks (CNN, introduced in Section 3.4) for these data terms.

Contrast measures

Point Process models in [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF][START_REF] Ortner | An adaptive simulated annealing cooling schedule for object detection in images[END_REF][START_REF] Crȃciun | Stochastic geometry for automatic multiple object detection and tracking in remotely sensed high resolution image sequences[END_REF][START_REF] Descombes | Multiple objects detection in biological images using a marked point process framework[END_REF] use contrast measures to build a data potential for each point. The idea is to compare statistical properties of pixels values of image X inside and around the object y.

The idea is to compare the distributions inside and outside the shape formed by y. For example [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF] use the T-test on the values inside shape i(y) and on its contour c(y):

V (y, X) ∝ |µ i(y) -µ c(y) | √ σ 2 i(y) n i(y) + σ 2 c(y) n c(y) , ( 4.4)
with µ i(y) , σ i(y) the mean and variance of values of image X in interior i(y), similarly for contour s(y).

Alternatively [START_REF] Kulikova | A Marked Point Process Model Including Strong Prior Shape Information Applied to Multiple Object Extraction From Images[END_REF] propose to check the image gradient ∇X against the normal vector -→ n y (t) along the contour line s y (t), t ∈ [0, 1]:

V (y, X) ∝ ∫ t∈[0,1] -→ n y (t) • ∇X[s y (t)] √ |∇X[s y (t)]| 2 + ϵ dt, ( 4.5)
where ∇X[s y (t)] is the image gradient at contour coordinate t.

While there are many other formulations for contrast measure, most rely on the assumption of simple backgrounds, limited object of interest variations, and having little to no high contrast background elements. Moreover, the definition of interior and outside contour of the object y is to be fine-tuned to the application; e.g. the gap and thickness of considered contours (see Figure 4.3). Similarly, in order to deal with heterogeneity inside the object [START_REF] Crȃciun | Stochastic geometry for automatic multiple object detection and tracking in remotely sensed high resolution image sequences[END_REF] measure the contrast between an interior contour and exterior one (red and blue in Figure 4.3); as the visual features inside the boats (e.g. windows, decks, etc...) increase the variance of inside pixels.

On the other hand CNN models have shown to be very efficient at extracting features from images for object detection and classification. In the following section, we will show how to interpret a CNN based object detector output to obtain an energy that measures the fitness of a configuration against an image. It allows us to go past the contrast measure design by utilizing a pre-trained CNN output. The next example shows a basic comparison between contrast measures and a CNN model.

Limitations of contrast measures: example

We compare the performances of contrast measures on both remote sensing images and synthetic data. We evaluate the performance of each by generating a set of object proposals from the image, 99% of which do not correspond to an object, 1% are actual objects of interest. For a contrast measure to be effective, it should be able to score high on objects of interest and low on negative proposals. For every contrast threshold we can assign positive and negative detections to each proposal, thus compute precision and recall values. Computing precision and recall for every threshold gives us the precision-recall curve.

We compare the two contrast measures given in (4.4) and (4.5) against our CNN model (see next sections) and a measure that always returns a constant value 1 .

This comparison is performed on both a synthetic image and a sample from the dataset DOTA [START_REF] Xia | DOTA: A Large-Scale Dataset for Object Detection in Aerial Images[END_REF]) (sampled at 0.5m resolution). The synthetic data consists of bright or dark rectangles on a gray background with some additive noise. The data and results are shows in Figure 4.4.

The contrast measure performs excellently on the synthetic data, but struggles on real images. Indeed, this image contains a lot of high contrast elements that are not objects of interest; airconditioning units, pavement features, trees etc... To build a data term from a CNN for the Point Process, we aim at inferring from a given image X (Figure 4.5 (a)) a map of potentials V pos (Figure 4.5 (b)), so that the data energy V pos of a point y can be computed by simply looking up the value of V pos at the corresponding position. Previous works on marked Point Processes using contrast measures have shown good results when applied to high contrast images where background clutter is limited. Thus, our first approach consists in applying a tailored contrast measure in an image transformed by a FCN model into a more simple representation. Ideally the image transformation model should produce simple contrasted patterns on object positions, while simplifying background to be as close as possible to uniformity. Applying the pattern matching would in turn provide potential map V pos . The following approach was published in (Mabon et al., 2021).

Inferring heatmaps. First we must determine the output of our CNN. We need to build target (or Ground Truth, GT) maps with desired properties on training data, for the model to be able to reproduce such maps on unseen data. Object detection models such as (T. [START_REF] Zhao | An Arbitrary-Oriented Object Detector Based on Variant Gaussian Label in Remote Sensing Images[END_REF][START_REF] Yi | Oriented Object Detection in Aerial Images with Box Boundary-Aware Vectors[END_REF][START_REF] Huang | A General Gaussian Heatmap Label Assignment for Arbitrary-Oriented Object Detection[END_REF] use heatmaps to identify object centers. Definition 4.2.1. For object detection, a heatmap refers to a map (often 2 dimensional) of probability-like values measuring the presence of an element or event; in our case object centers. [START_REF] Guo | On Calibration of Modern Neural Networks[END_REF] shows, these measures are often poorly calibrated and should rather be considered as scores. We will also refer to those as probability estimates.

Training the CNN model to infer an exact map such as it returns 1 at object centers and zeros elsewhere (see Figure 4.6) proves to be impractical; the image data and annotations are sometimes off by a few pixels (making the sharp target map incorrect) and the labels are too imbalanced to properly train the model (i.e. the ratio of center pixels to background is close to zero).

Relaxing the heatmap. We propose relaxing the training target, in learning to infer a map H such that for every pixel ρ in the image:

H[ρ] = max { 1 - d(ρ, y GT ) r max , 0 } , ( 4.6)
where d(ρ, y GT ) is the distance of pixel ρ to the nearest point in y GT and r max a threshold term (see Figure 4.6). This relaxing produces a relaxed heatmap target to train out model on.

Balancing. The works in (M. Bai & Urtasun, 2017) on inference of discrete watershed map for instance segmentation propose a solution to alleviate class imbalance issues. Their approach requires them to infer a distance map in the likes of (4.6); they propose to discretize the distance values into n c distinct bins or classes (see Figure 4.6). We now have to use the CNN to estimate the most likely class (or value bin). This allows to re-weight the loss over each class to deal with the imbalance.

Definition 4.2.2. In image processing, a watershed [START_REF] Beucher | Use of watersheds in contour detection[END_REF][START_REF] Vincent | Watersheds in digital spaces: An efficient algorithm based on immersion simulations[END_REF]) is a transform defined on a grayscale image, where lower levels (lower pixel brightness) of the image get progressively flooded to form a segmentation of the image. By extension the watershed map refers to the image the watershed method is applied on.

The ground truth (or training target) map is then, for every pixel ρ: for ease of reading.

H dscr [ρ, c] = ⎧ ⎨ ⎩ 1 if H[ρ] ∈ [ c-1 nc , c nc ] , 0 else
Estimation and training. The Unet estimates scores ˆ H dscr from image X, for each discrete class of values c, by the application of a Softmax (Goodfellow et al., 2016, 6.2.2.3 Softmax Units for Multinoulli Output Distributions) on the Unet output Z of size (H, W, n c ) as :

Softmax(Z[ρ]) c = exp (Z[ρ, c]) ∑ nc c ′ =1 exp (Z[ρ, c ′ ]) . (4.8)
The CNN model is trained to minimize the cross entropy loss [START_REF] Kline | Revisiting squared-error and cross-entropy functions for training neural network classifiers[END_REF]:

L(H dscr , ˆ H dscr ) = - ∑ ρ nc ∑ c=1 w c H dscr [ρ, c] log ( ˆ H dscr [ρ, c] ) , ( 4.9)
with c] a weighting factor introduced to amplify the loss on the less frequent classes (i.e. foreground classes in our case).

w c = 1 -1 nc ∑ ρ H dscr [ρ,
Building the data potential. Since we know by construction of the target map the aspect of the map around object centers, we devise a filter K as illustrated in Figure 4.7 (a) to look for object centers when cross-correlated with the output map ˆ H dscr . When the filter pattern matches the pattern in the output map we get a high response, otherwise it is low, as shown in Figure 4.7 (b). This approach is akin to using a contrast measure as the filter is looking at a specifically shaped pattern of contrast in the image. Using cross-correlation to look for this pattern allows to pre-compute the contrast for every pixel of the image easily. We denote

F[ρ] = (K * ˆ H dscr )[ρ]
the cross-correlation response at pixel ρ, and thus define the data potential2 as :

V pos (y, X) = q pos (F[ρ y ]) , (4.10)

q pos (x) = ⎧ ⎨ ⎩ -2(x-d 0 ) d 1 -d 0 + 1 if x < d 1 exp ( -2(x-d 0 ) d 1 -d 0 + 2 ) -2 else, (4.11)
where ρ y is the pixel containing y. The quality function and its two parameters d 0 , d 1 is illustrated in Figure 4.7. With the data potential defined as such, the map F can be pre-computed for a given image X from the Unet output, to be then used for any computation of U (y, X), simplifying the contrast computation to the sampling of an array. Limitations. While this first approach allows to extract an energy map, we notice that nearby blobs formed by objects in the energy map can sometimes clump together. We illustrate this in Figure 4.5; the clumping of nearby detections into one long blob is quite pronounced for cars parked side-to-side (as opposed to front-to-back). Moreover, the energy pits around detected centers are quite wide -as a side effect of the relaxation -making the positioning of object imprecise. Remark 4.2.2 -The focal loss [START_REF] Lin | Focal loss for dense object detection[END_REF] propose to weight the loss for imbalanced data such as heatmap without any discretization. In our case, learning to infer the heatmap H directly with focal loss does not improve significantly on the clumping issue.

d 0 d 1 x 2 0 q(x) (c) q pos (x)

Inferring vector fields

The previous approach to build a data energy is limited in how it can efficiently separate nearby instances. Here we propose an approach based on vectors and their divergence to better pinpoint the location of objects. The following approach has been published in (Mabon et al., 2022b) and (Mabon, Ortner, & Zerubia, 2022c). Vectors point to object center. To do instance separation [START_REF] Neven | Instance Segmentation by Jointly Optimizing Spatial Embeddings and Clustering Bandwidth[END_REF] propose to infer for each pixel of the image the vector from the pixel location to the center of the object the pixel belongs to. Outside the objects the vectors are not defined.

With our small objects, we rather choose to extend the range at which these vectors are inferred outside the area of the object; within a radius r max the model has to infer the unit vector that points towards the nearest object center as illustrated in Figure 4.8. In short, we look to infer the map:

H ∇ [ρ] = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1 ∥y-ρ∥ 2 [ y i -ρ i y j -ρ j ] , y = arg min y ′ ∈y GT {∥y ′ -ρ∥ 2 } if H M [ρ] = 1 [ 0 0 ] otherwise (4.12)
The map H M is a binary map, equal to one if less than r max away from an object center, else it is zero:

H M [ρ] = ⎧ ⎪ ⎨ ⎪ ⎩ 1 if arg min y ′ ∈y GT {∥y ′ -ρ∥ 2 } ≤ r max 0 otherwise. (4.13)
Remark 4.2.3 -Note that H ∇ can be obtained from the image gradient of H defined in (4.6):

H ∇ [ρ] = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∇H[ρ] ∥∇H[ρ]∥ 2 if H M [ρ] = 1 [ 0 0 ] otherwise , ∇H[ρ] = [ ∂H ∂i [ρ] ∂H ∂j [ρ] ] . ( 4

.14)

The gradient is defined over continuous space, but can be easily approximated in discrete pixel space by the following cross correlation:

∇H ≃ ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ [ -1 0 +1 ] * H ⎡ ⎢ ⎣ -1 0 +1 ⎤ ⎥ ⎦ * H ⎤ ⎥ ⎥ ⎥ ⎥ ⎦
Inference and training. To produce a field of vectors ˆ H ∇ from image X, we modify the Unet to output two tensors Z ∇ and Z M of respective sizes (H, W, 2) and (H, W, 1). The output Z M allows estimating H M as such:

ˆ H M [ρ] = σ (Z M [ρ]) . (4.15)
Thus, the estimate of H ∇ can be obtained as:

 H ∇ [ρ] = ˆ H M [ρ] Z ∇ [ρ] ∥Z ∇ [ρ]∥ 2 , (4.16) with ∥Z ∇ [ρ]∥ 2 the norm of 2-dimensional vector Z ∇ [ρ].
The normalizing allows to consistently maintain unit vectors in the neighborhood of objects, while the product with ˆ H M [ρ] sets to zero all out-of-range pixels. The model is trained to minimize the Mean Square Error (MSE) between the estimated and ground truth vectors (inside the r max -neighborhood of each point) while minimizing the Binary Cross Entropy (BCE) between the mask and its estimation:

L( ˆ H ∇ , ˆ H M , H ∇ , H M ) = MSE( ˆ H ∇ ,H ∇ ) where H M [ρ]=1 ⏟ ⏞⏞ ⏟ ∑ ρ H M [ρ] ⃦ ⃦ ⃦ ˆ H ∇ [ρ] -H ∇ [ρ] ⃦ ⃦ ⃦ 2 2 - ∑ ρ H M [ρ] ln( ˆ H M [ρ]) ⏞ ⏟⏟ ⏞ BCE( ˆ H M ,H M )
.

(4.17)

We provide some details about the training procedures in Appendix B.

Remark 4.2.4 -In practice, we infer a single tensor Z of size (H, W, 3) and then extract the relevant features/channels to obtain Z ∇ and Z M of respective sizes (H, W, 2) and (H, W, 1).

Divergence as data potential. Given the inferred ˆ H ∇ for image X, we need to build a scalar data potential for each element of y. By construction the center of each object corresponds to a point where vectors converge. It can be seen clearly in Figure 4.8 (c) and (d). Thus, we use the divergence operator to extract object locations; the object centers get negative values, while separation lines in between get positive values (see Figure 4.8 (d)). Denoting ˆ H div = div ˆ H ∇ :

V pos (y, X) = ˆ H div [ρ y ]. (4.18)
We show the potential map computed over a real unseen (i.e. not trained on) image in Figure 4.9.

Computing the divergence. The divergence in continuous space over a vector field F in 2D is defined as:

div F = ∂F i ∂i + ∂F j ∂j . ( 4

.19)

The divergence field over in raster grid coordinates can be approximated using the finite difference approximation for derivatives:

∂F i ∂i [y i ] ≃ F i [y i + 1] -F i [y i -1] 2 . ( 4

.20)

The approximate derivate in an image can be obtained by convolution using a Prewitt [START_REF] Prewitt | Object enhancement and extraction[END_REF] filter:

∂F i ∂i ≃ ⎡ ⎢ ⎣ +1 +1 +1 0 0 0 -1 -1 -1 ⎤ ⎥ ⎦ * F i , ∂F j ∂j ≃ ⎡ ⎢ ⎣ -1 0 +1 -1 0 +1 -1 0 +1 ⎤ ⎥ ⎦ * F j . (4.21)
The computation of the divergence for a vector field defined in a raster grid is detailed in Appendix I.3.

Interpolation. The energy described in (4.18) only samples the value of the ˆ H div at integer locations (exact pixel locations). However, points y in configuration y have real-valued locations (y i , y j ) ∈ R 2 . Thus, we propose interpolating values in between exact pixel locations, defining:

V pos (y, X) = ˆ H div [y i , y j ], (4.22)
with ˆ H div [y i , y j ] the bilinear interpolation of ˆ H div at location (y i , y j ). In this section we contributed in proposing an energy map built from the divergence of a field of vector inferred by a CNN, allowing for better separability of small tightly-packed instances.

Interpreting CNN outputs as energy

Until now, we proposed several network architectures to produce a map of potentials for the data term of our Point Process. In this section we propose a general framework to use the output of any CNN model (within some specifications) to produce the required potential. This was published in (Mabon et al., 2023a) and (Mabon, Ortner, & Zerubia, 2023b) and to be submitted in (Mabon, Ortner, & Zerubia, 2023). The core idea stems from [START_REF] Grathwohl | Your classifier is secretly an energy based model and you should treat it like one[END_REF] where they propose to interpret classifier model as Energy Based Models (EBM).

Here we consider a class of CNN model that infers a map of scores of centers for every pixel of the image such as the heatmap based models [START_REF] Huang | A General Gaussian Heatmap Label Assignment for Arbitrary-Oriented Object Detection[END_REF]. The model described in the previous section can also be adapted to produce probability estimates (scores) as explained in Remark 4.2.5. This class of models produce a probability estimates of center for pixel ρ from a tensor ˆ Z pos inferred by the CNN. The probability estimate p ˆ(ρ) is obtained through the sigmoid function as: This yields:

p ˆ(ρ) = σ ( ˆ Z pos [ρ] ) = 1 1 + exp ( -ˆ Z pos [ρ] ) . ( 4
V pos (y, X) = q pos ( ˆ Z pos [y i , y j ]
) (4.26)

q pos (x) = ln (1 + exp (-x + t pos )) , ∀x ∈ R, (4.27)
where t pos is a scalar parameter we introduce to offset the inflection point of function x ↦ → ln(1 + exp(-x)) as shown in Figure 4.10. Once again we can perform bilinear interpolation to have the potential be defined for any real-valued coordinate in S instead of discrete pixels.

Data energy on marks

For now, we only considered potentials on the position of the objects, in this section we propose a potential on the marks of the Point Process. First we discretize each mark κ (κ ∈ {a, b, α} in the case of a Point Process (PP) of rectangles), into n κ classes (or value bins) in the range [κ min , κ max ]. We define the integer class of a value v for mark κ as: with ⌈c κ (v)⌉ the corresponding integer class. We illustrate the discretization in Figure 4.11: here we pick n κ = 8 for κ ∈ {a, b, α}, a low value for illustration purposes (in practice we pick n κ = 32). For instance, in this illustration we have ⌈c α (y α )⌉ = c 1 α , ⌈c a (y a )⌉ = c 3 a etc... Now, supposing we have a model trained to classify the mark of an object at a given position (i, j) ∈ S, such model will produce a probability estimate of mark κ to be in class c ∈ 1, n κ , from model output ˆ Z i,j κ ∈ R nκ as :

c κ (v) = n κ v -κ min κ max -κ min , ∀v ∈ [κ min , κ max ], (4
Z κ i = y i j = y j H W n κ c 1 κ c 2 κ c 3 κ c 4 κ c 5 κ c 6 κ c 7 κ c 8 κ Z i,j κ Z i,j κ [ c κ (y κ ) ] - + ln c exp V κ (y, X) c κ (y κ ) y j y i Z κ A → -A + ln c exp(A) V κ (y, X)
p ˆ(c|i, j, X) = Softmax( ˆ Z i,j κ ) c = exp ( ˆ Z i,j κ [c] ) ∑ nκ c ′ =1 exp ( ˆ Z i,j κ [c ′ ] ) . (4.29)
As in [START_REF] Grathwohl | Your classifier is secretly an energy based model and you should treat it like one[END_REF], we can reinterpret the model output into a potential, giving:

V κ (y, X) = -ˆ Z i,j κ [⌈c κ (y κ )⌉] + ln ( nκ ∑ c=1 exp ( ˆ Z i,j κ [c] ) ) , (4.30)
we justify this formulation in Appendix C. As the CNN outputs a tensor ˆ Z κ of shape (H, W, n κ ), we get vector ˆ Z i,j κ by sampling tensor ˆ Z κ at location (⌈y i ⌉, ⌈y j ⌉), as illustrated in Figure 4.12 (a).

Interpolation. Same as the first proposed position energies, the energy on marks defined in 4.30 is defined over the discrete values of the mark. As the marks live in a continuous space, we 4.3 -Priors on configurations 83 propose using interpolation to defined V κ in between discrete mark values:

V κ (y, X) = ˆ Z ⋆ κ [y i , y j , c κ (y κ )] (4.31)  Z ⋆ κ [i, j] = -ˆ Z κ [i, j] + ln nκ ∑ c=1 exp ( ˆ Z κ [i, j, c] ) , ∀i, j ∈ S d . (4.32)
The above, illustrated in Figure 4.12 (b), has two benefits: first (4.31) defines the energy for continuous values of the marks using bilinear interpolation, thus ensuring Lipschitz continuity of the energy. Second we can pre-compute (4.32) once for a given image; making the mark energy computation a simple value lookup and interpolation for each point. This pre-computation is quite fast, as it is defined as an operation over a tensor, which is implicitly parallelized by tensor computation libraries such as (PyTorch, n.d.).

Remark 4.2.6 - [START_REF] Grathwohl | Your classifier is secretly an energy based model and you should treat it like one[END_REF] propose to use the extra degree of freedom of the model output to measure the probability of the observation itself (not its class); indeed the Softmax function is invariant to any offset on (Softmax(A

+ b) = Softmax(A), b ∈ R)
which provides a degree of freedom in A.

In our case we can use ln

( ∑ nκ c=1 exp ( ˆ Z i,j κ [c] ))
as a position potential; measuring the probability of an object at that location regardless of the class. While this would allow to join the position and mark potentials, our experiments have not been conclusive for the moment.

Our contribution in the last two subsections, is to propose to use any classical CNN for detection as a data term for the Point Process (for both position and marks), based on the suggestion of [START_REF] Grathwohl | Your classifier is secretly an energy based model and you should treat it like one[END_REF] to interpret classifiers as energy models.

Priors on configurations

With the data potentials defined, we focus now on the design of priors to regularize the configurations of objects. We first detail the potentials that apply to single points, then those that measure the interactions. Some priors introduced in this section are illustrated in Figure 4.13.

Point priors 4.3.1.1 Area and ratio priors

As our objects of interest have quite consistent sizes, we favor configurations where shape parameters fit within some expected range. The ratio prior, compares the point y width to length ratio against a parameter µ ratio

V ratio (y) = -exp ⎛ ⎜ ⎝- ( ya y b -µ ratio ) 2

Joint area and ratio prior

The previous two energy terms V ratio and V area have a limitation: Let us consider objects that fall into two groups with typical ratio and area parameters (µ area,1 , µ ratio,1 ) for the first group and (µ area,2 , µ ratio,2 ) for the second. With previously defined area and ratio potentials, we could define two of each. However, this would create four attraction points, favoring points of geometries close to (µ area,1 , µ ratio,1 ) and (µ area,2 , µ ratio,2 ) as needed, but also shapes close to (µ area,1 , µ ratio, 2) and (µ area,2 , µ ratio,1 ) that may not exist in reality. This is illustrated in Figure 4.14 (d).

To deal with this issue we propose a joint area and prior potential, that combine the two to favor points around (µ area , µ ratio ) in the area-ratio subspace as illustrated in Figure 4.14 (e). This joint area and ratio prior is defined as:

V jntAR (y) = -exp ⎛ ⎜ ⎝- ( ya y b -µ ratio ) 2 2σ 2 ratio - (y a y b -µ area ) 2 2σ 2 area ⎞ ⎟ ⎠ . (4.34)
Measured joint distribution. Plotting the effective distribution of area and prior over the data, two modes can be identified for each (see Figure 4.14 (a) and (b)). However, using the V area and V ratio priors would create four attractive points. As can be seen in Figure 4.14 (c), the data exhibits only two modes (cars and trucks); justifying the use of two joint priors V jntAR . One approach is to set parameters (µ area , µ ratio ) manually from looking at the joint distribution in Figure 4.14 (c); however we propose in Section 5.2.2, to automatically estimate these parameters from the data, in a unified estimation procedure. We show the results of a model using these automatically estimated parameters in Section 6.2.2.

Interaction priors 4.3.2.1 Aggregation operators

Some energy models for Point Processes sum all the interactions over the whole configuration as seen in (3.23). However, we find this creates some problems with attractive potentials: Objects sometimes accumulate the negative potentials from the interactions, as the number of interactions grows quadratically relative to the number of points. Meanwhile, isolated points remain at high energy.

For interaction priors (V e , e ∈ ξ prior-interact. ), we propose using a carefully chosen aggregation operator F e along a pair-interaction potential v as such: To avoid the quadratic increase of interaction potentials, we choose F e such that:

∀A > 0, F e : P([-A, A]) ↦ → [-A, A], ( 4.36)
where P([-A, A]) is the power set (set of all subsets) of [-A, A]; i.e. if v e is bounded by A, V e is too. For instance if F e is a maximum, minimum, or average the above condition is met (as in Example 4.1.1). If it is a sum, then the model is equivalent to the above-mentioned (3.23) (up to a factor 2 on the interaction potentials, as each interaction would be counted twice).

In practice this helps to build some versatile priors; e.g. "the object needs to be aligned with at least one object", is easily enforceable with a min over all alignment potentials. Meanwhile, it helps avoid accumulation of points in attractive neighborhoods. We further discuss the choice of aggregation function F e in Section G.

Non overlapping prior

The non-Bayesian approach mentioned in Section 3.2, may lead to the accumulation of objects in areas with low data potential without a strong prior on object not overlapping. For the objects we study -in most situations -no overlap is permitted. We propose a simple prior to penalize overlaps:

V ovrlp (y, N y {y} ) = max u∈N y {y} {v ovrlp (y, u)} v ovrlp (y, u) = max { 0, Area(u ∩ y) min{Area(u), Area(y)} -t ovrlp } , (4.37)
where t ovrlp is a positive threshold that may allow some permeability between objects. We use a maximum to aggregate the interaction potentials v ovrlp as we aim at penalizing overlap with any neighboring objects. We find there is no simple closed-form formula for the area of intersection of two oriented rectangles. Thus, in practice, we approximate the intersection measure as described in Appendix D.

-4.3.2 Interaction priors 87

Alignment prior

The objects of interest (small vehicles), tend to be aligned: either parked with the same heading, or circulating in parallel lanes. We introduce a prior that favors aligned configurations where t align ∈ {0, π/2} is a parameter picked to either favor parallel objects (if 0) or orthogonal objects (if π/2). The aggregation operation over all the interaction potentials v align is picked to be a minimum to favor being aligned with at least one object.

Repulsive and attractive priors

Vehicles often maintain a certain distance, especially when parked close. We introduce two concurring priors based on the gap between two rectangles (see (D.16) in Appendix D), one repulsive, the other attractive; with carefully picked thresholds t repls , t attrc we can favor objects within a specific distance range (see Figure 4.15): with Gap(y, u), measuring the smallest distance between the two shapes, as defined in (D.16).

V repls ( y|N y {y} ) = max u∈N y {y} {v repls (u, v)} v repls (u, v) = max { 0, 1 - Gap(y, u) -t repls d max -t repls } , ( 4
Gap(y, u) vrepls (y, u)

0 t repls Gap(y, u) vattrc(y, u) 0 tattrc
Gap(y, u)

0 t repls tattrc + =

Zero neighborhood prior

For all the interaction priors defined above, the negative of positive potentials added to interacting points might offset the overall energy values of interacting points over points that have no interactions. For that matter we introduce a simple prior aimed at compensating this effect; the potential is 1 when no neighbor is present, else 0:

V zrNbr (y, N y {y} ) = 1 |N y {y} |=0 (y). (4.41)
Here we pick a value of 1 as it we will be able to scale it as needed using the energy weight w zrNbr (see (4.3)). This extra term acts similarly as the offset term w 0 , but applies only to objects without neighbors. In (Mabon et al., 2021) we propose a prior on triplets of points. In this work we consider a non-marked Point Process (i.e. no orientation for the detected objects); there is no possibility to define alignment with only two objects. Thus, we look at the angle formed by triplets of vehicles within the neighborhood N y {y} of each object y in configuration y. We observe that vehicles are most often stopped or circulating in lines; that is why we choose to enforce an angle of 0 mod π. This prior is no longer needed once dealing with marked objects, as the simple alignment prior already favors this type of configurations.

Triplet priors

y u v δ α = π 2 y u v δ α = 0 y u v δ α 0 (a) Example triplet alignments 0 max /2 1 0 q( ) (b) Plot of q(δ α )
Denoting ∡uyv the value (in radiant) of the angle formed by the triplet of points u, y, v (as shown in Figure 4.16 (a)), we define δ α (y, u, v) the absolute angle, and the triplet alignment potential v triAl as follows: (y, u, v)) . (4.43) To turn the absolute difference δ α into a potential, we use the following quality function, parametrized by δ α max , which sets a margin for the angle as illustrated in Figure 4.16 (b):

δ α (y, u, v) = min{|∡uyv|, |π -∡uyv|}, (4.42) v triAl (y, u, v) = q triAl (δ α
q triAl (δ) = - 1 δ 2 α max ( 1 + δ 2 α max 1 + min{δ, δ α max } 2 -1 ) (4.44)
Finally, the potential for point y and its neighborhood N y {y} takes the minimum over all potentials; this favors points that are aligned within at least one triplet: With all potentials defined above, we can compute the energy for any configuration y. We compute the energy as follows:

V triAl (y, N y {y} ) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ min u,v∈N y {y} u̸ =v {v triAl (y, u, v)} if |N y {y} | ≥ 2 0 else. ( 4
1. For a given image X, pre-compute the potential maps once: Within a MCMC framework, in which the energy has to be computed for different configurations with the same image, only steps 2 (a) to (c) have to be repeated as the image and inferred tensors are unchanged (no need to redo step 1 each time). We illustrate this pipeline in Figure 4.17: the pre-computed block corresponds to the results from step 1 that needs to be computed only once per image. The blocks κ ∈ {a, b, α}, e ∈ ξ prior-interact. , e ∈ ξ prior-point are repeated for each respective mark, point prior and interaction prior. For the interaction prior the distance matrix is only computed once for all interaction priors. For more precision on the actual implementation of this pipeline we refer the reader to Appendix I.

Discussions

Continuity

In Section 3.3.3, we mention that the energy function has to be Lipschitz-continuous so that diffusion can be used to sample configurations. The data energy terms are made continuous through interpolation, while the local prior terms are continuous by construction. The interaction energies v e are continuous too. However, object y moving continuously into (or out of) the neighborhood of another point u creates a discontinuity as it crosses the boundary of the neighborhood N y {u} . It is made obvious when simplifying to only two points, we get:

V e (y, {u}) = 1 d(u,y)≤dmax (y)v e (y, u).
y 1 i , y 1 j y 2 i , y 2 j y 3 i , y 3 j Z pos [y i , y j ]
(n(y))

V pos (•)

Z κ (H, W, nκ) y 1 i , y 1 j , y 1 κ y 2 i , y 2 j , y 2 κ y 3 i , y 3 j , y 3 κ Z κ [y i , y j , c κ (y κ )]
(n(y))

V κ (•)

•

y 1 • y 2 • y 3 y (n(y), d) V e (•)
•

y 1 • y 2 • y 3 y (n(y), d) V e (•) d(y, y ) (n(y), n(y)) ⊕ ×w y∈y U (y, X) y = {y 1 , y 2 , y 3 } f (•) Element-wise operation ⊕ Concatenation yi, yj, yκ Value of at (y i , y j , y κ )

Mixing of energy terms

The classical approach for mixing energy terms via a weighted sum (see (4.3)) has its limitations as illustrated in Section 4. 3.1.2. In (Mabon et al., 2022a) we propose using a hierarchical model, where the priors are only taken into account if the position potential reaches a certain threshold:

V ( y, X, N y {y} ) = w pos V pos (y) + 1 Vpos(y)<tpos ⎛ ⎝ ∑ e∈ξ, e̸ =pos w e V e ( y, X, N y {y} ) ⎞ ⎠ .
The latter induces some more discontinuities in the energy model that do not go well with the jump-diffusion sampling [START_REF] Geman | Diffusions for Global Optimization[END_REF]) (see Section 3.3.3). Alternatively, we can see the function V can be any function that transforms the vector of all energy terms into a scalar energy:

V : R |ξ| → R v y ↦ → V (v y ), v y = [ V e 1 (y) . . . V e |ξ| (y) ] ⊤ (4.48)

CHAPTER 4 -Building an Energy Based Model for object detection

For instance one could use a Multi Layer Perceptron (MLP) parametrized by θ, such that V (v y ) = MLP θ (v y ). As the MLP is a universal approximator (given enough hidden layers) [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF] it could better approach the optimal mixing function. The criterion for optimality will be discussed in Section 5.2.

Towards generic (learnable) energy terms

In this chapter we proposed a series of priors with parameters to either set manually of estimate from the data. The selection and design of the energy prior functions to put in the total energy U is by itself a prior on the model, as it restricts the space of energies U (. . . , θ) generated by parameters θ. As prospective work, we propose in Appendix G to push further the parametrization of the energy priors in order to bypass the function-design-induced constraints, while allowing to build energy functions that better fit the data. We first go through a design of per-point priors with Multi Layer Perceptron (MLP), and then show that attention mechanisms can be used to learn a generic aggregation function for interaction priors.

Conclusion

In this Chapter we built the Point Process through its energy function, using a CNN output for the data term and introducing multiple priors. We summarize our contributions as follows:

• A first data model for non-marked Point Processes, using a heatmap inferred by a CNN (Section 4.2.2.1).

• We further refine this approach by inferring a map of vectors and computing the divergence; it achieves better instance separation (Section 4.2.2.2).

• Then we propose a way to interpret most CNN outputs as energies that we can use on our model, both for object positions and marks (Section 4.2.2.3).

• All the above data terms based on a CNN output benefit from replacing the contrast measures computation by a simple value lookup and interpolation (Section 4.2.2.2 and 4.2.3), thanks to the pre-computation of energy maps (Section 4.4.1).

• Joint priors on object parameters to avoid favoring non-existing shapes when dealing with multiple shape modes (e.g. cars and trucks) (Section 4.3.1.2).

• Carefully chosen aggregation functions on interaction priors to avoid energy explosion and favor specific configurations (e.g. alignment with at least one other object, no overlap with any other object) (Section 4.3.2.1).

• As prospective work, we propose a way to solve discontinuities in object interactions (Section 4.4.2.1) and new models to mix the energy terms (Section 4.4.2.2).

CHAPTER 5

Point Processes as Energy Based Models

Dans ce chapitre, nous passons en revue les procédures d'échantillonnage et d'estimation des paramètres pour le Processus Ponctuel défini au chapitre 4. Nous proposons d'exploiter les cartes de potentiels fournies par le CNN, en élaborant un noyau de perturbation local et un noyau de naissance, tous deux construits à partir des données. De plus, nous proposons un schéma de calcul parallèle modifié pour échantillonner de nouveaux points et perturbations qui se concentre sur les zones de l'image à forte densité d'objets. Pour l'estimation des paramètres du Processus Ponctuel, nous proposons d'abord une approche fondée sur les séparateurs à vaste marge pour résoudre les limitations des méthodes précédentes. La deuxième approche et contribution majeure consiste à appliquer des méthodes d'apprentissage de modèles basés sur l'énergie pour estimer tous les paramètres du Processus Ponctuel au sein d'un algorithme unifié par divergence contrastive. Enfin, nous proposons un score d'objet pour le calcul de métriques utilisant l'intensité de Papangelou.

In this chapter we go through the sampling and parameter estimation procedures for the Point Process defined in Chapter 4. We propose to leverage the potential maps provided by the CNN, by building a data driven local perturbation kernel and birth kernel. Moreover, we propose a modified parallel computation scheme to sample new points and perturbation that focuses on areas of the image with high object density. For the estimation of the Point Process parameters, we first propose an approach based on Support Vector Machines to alleviate limitations of previous methods. The second approach and major contribution of this thesis, consists in applying Energy Based Model learning methods to estimate all parameters of the Point Process within one unified algorithm through contrastive divergence. Finally, we propose an object score for metrics computation based on the Papangelou intensity.

Sampling

In Chapter 4 we built our Point Process model through the definition of its energy function U . This section presents the improvements over the sampling procedures introduced in Section 3.3. We leverage the pre-computed energy potential maps provided by the CNN model to build more efficient sampling. First by using local perturbations based on the potential map adapted from [START_REF] Ortner | Processus ponctuels marqués pour l'extraction automatique de caricatures de bâtiments à partir de modèles numériques d'élévation (These de doctorat[END_REF][START_REF] Descombes | Stochastic geometry for image analysis[END_REF], which we published in (Mabon et al., 2021). We adapt the birth map from [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF] to our model, and propose using a truncated energy model to sample new points (Mabon et al., 2023a). We then adapt parallel Point Process sampling from [START_REF] Verdié | Efficient Monte Carlo sampler for detecting parametric objects in large scenes[END_REF] with cell picking based on the potential maps, as published in (Mabon et al., 2023a) and to be submitted in (Mabon et al., 2023). Finally, the resulting sampling procedure is presented in Algorithm 5.1, based on Jump Diffusion while leveraging automatic gradient computation.

Data driven kernels

Local perturbation based on data

The local perturbation kernel introduced in Section 3.3.2.3 usually uses a zero centered Gaussian distribution to generate the perturbation δ on a point y on configuration y. When computing the total energy we first infer a position potential map V pos so that we simply have to sample values from it to get V pos . It is defined as: (5.1) for any pixel ρ ∈ 0, H × 0, W . Similarly, we can define V κ for each mark κ ∈ {a, b, α}.

V pos [ρ] = V pos (ρ, X),
Given that we have access to the pre-computed position potential map, we could intuitively use it to make better move propositions within the Markov chain.

In order to propose more relevant perturbations, we propose to use this potential map. From (3.64) we have:

α(y, x) = min { 1, h(x)g ′ (δ ′ ) h(y)g(δ) } ,
Where g(δ) is the density for picking perturbation vector δ and g ′ (δ ′ ) the density of the reverse move.

In (Mabon et al., 2021) we propose a data-based translation kernel adapted from [START_REF] Ortner | Processus ponctuels marqués pour l'extraction automatique de caricatures de bâtiments à partir de modèles numériques d'élévation (These de doctorat[END_REF][START_REF] Descombes | Stochastic geometry for image analysis[END_REF] using the position potential map V pos . This proposal is applied to a non-marked Point Process (thus only performing translations) but can be generalized to a marked Point Process too.

For a point y in configuration y, y ∈ S, we generate a translation vector δ ∈ -δ max , δ max 2 , with δ max ∈ N + the maximum move distance.

Vector δ is sampled from δ ∼ g(δ), with:

g(δ) = exp(-w pos V pos [⌊y⌋ + δ]) ∑ δ ′ ∈ -δmax ,δmax 2 exp(-w pos V pos [⌊y⌋ + δ ′ ]) , (5.2)

CHAPTER 5 -Point Processes as Energy Based Models

With ⌊y⌋ the integer coordinates vector of point y. In practice, we build a tensor D y,X of size (2δ max + 1, 2δ max + 1) such that:

D y,X [i, j] = 1 Z exp ( -w pos V pos [ ⌊y⌋ + [ i j ]]) , i, j ∈ -δ max , δ max 2 , Z = ∑ δ ′ ∈ -δmax ,δmax 2 exp(-w pos V pos [ ⌊y⌋ + δ ′ ]
).

(5.3) 2D tensor D y,X defines a density over the raster space that approximates g. Note that for the return move from y ′ = y + δ to δ, the density is

g ′ (δ ′ ) = D y ′ ,X [-δ i , -δ j ]
. 1 To ensure this kernel can access all real-valued positions within δ max , δ is drawn as δ ∼ D y,X , then we set: ]).

δ ← δ + [ U([0, 1]) U([0, 1
]

( 5.4) The change in densities g and g ′ gets balanced out in α(y, x). The resulting kernel has the following acceptance ratio and is illustrated in Figure 5.1:

α(y, x) = min { 1, h(x)D y ′ ,X [-⌊δ i ⌋ , -⌊δ j ⌋] h(y)D y,X [⌊δ i ⌋ , ⌊δ j ⌋] } . y y δ ∼ Dy,X δ ∼ D y ,X D y,X D y ,X

Birth with density

The uniform birth kernel introduced in Section 3.3.2.3 proposes points uniformly in S × M. However, the density and scattering of objects in the image might not be uniform, making a lot of the object proposals superfluous. In their work, [START_REF] Crȃciun | Stochastic geometry for automatic multiple object detection and tracking in remotely sensed high resolution image sequences[END_REF] use a binary map to restrict the space S into a smaller subspace; boats being most likely on water, the binary map corresponds to a detection of water bodies. In [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF], the authors use a map of potentials to build a proposal density for the birth kernel. This proposal density is built from a pre-computation of the contrast measure for numerous possible object positions and parameters (i.e. data potentials in that case).

A density birth kernel Q B that proposes a new point u to configuration y with density d(u), has the following Green ratio:

r(y, y ∪ {u}) = p D p B Q D (y ∪ {u} → y) Q B (y → y ∪ {u}) h(y ∪ {u}) h(y) = p D p B λ d(u)(n(y) + 1)
h(y ∪ {u}) h(y) .

(5.5)

As in [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF], we propose using the data energy to drive this density d using the pre-computed position and mark tensors V pos and V κ .

Discrete points space. The position and mark energy maps V pos and V κ , are raster maps (i.e. defined over integer positions), meanwhile the density has to be defined on S × M. As it is simple to sample in a discrete finite space, we define a discrete equivalent of the space S × M to sample in. The discretization stems from the pixel space in the image, and mark discretization performed in Section 4.2.3:

S d = 0, H × 0, W , M d = ∏ κ∈{a,b,α}
{c(κ maxκ min ) + κ min , c = 0, . . . , n κ } .

(5.6)

Sampling from truncated energy. Ideally we would propose a point u to the current configuration y, sampled with the marginal density p(u|y) as is done within Gibbs sampling. While we cannot easily compute the marginal density of u knowing y we can approximate it through the position potential of point y, with p(u|y) ∝ h(y ∪ {u})/h(y):

h(y ∪ {u}) h(y) = exp (U (y, X) -U (y ∪ {u}, X)) (5.7) = exp ⎛ ⎜ ⎜ ⎝ -V (u, N y {u} , X) ⏞ ⏟⏟ ⏞ added point energy + ∑ y∈N y {u} V (y, N y {y} , X) -V (y, N y∪{u} {y} , X) ⏞ ⏟⏟ ⏞ energy variation of existing points ⎞ ⎟ ⎟ ⎠ (5.8) ≃ exp ( -V (u, N y {u} , X)
)

(5.9)

≃ exp ⎛ ⎝ -w pos V pos (u, X) - ∑ κ∈{a,b,α} w κ V κ (u, X) ⎞ ⎠
(5.10)

In ( 5.9) we first neglect the energy variations of the existing points, moreover the cardinality of N y {u} is limited by the number of object in a d max radius. In (5.10) we simplify the energy of point u to the easy-to-sample energy terms: the position and marks potentials. While for instance the non-overlap prior is of great importance, it makes the sampling of a new point u more computationally expensive. That way our truncated law for point proposal only depends on the image X.

We derive the following density for sampling new points ρ in discrete space

S d × M d :  d(ρ) = 1 Z exp ( - w pos V pos [ρ] + ∑ κ∈{a,b,α} w κ V κ [ρ] T d )
, (5.11) where Z is a normalizing constant (which we can compute easily as ˜ d is defined over a discrete finite space). Density ˜ d(ρ) can be computed for all pixels ρ in S d ×M d , as the normalizing constant is simply a sum over HW n 3 κ elements; thus this density can be easily sampled from. Temperature parameter T d > 0 is used to influence the density: high T d tends towards a uniform distribution, (i.e. getting closer to the uniform birth kernel); low T d makes the birth kernel only propose low data potential points (i.e. likely higher acceptance rate at the cost of less space exploration). Eventually, to sample a point u in continuous space S × M we proceed as follows:

1. Sample a pixel ρ in S d × M d with density ˜ d.

Sample u in ρ uniformly (u ∈ S × M).

This resulting sampling density d of u is:

d(u) = 1 |ρ u | ˜ d(ρ u ),
(5.12)

with ρ u the pixel in S d × M d containing u, and |ρ u | the measure of pixel ρ u (= 1 if the unit of measure is set to a pixel).

Jump diffusion

For sampling the Point Process we use the Jump Diffusion mechanism introduced in Section 3.3.3. It alternates diffusion with birth and death kernels to explore the whole configuration space Y. As a reminder the discrete diffusion as introduced (3.68) gives:

y t+1 = y t -γ∇U (y t ) + √ 2T (t)w t , w t ∼ N (0, γ).
The gradient of the energy ∇U is computed thanks to an automatic differentiation engine (see following paragraph on implementation 5.1.2), lifting the burden of manual gradient derivation, and allowing to implement and test new energy functions easily (as long as those are defined using the set of operations for which the automatic differentiation is available). Diffusion can be seen as a continuation of the local perturbation based on data. While with the latter we approximate the energy locally with the position potential to sample a new position for a point, the diffusion approximates the energy locally as U (y + δ) ≃ U (y) + ∂U (y) ∂y δ and samples a small perturbation in the vicinity of y.

For most of the energy model, the gradient computation is well-defined: for both single point priors (V e , e ∈ ξ prior-point ) and interaction priors (V e , e ∈ ξ prior-interact. ) minus the discontinuities discussed in Section 4.4.2.1. For the data terms (V e , e ∈ ξ data ), the computation of the energy corresponds to looking up a value at a position in a tensor Z. Thanks to the interpolation used to obtain values for float pixel coordinates, the operation is made differentiable, both relative to the configuration y and relative to the potential map Z2 .

In practice, as parallelization requires some maximum spatial displacement on points (see next Section 5.1.2), the step ∆C c (x t ) from x t to x t+1 (Alg. 5.1 line 9) is clipped; For each point y ∈ x t updated to y ′ , we bound the i and j components of ∆C c (x t ) within [-δ max , δ max ]:

∆C c (x t ) i ← min (δ max , max (-δ max , ∆C c (x t ) i )) , ensuring |y i -y ′ i | ≤ δ max ,
and similarly for j.

Parallel sampling

In Section 3.3.2.4 we define the basics of parallelization of the Point Process samples, which allows performing perturbations over multiple independent cells at once. With parallelization, we aim to take advantage of the spatial Markovianity of the Point Process; i.e., points further than a certain distance have no effect on each other's energy. Contrary to [START_REF] Verdié | Detecting parametric objects in large scenes by Monte Carlo sampling[END_REF] we use a constant cell size instead of quadtrees.

We build each cell to be a square of size d c , each cell is assigned to a mic-set s, such that no two cells in set s are neighbors (considering the 8-connectivity). To build such cells on a 2D plane requires four colors as illustrated in Figure 5.2. We demonstrate in Appendix E.1 that the minimum cell size necessary to maintain independence of the moves in cells of the same mic-sec needs to be such that:

d c ≥ 2d max + 2δ max .
(5.13)

To simulate the MCMC, one could pick a mic-set s uniformly and run the perturbations in parallel over each cell c in set s. However, we might do more perturbations than necessary in some low density areas of the image. Thus, we propose a selection scheme for the cells to compute, that aims at reproducing sampling density d introduced in Section 5.1.1.2. For a given configuration y we do:

1. Pick a kernel Q m with probability Q m (y, Y),
2. pick one mic-set s with probability p(s), 3. keep each cell in mic-set s with probability p(c|s) to form s 4.

for every cell c in set s

˜run kernel Q m restricted to c (denoted Q c,m ).
The process of cell selection (step 3), allows to limit the number of cells processed at once on big images, hence limiting the computational cost while maintaining the desired sampling density d. Probabilities of picking cells. We now build the probabilities p(s) and p(c|s) with which cells are picked, with the aim to fit the sampling density d introduced in (5.12). For a density birth kernel Q c,B , the kernel density over the whole image Q B is now:

d(u) λ ⏞ ⏟⏟ ⏞ Q B (y→y∪{u}) = p(s u )p(c u |s u )Q cu,B (y → y ∪ {u}), (5.14)
with c u and s u the cell and mic-set containing u respectively. For the above defined procedure we have the following conditions stemming from it:

• For step 2, we need p(s) ∈ [0, 1] and

∑ s p(s) = 1.
• For step 3, we need p(c|s)

∈ [0, 1]. • Q cu,B is of the form Q cu,B (y → y ∪ {u}) = dc(u) λ , with d c (u) density we can sample from, with ∫ u∈c d c (u)du = 1.
We propose the following solution that verifies (5.14) and the conditions mentioned above:

p(s) = ∫ v∈s d(v)dv = d(s),
(5.15)

p(c|s) = ∫ v∈c d(v)dv ∫ v∈s d(v)dv = d(c) d(s) , ( 5.16
) Skewing distributions for more parallelism. While the procedure to pick cells s ˜allows sampling points with density d, by construction the expected number of cells in s ˜is 1. In practice, we skew probability p(c|s) to ensure more cells are kept for simulation, by replacing p(c|s) by p ′ (c|s) = min (1, np(c|s)), with n > 1. It increases the amount of parallel computations of cells (thus reducing simulation time) at the cost of deteriorating the approximation of global sampling density d.

d c (u) = d(u) ∫ v∈c d(v)dv = d(u) d(c) , ( 5.17
Acceptance ratio in cells. From the mutual independence of cells we get the following property: the acceptance ratio for a move in cell c can be simplified (at least) as such:

α(y, y ′ ) = α(C c ¯(y) , C c ¯(y ′ ) ), (5.18)
where C c (y) is the set of points from configuration y in cell c, and c ¯the 8 cells neighboring c and c itself (thus by extension C c ¯(y) is points in y inside c or its neighbors). More details in Appendix E.2. This allows making the computation of acceptance ratio local, thus easy to perform in parallel.

Implementation. Instead of using multiple CPU threads as in [START_REF] Crȃciun | Stochastic geometry for automatic multiple object detection and tracking in remotely sensed high resolution image sequences[END_REF], we defer the parallelization to the implicit multithreading capabilities of GPU (Graphical Processing Unit) when performing matrix/tensor operations. In short, we write the whole energy model in terms of tensor operations, using an extra batch dimension to process multiple cells at once. Libraries such as (PyTorch, n.d.) combined with CUDA (CUDA, 2017), allow to transparently run tensor operations over multiple threads within a GPU. On top of it, granted operations are defined properly, PyTorch allows to automatically compute the gradient of the energy model; either relative to parameters θ -which we use for parameters estimation in Section 5.2.2 -or to the configuration y, which are used for diffusion in Section 3.3.3. We go further into details about the data structure and implementation of our approach in Appendix I.

Sampling method

The sampling method is summarized into Algorithm 5.1. The input for this algorithm are:

• x 0 : initial configuration;
• X: image;

• θ: energy model parameters;

• T 0 : initial temperature;

• n s : number of samples;

• α T : temperature decay rate (see annealing in Section 3.3.4).

CHAPTER 5 -Point Processes as Energy Based Models

Input: x 0 , X, θ, T 0 , n s , α T 1: for t = 0, . . . , n s -1 do 2:

Pick diffusion with probability 0.8, else jump 3:

Pick mic-set s with probability p(s)

4:
Keep each c in s with probability p(c|s) to make s 5: dw ∼ N (0, γ)

x t+1 ← x t /*
9: ∆C c (x t ) = -γ∇ Cc(xt) U (x t , X, θ) + dw √ 2T t 10: C c (x t+1 ) ← C c (x t ) + ∆C c (x t ) 11: else 12: Q c ← Q c,B with probability 0.5 else Q c,D
/* pick birth or death */ 13:

x ′ ∼ Q c (x → •) 14: r ← Qc(x ′ →x) Qc(x→x ′ ) exp ( -U (x ′ ,X,θ)-U (x,X,θ) Tt ) /* compute Green ratio */ 15: C c (x t+1 ) ← C c (x ′ ) with probability min(1, r) 16: end if 17:
end for 18:

T t+1 ← α T T t 19: end for Output: x ns Algorithm 5.1: Sampling method Sample(x 0 , X, θ, T 0 , n s , α T ).

Alternatively, in Algorithm 5.1 we can run the diffusion kernel interrupted by jumps (birth or death) after a wait time ω, with ω following a Poisson distribution:

ω ∼ p(ω) = τ ω ω! exp(-τ ), (5.19)
where τ = E(ω) represents the expected waiting time thus controls the frequency of the jumps [START_REF] Descombes | Stochastic geometry for image analysis[END_REF].

In this section we first introduced using the potential map provided by the CNN to first build a local transform kernel based on the position potential map V pos . Then we improved on [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF], by building a birth map from truncated model given by the CNN output instead of precomputing contrast for numerous objects in the image. In the continuation of data driven kernels, we implement the diffusion mechanism, by proposing to leverage the modern automatic differentiation engines. Finally, we adapted the parallel approach from [START_REF] Verdié | Efficient Monte Carlo sampler for detecting parametric objects in large scenes[END_REF], proposing a cell selection scheme that simulates the Point Process in parallel given the density derived from the potential map.

Parameters estimation

Our model introduced throughout Chapter 4 contains several parameters that require to be set before sampling configurations. The most important of which being the relative weights on the 5.2 -5.2.1 Learning energy weights with local perturbations 103 energies introduced in (4.3). We also have a few parameters in the various priors of the model. We denote by θ the set of parameters and by U (y, X, θ) the resulting energy model. The ideal parameters θ * should be such that the ground truth configuration y GT corresponds to the minimum of the energy U for an image X:

y GT = arg min y∈Y U (y, X, θ * ) .
(5.20)

We first present some previous approaches for Point Process parameter estimation, for instance using linear constraints. Then, we propose a method to learn the weights of the energies θ = {w 0 } ∪ {w e , e ∈ ξ} based on margin maximization, solving the over-constraint issue of the previous approach from (Q. [START_REF] Yu | Multiple-Target Tracking by Spatiotemporal Monte Carlo Markov Chain Data Association[END_REF][START_REF] Craciun | Joint Detection and Tracking of Moving Objects Using Spatio-temporal Marked Point Processes[END_REF]. This was published in (Mabon et al., 2021). Eventually we adapt learning methods from energy based models, using contrastive divergence [START_REF] Hinton | Training Products of Experts by Minimizing Contrastive Divergence[END_REF][START_REF] Teh | Energy-based models for sparse overcomplete representations[END_REF] with a replay buffer [START_REF] Du | Implicit generation and modeling with energy based models[END_REF] to estimate the Point Process parameters with gradient descent. This approach was first published in (Mabon et al., 2022a) andrefined in (Mabon et al., 2023a) and (Mabon et al., 2023b).

Learning energy weights with local perturbations

Linear Programming

In their work [START_REF] Crȃciun | Stochastic geometry for automatic multiple object detection and tracking in remotely sensed high resolution image sequences[END_REF] and (Q. [START_REF] Yu | Multiple-Target Tracking by Spatiotemporal Monte Carlo Markov Chain Data Association[END_REF] use random perturbations on the Ground Truth y GT to generate non-valid configurations and build a set of linear constraints to estimate the weights w e . As a matter of fact, re-formulating (5.20) locally yields:

U ( y GT , X, θ * ) < U ( y -, X, θ * ) y -∼ Q - ( y GT , • ) , (5.21)
where y -is a negative configuration, generated from perturbation kernel Q -. The kernel Q -is build such that the perturbation on y GT results in a perceptually worse configuration. For instance [START_REF] Crȃciun | Stochastic geometry for automatic multiple object detection and tracking in remotely sensed high resolution image sequences[END_REF] use random translation, rotations and scaling (with a great amplitude) along with random addition and removal of points.

Once a certain number of linear constraints are built, [START_REF] Crȃciun | Stochastic geometry for automatic multiple object detection and tracking in remotely sensed high resolution image sequences[END_REF] and (Q. [START_REF] Yu | Multiple-Target Tracking by Spatiotemporal Monte Carlo Markov Chain Data Association[END_REF] use Linear Programming to solve for the weights w e . However, this approach has several limitations:

First, the number of constraints are to be limited, otherwise this leads to an over-constrained unsolvable problem. The latter couples badly with the imprecise nature of the ground truth; with noisy ground truth, some linear constraints can become contradictory, thus making the problem unsolvable. With noisy input data, one would resort to use more data to compensate the variance of the input; the latter Linear Programming approach does not offer this possibility. We propose to address this issue in the following Section 5.2.1.2.

Finally, how the negative sample y -is built from to the ground truth depends on the design of this kernel Q -. In practice, we find we need to iterate between crafting the perturbation Q - and solving for w e . For instance, given a set of negative samples y -, we obtain weights w e and find out that, at inference, it exhibits some bad characteristics (e.g. too much overlap). We then proceeded to build a new set of negative configurations y -via a new Q -, with more examples to Figures are adapted from [START_REF] Cortes | Support-vector networks[END_REF] and [START_REF] Larhmam | Maximum-margin hyperplane and margin for an SVM trained on two classes[END_REF]. Support Vector Machine. With the Support Vector Machine (SVM) [START_REF] Cortes | Support-vector networks[END_REF] propose to find the hyperplane that maximizes the gap in between the two categories within a classification problem. We consider a set of samples y k , k = 1, . . . , n with associated vectors v k , with labels l k = -1 for positive samples and l k = 1 for negative samples 3 .

For linearly separable data, the Hard-margin criterion minimize ∥w∥ 2 2 (i.e. maximizes the gap 2/ ∥w∥ 2 ) with constraint l k (w • v kb) ≥ 1 for all k = 1, . . . , n; i.e. all configurations lie in the correct side of the hyperplane (see Figure 5.3 (c)).

However, our data might not be separable: [START_REF] Cortes | Support-vector networks[END_REF] propose to use the hinge loss and minimize the following instead:

γ ∥w∥ 2 2 + 1 n n ∑ k=1 max (0, 1 -l k (w • v k -b)) , ( 5.23)
with γ > 0 a scalar setting the balance between increasing the margin and having points in the correct side.

Remark 5.2.1 -SVM are often used with nonlinear kernels, allowing for nonlinear separation functions. In our case the energy model is linear; we discuss the addition of nonlinearities in the energy mixing function in Section 4.4.2.2.

Application and limitations. Minimizing this loss on the training data gives us the weights w and b. In our case, parameter b in (5.23) is superfluous as the energy needs to be known up a to a constant. This method allows learning weights w e even from non-separable set of positive and negative samples. Moreover, we can consider many samples without risking to over-constraint the optimization. However, this method is still reliant on the two perturbation kernels Q -and Q + that respectively generate the negative and positive samples. Finally, this does not allow to tune the energy term parameters such as t pos , µ area , t ovrlp , . . .. Those have to be estimated separately or set manually.

Maximum likelihood learning

In this section we approach our parameter estimation problem from the angle of Energy Based Models (EBM). Definition 5.2.1. Energy Based Model (EBM) is a form of generative model; it captures the dependencies between variables by associating a scalar energy to each configuration of variables. Importantly it allows for an implicit sample generation procedure, where sample x is found from x ∼ exp(-U (x)). Learning the EBM consists in finding the energy function that associates low energies to correct values of the variables. [START_REF] Lecun | A Tutorial on Energy-Based Learning[END_REF][START_REF] Du | Implicit generation and modeling with energy based models[END_REF].

CHAPTER 5 -Point Processes as Energy Based Models

Formally, a Gibbs Point Process model is an Energy Based Model. The works of [START_REF] Lecun | A Tutorial on Energy-Based Learning[END_REF] propose to learn EBMs by maximizing the likelihood for the data D:

p ( y GT 1 , . . . , y GT n D |X 1 , . . . , X n D , θ ) = n D ∏ k=1 p ( y GT k |X k , θ )
.

(5.24)

The latter can be derived into the Negative Log-Likelihood (NLL) loss:

L NLL (θ, D) = 1 n D n D ∑ k=1 ( U ( y GT k , X k , θ ) + β log ∫ y∈Y exp ( -β -1 U (y, X k , θ) ) )
, (5.25) with β a temperature parameter (from the Gibbs distribution), that has no effect on the position of the minimum. While the left part U (y GT k , X k , θ) pulls down the energies of ground truth configurations, the right part U (y, X k , θ) pulls up the energies of all configurations in Y.

When minimizing the loss via gradient descent, the gradient at step n for parameters θ n and data sample k = 1, . . . , n D , is given as:

∂L NLL ( θ n , y GT k , X k ) ∂θ n = ∂U ( y GT k , X k , θ n ) ∂θ n - ∫ y∈Y ∂U (y, X k , θ n ) ∂θ n p (y|X k , θ n ) ⏞ ⏟⏟ ⏞ I k,n
, (5.26) where p(y|X i , θ n ) is given from the Gibbs distribution:

p(y|X i , θ n ) = exp(-β -1 U (y, X k , θ n )) ∫ y ′ ∈Y exp(-β -1 U (y ′ , X k , θ n ))
.

(5.27)

While the integral I k,n remains intractable, it can be approximated through Monte Carlo sampling, where ˜ y 1 , . . . , ˜ y ns are n s samples drawn from the law defined by p(•|X k , θ n ), yielding:

I k,n ≃ 1 n s ns ∑ s=1 ∂U (˜ y s , X k , θ n ) ∂θ n .
(5.28)

As the law defined by p(•|X k , θ n ) is simply the Gibbs Point Process with energy U (y, X k , θ n ) at temperature β, we can use the Point Process sampling procedure we introduced previously.

Contrastive divergence

Hinton et al. in [START_REF] Hinton | Training Products of Experts by Minimizing Contrastive Divergence[END_REF][START_REF] Teh | Energy-based models for sparse overcomplete representations[END_REF] propose to use a single sample in his contrastive divergence method. This method also uses few simulation steps for the Monte Carlo Markov chain (MCMC) to generate y -, starting from the desired answer y GT . The general idea is to generate contrastive samples y -that follow the density derived from U (•, X, θ n ) at step n of the optimization. Then we proceed to update θ n to θ n+1 , by gradient descent, with the goal to minimize the energy of the valid sample y GT , while maximizing the energy of the contrastive sample y -(see Figure 5.4). Alternatively we can augment the data and use positive samples y + = y GT + N (0, σ + ) to replace y GT as in [START_REF] Du | Implicit generation and modeling with energy based models[END_REF]. The loss to minimize is then: with γ > 0 the weight of regularization term R. We introduce the regularization term to avoid an explosion of the per-point energy V (y, X, N y {y} , θ n ). To ensure for a sparse weighting of energies (i.e. minimize the number of non-zero weights w e ) we can introduce a new regularization term as an L 1 norm on the vectors of weights [START_REF] Goodfellow | Deep learning[END_REF]:

L ( θ n , y + , y -, X ) = U (y + , X, θ n ) -U (y -, X, θ n ) + γR V R V = ∑ y∈{y + ,y -} 1 |Y | ∑ y∈y ⃓ ⃓ ⃓V ( y, X, N y {y} , θ n )⃓ ⃓ ⃓ , ( 5
R 1 = ∑ e∈ξ |w e |
(5.30)

The broad strokes of the estimation procedure for parameters θ are as follows:

1. Pick a pair (y GT , X) from data D.

2. Generate positive sample y + = y GT + N (0, σ + ).

Generate negative sample y

-∼ exp(-β -1 U (y -, X, θ n )).

Compute the loss

L ( θ n , y + , y -, X )
(see (5.29)).

5. Update θ n to θ n+1 according to the gradient ∇L, with the Stochastic Gradient Descent (SGD) method [START_REF] Bottou | Stochastic gradient descent tricks[END_REF]. While the local perturbation method (Section 5.2.1) only learns from constraints with states in the neighborhood of the Ground Truth y GT (or positive samples y + ), the contrastive divergence method allows learning by contrast against any state in Y. It also focuses the contrastive samples on low energy states (instead of uniformly sampling in Y which would be inefficient).

Replay buffer

As sampling contrastive samples y -can be time-consuming, [START_REF] Du | Implicit generation and modeling with energy based models[END_REF] propose an adaptation of Hinton's method [START_REF] Hinton | Training Products of Experts by Minimizing Contrastive Divergence[END_REF][START_REF] Teh | Energy-based models for sparse overcomplete representations[END_REF], making use of a replay buffer. The replay buffer saves Markov chain results at current optimization step, to use for initialization in the next steps, thus saving computing time. This allows reducing the long simulation time necessary to pass the burn-in period of the Markov Chain [START_REF] Robert | Diagnosing Convergence[END_REF]) and get samples y -. We use a sample from the law derived from U (•, X, θ n-1 ) to initialize the chain that samples the law derived from U (•, X, θ n ).

The use of the replay buffer within one step n of the optimization loop goes as follows:

1. With probability p B set configuration x 0 as a value of the replay buffer B picked at random. With probability 1p B (or if the buffer is empty) set configuration x 0 as a random configuration in Y.

2. Run the Markov Chain (x t ) t=0,...,ns to simulate model of energy U (•, X, θ n ).

3. Use y -= x ns as a negative sample for the loss computation and update of θ n to θ n+1 .

4. Update buffer B ← B ∪ {y -}.

The buffer allows to virtually run longer chains across the epochs 4 as shown in Figure 5.5 (a). At each epoch we pick the result from a previous chain and continue it for a n s steps before saving the result to the buffer. We find it is useful to limit the size of the buffer to n B , and update it on a first in first out fashion; with lower size n B , the Markov chain starts more often from recent contrastive samples, i.e. likely closer to an energy minimum. At the limit, with n B = 1, the contrastive sample generations pick up the last sample y -, and updates it with n s MCMC steps (see Figure 5.5 (b)). This yields an expect cross-epoch chain length of nsp B 1-p B , allowing to run longer chains with fewer steps. This corresponds to the persistent contrastive divergence proposed by [START_REF] Tieleman | Training restricted Boltzmann machines using approximations to the likelihood gradient[END_REF]. We discuss the effects of the choice of p B and n s on the cross-epoch chain in Appendix F. Remark 5.2.2 -Contrary to [START_REF] Du | Implicit generation and modeling with energy based models[END_REF] -where the authors perform image generation -our generated configurations y -depend on the current image X k . Thus, we define the buffer B as a set of sub-buffers B X k , each associated to an image X k . A configuration y -generated from image X k get added to sub-buffer B X k .

Effect of temperature on samples

In the contrastive divergence method (Section 5.2.2.1), we introduce a temperature parameter β to generate the negative samples y -. This temperature parameter is the same as the one introduced for simulated annealing (Section 3.3.4) and has the same effects on the generated samples. At high temperature the samples get closer to a uniform distribution in Y, and as temperature lowers samples get closer to local or global minima. Within the contrastive divergence optimization scheme:

• Picking a high β allows to increase the energy of configurations far away from y GT . Samples are not likely to be drawn close to y GT , thus the learned model will not discriminate efficiently around y GT .

• At low β sample are drawn consistently in local minima around y GT , the ability to discriminate against configurations that differ greatly from y GT will be diminished. We propose using a set of β, {β 1 , β 2 , β 3 } corresponding to low, medium and high temperatures, allowing to mix samples explore Y while drawing some samples near local minima or U (y, X, θ n ).

epoch n = 4 B X1 at n B X2 at n epochs image samples X 1 : y - 1,1 U(Y) K ns θ1 y - 1,2 U(Y) K ns θ2 y - 1,3 K ns θ3 y - 1,4 K ns θ4 y - 1,5 K ns θ5 X 2 : y - 2,1 U(Y) K ns θ1 y - 2,2 K ns θ2 y - 2,3 K ns θ3 y - 2,4 K ns θ4 y - 2,5 K ns θ5 (a) Unlimited buffer size epoch n = 4 B X1 at n B X2 at n epochs image samples X 1 : y - 1,1 U(Y) K ns θ1 y - 1,2 U(Y) K ns θ2 y - 1,3 K ns θ3 y - 1,4 K ns θ4 y - 1,5 K ns θ5 X 2 : y - 2,1 U(Y) K ns θ1 y - 2,2 K ns θ2 y - 2,3 K ns θ3 y - 2,4 K ns θ4 y - 2
Remark 5.2.3 -As for the previous remark, we further separate the buffer to associate each sub-buffer to a temperature β. Now the buffer corresponding to image X and temperature β is B X,β .

Remark 5.2.4 -When estimating weights w 0 , w e we notice that scaling the weights by a factor ς, is equivalent to scaling the temperature by ς -1 :

U (y, X, {w 0 } ∪ {w e , e ∈ ξ}) T ς -1 = U (y, X, {ςw 0 } ∪ {ςw e , e ∈ ξ}) T .
(5.31)

Thus picking, a temperature might have different effect depending on the weights. In practice, we find that scaling the energy model to unit variance solves this consistency issue; we compute the energy model variance given its parameters θ:

Var(U T =∞ ) = ⟨U 2 T =∞ ⟩ -⟨U T =∞ ⟩ 2 , ( 5.32)
with ⟨U T =∞ ⟩ the mean of energy U over multiple samples drawn at infinite temperature (i.e. drawn from U(Y)).

We then simulate the energy model

U ¯= U/ √ Var(U T =∞ ) (so that Var(U ¯T =∞ ) = 1
). This is equivalent to specifying the temperature for a unit variance energy model and scaling it to the current energy model it is being applied on. The simple procedure to compute the scaling factor is presented in Algorithm 5.2, where we pick 10000 as an arbitrary large number of samples to compute the variance on.

Input: θ, X

1: Values ← ∅ 2: for all k = 1, . . . , 10000 do 3:

x ∼ U(Y)

/* sample random configuration */

4:

Values ← Values ∪ {U (x, X, θ)} 5: end for 6: ς ← √ Var(Values) Output: ς Algorithm 5.2: Computation of the energy scaling factor, ComputeScaling.

Algorithm

Resulting from the above, we have the following algorithm for inferring parameters ˆ θ: for all (y GT , X) in D do y -← Sample (x 0 , X, θ n , ςβ, K, 1) /* see Sample(x 0 , X, θ, T 0 , n s , α T ) in Alg.

*/

8:

y + ← y GT + N (0, σ + ) 9: ∆θ n ← ∇ θn L ( θ n , y + , y -, X ) 10:
Update θ n+1 with ∆θ n using SGD 11:

B X,β ← {y -} ∪ B X,β
/* update the buffer */ 12:

n ← n + 1 13:
end for 14: end while Algorithm 5.3: Contrastive divergence for parameters estimation.

The sampling method at line 7 is described in Algorithm 5.1. In practice one iteration of the loop on data at line 3 is done on a mini batch of ground truth and image pair instead of a single pair; it allows smoothing the gradient between iterations by averaging on data samples.

In this section we contribute in proposing to use contrastive divergence (adapted from energy based models) to estimate the parameters of the Point Process model. The use of replay buffer (or persistent contrastive divergence) allows cutting on sampling time when sampling negative configurations. While previous approaches using linear constraints and linear programming would often be limited to learning energy weights, and risk over-constraint with too noisy or too many data points, our method does not suffer from over-constraint, and allows solving for all parameters of the model at once.

Discussions

Relating contrastive divergence to manual methods. In practice, the contrastive divergence method relates to the empirical trial and error procedure used in previous works to set ˆ θ. Indeed, oftentimes the manual procedure to find weights w e goes as follows:

1. the user sets some approximate values to weights w ˜e within ˜ θ,

the user simulates configurations y

˜∼ U (•|X, ˜ θ),
3. if the inference is of good quality, the procedure ends. Otherwise, the user updates values for weights w ˜e, picked to counteract the observed bad configurations; (e.g., the configuration has too much overlap → increase the weight of the overlap energy). Then go back to step 2.

In a broad sense, the repeat of steps 2 and 3 maps to the alternated sampling (Alg. 5.3 line 7) and loss minimization (Alg. 5.3 lines 9 to 10). As step 3 in the manual procedure, attempts to maximize the energy of the (bad) sampled configuration, minimizing the loss in Algorithm 5.3 tends to increase U (y -, X, θ n ) and decrease U (y + , X, θ n ).

Learning more than just weights. For previous parameter inference method, we could only infer the weights of the linear combination of energy terms. Here, the only requirement to be able to learn a parameter is to be able to back propagate the energy model gradient up to that parameter. This allows to learn not only weights but also energy term parameters so that: θ = {w 0 } ∪ {w e , e ∈ ξ} ∪ {t pos , µ area , σ area , µ ratio , µ ratio , t ovrlp , t repls , t attrc }.

We detail in Appendix I the implementation allowing such back propagation of the gradient.

It would be possible to infer the parameters of the whole data energy model (i.e., the backbone CNN), however this proves quite computationally expensive; training the CNN requires many epochs as parameters are numerous (> 10 6 parameters). Meanwhile, the contrastive divergence method (Algorithm 5.3) requires simulating a few Markov Chain steps in between epochs to generate contrastive samples, which increases the computational complexity of each epoch.

Papangelou intensity as a score

Computing the detection score

Classical CNN models for object detection yield a confidence score s(y) ∈ R for each proposed object y in the image. This confidence score is often interpreted, for each detection, as proportional to the probability of proposed element y to be a true object (true positive), s(y) ∝ p(y|X). Applying a score threshold t s gives a set of detections, for which metrics such as precision and recall can be computed by matching the detections with the ground truth. As the threshold has to be adapted according to the need of the application of the model; i.e. some applications may require few false positive (high precision) while other need less missed detections (high recall). To assess the performance independently of the threshold selection, the Average Precision (AP) metrics sums up the performance as the area under the precision-recall curve.

With the PP approach, the probability of one proposed point being an object of interest depends on the rest of the inferred configuration y ˆ, thus the scoring function should too s(y|y ˆ\ {y}) ∝ p(y|y ˆ\ {y}, X). From Definition 3.1.7 we have that the Papangelou intensity is proportional to the probability of finding a point y ∈ y in a small neighborhood dy knowing the rest of the configuration y \ {y}. We propose to use the Papangelou intensity as a score :

λ(y; y \ {y}) = h(y) h(y \ {y}) = exp ( U (y \ {y}, X, ˆ θ) -U (y, X, ˆ θ) )
.

(5.33)

Pruning sequence. However, the dependency of the score on the current configuration yields a complication while computing the Average Precision: when applying a threshold t s to prune the configuration y into y ′ ⊂ y, for any y ∈ y ′ , the score s(y|y \ {y}) may differ from s(y|y ′ \ {y}).

With a score of the form s(y), that only depends on y and the image X -such as those from classical CNN models -the score from one object after pruning is unchanged.

In the PP case, we compute the scores by sequentially removing the lowest scoring point until none is left; i.e., we build a sequence of configurations y 1 ⊃ y 2 . . . y n(y ˆ)-1 ⊃ y n(y ˆ) ⊃ ∅, with Case 3: Arbitrary sequence of moves. We now consider the generic case where y -is generated from a sequence of moves from y + . Note that a local transformation (translation, rotation, scaling etc...) can be formulated as a death followed by a birth. We built a sequence (y k ) k=0,...,n of n configurations as:

∀k = 1, . . . , n, y k = { y k-1 \ {y k } if y k ∈ y + y k-1 ∪ {y k } otherwise,
with y 0 = y + , y -= y n , and y k elements of either S × M or y +5 . Without loss of generality we can reorder the sequence to match the pruning order defined in (5.34). The energy change for one move is given as:

∆U (y k-1 → y k ) = { log(λ(y k ; y k-1 \ {y k })) if y k ∈ y + -log(λ(y k ; y k-1 ))
otherwise.

As we have (by definition) ∆U (x → x ′′ ) = ∆U (x → x ′ ) + ∆U (x ′ → x ′′ ), the loss is given as:

L(θ, y + , y -, X) = - n ∑ k=1 ∆U (y k-1 → y k ) = ∑ y k / ∈y + log(λ(y k ; y k-1 ) ⏞ ⏟⏟ ⏞ (a)
) -

∑ y k ∈y + log(λ(y k ; y k-1 \ {y k }) ⏞ ⏟⏟ ⏞ (b)
).

With ordering of the y k , y k matching the pruning order in (5.34) each λ(y k ; . . . ) can be matched to their respective score:

• part (a) corresponds to non-valid points added to y + , their score is minimized as the loss is decreased;

• part (b) corresponds to valid points removed from y + , their score is increased as the loss is minimized.

Hereby we show that minimization of the loss at a configuration level leads to the expected results on object scores.

Conclusion

In this chapter we presented our methods to sample and estimate parameters of the Point Process model. We summarize our contributions as follows:

• Adapting data-based perturbation kernels to our pre-computed energy maps: both for local perturbation kernels (Section 5.1.1.1) and Birth kernels by sampling from a truncated energy function (Section 5.1.1.2).

• Adapting parallel computation to focus sampling on cells with estimated high density of object (Section 5.1.2).

• A simple procedure to estimate energy weight parameters, based on SVM, that avoids the over-constraint issue of previous approaches (Section 5.2.1.2).

• Another parameter estimation method based on contrastive divergence initially introduced for EBM. It allows learning internal energy parameters along with energy weights in a unified procedure. We adapt the replay buffer to the Point Process model, in order to shorten sampling times during parameter estimation (Section 5.2.2).

• Using the Papangelou intensity as scoring function for the inferred points, to compute benchmark metrics (Section 5.3).

CHAPTER 6

Experimental results In this chapter we apply on both synthetic and real data the models and algorithms presented in the thesis. First we present the application of the contrastive divergence on a simple example using synthetic data to demonstrate its capabilities. We then show our first results with non-marked Point Processes on a first dataset: Cars Overhead With Context. A more advanced model is applied on DOTA, which shows great results compared to other CNN based approaches; inferred configurations are more regularized and our models are more resilient to additive noise. In the following part we demonstrate a use of the contrastive divergence method for parameter estimation in a simple case, and how it is able to properly weight the energy terms within the energy model; i.e. the estimated weights reflect the relative importance of energy terms within the generative model. We proceed as follows:

Dans

1. Define an energy model.

2. Sample some configurations given a set of energy weight parameters.

3. Estimate the energy weight parameters from the generated data in order to simulate similar data.

We consider a Point Process of circles, where each point y ∈ y is such that y ∈ S × M = [0, H] × [0, W ] × [r min , r max ], with y = (y i , y j , y r ). The energy model is composed of the following energy terms:

V pos (y, X) = X[y i , y j ] (6.1)

V ovrlp (y, N y {y} ) = max u∈N y {y} { max ( 0, d(u, y) -y r -u r y r + u r )} (6.2)
V expnd (y) = y rr min r maxr min (6.3) V noise (y) = v y , v y ∼ N (0, 1) (6.4)

The position potential (6.1) samples the value of the image at the object center location, overlap potential (6.2) assigns high energy to overlapping circles, the expansion potential (6.3) favors circles of a larger size, and (6.4) gives a random energy to each point of the configuration (and as such is not a useful potential). To each we assign a corresponding weight w pos , w ovrlp , w expnd and w noise . We generate training samples by sampling the energy model with weights set as w noise = 0, w pos = -2, w ovrlp = 8, w expnd = -1, picked semi-arbitrarily to limit overlap while favoring big circles: some samples are shown in Figure 6.1 (a) and (b). With these image and configuration pairs, we now use the contrastive divergence method (Algorithm 5.3) to estimate the weights used for generating the configurations. We plot the energy for positive samples (y + ) and negative samples (y -) through the iterations of the training procedure (epochs) in Figure 6. 1 (d) and the inferred weights in (e). We see that the noise priorthat is irrelevant to the Point Process energy mixture -remains with an almost zero weight (even through multiple iterations of the inference procedure, see Figure 6.1 (f)). Although the resulting weights are not exactly the ones used to generate the training data, both models (the data generator model, and the model with inferred weights) produce similar results as shown in Figure 6.1 (c); both models seemingly have similar distributions at low temperature. Remark 6.1.1 -For usual gradient descent methods, the loss should decrease monotonically through epochs. In this case, the loss is computed against a set of contrastive samples that is constantly updated; at each epoch new contrastive samples are generated, hopefully closer to the ground truth configurations. Here, as the model converges, we expect the positive (U (y + )) c) shows a configuration inferred from the model using the estimated parameters. In (d) we plot the loss and its components (see (5.29)) through the estimation procedure steps. (e) shows the weights value through the estimation procedure, as in (f) where we compare over multiple random initialization.

6.2 -6.1.2 Stable perturbation kernel and negative (U (y -)) samples energies to progressively get closer. However, in Figure 6.1 (d) we observe a consistent offset between U (y + ) and U (y -) after epoch 30; this is due to some high temperatures β in the mix used to generated y -(see Section 5.2.2.3). At high temperatures the samples y -are drawn closer to a uniform distribution in Y, thus get rated with a higher energy.

Stable perturbation kernel

Within the contrastive divergence training procedure in Algorithm 5.3 the sampling of negative sample y -is very costly. One can wonder whether we could instead generate negative samples y - with a static kernel Q -instead; reducing computational time and complexity while maintaining the advantages of the gradient descent approach to estimate parameters θ.

In this ablation study, we place ourselves in a simple setting similar to the energy model and training data described in previous Section 6. 1.1. We propose three different experimental settings:

• Contrastive: using Algorithm 5.3 as such.

• Static uniform: replacing the costly line 7 in Algorithm 5.3 by y -∼ Q -(y + , •), where perturbation kernel Q -shifts location and marks of points with a random Gaussian, removes points at random, and proposes new points uniformly in S × M.

• Static density: same as previously, but new points are proposed with birth density d defined in (5.12).

We show results in Figure 6.2: from this simple text we see that Algorithm 5.3 fails to learn an energy that can generate proper point configurations when the negative samples are not generated with y -∼ exp(-U (•, X, θ)).

Results on optical data

Our application goal is the detection of small objects in images from satellites such as Pléiades-HR (Pleiades-HR, 2012) (0.7m resolution), Pléiades Neo (Pléiades Neo, n.d.) (0.3m resolution) or CO3D (CO3D Constellation, 2023) (0.5m resolution). Within the LiChIE project and our collaboration with Airbus Defense and Space, we set the resolution for our application to 0.5m. [START_REF] Mundhenk | A Large Contextual Dataset for Classification, Detection and Counting of Cars with Deep Learning[END_REF].

The dataset contains images of several urban areas, for which the centers of vehicles are labeled in their center. This dataset is quite diversified as it shows dense urban areas as well as forested areas or fields.

As we aim to develop detection model for satellite imagery, we sub-sample images to as resolution of 0.50m and refer to the resulting dataset as COWC 0.5 . We use an anti-aliasing filter when 6.2 -6.2.1 COWC with non-marked Point Process 123 sub-sampling in order to avoid artifacts on object edges and better approximate the acquisition of a lower resolution image. We only use the RGB images of the dataset (no grayscale). For each urban area we split it into two thirds for training and one third for testing.

Point Process model

As the dataset only contains object location and no geometric data, we consider a Point Process without marks y ∈ S.

The energy model contains the following energies, with corresponding weights:

• Position potential V pos as defined in (4.10).

• Overlap prior V ovrlp from (4.37), where each point is considered as a disk of constant diameter a set as the typical car width in pixels at resolution 0.5m.

• Triplet alignment potential V triAl defined in (4.45).

The total energy is then given as:

U (y, X) = ∑ y∈y w pos V pos (y, X) + w ovrlp V ovrlp (y, N y {y} ) + w triAl V triAl (y, N y {y} ). (6.5)

Parameter estimation and inference

To estimate the weights w pos , w ovrlp , w triAl , we use the linear SVM methods as described in Section 5.2.1.2. We generate positive samples y + from e ground truth configuration y GT for each element of the training data with random Gaussian perturbation of low amplitude. Negative samples y -are generated from uniform sampling of points or high amplitude Gaussian perturbation or random point removal from y GT . The estimate y ˆof optimal configuration y * = arg min y∈Y U (y, X) is obtained with the RJMCMC and simulated annealing method using the following kernels:

• Birth and Death kernels as introduced in 3.3.2.3 with non-uniform proposal density based on V pos .

• A transform kernel based on data as described in Section 5.1.1, using the energy map from V pos .

Results

Metrics. To measure the validity of an estimated point configuration y ˆagainst the ground truth y GT we try to match every ground truth center with estimated points within a radius r = a/2 as in [START_REF] Mundhenk | A Large Contextual Dataset for Classification, Detection and Counting of Cars with Deep Learning[END_REF]. Thus, an object is considered a detection if it matches to at least one point from y ˆ. DOTA. We compile DOTA 0.5 , a 0.5 meter per pixel vehicles in remote sensing dataset, from the DOTA dataset [START_REF] Xia | DOTA: A Large-Scale Dataset for Object Detection in Aerial Images[END_REF], and use it to train the CNN backbone, and infer the model parameters ˆ θ. This dataset is built by sub-sampling images from DOTA to the desired spatial resolution.

-6.2.2 Point Process of rectangles on DOTA and ADS data 125

To test the resilience of our models against noise, we test the methods on the same data adding Gaussian noise (with σ = 0.3 on input image X at inference, see Figure 6.4). We denote this dataset DOTA N 0.5 .

ADS data. Moreover, we evaluate models on data provided by Airbus Defense and Space (ADS) of two areas around Lyon, this aerial data is sub-sampled to the desired resolution (0.5m ), emulating the CO3D satellite characteristics. This dataset is unlabeled, thus will only serve to evaluate the models qualitatively, and was not used for training/parameters estimation.

Models

In this section we show results both on CNN based models, and PP models on the DOTA 0.5 dataset.

The PP models use the following energies from Table 6 To each energy term e we assign a respective weight w e as in (4.3). For some above-mentioned energies we use multiple instances; e.g. we have one mark energy for each mark κ ∈ {a, b, α}. Similarly, we have one alignment ratio for parallel objects (t align = 0) and one for perpendicular objects t align = π/2.

The models tested on the test dataset are as follows:

• CNN-LocalMax.: naive detection from the CNN backbone (used in the PP models); we find objects through local maxima in the output probability maps (or local minima in the energy maps).

• CNN-PP ♦ : PP model with minimal inferred parameters: inferred parameters θ are the energy weights {w e , e ∈ ξ}. The energy term parameters in the table above are set manually.

• CNN-PP ⋆ : PP model with parametrized priors: inferred parameters θ include both the energy weights {w e , e ∈ ξ} and the energy term parameters in the table above.

• BBA-Vec. and YOLOV5-OBB: lastly, we compare all of our models above with BBA-Vec. from [START_REF] Yi | Oriented Object Detection in Aerial Images with Box Boundary-Aware Vectors[END_REF] and YOLOV5-OBB from [START_REF] Yang | On the Arbitrary-Oriented Object Detection: Classification Based Approaches Revisited[END_REF][START_REF] Jocher | Ultralytics/yolov5: V7.0 -YOLOv5 SOTA realtime instance segmentation[END_REF]. These models are trained on the same data as the above-mentioned models.

Parameter estimation and inference

Estimation of energy parameters. The parameters θ for both PP methods CNN-PP ♦ and CNN-PP ⋆ are estimated with the Algorithm 5.3. The procedure is applied over 512 epochs (iterations on the train data which consists of 2048 ground truth patches of size (128, 128)). Temperatures for sampling negative/contrastive configurations are picked as β 1 = 1.0, β 2 = 10 -2 and β 3 = 10 -4 . We set the buffer size n B = 1 and its probability p B = 0.99. The regularization weight is set to γ = 1 for the loss in (5.29).

Inferring configurations from images. At inference, the configuration y ˆis obtained through Algorithm 5.1. We set a number of iterations per pixel n ipp = 2.0, with the expected number of cells being process in parallel n // = E[|s ˜|], we get the number of iterations for an image X of size (H, W ):

n s = n ipp HW n // . (6.9)
The initial temperature is set to T 0 = 0.1ς; i.e. the temperature is set to 0.1 for a unit variance model, and then scaled with ς computed from Algorithm 5.2. The annealing rate α T is set so that T ns = 10 -10 as:

α T = ( T ns T 0
) 1 ns (6.10) The above models are evaluated on a test split of the DOTA 0.5 dataset. Metrics are computed by matching the object proposals and Ground Truth using the IOU (Intersection Over Union), with a 0.25 threshold. We show detection results in Figure 6.4, with a score threshold set for each model to the one maximizing the F1 score. The Average Precision (AP) metrics (which are threshold independent), are shown in Table 6.3.

Results

Method

We notice in Figure 6.4, that although Ground Truth labels are noisy -which could have caused over-constraining problem with linear programming (see Section 5.2.1) -our models infer more regular configurations, that better fit the objects and their physical constraints.

Testing the models on noisy data in Figure 6.5 and on the ADS data in Figure 6.6 (we show highlights in Figure 6.7), shows that the added priors allow the PP based models to be more Lastly, CNN-PP ♦ and CNN-PP ⋆ show comparable performances (qualitatively and quantitatively), demonstrating the capability of the contrastive divergence method to infer the prior/interaction model (via the prior parameters), thus limiting the manual setting of parameters via trial and error.

Computational complexity

Our two MPP based methods -running on an Nvidia Quadro RTX 8000 Graphics Processing Unit -execute in 300s on average, for 16k parallel iterations (equivalent to 77k sequential iterations) for one image. With an efficient implementation, considering a density of 1.9 × 10 -3 (95 th percentile of the observed object densities2 ), we estimate the cost of one iteration to be around 5 operations per pixel per iteration, thus around 4 × 10 5 operations per pixel in total. We show the derivation of those values in Appendix H. As a point of comparison transformer models for image classification range from 5 × 10 6 to 1.8 × 10 7 operations (Y. Zhao et al., 2021). The complexity could be greatly reduced by reducing the number of iterations of the Point Process sampler, at the cost of straying away from the proper convergence properties.

Conclusion

In this chapter we applied our methods on remote sensing data sets; both on COWC [START_REF] Mundhenk | A Large Contextual Dataset for Classification, Detection and Counting of Cars with Deep Learning[END_REF] and DOTA [START_REF] Xia | DOTA: A Large-Scale Dataset for Object Detection in Aerial Images[END_REF] that are publically available, and on some data provided by Airbus Defense and Space. Even with the limited complexity of our models, we achieved good qualitative results by ensuring regularity in the inferred geometry. Quantitatively our models show better results than other CNN based methods. The latter is more pronounced with noisy input data: our models maintain good qualitative and quantitative results even in adverse noise conditions. We also demonstrate the proposed parameter estimation method consistently matches or outperforms slightly the manual setting of model parameters, qualitatively and quantitatively, with the intrinsic advantage of limiting tedious manual trial and error. Lastly we show our model can generalize to new data by looking at qualitative results on ADS data after training on DOTA.

CHAPTER 7

Conclusion and Perspectives

Conclusion

In this PhD thesis we address the incorporation of interaction models in object detection methods, while leveraging the capabilities of deep convolutional neural networks. Satellite images are inherently flawed due to the atmospheric perturbations, partial occlusions and limited spatial resolution. To compensate for this lack of visual information, it becomes essential to incorporate prior knowledge about the arrangement of objects of interest. On one hand CNN based methods are great at extracting patterns from images, but struggle to learn object to object interaction models without introducing attention mechanisms such as Transformers that drastically increase the complexity of the model: for instance, [START_REF] Lu | Few-Shot Object Detection in Aerial Imagery Guided by Text-Modal Knowledge[END_REF] proposed to include priors as text-model descriptors, while [START_REF] Zeng | Dynamic Cascade Query Selection for Oriented Object Detection[END_REF] made use of cascades of attention modules.

On the other hand, Point Processes propose to jointly solve for both the object likelihoods relative to the image (data term), and the coherence of the arrangement of objects itself (prior term). Firstly, Point Process approaches model the configurations as vectorial geometry, constraining the state space by construction. Secondly, these models allow formulating explicit priors on configurations as energy functions. However, in most of the literature [START_REF] Verdié | Detecting parametric objects in large scenes by Monte Carlo sampling[END_REF][START_REF] Schmidt | Forest point processes for the automatic extraction of networks in raster data[END_REF], the data terms relied on contrast measures between the foreground and background. We showed the limitations of those on satellite data in this manuscript. Although we should note (T. [START_REF] Li | Feature Extraction and Tracking of CNN Segmentations for Improved Road Detection from Satellite Imagery[END_REF] proposed transforming the original image with a CNN trained for classification, to then compute the contrast measures on that simplified result.

Instead of going the Transformers/complexity route, we proposed in this thesis to combine the CNN pattern extraction with the Point Process framework. The basic idea of this approach is to use the CNN output as the data term for a Point Process model. At first, we set up an original contrastlike method (pattern matching using cross correlation) on the output of CNN trained to map blobs at object centers. We then proposed a novel method using a CNN trained to infer a field of vectors that point towards the nearest object center. To obtain a heatmap we simply need to compute the divergence of this field. This allows to easily separate close objects at a low computational cost. Finally, we showed how we can reinterpret any heatmap-based CNN output as an energy that can be used for our model, both for the position and marks of each object. With our formulation of the data energy as the interpolation of a potential map, we transformed the contrast measure computation, into a simple lookup and interpolation of values in a 2D or 3D tensor.

CHAPTER 7 -Conclusion and Perspectives

With these pre-computed potential maps, we proceeded in adapting the Point Process sampling methods to make use of that easily available information. Indeed, these tensors define densities over a discrete space that is easy to sample from, and that approximates a truncated version of the energy model (which is defined in continuous space). First, we applied this data-guided sampling with the local transform perturbation kernel, and then built a birth map. These spatial densities also allow us to focus the parallel sampling of the Point Process to areas where the number of point is likely to be higher and thus requires more proposals. Finally, we leveraged the easily available automatic differentiation engines, to compute the energy gradient and perform diffusion for a more efficient exploration of the state space at a fixed number of points.

Previous methods would use Linear Programming to estimate the relative weights of the energies of the model (Q. [START_REF] Yu | Multiple-Target Tracking by Spatiotemporal Monte Carlo Markov Chain Data Association[END_REF][START_REF] Craciun | Joint Detection and Tracking of Moving Objects Using Spatio-temporal Marked Point Processes[END_REF] from a set of constraints. However, this is prone to over-constraining and only takes into account a limited amount of data points. Our first parameter estimation resulted in an SVM based method, that looks for the best possible separation between positive and negative samples. Thanks to its permeable boundary, this method allows to consider many more data points, even noisy ones that would contradict with hard constraints. Ultimately, we proposed adapting a method used for Energy Based Models: the contrastive divergence [START_REF] Hinton | Training Products of Experts by Minimizing Contrastive Divergence[END_REF][START_REF] Du | Implicit generation and modeling with energy based models[END_REF]. Within a gradient descent scheme, it alternates between generating contrastive/negative samples from the current energy, and maximizing the energy difference between positive samples (similar to the ground truth) and the contrastive samples. The procedure we presented allows estimating not only energy weights, but also the various internal parameters of the energy terms. Otherwise, those would have to be set manually (by trial and error) or with a separate procedure for each parameter of each energy term.

To be able to compare our model with other CNN based object detection models we introduced a scoring function for the objects inferred by the Point Process model. This score is computed from the Papangelou conditional intensity, which measures the likelihood of a point given the rest of the configuration. This is a key difference from classical scoring functions that only consider individual points, regardless of their neighborhood.

We apply our methods on remote sensing data sets; both on COWC [START_REF] Mundhenk | A Large Contextual Dataset for Classification, Detection and Counting of Cars with Deep Learning[END_REF]) and DOTA [START_REF] Xia | DOTA: A Large-Scale Dataset for Object Detection in Aerial Images[END_REF] that are publically available, and some data provided by Airbus Defense and Space. Even with the limited complexity of our models, we achieved good qualitative results by ensuring regularity in the inferred geometry. Quantitatively our models are better than some other recent CNN based methods. The improvement of results is more pronounced with noisy input data: our models maintained good qualitative and quantitative results even in adverse noise conditions. We also demonstrated the proposed parameter estimation method consistently matches or outperforms slightly the manual setting of model parameters, qualitatively and quantitatively, with the intrinsic advantage of limiting tedious manual trial and error procedures. Lastly we show our model can generalize to new data by looking at qualitative results on ADS data after training on DOTA.

Perspectives

Due to the iterative nature of MCMC sampling methods, and the limited optimization of our model (mainly to allow to modularity while experimenting), the inference time is significantly longer that state-of-the-art CNN based object detection methods. Speeding up the inference could be achieved by using fewer iterations, initializing the chain with an approximate inference (e.g. the CNN-LocalMax. output shown in Section 6.2.2) along with a faster simulated annealing. This would yield faster results at the cost of straying further from the optimal convergence conditions. Other solutions exist such as discretizing the state space and solving a binary variable optimization problem (T. T. [START_REF] Pham | Efficient Point Process Inference for Large-Scale Object Detection[END_REF].

While in this thesis we focused on small vehicle detection, the framework we proposed can be easily adapted to any other modality of objects where interactions between objects of interest are a key element. For instance road networks can be modeled with Point Processes of segments as in [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF], and have important priors on their interactions (e.g. limited bend radius). Point Processes can also be applied to temporal data; [START_REF] Craciun | Joint Detection and Tracking of Moving Objects Using Spatio-temporal Marked Point Processes[END_REF] performed tracking in remote sensing data. In this case, priors on dynamics are strong and could be able to complement the lack of visual information and even complete occlusions (e.g. car going under a bridge). Also, as our model is quite resilient to noise, we expect it would perform well on SAR data, as the input noise is more important with this type of data.

Future work could focus on the composition capabilities of such energy models: in short, it would be feasible to estimate the data energy functions on one set of data, and the prior energy model on another one, to then compose the two and form a complete model. In practice this would allow to learn a data model on imperfect annotations, and learn a prior model on perfect synthetic data that reproduces the expected structures; here we would only need configurations of objects, no difficult to produce images.

Finally, the proposed parameter estimation method is not bound to object detection: we noted earlier the generative nature of our Point Process model, thus we could apply it to learning patterns and structures of points in order to produce a generative model that will be capable of imitating patterns as done in [START_REF] Hurtut | Appearance-guided synthesis of element arrangements by example[END_REF]. The latter could benefit from non-application-specific prior model that would be learned on the data, such as we propose in Appendix G: using a universal approximator for the per-object prior, and learning aggregation operators over the interactions using attention mechanisms.

A Jacobian for the local transform kernel

In Section 3.3.2.3 we consider a Marked Point Process with dim(M) marks and denote d = dim(M) + dim(S) the dimension of each point.

For a move from y to x, with perturbation δ We denote y = {y 1 , . . . , y n } and x = {x 1 , . . . , x n }, with n = n(y) = n(x). We have

x i = y i , i ̸ = k and x k = y k + δ, also δ ′ = δ.
Here we illustrate a case where n ̸ = k for sake of readability, but the reader will easily picture how the result is similar if n = k. I d and O d respectively denote the identity matrix and zero matrix of size d × d.

First we have the following partial derivatives:

i = j ̸ = k i = j = k i ̸ = j ∂x i ∂y j I d I d O d ∂x i ∂δ O d I d - ∂δ ′ ∂y j O d O d - ∂δ ′ ∂δ I d
The Jacobian is then computed as:

⃓ ⃓ ⃓ ⃓ ∂(x, δ ′ ) ∂(y, δ) ⃓ ⃓ ⃓ ⃓ = ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ∂x 1 ∂y 1 . . .
. . . ∂δ ′ ∂y n ∂δ ′ ∂δ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ (A.1) = ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ I d . . .
⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ ⃓ (A.2) = det(I nd ) det(I d ) (A.3) = 1 (A.4)
B Training a CNN for small object detection in remote sensing data

In this section we present a few methods used to train the CNN models to build the data term regardless of the loss function.

Sampling on objects. When training on large satellite images such as in DOTA 0.5 , the density of objects can vary greatly. As we train the CNN with constant size patches (128, 128, 3), we sample those from the large images of the training dataset. However, sampling those patches uniformly produces many empty patches (no object of interest). Thus, we sample patches 1/3 of patches uniformly in space, and 2/3 centered on objects (with a random Gaussian perturbation to avoid consistently centered objects). This allows maintaining a good number of examples of objects, while still avoiding over-detection in empty areas. These patches are resampled in the training dataset every 8 epochs to avoid overfitting.

Augmentations. We perform image augmentations on the training patches to add synthetic diversity to the data. We implement two sets of augmentations: a realistic one, that produces images that perceptually fit the aspect of images we get from satellite imagery; and a strong one that produce images that could not realistically be produced by the sensor (if not for a technical issue). The strong augmented image remains legible though, for the vast majority (except some edge cases), the objects can be still be seen in the image. The augmentations are illustrated in The augmentations are implemented using Albumentations [START_REF] Buslaev | Albumentations: Fast and Flexible Image Augmentations[END_REF]. We find that training the CNN model on the strong augmentation gives better results on the validation dataset.

We believe this makes the model more robust thanks to the greater (artificial) variety of data seen at training.

Sampling hard patches.

As training progresses, we focus the training on some hard patches.

After some number of epochs, we compute the loss (or a part of it) on all the training images; for instance we compute BCE( ˆ H M , H M ) from (4.17) on every image (not patches). We then define a per pixel density from the error by normalizing it. Difficult patches are sampled using this density. We find that blurring the density map with a Gaussian blur allows sampling patches with a higher average error. These difficult patches are mixed in the previously described sampling procedure, with 1/3 difficult patches, and 2/3 from previous method. Inferring on large images. The Fully Convolutional aspect of the Unet allows it to be translation invariant. When processing a very large image, we cannot always load the full image into the model at once. We can cut the image into patches, infer on each patch, and stitch those back together. Given the translation invariant property we expect the same results as if we inferred on the whole image at once.

However, one has to consider the padding used within the model. When performing a convolution with a 3 × 3 kernel, the output will be 2 pixels smaller on both width and height if not padding is applied, as the convolution kernel reaches the edge. To maintain the same output size, we apply mirror padding to add 1 pixel on each side on each dimension to ensure the convolution output is of the same size. This, in turns, introduces some edge artifacts, as the mirror padding tries to mimic data that is not there, thus losing the translation invariance of the result. This padding is however useful when processing image, in order to be able to get results on the borders of the image.

As we can compute the size of this edge effect, we can take it into account when cutting the image into patches, so that we only keep the non-affected center part. We illustrate the cutting scheme in Figure B.2. The size of the edge effect s is approximated by:

2s = 2 d + 2 d+1 l + 8(2 d-1 -1) + 4, (B.5)
with d the depth of the Unet (number of pooling operations), and l the number of convolution in a bloc (l = 2 in the classical implementation of Unet we use). Knowing this edge size, we pad the image with this amount on each side, then split it into patches so that the inner part of the patches that have no edge effect connect together.

C Marks data energy derivation

This Section justifies the formulation of (4.30) in Section 4.2.3. First we consider a simple energy model with no interactions such that:

U (y, X) = ∑ y∈y V pos (y, X) + V κ (y, X) (C.6)
For any configuration of points y we have:

p(y|y, X) = p({y} ∪ y|X) p(y|X) = exp ( - ∑ y ′ ∈{y}∪y V pos (y ′ , X) + V κ (y ′ , X) ) exp ( - ∑ y ′ ∈y V pos (y ′ , X) + V κ (y ′ , X) ) = exp (-V pos (y, X) -V κ (y, X)) . (C.7)
Here p(y|y, X) gives the probability of y being a point, given points in the configuration y. Without interactions between points in the model, we have p(y|y, X) = p(y|X):

p(y|X) = 1 Z exp (-V pos (yX) -V κ (y, X)) p(y i , y j |X)p(y κ |y i , y j , X) = 1 Z exp (-V pos (y, X)) exp (-V κ (y, X)) (C.8)
From the position energy definition, we identify p(y i , y j |X) ∝ exp(-V pos (y, X)), thus:

p(y κ |y i , y j , X) = 1 Z ′ exp (-V κ (y, X)) . (C.9)
Meanwhile, the CNN classifier gives us:

p ˆ(y κ ∈ c κ |y i , y j , X) = Softmax( ˆ Z i,j κ ) c = exp ( ˆ Z i,j κ [c κ ] ) ∑ nκ c ′ =1 exp ( ˆ Z i,j κ [c ′ κ ] ) . (C.10)
We approximate probability p(y κ |y i , y j , X) over continuous values of y κ , as its discrete counterpart p ˆ(y κ ∈ ⌈c κ (y κ )⌉|y i , y j , X) of y κ belonging to the ⌈c κ (y κ )⌉ th bin of values (see (4.28)). Thus we have:

V κ (y, X) = -ln (p ˆ(y κ ∈ ⌈c κ (y κ )⌉|y i , y j , X)) (C.11) V κ (y, X) = -ˆ Z i,j κ [⌈c κ (y κ )⌉] + ln ( nκ ∑ c=1 exp ( ˆ Z i,j κ [c] ) ) (C.12)

D Fast computation of overlap

We find there is no simple closed-form formula for the area of intersection of two oriented rectangles: Intersec(u, y) = Area(u∩y) min{Area(u),Area(y)} . In order to be able to compute overlap values over numerous rectangles without needing to branch the program for each pair of objects, we use a faster approximation to replace the precise overlap value. This makes use of the hyperplane separation theorem applied to convex polygons (Polygon Collision -GPWiki, 2012): if two convex polygons are not intersecting there exists a line that passes between them; such a line exists only if one of the sides of one of the polygons forms such a line.

For a given pair of rectangles, we project both into one dimension four times; one for each long and short axis of each of the two rectangles. The 1D overlap can then be easily computed in each projection. The smallest intersection value approximates the intersection of the two objects. While is this not exactly the intersection value Intersec(u, y) defined in (4.37), this approximation returns 0 only when no overlap is present, and increases (continuously) to 1 when a full overlap happens. (i.e. a 2d max radius). Also, we built all perturbations to only translate objects up to δ max . Moreover, any complex perturbation can be simplified to a sequence of births and deaths (e.g. moving a point is one death and one birth). Thus, any perturbation is a combination of the following:

Intersec(u, y) ≃ max ⎧ ⎪ ⎨ ⎪ ⎩ 0, - 1 c max v∈{u,y} k∈{a,b} { Gap v,k (u, y) } ⎫ ⎪ ⎬ ⎪ ⎭ (D.13) c = min v∈{u,y} k∈{a,b} { min ( Proj v,k (u), Proj v,k (y)
• A death in cell c, given the Markovianity relation the acceptance rate α depends only on objects less than 2d max away from c.

• A birth in cell c or within distance δ max of cell c, thus the acceptance ratio if only function of points within 2d max + δ max . Note that other cells of the same mic-set as the current cell might add objects in their δ max -neighborhood.

We illustrate the above in

E.2 Acceptance ratio for a move in a cell

Here we aim to show that for a move from y to y ′ in cell c we have:

α(y, y ′ ) = α(C c ¯(y) , C c ¯(y ′ )
), (E.17)

This derives from:

∆U (y → y ′ ) = ∆U (C c ¯(y) → C c ¯(y ′ )
), (E.18)

with ∆U (y → y ′ ) = U (y ′ , X, θ) -U (y, X, θ). Any move in cell c can be decomposed into the removal of points from y in cell c, and the addition of points in the δ max -neighborhood of cell c.

Considering the series of configurations induced by this decomposition into N simple additions or removals we have:

∆U (y → y ′ ) = ∆U (y → y 1 . . . y N -1 → y ′ ) (E.19) = N -1 ∑ k=0 ∆U (y k → y k+1 ), (E.20)
with y = y 0 and y ′ = y N .

From here we only have to show (E.18) for a birth of a single point, as the energy change of a death is minus the energy of the opposed birth move. So we now need to show (E.18) for y ′ = y ∪ {u}. Then we have y = P 0 ∪ P 1 ∪ P 2 ∪ P 3 . We defined similarly P ′ 0 , . . . , P ′ 3 for y ′ . Here we suppose y ′ = y ∪ {u}. Since u is placed in cell c or δ max away from c we have P 1 = P ′ 1 , P 2 = P ′ 2 , P 3 = P ′ 3 , but P 0 ̸ = P ′ 0 .

We split the energy difference into our partitions: • The difference cancels out between P 3 and P ′ 3 as P 3 = P ′ 3 and ∀y ∈ P 3 , N y {y} = N y ′ {y} ⊆ P 3 ∪ P 2 by construction (i.e. the d max -neighborhood of points in P 3 in unchanged too).

∆U (y → y ′ ) = ∑ y∈Y V (y|N y ′ {y} ) -
• The same holds for P 2 and P ′ 2 , as P 2 = P ′ 2 and ∀y ∈ P 2 , N y {y} = N y ′ {y} ⊆ P 3 ∪ P 2 ∪ P 1 .

• For P 1 , P ′ 1 however we have P 1 = P ′ 1 but ∀y ∈ P 1 , N y {y} ⊆ P 2 ∪ P 1 ∪ P 0 and P 0 ̸ = P ′ 0 . Still we can write the following (also valid for P ′ Remark E.1 -The partition P 2 could be reduced further. Indeed, if we set P 2 = {y ∈ Y | ∃y ˜∈ c × M, d(y, y ˜) ≤ 2d max + δ max } \ P 1 , the above demonstration still hold, and we can compute the ∆U on a smaller partition of y. In practice, we use the larger P 2 , as points are stored in memory by chunks corresponding to their cell, making it easier to simply pull 9 cells from memory than checking for distances of each point against point u.

F Buffer: size and chain length Notations. One step of the contrastive divergence performs n s Markov chain iterations (with transition kernel K θ (•, •), resulting of parametrized energy model U (•, θ)) starting from a configuration x resulting in x ′ . We denote the latter by:

x ns×K θ ----→ x ′ (F.30)
We also denote y - n the negative sample at step n.

Basic contrastive divergence. Without buffer, at epoch n, the negative samples y -is generated as :

x 0 ns×K θn -----→ y - n , (F.31)

with x 0 = y GT [START_REF] Hinton | Training Products of Experts by Minimizing Contrastive Divergence[END_REF] or x 0 ∼ U(Y) [START_REF] Du | Implicit generation and modeling with energy based models[END_REF]. To increase chain length (and have negative samples of better quality), one has to increase n s , at a great computational cost. As we alternate between sampling negative configurations and upgrading the parameters, this greatly increases the inference time.

Introducing persistent chains. Here we consider a one-sized buffer (or persistent contrastive divergence [START_REF] Tieleman | Training restricted Boltzmann machines using approximations to the likelihood gradient[END_REF]) as described in Section 5.2.2.2. The resulting procedure for step n is as such: For this cross-epoch chain to converge, the updates from each θ n to θ n+1 has to be small enough so that the stationary density of the Markov chain at n + 1 is close enough to the density at n. Otherwise, the burn-in period for the chain to stabilize to its new stationary density might be longer than the number of Markov chain iterations n s . In practice θ n is updated via small gradient descent steps, and as the energy function U is Lipschitz-continuous w.r.t. to its parameters θ, it guarantees limited variations in the stationary density.

The size of the chain from an initial state x 0 to y - n depends on the random variable L representing the consecutive times the memory was picked in a row after a random initialization. The low of L is simply:

p(L = l) = (1 -p B )p l B .
(F.34)

In turn, the chain presented in (F.33) from x 0 to y - n , has an expected length given by:

E[Ln s ] = n s p B 1 -p B . (F.35)
In practice using p B = 0.95 yields an expected length of 19n s . This allows using a lower n s thus reducing the computational cost, while maintaining qualitative samples from the density derived from U (•, X, θ n ).

Increasing the buffer size, inevitably leads to a lower cross-epoch chain length as it introduces the possibility at each step to resume from a shorter cross-epoch length.

G Towards generic (learnable) energy terms

In Chapter 4 we proposed a series of priors with parameters to either set manually of estimate from the data. The selection and design of the energy prior functions to put in the total energy U is by itself a prior on the model, as it restricts the space of energies U (. . . , θ) generated by parameters θ. In this section, as a perspective, we propose to push further the parametrization of the energy priors in order to go beyond the function-design-induced constraints, while allowing to build energy function that better fit the data. We first go through a design of per-point priors with MLP, and then show that attention mechanisms can be used to learn a generic aggregation function for interaction priors.

Generic single-point priors. The above-defined priors on single points (e ∈ ξ prior-point ), are of the form V e (y). As the MLP is a universal approximator, a sufficiently generic prior function would be:

V sgl (y) = MLP θ sgl (y), (G.36) this single point prior could replace all single point priors defined previously (e.g. ( 4.33) or (4.34)) as the MLP can learn an approximation of those functions from the data.

Generic interaction priors. Similarly, we can define a generic form of prior for interaction energies: where MLP θ int (y, u) gives a potential per interaction, which are aggregated with the operator F int . Here, the choice of the aggregation operator F int remains. One could set if to a preselected operator min, max or average. But in the aim to generalize, we instead would need a family of operators parametrized by θ that can approximate these two operators, preferably within a continuous parameter space so that we can estimate the operators by gradient descent.

Attention as aggregation operator. In broad strokes, the operator F int needs to aggregate a set of point to point interactions into a single value by considering only the relevant elements. This reminds us of attention mechanisms in [START_REF] Bahdanau | Neural Machine Translation by Jointly Learning to Align and Translate[END_REF] [START_REF] Vaswani | Attention is all you need[END_REF], where the goal is to aggregate the query the relevant relations of a word in a sentence with its neighbors.

In that paragraph we will see that the attention mechanism fit our criterion for the design of an interaction energy prior. Attention is usually defined as:

Attention(Q, K, v) = Softmax ( QK ⊤ √ E ) v, (G.38)
with v a vector of N values K the corresponding keys of shape (N, E) and Q a query matrix of shape (L, E). The result is of size L, the number of query elements. In short, the keys K represent each element of v in an embedding space of size E. Each query vector Q get compared to each key in K in the matrix product QK ⊤ . After normalizing with a Softmax, this result serves to weight each element on v according to its dot product in the embedding space with the query elements. Since keys K, are supposed to represent each value in v, we can have them generated from the same input matrix V of shape (N, F ). Matrix V represents each of the N elements in a feature space of size F . Equation G.38 then becomes:

Attention(Q, K, v) = Attention(qQ w , VK w , Vv w ), (G.39) with Q w , K w , v w matrices of weights of respective sizes (F, E), (F, E) and (F, 1), and q a vector of L queries in feature space (shape (L, F )).

In our case, the elements to apply attention on are the neighbors N y {y} of y (N = |N y {y} |), queried by the current point y (L = 1). To transform each point into feature space if size F we propose using an MLP. Eventually we have:

V att (y, N y {y} ) = Attention(q y Q w , V y,y K w , V y,y v w ), (G.40) Similarly, for the aggregation operator we have: (G.44)

V
F att (v) = ∑ v∈v ω v v,
ω v = Softmax ( q y Q w (V y,y K w ) ⊤ √ E ) , (G.45)
where the construction of the Softmax operator ensure that have Remark G.1 -A careful reader might have noticed the shift from an unordered set representation to a vector representation in (G.41); it is crucial that the aggregation of interactions be permutation invariant as there is no inherent ordering of elements in the set of neighbors N y {y} . Hopefully the attention as defined above (i.e. without positional encodings) in permutation invariant [START_REF] Lee | Set Transformer: A Framework for Attention-based Permutation-Invariant Neural Networks[END_REF].

The greater number of parameters in this proposed energy model combining learned priors and attention makes the parameter estimation procedure harder to tune. We are still tuning the details of that approach, and hope to produce meaningful results soon.

Remark G.2 -From going through the proposed attention on neighbors approach, we notice striking similarities of the resulting model with the graph attention networks proposed in [START_REF] Veličković | Graph Attention Networks[END_REF]. Indeed, the energy computation over a configuration can be seen as a node level inference on a graph; each node of the graph is an object y in configuration y and the links are the neighborhood relationship. Within a message passing framework, [START_REF] Veličković | Graph Attention Networks[END_REF] (G.46) where -→ h i is the feature vector of node i, K the number of attention heads2 and α k i,j the weights resulting from the attention head k as illustrated in Figure G.6.

- → h ′ i = σ ⎛ ⎝ 1 K K ∑ k=1 ∑ j∈N i α k i,j W k - → h ′ j ⎞ ⎠ ,
Considering our current proposed energy model and the similarity it bears with graph models, future works could look into approaching the energy inference through the lens of Graph Neural Networks. From [START_REF] Veličković | Graph Attention Networks[END_REF].

H Algorithmic complexity

In this appendix we estimate the number of CPU cycles to perform an inference. Here we do not consider the parallelization of operations which can save a lot of computation time. Table 1: Cycles cost per operation, in number of cycles (approximate).

For F(Gap) we consider the procedure given in Appendix D: One has to compute the maximum over 4 Gap v,k (u, y) (v =∈ {u, y}, k =∈ {a, b}). Computing one Gap v,k (u, y) amounts to projecting the four corners of a rectangle into the axis of the other; i.e. performing a change of reference axis. Once the corners projected, one has to find the minimum and maximum of those four values. With the projection bounds found, it is only a matter to subtract values to find the gap. Considering the heading and orthogonal-heading vectors (long and short axis) of each shape precomputed, the cost amounts to: The pre-computation of the heading ( -→ y t ) and orthogonal-heading ( -→ y n ) vectors (for projection) necessitates computing cos(y α ) and sin(y α ), thus we approximate a cost of F( -→ y t , -→ y n ) ≃ 200 for one shape. In turn the total cost of Gap becomes:

F (Gap) ≃ 4F ( Gap v,k (u, y) ) + 3F(a, b) ⏞ ⏟⏟ ⏞ max Gap v,k = 73
(H.49)

The vectors -→ y t and -→ y n precomputed (to avoid computing those fore each pair of shapes), thus accounted for in a separate line in , V e (y, X) does not depend on u if y ̸ = u. To further simplify the computation we look into the interaction energies in the right-hand part of (H.71). First we specify the form of F e : With this we can approximate the complexity of computing ∂U (y) ∂u . We approximate F( ∂Ve(u,X) ∂u ) ≃ F(V e (u, X)) and F( ∂ve(y,y ′ ) ∂u ) ≃ F(v e (y, y ′ )): To get the operations per pixel (ops/px/iter) we compute F(K c )/d 2 c . In practice, we use n s = 77000 iterations of the Markov chain, so the total number of operations (ops/px) is n s F(K c )/d 2 c . In the testing data we observe a variety of object densities, we report complexity values for the minimum and maximum, average and 95 th percentile (q95) and report those values in Table 3 (a). For these density values, we compute the cost per pixel with diffusion on whole cells, on single object and without diffusion in Table 3

I Implementation: data structures and parallelism

We define all the energy terms in terms of tensor computation with an extra batch dimension so that we can compute energies over multiple cells at once, and leave the tensor computation library deal with the parallelization for us.

Example I.1 -Consider a simple sum operation: given a tensor a of size (N ), we compute the sum to return a single scalar

I.1 Batch computation of cells

When performing the RJMCMC in parallel, we have to compute the acceptance ratio for a set s ˜of B cells in parallel. We have the following property: We complement Y with a binary mask M Y to encode actual values against undefined, as in practice the ∅ contain some filler values. The mask is of size (B, N ). For our example:

M Y = ⎡ ⎢ ⎢ ⎢ ⎣ 1 1 1 0 0 1 1 1 1 0 1 1 0 0 0 1 1 1 1 1 ⎤ ⎥ ⎥ ⎥ ⎦
interpolate function interpolates values over the map Z in a way that the gradient can be back propagated to the configuration Y. This allows to perform diffusion easily (see Section 3.3.3).

I.2 Generic energy models

The above defined modules all join into the GenericEnergyModel module, that computes the total energy of a configuration. Here we present it stripped of the many debugging and logging functions our code actually contains in order to convey the main ideas of the implementation:

  = X|Φ = y)P(Φ = y), (0.0.2) où l'image X et la configuration y sont respectivement les réalisations de X et Φ. Ici P(X = X) correspond à la loi marginale des observations, construite à partir de la loi des observations P(X = X|Φ = y) et de la loi a priori P(Φ = y) sur le processus sous-jacent. Avec la formule de Bayes, nous pouvons inverser le conditionnement pour étudier la loi de Φ conditionnellement à l'observation X :P(Φ = y|X = X) ∝ P(X = X|Φ = y)P (Φ = y).(0.0.3)Pour minimiser le coût bayésien L(Φ * , Φ) = 1(Φ ̸ = Φ * ), l'estimateur optimal est le Maximum

  Construction du modèle d'énergie. Nous construisons le Processus Ponctuel de densité h à travers sa fonction d'énergie U , pour une configuration y et une image X : h(y) ∝ exp(-U (y, X)) y, X, N y {y} ) est l'énergie du point y, sachant l'image X et le voisinage N y {y} de y dans y. Cette énergie est une combinaison linéaire de termes d'attache aux données et de terme a priori.

Figure 1 . 1 :

 11 Figure 1.1: Nadar Élevant la Photographie à la Hauteur de l'Art (Honoré Daumier, 1862).

Figure 1 . 2 :

 12 Figure 1.2: Ground track of the Pleiades NEO-4 satellite orbited in August 2021. From (OrbTrack.Org, 2023).

Figure 1 . 3 :

 13 Figure 1.3: Optical pushbroom: the CCD linear detector acquires at a period ∆t. The interval between two rows of pixels ∆y is given by the time interval ∆t multiplied the satellite ground speed V S . Adapted from[START_REF] Tupin | Remote sensing imagery[END_REF].

  prior law P(Φ = y) on the underlying process. With Bayes formula we can revert the conditioning to study the law of Φ conditional to the observation X: P(Φ = y|X = X) ∝ P(X = X|Φ = y)P (Φ = y).

Figure 1 . 4 :

 14 Figure 1.4: Limitation of the Bayesian approach: (a) synthetic image; (b) Bayesian model; (c) detector model using contrast measures. From (Ben Hadj et al., 2010).

Figure 2 . 1 :

 21 Figure 2.1: A graph of the bibliographic references, papers are linked together through (most of) their references and to the themes (or keywords) of these papers. We highlight the main clusters in this graph.

Figure 2 . 2 :

 22 Figure 2.2: Example rectangle features shown relative to the enclosing detection window. The sum of the pixels which lie within the white rectangles are subtracted from the sum of pixels in the gray rectangles. From (Viola & Jones, 2001).

Figure 2 . 4 :

 24 Figure 2.4: Heatmap for center detection. Here the heatmap resolution is coarser than the image resolution, the local offset allows compensating for the discretization error. From (X. Zhou et al., 2019).

Figure 2 . 5 :

 25 Figure 2.5: Evolution of remote sensing datasets dedicated to machine learning tasks. The vertical axis represents data volume, while the circle size maps to the number of spatial pixels covered.From[START_REF] Schmitt | There is No Data Like More Data -Current Status of Machine Learning Datasets in Remote Sensing[END_REF].

Figure 2 . 6 :

 26 Figure 2.6: Diffusion model for image generation. The model learns the transition kernel p θ via the gradient of an energy within Langevin dynamics. From (Ho et al., 2020).

3. 3 . 2 . 4

 324 Parallel sampling . . . . . . . . . . . . . . . . . 3.3.3 Jump diffusion . . . . . . . . . . . . . . . . . . . . . . . . 3.3.4 Simulated annealing . . . . . . . . . . . . . . . . . . . . . 3.3.5 Stopping conditions . . . . . . . . . . . . . . . . . . . . . 3.4 CNN fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 1 :

 31 Figure 3.1: Poisson Point Processes on S = [0, 1] 2 : left is homogeneous with λ = 400; right has intensity measure ν(i, j) = λ(cos(2πi) + 1)/2.

  Definition 3.1.8. Let (S, d) and (M, d ′ ) complete, separable metric spaces, ν a finite Borel measure on S, η a probability distribution on the sigma-algebra of M, and f ν×η the density of a Poisson process on S × M.

Figure 3 . 2 :

 32 Figure 3.2: Realizations of the Strauss Point Process for several values of γ with β = 1. The underlying homogeneous Poisson process intensity over S = [0, 1] 2 is λ = 100. The radius of the circles represent half the interaction distance r.

Condition 3 . 1 . 1 .Figure 3 . 3 :

 31133 Figure 3.3: Realizations of the marked pairwise interaction process from (3.15) on S × M = [0, 1] 2 × [0.02, 0.2], with β = 1. The mark corresponds to the radius of each circle.

  .17) Secondly, the Monte Carlo Markov Chain sampling (see Section 3.3.2) requires the following condition: Condition 3.1.2. A Point Process derived from an unnormalized density h w.r.t. the measure µ of the Poisson process is locally stable if there exist M ∈ R such that h(y ∪ {y}) ≤ M h(y), ∀y ∈ Y, ∀u ∈ S.

n

  the interaction relation y ∼ n y ′ ⇔ ||yy ′ || < d max resulting in neighborhood N y {y} . We need to determine the relation ∼ m (with associated neighborhood ∂ y {y} ) for which the Point Process in Markovian. Given Definition 3.1.6, the dependency of ratio h(y ∪ {u})/h(y) will give us the minimal requirement for ∼ m (i.e. the most restrictive relation that still ensures Markovianity of the process):

  y∪{u} {y} (as ∼ n is symmetric). This shows the density ratio depends on u, and the first and second degree neighbors of u. Defining the relation for Markovianity such as y ∼ m y ′ ⇔ ∃y ′′ ∈ y, y ∼ n y ′′ ∼ n y ′ is impractical and depends on y; with ∼ n derived from d max , the following is sufficient to have Markovianity w.r.t. relation ∼ m : y ∼ m y ′ ⇔ ||yy ′ || < 2d max .

Figure 3 . 4 :

 34 Figure 3.4: Edge case for maximum number of neighbors. Current point is labeled y with its neighborhood represented as a large dotted circle. Other points in y are represented as small hollow circles. In this illustration n c,max = 4 for all cells c.

Figure 3 . 5 :

 35 Figure 3.5: Examples of local perturbations. "All" combines all transforms at once.

Figure 3 . 6 :

 36 Figure 3.6: Independence of cells. On the left the two cells are not independent as c 2 is not wide enough. On the right cells c 1 and c 3 are independent. From (Verdié & Lafarge, 2014) (Here ϵ corresponds to 2d max ).

Figure 3 . 7 :

 37 Figure 3.7: Unet architecture. From[START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF].

Contributions CHAPTER 4

 4 Building an Energy Based Model for object detection Dans ce chapitre, nous construisons le Processus Ponctuel à travers sa fonction d'énergie. Nous montrons d'abord les limites des mesures basées sur le contraste pour les images satellitaires, puis procédons à la construction d'un terme d'attache aux données à partir de la sortie d'un CNN. Nous proposons plusieurs manières de construire une carte de potentiel ; nous constatons qu'une approche à base de vecteurs pointant vers les centres, puis divergence donne les meilleurs résultats. Nous montrons également comment réinterpréter une sortie de pseudo-probabilité classique d'un CNN en une énergie pour le Processus Ponctuel. De plus, nous proposons un ensemble d'a priori pour régulariser les configurations de véhicules, tout en étant prudent sur le choix de l'opérateur d'agrégation des multiples interactions d'un objet. Enfin, nous discutons de la façon dont l'inférence du CNN est factorisée dans le pipeline de calcul d'énergie pour un coût de calcul limité.

Figure 4 . 1 .

 41 The marks are bounded, and the angles are defined modulo π such that M = [a min , a max ] × [b min , b max ] × [0, π]. Describing rectangles using sizes (y a y b ) and ratios (y a /y b ) is equivalent.

Figure 4 . 1 :

 41 Figure 4.1: Parametrization of points for a marked Point Process of rectangles.

Figure 4 . 2 :

 42 Figure 4.2: Realizations of the processes described in Example 4.1.1. The radius of each circle corresponds to the radius marks y r .

  interior

Figure 4 . 3 :

 43 Figure 4.3: Sampling areas to compute contrast for an ellipse. Parameters ρ int , ρ ext , δ int , δ ext drive the size and locations of these surfaces.

Figure 4 . 4 :

 44 Figure 4.4: Comparing contrast measures on synthetic (first line) and real data (second line). Positive and negative samples, represent respectively 1% and 99% of the samples for which the metric is computed. For precision-recall; blue: T-test; orange: image-gradient; red: CNN output; green: constant value. On the top right plot, blue and orange plots are overlapping.

Figure 4 . 5 :

 45 Figure 4.5: Example of energy map: (a) input image X from COWC (Mundhenk et al., 2016); (b) results of F passed through function q pos , red is high positive potential, blue is negative. We show a zoom on part of the image where nearby blobs clump together.

Figure 4 . 6 :

 46 Figure 4.6: Relaxation and discretization of the target map illustrated for a single object in 1D. For the discrete plot, values are scaled from range [0, n c ] to [0, 1] for ease of reading.

Figure 4 . 7 :

 47 Figure 4.7: Filter and response for energy potential

Figure 4

 4 Figure 4.8: (a) Closeup image with centers overlay; (b) centers heatmap H; (c) vector field H ∇ ; (d) divergence. For (b) and (d): red> 0, blue< 0.

Remark 4 . 2 . 5 -

 425 This vector inference approach can be adapted to produce score of object centers. Given tensor ˆ H div , we can produce an approximated probability p ˆ(y|X) = σ(a ˆ H div + b), where sigma is the sigmoid function and a, b two scalar parameters (e.g. from a linear block in the neural network). As the divergence operator is differentiable, the parameters a, b can be learned during the training of the whole model.

Figure 4 . 9 :

 49 Figure 4.9: Example of inferred energy map. Image from DOTA (Xia et al., 2018) with corresponding inferred energy map. Blue areas corresponds to negative potentials while red to positive potentials.

  .23) Meanwhile, if we consider the simplest Point Process model derived from U (y, X) = ∑ y∈y V pos (y, X), Definition 3.1.7 and Theorem 3.1.2 give us: p (N Φ (dy) = 1|Φ ∩ (dy) c = y ∩ (dy) c ) = h({u} ∪ y) over the whole pixel ρ, and relating it the probability estimate in (4.23), we propose to define V pos such as: exp(-V pos (y, X)) = σ

Figure 4 .Figure 4 .

 44 Figure 4.10: Illustration of function from (4.27) and its inflection point t pos .

Figure 4 .

 4 Figure 4.12: Energy over discrete marks (a) and continuous marks (b).

CHAPTER 4 -Figure 4 .

 44 Figure 4.13: Illustration of the priors.

Figure 4 .

 4 Figure 4.14: Distribution of area and ratio over the training data. (a) and (b) show distributions of areas and ratios over the data, and (c) their joint distribution. In (d) we propose a setup with two ratio priors plus two area priors, while (e) uses two joint priors (green is high energy, pink is low).

  {v align (y, u)} v align (y, u) =cos (|y αu α |t align ) ,(4.38)

Figure 4 .

 4 Figure 4.15: Combining attractive and repulsive priors.

Figure 4 .

 4 Figure 4.16: Illustration of triplet alignment prior.

  (a) Pass the image X through the CNN model to produce ˆ Z pos and ˆ Z κ . 2. For any configuration y ∈ Y relative to image X: (a) Compute data terms for each y ∈ y given ˆ Z pos and ˆ Z κ . (b) Compute each prior term for each y ∈ y; each depends at most on N y {y} . (c) Combine all terms into U (y, X) as in (4.2) (4.3).

CHAPTER 4 -

 4 v e is Lipschitz-continuous, the above is made discontinuous at d(u, y) = d max by the step function induced by the indicator 1. As such our energy model is Lipschitz-continuous for one object y everywhere in S × M except in set ∪ u∈y\{y} {x, x ∈ S, d(x, u) = d max }, i.e. the union of boundaries of each neighborhood, as shown in mauve in Figure 4.18 3 . 90 Building an Energy Based Model for object detection κ ∈ {a, b, α} e ∈ ξ prior -interact. e ∈ ξ prior -point

Figure 4 .

 4 Figure 4.17: Energy model pipeline. The pre-computed section can be reused for computing different y ∈ Y for a given image X. d(y, y ′ ) is the distance matrix.

Figure 4 .

 4 Figure 4.18: Energy continuity illustrated for a simple example: on the left, we move continuously one object through space and marks (trajectory from violet to green). The energy terms V e are plotted on the right. The total energy of a point (bottom plot on the right) is V (y) = V jntAR + V align +8V ovrlp -0.5V zrN br . In mauve, we represent the discontinuity surface in S, with t 1 , . . . , t 4 the points where the path of the sliding object intersects. The interaction radius d max is set at a low value for illustration purposes.

Figure 5 . 1 :

 51 Figure 5.1: Data-driven translation kernel applied to y. Perturbation δ is generated from the 2D density map D y,X . The reverse mode follows density D y ′ ,X which is centered on y ′ . The red and blue centered squares represent the integer positions ⌊y⌋ and ⌊y ′ ⌋ for y and y ′ .

Figure 5 . 2 :

 52 Figure 5.2: Mic sets and cells for a Point Process. Each point is represented in the color of the cell it belongs to. Cell c k s is the k th cell of mic-set s.We highlight cells of mic-set 0, for which we represent the radius in which each point can be moved δ max (colored dashed circle). Gray semitransparent circles represent each point interaction radius d max (neighborhood N y {y} ).

  extension) d(s) and d(c) the density d integrated over all cells of mic-set s and cell c respectively.

6 . 1 Figure 5 . 4 :

 6154 Figure 5.4: Effect of training the energy with contrastive samples generated from the current θ n .

Figure

  Figure adapted from[START_REF] Lecun | A Tutorial on Energy-Based Learning[END_REF].

1 Figure 5 . 5 :CHAPTER 5 -

 1555 Figure 5.5: Buffer size effect on contrastive samples. Buffers are represented in blue. In this example the Markov chain starts from y - •,0 ∼ U(Y), thus the unlimited buffer is only of size 3 at epoch n = 4. Configuration y - k,n corresponds to contrastive sample for image X k at epoch n, and K ns θn is the application of the transition kernel of the RJMCMC n s times with parameters θ n .

  U(B X,β ) with probability p B , else U(Y)

1

 1 Training example: Filtering out irrelevant terms

Figure 6 . 1 :

 61 Figure 6.1: Example of contrastive divergence for training. (a) and (b) show two training samples, a synthetic image with configuration of circles in lilac. (c) shows a configuration inferred from the model using the estimated parameters. In (d) we plot the loss and its components (see (5.29)) through the estimation procedure steps. (e) shows the weights value through the estimation procedure, as in (f) where we compare over multiple random initialization.

Figure 6 . 2 :

 62 Figure 6.2: Inference results on a test image (a to c) and weights during training (d to f) for each experimental setup.

Figure 6 . 3 :

 63 Figure 6.3: Detection on a sample of COWC 0.5 : (a) ground truth, (yellow : detected vehicle, orange: false negative) (b) estimated configuration y ˆ(green: true positive, red: false positive).Table shows metrics on this patch.

Table 6 . 3 :

 63 Average Precision (AP) values, respectively for default DOTA 0.5 and noisy DOTA N 0.5 data. Right : Precision Recall (PR) curves, color correspondence in left table.

Figure 6 . 4 :Figure 6 . 5 :

 6465 Figure 6.4: Samples of detection on the test dataset for various models. The score threshold (only used for display purposes to prune low score objects) is set to maximize the F 1 score for each model.

Figure 6 . 6 :

 66 Figure 6.6: Models applied to ADS data. [© Airbus Defense and Space]

Figure

  or fog channel shuffle or dropout of channels random brightness and contrast

  Figure B.1: Image augmentation examples

Figure B. 2 :

 2 Figure B.2: Cutting a large image for inference. Right, input image; left, patches overlay, plain lines are patches boundaries while dotted lines correspond to the center parch with no edge effect.

WhereFigure D. 3 :

 3 Figure D.3: Fast overlap computation illustrated: left has no overlap (gap>0); right has overlap (gap<0). In dashed orange the current projection axis, in gray the other three axes; the final result is the maximum over the four projection axes.

FigureFigure E. 4 :

 4 Figure E.4: Determining cell size from basic moves. A point is moved starting from a cell in the current mic-set. The point can be moved in the δ max -neighborhood. This move might change the energy of points in the (δ max + d max )-neighborhood. In short, the Markovianity ensures this move depends at most on the dashed area.

Figure E. 5 :¯

 5 Figure E.5: Space partitioning around cells

•0

  1 and C c ¯(y ′ )):∀y ∈ P 1 , N y {y} = N P 0 ∪P 1 ∪P 2The above equality is also valid for y ∈ P 0 and y ∈ P ′We are now left with what we needed to show (we highlight the changes from (E.26) in red):

  -1 pulled from the single element buffer. With this use of the buffer, we virtually run a Markov chain across the epochs/steps, for which the transition kernel is updated every n s steps:

  V int (y, N y {y} ) = F int { MLP θ int (y, u), u ∈ N y

  y,y [k] = MLP θ att (u k ), k = 1, . . . , N, {u 1 , . . . , u N } = N y {y} (G.41) q y = MLP θ att (y) (G.42)were Q w , K w , v w are parameters of the model (i.e. part of θ).The above definition fits within the generic formulation defined in(4.35) for the interaction potential v att :v att (y, u k ) = MLP θ att (u k )v w .(G.43)

∑ v∈v ω v = 1 ,

 1 thus ensuring condition (4.36) is met.

Figure G. 6 :

 6 Figure G.6: Left: the attention mechanism a(W -→ h i , W -→ h j ) employed by our model, parametrized by a weight vector -→ a ∈ R 2F ′ , applying a LeakyReLU activation. Right: an illustration of multi-head attention (with K = 3 heads) by node 1 on its neighborhood. Different arrow styles and colors denote independent attention computations. The aggregated features from each head are concatenated or averaged to obtain -→ h ′ 1 . From (Veličković et al., 2018).

F 1 N•

 1 e (v 1 , . . . , v N ) = N ∑ k=1 ω k (v 1 , . . . , v N )v k . (H.72)This form remains quite general as it encompasses the following aggregation operations (amongst many others):• min: with ω k (v 1 , . . . , v N ) = 1 k= arg min k ′ =1,...,N v k (k) • max: with ω k (v 1 , . . . , v N ) = 1 k= arg max k ′ =1,...,N v k (k) • Average: with ω k (v 1 , . . . , v N ) = Attention mechanisms as described in Appendix I.2Moreover we consider symmetrical interaction energies such that v e (y, u) = v e (u, y).The right-hand part of (H.71) can be then derived as: v e (y, y ′ 1 ), . . .

F≃

  |C c (y) |(399N y + 2480). (H.80) To compute the diffusion update, we sample |C c (y) | Gaussian distributed values. We grossly approximate this sampling to 100 (this is negligible compared to the cost with a N y |C c (y) | factor):F(Diffusion c ) ≃ 2580|C c (y) | + 399N y |C c (y) |. (H.81)Alternatively the diffusion applied to a single element y ∈ y costs:F(Diffusion y ) ≃ 2580 + 399N y . (H.82)H.3 Total costWe can now compute the cost of one Markov chain step K c in cell c. In practice, we have d max = 16 and δ max = 8 thus d c = 48 and n κ = 32. From these and the object density λ we can derive approximate values for the parameters in the equations above: c ) = 0.8F(Diffusion c ) + 0.1F(K B ) + 0.1F(K D ) ≃ 3.8 × 10 2 + 4.9 × 10 6 λ + 5.9 × 10 8 λ 2 . (H.87)

  (b).

  (a) Density values λ value λ min 1.4 × 10 -5 λ avg 7.9 × 10 -4 λ q95 1.9 × 10 -3 λ max 5.2 × 10 -3 (b) Cost per pixel

  i]. If we want to perform this operation over a batch of B tensors a 1 , . . . , a B , we concatenate those into A of shape (B, N ) such that A[j] = a j , now the sum operation over the batched tensor is such that the resulting tensor R of size(B) is R[j] = ∑ i A[j, i].The sum operation is parallelized across the batch dimension transparently, and results for each input element a j is stored in R[j].

  Denoting the point dimension D = dim(S) + dim(M), we can represent the current cells s content as a tensor Y of shape (B, N, D), where N is the maximum number of points in any cells and its neighbors:N = max c∈s ˜n(C c ¯(y)).Each element Y[k, l] is a D-dimensional vector, where each value represents (in that order) the two coordinates, width, length and angle of y c ¯k,l . For instance, considering a set of cells c 1 , . . . , c 4 with y c ¯1,1 , y c ¯1,2 , . . . the content of c ¯1 and similarly for the others. Denoting by ∅ the absence of point (since all c ¯do not have the same number of points necessarily), the following shows an example of Y for B = 4, N = 5, D = 5 (if we consider rectangles, each y is a vector of size D = 5):

  from torch import Tensor, nn from torch.nn import Module from torch.nn.functional import grid_sample c l a s s PositionEnergy(Module): d e f __init__( s e l f ,threshold: f l o a t ): s u p e r (PositionEnergy, s e l f ).__init__() s e l f .threshold = nn.Parameter(torch.tensor(threshold)) d e f q_fun( s e l f , Z_values:Tensor) -> Tensor: r e t u r n torch.log(1 + torch.exp( -Z_values + s e l f .threshold)) d e f forward( s e l f ,Y:Tensor,Z:Tensor, ** kwargs) -> Tensor: """ :param Y: batched configuration Tensor of shape (B,N,D) :param Z: energy map Tensor of shape (1,1,H,W) :return: position energy Tensor of shape (B,N) """ Z_values = interpolate( # interpolate map values at points positions positions = Y[:,:,:2].view(1, n_cells, n_points, 2), image = Z ) r e t u r n s e l f .q_fun(Z_values) # apply quality function d e f interpolate(positions:Tensor,image:Tensor) -> Tensor: """ :param positions: Tensor of coordinates (1,B,N,2) :param image: image Tensor (1,C,H,W) :return: tensor of values at positions (1,1,N,C) """ a s s e r t l e n (image.shape) == 4 a s s e r t l e n (positions.shape) == 4 h, w = image.shape[2:] samples = torch.stack([ # remap sample coordinates to range [-1,1] positions[..., 1] / (w -1), positions[..., 0] / (h -1) ], dim=-1) * 2 -1 r e t u r n grid_sample(image, samples, mode='bilinear', align_corners=True, padding_mode='border') Algorithm I.3: Computation of the position energy over batched cells. Interactions. For interaction priors such as the alignment prior (4.38), we compute energies over interactions as shown in Algorithm I.4. That means handling interaction energy tensors of shape (B, N, N ). The quadratic computation cost is kept in check as N is bounded by the maximum number of objects in a cell and its neighbors (9n c,max ). The tensor deltas represents the element-wise differences deltas[k, l, l ′ , d] = Y[k, l, d] -Y[k, l ′ , d]. Parameter inbound is a binary tensor encoding if two points Y[k, l], Y[k, l ′ ] are neighboring each other. Those are computed with the procedure shown in Algorithm I.5. Computing these interaction tensors outside each energy term module allows those to be re-used for each energy term.

from

  torch import Tensor, nn from torch.nn import Module import numpy as np d e f compute_distances_and_deltas(Y:Tensor, M:Tensor, dmax: f l o a t )-> Tensor: """ :param Y: batched configuration Tensor of shape (B,N,D) :param M: batched mask Tensor of shape (B,N) :param dmax: maximum interaction distance :return: deltas (B,N,N,D), distances (B,N,N), inbound (B,N,N) """ b = Y.shape[0] # =B n = Y.shape[1] # =N d = Y.shape[2] # =D # use torch boradcasting rules to form a difference tensor of shape (B,N,N, D) deltas = Y.reshape((b, n, 1, d)) -Y.reshape((b, 1, n, d)) distances = torch.sqrt( # compute distances from coordinates difference torch.sum(torch.square(deltas[..., :2]), dim=-1) + 1e-8 ) mask_1 = ~M.reshape((b, n, 1)) # mask out missing points in dim=1 mask_2 = ~M.reshape((b, 1, n)) # mask out missing points in dim=2 self_mask = ~torch.eye(n,dtype= b o o l ) # mask out distance to self dist_mask = mask_1 | mask_2 | self_mask # set missing points (or self) to be out of bounds distances = distances + dist_mask * (dmax + 1) inbound = distances < dmax # check who is inbound r e t u r n deltas, distances, inbound Algorithm I.5: Computation of the distances and component differences for points in batched cells.

  Proposed approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Thesis structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Point Processes for object detection . . . . . . . . . . . . . . . . . . . . . . . . .3 Energy Based Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Energy Based Models as generative models . . . . . . . . . . . . . . . . 2.3.2 Learning EBMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 EBM applications to computer vision . . . . . . . . . . . . . . . . . . . Point Process fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Stability conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Point Process for object detection . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Energy model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Sampling the Point Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Markov chains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Reversible Jump Monte Carlo Markov chain . . . . . . . . . . . . . . . 3.3.3 Jump diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.4 Simulated annealing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.5 Stopping conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 CNN fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Parameters estimation with Contrastive Divergence . . . . . . . . . . . . . . . . 6.1.1 Training example: Filtering out irrelevant terms . . . . . . . . . . . . . . 6.1.2 Stable perturbation kernel . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Results on optical data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.1 COWC with non-marked Point Process . . . . . . . . . . . . . . . . . . 6.2.2 Point Process of rectangles on DOTA and ADS data . . . . . . . . . . . 6.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Conclusion and Perspectives 7.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	xii	TABLE OF CONTENTS	CONTENTS	xiii
		Model foundations 3 Foundations for Point Processes and Convolutional Neural Networks 3.1 xi 6 Experimental results 6.1 Bibliography List of Figures List of Tables Notations & Acronyms Notations Point Process S Image space Sec. 3.1.1 M Mark space Def. 3.1.8 S d , M d Discretized image and mark spaces Eq. (5.6) y Point configuration Eq. (3.1) n(y) Number of points in y Def. 3.1.2 Y, Y n Set of all configurations with any number of points / n points Eq. (3.2) N lf Set of all locally finite configurations Sec. 3.1.1 µ, ν Density/intensity measure of the Poisson Point Process Def. 3.1.4 f, h Normalized / unnormalized density of the Point Process Eq. (3.5), (3.8) ∂ A B , N A B Neighborhood of B in A for a relation ∼ Def. 3.1.5 λ(y; y) Papangelou intensity Def. 3.1.7 Sampling π Stationary measure of the Markov chain Eq. (3.32) K transition kernel for the Markov chain Def. 3.3.3 Q, Q m perturbation kernel / sub-kernel for the RJMCMC Sec. 3.3.2 α(y, y Sec. 3.3.2 T t Temperature at step t Sec. 3.3.4 Energy model U (y, X) Energy function defining density h Sec. 3.2.1 V e Energy term Eq. (4.3) V Energy of a point: sum of energy terms Eq. (4.3) v e Interaction potential Eq. (4.35) 3.1.4 xii Aggregation function for interaction potential v e Eq. (4.35) F e ξ Set of energy terms Sec. 3.2.1

Introduction 1 Introduction 1.1 Computer vision in remote sensing . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Satellite image acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Object detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Image probability model . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Point Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Convolutional Neural Networks . . . . . . . . . . . . . . . . . . . . . . 1.4 2 State of the art 2.1 Convolutional Neural Networks for object detection . . . . . . . . . . . . . . . . 2.1.1 The rise of CNN models for object detection . . . . . . . . . . . . . . . 2.1.2 CNN for detection in remote sensing . . . . . . . . . . . . . . . . . . . . 2.1.3 Remote sensing datasets . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 2.2.1 Point Process models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Sampling the Point Process . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3 Point Process parameter estimation . . . . . . . . . . . . . . . . . . . . 23.1.1 Poisson Point Processes . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Markov Point Processes . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3 Markov marked Point Processes . . . . . . . . . . . . . . . . . . . . . . ′ ), r(y, y ′ ) Acceptance probability / Green ratio for move y → y ′

  Les travaux futurs pourraient se concentrer sur les capacités de composition de tels modèles énergétiques : en bref, il est possible d'estimer les fonctions énergétiques d'attache aux données sur un premier ensemble de données, et le modèle énergétique a priori sur un autre, pour ensuite composer les deux fonctions d'énergie et former un modèle complet (attache aux données et a priori). En pratique, cela permettrait d'apprendre un modèle de données sur des annotations imparfaites, et d'apprendre un modèle a priori sur des données synthétiques parfaites qui reproduisent les structures attendues ; pour ces données synthétiques nous n'aurions besoin que de configurations d'objets, pas d'images.

effectue un suivi des données de télédétection. Dans ce cas, les a priori sur la dynamique sont forts et pourraient être capables de compléter le manque d'information visuelle et même d'occlusions complètes (par exemple une voiture passant sous un pont). De plus, comme notre modèle est assez résistant au bruit, nous nous attendons à ce qu'il fonctionne bien sur les données SAR, car le bruit d'entrée (bruit de chatoiement) est plus important avec ce type de données. Enfin, la méthode d'estimation des paramètres proposée n'est pas nécessairement liée à la détection d'objets. Comme nous avons noté précédemment la nature générative de notre modèle

-CNN: Regions with CNN features Figure

		warped region	aeroplane? no. . . .
				person? yes. . . . tvmonitor? no.
	1. Input image	2. Extract region proposals (~2k)	3. Compute CNN features	4. Classify regions

CNN R2.3: Region with CNN features method overview. From (

  Stability conditions . . . . . . . . . . . . . . . . . . . . . . 3.2 Point Process for object detection . . . . . . . . . . . . . . . . . 3.2.1 Energy model . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1.1 Markovianity relation . . . . . . . . . . . . . . . 3.2.1.2 Local stability . . . . . . . . . . . . . . . . . . . 3.3 Sampling the Point Process . . . . . . . . . . . . . . . . . . . . . 3.3.1 Markov chains . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Reversible Jump Monte Carlo Markov chain . . . . . . . . 3.3.2.1 Detailed Balance . . . . . . . . . . . . . . . . . . 3.3.2.2 Balancing in terms of random numbers . . . . . . 3.3.2.3 Perturbation kernels . . . . . . . . . . . . . . . .

3.1 Point Process fundamentals . . . . . . . . . . . . . . . . . . . . . 3.1.1 Poisson Point Processes . . . . . . . . . . . . . . . . . . . 3.1.2 Markov Point Processes . . . . . . . . . . . . . . . . . . . 3.1.3 Markov marked Point Processes . . . . . . . . . . . . . . .

3.1.4

  4. There exists a symmetric measure ψ m (dy, dy ′ ), such that π(dy)Q m (y, dy ′ ) is absolutely continuous w.r.t. ψ m (dy, dy ′ ). The associated Radon-Nikodym derivative is denoted D m ;

	D m (y, y ′ ) =	π(dy)Q m (y, dy ′ ) ψ m (dy, dy ′ )	.	(3.40)
	3.3.2.1 Detailed Balance

  .53) We now have to find the symmetric measure on Y × Y. Note that both parts of the kernel only move between Y n to Y n+1 , thus we can restrict the study to subsets A n ⊆ Y n and B n+1 ⊆ Y n+1 . We can then write measures πQ(A n , B n+1 ) and πQ(B n+1 , A n ), with µ the density of the reference Poisson Point Process

	∫
	An

  unvisited cells will keep previous state */

	6: 7:	s do for all c ∈ ˜ if diffusion then
	8:

  ce chapitre, nous appliquons sur des données synthétiques et réelles les modèles et algorithmes proposés dans la thèse. Nous présentons d'abord l'application de la divergence contrastive sur un exemple simple utilisant des données synthétiques pour démontrer ses capacités. Nous présentons ensuite nos résultats avec des Processus Ponctuels non marqués sur un premier jeu de données : Cars Overhead With Context. Un modèle plus avancé est appliqué sur le jeu de données DOTA, qui montre d'excellents résultats par rapport à d'autres approches fondées sur les CNN ; les configurations inferrées sont plus régularisées et nos modèles sont plus résilients au bruit additif.

  Stable perturbation kernel . . . . . . . . . . . . . . . . . . 121 6.2 Results on optical data . . . . . . . . . . . . . . . . . . . . . . . . 121 6.2.1 COWC with non-marked Point Process . . . . . . . . . . . 121 6.2.1.1 Data . . . . . . . . . . . . . . . . . . . . . . . . 121 6.2.1.2 Point Process model . . . . . . . . . . . . . . . . 123 6.2.1.3 Parameter estimation and inference . . . . . . . . 123 6.2.1.4 Results . . . . . . . . . . . . . . . . . . . . . . . 123 6.2.2 Point Process of rectangles on DOTA and ADS data . . . . 124 6.2.2.1 Data . . . . . . . . . . . . . . . . . . . . . . . . 124 6.2.2.2 Models . . . . . . . . . . . . . . . . . . . . . . . 125 6.2.2.3 Parameter estimation and inference . . . . . . . . 126 6.2.2.4 Results . . . . . . . . . . . . . . . . . . . . . . . 126 6.2.2.5 Computational complexity . . . . . . . . . . . . 130 6.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

	6.1 Parameters estimation with Contrastive Divergence . . . . . . . 119
	6.1.1 Training example: Filtering out irrelevant terms . . . . . . . 119
	6.1.2

  An object of y GT is a False Negative (FN) if it corresponds to no point in y ˆ. A point of y ˆis a False Positive (FP) if there is no matching object in y GT . If multiple point in y ˆmatch the same element of y GT , we count one True Positive (TP) and every extra matching point as FP. Lastly if one detected element in y ˆmatches to multiple ground truth elements, we set the match to the closest ground truth object, if the remaining ground truth objects have no other matching detections they are counted as FN. From the TP, FP and FN we can compute the precision, recall and F1 score for an estimated configuration y ˆas shown in Table6.1.

	124			CHAPTER 6 -Experimental results
	Precision (Pr), Recall (Rc) and F1 scores for a given score threshold t are given as:
		Pr(t) =	TP(t) TP(t) + FP(t)		(6.6)
		Rc(t) = F1 (t) = 2 TP(t) |y GT | Pr(t)Rc(t) Pr(t) + Rc(t)		(6.7) (6.8)
	(a) ground truth		(b) estimated configuration
	nb. objects TP FP FN Precision Recall F1
	74	48 17 12	0.73	0.64	0.68

Table 6 .

 6 Table shows metrics on this patch. 1: Precision, recall and F1 on each test image of COWC 0.5 .

	Location No. obj. Prec. Rec. F1	File
	Potsdam	425	0.74 0.66 0.70 top_potsdam_6_9_RGB
	Potsdam	262	0.56 0.48 0.52 top_potsdam_6_8_RGB
	Potsdam	110	0.76 0.68 0.72 top_potsdam_6_7_RGB
	Selwyn	121	0.87 0.50 0.75 Selwyn_BX22_Tile_RIGHT_15cm_0003
	Toronto	3969	0.87 0.63 0.74 03747
	6.2.2 Point Process of rectangles on DOTA and ADS data
	6.2.2.1 Data

Table 6 . 2 :

 62 .2: Energy terms of the model.

		Energy term	Notation	Instances	Parameters	Eq. No.
	DT a Prior Terms	position mark joint area & ratio overlap alignment repulsive attractive	V pos V κ V jntAR V ovrlp V align V repls V attrc	-κ ∈ {a, b, α} k ∈ {car, truck} -t align ∈ {0, π 2 } --	t pos -µ k,ratio , µ k,area , σ k,ratio , σ k,area t ovrlp -t repls t attrc	(4.26) (4.30) (4.34) (4.37) (4.38) (4.39) (4.40)
		zero neighbor	V zrNbr	-	-	(4.41)
	a Data Terms

  O d . . . O d O d d . . . O d . . . I d O d O d . . . O d . . . O d I d

	. . .	. . .	. . .	. . .	. . .	. . .
	O d . . . I d . . . O d I d
	. . .	. . .	. . .	. . .	. . .	. . .

O

Table 2

 2 To pass from H.68 to H.69, we note that V e (y, X, N y {y} ) does not depend on u if y ̸ = u and y / ∈ N y {u} . From H.69 to H.70 we note that for e ∈ ξ loc.

	∂U (y) ∂u	=	∑ y∈y	∂V	(	{y} y, X, N y ∂u	)

  , v e (y, y ′ {y} = {y ′ 1 , . . . , y ′ Ny }, denoting the number of neighbors N y = |N y {y} | for simplicity. If we consider F e to be a min (or max), as our current model is done, the computation simplifies to:

	with N y
	∂V e	(	y, X, N y {y} ∂u	)	=	Ny ∑ k=1	1 k= arg min l=1,...,Ny	ve(y,y ′ l ) (k)	∂v e (y, y ′ k ) ∂u	(H.75)
									=	∂v e (y, y ′ ky ) ∂u	, with k y = arg min l=1,...,Ny	v e (y, y ′ l ).	(H.76)
	Equation (H.71) then becomes:
	∂U (y) ∂u	=	∑ e∈ξ loc.	w e	∂V e (u, X) ∂u	+	∑ e∈ξ inter.	∑ y∈{u}∪N y {u}	w e	∂v e (y, y ′ ky ) ∂u	(H.77)
												∂u	Ny ))v e (y, y ′ k )	(H.73)
					=	Ny ∑ k=1	∂ω k (. . . ) ∂u	v e (y, y ′ k ) +	∂v e (y, y ′ k ) ∂u	ω k (. . . )	(H.74)

Table 3 :

 3 Number of operations depending on object density.

	λ	ops/px/iter	ops/px
		With diffusion on c
	λ min λ avg λ q95 λ max	0.2 2.0 5.2 18.3 With diffusion on y 1.5 × 1.5 × 4.0 × 1.4 ×
	λ min λ avg λ q95 λ max	1.1 1.2 1.4 1.9 Without diffusion 8.2 × 9.2 × 1.1 × 1.5 ×
	λ min λ avg λ q95 λ max	0.8 1.0 1.4 2.3	6.3 × 8.0 × 1.0 × 1.8 ×

Previous approaches in NLP would build models that would run through each work of the sentence one by one, while updating a hidden state. This sequential approach limits parallelization and is limited by the memory that the hidden state represents(Vaswani et al.,

2017).

One can easily combine the densities of two EBMs by summing their energies.

tivates us to use CNN methods to transform the image into a representation that is easier to then factor into the energy function of the Point Process, in a way reducing this pixel-to-object gap.

For heavily imbalance datasets such a classifier can get decent results (e.g. if we invert the positive to negative proportions); we introduce it as a baseline check.

Denoting pos for position, as the model can include a potential for marks too, that will be discussed in a further Section.

Within a neighborhood (or intersection of neighborhoods) there are no discontinuities.

We use negative array indexing for simplicity, one can either offset all indexing by δmax for all indexing to be positive, or use the Python negative indexing convention A[-i] = A[len(A)i].

The heavy lifting of interpolation and differentiation is handled by the Pytorch library, here the grid_sample operation (Torch.Nn.Functional.Grid_sample -PyTorch 2.0 Documentation, n.d.)

We swap positive and negative labels due to the reverse ordering on energies; better configurations get lower energy. That way the equations fall back to the canonical formulations for SVM.

We shall add the condition that (y k ) k=1,...,n contains no more than one instance of each element of y + so that the set subtraction is properly defined.

Available at: https://gdo152.llnl.gov/cowc/

i.e. 95% of observed object densities are below this value.

Instead of one attention mechanism they use multiple attention heads at once and average the result.

Remerciements

CHAPTER 5 -Point Processes as Energy Based Models the bad characteristic (e.g. more overlaps). We obtain a new set of weights w e that would hopefully produce less of the negative characteristic. Section 5.2.2 proposes to address the Q -kernel design issues.

Margin maximization

With this first approach proposed in (Mabon et al., 2021), we try to find the best linear separation between ground truth and negative configurations/samples.

First we consider a set of negative samples y -generated -as with the previous method -from a kernel Q -; y -∼ Q -(y GT , •). We also introduce positive samples y + that are valid configurations, obtained by small perturbations of the ground truth with a perturbation kernel Q + : y + ∼ Q + (y GT , •). On one hand this models the uncertainty on the ground truth data due to its noisy nature. This also help to balance the separation problem introduced later on and artificially augments the dataset. Lastly, introducing energy term V 0 (y) = 1 and thus ξ ′ = {0} ∪ ξ we can rewrite (4.3) (5.22)

In a sense, we describe any configuration y as a point with coordinates v in the energy terms space; configuration energy U is a linear function from this space to R, defined as the dot product of a weight vector w and the coordinates vector v of y in this space. Given a set of negative and positive samples, we look for a hyperplane defined by w•v-b = 0, that best separates negative and positive samples. As shown in Figure 5.3, there might be more than one hyperplane (H 1 and H 2 in (a)) or no hyperplane that can separate all the data (H 3 in (b)). Equation ( 5.34) provides a pruning order y 1 , . . . , y |y ˆ| of points in y ˆ. This ordering allows to plot the precision and recall curve. Indeed, to trace a precision recall-curve, one only requires the sequence of (Recall(t s ), Precision(t s )) pairs, which are obtained by sequentially pruning the lowest scoring points. Equation (5.35) provides a score to each point y n .

Contrastive divergence loss and Papangelou intensity

On one hand we minimize the loss in (5.29) derived from the likelihood maximization, on the other we evaluate the performance of our detections with the scoring method in (5.35) sourced from the Papangelou intensity. In this part we show that while the two are derived differently, minimization of the loss function leads to good properties on the score function. Here we consider the simplified loss, as γ in (5.29) is small:

denoting ∆U (y → x) = U (x) -U (y). Theorem 3.1.2 gives us the following expression of the Papangelou intensity:

= exp(∆U (y ∪ {u} → y))

(5.37)

Case 1: single point addition. Considering a basic pair of positive contrastive samples, in which we add a non-valid point u to y + :

Now we have:

We get the expected behavior: minimizing the loss L leads to minimizing the score of point u.

Case 2: single point removal. Now considering the removal of a valid point y from y + to form y -: y -= y + \ {y}, u ∈ y + .

Then:

=log(λ(y; y + \ {y})).

Here we get the expected result that minimizing the loss lead to increasing the score of y.

List of Tables

List of Definitions

List of Algorithms Notations.

• n κ : number of classes for mark κ

We assume the distance matrix is precomputed, we only have to update parts of it when new points are added, or when points are updated.

H.1 Cost per energy term

We give the costs of each energy term (and the precomputing of -→ y t , -→ y n ) in Table 2. For interaction energies, we only specify the cost of computing on one interaction (not computing the reduction operator yet). With symbol ↑ we signify we use the results from the previous computation thus only counting it once. Interpolation (trilinear or bilinear), requires several lookups to interpolate (4 for bilinear, 8 for trilinear).

Term

Eq. We decompose the cost of computing V (y) into the cost of computing the local energies V e (y), the cost of computing the N y = |N y {y} | interactions v e (y, u), u ∈ N y {y} , and lastly the cost of aggregating those interactions with F e : The birth kernel application is split into tree phases : sampling a new point (K sample ), adding the point to the configuration (K add ), computing the Green ratio (r) to accept or reject the addition:

Sampling. To sample a new point, we pick one new point in the cell. The new point is picked in cell c (of size d c ) given a density map of size d 2 c , the marks are sampled with densities of size n κ each (3 marks). We approximate the cost of sampling an index in an array of size n to be around ln nF(a, b): we draw a random number with U([0, 1]) (complexity O( 1)) then look in the array representing the cumulative distribution (that can be seen as a sorted array) where the random number fits (O(ln n) comparisons).

Considering the cost to sample an array of size n to be n comparisons (3 cycles):

Adding the point to the configuration. Adding a new point requires computing its distance with the possible neighbors. Computing the distance itself is

We also compute the difference of each mark u κy κ ; the total cost of computing the distance and marks difference is F(dist, diff ) ≃ 25 When inserting a point we need to compute the distance with all points in c ¯, i.e. |C c ¯(y) | elements:

Computing the Green ratio. Computing the Green ratio for the addition of one point y requires computing the energy change induced by the addition of that point.

Sampling. The point to delete is sampled uniformly, we thus approximate:

H.2.3 Diffusion

To perform Diffusion on the configuration we have to compute ∂U (y) ∂u for u ∈ y. In our current code we compute directly ∂U (y) ∂u to facilitate implementation. However, we can simplify the whole computation as follows. Prior to that, we define two non-overlapping sets of energy terms such that ξ loc. ∪ ξ inter. = ξ:

• Local energy terms ξ loc. , such as V e , e ∈ ξ loc. can be written as a function of y, X and θ only; such as V pos , V κ , V zrNbr ...

• Interaction energy terms ξ inter. , such as V e , e ∈ ξ inter. can be written as a function of y, N y {y} , X and θ only; such as V align , V ovrlp ...

APPENDICES

Total energy. Denoting V(•) the equivalent in batched tensor operations of per point energy V (see (4.3)) we get:

)

[k],

with A ⊙ B the element wise multiplication of matrices A and B. In Python, we get Algorithm I.1. Point priors. For point-wise priors such as V ratio defined in (4.33), tensor operation V ratio is defined as a per-point operation (only acts on last dimension D of shape (B, N, D)). We show the resulting simplified code in Algorithm I.2. Data terms. Data terms are akin to point priors as they are computed for each point independently. We show an example for the position energy term defined in (4.26). Importantly, the APPENDICES Remark I.1 -Some optimizations remain to be implemented here: for instance the deltas tensor could be computed only on its upper triangular part as each deltas[k] is an antisymmetric matrix (A ⊤ = -A). d e f __init__( s e l f ,config): s u p e r (GenericEnergyModel, s e l f ).__init__() """ [...] """ s e l f .energy_modules: ModuleList = load_modules(config) # for instance self.energy_modules is a list of Modules such as PositionEnergy, AlignEnergy, RatioEnergy etc defined in the configuration file. s e l f .max_dist = s e l f .config['maximum_distance'] s e l f .energy_names = [p.name f o r p i n s e l f .energy_modules] s e l f .energy_combination_module = LinearCombinator( energy_names= s e l f .energy_names) d e f forward( s e l f , Y:Tensor, M:Tensor, Z:Tensor):

deltas, distances, inbound = compute_distances_and_deltas( Y=Y, M=M, dmax= s e l f .max_dist) energy_dict = {} f o r e i n s e l f .energy_modules:

energy_dict[e.name] = e.forward( Y=Y, M=M, Z=Z, deltas=deltas, inbound=inbound) energy_vector = torch.stack ([energy_dict[k] f o r k i n s e l f .energy_names ], dim=-1) energy_per_point = s e l f .energy_combination_module.forward(energy_vector ) * M #mask out the energy of non existing points energy_per_cell = torch.sum(energy_per_point, dim=-1) total_energy = torch.sum(energy_per_cell) r e t u r n {'energy_per_point': energy_per_point, 'energy_per_subset': energy_per_cell, 'total_energy': total_energy, ** energy_dict} c l a s s LinearCombinator(Module): d e f __init__( s e l f ,energy_names:List[ s t r ]): s e l f .f = nn.Linear(in_features= l e n (energy_names),out_features=1,bias= True) d e f forward(x:Tensor):

""" :param x: tensor (B,N,V) with V the number of energy terms :return: tensor (B,N), energy per point """ r e t u r n s e l f .f(x).squeeze(dim=-1) #squeeze to remove the last dim of size 1

Algorithm I.6: GenericEnergyModel, computes the energies of a batch of cells. APPENDICES 189

I.3 Miscellaneous methods

Divergence computation. The divergence (needed for computing position potential (4.22)) is computed over a field of 2D vectors defined in the image raster space. In Algorithm I.7 we show how it is computed using Pytorch: using such a library allows to leverage the automatric gradient computation.

from torch import Tensor d e f divergence(f: Tensor, indexing='ij'):

""" :param f: vector field compunents of shape (B,C,H,W) where C is the number of dims :param indexing: indexing type, "xy" or "ij" :return : Divergence, Tensor of shape (B,1,H,W) """ b = f.shape[0] n_dims = f.shape Apprentissage de modèles de géométrie stochastique et réseaux de neurones convolutifs. Application à la détection d'objets multiples dans des jeux de données aérospatiales.

Jules MABON

Résumé

Les drones et les satellites en orbite basse, dont les CubeSats, sont de plus en plus utilisés pour la surveillance, générant d'importantes masses de données à traiter. L'acquisition d'images satellitaires est sujette aux perturbations atmosphériques, aux occlusions et à une résolution limitée. Pour détecter de petits objets, l'information visuelle est limitée. Cependant, les objets d'intérêt (comme les petits véhicules) ne sont pas uniformément répartis dans l'image, présentant des configurations spécifiques. Ces dernières années, les Réseaux de Neurones Convolutifs (CNN) ont montré des compétences remarquables pour extraire des informations, en particulier les textures. Cependant, modéliser les interactions entre objets nécessite une complexité accrue. Les CNN considèrent généralement les interactions lors d'une étape de post-traitement. En revanche, les Processus Ponctuels permettent de modéliser la vraisemblance des points par rapport à l'image et leurs interactions simultanément. La plupart des modèles stochastiques utilisent des mesures de contraste pour la correspondance à l'image ; elles sont adaptées aux objets à contraste fort et faible complexité du fond. Cependant, les petits véhicules sur les images satellitaires présentent divers niveaux de contraste et une grande variété d'objets de fond et de fausses alarmes. Cette thèse de doctorat propose d'utiliser les CNN pour l'extraction d'informations, combinées aux Processus Ponctuels pour modéliser les interactions, en utilisant les sorties CNN comme données. De plus, nous introduisons une méthode unifiée pour estimer les paramètres du modèle de Processus Ponctuel. Nos résultats montrent l'efficacité de ce modèle sur plusieurs jeux de données de télédétection, avec régularisation géométrique et robustesse accrue pour un nombre limité de paramètres.

Abstract

Unmanned aerial vehicles and low-orbit satellites, including CubeSats, are increasingly used for wide-area surveillance, generating substantial data for processing. Satellite imagery acquisition is susceptible to atmospheric disruptions, occlusions, and limited resolution, resulting in limited visual data for small object detection. However, the objects of interest (e.g., small vehicles) are unevenly distributed in the image: there are some priors on the structure of the configurations. In recent years, convolutional neural network (CNN) models have excelled at extracting information from images, especially texture details. Yet, modeling object interactions requires a significant increase in model complexity and parameters. CNN models generally treat interaction as a post-processing step. In contrast, Point Processes aim to simultaneously model each point's likelihood in relation to the image (data term) and their interactions (prior term). Most Point Process models rely on contrast measures (foreground vs. background) for their data terms, which work well with clearly contrasted objects and minimal background clutter. However, small vehicles in satellite images exhibit varying contrast levels and a diverse range of background and false alarm objects. In this PhD thesis, we propose harnessing CNN models information extraction abilities in combination with Point Process interaction models, using CNN outputs as data terms. Additionally, we introduce a unified method for estimating Point Process model parameters. Our model demonstrates excellent performance on multiple remote sensing datasets, providing geometric regularization and enhanced noise robustness, all with a minimal parameter footprint.