<figure><img src="image_1.png" /></figure>

HAL
open science

# Structured Learning of Geospatial Data

Loic Landrieu

► To cite this version:

Loic Landrieu. Structured Learning of Geospatial Data. Machine Learning [stat.ML]. Paris-Est Sup, 2023. tel-04095452

HAL Id: tel-04095452

https://hal.science/tel-04095452

Submitted on 11 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

École des Ponts
ParisTech

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

IGN
INSTITUT NATIONAL
DE L'INFORMATION
GÉOGRAPHIQUE
ET FORESTIÈRE

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

Université
Gustave Eiffel

# HABILITATION À DIRIGER DES RECHERCHES DE L'UNIVERSITÉ GUSTAVE EIFFEL

**Spécialité:** Sciences et Technologies de l'Information Géographique

Loïc LANDRIEU

LIGM, Ecole des Ponts, CNRS, Univ. Gustave Eiffel
LASTIG, IGN/ENSG

## STRUCTURED LEARNING OF GEOSPATIAL DATA

Defended on May 9 2023 in front of a jury composed of:

**Reviewers :**

Pierre ALLIEZ

Nicolas COURTY

Maks OVSJANIKOV

Inria Sophia-Antipolis (TITANE)

Univ. Bretagne Sud (IRISA/Obelix)

Ecole Polytechnique (GeoVIC/LIX)

**Examinators :**

Camille COUPRIE

Christian HEIPKE

Konrad SCHINDLER

Martin BURGER

Carola-Bibiane SCHÖNLIEB

Meta (FAIR)

Univ. of Hannover (IPG)

ETH Zürich (IGP)

Univ. of Erlangen (Dep. of Mathematics)

Univ. of Cambridge (DAMTP, CIA)
<figure><img src="image_8.png" /></figure>
Acknowledgments

I can't thank my students enough for their hard work, enthusiasm, and creativity. It has been a privilege to guide such motivated future researchers through their first steps in academia.

I am thankful to the LASTIG lab, IGN, and the ISPRS for making me feel welcome in their community. In particular, I appreciated the trust and continuous support of Clément Mallet, Mathieu Brédif, Bruno Vallet, and Nicolas Paparoditis.

Lastly, I want to thank Alix for her world-class editing skills and for making me take a much-needed break once in a while.
# Preamble

I started my academic career with a Ph.D. in so-called “pure applied maths” [Loi16], working on functional optimization[RL15, LO16] and probabilistic graphical models [LO14]. While I immensely enjoyed the powerful tools I got to study and later develop, I also regretted not tackling concrete issues. I decided that after my Ph.D., I would start with real problems and let their characteristics and constraints guide my research. Circumstances and opportunities led me to study geospatial machine learning at IGN, and I could not have asked for a more exciting, impactful, and fertile playing field.

The message of this manuscript is simple: geospatial data has distinctive characteristics leading to unique machine-learning challenges and impactful applications. I draw two conclusions from this observation: (i) the analysis of geospatial data should rely on bespoke methods to leverage its unique structure; (ii) geospatial tasks are great candidates for evaluating and motivating new machine learning approaches. This habilitation showcases various methods involving abstract mathematical reasoning, sensor-specific considerations, and everything in between. The common motto is always to leverage the specificities of the considered problem into the architecture of the algorithm for added precision, speed, and parsimony.

In the title of this manuscript, I use the phrase “structured learning”, a reference to “structured optimization”. This sub-field of functional optimization consists in designing efficient algorithms that exploit the structure of the functional to minimize. Although most of the work here is more applied than pure mathematical optimization, I have applied this principle throughout all my research.

The term “geospatial data” refers to acquisitions with a large-scale spatial structure, without necessarily referencing their absolute position on Earth. This encompasses remote sensing data, such as satellite or aerial acquisitions, but also 3D scans or images taken from large scenes.
# Abstract

## Abstract

This manuscript presents an overview of my work in the field of geospatial machine learning, a rapidly growing interdisciplinary field that poses many methodological challenges and has a wide range of impactful applications. Throughout my research, I have focused on developing bespoke approaches that leverage the unique properties of geospatial data to create more efficient, precise, and parsimonious models.

This manuscript is divided into four main chapters, each covering a different property of geospatial data structures that can be leveraged algorithmically. The first chapter presents a versatile mathematical framework formalizing the concept of spatial regularity with graphs. We propose an efficient algorithm that tackles a broad family of spatial problems and provides novel convergence guarantees and significant speed-ups compared to generic approaches.

The second chapter introduces a deep learning method that extends the idea of exploiting graph regularity to the case of massive 3D point clouds. We simplify the task of large-scale semantic segmentation by formulating it as a small graph labelling problem. Our compact models reach high precision at a fraction of the computational cost of other approaches.

In the third chapter, we present a collection of methods designed to take advantage of the data structure inherited from 3D sensors. By considering the sensors' structure, we develop powerful networks with state-of-the-art accuracy, latency, and robustness for various applications and data types.

The last chapter dives into the real-life challenge of automated satellite time series analysis for crop mapping. Recognizing the difference between such data and standard formats used in computer vision, we propose novel and streamlined architectures that achieve unprecedented precision while remaining efficient and economical in memory and preprocessing. We also introduce the task of panoptic segmentation for satellite time series and an efficient architecture to solve this problem at scale.

In summary, this manuscript argues that geospatial problems represent a challenging and impactful venue for evaluating the newest machine learning and vision methods and a fertile source of inspiration for designing novel approaches.

An extended abstract in French is available at the end of the manuscript.
# Contents

<table><tr><td>Acknowledgments</td><td>i</td></tr><tr><td>Preamble</td><td>iii</td></tr><tr><td>Abstract</td><td>v</td></tr><tr><td>1 Introduction</td><td>3</td></tr><tr><td>1 The Unique Structure of Geospatial Data</td><td>3</td></tr><tr><td>2 Impactful Applications of Geospatial Analysis</td><td>7</td></tr><tr><td>3 Exploiting the Structure of Geospatial Data</td><td>8</td></tr><tr><td>4 Organization of the Manuscript</td><td>11</td></tr><tr><td>2 Exploiting Graph Regularity</td><td>13</td></tr><tr><td>1 The Cut Pursuit Algorithm</td><td>14</td></tr><tr><td>2 Cut Pursuit for Graph-Total Variation</td><td>18</td></tr><tr><td>3 l<sub>0</sub>-Cut Pursuit for Contour Regularizing</td><td>31</td></tr><tr><td>4 The Plane-Pursuit Algorithm</td><td>35</td></tr><tr><td>3 Exploiting the Spatial Regularity of 3D Data</td><td>39</td></tr><tr><td>1 Regularizing 3D Point Clouds Segmentation</td><td>40</td></tr><tr><td>2 The Superpoint Approach</td><td>43</td></tr><tr><td>3 Large-Scale Surface Reconstruction</td><td>54</td></tr><tr><td>4 Exploiting the Structure of 3D Sensors</td><td>59</td></tr><tr><td>1 Online LiDAR Segmentation</td><td>60</td></tr><tr><td>2 Image and LiDAR Fusion</td><td>65</td></tr><tr><td>3 Surface Reconstruction with Visibility Information</td><td>71</td></tr><tr><td>4 Forestry Analysis from Aerial LiDAR</td><td>78</td></tr><tr><td>5 Exploiting the Structure of Satellite Time Series</td><td>83</td></tr><tr><td>1 Stakes and Challenges of SITS Analysis</td><td>84</td></tr><tr><td>2 Parcel Classification</td><td>86</td></tr><tr><td>3 Panoptic Segmentation of SITS</td><td>94</td></tr><tr><td>4 Multi-Modal TAEs</td><td>101</td></tr><tr><td>5 Modelling Crop Rotations</td><td>107</td></tr><tr><td>6 Leveraging Class Hierarchies</td><td>109</td></tr><tr><td>6 Perspectives</td><td>115</td></tr><tr><td>1 Efficient Learning with Hierarchical Partitions</td><td>115</td></tr><tr><td>2 Cross-Modal Reciprocal Learning</td><td>116</td></tr><tr><td>3 Other Works</td><td>117</td></tr></table>
CONTENTS

1

7 Curriculum Vitae 119

Synthèse en Français 123

1 Contributions Scientifiques 123

2 Spécificités des Données Géospatiales 124

3 Structure de Régularité sur Graphe 125

4 Structure de Régularité des Données 3D 128

5 Structure des Capteurs de Données 3D 129

6 Structure des Séquences Temporelles Satellite 132

Bibliography 135
CHAPTER 1

# Introduction

Geospatial data analysis has obvious assets as a scientific field: virtually unlimited data, abundant annotations, challenging structure, unparalleled scale, and exciting applications. This combination makes it a promising and exciting field for evaluating and designing novel machine learning algorithms. This chapter details the specificities of geospatial data and presents some of its most motivating applications. We then argue for the benefits of using dedicated architectures, and present an overview of how the works in this manuscript leverage the unique structures of geospatial data.

# 1 The Unique Structure of Geospatial Data

Geospatial data refers to data augmented with their location on Earth. By extension, we use this term to designate data with a large-scale spatial extent without necessarily referring to an absolute positioning system. Such data are typically collected with various remote sensing sensors and exhibit particular properties not commonly encountered in natural images and videos. In particular, the diversity of complex sensors requires special considerations. Geospatial problems are typically large-scale, with entangled spatial, temporal, and spectral dimensions. Lastly, the absolute temporal and spatial frame of reference is a crucial characteristic that can be exploited for added precision and speed.

**Multi-Sensor.** As for photography, geospatial data inherits the structure of the sensors used for acquisition. While intrinsic and extrinsic parameters are often sufficient to characterize the acquisition of a natural image, remote sensing sensors require specific considerations. Geospatial data can be collected by various sensors: active or passive, mobile or static, terrestrial, aerial or spatial, within or beyond the visible light spectrum. In this work, we focus mainly on two types of data:

* *LiDAR Scans.* This active sensor uses lasers whose return times can be used to infer the location of reflective surfaces. LiDARs can be mounted on a fixed frame [HSL$^{+}$17, ASZ$^{+}$16], a mobile platform [GLSU13, LAL22], manned or unmanned aircrafts [KLMC22a, VAG20, IGN], and even satellite [DBG$^{+}$20]. The high precision and
<figure><img src="image_9.png" /><figcaption>Figure 1.1: Example of geospatial data types considered in this manuscript.</figcaption></figure>

low latency of LiDARs have contributed to their wide adoption for in-
dustrial applications. Spinning LiDAR, in particular, have a complex
acquisition geometry in which the sensor spins while being mounted
on a mobile platform in a dynamic environment.

* *Satellite Image Time Series (SITS)*. Satellites offer a privileged perspective on our planet in terms of extent and viewpoints. Some satellites provide optical information across a broad spectrum of frequencies thanks to sophisticated sensors [DDBC$^{+12}$], while others use radar waves to gather rich information about objects reflecting specific bandwidths [PRG$^{+16}$]. With carefully planned orbits, the satellites fly over certain areas at a periodical rate of a few days or less, producing temporal sequences of images.

**Large-Scale.** Earth Observation (EO) is one of the primary sources of geospatial data. Global coverage is standard for satellite [DDBC$^{+12}$, PRG$^{+16}$, DBG$^{+20}$], and country-scale acquisition campaigns are frequent for aerial imagery [IGN]. Geospatial data analysis tasks often require processing data at a near-global scale, e.g. biomass estimation for carbon capture assessment [LKA$^{+22}$] or agricultural monitoring for worldwide food security [MRB$^{+19}$, PI18].

Geospatial processes display influences at different spatial and temporal scales simultaneously: the crop cultivated in a field is influenced by the parcel's slopes, but also its bioclimatic regions; the fields' appearance dur-
ing an acquisition will be influenced by the position of the sun and the at-
mospheric conditions [RMM20], but also the cumulative weather patterns
of the last few months [NPA22] and local agricultural practices. Some pro-
cesses require considering large areas simultaneously: some urban land use
classes can only be decided based on their position within complex urban
configurations. The extent of EO is typically arbitrary and only depends on
the memory and computation available. Consequently, modelling large-
scale geospatial processes necessitates efficient and parsimonious methods.
In contrast, computer vision considers images with a fixed and predictable
extent and prioritizes precision over speed and scalability.

Geospatial processes also exhibit multiple and intricate data shifts. We
observe a combination of *prior-shift* (more vines in Italy than Norway),
*concept-shift* (roofs look different in the South or North of France), and
*covariate-shift* (meteorological conditions vary year to year) [KL18].

**Multidimensional.** The influence of the spatial dimensions on geospatial processes is evident; *"everything is related to everything else. But near things are more related than distant things"* is often referred to as the first law of geography [Tob70]. However, the temporal and spectral dimensions also play a crucial role in remote sensing data analysis.

The multi-temporality of geospatial data can be a consequence of the
orbit of satellites or the very nature of the sensor, such as spinning LiDARs.
For many applications, the temporal dimension of the acquisitions is cru-
cial for differentiating different classes, such as tree species or crop types
[GLGC19], which may appear similar at specific dates and only differs in
their evolution. For other applications, time is a structuring constraint,
e.g. real-time 3D perception for autonomous driving [LAL22].

Satellite images often have a greater spectral extent than natural images. This goes from a fourth *near-infrared* channel or the 13 bands of Sentinel-2, to hyperspectral imagery, which uses hundreds [MPH+22]. This information proves decisive when analyzing materials with specific spectral responses, such as vegetation [CRC05, TDH+13], minerals [CS95], water resources [GCB07], and anthropogenic structures [WHL08].

**Geo-referenced.** A key difference between natural images and geospatial data is the absolute spatio-temporal coordinates of the latter. One of video analysis' most challenging aspects is tracking the objects of interest across frames. The positions of static geospatial objects such as roads or agricultural parcels are given in an absolute referential (typically WGS84 [SM98], or Lambert93 [Tho52]), greatly simplifying the co-registration between acquisitions taken at different dates or with different sensors.

The spatio-temporal coordinates are also structuring information.
Knowing that a video frame is the 13-th of the sequence does not give much
information about its content nor impacts its interpretation. In contrast, the
time of year when a satellite image is taken drastically influences its anal-
ysis, and the exact time of day also impacts the radiometry through the
sun's position. Weather patterns also significantly affect geospatial pro-
cesses such as plant growth [NPA22], and dictate cloud cover, which com-
pletely blocks optical observations. Likewise, the exact coordinate of a pixel
in an image is usually not informative. In contrast, the absolute position
of a geospatial acquisition may strongly influence its interpretation due to
spatial domain shift.

**Richly Annotated.** Language models have access to a near-infinite corpus by training on the entire Internet with efficient self-supervised tasks [RWC+19, GBB+20]. Many aligned text-image pairs [RKH+21] are also available online as captioned images. However, direct annotations of text and images often require costly and time-consuming human intervention.

Geospatial data benefit from extensively available annotations. In-
deed, Geographical Information Systems (GIS) [Goo10, Cha08] are widely
used by public and private actors alike for numerous tasks such as ur-
ban and terrain modelling [For08], driving guidance systems [LWZ20],
or econometrics [BJL+02]. Thanks to increasingly open public data poli-
cies [geo], crowd-sourced initiatives [Map], and academic efforts [GL21b,
LAL22, GLSU13, SBK21], the amount of open-access annotated geospatial
data has starkly increased. For example, countries such as France and Den-
mark release yearly the cultivated crops in each agricultural parcel of the
country [RPG], amounting to millions of labels every year. Many govern-
ments publish publicly the cadastre (extent of properties) of their entire ter-
ritory [BDT]. **The wide availability of annotated geospatial data results
in a uniquely favourable setting for developing and evaluating machine
learning algorithms.**

Despite a surface-level resemblance, the structure of geospatial data
differs vastly from natural images and videos. In particular, the data
is gathered with a variety of active or passive sensors with unique
properties. The data also displays a complex multidimensional
structure with entangled spatial, temporal, and spectral dimensions.
The scale of data is also more extensive than typically encountered in
computer vision, requiring scalable algorithms for training and infer-
ence. However, remote sensing data benefits from absolute spatial
and temporal referencing, making tracking more accessible and con-
ditioning information extraction. Contrary to many machine learning
applications, remote sensing data can be paired with a large amount
of open-access annotations at a low cost.
## 2 Impactful Applications of Geospatial Analysis

After years of unbridled enthusiasm from the industry as well as the general public, some voices have called the true potential of computer vision for real-life applications to be over-hyped [Man18], unsustainable [TGLM21], and even unethical [VN20]. Such declarations can appear inflammatory and provocative. Notwithstanding, we argue that automated geospatial analysis assuredly leads to impactful and beneficial applications. Below is a non-exhaustive list of applications that motivated this manuscript's work.

**Crop Mapping.** The ability to classify and predict the yield of crops at a large scale is a crucial part of the growing concern for global food security [MRB+19, PI18]. Through its Common Agricultural Policy, Europe distributes over 57B euros each year in subsidies [CAP] (25B dollars in the US). The fair allocation of this assistance implies at least partial automation, as France alone counts close to 10 million individual parcels [RPG]. Moreover, crop type classification can also ensure that the best agricultural practices are employed to preserve soil [AAM+11, Bul92] and to maximize yield without relying on pesticides [KKN+15, SL08].

**Forest Inventory.** Automated vegetation analysis from aerial observation is an essential step for many forestry applications [BLBK17], such as forest management [JSM01], biomass estimation [FSM+16, Lu06], or forest fire modelling and prevention [MR11, Mac96, SOC01]. In turn, this allows for further ecological studies of the forest biome, such as estimating habitat suitability [HDW+06, MVG+09] or deriving biodiversity indicators [BD05, Nag01, IK98].

LiDAR-equipped satellite [DBG+20] offers the perspective of global monitoring of forest resources [LKA+22]. However, aerial or space-borne remote sensing of natural forests is often limited to canopy observation and misses its multilayer structure of [KLMC22a, FMJ+15]. Furthermore, annotations often require costly in-situ observations [KLMC22b].

**Land Use/Cover Mapping.** Thanks to aerial/spatial points of view, remote sensing allows states to oversee their territory's land usage extensively. This has proven especially important to monitor the increase in impervious surface [AJG96, SBT+05] and its harmful environmental effects. Automated and large-scale land-use classification is also helpful in monitoring urban sprawl [HL03, VOOR19] and is a step towards sustainable urban development [TL16].

Another use of automated remote sensing analysis is the possibility of emergency mapping [BGT15] to better facilitate disaster response and
damage assessment, for floods [ASM<sup>+</sup>22], earthquakes [DG12, DBC<sup>+</sup>11], tsunamis [KMMB20], or forest fires [CC89], and even conflicts [LS19].

**Digital Twins.** Extensive acquisition capabilities combined with efficient processing of 3D data allow for the perspective of digitizing entire cities [Bat18, KNL<sup>+</sup>20] as well as extensive industrial facilities [NSGW21, KFC<sup>+</sup>19]. Both prospects represent significant environmental and economic stakes. For example, digital cities allow to better design and monitor urban ecological corridors [PZL17]. Modelling of industrial infrastructures can increase productivity [VBR<sup>+</sup>17, MLL<sup>+</sup>19] or facilitate their decommisioning [PTB16].

**Autonomous Driving.** Autonomous driving is predicted to become a norm-shattering reality in the near future. Nevertheless, there are still several significant hurdles to overcome for large-scale deployment, not only regarding social, ethical, and engineering issues but also in terms of machine learning, computer vision, and remote sensing research. As the debate rages on over whether or not cameras are the only necessary sensor for fully autonomous vehicles, several major players have opted for a highly diversified array of sensors on their prototypes: thermal cameras, radar, LiDAR in various configurations, ultrasounds sensors, audio sensors, and so on [Zoo, GKM<sup>+</sup>20, Val, Lyf]. This diversity and complementary of sensors, combined with the entangled spatial and temporal dimensions of the acquisitions, make some aspects of this problem a flagship remote sensing / geospatial application.

**Cultural Heritage / Archaeology.** Remote sensing technologies have been used successfully for monitoring heritage sites and discovering hidden archaeological sites for several decades [LWG<sup>+</sup>19]. Aerial LiDAR has allowed archaeologists to find remnants of forgotten civilizations hidden under heavy forest cover in the Amazon [PBI<sup>+</sup>22] and the jungles of Angkor [EFP<sup>+</sup>13]. LiDAR scanning has also been used to save endangered artifacts [Ico, LL19] or help their protection [PBM<sup>+</sup>18], and more generally, backing up our shared cultural and ecological heritage [FLE<sup>+</sup>22].

## 3 Exploiting the Structure of Geospatial Data

Due to their unique properties, geospatial problems present unique challenges, making them a compelling opportunity for evaluating and designing new approaches in machine learning and computer vision.
**3.1 Geospatial Analysis as an Evaluation**

Machine learning methods are typically evaluated on problems taken from vision, language, or medicine. Geospatial tasks constitute an excellent standard for measuring models' performance, scalability, and expressiveness. We provide here some pointers on how to increase the appeal of geospatial problems as a standard evaluation task for learning algorithms.

**Accessible Datasets:** Accessibility and ease of use are central to the success of image datasets: images are easy to open and manipulate with most modern programming frameworks. This is far from the case for geospatial data, which tends to use dedicated formats which are typically harder to manipulate, e.g. .tiff for georeferenced images or .las for 3D scans. Using generic formats such as .npy or .h5 format can significantly lower the barrier of entry for researchers interested in applying their method to geospatial data. Datasets should also, of course, be freely distributed and accessible. Platforms like Zenodo offer free hosting with unlimited downloads for academic datasets.

**Cleaned and Focused Problems:** Large-scale spatial data exhibits a variety of complex domain shifts. For example, the composition of acquisitions at different dates, times, or days, translates into radiometric discrepancies. Rare classes can only appear in select regions, complicating the creation of meaningful train/validation/test sets. The class nomenclatures useful for practitioners can often be problematic for learning because of their ambiguity, excessively rare classes, or significant inter-class similarities.

When constructing datasets to assess the performance of learning algorithms, it is essential to curate the data to insulate the precise task evaluated. Indeed, spatial domain shifts or rare classes may bias the metrics in favour of specific architectures for reasons unrelated to the evaluated tasks. Adaptive sampling or reduced nomenclature can help mitigate these issues. Keeping adequate metadata can allow specific evaluation settings focused on the messiness and complexity of real-world data.

Geospatial benchmarks need to consider the spatial auto-correlation of spatial processes in their train/test split. There should be a clear demarcation and buffer between these parts of the dataset, as failing to do so may overtly favour large models overfitting the training set. However, unless the goal is to evaluate low-shot learning or domain adaptation, the data distribution of the training set should be relatively similar to the test sets. It is a tricky balance to find, and these questions should be attentively considered when designing geospatial benchmarks.
**Reproducible Code:** The geospatial community should systematize the release of open-access code for every publication. This will encourage authors to compare their work with recent baselines and allow for more thorough studies of the impact of a paper's contributions. The (almost) systematic release of research code played a significant part in computer vision's rapid evolution and success. In contrast, many geospatial articles still present methods evaluated on closed datasets without meaningful comparison to the state-of-the-art, open-access implementation or enough detail to re-implement them. Regardless of the brilliance of the ideas they explore, such papers have nearly no scientific value.

## 3.2 The Need for Dedicated Architectures

Remote sensing as a field shares the same goal as computer vision: analyze an object or environment from information gathered by a sensor situated at a distance. As detailed in Section 1, geospatial data and remote sensing differ from natural images in several key ways. However, a substantial part of geospatial analysis research directly applies computer vision methods to remote sensing data methods with little adaptation, limiting their efficiency [RCVS+19, TRW+21]. Developing algorithms and architectures that exploit the uniqueness of geospatial data leads to several benefits, which we list below.

**Adaptive.** Geospatial processes are typically large-scale and require networks with large receptive fields to capture. Likewise, as remote sensing sensors are typically further away than cameras, they have lower resolution and density, necessitating adapted processing. Computer vision focuses on analyzing objects or indoor scenes from dense acquisitions, which does not benefit from modelling long-range interactions to the same extent. For example, the nadir or near-nadir point of view of aerial scans differs from acquisitions at ground level. This results in distinct sampling density and a unique perspective on the world. Terrestrial LiDAR in urban areas captures building facades with centimetric precision, whose geometry is well characterized by grid-based approaches [CGS19] or dense convolutions. On the contrary, aerial scans have a density of only a handful (10-20) echos per square meter and are mostly limited to roofs, resulting in significant gaps between related entities. This motivates the use of networks that are robust to stark changes in density and able to connect distant entities.

**Efficient.** Computer vision favours large models with a high-level understanding of natural images. Remote sensing tasks are typically more grounded, less abstract, and require smaller models. Furthermore, the size of remote sensing acquisition favours algorithms that can scale to large in-
puts. Using gigantic pretrained computer vision models such as Founda-
tion Models [BHA+21] may give satisfying results out of the box, but incurs
superfluous computation and have unnecessarily high hardware require-
ments. Streamlining architectures to the bare essential results in faster and
more data-efficient approaches, at virtually no cost in precision.

**Multimodal.** While multimodality is explored in computer vision, it is an essential aspect of geospatial data frequently captured by many different sensors with complementary characteristics. Furthermore, georeferencing leads to a natural alignment across modalities, which is often difficult or impossible with natural images. While this alignment facilitates the fusion of various sources of information, this remains a nontrivial task due to their differing nature and resolution.

**Attractive.** The use of unaltered vision models for geospatial tasks may wrongly lead students and researchers outside of the community to believe that geospatial analysis is purely an applied field without methodological challenges. On the contrary, the successful use of novel and dedicated architectures for spatial problems promotes the field as a vector of innovation. This is a key endeavour to present geospatial tasks as a challenging and impactful source of evaluation and inspiration for novel methods.

# 4 Organization of the Manuscript

We present approaches to leveraging geospatial data's structure into more
efficient, parsimonious, and precise algorithms. The works are organized
into chapters covering different characteristics of geospatial data that can
be exploited.

**Chapter 2: Exploiting Graph Regularity.** We first present an abstract formalization of spatial regularity into a generic graph framework. We introduce the cut pursuit algorithm to minimize a large family of functionals commonly encountered in geospatial learning and other applications. We present a unified analysis for several versions of the method, including strong and unique theoretical guarantees and speed improvements of several orders of magnitude. This chapter contains this manuscript's most theoretical work, which will serve as the basis or inspiration for much of the work presented thereafter.

**Chapter 3: Exploiting Spatial Regularity of 3D Data.** This chapter presents a deep learning implementation of the ideas of the last chapter
for the particular case of large-scale 3D data analysis. By computing adaptive partitions of point clouds or 3D space, we transform complex problems of semantic segmentation or surface reconstruction into small-scale graph analysis problems. This approach allows us to reach state-of-the-art results with compact and efficient methods.

**Chapter 4: Exploiting Sensor Structure of 3D Data.** This chapter describes three distinct approaches that exploit the specific structure of 3D acquisition sensors. Our methods apply to LiDAR time sequences, joint image and point cloud scans, and large point clouds in the wild. Our approaches result in significant improvements in terms of precision, latency, and model size for semantic segmentation and surface reconstruction tasks.

**Chapter 5: Exploiting the Structure of Satellite Time Series.** This final chapter presents our work on satellite image time series (SITS) for automated crop mapping. By identifying the profound differences between SITS and videos, we design new dedicated approaches to exploit their specificities, which results in significant gains in precision and efficiency. We also present the first state-of-the-art for SITS panoptic segmentation and a unique dedicated dataset.
# Exploiting Graph Regularity

Many geospatial analysis problems can be formulated as the minimization of a functional defined on a graph and whose solutions are spatially regular. We first formalize this property as a form of graph-structured sparsity, then introduce the cut pursuit algorithm, which exploits this property for computational efficiency. Our approach can be adapted to functionals with various degrees of smoothness and continuity, covering a wide array of spatially structured problems such as inverse brain imaging or large-scale surface reconstruction. When applied to problems whose solution is spatially regular, the cut pursuit algorithm provides a considerable acceleration compared to other widely used optimization algorithms. In some settings, we also prove unique convergence guarantees that do not require smoothness or convexity of the minimized functional.

This chapter is organized around the following publications:

† [LO17] Loic Landrieu, Guillaume Obozinski, "Cut pursuit: Fast Algorithms to Learn Piecewise Constant Functions on General Weighted Graphs", *SIAM Journal of Imaging Science*, 2017

[RL18] Raguet Hugo, Loic Landrieu, "Cut-Pursuit Algorithm for Regularizing Nonsmooth Functionals with Graph Total Variation", *ICML*, 2018

[RL19] Raguet Hugo, Loic Landrieu, "Parallel Cut Pursuit For Minimization of the Graph Total Variation", *ICML Workshop on Graph Reasoning*, 2019

[GLCV19] Stéphane Guinard, Loic Landrieu, Laurent Caraffa, Bruno Vallet, "Piecewise-planar Approximation Of Large 3D Data As Graph-Structured Optimization", *ISPRS Annals of Photogrammetry, Remote Sensing & Spatial Information Sciences*, 2019

<footnote>† Work presented in the PhD thesis of the defendant [Loi16].</footnote>
# 1 The Cut Pursuit Algorithm

Graphs are commonly used to describe the structure of spatial data, such as point clouds or vector maps, by encoding the adjacency between individual elements. As we will exemplify in our numerical experiments, many real-life geospatial problems can be expressed as the minimization of a well-chosen functional whose variables are associated with the vertices of such a graph. We can express the spatial regularity of the minimizers as a form of graph-structured sparsity. We first propose formalizing these concepts mathematically and then introducing an efficient optimization algorithm capable of leveraging this graph sparsity.

We define an unoriented graph $G = (V, E, w)$ with edge weights $w \in \mathbb{R}_+^E$. We consider a functional $F: \Omega^V \to \mathbb{R}$ whose input variable $x$ can be decomposed with respect to the vertices of $G$: each $x_v$ belongs to the set $\Omega$. We seek $x^*$ a minimizer of $F$:<fnref n="1" />

$$x^* \text{ minimizes } F(x). \qquad (2.1)$$

We are interested in functionals $F$ whose minimizers $x^*$ are constant with respect to a partition of $V$ into a small number of components w.r.t $G$. This property corresponds to a regularity prior on spatial processes and can be exploited algorithmically for faster solving.

## 1.1 Graph-Structured Sparsity

We define two notions of regularity, which are related but not equivalent.

**Graph-Coarseness.** Mumford and Shah describe an image as *simple* if it is piecewise smooth, i.e. can be decomposed into a small number of regions with short contours and smooth variations [MS89]. This chapter focuses on a stricter form of simplicity: the signal is constant within each region. We can translate this concept for a signal $x$ defined on the vertices of a graph $G$: $x$ is simple if it is constant within the regions of a partition $V$ of $V$ with much fewer connected components than vertices, i.e. $|V| \ll |V|$. This property is referred to as *graph-coarseness* and is a widely used prior for natural images [ROF92] and medical imagery [EZE07] images. It is also commonly encountered when the vertices of $V$ correspond to elements positioned in space, such as pixels or 3D points. Since the concept of connectedness in a graph is inherently discontinuous, directly encouraging graph-coarseness makes the optimization problem significantly more complicated.

<footnote n="1">¹Depending on its properties, *minimizing* $F$ can mean finding a global or local optimum, a critical point, or simply trying to achieve low values of $F$.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 2.1: **Graph-Sparsity and Graph-Coarseness.** Example of two vertex-valued signals which are both graph-coarse: they only have two constant connected components characterized by white or black vertices. We represent the non-zeros values of the graph-gradient as red edges. The signal in (a) is also graph-sparse, while the signal in (b) is not.</figcaption></figure>

**Graph-Sparsity.** We qualify $x$ as graph-sparse if most vertices linked by an edge share the same value: $|\{(u, v) \in E | x_u \neq x_v\}| \ll |E|$. When $\Omega$ is equipped with a vector structure, this property translates into the sparsity of the graph-gradient, i.e. $\{x_u - x_v | (u, v) \in E\}$ is mostly zero. This corresponds to a form of structured sparsity as described by Bach et. al [BJMO12], and hints at the possibility of more parsimonious and efficient solving methods.

**Encouraging Coarseness.** As illustrated in Figure 2.1, graph-sparsity implies graph-coarseness. Indeed, if $G$ itself is connected, the number of constant connected components of $x$ is bounded by the number of non-zero values of its graph-gradient plus one. However, the converse is false: a graph-coarse signal is not necessarily graph-sparse.

Motivated by this observation, a straightforward way to impose a coarseness prior to a graph signal $x$ is to use a sparsity-inducing regularizer on the graph gradient. This observation leads to the following family of functionals, first characterized by Geman and Reynolds [GR92]:

$$F(x) := f(x) + \sum_{(u,v) \in E} w_{u,v}h(x_u, x_v), \qquad (2.2)$$

with $f: \Omega^V \to \mathbf{R}$ and $h: \Omega \to \mathbf{R}$ a function that reaches its minimum only when $x_u = x_v$. More precisely, we are interested in functions $h$ that encourage strict equality between $x_u$ and $x_v$.<fnref n="2" /> Consequently, we can expect

<footnote n="2">²When $\Omega$ is continuous and equipped with a vector structure, we typically have $h(x_u, x_v) = \rho(x_u - x_v)$ with $\rho$ a sparsity-inducing function, e.g. a non-smooth norm [BJMO11a].</footnote>
a minimizer $x^*$ of such functional $F$ to be graph-sparse and hence graph-
coarse. In particular, with $h(x_u, x_v) = \|x_u - x_v\|$ we obtain the well-studied
regularity-inducing graph total variation (TV).

The functional optimization literature explores many approaches to exploit the sparsity of the solution of some large-scale problems with fast and parsimonious algorithms [BJMO11b]. Analogously, we propose to exploit the graph-coarseness of functional minimizers of the form Equation 2.2 with a working set strategy. <fnref n="3" />

**1.2 Principle of Cut Pursuit**

The cut pursuit algorithm maintains a current partition $\mathcal{V}$ of $V$ and alternates between two steps: the *reduction step* and the *refinement step*, see Figure 2.2. The reduction step finds minimizers of $F$ that are also piecewise constant $w.r.t \ \mathcal{V}$. The refinement step splits $\mathcal{V}$ into a finer partition. Both steps can be performed efficiently by exploiting the structure of their corresponding optimization problems.

<figure><img src="image_11.png" /><figcaption>Figure 2.2: **Steps of the Cut Pursuit Algorithm.** At first, the vertices of the graph $G$ are combined into a single component, defining a single-vertex reduced graph $\mathcal{G}$, and the resulting reduced problem is solved (a). The partition is then iteratively refined (b,d) and the corresponding reduced functionals minimized (c,e). The reduction steps are performed $w.r.t$ the smaller reduced graph $\mathcal{G}$, while the refinement step involves the full graph $G$. The vertex colours represent their values, from blue to red.</figcaption></figure>

**Problem Reduction.** From a partition $\mathcal{V}$, we can define a *reduced graph*
$\mathcal{G} = (\mathcal{V}, \mathcal{E}, \mathcal{W})$ whose vertices are the components of $\mathcal{V}$ itself, the edges

<footnote n="3"><sup>3</sup>By analogy with the literature on sparse optimization, we distinguish between *working set* and *active set* algorithms. The former, also called *column generation algorithms*, maintains a set of constraints that can be active or not, while the latter considers the exact set of active constraints.</footnote>
denote the adjacent components, and W their weights:

$$ \mathcal{E} := \{(U, U') \in \mathcal{V}^2 | (U \times U') \cap E \neq \emptyset\} \quad (2.3) $$

$$ W_{U,U'} := \sum_{(u,v) \in U \times U'} w_{u,v}. \quad (2.4) $$

We can use this reduced graph to advantageously express the target functional $F$ for iterates $x$ that are piecewise constant $w.r.t \ \mathcal{V}$. We first notice that such $x$ can be written as $\mathbb{1}(r, \mathcal{V})$ with $r \in \Omega^{\mathcal{V}}$ such that:

$$ [\mathbb{1}(r, \mathcal{V})]_v = \begin{cases} r_U & \text{if } v \in U \\ 0 & \text{otherwise.} \end{cases} \quad (2.5) $$

In other words, a signal that is piecewise constant $w.r.t \ \mathcal{V}$ is completely characterized by only $|\mathcal{V}|$ values. In this case, we call $r$ the reduced variable, indicating the value of $x$ for each component of $\mathcal{V}$.

Without loss of generality, we can assume that $h(x, x) = 0$ in Eq. (2.2). In this case, $h(x_u, x_v)$ is zero for all edges $(u, v)$ linking vertices in the same constant component. Consequently, the regularization cancels out inside constant components. We introduce the reduced functional $F^{(\mathcal{V})} : \Omega^{\mathcal{V}} \to \mathbb{R}$, defined as the restriction of $F$ to variables that are piecewise constant $w.r.t \ \mathcal{V}$, and which can be written as such:

$$ F^{(\mathcal{V})}(r) := F(\mathbb{1}(r, \mathcal{V})) = f(\mathbb{1}(r, \mathcal{V})) + \sum_{(u,U') \in \mathcal{E}} W_{U,U'}h(r_U, r_{U'}). \quad (2.6) $$

We remark that $F^{(\mathcal{V})}$ has a very similar structure to the target functional $F$ and inherits its regularity. Furthermore, $F^{(\mathcal{V})}$ is defined $w.r.t$ the reduced graph $\mathcal{G}$, which is much smaller than $\mathcal{G}$ as long as $\mathcal{V}$ is coarse. Under this hypothesis, minimizing $F^{(\mathcal{V})}$ is significantly easier than minimizing $F$. If we had access to a constant partition of the sought solution $x^*$, solving Eq. (2.6) would also solve the full problem defined in Eq. (2.1), but at a drastically reduced cost. The next step consists in approximating this partition.

**Refining $\mathcal{V}$**. While the reduction step allows us to solve a simplified version of the problem, the resulting minimizer is only as good as the current partition $\mathcal{V}$. The goal of the refinement step is to add new degrees of freedom to Eq. (2.6) so that its minimization yields a lower value for $F$. In practice, we seek a refining partition $\mathcal{V}'$ that can be combined with $\mathcal{V}$ to obtain a new finer partition through cross-partitioning:

$$ \text{refine}(\mathcal{V}, \mathcal{V}') := \{U \cap U' | U \in \mathcal{V}, U' \in \mathcal{V}'\}. \quad (2.7) $$

In cases where $F$ is not convex, it can be useful to define a more complex refining operator refine($\mathcal{V}, \mathcal{V}'$) to mitigate the influence of suboptimal early decisions [SIBD11], see Section 3.3.
We consider a class of graph-functionals whose minimizers are piecewise constant with respect to a coarse partition and hence only defined by a few values. The cut pursuit algorithm exploits this property by iteratively refining a graph partition and solving a constrained version of the objective functionals. Both of these steps have specific properties allowing for their efficient computation.

**Algorithmic Structure.** Due to the structure of $F$ defined in Eq. (2.2), there exist minimizers $x^*$ of $F$ that are piecewise constant with respect to a partition with few components. Hence, we iteratively refine a partition $\mathcal{V}$ with a top-down approach until we find a partition that lets us compute a satisfactory minimizer of $F$. The different steps of the cut pursuit algorithm are represented in Figure 2.2, and can be summarized as follows:

* **Init.** Set $\mathcal{V}$ as the connected components of $V$, i.e. the coarsest possible partition of $V$.

* **Reduction.** Minimize the reduced functional $F^{(\mathcal{V})}$ with respect to the current partition $\mathcal{V}$:

$$r^{(\mathcal{V})} \leftarrow \min_{r \in \Omega^{\mathcal{V}}} F^{(\mathcal{V})}(r). \qquad (2.8)$$

The complexity of this step depends on the characteristics of $F$ and the coarseness of $\mathcal{V}$.

* **Refinement.** We now split the current partition $\mathcal{V}$ into finer components by computing a *refining partition* $\mathcal{V}'$:

$$\mathcal{V} \leftarrow \text{refine}(\mathcal{V}, \mathcal{V}') . \qquad (2.9)$$

The rationale for finding $\mathcal{V}'$ depends on $F$ and $\Omega$. In favourable cases, this amounts to finding a minimum cut in a well-chosen flow graph.

* **Stopping Criterion.** If the chosen stopping criterion is not reached, we return to the refinement step.

In the following sections, we explain how this approach can be implemented for different choices of $f, h$, and $\Omega$.

## 2 Cut Pursuit for Graph-Total Variation

With $h(x, y) = \|x - y\|$, finding the stationary points of $F$ as defined in Eq. (2.2) amounts to regularize $f$ with the graph total variation. This setting is well-studied and understood, and we are able to prove the convergence of our proposed scheme. We first provide a quick survey of the functional optimization literature on total variation. Then, we explain how to implement the cut pursuit algorithm for $\Omega = \mathbf{R}$ when $f$ is differentiable, and
to what extent these hypotheses can be relaxed. We then present heuristic approaches to deal with data from multivariate vertices, e.g. $\Omega = \mathbf{R}^n$, and to parallelize the graph cut-based refinement stage of cut pursuit. Finally, we provide several numerical experiments showcasing the considerable acceleration provided by the cut pursuit algorithm for problems with coarse solutions.

### 2.1 Background on Total Variation

Rudin, Osher and Fatemi [ROF92] first introduced the total variation regularizer as a tractable surrogate of the Mumford-Shah functional [MS89]. It has been at the center of sustained interest in the image processing community as a simple and versatile prior on natural and medical images [EZE07]. The total variation is also the subject of numerous works in the functional optimization community. Large-scale problems regularized by the total variation are typically solved by first-order proximal splitting algorithms [CP08, RFP13, CGN$^+$13] with or without pre-conditioning [PC11, RL15, MYWC18]. While these iterative methods enjoy convergence guarantees in the convex case, they require many iterations over thousands or even millions of variables. Luckily, most approaches can be parallelized in a straightforward fashion, which can mitigate this drawback [BS14, KBJ$^+$15]. Furthermore, a large array of functions can be regularized with such schemes [LP15].

An alternative approach exploits a deep connection with graph cuts: the total variation is the Lovász extension of the submodular cut function in a well-chosen flow graph [Bac11]. This link has led to several fast algorithms, among which the one of Chambole & Darbon [CD09] who reformulate finding the proximal of the graph total variation as a parametric max-flow problem. Such methods are typically very fast thanks to efficient graph cut solvers [BK04]. However, these approaches can only be applied to regularize functions $f$ with restrictive properties (convex, separable, smooth) and cannot be easily parallelized.

The cut pursuit algorithm for graph-total variation regularizing combines all the advantages of the previously mentioned methods:

(i) the speed of graph cut-based approaches,

(ii) the parsimony of sparsity-aware optimizers,

(iii) the ease of parallelization and

(iv) versatility of iterative proximal-based algorithms.
## 2.2 Regularizing Smooth Functionals

We first describe the cut pursuit algorithm in its most straightforward setting: $\Omega = \mathbf{R}$ and the function $f: \Omega^V \to \mathbf{R}$ is differentiable everywhere, but not necessarily convex.

**Reduction Step.** This step consists in retrieving a critical point $r^{(V)}$ of the reduced functional $F^{(V)}$ defined w.r.t a current partition $\mathcal{V}$, see Eq. (2.6). The difficulty in finding $r^{(V)}$ depends on the structure of $f$. If $f$ corresponds to an inverse problem, e.g. $f(x) = \|Hx - y\|^2$ with $H$ a linear operator, then the minimization of $F^{(V)}$ has the exact same form as $F$ only with $|\mathcal{V}|$ vertices instead of $|\mathcal{V}|$. For other forms of $f$ involving, for example, nonlinear operators across multiple coordinates of $x$, the reduced functional may not necessarily be easier to minimize than $F$. In this case, one should not use the cut pursuit algorithm.

The reduced functional can often be efficiently minimized with a first-order proximal-based algorithm, such as Generalized Forward Backward Splitting [RFP13]. Given the small number of variables, we can run many iterations and provide high-precision critical points. However, the reduced functional is typically poorly conditioned<fnref n="4" /> since the size of the components of $\mathcal{V}$ and the length of their interface can vary significantly. We can address this difficulty with a preconditioning strategy such as the one proposed by Raguet and Landrieu [RL15].

**Refinement Step.** The goal of the refinement step is to split the components of the current partition $\mathcal{V}$ such that the following reduction step can decrease the objective functional $F$ as much as possible. Intuitively, this can be achieved by refining the components of $\mathcal{V}$ to unlock new steep descent directions. The directional derivative of $F$ at the current iterate $x$ in direction $d \in \Omega^V$ can be written as follows:

$$F'(x, d) = \sum_{v \in V} \delta_v(x) d_v + \sum_{(u,v) \in E^{(x)}} w_{u,v} \|d_u - d_v\|, \quad \text{with} \quad (2.10)$$

$$\delta_v(x) := \nabla_v f(x) + \sum_{v \in V} \sum_{u \in \mathcal{N}(v)} w_e \operatorname{sign}(x_v - x_u), \quad (2.11)$$

where $\mathcal{N}(v) = \{u \in V | (u, v) \text{ or } (v, u) \in E\}$, $\operatorname{sign}(t) = -1, 0, \text{ or } 1$ depending on $t < 0, t = 0, \text{ or } t > 0$, and $E^{(x)} = \{(u, v) \in E | x_u = x_v\}$. We aim to find a steep descent direction $d$, i.e. leading to a small directional derivative $F'(x, d)$. The first term of Eq. (2.10) indicates the sensitivity of each vertex $v$ to an increase of their iterate $x_v$. This value depends on the gradient of

<footnote n="4">⁴Some directions are much more sensitive than others, which complicates finding a fitting gradient step size.</footnote>
<figure><img src="image_12.png" /><figcaption>Figure 2.3: **Illustration on an Image.** We construct a graph from the pixels' adjacency, and use the pixels' intensity to assign a value $y_v$ to each vertex $v$. We visualize the first two steps of the cut pursuit algorithm for a simple TV denoising problem: $\Omega = \mathbf{R}$, $f(x, y) = \|x - y\|^2$ and $h(x_u, x_v) = \|x_u - x_v\|$. We represent the graph cuts from the refinement stage as black (source-size) and white (sink-size) pixels separated by a red line representing the resulting minimal cut.</figcaption></figure>

$f$ and the values of adjacent vertices, see Eq. (2.11). The second term of Eq. (2.10) favours descent directions with identical amplitude for adjacent vertices with equal iterates. In other words, we must make a compromise between the individual tendency of each vertex to increase or decrease and avoid the creation of new discontinuities between adjacent vertices.

In general, finding the steepest direction at $x$ is a hard problem, potentially as hard as minimizing $F$ itself. However, we need not necessarily find the steepest direction but add relevant degrees of freedom to $F^{(V)}$. To make the problem more tractable, we restrict the space of potential directions to $D = \{-1, 1\}$: increase or decrease. Now, we can write the search for the steepest binary descent direction $d^{(x)}$ as a combinatorial optimization problem:

$$d^{(x)} \in \underset{d \in D^V}{\operatorname{argmin}} F'(x, d). \qquad (2.12)$$
<figure><img src="image_13.png" /><figcaption>Figure 2.4: **Flow Graphs.** Finding a minimal cut in the flow graphs above provides the steepest binary/ternary direction, and thus a refining partition. In this illustration, $x_u = x_v \neq x_z$. The *vertical* edges' capacity of the flow graphs are derived from the $\delta_v$ in the smooth setting, and from $\delta_v^+$, $\delta_v^-$ as well $m_u := \max(0, \delta_v^-(x), -\delta_v^+(x))$ in the nonsmooth setting. The *horizontal* edges are defined by the edge weights of $G$.</figcaption></figure>

Fortunately, this problem has a favourable structure and can be efficiently solved by finding a minimum cut in a well-chosen *flow graph*, see Figure 2.4a. Such a cut can be retrieved in quasi-linear time with efficient solvers such as the one of Boykov and Kolmogorov [BK04]. The refining partition $V'$ is then defined as the maximal constant connected components of the steepest binary descent direction $d^{(x)}$.

**Convergence Guarantees.** Despite restricting the search directions to $\{-1, 1\}^V$, the cut pursuit algorithm is guaranteed to find a stationary point of $F$ in a finite number of steps. If $x$ is not a stationary point of $F$, then $d^{(x)}$ is a *strict* descent direction. In other words, the following reduction step will decrease $F$. On the contrary, if $d_x$ does not refine $V$, then $x$ is provably a stationary point of $F$. In practice, since the partition induced by $x^*$ is coarse and the size of the current partition increases exponentially with each refinement step, we only need a few iterations to reach convergence at machine precision.

In [LO17], Guillaume Obozinski and I provided a first convergence proof by identifying the flows in the flow graph with the sub-gradient of
F. The idea of this proof is that if no non-trivial cut<fnref n="5" /> of the flow graph can be found, no descent direction of $f$ at $x$ is steep enough to justify creating a new discontinuity—more formally: $-\nabla f(x) \in \partial TV(x)$. We also provided additional motivation for our method by viewing the total variation as an atomic gauge [CRPW12] and the steepest binary direction as the Frank-Wolf direction of the regularization problem [Jag13]. In [RL18], Hugo Raguet and I developed a more elementary proof based on a systematic analysis of the directional derivative of $F$. The gist of the proof is that for a wide variety of functionals, if a strict descent direction exists, there also exists another strict descent direction whose vector representation is only composed of 1s and -1s. A major benefit of this last proof is that it does not assume the convexity of $F$.

The cut pursuit algorithm can efficiently recover critical points of a non-convex functional regularized by the total variation with no more than a few graph cuts. Furthermore, the graph gradient of our solutions is *exactly* sparse since all considered iterates are piecewise constant. This is in contrast to iterative algorithms, which are only sparse after thresholding procedures which often yield artifacts. Finally, our top-down strategy favours retrieving partitions with few components, offering a solution to the undesirable staircasing effects often encountered using iterative methods.

> The cut pursuit algorithm provably minimizes $F$ in just a few graph cuts. It provides convergence guarantees to machine precision, and its solutions are coarse without post-processing.

## 2.3 Extension to Non-Smooth Functionals

Proximal-based methods can be used to regularize a large array of functionals with the graph total variation [RFP13, LP15], including nonsmooth ones. In [RL18], we show that when $\Omega = \mathbf{R}$, the cut pursuit algorithm can be extended for $f$ consisting of a smooth part $s$ and a non-differentiable, vertex-separable part. In other words, $f$ can be written as follows:

$$f(x) = s(x) + \sum_{v \in V} g_v(x_v), \qquad (2.13)$$

with $s: \Omega^V \to \mathbf{R}$ differentiable everywhere and $g_v: \Omega \to ]-\infty, \infty]$ for all $v \in V$. We only assume that all $g_v$ are directionally differentiable in the sense that $\lim_{t\downarrow 0} \frac{1}{t}(g_v(x+td) - g_v(x))$ exists in $]-\infty, \infty]$ for all $x$ such that $g_v(x) < \infty$ and all $d \in \Omega^V$. Note that this definition is different from regular differentiability. In fact, any convex function is directionally differentiable regardless of smoothness [HUL04, part D], while the converse does not

<footnote n="5">⁵A cut is trivial if all vertices are on the same side: source or sink.</footnote>
hold. In practice, we choose $g_v$ as set indicators of convex subsets of $\Omega$, or $g_v = \| \cdot \|$ to obtain the LASSO regularizer [EHJT04].

**Reduction Step.** The reduction step is unchanged from the differentiable case as long as a sufficiently versatile solver is used to find a stationary point of the reduced function $F^{(V)}$, such as the Generalized Forward Backward [RFP13].

**Refinement Step.** The rationale of this step is also unchanged: we add degrees of freedom to the current partition by looking for a steep descent direction. However, since the functions $g_v$ are not differentiable, we must adapt Eq. (2.10). In particular, the directional derivative of $F$ at $x \in \Omega^V$ in direction $d \in \Omega^V$ now writes:

$$F'(x, d) = \sum_{\substack{v \in V \\ d_v > 0}} \delta_v^+(x) d_v + \sum_{\substack{v \in V \\ d_v < 0}} \delta_v^-(x) d_v + \sum_{(u,v) \in E^{(x)}} w_{u,v} \|d_u - d_v\|, \quad (2.14)$$

with

$$\delta_v^+(x) := \nabla_v f(x) + g_v'(x_v, +1) + \sum_{v \in V} \sum_{u \in \mathcal{N}(v)} w_e \operatorname{sign}(x_v - x_u) \quad (2.15)$$

$$\delta_v^-(x) := \nabla_v f(x) - g_v'(x_v, -1) + \sum_{v \in V} \sum_{u \in \mathcal{N}(v)} w_e \operatorname{sign}(x_v - x_u). \quad (2.16)$$

The main difference between Eq. (2.14) and Eq. (2.10) is that non-smooth functions $g_v$ may have different directional derivatives in the increasing (+1) or decreasing directions (-1). Consequently, it is now possible for a given vertex to not favour either the increasing or decreasing directions. This leads us to add the null value to the discrete set of possible directions $D$ when simplifying the search for a steep descent direction Eq. (2.12): $D = \{-1, 0, 1\}$. Solving this steepest ternary descent direction problem corresponds to finding a minimum cut in a suitable flow graph, represented in Figure 2.4b.

**Convergence Guarantees.** The non-smooth setting benefits from the same guarantees as the smooth setting: if $x$ is not a stationary point, then a strict ternary descent direction exists, i.e. $d^{(x)} \in D^V$ such that $F'(x, d^{(x)}) < 0$. Consequently, if solving Eq. (2.12) returns a trivial cut, then $x$ is a stationary point. In practice, we only need a few iterations (less than 10) to achieve machine precision for most standard total variation parameterizations.

## 2.4 Extension to Multidimensional Values

We can extend the cut pursuit algorithm to the case where each vertex is associated with a multidimensional variable. Although our convergence
guarantees do not hold, this approach efficiently leads to satisfactory approximate solutions.

**Choosing a Direction Set.** In the previous setting, restricting the set of candidate directions $D$ to increasing $(+1)$, decreasing $(-1)$, or null $(0)$ is enough to obtain our optimality certificate. This allows us to efficiently find the steepest binary/ternary descent directions with a single graph cut. In the setting where $\Omega$ is multidimensional, i.e. when the signal $x$ has several values for each vertex, such restricted directions may not be sufficient. To preserve optimality guarantees, $D$ would have to be large or even infinite, preventing us from efficiently finding a steep descent direction. However, we can choose $D$ heuristically to achieve empirical efficiency without a theoretical convergence guarantee.

A key point to note is that the set of potential directions $D$ can be chosen for each vertex independently and change according to the current iterate. If each vertex $v$ has its own set of candidate directions $D_v$, then the total set of directions in Eq. (2.12) is the Cartesian product of all directions $\times_{v \in V} D_v$.

**Examples.** In the case of $\Omega = \mathbf{R}^n$, we can define a heuristic set of directions for each component $U$ of $\mathcal{V}$ by computing $e_U$ the first eigenvector [Pea01] of the centered values of the iterate for the component: $\{x_v\}_{v \in U}$. We can then define $D_v = \{-e_U, +e_U\}$ for all $v \in U$, or add even more directions corresponding to other eigenvalues. Each component can adapt its descent directions according to its content, resulting in better splits and fewer computations than generic candidate directions.

We also propose a meaningful heuristic direction set when $\Omega = \Delta_K$, the $K$-dimensional simplex. For each vertex $v \in V$, we denote by $k_v := \operatorname{argmax}_{k \in [1,K]} x_{v,k}$ the corner of the simplex which is closest to $x_v$, i.e. the label for which $x_v$ associates the highest predicted probability. We then define $D_v = \{1_k - \mathbf{1}_{k_v} \mid k \in [1,K]\}$, with $[1_k]_l = 1$ if $k = l$ and 0 otherwise. In other words, each vertex can decide to transfer probability from a label to its most likely label or remain unchanged. Our split can factor in various configurations by adapting the descent direction at the vertex level.

## 2.5 Parallel Refinement

In practice, the refinement step is often the computational bottleneck of the cut-pursuit algorithm. The reduction step only involves $|\mathcal{V}|$ variables and relies on easily parallelizable solvers. In contrast, finding a steep descent direction implies finding minimal cuts in a flow graph with at least $V + 2$ vertices. Although the Ford-Fulkerson graph cut solver can be surprisingly fast [FF56], it is not easily parallelizable.
**Parallelizing Graph Cuts.** All binary terms in Eq. (2.10) and (2.14) correspond to edges within the constant components of the current partition $V$. There are no binary terms connecting vertices in different components. Consequently, $F'(x,d)$ is separable along the components of $V$. For the binary case, this means that we can rewrite $F'(x,d)$ as follows:

$$F'(x, d) = \sum_{U \in \mathcal{V}} F'_{U}(x_{U}, d_{U}), \quad \text{with for all } U \in \mathcal{V}, \qquad (2.17)$$

$$F'_{U}(x_{U}, d_{U}) = \sum_{v \in U} \delta_{v}(x) d_{v} + \sum_{(u,v) \in E^{(x)} \cap U^{2}} w_{u,v} \|d_{u} - d_{v}\|, \qquad (2.18)$$

In terms of graph cuts, no flow can be exchanged between the components of $V$, allowing us to compute the cuts in each component *independently*. We implemented a *multi-threaded* version of the already efficient implementation of the Ford-Fulkerson algorithm [FF56] by Boykov and Kolmogorov [BK04].

> Taking advantage of the particular structure of the refinement problem, we propose the first multi-threaded application of the Ford-Fulkerson algorithm.

**Balancing the Parallel Workload Distribution.** While our parallelizing strategy already considerably accelerates the cut pursuit algorithm, we can improve thread utilization even further. Indeed, when performing the refinement step in parallel, each component is assigned to a single thread. Consequently, the computation takes at least as long as the time required to refine the most challenging component, typically the largest one. If the partition is unbalanced, this can lead to poor threads usage. Consider the first iteration, where the partition only has one component, preventing parallelization.

For this reason, we advocate decomposing the largest components into smaller *balancing* components with a maximum size defined by the problem size and the number of available threads. This can be done efficiently with a breadth-first search from random starting vertices. In the refinement step, we set the border edge capacities between balancing components to 0 in Eq. (2.17). In the reduction step, we can ignore the balancing components. In theory, this partitioning strategy may lead to worse local minima or reduce the accuracy of the solution. In practice, the acceleration provided might offset the disadvantage of slightly less optimal solutions depending on the application.

Let us finally note that there is still room for improvement. For instance, we could replace the greedy breadth-first construction of the balancing components with better domain-specific heuristics. This is especially true for spatially embedded graphs as many fast and data-adaptive
partition algorithms exist for 3D point clouds [MD03, Mea82, YL22]. Furthermore, the compelling time and space complexity of the Ford-Fulkerson algorithm does not simply depend on the size of the flow graphs but also on the ratio between *vertical* and *horizontal* capacities: large horizontal capacities tend to incur longer computation time. This observation may lead to better parallel scheduling. However, an important limitation of this line of improvement is that the balancing step must be efficient or risk negating any computational benefit otherwise provided.

<figure><img src="image_14.png" /><figcaption>**Figure 2.5: Qualitative Results.** Illustration of the solutions given by the cut-pursuit algorithms. (a) noisy and regularized point cloud reflectance (detail of a much larger scene); (b) noisy prediction and regularized point cloud classification. (c) Segmentation Regularization</figcaption></figure>

## 2.6 Numerical Experiments

We showcase the performance of the cut pursuit algorithm across three different experiments. See Figure 2.5 for a visual representation of these tasks.

**Point Cloud Reflectance Regularization.** We spatially regularize the reflectance of 3D points acquired with a mobile LiDAR scanning vehicle as described by [PPC+12]. This corresponds to a denoising problem $f: x \rightarrow x$
$\|x - y\|^2$ with $y$ the observed noisy reflectance. $G$ is the 5-nearest neighbours graph of the point cloud ($|V| = 109$m and $|E| = 317$m). We expect the solution to be simple as anthropic scenes typically contain parts with uniform materials and reflectance.

**Brain Source Identification in Electroencephalography.** We consider the inverse problem of brain source identification from electroencephalography. The brain of a patient is mapped to a triangular mesh whose adjacency structure is described by a graph $G = (V, E)$ with $|V| = 20$k and $|E| = 29$k. A set of $N = 91$ electrodes records the patient's brain activity $y \in \mathbb{R}^N$, and the goal is to retrieve the neuron activation on the detailed mesh. The relationship between the electrodes' output and the brain activity is given by a *lead-field* linear operator $\phi: \mathbb{R}^V \to \mathbb{R}^N$, which is derived from physical and physiological considerations by experts of the domain [BAC$^{+}$14]. To model the regularity, sparsity, and positivity of brain signals, we set

$$f: x \mapsto \frac{1}{2} \|y - \phi x\|^2 + \sum_{v \in V} (\lambda_v |x_v| + \iota_{\mathbb{R}_+} (x_v)), \quad (2.19)$$

with $\lambda_v$ the parameters of the weighted LASSO regularization and $\iota_{\mathbb{R}_+}$ the set indicator function of $\mathbb{R}_+$.

**Point Cloud Classification Regularization.** We consider the problem of spatially regularizing noisy semantic labellings of point clouds among a class set $K$ [LRV$^{+}$17]. $G$ is the 10-nearest-neighbours graph of the point cloud, ($|V| = 3000111$ and $|E| = 17206938$). A noisy classification $y \in \mathbb{R}^{V \times K}$ is given by a random forest classifier operating on handcrafted vertex geometric features as described in [GL17]. Noting $\iota_{\Delta_K}$ the convex indicator of the $|K|$-dimensional simplex, and KL $(r,s) = \sum_{k \in K} r_k \log(r_k/s_k)$ the *Kullback-Leibler divergence*, we choose

$$f: x \mapsto \text{KL}(y_v, x_v) + \sum_{v \in V} \iota_{\Delta_K}(x_v). \quad (2.20)$$

In this multidimensional setting, we use the heuristic direction set described in Section 2.4 for simplices.

**Competing Methods.** We compare the performance of the cut pursuit algorithm against highly specialized graph cut-based approaches and flexible proximal algorithms. The competing methods—all implemented in C++—are as follows:

* **PMF:** the parametric max flow-based algorithm of [CD09] for the proximity operator of the graph total variation;
* **PFDR:** the Forward-Douglas-Rachford splitting algorithm of [Rag18] preconditioned according to [RL15], with all proximal steps parallelized;

* **CP:** cut pursuit algorithm with PFDR to solve the reduced problem;

* **PCP:** CP with parallel refinement, without balancing the parallel workload distribution;

* **PCP-balanced:** PCP with component balancing.

**Analysis.** We report the performance across methods and tasks in Figure 2.6. We first compute a high-precision estimate of the optimal value for these convex problems by running PFDR for $1e5$ iterations, which leads to an iterate evolution under $10^{-8}$ as measured by the largest coordinate change. This allows us to approximate and plot the distance towards optimality.

> When the sought solution is coarse, the cut pursuit algorithm can accelerate by several orders of magnitude compared to state-of-the-art solvers. In simple settings such as TV-denoising, our approach outperforms the highly specialized PMF algorithm [CD09] by a significant margin. We can also employ cut pursuit to more complex settings, such as the regularization of nonconvex and nondifferentiable functions, and observe comparable improvements.

In all experiments, the cut pursuit algorithms converge in less than 10 iterations, translating to faster recovery than iterative methods as detailed in Table 2.1. In the point cloud smoothing experiment, the cut pursuit algorithm performs as well as the highly specialized Parametric Max-Flow formulation. Parallelizing the refinement steps brings a significant speed-up, and is further improved when using the component balancing strategy.

In the brain source identification, the final partition contains only 20 constant component. This setting advantages cut pursuit, which can solve this complex, poorly conditioned, and medium-sized problem several orders of magnitude faster than iterative methods. However, parallelization does not bring significant improvements in this setting, and the balancing strategy leads to a suboptimal solution.

Cut pursuit also brings significant acceleration to the semantic segmentation regularization problem. Despite the non-applicability of our convergence guarantees, our method can retrieve high-quality solutions. The balancing strategy further accelerates the recovery.

**Warm Restart Strategy.** In contrast to methods using dual [CD09] or auxiliary variables [Rag18, RFP13], the cut pursuit algorithm can leverage a
<table><caption>Table 2.1: **Recovery Speed.** Time for the method to reach a solution within $10^{-2}$ of the approximate optimal. PMF cannot be used to minimize the functionals in the second and third experiments. The cut pursuit algorithm offers acceleration of several orders of magnitude.</caption><thead><tr><th></th><th>Reflectance<br>Regularization</th><th>Brain Source<br>Identification</th><th>Segmentation<br>Regularization</th></tr></thead><tbody><tr><td>PFDR</td><td>1123</td><td>1.40</td><td>838</td></tr><tr><td>PMF</td><td>732</td><td>-</td><td>-</td></tr><tr><td>CP</td><td>511</td><td><b>0.04</b></td><td>271</td></tr><tr><td>PCP</td><td>255</td><td><b>0.04</b></td><td>145</td></tr><tr><td>PCP-Bal.</td><td><b>150</b></td><td>0.11</td><td><b>39</b></td></tr></tbody></table>

candidate iterate believed to be close to the solution with a simple warm restart strategy. Indeed, instead of starting the current partition $V$ with $\{V\}$, we can use the maximal constant connected components of the candidate iterate. When computing a sequence of problems with decreasing regularization strength, we can use the previous partition as initialization, see [LO17][2.8]. This strategy allows us to approximate the entire regularization path of TV for a cost comparable to solving the lowest regularization setting.

**Conclusion.** We have provided a theoretical and practical framework for harnessing the speed of efficient graph-cut algorithms for a large class of graph-structured problems involving nondifferentiable terms alongside the total variation. Cut pursuit is not meant to be an all-purpose algorithm. It should only be used if the solution is expected to be piecewise constant w.r.t a coarse partition, and this is not generally the case with natural images. However, we believe that cut pursuit addresses some of the limitations of common solvers for such problems:

(i) Cut pursuit is both fast and versatile as does not require convexity or differentiability. In contrast, competing methods are either fast but very specialized (PMF) [CD09], or generic but slow [Rag18].

(ii) Cut pursuit retrieves a solution that is *exactly* piecewise constant. In contrast, the iterative method requires a thresholding step which can engender artifacts.

(iii) By construction, the cut pursuit algorithm retrieves a coarse solution. In contrast, iterative methods generally find solutions with many small constant components due to the *staircasing effect* of total variation [Ja16].

(iv) Cut pursuit allows us to run the fast Ford-Fulkerson graph cut algorithm in parallel to solve continuous optimization problems.
<figure><img src="image_15.png" /><figcaption>Figure 2.6: Relative distance to the approximate optimal against the running time of the algorithms; optimal values are estimated by longer, high-precision runs.</figcaption></figure>

## 3 $\ell_0$-Cut Pursuit for Contour Regularizing

We now consider the case where $h(x_u, x_v) = [x_u - x_v]$ with $[\cdot] : \Omega \to \mathbf{R}$ equals to 0 at $0_\Omega$ and 1 elsewhere. In other words, the graph-structured regularization in Eq. (2.2) becomes the cut between constant components, i.e. the total contour length of the induced partition. This regularization can be linked to several other commonly used regularizers and problems. We can see this penalizer as a continuous-space version of the Potts penalty [Pot52], a generalized version of the minimal partition problem in which the number of components is not predefined [Lec89], a discrete (graph) version of *Caccioppoli partitions* [TC96], or an $\ell_0$ version of the total variation
following the terminology of Candes et. al for the LASSO penalty [CWB08].

In this setting, the functional $F$ is neither continuous nor convex. Consequently, finding a global optimal is generally not an option regardless of the convexity of $f$. However, we show that our approach can be adapted into a greedy formulation, dubbed $l_0$-cut pursuit, which can find quickly and reliably good local minima of $F$. We restrict ourselves to the setting in which $f$ is separable into continuous functions: $f(x) = \sum_{v \in V} f_v(x_v)$ with $f_v: \Omega \to \mathbb{R}$. This assumption is not necessary but greatly simplifies the formulation of the following sub-problems.

## 3.1 Reduction.

For a given partition $\mathcal{V}$ with reduced graph $\mathcal{G} = (\mathcal{V}, \mathcal{E}, w)$, we define the reduced functional $F^{(\mathcal{V})}$ for the reduced variable $r \in \Omega^{\mathcal{V}}$ as follows:

$$F^{(\mathcal{V})}(r) = \sum_{U \in \mathcal{V}} \sum_{v \in U} f_v(r_U). \qquad (2.21)$$

Note that we drop the edge-wise regularization term altogether. The ratio-nale is that during the reduction step, we only change the values of components and do not add nor remove existing *borders*. Since the $l_0$ norm penalizes all non-zeros values identically, we can completely ignore the penalty term in the reduction phase. We can now solve the reduced problem for each component independently and in parallel. A caveat of this simplification is that adjacent components with very close values may benefit from adjusting their iterates to remove their mutual borders. This case is handled by a backward step defined below.

## 3.2 Refinement.

In the original cut pursuit formulation, we perform the refinement step by minimizing an energy determined by the directional derivative of $F$. Since the functional $F$ is non-continuous in the $l_0$ setting, we do not have access to such values. Instead, we define a greedy counterpart of the steepest binary partition that we dub the *optimal binary partition*, which we use to split the current partition $\mathcal{V}$. The optimal binary partition $B^*$ is a solution to the following bilevel optimization problem:

$$B^* \in \underset{B \subset \mathcal{V}}{\operatorname{argmin}} \min_{h, h' \in \Omega^{\mathcal{V}}} \sum_{U \in \mathcal{V}} \left[ \sum_{v \in U \cap B} f_v(h_U) + \sum_{v \in U \cap B^c} f_v(h'_U) \right] \\ + \sum_{U \in \mathcal{V}} w(U \cap B, U \cap B^c), \qquad (2.22)$$

with $B^c$ the complementary of $B w.r.t \mathcal{V}$. The rationale of this step is to split each component $U$ of the current partition $\mathcal{V}$ along a binary partition $(B \cap$
$U/B^c \cap U$), associated with the constant values $h_U, h'_U \in \Omega^2$ respectively. The second term quantifies the additional boundary incurred by the new partition.

The optimization problem defined in Eq. (2.22) is difficult as the second term is still nonconvex and noncontinuous. We can, however, find an approximate solution by alternating between solving for $(h, h')$ and $B$. The partition retrieved is provably a local minimum of the nonconvex optimization problem defined in Eq. (2.22). In practice, only a few (two or three) iterations are necessary.

## 3.3 Backward-Step.

It is common for greedy working set methods to implement a backward step in which previous decisions can be reversed, see the Single Best Replacement method of Soussen et. al [SIBD11]. In our setting, this amounts to removing existing borders by merging adjacent components. While we explored different strategies, the simplest one is, in most cases, sufficient.

We perform this step after the refinement step. We start by computing for each adjacent component $(U, U') \in \mathcal{E}$ in the reduced graph the change in the objective functional $F$ incurred by their merging. If this value is negative, the reduction in contours length outweighs the removal of one degree of liberty incurred by forcing the vertices in $U \cup U'$ to share the same value. Such cases typically happen for borders added at early iterations that are no longer relevant after further refinement. We then greedily select the most profitable move, i.e. which most decreases the objective functional and update the potential gains of all adjacent components. This step is repeated until there are no more profitable fusions. With $f$ vertex-separable, this step is performed by considering the reduced graph only, and is hence very fast.

> While the contour penalty defines a difficult nonconvex, noncontinuous problem, the $\ell_0$-cut pursuit algorithm can efficiently find approximate solutions, outperforming celebrated approaches such as $\alpha$-expansion [BVZ01] with observed accelerations up to 100-fold.

## 3.4 Numerical Experiments.

We evaluate the performance of the $\ell_0$-cut pursuit algorithm for a 3D point cloud segmentation. This large-scale optimization problem consists of grouping points into clusters with homogeneous properties and simple contours, called superpoints. We consider 271 point clouds from the S3DIS dataset [ASZ$^+$16], averaging 116 000 points each after subsampling. In the next chapter, we give more details about this experiment and its motivation.
For a 3D point cloud, we define a graph $G = (V, E, w)$ such that $V$ is the set of 3D points, $E$ encodes their adjacency, and $w$ their proximity. We associate with each vertex $v$ a value $y_v \in \mathbb{R}^7$, encoding its radiometry and local geometry. We aim to compute a piecewise constant approximation of this signal into few components. This amounts to minimizing the following functional for $x \in \Omega^V$:

$$F(x) := \sum_{v \in V} \|x_v - y_v\|^2 + \sum_{(u,v) \in E} w_{u,v}[x_u - x_v]. \quad (2.23)$$

We can approximately minimize this functional with the $l_0$-cut pursuit algorithm. As a baseline, we define the following approach:

* We set the number of classes $K$.

* We run the K-means algorithm on the observations $y \in \mathbb{R}^7 \times V$ to obtain class values $k \in \mathbb{R}^{7 \times K}$. The resulting centroids define a discrete set of values that each vertex can take.

* We compute the combination of classes minimizing $F$. This corresponds to a classical energy minimization problem which can be solved approximately with $\alpha$-expansion.

* We update the class values by computing their feature average. This step is homologous to the debiasing post-processing often applied to the solution of regularized problems [DPS15].

Note that this baseline is only partially comparable to cut pursuit: the former operates on a discrete set of values while the latter works in a continuous space. Since the problem itself is defined w.r.t a continuous domain, we can expect the cut pursuit approach to reach lower energies. However, we can compare the execution time of the two approaches. Indeed, both heavily rely on graph cuts to minimize the energy, albeit in a different manner: $\alpha$-expansion solves a graph cut problem for each class at each iteration, while the cut pursuit algorithm uses the cuts to split existing components. Comparing the execution speed between these methods may inform us which graph cut formulation is more efficient for piecewise constant approximation.

In Figure 2.7, we represent the evolution of the performance of $l_0$-cut pursuit, with or without parallelization, as well as two versions of our baseline with $K = 10$ and $20$. We observe that $l_0$-cut pursuit reaches lower energies, which we expected. More interestingly, our approach provides an acceleration of over 100 fold compared to $\alpha$-expansion. This indicates that the component-wise splitting of $l_0$-cut pursuit might be a much better strategy than the class-wise splitting of $\alpha$-expansion when the sought signal is coarse. Furthermore, we can bring further speed-ups with parallelization.
<figure><img src="image_16.png" /><figcaption>Figure 2.7: **Numerical Experiment for $l_0$-cut pursuit.** We represent the evolution of $F/|V|$ with respect to computation time in seconds. The cut pursuit algorithm produces its approximation in little more than 2 minutes with parallelization, and under 7 minutes without. In contrast, $\alpha$-expansion takes up between 2 and 8 hours depending on the chosen number of labels and desired precision.</figcaption></figure>

# 4 The Plane-Pursuit Algorithm

Our approach is particularly well-suited for data exhibiting spatial regularity and well-represented with graphs, such as 3D point clouds. This setting has attracted most of the interest in cut pursuit in terms of applications: [TGB19, GL17, LRV$^+$17, LS18, MSF20]. In the next chapter, we will detail an approach to scale the analysis of deep networks to point clouds containing millions of points with cut pursuit. This section presents a specialized adaptation of cut pursuit for computing piecewise-planar approximations of large 3D point clouds.

We derive the graph $G = (V, E, w)$ from an unstructured 3D point cloud: $V$ is the set of 3D points, $E$ captures local neighbourhoods, and $w \in \mathbb{R}_+^E$ the proximity of points. We aim to approximate a point cloud with a small number of 2D planes embedded in $\mathbb{R}^3$. We propose to formulate this task as the recovery of a graph-coarse signal on graph $G$.

Let $\mathcal{P}$ be the set of all planes of $\mathbb{R}^3$. We denote by $d(v, \pi)$ the distance between a point $v$ and a plane $\pi \in \mathcal{P}$. We define $\Pi: V \mapsto \mathcal{P}$ as the function that associates each vertex of $V$ with a plane $\Pi_v$ in $\mathcal{P}$. We want $\Pi$ to define as few planes as possible while remaining a good approximation of the point cloud. This task amounts to recovering a signal $\Pi$ on the set of vertices $V$ for which (i) the points are well approximated by their pro-
jection in their respective plane: $d(v, \Pi_v)$ is small and (ii) the signal $\Pi$ is graph-coarse. We can translate these objectives into a unique functional $F$ to minimize:

$$F(\Pi) := \sum_{v \in V} d(v, \Pi_v)^2 + \mu \sum_{(u,v) \in E} w_{u,v} [\Pi_u \neq \Pi_v], \quad (2.24)$$

where $[\pi \neq \pi']$ is the function of $\mathcal{P}^2 \to \{0,1\}$ equal to 0 when $\pi$ and $\pi'$ are identical planes and 1 otherwise. The parameter $\mu \in \mathbb{R}_+$ is the regularization strength. The first term of Eq. (2.24) measures the fidelity of the planar reconstruction defined by $\Pi$. The second part of Eq. (2.24) encourages planes associated with adjacent vertices to be identical.

The functional defined in Eq. (2.24) is similar in structure to Eq. (2.2), and we can apply the same algorithmic scheme as in the previous sections with minor adjustments: the reduction step now perform plane-fitting for each component, and the refinement step use RANSAC [FB81] to initialize the splitting procedure described in Section 3.2.

**Numerical Experiments.** We compare our results with a traditional planar segmentation approach based on a region-growing clustering by Cohen et. al [CSAD04]. For a given number of components/regions, plane pursuit improves the geometric error by an entire order of magnitude and speed by more than 50-fold. In Figure 2.8, we represent qualitative results illustrating the reason behind this performance: plane pursuit fits large planes where the geometry is simple, adapting the size of planes to the local geometry's complexity.

The cut pursuit algorithm can be adapted for the efficient planar reconstruction of unstructured 3D point clouds. The resulting plane-pursuit algorithm takes advantage of the regularity of the solution to provide both a significant acceleration and a better approximation.

## Conclusion

Cut pursuit is a meta-algorithm, which can take different forms depending on the property of the functional to minimize. We provide convergence proofs for a wide range of data functionals when using the total variation as graph regularizer. The cut pursuit algorithm can be seen as an efficient heuristic strategy in more complex settings. In all cases, we provide considerable acceleration when the sought signal exhibits graph-coarseness. This results from the algorithmic design that exploits such structured sparsity but also from an efficient parallel implementation.
<figure><img src="image_17.png" /><figcaption>Figure 2.8: **Qualitative Results for Plane-Pursuit.** We represent the projection of each point on their supporting plane, and each colour represents a different component/region. Our method yields planes of adaptive size: large in geometrically simple regions such as roads, and smaller for complex areas such as facades and window panes.</figcaption></figure>
Exploiting the Spatial
Regularity of 3D Data

This chapter presents different methods for exploiting the spatial regular-
ity of large 3D point clouds to accelerate their automated analysis. We first
propose a versatile family of spatial regularizers for semantic labelling. We
then introduce superpoint graph, a compact yet rich representation of large
point clouds. This structure can be computed efficiently and leads to signif-
icant improvement in terms of precision and scalability. We also present a
hybrid approach combining deep learning and combinatorial optimization
to mesh large point clouds.

This chapter is organized around the following publications:

[LS18]: Loic Landrieu, Martin Simonovsky, "Large-scale Point Cloud Semantic Segmentation with Superpoint Graphs", CVPR, 2018

[LB19]: Loic Landrieu, Mohamed Boussaha, "Point Cloud Oversegmentation with Graph-Structured Deep Metric Learning", CVPR, 2019

[CCHL20]: Chaton, Thomas, Nicolas Chaulet, Sofiane Horache, and Loic Landrieu, "Torch-Points3D: A modular multi-task framework for reproducible deep learning on 3D point clouds", 3DV, 2020

[SLMV21]: Raphael Sulzer, Loic Landrieu, Renaud Marlet, Bruno Vallet, "Scalable Surface Reconstruction with Delaunay-Graph Neural Networks", Symposium on Geometry Processing (SGP), 2021

[LRV<sup>+</sup>17]: Loic Landrieu, Hugo Raguet, Bruno Vallet, Clément Mallet, Martin Weinmann, "A structured regularization framework for spatially smoothing semantic labelings of 3D point clouds", ISPRS Journal of Photogrammetry and Remote Sensing, 2017

[GL17]: Stéphane Guinard, Loic Landrieu, "Weakly supervised segmentation-aided classification of urban scenes from 3D LiDAR point clouds", ISPRS - *International Archives of the Photogrammetry and Remote Sensing and Spatial Information Sciences*, 2017
# 1 Regularizing 3D Point Clouds Segmentation

We present a structured-optimization approach to impose a spatial regularity prior to probabilistic 3D point classifications. Our approach increases performance while retaining the probabilistic nature of the prediction.

**Classical Approaches to 3D Semantic Segmentation.** Steady improvements in sensor technology allow the acquisition of 3D point clouds with high resolution and precision. A straightforward consequence of the point density is that the resulting point clouds exhibit a strong spatial regularity. This property manifests itself in the semantic labels associated with each 3D point: a point *surrounded* by points belonging to the class "car" most likely has the "car" label as well, as seen in Figure 3.1a.

Until the recent rise of 3D deep learning [QSMG17a, GWH$^{+}$20], most 3D point cloud semantic segmentation pipelines classified individual points based on hand-crafted features [WJHM15]. These features are derived from spatial neighbourhoods of each point, which can be fixed [DMDV11], adaptive [WSM$^{+}$15], or multiscale [TGD$^{+}$18]. However, the label associated with each point is ultimately decided independently for each point, leading to the low spatial regularity of the labels as observed in Figure 3.1b.

Many works use graphical models such as Markov Random Fields (MRF) [MBVH09, SVB10, LR12, NNSP14] or their discriminative counterpart, Conditional Random Fields (CRF) [SNRS14] to increase the spatial regularity of the prediction. Computing marginal inference on a large scale is costly and can lead to poor results. Maximum a Posteriori (MAP) inference can be efficiently and reliably approximated with graph cut-based algorithms [BVZ01]. However, this type of inference produces a single label for each point and discards any confidence information from the prediction. Calibrated predictions, as seen in Figure 3.1d, are crucial in 3D analysis, particularly for applications such as autonomous driving [HWB$^{+}$13] or surface reconstruction [CBV17, CKR04].

**Spatial Smoothing as Structured Optimization.** Let $G = (V, E, w)$ be an edge-weighted graph characterizing the adjacency structure of a 3D point cloud indexed by $V$. We consider a probabilistic labelling $\delta \in \Delta^V$ of the points of $V$, with $\Delta_K$ the $K$-simplex when $K$ is the number of classes. We denote $\delta_{v,k}$ the probabilistic prediction associated with the vertex $v$ and class $k$. $\delta$ is typically obtained by classifying each point individually as proposed in [WJHM15]. We propose to advantageously formulate the problem of increasing the spatial regularity of $\delta$ as a structured optimization problem in
<figure><img src="image_18.png" /><figcaption>Figure 3.1: **Smoothing a Pointwise Prediction.**: A 3D point cloud taken from the Paris-rue-Cassette Dataset [VBS+15] and with true annotation (a). A pointwise classifier first gives a noisy semantic segmentation (b) which we smooth to improve the spatial regularity (c). Our approach allows to retain the probabilistic nature of the segmentation (d) : confident <color>: uncertain. Note that misclassifications in (c) correspond to the least confident area in (d).</figcaption></figure>
the domain $\Omega^V$:

$$q^* \in \underset{q \in \Omega^V}{\operatorname{argmin}} \sum_{v \in V} f(\delta_v, q_v) + \sum_{(u,v) \in E} w_{u,v}h(q_u - q_v). \quad (3.1)$$

For different choices for $\Omega$, $f$, and $h$, this versatile formulation allows us to define both classic and new smoothing problems:

* **MAP Inference in a CRF. (MAP)** We set $\Omega$ to $\Delta_K$ the corner of the $K$-simplex, $f(\delta_v, q_v) = \sum_{k=1}^K q_{v,k} \log(\delta_{v,k})$, and $h = [\cdot = 0_\Omega]$ equal to 1 at $0_\Omega$ and 0 otherwise. This formulation corresponds to the MAP problem in a CRF defined by $G$. This combinatorial problem can be efficiently approximated with $\alpha$-expansion, but the solution is a simple labeling without confidence values.

* **TV-regularized KL-Divergence. (KL-TV)** We set $\Omega = \Delta_K$, $f(\delta_v, q_v) = \sum_{k=1}^K \delta_{v,k} \log(q_{v,k})$, and $h(x) = \|x\|_{1,1} = \sum_{k=1}^K |x_k|$. The choice of $f$ corresponds to the (variable part of the) Kullback-Leibler divergence, a similarity measure between distributions based on information theory [KL51]. The choice of $h$ leads to the graph-total variation presented in Section 2.1. The resulting convex problem is typically solved using first-order proximal methods such as the Generalized Forward-Backward Splitting algorithm [RFP13].

* **Continuous Potts with KL-Divergence. (KL-Potts)** In this formulation, we change $h$ to $[\cdot = 0_\Omega]$. This results in a continuous-space, non-convex, and noncontinuous optimization problem, which we have extensively described in Chapter 2. This functional can be efficiently minimized with the $l_0$-cut pursuit algorithm.

As a baseline, we also consider the marginal inference in the CRF defined by $G$ (Marginal), solved with Loopy-Belief Propagation [Pea82].

**Numerical Experiments.** We evaluate these four approaches on three open-access medium-scale 3D scenes containing between 1.3 and 12 million points and between 3 and 7 distinct semantic classes. We first perform adaptive-scale geometric feature extraction [WSM+15] and train a random forest classifier on a small subset of 1000 points to obtain noisy pointwise predictions. We observe in Figure 3.2 that the performance of marginal inference provides significantly lower results than the other methods in terms of the F1-score. MAP inference leads to results of comparable quality to KL-TV and KL-Potts, but loses the probabilistic nature of the prediction. In contrast, we can associate confidence with each prediction of KL-TV and KL-Potts by computing the entropy of the class prediction. As seen in Figure 3.2 that almost all errors occur for the 30% least confident points, illustrating the interest of this approach when precision is critical.
<figure><img src="image_19.png" /><figcaption>Figure 3.2: **Quantitative Results for Label Smoothing.** Accuracy/coverage plot: we sort the points by decreasing confidence for the Oakland-3C dataset (top left), the Oakland-5C dataset (top right) and the Paris-rue-Cassette database (bottom left). The precision reaches nearly 100% for probabilistic methods at coverage 70%, i.e. once removing the 30% least confident points</figcaption></figure>

> We propose a versatile optimization-based framework for spatially regularizing pointwise prediction on 3D point clouds. In particular, our approach increases the precision of the prediction while retaining the probabilistic nature of the pointwise prediction, which is critical for many applications.

## 2 The Superpoint Approach

Following the steps of image analysis, 3D point cloud analysis has now fully adopted the deep learning paradigm [GWH$^{+}$20]. However, neural networks are limited in the size of the inputs they can handle simultane-
ously. LiDAR scans often reach several million points and require sliding
window strategies [QYSG17, CGS19, TQD<sup>+</sup>19, CCHL20]. Although this
is sufficient for problems such as indoor segmentation, the loss of global
structure can be a limiting factor for large-scale geospatial data.

The semantic labels of 3D point clouds are generally regular, as ex-
plored in the previous section. Furthermore, the transitions between ob-
jects typically occur in regions with geometric or radiometric discontinu-
ities. In other words, the interface between objects of different natures is
often characterized by a sharp shift in shape and/or colour—or a learned
combination of such features. Consequently, by grouping adjacent points
with similar local geometry and radiometry, the resulting partition should
also be semantically homogeneous: each component contains mainly the
same semantic or instance labels. Instead of using this prior *ex post* to
improve a noisy prediction, we propose to exploit this property algorithm-
mically to improve processing efficiency. The resulting deep learning ap-
proach we propose can handle millions of points simultaneously.

**Principle of SuperPoint Graph.** We represent large 3D point clouds as collections of interconnected shapes called superpoints, in the spirit of superpixel methods for images [ASS$^{+}$12]. As illustrated in Figure 3.3, we represent this structure with an attributed directed graph called the SuperPoint Graph (SPG). The vertices are the shapes, and the edges represent their adjacency relationship in 3D space. We also equip each edge with a descriptor of the nature of the adjacency between superpoints (e.g. max and average distance, interface size, size ratio).

Instead of considering individual points or voxels, we directly classify
superpoints. Since the number of superpoints is typically several orders of
magnitude smaller than the number of points, this allows us to consider
large point clouds and model long-range interactions. The proposed SPG
representation divides the semantic segmentation problem into three dis-
tinct problems of different complexity:

1 **Semantically homogeneous partition:** We first partition the point
cloud into geometrically simple superpoints. We can easily derive
the SPG from this partition.

2 **Superpoint embedding:** Since superpoints are geometrically simple, we can embed them with a lightweight point set network.

3 **Contextual segmentation:** The SPG typically contains a few hundred vertices and a few thousand edges at most, even for complex scenes. This small scale allows us to employ powerful graph convolution networks to leverage the context of superpoints.

Note that only the partitioning step considers the entire point cloud; all
other steps operate directly on the superpoint graph. An essential hypothe-
sis is that geometrically simple superpoints are semantically homogeneous,
<figure><img src="image_20.png" /><figcaption>**Figure 3.3: Superpoint Graph.** An input point cloud (a) is divided into geometrically simple shapes, called superpoints (b). We construct a graph by linking nearby superpoints with superedges (c). We then use graph convolutions to classify the superpoints according to their shape and context.</figcaption></figure>

i.e. contain mostly points with the same label. In Section 2.1, we first assume that geometric and radiometric homogeneity implies semantic purity. We then present in Section 2.2 a method to directly learn to partition the point cloud into semantically homogeneous regions.

> By grouping points into geometrically and semantically homogeneous regions, we cast the semantic segmentation of a large 3D point cloud as a vertex classification problem on a small graph.

### 2.1 Geometrically Homogeneous Point Cloud Partition

This section details the partition of an input point cloud V into geometrically and radiometrically simple and contiguous regions. Our objective is not to retrieve individual instances, such as cars or chairs, but to break down objects into geometrically simple parts, as seen in Figure 3.3.

**Geometrically Homogeneous Partition.** We associate with each point a set of geometric features based on hand-made dimensionality features
<figure><img src="image_21.png" /><figcaption>Figure 3.4: **Pipeline.** We illustrate our approach on a toy scene with a table and a chair. We partition the point cloud (a), which allows us to build the superpoint graph (b). Each superpoint is embedded with a PointNet network. The embeddings are then refined with message passing along superedges to produce the final labelling (c).</figcaption></figure>

[DMDV11, GL17] and, if available, radiometric information (e.g. colour, intensity). We denote these features by $f \in \mathbb{R}^{d_v \times V}$. We also compute $G_{nn} = (V, E_{nn}, w)$ the k-nearest neighbours *adjacency graph* of the point cloud in 3D space (this is *not* the SPG). Each edge $(u, v)$ in $E_{nn}$ is attributed with a non-negative weight $w_{(u,v)} \in \mathbb{R}^+$ encoding the proximity of the points.

Partitioning amounts to computing a piecewise constant approximation of $f$ with respect to the graph $G_{nn}$. We define the superpoints as the constant connected components of the solution of the following optimization problem:

$$ \underset{g \in \mathbb{R}^{d_v \times V}}{\operatorname{argmin}} \sum_{v \in V} \|g_v - f_v\|^2 + \lambda \sum_{(u,v) \in E_{nn}} w_{u,v}[x_u - x_v], \quad (3.2) $$

with $[x_u - x_v] = 0$ if $x_u = x_v$ and 1 else, and $\lambda$ the regularization strength that determines the coarseness of the resulting partition. By construction, the recovered superpoints $S = \{S_1, \dots, S_k\}$ are geometrically and radiometrically homogeneous regions of the point cloud. The defined functional is neither convex nor continuous, making its exact minimization unrealistic for large point clouds. However, the $l_0$-cut pursuit algorithm introduced in Chapter 2 can quickly find an approximate solution with a few graph cuts.

**Computing SuperEdges.** The SPG is an oriented attributed graph $\mathcal{G} = (S, \mathcal{E}, F)$ whose vertices are the superpoints $S$ and whose edges $\mathcal{E}$ (referred to as *superedges*) represent the adjacency between superpoints. We link superpoints whose points are connected in the Delaunay tetrahedralization [JT92, WH94] of the point cloud. Using this adjacency definition instead of
nearest neighbours allows us to model long-range relationships not typi-
cally captured by nearest neighbour graphs. We attribute each superedge
with features that describe the adjacency relationships between the con-
nected superpoints, such as the offset vector between their centroids and
comparisons between the superpoint size and shapes. Note that the asym-
metry in the superedge features makes the SPG a directed graph.

## 2.2 Learning to Partition 3D Point Clouds

<figure><img src="image_22.png" /><figcaption>Figure 3.5: Learning to Partition. We process a point cloud (a) with our point embedding network (b). We compute a piecewise constant partition of these embeddings (c), which closely follows the ground truth objects (d).</figcaption></figure>

Not all interfaces between objects are characterized by a change in ge-
ometry or colour. For example, whiteboards and walls are flat, white, and
vertical. The transition between these objects is subtle yet immediately rec-
ognizable to the human eye. Additionally, certain objects present heteroge-
neous radiometry and geometry: trees are composed of trunks (flat, brown,
and vertical) and foliage (scattered, green, and volumic). In this section, we
propose learning point descriptors with high contrast along object transi-
tions; see Figure 3.5.

**Graph-Structured Contrastive Loss.** We denote by $\sigma: V \rightarrow S_m^V$ a neural network which extracts local geometric and radiometric features of 3D points, with $S$ the $m$-dimensional hypersphere. This model is typically
<figure><img src="image_23.png" /><figcaption>Figure 3.6: **Contrastive Functionals.** The homogeneity-promoting functional $\phi$ (in blue) and the contrastive functional $\psi$ (in red) used in the graph-structured contrastive loss.</figcaption></figure>

straightforward and operates with a limited receptive field (e.g. 50 nearest neighbours of each point). Our goal is for the vertex-valued function $f = \sigma(V)$ to be homogeneous within objects and to present high contrast at their interface. To this end, we define the set of *intra-edges* $E_{\text{intra}}$ as the edges of $E_{\text{nn}}$ connecting points within the same object. Equivalently, $E_{\text{inter}}$ is the set of *inter-edges* and represents the transition between distinct objects.

We propose to learn $f$ by minimizing the following graph-structured contrastive loss:

$$ \ell(f)|E| := \sum_{(u,v) \in E_{\text{intra}}} \phi(f_u - f_v) + \sum_{(u,v) \in E_{\text{inter}}} \mu_{u,v} \psi(f_u - f_v), \quad (3.3) $$

with $\phi$ (resp. $\psi$) a function minimal (resp. maximal) in 0, and $\mu_{u,v} \in \mathbb{R}^{E_{\text{inter}}}$ a weight on inter-edges. A point embedding function that minimizes this loss should be uniform within objects and have sharp contrast at their interface. Consequently, the components of the piecewise constant approximation of Eq. (3.2) should follow the objects' borders.

We chose $\phi$, the function that promotes intra-object homogeneity as $\phi(x) := \delta(\sqrt{\|x\|^2/\delta^2 + 1} - 1)$ with $\delta = 0.3$. The first term of $\ell$ becomes the pseudo-Huber graph-total variation on the $E_{\text{intra}}$ edges [H$^+$73, CBFAB97], and promotes the homogeneity of embeddings within the same object. With $\psi(x) = \max(1 - \|x\|, 0)$, the second part of $\ell$ penalizes vertices connected across objects' borders and that share similar embeddings. A point embedding constant within objects and with *enough* contrast between adjacent objects will have a loss of 0, which implies perfect segmentation with a well-chosen parameterization of Eq. (3.2). The four-colour theorem
[Gon08] tells us that such embedding can be achieved for hyperspheres with a dimension of 3 or more. However, because $\sigma$ operates on point features, borders that do not present sufficiently distinct geometric or radiometric signatures will remain undetected.

**Cross-Partition Weighting.** Equation Eq. (3.3) views the problem of learning to partition as an edge classification problem. However, a single missed edge can erroneously fuse two large superpoints that cover different objects and drastically impact the quality of the partition. We propose to translate the influence of each edge on the resulting partition into a set of inter-edge weights $\mu_{u,v}$, see Figure 3.7.

We compute the partition $S$ obtained by the piecewise constant approximation of the point features $f$ as defined by Eq. (3.2). We define the cross-partition $V_{cross}$ between $S$ and the true partition $O$ of $V$ into objects: $V_{cross} := \{O \cap S | O \in O, S \in S\}$. For all inter-edge $(u, v) \in E_{inter}$, we denote by $U$ (resp. $V$) the element of $V_{cross}$ containing $u$ (resp. $v$) and associate the following edge weight:

$$ \mu_{u,v} := \mu_0 \frac{\min (|U|, |V|)}{|U \times V \cap E_{inter}|}, \qquad (3.4) $$

with $\mu_0$ a normalizing parameter. The rationale of this formula is the following: since $(u, v)$ is an inter-edge, $u$ and $v$ belong to different objects. If $(u, v)$ is not detected as a transition, then $U$ and $V$ will be erroneously merged Since $U$ and $V$ cover different objects, this would incur at least $\min (|U|, |V|)$ *trespassing* vertices, i.e. vertices that do not belong to the same object as the majority of the vertices of their superpoint. The weights are also divided by the size of the interface between $U$ and $V$ to evenly distribute the penalty over the number of edges that make up an interface. This prevents long borders from being overrepresented in the loss and allows the network to handle *choke points*.

We partition large point clouds into semantically homogeneous regions by computing a piecewise constant approximation of local point descriptors on an adjacency graph. These descriptors can be either hand-crafted or learned to present high contrast at the interface between objects.

## 2.3 SuperPoint Classification

After constructing the superpoint graph, our goal is to assign a semantic label to each superpoint based on its shape and the global context of the scene, see Figure 3.4.
<figure><img src="image_24.png" /><figcaption>**Figure 3.7: Cross-Partition Weighting.** Scene comprised of a door (D) and a wall (W). Two superpoints L (left) and R (right) overlap the door. The edges of (LW, LD) (resp. (RW, RD)) connect the vertices of the left (resp. right) superpoint that belong to the wall and the door. With fewer trespassing vertices and a longer interface than (RW, RD), the weights of (LW, LD) are smaller.</figcaption></figure>

**Superpoint Embedding.** Since superpoints are geometrically simple by construction, we sub-sample them on-the-fly to $n_p = 128$ points and normalize them to the unit cube. This adaptive sampling increases the computational and memory efficiency, and doubles as a powerful data augmentation. We then embed each superpoint $U \in S$ independently and in parallel into a vector $z_U$ of size $n_z = 32$ with a simple shape-embedding network $\pi: \mathbf{R}^{3 \times n_p} \to \mathbf{R}^{n_z}$ such as PointNet [QSMG17b]. To take the superpoint size into account, we concatenate its diameter to the maxpooled shape vector.

**Contextual Segmentation.** In order for each superpoint to refine its embedding based on its local context, we employ a message-passing scheme inspired by Edge-Conditioned Convolutions (ECC) [SK17]:

(i) We initialize the state of each superpoint $U$ with the shape embedding $z_U^{(0)} = z_U$.

(ii) Each superedge $(U, V)$ is assigned a filter $\Theta_{U,V} \in \mathbf{R}^{n_z}$ mapped from its features by an MLP.

(iii) Each superpoint $U$ sends a message $\Theta_{U,V} \odot z_V^{(t)}$ to each of its neighbours $V \in N_U$ in the SPG, consisting of its current state modulated by the edge filter.

(iv) Once all messages are sent, each superpoint $U$ combines the incoming messages and updates its current state accordingly.

(v) After a set number of (iii) $\mapsto$ (iv) steps, the last current state is mapped to a vector of class scores by another MLP.
The following equation summarizes the main mechanism behind this scheme:

$$z_{U}^{(t+1)} = \text{update} \left( z_{U}^{(t)}, \text{pool} \left( \left\{ \Theta_{U,V} \odot z_{V}^{(t)} \right\}_{V \in N_{U}} \right) \right), \quad (3.5)$$

with pool the averaging operator, $\odot$ the Hadamard product, and update a modified Gated Recurrent Unit (GRU) [CvMG$^{+}$14]. We added a gating mechanism allowing the recurrent network to block parts of the input vector depending on its current state. The GRU can ignore specific channels in straightforward cases, freeing their usage for ambiguous situations. Note that the superpoint and superedge embedding networks, as well as the GRU, are all trained end-to-end simultaneously.

This message passing scheme is similar to some inference methods used in CRFs, such as Loopy Belief Propagation [MWJ13], or its numerous deep learning-based reformulations [ZJR$^{+}$15, SU15, LSvdHR16, CK16, LKZ$^{+}$17]. The main difference is that the messages operate directly on learned representations instead of relating to class compatibility. This allows the network to leverage weak unaries or to discover latent sub-classes corresponding to object parts, e.g. the legs of tables and chairs.

## 2.4 Numerical Experiments

**Datasets.** We evaluate the performance of our 3D point cloud analysis algorithm on three datasets of different natures:

* **S3DIS [ASZ$^{+}$16]**. This indoor dataset of office buildings contains more than 278 million semantically annotated 3D coloured points in 6 building areas and annotated with 13 classes and individual instances.

* **Semantic3D [HSL$^{+}$17]**. This large outdoor dense RGB LiDAR dataset consists of over 3 billion points from various urban scenes. The dataset consists of 15 training scans and 15 test scans with withheld labels and a reduced test set of 4 subsampled scans.

* **vKITTI3D [GWCV16, EKHL17]**. This autonomous driving dataset is made up of virtual LiDAR scans [GWCV16] with point annotations [EKHL17].

**Point Cloud Oversegmentation.** We first evaluate the partition quality obtained from hand-crafted and learned features. We compare our results with two state-of-the-art superpoint methods: VCCS [PASW13] and the method of Lin et. al [LWZ$^{+}$18]. We also implemented a 3D version of the deep superpixel algorithms proposed by [JSL$^{+}$18] and inspired by SLIC superpixels [ASS$^{+}$12]. We propose three metrics to evaluate the quality of a partition. Oracle Overall Accuracy (OOA), which is the accuracy obtained when associating each superpoint to its majority label; Border Recall and
<figure><img src="image_25.png" /><figcaption>Figure 3.8: Quantitative Oversegmentation Results. Performance of the different algorithms on the 6-fold S3DIS dataset (a, b, c), and the 6-fold vKITTI dataset (d, e, f). SSP-Cluster and VCCS are not represented for vKITTI for the sake of legibility as their performance is too low.</figcaption></figure>

Precision (BR and BP), assessing how well the contour of objects is correctly recovered (with a set tolerance).

In Figure 3.8, we present the results of the different oversegmentation approaches. Using learned features leads to partitions requiring five times fewer superpoints than competing methods to reach the same purity and border quality. We can observe in the qualitative results displayed in Figure 3.9 that networks learn features with high contrast around objects, even ones with subtle interfaces such as whiteboards and white walls. Furthermore, the features are homogeneous within objects even when the radiometry and geometry change due to light reflections or composite objects such as trees.

**Point Cloud Semantic Segmentation.** We report in Table 3.1 the performance of SPG, with and without learned partition and for different datasets. Since our approach classifies superpoints and not points, its precision is bound by the quality of the partition. However, with only a fraction of the parameters of large and complex neural networks such as KP-Conv [TQD$^{+}$19], or MinkowskiNet [CGS19], our method achieves comparable performance. Learning the partition provides a consistent improvement.
<figure><img src="image_26.png" /><figcaption>(b) vKITTI scene with true 233 objects. Superpoints: SSP 420, VCCS 422, Lin 425.</figcaption></figure>

**Figure 3.9: Qualitative Oversegmentation Results.** Our learned point features are homogeneous within objects and present high contrast at their border, allowing for their precise segmentation.

ment, with few additional parameters (under 15k). SPG training epochs are at least an order of magnitude faster than its competitors, and the entire model can be trained in under 2h on a single commercial GPU, compared
to days on GPU clusters for the most demanding networks.

<table><caption>Table 3.1: **Quantitative Semantic Segmentation Results.** Results given in % of mIoU over several datasets.</caption><thead><tr><th>Model</th><th>Size<br>×10<sup>6</sup></th><th>S3DIS<br>6Fold</th><th>Semantic3D<br>Red/Full</th><th>vKITTI3D</th></tr></thead><tbody><tr><td>PointNet [QSMG17a]</td><td>3.5</td><td>47.6</td><td>-</td><td>34.4</td></tr><tr><td>PointNet++ [QYSG17]</td><td>12.4</td><td>54.5</td><td>-</td><td>-</td></tr><tr><td><b>SPG [LS18]</b></td><td>0.25</td><td>62.1</td><td>73.2/76.2</td><td>55.4</td></tr><tr><td>MinkoNet [CGS19]</td><td>31</td><td>65.9</td><td>-</td><td>-</td></tr><tr><td>ConvPoint [Bou19]</td><td>2.8</td><td>68.2</td><td>-/76.5</td><td>-</td></tr><tr><td><b>SSP + SPG [LB19]</b></td><td>0.25</td><td>68.4</td><td>-</td><td>57.0</td></tr><tr><td>SPNet [HYC<sup>+21</sup>]</td><td>-</td><td>68.7</td><td>-</td><td>57.0</td></tr><tr><td>RandLA-Net [HYX<sup>+20a</sup>]</td><td>1.2</td><td>70.0</td><td>77.4/-</td><td>-</td></tr><tr><td>KPConv [TQD<sup>+19</sup>]</td><td>14.9</td><td>70.6</td><td>-</td><td>-</td></tr></tbody></table>

# 3 Large-Scale Surface Reconstruction

Most surface reconstruction methods rely either on visibility-based energy formulations [BRV16, CBV17, JP11, JP14, VLPK12, ZSH19, LPK09] or shape-based deep learning approaches [GFK<sup>+18</sup>, YFST18, SO20, LZS20, DGF<sup>+19</sup>]. Energy-based methods are typically scalable and robust, and provide useful topological guarantees. Deep learning approaches can learn surface priors directly from training data leading to high-performance in-distribution. This section presents a hybrid method in which a graph neu-ral network predicts the parameter of a graph-cut problem, combining the advantages of both approaches, see Figure 3.10.

## 3.1 Tetrahedralization-based Surface Reconstruction

The point density of a 3D point cloud is highest near the surface of the scanned objects. A standard approach to exploit this regularity is to discretize the 3D space with a Delaunay tetrahedralization [CG<sup>+04</sup>] and directly classify each cell as inside or outside. Indeed, the size of the resulting cells is inversely proportional to the density of points, simplifying the decision for large empty regions of space. The interface between the inside and outside cells determines the predicted surface.

We denote by $\mathcal{T}$ the set of tetrahedron partitioning the 3D space, and $\mathcal{E} \subset \mathcal{T} \times \mathcal{T}$ the edges connecting cells with a shared facet. We aim to find a labelling $l_t$ in $\{0,1\}$ for each tetrahedron $t$, with 0 meaning outside and 1 inside. Each cell is attributed a potential $U_t$ quantifying its likelihood of
<figure><img src="image_27.png" /><figcaption>Figure 3.10: Pipeline. We discretize the input point cloud (a) and compute tetrahedron visibility information from lines of sight —► and from camera positions ⌖. We use a graph neural network (b) to predict an occupancy score for each tetrahedron, which in turn define an energy (c) leading to a minimal cut --- in an adapted flow graph. The reconstructed surface (d) — is the interface between cells with different labels.</figcaption></figure>

being inside based on visibility considerations, i.e. cameras are outside, and
the immediate continuation of a line-of-sight after impact is inside. We also
define binary potential $B_{s,t}$, whose value indicates the propensity of a facet
to be part of the sought surface. The optimal label assignment minimizes
the following energy:

$$
E(l) = \sum_{t \in \mathcal{T}} U_t(l_t) + \lambda \sum_{(s,t) \in \mathcal{E}} B_{s,t}(l_s, l_t), \quad (3.6)
$$

where $\lambda \ge 0$ is the regularization strength. We can find a global minimizer
of the energy $E$ by computing a minimum cut in an appropriate flow graph
$[BVZ01]$ or using a linear programming approach $[BdLGM14]$. Of course,
the quality of the resulting surface highly depends on the relevance of the
unary and binary terms. Handcrafted potentials are usually quick to com-
pute but may be insufficient in complex cases.

We can exploit the irregular density of 3D scans with a partition of
space that is adaptive to the point density. We then use a hybrid ap-
proach to classify each cell as inside or outside. Our method benefits
from the scalability and robustness of energy-based models and the
adaptability of deep learning.

3.2 Visibility-Based Occupancy Prediction

We propose to use a graph neural network to predict unary potentials in-
stead of using handcrafted values. We first associate each cell with a set
of visibility features based on ray traversal from the sensor positions. We
then use the scalable Graph-SAGE [HYL17] graph-convolution scheme to
<table><caption>Table 3.2: **Object Reconstruction.** Evaluated on shapes from Berger et. al [BLN<sup>+</sup>13] with different scanning settings, our approach outperforms both traditional and deep approaches despite only using a fraction of the training set of ConvONet [PNM<sup>+</sup>20].</caption><thead><tr><th>Method</th><th>Chamfer [↓]<br>distance<br>(point ave. %)</th><th>Volumetric [↑]<br>IoU (%)</th><th>Number [↓]<br>of<br>components</th></tr></thead><tbody><tr><td>ConvONet [PNM<sup>+</sup>20]</td><td>2.53</td><td>64.1</td><td>7.6</td></tr><tr><td>IGR [GYH<sup>+</sup>20]</td><td>5.13</td><td>62.6</td><td>38.3</td></tr><tr><td>Poisson [KH13]</td><td>0.74</td><td>86.1</td><td>7.8</td></tr><tr><td>Labatut et al. [LPK09]</td><td>0.72</td><td>86.4</td><td>2.0</td></tr><tr><td><strong>Ours</strong></td><td><strong>0.65</strong></td><td><strong>88.5</strong></td><td><strong>1.1</strong></td></tr></tbody></table>

leverage information from adjacent cells and predict for each *t* a pair of class scores $i_t, o_t$. These scores can be combined to form a soft occupancy prediction with a softmax function. This network is supervised with the Kullback-Leibler divergence between the predicted and true occupancy ratio of each cell. Finally, we define the unary terms according to the occupancy scores:

$$U(l_t) = i_t[l_t = 0] + o_t[l_t = 1], \qquad (3.7)$$

with $[x = y]$ the Iverson bracket, equal to 1 if $x = y$ and 0 otherwise.

## 3.3 Numerical Experiments

We present numerical experiments to show the performance of our reconstruction method for both objects and large-scale scenes. In both settings, our method is only trained on a small synthetic dataset (130 shapes from ShapeNet [CFG<sup>+</sup>15], artificially scanned) and yet outperforms state-of-the-art learning and non-learning-based methods, highlighting its capacity for generalization.

**Object Reconstruction.** We scan virtually 5 shapes from Berger et. al [BLN<sup>+</sup>13], each with 5 noise, resolution, and outlier ratio settings. As seen in Table 3.2, our approach outperforms both traditional and deep learning approaches trained on the entirety of ShapeNet (50k objects).

**Scene Reconstruction.** Even though we only train our model on a small set of synthetic shapes, it can generalize to large real-life scenes, see Figure 3.11. Our method consistently outperforms traditional methods in terms of precision. It is also faster and more memory efficient than deep learning methods such as ConvONet.
<figure><img src="image_28.png" /><figcaption>Figure 3.11: **Qualitative Results on ETH3D.** Our mesh reconstruction method takes as input a dense MVS point cloud (a) and produces a mesh (b), simultaneously preserving fine details and completing missing parts (here textured with [WMG14]). We represent: in (c), a cropped image of a detail from the *terrace* scene of the ETH3D benchmark [SSG+17b]; in (d), the reconstruction by Jancosek and *et. al* [JP14]; and in (e), our reconstruction. Notice the missing staircase and spurious vertical pattern on the concrete wall in (d). In contrast, our method (e) reconstructs part of the staircase as well as the fine-grained wall textures.</figcaption></figure>

Conclusion

We presented three methods to leverage the underlying structure of 3D
point clouds for three tasks: (i) improving pointwise predictions, (ii) an-
alyzing large point clouds, (iii) large-scale surface reconstruction. By ex-
ploiting spatial regularity, our approaches resulted in higher precision, re-
duction in model size, faster training, and higher generalization.
# Exploiting the Structure of 3D Sensors

This chapter presents four methods that exploit the acquisition structure of 3D data for added speed and precision. We study the structure of LiDAR time sequences and show how a sensor-aware method can reach real-time speed without sacrificing accuracy. We then propose an end-to-end approach for merging images and raw 3D point clouds which defines a new state-of-the-art with minimalistic pre-processing. We introduce a simple way to consider the acquisition geometry of 3D scans to infuse visibility information into any deep surface reconstruction methods. Finally, we present a deep learning approach for automated forest inventory that use the ability of LiDAR to penetrate the tree canopy to predict a multi-layer vegetation structure.

This chapter is based on the following publications:

[LAL22]: Romain Loiseau, Mathieu Aubry, Loic Landrieu, "Online Segmentation of LiDAR Sequences: Dataset and Algorithm", ECCV, 2022

[RVL22]: Damien Robert, Bruno Vallet, Loic Landrieu, "Learning Multi-View Aggregation In the Wild for Large-Scale 3D Semantic Segmentation", CVPR, 2022

[SLB+22]: Raphael Sulzer, Loic Landrieu, Alexandre Boulch, Renaud Marlet, Bruno Vallet, "Deep Surface Reconstruction from Point Clouds with Visibility Information", ICPR, 2022

[KLMC22a]: Ekaterina Kalinicheva, Loic Landrieu, Clément Mallet, Nesrine Chehata, "Multi-Layer Modeling of Dense Vegetation from Aerial LiDAR Scans", CVPR Workshops, 2022

[KLMC22b]: Ekaterina Kalinicheva, Loic Landrieu, Clément Mallet, Nesrine Chehata, "Predicting Vegetation Stratum Occupancy from Airborne LiDAR Data with Deep Learning", Journal of Applied Earth Observation and Geoinformation, 2022
# 1 Online LiDAR Segmentation

Roof-mounted spinning LiDAR sensors are widely used by autonomous vehicles [RBG19]. Most semantic datasets [BGM+19, JOWS21, LXG21] and algorithms used for LiDAR sequence segmentation operate on frames corresponding to a 360° degree arc around the sensor. This incurs an acquisition latency of over 100ms, which is incompatible with real-time applications. To remedy this limitation, we introduce HelixNet, the largest dataset of LiDAR time sequences with fine-grained point information, allowing for precise latency estimation. We also introduce Helix4D, a spatio-temporal transformer architecture operating on *slices* of acquisition corresponding to a fraction of a complete sensor rotation, significantly reducing the latency without hampering precision, see Figure 4.1.

<figure><img src="image_29.png" /><figcaption>Figure 4.1: **Online LiDAR Segmentation.** The 3D point sequences of rotating LiDAR follow a complex helix-like structure in space and time, represented in (a). We propose an efficient spatio-temporal transformer operating on angular slices centred on the sensor's position (b). The slices are partitioned into voxels, each gathering information with voxels from past slices to build a large spatio-temporal receptive field. Our proposed model can segment the LiDAR point stream (c) with state-of-the-art accuracy and in real-time.</figcaption></figure>

## 1.1 HelixNet: A Dataset for Online LiDAR Segmentation

We introduce HelixNet, a large-scale, open-access LiDAR benchmark for real-time semantic segmentation algorithms. HelixNet contains 20 3D sequences from the STEREOPOLIS II dataset [PPC+12], corresponding to 10 billion 3D points across 78 800 frames and 8.85 billion individual labels. As shown in Table 4.1, HelixNet is the largest densely annotated open-access rotating LiDAR dataset by a significant margin.
<table><caption>**Table 4.1: Embarked LiDAR Datasets with Semantic Point Annotations.**
With over 8.8B annotated 3D points, HelixNet is 70% larger than SemanticKITTI, and includes diverse scenes spanning 6 different French cities. HelixNet arranges points with respect to the sensor rotation and contains fine-grained information about their release time.</caption><thead><tr><th>Dataset</th><th>labels</th><th>frames</th><th>classes</th><th>span</th><th>format</th></tr></thead><tbody><tr><td>HelixNet (Ours)</td><td>8.85B</td><td>78k</td><td>9</td><td>6 cities</td><td>sensor rotation</td></tr><tr><td>SemanticKITTI [BGM<sup>+</sup>19]5.2B</td><td>5.2B</td><td>43k</td><td>19</td><td>1 city</td><td>frame</td></tr><tr><td>Rellis3D [JOWS21]</td><td>1.5B</td><td>13k</td><td>16</td><td>1 city</td><td>frame</td></tr><tr><td>KITTI-360 [LXG21]</td><td>1.0B</td><td>81k</td><td>37</td><td>1 city</td><td>frame</td></tr><tr><td>A2D2 [GKM<sup>+</sup>20]</td><td>387M</td><td>41k</td><td>38</td><td>3 cities</td><td>frames</td></tr><tr><td>Paris-Lille-3D [RDG18]</td><td>143M</td><td>N/A</td><td>50</td><td>2 cities</td><td>multi-frame</td></tr><tr><td>Toronto3D [TQM<sup>+</sup>20]</td><td>78M</td><td>N/A</td><td>8</td><td>1 city</td><td>multi-frame</td></tr></tbody></table>

We format the sequences to closely follow the LiDAR data stream released by the sensor. More specifically, points are grouped by packets sharing the same *release time*, i.e. the moment they are made available to the segmentation algorithm. We also associate each point with the position and angle of the sensor, rather than only once per frame for SemanticKITTI [BGM$^{+}$19]. This fine-grained information is critical to measuring the acquisition latency (acquisition to release) and inference latency (release to classification). The fine-grained information about the sensor position allows us to group the points into meaningful acquisition slices, a critical step for the approach detailed in the next section.

## 1.2 Helix4D: Fast LiDAR Segmentation with Transformers

We consider a sequence of 3D points acquired by a rotating LiDAR on a mobile platform, which we split into chronologically ordered slices of acquisition. As represented in Figure 4.3, we process each slice with a U-Net architecture [RFB15] with cylindrical convolutions [ZZW$^{+}$21]. At the lowest resolution, a spatio-temporal transformer network connects neighbouring voxels in space and time, resulting in a large receptive field.

**Temporal Slicing** Instead of processing the data frame-by-frame, we propose to split the sequence into slices covering a fixed portion $\Delta\theta \in [0, 2\pi]$ of the sensor rotation, resulting in a shorter acquisition time and lower latency. Choosing $\Delta\theta = 2\pi$ corresponds to the classic frame-by-frame setting and implies an acquisition latency of 104ms in HelixNet or SemanticKITTI [BGM$^{+}$19]. A slice size of $\Delta\theta = 2\pi/5$ leads to an acquisition latency of
<figure><img src="image_30.png" /><figcaption>Figure 4.2: Extracts from HelixNet. Our proposed dataset contains various urban scenes from motorways to pedestrian plazas and historical centers.</figcaption></figure>

21ms, which is more conducive to the real-time processing of driving data.

**Spatio-Temporal Transformer** Inspired by the Cylinder3D model [ZZW+21], we discretize into grids of increasing coarseness. We use a cylindrical convolutional encoder to produce feature maps at the lowest resolution. These features are processed by a spatio-temporal
<figure><img src="image_31.png" /><figcaption>Figure 4.3: **Helix4D Architecture.** A point sequence is split into angular slices and encoded with a cylindrical convolutional encoder $\mathcal{E}^{\text{grid}}$. At the lowest resolution, we apply $W$ consecutive spatio-temporal transformer blocks $\mathcal{T}^w$ with attention spanning current and past slices. The resulting features are up-sampled to full resolution with a convolutional decoder $\mathcal{D}^{\text{grid}}$ and used to classify the individual points.</figcaption></figure>

transformer and, in turn, upsampled to full resolution with a convolutional decoder. The points are then individually classified to produce a semantic segmentation of the sequence. Our simplified architecture results in a lighter computational and memory load but can still learn rich spatio-temporal features thanks to the addition of the transformer module described below.

We consider all non-empty voxels at the lowest resolution, forming a non-strictly ordered time sequence. We apply to each voxel $W$ independent transformer blocks $\mathcal{T}^1, \cdots, \mathcal{T}^W$ successively. For added efficiency, the transformers use a sparse attention scheme by only comparing each voxel $v$ with other voxels within a mask determined by a spatial radius $R$ and a set of rotation offsets $P \subset \mathbb{N}$. In the context of autonomous driving, we choose $R = 6m$ and $P = \{0, 5, 10\}$, corresponding to slices 0.5 and 1 seconds in the past along with the current one. We encode the relative position between voxels based on their spatio-temporal offset in the manner of Wu et. al [WPC$^{+}$21].

We design a simplified transformer architecture [VSP$^{+}$17] with only two learnable modules: a linear block generating the keys and values and the relative positional encoding. We save further computation at inference time by storing in memory the keys, values, and absolute positions of the voxels in past slices with a fixed buffer of $\max(P)$ rotations. This allows us to allocate a large spatio-temporal receptive field to each voxel without supplementary computations.
<figure><img src="image_32.png" /><figcaption>**Figure 4.4: Influence of Slice Size.** We plot the processing time (left, in ms) and precision (right, in mIoU) of different methods for the considered slice sizes, estimated on the validation set of SemanticKITTI [BGM+19]. Methods whose inference time is slower than the acquisition time of the slice (red shaded area) do not meet the real-time requirement.</figcaption></figure>

## 1.3 Numerical Experiments

**Evaluating Online Semantic Segmentation** We evaluate the performance and inference time for our approach and other state-of-the-art methods on our proposed dataset HelixNet and the standard SemanticKITTI dataset. We consider different slice sizes, from $\Delta\theta = 2\pi/5$ to $\Delta\theta = 2\pi$, i.e. frame-by-frame. We measure the inference latency of a segmentation method as the average time between the release of the last point of a slice and its segmentation. To meet the real-time requirement, the classification of a slice must be faster than its acquisition, as slower processing would cause the prediction module to continuously fall behind. Although thinner slices directly reduce acquisition latency, they also make the real-time requirement more strict: a complete turn must be processed in less than 104ms, and a fifth turn must be in at most 21ms.

Because SemanticKITTI [BGM+19, GLSU13] lacks the pointwise information we provide in HelixNet, we cannot directly run semantic segmentation algorithms in the online setting. This forces us to make several approximations concerning the laser alignments, the movement of the vehicle, the rotation speed of the sensor, and the release time of the points. We must also adapt existing frame-by-frame methods to process slice and not full frames. We propose to process the point clouds corresponding to each slice independently and sequentially. This approach restricts the receptive field to the extent of the slices. However, since the sensor moves, the relative
positions given in past slices whose may no longer be valid. By explic-
itly modelling the spatio-temporal offset between voxels, Helix4D does not
suffer from this limitation.

**Analysis.** In Figure 4.4, we report the inference time and mIoU for dif-
ferent slice sizes. Due to its design, the performance of Helix4D is not af-
fected by the slice size. In contrast, competing methods perform worse with
smaller slices. Helix4D yields state-of-the-art accuracy in the frame-by-
frame setting, with mIoU scores only matched by Cylinder3D [ZZW$^+$21].
However, Cylinder3D has 50 times more parameters and is twice slower,
failing to meet the real-time requirement.

Only two approaches can perform real-time classification of slices with $\Delta\theta \le 2\pi/3$: SalsaNeXt [CTA20] and Helix4D. Our approach outperforms SalsaNeXt by over 10 mIoU points for both frame and slices. Our model's total latency (acquisition plus inference time) in the online setting is 40ms (21 + 19ms). However, it reaches the same performance as Cylinder3D evaluated on full frame with a latency of 212ms (104 + 108ms), an acceleration of more than 5 folds. In short, Helix4D is as accurate as the largest and slowest models, with an inference speed comparable to that of the fastest and least accurate models.

> We exploit the helix-like structure of point cloud sequences acquired with a rotating LiDAR mounted on a mobile platform. Instead of waiting for a complete sensor rotation to start inference, we process slices of rotation, drastically decreasing the acquisition and processing latency. Evaluated on our proposed dataset and a classic benchmark, our approach reaches the performance of state-of-the-art methods with a latency reduced by 5× and a model size reduced by 50×.

# 2 Image and LiDAR Fusion

LiDAR scans capture geometric information with high precision; images provide rich textural and contextual cues. We can exploit this complementarity by processing each modality with a dedicated network and projecting onto 3D points the 2D features learned from real [DN18, HZJ$^+$21, JGS19] or virtual images [KYF$^+$20, CLLH19]. However, merging large-scale point clouds and images raises challenges such as mapping pixels with points and aggregating features between multiple views. Current methods rely on a costly mesh reconstruction or specialized sensors to recover occlusions, and use heuristics to merge available image information. In contrast, we propose an end-to-end trainable multiview aggregation model to combine features from images taken at arbitrary positions.
<figure><img src="image_33.png" /><figcaption>Figure 4.5: **Multiview Information.** A 3D point ◊ is seen in several images with different insights. Here, the **green** image contains contextual information, while the **pink** image captures the local texture. The **orange** image sees the point at a slanted angle and may not contain additional relevant information.</figcaption></figure>

## 2.1 Deep Multi-View Aggregation

We consider *P* a set of 3D points and *I* a collection of co-registered images from the same scene. We seek to exploit the correspondence between 3D points and 2D pixels to perform 3D point cloud semantic segmentation with features learned from both modalities. As illustrated in Figure 4.5, the conditions under which a pixel sees a 3D point heavily influence the quality of the image information. Intuitively, closer images contain textural information, while far-away images are more informative about the context. On the other hand, a slanted viewing angle may result in unreliable or irrelevant 2D features. We propose a method to learn to aggregate relevant image features onto 3D point clouds based on such viewing conditions.

**Point-Pixel Mapping.** Combining point clouds and images involves computing point-pixel mappings that takes occlusions into account. This generally necessitates accurate depth maps from specialized sensors [VM98, CHLS17] or a costly meshing step [BGLSA18]. We propose a GPU implementation of the straightforward Z-buffering method [Str82], which efficiently computes the sought mapping directly from point clouds, images, and poses. This algorithm defines for each point *p* the set $v(p) \subset I$ of images in which it is visible. If $i \in v(p)$, we say that image *i* and point *p* are compatible, and denote by $pix(p, i)$ the index of the pixel of *i* which sees *p*.

We associate with each compatible point-pixel pair $(p, i)$ a vector $o(p, i)$ describing the conditions under which the point $p$ is seen in $i$. We use a
set of 8 handcrafted features describing the point-pixel distance, local geometry, viewing angle, distortion, local density, and occlusion rate, i.e. the ratio of neighbouring points also seen in the image. While we could learn to describe the viewing conditions based on learned 2D and 3D features in an end-to-end fashion, this incurs a significant increase in memory load and did not improve performance in our experiments. Our proposed features contain enough information to learn point-image compatibility and are easy to compute.

<figure><img src="image_34.png" /><figcaption>Figure 4.6: **Bimodal 2D/3D Architecture.** We combine a 2D convolutional encoder $\mathcal{E}^{2D}$ and a 3D network composed of an encoder $\mathcal{E}^{3D}$, a decoder $\mathcal{D}^{3D}$, and a classifier $\mathcal{C}^{3D}$ with three different 2D/3D fusion strategies: early (our choice in the experiments), intermediate, and late fusion.</figcaption></figure>

**Learning Multi-View Aggregation.** We denote by $\{f_i^{2D}\}_{i \in I}$ a set of 2D feature maps of width $C$ associated with the images $I$, typically obtained with a convolutional neural network. We aim to transfer these features to the 3D points by exploiting the correspondence between points and images. However, not all viewing images contain equally relevant information for a given 3D point. For example, an image viewing a point from a distance may give important contextual cues, while an image taken close and at a straight angle may give detailed textural information. In contrast, an image in which a point is seen from a slanted angle or under high distortion may not contain relevant information and may need to be discarded. Note that
the relevance of an image depends on the other views: while an image may see a point with a less-than-ideal viewing angle, it may be the only one with contextual cues and therefore should be kept.

To model these complex dependencies, we predict for each compatible point-image pair $(p, i)$ a *quality score* $x(p, i) \in \mathbb{R}$ describing the relevance of image $i$ for point $p$. The scores $x(p, i)$ are defined as a deep set function $[ZKR^+17]$ of the viewing conditions $o(p, i)$ for all images $i$ seeing $p$:

$$ \{x(p, i)\}_{i \in v(p)} = \text{DeepSet}(\{o(p, i)\}_{i \in v(p)}). \quad (4.1) $$

We use a softmax function to compute attention scores $a(p, i)$ in $[0, 1]$ representing the relative relevance of image $i$ for point $p$ among the image set $v(p)$. We associate to the point $p$ the 2D features $\mathcal{P}(f^{2D}, p)$ defined as the sum the 2D features for all pixels seeing $p$ weighted by their respective attention scores $a(p, i)$:

$$ \mathcal{P}(f^{2D}, p) = \sum_{i \in v(p)} a(p, i) \phi(f_i^{2D}[\text{pix}(p, i)]), \quad (4.2) $$

with $\phi$ a learned linear function. The vector $\mathcal{P}(f^{2D}, p)$ corresponds to the image information gathered from all views of $p$ according to their relevance. In Figure 4.7, we observe the influence of several viewing conditions on the attention scores.

We propose to exploit the synergy between 3D point clouds and images by learning a multiview aggregation scheme based on the viewing conditions of points in images: distance, viewing angle, etc. Combined with standard 2D and 3D networks, our methods define a new state-of-the-art while only requiring raw point clouds, images, and poses. In contrast, all other methods operate on colourized point clouds, and 2D/3D fusion approaches require either a meshing step or specialized depth sensors.

**Extensions.** We propose two extensions of this approach: view gating and feature grouping. View gating is a mechanism that allows us to completely block the image features for a point if no viewing conditions are satisfactory. The gating parameter $g(p) \in [0, 1]$ for a point $p$ is a non-decreasing function of the maximum class score $x_{(p,i)}$ for all viewing images $i \in v(p)$. We then multiply $\mathcal{P}(f^{2D}, p)$ by this parameter, allowing the module to block possibly corrupted information from poor scanning conditions.

The motivation for feature grouping is that there may be several manners for an image to be relevant or not for a point $p$, such as contextual cues, textural information, and colour. We group the channels of
<figure><img src="image_35.png" /><figcaption>Figure 4.7: **Influence of Viewing Descriptors.** Given a point ✧ seen in an image (top left), we compute the quality scores when varying the viewing conditions from their initial values •. We observe some interesting behaviours for the aggregation module which can combine multiple viewing conditions in a nontrivial way. Top right, images with a high occlusion rate despite a straight viewing angle are blocked from transferring information. Bottom left, the module selects close images, but not too close—except if the surface is planar. Bottom right, images with a high occlusion rate despite a low point density are discarded.</figcaption></figure>

$\phi(f_i^{2D}[\text{pix}(p, i)])$ into $K$ contiguous groups. We now associate for each point-image pair $(p, i)$ a set of $K$ quality values, one for each feature group. We then define $K$ attention scores, each measuring the relative importance of images for its feature group. After their weighted summation, we con-
catenate the features channelwise.

**Modality Fusion Network.** We represent in Figure 4.6 the different fusion strategy that we have explored. In practice, we observed that replacing the RGB values of point clouds in colorized point clouds with the pooled image features leads to the best performance. We use the implementation of the Minkowski Engine [CGS19] from torch-Points3D [CCHL20] as 3D backbone, and a pretrained ResNet18 network for embedding images.

## 2.2 Numerical Experiments

<figure><img src="image_36.png" /><figcaption>Figure 4.8: **Qualitative illustration.** Scenes from our considered datasets (top: S3DIS, middle: ScanNet, bottom: KITTI-360) with (a) colorized point clouds, (b) ground truth point annotations, (c) prediction of the backbone network operating on the colorized point cloud, and (d) our method operating on raw uncolored point clouds and images. Our approach can use images to resolve cases in which the geometry is ambiguous or unusual, such as a large amphitheatre with tiered rows of seats (top row).</figcaption></figure>

We present in Table 4.2 the performance of our approach compared to state-of-the-art semantic segmentation methods. Using our image fusion module consistently improves the performance of our 3D backbone by a large margin of 3 to 5 points. Remarkably, we outperform the recent Point-Transformer network even though its performance is much higher than our backbone. We also perform better than all 2D/3D fusion methods. Note that all evaluated methods use colorized point clouds except ours, which only uses raw point clouds and images. See Figure 4.8 for qualitative illustrations. The gating mechanism accounts for 0.4 to 3 mIoU points depending on the considered dataset, and channel grouping for 0.4 to 4.8 points.
## 3. SURFACE RECONSTRUCTION WITH VISIBILITY INFORMATION 71

<table><caption>Table 4.2: **Quantitative Evaluation.** Mean Intersection-over-Union of different state-of-the-art methods on S3DIS's Fold 5 and 6-fold, and KITTI-360 Test. All methods except the last line are trained on colorized point clouds. **State-of-the-art**, second highest.</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">S3DIS</th><th colspan="2">KITTI</th></tr><tr><th>Fold 5</th><th>6-Fold</th><th>Val</th><th>360 Test</th></tr></thead><tbody><tr><td colspan="5"><em>Methods operating on colorized point clouds</em></td></tr><tr><td>PointNet++ [QYSG17]</td><td>-</td><td>56.7</td><td>67.6</td><td>35.7</td></tr><tr><td>SPG+SSP [LS18, LB19]</td><td>61.7</td><td>68.4</td><td>-</td><td>-</td></tr><tr><td>MinkowskiNet [CGS19]</td><td>65.4</td><td>65.9</td><td>72.4</td><td>-</td></tr><tr><td>KPConv [TQD<sup>+19</sup>]</td><td>67.1</td><td>70.6</td><td>69.3</td><td>-</td></tr><tr><td>RandLANet [HYX<sup>+20b</sup>]</td><td>-</td><td>70.0</td><td>-</td><td>-</td></tr><tr><td>PointTrans.[EBD21]</td><td><strong>70.4</strong></td><td><u>73.5</u></td><td>-</td><td>-</td></tr><tr><td>Our 3D Backbone</td><td>64.7</td><td>69.5</td><td>69.0</td><td><u>53.9</u></td></tr><tr><td colspan="5"><em>Methods operating on point clouds and images</em></td></tr><tr><td>MVPNet [JGS19]</td><td>62.4</td><td>-</td><td>68.3</td><td>-</td></tr><tr><td>VMVF [KYF<sup>+20</sup>]</td><td>65.4</td><td>-</td><td><strong>76.4</strong></td><td>-</td></tr><tr><td>BPNet [HZJ<sup>+21</sup>]</td><td>-</td><td>-</td><td>69.7<sup>1</sup></td><td>-</td></tr><tr><td>3D Backbone+</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>DeepViewAgg (ours)</td><td><u>67.2</u></td><td><strong>74.7</strong></td><td>71.0</td><td><strong>58.3</strong></td></tr></tbody></table>

# 3 Surface Reconstruction with Visibility Information

Most deep surface reconstruction methods from point clouds ignore pose information and only operate on point location. However, sensor visibility holds valuable information regarding space occupancy and surface orientation. This section presents two simple ways to augment point clouds with visibility information that surface reconstruction networks can leverage with minimal adaptation. Our proposed modifications consistently improve the accuracy of generated surfaces as well as the generalization capability of the networks to unseen domains.

## 3.1 Visibility-Augmented Point Clouds

By incorporating visibility information into common Deep Surface Reconstruction (DSR) architectures, we aim to improve the occupancy [LPK09, VLPK12, JP11] and orientation [SSG17a] of the predicted surface, see Figure 4.9. While DSR methods have achieved impressive results due to their ability to learn shape priors [PFS<sup>+19</sup>, MON<sup>+19</sup>] and shape similarities [HMGCO20, DDN20], they tend to generalize poorly to unseen categories or settings. We aim to increase their generalization capacities by learning a more grounded and *physical* visibility model.

Real-world point cloud acquisition techniques, such as LiDAR scanning or multi-view stereo (MVS), naturally provide the sensor position. We can
<figure><img src="image_37.png" /><figcaption>Figure 4.9: **Reconstruction with Visibility Information.** We augment each 3D point ● with a sightline vector → and two auxiliary points are placed before ● and after ● the observed point. We can then easily adapt surface reconstruction algorithm to reconstruct significantly more accurate surfaces.</figcaption></figure>

then form lines of sight between the sensor and the acquired points. For each point, we associate (i) a unit *line-of-sight vector* pointing in the direction of the sensor; (ii) two auxiliary points situated before and after the point along the line of sight, see Figure 4.10 for an illustration. By construction, the *before* is likely outside the scanned object, and the *after point* is likely inside. We position these auxiliary points at a distance *d* on each side of the real point, where *d* is a typical level-of-detail value, e.g. the scanning resolution at a set distance.

We can adapt most DSR networks to handle visibility-augmented point clouds with only two simple modifications:

* We concatenate the 3 coordinates of the line-of-sight vectors to the point location and features.

* We add auxiliary points to the point cloud, thus tripling the num-
<figure><img src="image_38.png" /><figcaption>Figure 4.10: **Visibility-Augmented Point Cloud.** Each observed point is associated with a sightline unit vector pointing toward its sensor. New points are added before and after each point to help disambiguate occupancy.</figcaption></figure>

ber of input points. We can perform a volumetric sampling for each
category of points (before, after, and real points) for efficiency.
Since we only add point and input features, most existing architectures can
be used as is. As most of the memory and computation is taken up by
the computation of wide feature maps for subsampled clouds in the inner
layers, our modification typically does not result in significant overheads.

**3.2 Numerical Experiments**

We report the results of experiments on object- and scene-level point clouds by evaluating 5 of the top-performing DSR approaches with and without our proposed improvements. Because sensor position is typically not given in current datasets and benchmarks, we perform a synthetic scan of all 3991 shapes from ModelNet10 to form a synthetic training set.

**Object-Level Reconstruction** We report in Table 4.3 the performance on ModelNet10 for various models, with and without sightline vectors or auxiliary points. We observe a general performance improvement across all methods, as illustrated in Figure 4.11. In practice, adding line-of-sight vectors and auxiliary points is comparable to directly adding the *true oriented normals* from the sought surface in terms of performance. However, these normals are not accessible without ground truth surface, while sensor positions are freely given.

The effect on the run time of adding sightline vectors and auxiliary
points depends on the method: under 2% increase for POCO and Con-
<figure><img src="image_39.png" /><figcaption>Figure 4.11: **Object-Level Reconstruction.** Shapes from the ModelNet10 test set reconstructed using four different DSR methods operating on point clouds with and without our proposed visibility-augmentation.</figcaption></figure>

vONet, 25% for Shape As Points, and 125% for Points2Surf whose bottle-neck is the nearest neighbour computation between points.
<table><caption>Table 4.3: **Object-Level Reconstruction.** We report the volumetric IoU (IoU), average Chamfer distance (CD) and normal consistency (NC) for DSR methods trained and tested on ModelNet10, with and without line-of-sight vectors (SV) or auxiliary points (AP). † Trained on ShapeNet.</caption><thead><tr><th>Model</th><th>SV</th><th>AP</th><th>IoU ↑</th><th>CD ↓</th><th>NC ↑</th></tr></thead><tbody><tr><td>ConvONet-2D [PNM<sup>+</sup>20]</td><td></td><td></td><td>0.853</td><td>0.618</td><td>0.934</td></tr><tr><td>ConvONet-2D [PNM<sup>+</sup>20]</td><td>✓</td><td></td><td>0.871</td><td>0.557</td><td>0.936</td></tr><tr><td>ConvONet-2D [PNM<sup>+</sup>20]</td><td>✓</td><td>✓</td><td><strong>0.889</strong></td><td><strong>0.508</strong></td><td><strong>0.944</strong></td></tr><tr><td>ConvONet-3D [PNM<sup>+</sup>20]</td><td></td><td></td><td>0.885</td><td>0.493</td><td>0.949</td></tr><tr><td>ConvONet-3D [PNM<sup>+</sup>20]</td><td>✓</td><td></td><td>0.911</td><td>0.424</td><td>0.956</td></tr><tr><td>ConvONet-3D [PNM<sup>+</sup>20]</td><td>✓</td><td>✓</td><td><strong>0.923</strong></td><td><strong>0.393</strong></td><td><strong>0.959</strong></td></tr><tr><td>Points2Surf [EOMW20]</td><td></td><td></td><td>0.842</td><td>0.590</td><td>0.890</td></tr><tr><td>Points2Surf [EOMW20]</td><td>✓</td><td></td><td><strong>0.859</strong></td><td><strong>0.544</strong></td><td>0.896</td></tr><tr><td>Points2Surf [EOMW20]</td><td>✓</td><td>✓</td><td>0.856</td><td>0.548</td><td><strong>0.897</strong></td></tr><tr><td>Shape As Points [PJL<sup>+</sup>21]</td><td></td><td></td><td>0.903</td><td>0.438</td><td>0.948</td></tr><tr><td>Shape As Points [PJL<sup>+</sup>21]</td><td>✓</td><td></td><td>0.907</td><td>0.430</td><td>0.950</td></tr><tr><td>Shape As Points [PJL<sup>+</sup>21]</td><td>✓</td><td>✓</td><td><strong>0.914</strong></td><td><strong>0.410</strong></td><td><strong>0.954</strong></td></tr><tr><td>POCO [BM22]</td><td></td><td></td><td>0.907</td><td>0.422</td><td>0.945</td></tr><tr><td>POCO [BM22]</td><td>✓</td><td></td><td>0.915</td><td>0.408</td><td><strong>0.950</strong></td></tr><tr><td>POCO [BM22]</td><td>✓</td><td>✓</td><td><strong>0.917</strong></td><td><strong>0.406</strong></td><td><strong>0.950</strong></td></tr><tr><td>† LIG [JSM<sup>+</sup>20]</td><td></td><td></td><td>–</td><td>0.974</td><td>0.849</td></tr><tr><td>† LIG [JSM<sup>+</sup>20]</td><td>✓</td><td></td><td>–</td><td><strong>0.880</strong></td><td><strong>0.882</strong></td></tr><tr><td>DGNN [SLMV21]</td><td>✓</td><td></td><td>0.866</td><td>0.543</td><td>0.884</td></tr></tbody></table>

**Out-of-Domain Reconstruction** We evaluate the impact of adding visibility information to existing networks in terms of generalization capacity. We trained POCO on ModelNet and tested it on large-scale scenes from ScanNet (real) and SceneNet (synthetic). In both cases, the visibility-augmented models produce smoother reconstructions with better completeness, see Figure 4.12. Quantitatively, the surface IoU increased by almost 3%.

We also evaluate out-of-category reconstructions by testing models trained on ModelNet on shapes from ShapeNet [CFG$^{+}$15]. We observe a significant performance increase for all evaluated methods, up to +44% surface IoU points. This illustrates the benefit of learning a visibility model instead of shapes. Lastly, we evaluate models trained on synthetic shapes from ModelNet10 on real-world scans acquired with LiDAR or MVS. In Figure 4.13, we show that all considered networks reconstruct more accurate surfaces using visibility information.
<figure><img src="image_40.png" /><figcaption>Figure 4.12: **Out-of-Domain Scene-Level Reconstruction.** POCO trained on ModelNet10, with and without visibility information, is run on scenes from SceneNet (synthetic RGB-D scan) and ScanNet (real RGB-D scan).</figcaption></figure>

We propose to incorporate visibility information into deep surface reconstruction methods by adding for each point a unit vector pointing towards the sensor and a pair of auxiliary points along this vector. This only requires minimal adaptation of the considered architecture and results in consistently improved predicted surfaces. More strikingly, the generalization capacity of the networks is significantly increased, and models trained on small synthetic objects can be applied to unseen classes, entire scenes, and even real scans.
<figure><img src="image_41.png" /><figcaption>Figure 4.13: **Out-of-Domain Reconstruction.** Reconstructed shapes from a LiDAR point cloud (left, *Ignatius* from Tanks And Temples) and a MVS point cloud (right, *TempleRing* from Middlebury) using four different DSR methods trained on ModelNet10, with and without visibility information.</figcaption></figure>
<figure><img src="image_42.png" /><figcaption>Figure 4.14: **Multi-Layer Forest Analysis.** We introduce WildForest3D, a novel dataset of annotated UAV-LS point clouds of dense forest (left), and a new model for the multi-layer analysis of vegetation. Our network performs 3D semantic segmentation (middle), and produces height maps and watertight meshes for three vegetation layers (right).</figcaption></figure>

# 4 Forestry Analysis from Aerial LiDAR

The automated analysis of vegetation is one of the key challenges of environment monitoring, allowing us to measure carbon stock and biodiversity [LKA+22, FSM+16]. Since cameras and low-resolution LiDAR cannot capture meaningful information underneath the tree tops [TVB+10], most methods focus on analyzing the canopy [HCZ17]. However, such analysis is limited: (i) tree height makes a poor substitute for biomass [Lu06], (ii) canopy-only analysis ignores the critical biodiversity and natural habitats from the understory [GAMP12], and (iii) ground vegetation is pivotal for forest fire modelling [MR11]. On the other hand, aerial LiDAR scanners with sufficient resolutions can perfectly capture precise geometric information below tree tops [FMJ+15]. We propose to directly model the multi-layer structure of natural forests, whose observation is unlocked by LIDAR.
**4.1 Multi-Layer Vegetation Modeling**

We release the first dataset with pointwise information for natural forests
and introduce a network to explicitly model the multi-layer structure of
natural forests.

**WildForest3D.** We introduce a first-of-its-kind dataset of 29 scans of dense natural forest, which contains 7 million 3D points and 2.1 million individual labels. This corresponds to over 2000 tree instances which were individually located and classified by *in situ* observations from forestry experts. The labels indicate the nature of the vegetation: deciduous canopy, coniferous canopy, understory, stems, and ground. We also produce high-resolution maps of the occupancy and thickness of three vegetation layers: ground vegetation, understory, and overstory.

<figure><img src="image_43.png" /><figcaption>Figure 4.15: **Multi-Layer Modeling.** Our network yields label predictions for each 3D point, which are projected onto rasters to obtain soft occupancy maps for 3 different vegetation layers. The network is supervised using 2D and 3D annotations. During inference, the predictions are used to derive the minimum and maximum elevation of each layer, which we convert in a watertight 3D mesh.</figcaption></figure>

**Layer-Wise Modeling.** We design a model for the automated analysis of multi-layer vegetation from aerial laser scans. Our network operates directly on the 3D points to perform semantic segmentation of the point clouds and generate layer occupancy rasters, see Figure 4.15. Both 2D and 3D tasks are supervised end-to-end and simultaneously. Once trained, our
<figure><img src="image_44.png" /><figcaption>Figure 4.16: **Qualitative Results.** In the top row, we show the ground truth for the point labels and layer occupancy maps, with un-annotated 3D points in grey. In the bottom row, we report our predictions with errors in red.</figcaption></figure>

model produces height maps for all layers, which we transform into watertight meshes. This surface-based representation is helpful for downstream applications such as biomass, carbon stock, and fire fuel estimation [FBJ+12, GCVY21], soil illumination [HvOW03], or vegetation parameter extraction for the forest inventory [AHK+21], see Figure 4.14.

Trained on our WildForest3D, our model can accurately predict the height and thickness of all vegetation layers at a resolution of 0.25m². The estimation error is divided tenfold compared to classical approaches such as regression trees and linear regressions operating on handcrafted descriptors, see Figure 4.16 for a qualitative illustration.

## 4.2 Weakly-Supervised Learning for Forestry Analysis

Pointwise annotation is exceptionally costly and tedious for forestry data. Contrary to most 3D modalities whose scans can be annotated by specialized companies, forest data require experts to take physical measurements and determine tree species *in-situ*. This task is particularly complex for natural forests, which may be poorly accessible, and whose trees are often interlocking. Experts produce tree-wise annotations that need to be converted to a useful format for training machine learning algorithms, such as point labels, which is a laborious and error-prone process. We propose a weakly-supervised scheme from coarse annotations to remediate this critical limitation. The *in-situ* forestry experts are asked to estimate the coverage ratio of three vegetation layers (overstory, understood, ground vegetation)
tion) in a 10m radius around a geo-localized position, see Figure 4.17. This task is much simpler than detailed tree inventory and can be performed in minutes for each plot of over 300m². We release a novel dataset of 200 cylindrical plots with dense 3D scans and cover ratios to evaluate if such annotations are sufficient for training deep models.

<figure><img src="image_45.png" /><figcaption>Figure 4.17: **Plot-Based Annotations.** We represent in (a) the 3D point clouds corresponding to two distinct plots. In (b), we show an aerial view of the plots, and in (c) the perspective of the *in situ* annotator. The top plots has a layer coverage of 0%, 10%, 50% for overstory, understory, ground vegetation, and the bottom plot 40%, 60%, 50%.</figcaption></figure>

We consider the problem of strata occupancy regression presented in the last section. We train a network that predicts fine-grained occupancy maps, only supervised with plot-aggregated coverage ratio. We combine this approach with two auxiliary losses: (i) pixel-wise negative cross-entropy to discourage low confidence occupancy (i.e. 50%) and produce *crisper* maps; (ii) we model the distribution of point elevation as a mixture of two Gamma distributions, which can be fitted with the expectation-conditional-maximization (ECM) algorithm [YCHNP19]. Models trained on purely plot-aggregated data can predict occupancy maps at 0.25m² with high precision.
We exploit the unique perspective offered by aerial LiDAR on the multi-layer structure of natural forests. We propose a model that performs jointly 3D semantic segmentation and layer occupancy regression. Our model yields a precise estimation of the density and height of different vegetation strata, which is helpful for a variety of forestry applications. We also propose a weakly-supervised training scheme that greatly simplifies the annotation process while yielding precise occupancy maps.

## Conclusion

We have shown through four different 3D models and applications the advantage of considering the specificities of the sensor when designing network architectures. Considering the viewing conditions of images in point clouds allows us to overtake state-of-the-art 2D/3D hybrid segmentation models while reducing costly preprocessing. Exploiting the rotating LiDAR structure motivates a novel architecture with high precision and significantly decreased latency and model size. We propose a simple alteration to recent deep surface reconstruction methods to incorporate visibility information. Our approach leads to across-the-board improvements and a significant increase in generalization capacity. Finally, we exploit the ability of LiDAR to penetrate tree canopy to learn the multi-layer structure of natural forests.
# Exploiting the Structure of
Satellite Time Series

This chapter presents our work on the automated analysis of Satellite Image Time Series, or SITS. Due to accessible data and plentiful annotations, SITS analysis is a prime example of a remote sensing task which benefits from the deep learning approach. However, SITS follow a complex multi-modal, spatial, temporal, and spectral structure. We propose a family of attention-based algorithms exploiting these properties to increase the speed and precision of crop mapping.

This chapter is organized around the following publications:

[GLGC20]: Vivien Sainte Fare Garnot, Loic Landrieu, Sebastien Giordano, Nesrine Chehata, "Satellite image time series classification with pixel-set encoders and temporal self-attention", CVPR, 2020

[GL21b]: Vivien Sainte Fare Garnot, Loic Landrieu, "Panoptic Segmentation of Satellite Image Time Series with Convolutional Temporal Attention Networks", ICCV, 2021

[GL20]: Vivien Sainte Fare Garnot, Loic Landrieu, "Lightweight Temporal Self-Attention for Classifying Satellite Image Time Series", ECML Workshop on Advanced Analysis and Learning on Temporal Data, 2020

[GL21a]: Vivien Sainte Fare Garnot, Loic Landrieu, "Leveraging Class Hierarchies with Metric-Guided Prototype Learning", BMVC, 2021

[GL22]: Vivien Sainte Fare Garnot, Loic Landrieu, "Multi-Modal Temporal Attention Models for Crop Mapping from Satellite Time Series", ISPRS Journal, 2021

[QL21]: Félix Quinton, Loic Landrieu, "Crop Rotation Modeling for Deep Learning-Based Parcel Classification from Satellite Time Series", Remote Sensing, 2021

[GBLC20]: Sébastien Giordano, Simon Bailly, Loic Landrieu, Nesrine Chehata, "Improved crop classification with rotation knowledge using Sentinel-1 and -2 time series", Photogrammetric Engineering & Remote Sensing, 2021

[GLGC19]: Vivien Sainte Fare Garnot, Loic Landrieu, Sebastien Giordano, Nesrine Chehata, "Time-space trade-off in deep learning models for crop classification on satellite multi-spectral image time series", IGARSS, 2020
# 1 Stakes and Challenges of SITS Analysis

<figure><img src="image_46.png" /><figcaption>Figure 5.1: **Sentinel-2 Image.** RGB values for the 31TFM tile, covering 110 × 110 km and containing over 200 000 agricultural parcels.</figcaption></figure>

Advances in space-based remote sensing, such as the launch of the Planet [BMKS14] and the open-access Sentinel constellations [DDBC$^{+12}$], have allowed for sustained improvements in the precision and availability of Earth observation. In particular, satellites with high revisit frequency are ideal for analyzing phenomena with complex temporal dynamics [GLGC19] such as crop mapping—the driving application of this chapter, see Figure 5.1.

The automated analysis of Earth observation creates impactful financial and environmental opportunities for private and public actors. Crop monitoring is necessary for the fair allocation of farming subsidies (57 and 22 billion euros per year in Europe and the US, respectively [CAP]). It can help enforce crop rotation best practices for maximizing yields [KKN$^{+15}$, SL08] and soil protection [AAM$^{+11}$, Bul92]. Automated SITS analysis is also used for other applications such as surveying urban development [TL16] and deforestation [RHV$^{+18}$].

Despite the inherent difficulty of differentiating between the complex growth patterns of cultivated plants, virtually limitless access to data and annotations have encouraged the development of algorithmic solutions for automated crop monitoring from satellites [KDB$^{+19}$]. Indeed, Sentinel-2 averages one multi-spectral observation every five days, which is beneficial for characterizing crop phenology [SBAK20, VMD$^{+18}$], see Figure 5.2. Moreover, farmers declare the kind of crop cultivated in their parcels yearly. This represents over 10 million yearly annotations for France alone [RPG], all openly accessible in the Land-Parcel Identification System.

SITS present a set of unique challenges and structures:

* **Large Scale.** The amount of data acquired by Sentinel-2 is considerable: more than 25 To / year in Europe for the optical modality alone. Training and applying deep models require computation- and
memory-efficient architectures.

* **Absolute Spatio-Temporal Frame of Reference.** The position of a pixel in a picture or the frame number in a video is typically considered arbitrary and uninformative information. In contrast, absolute spatial and temporal coordinates are crucial information in SITS as they significantly influence observations: wheat crops have different spectral signatures in Denmark or Italy, and in October or March.

* **Cloud Occlusion.** Agricultural parcels generally do not occlude each other from the near-nadir perspective of a satellite.<fnref n="1" /> However, cloud cover impacts acquisitions in several ways, from total occlusion to direct and indirect (shadow) radiometric corruption [LSL<sup>+</sup>17, GYL<sup>+</sup>20].

* **Low Spatio-Temporal Resolution.** The spatial resolution and revisit time of the Sentinel-2 satellites are nothing short of a technological marvel and revolutionized Earth observation. However, agricultural parcels are structured by elements such as hedges and furrows, which are typically smaller than the size of Sentinel-2 pixel (10m). Likewise, the evolution of agricultural land is organized around technical acts such as harvest or mowing, which are typically shorter than the satellite's revisit time of five days. This differs from the analysis of natural images or videos, in which the target objects or events are typically covered by multiple pixels and frames.

* **Multi-Modality.** Earth's surface is constantly monitored from space by an array of public, private, and military sensors that differ in nature and resolution. The Sentinel constellation itself comprises synthetic-aperture radar (S1), optical, microwave, and infrared radiometers (S2,S3), spectrometers (S3,S4,S5), and altimeters (S3,S6). In particular, C-band radar and visible/near-infrared radiometers collect complementary information [VTGGP18], which is useful for crop and land mapping. Furthermore, while cloud cover is highly disruptive for optical imagery, radar waves are essentially unaffected [SWNW18].

These characteristics make the analysis of SITS for crop monitoring an exciting and challenging machine learning task, deserving of uniquely tailored solutions. Throughout this chapter, we define SITS as 4-dimensional tensors of size $T \times C \times H \times W$, with $T$ the length of the sequence, $C$ the number of channels per pixel, and $H, W$ the spatial extent of the images. Note that in the general case, the acquisition dates depend on the region

<footnote n="1">¹Note that images with oblique viewing conditions also exist, causing self-occlusions in mountainous areas.</footnote>

<footnote n="2">²Clouds are surprisingly pervasive: up to 65% of the Earth's surface is occluded from space at any given time, and around 35% for land masses [JR08]. This makes crop monitoring [ESR<sup>+</sup>16] and emergency mapping [RRW19] particularly difficult in humid tropical regions.</footnote>
<figure><img src="image_47.png" /><figcaption>Figure 5.2: **Sentinel-2 Time Series.** Images (RGB bands, 10m per pixel) for parcels of *Winter cereal* and *Spring cereal*. The dots on the horizontal axis represent the unevenly distributed acquisition dates throughout the period of interest. Note the importance of the temporal evolution of the parcels to discriminate between classes.</figcaption></figure>

(tile) considered, a complication that our proposed solutions handle well.

## 2 Parcel Classification

In this setting, we know the exact extent of each parcel and try to predict the nature of the crop from a time sequence of observations. This is the operational setting in countries with an open-access Land Parcel Identification System (LPIS), such as France or Denmark.

### 2.1 Spatial Encoding

<figure><img src="image_48.png" /><figcaption>Figure 5.3: **Architecture of the PSE.** Each image of the considered sequence is processed as an unordered set of pixels.</figcaption></figure>
The spatial resolution of 10m per pixel of the Sentinel-2 optical sensor limits the nature of the spatial information that can be captured. In particular, textural information, such as furrows and hedges, is typically lost. This coarseness questions the widespread use of CNNs [KLSS17, RK18, GLGC19] as spatial encoders for SITS. In fact, CNNs are known to rely mainly on texture to extract spatial features [GRM+19].

Instead, we propose to consider each image as an unordered set of pixel, i.e. a point cloud in spectral space. Our proposed *Pixel-Set Encoder* (PSE) network randomly samples pixels for each date, and uses a modified set-based architecture [QSMG17a, ZKR+17] to learn statistical descriptors of the parcels' spectral distribution. While this approach destroys spatial information, the randomness of sampling acts as a powerful augmentation strategy to fight overfitting. We also incorporate a set of handcrafted spatial features (e.g. perimeter, surface), which are easy to compute and allow the network to retain high-level information about the parcel shape. Another advantage of this approach compared to image-based approaches is that parcels can be processed in parallel, regardless of their size. This allows us to forego the memory-intensive and information-altering resizing and padding steps required to embed images of various shapes with CNNs.

> Recognizing that the spatial resolution of Sentinel-2 images may be too low to represent texture, we view parcels as unordered sets of pixels and learn statistical descriptors of the parcels' spectra with a simple set-based network.

## 2.2 Attention-Based Temporal Encoding

The temporal dimension of remote sensing time series can be handled in a number of ways, such as temporal concatenation [KLG<sup>+</sup>16], temporal statistics [PVI<sup>+</sup>16], histograms [BMT<sup>+</sup>15], time kernels [TMC<sup>+</sup>17], shapelets [YK09], or probabilistic graphical models [GBLC20]. However, attention-based networks [VSP<sup>+</sup>17] have initiated a new era for sequential information analysis. Initially designed for Natural-Language Processing (NP), Transformers have proved more expressive and faster to train than RNNs for various tasks. We propose to adapt this concept to satellite image time series. Given the difference between the analysis of NLP and SITS, this requires several adjustments.

**Transformer Network.** In the original formulation, a *query-key-value* triplet $q^{(t)}, k^{(t)}, v^{(t)}$ is computed simultaneously for each element $x^{(t)}$ of the input sequence. The key $k^{(t)}$ conveys information about the nature of the content, the value $v^{(t)}$ encodes the content itself, and the query $q^{(t)}$ the type
<figure><img src="image_49.png" /><figcaption>Figure 5.4: **Architecture of the Temporal Attention Encoder.** Our simplified transformer architecture can efficiently embed sequences of image embeddings. Variables in bold are tensors concatenated along the temporal dimension, e.g. $\mathbf{e} = [e^{(0)}, \dots, e^{(T)}]$.</figcaption></figure>

of contextual information needed to embed $x^{(t)}$. Each element of the sequence associates an attention score to the preceding elements by computing the dot product compatibility between $k^{(t<T)}$ and $q^{(t)}$, rescaled with a modified softmax layer.

This procedure can be computed several times in parallel with different sets of independent parameters, or *heads*. This approach, called *multi-head attention*, allows for the specialization of different query-key compatibility. Since all computations are parallel, the Transformer takes full advantage of modern GPU architecture and yields a significant speed increase compared to recurrent architectures.

**Positional Encoding.** In their paper on text translation, Vaswani et. al encode the position of elements in the sequence with discrete Fourier coefficients, as defined in Equation 5.1. Instead of the position in the sequence, we use the number of days since the first observation to help to account for inconsistent temporal sampling (see Figure 5.2).

$$[p^{(t)}]_{i=1}^{E} = \sin \left( \frac{\text{day}(t)}{\tau^{2i}} + \frac{\pi}{2} \mod(i, 2) \right), \quad (5.1)$$

with $E$ the size of the spatial embedding.

**End-to-End Encoding.** Since we train our spatial and temporal encoders simultaneously, we can directly use spatial features as values. This removes
needless computation and avoids a potential information bottleneck. Using a single linear layer $FC_1^{(h)}$, we generate the keys $k_h^{(t)}$ and queries $k_h^{(t)}$ for each element $t$ of the sequence and each head $h$ in parallel :

$$k_h^{(t)}, q_h^{(t)} = FC_1^{(h)} (e^{(t)} + p^{(t)}) \qquad (5.2)$$

**Master Query Generation.** Our goal is to encode the sequence corresponding to one year of observation for each parcel. In NLP, this corresponds to the sentence classification setting, e.g. for sentiment analysis, and is typically done using a special “beginning of sentence” token [DCLT18]. Here, we generate a single master-query $\hat{q}_h$ per sequence and head $h$ as the temporal average of all queries of the sequence, processed by a single fully-connected layer $FC_2^{(h)}$:

$$\hat{q}_h = FC_2^{(h)} \left( \text{mean} \left( \{q_h^{(t)}\}_{t=1}^T \right) \right) \qquad (5.3)$$

Two benefits of this approach are that the query can adapt to the content of the sequence and is linear in time instead of quadratic for token embedding tasks.

**Multi-Head Temporal Attention.** We multiply the master query with the keys to produce attention weights $a^{(h)} \in [0,1]^T$ for each element of the sequence, determining which dates contain the most helpful information. We use the attention weight to obtain the output $o_h$ as the temporal average of the input embeddings. Finally, we concatenate the output of all heads and process the results with MLP3:

$$a_h = \text{softmax} \left( \frac{1}{\sqrt{d_k}} \left[ \hat{q}_h \cdot k_h^{(t)} \right]_{t=1}^T \right) \qquad (5.4)$$

$$o_h = \sum_{t=1}^{T} a_h[t] (e^{(t)} + p^{(t)}) \qquad (5.5)$$

$$\hat{o} = \text{MLP}_3 ([o_1, \dots, o_H]) . \qquad (5.6)$$

In Figure 5.5, we illustrate head specialization by plotting the average attention masks for two types of cereals. We can see that each of the four heads specializes in a specific portion of the time series. Additionally, the dependency of the attention mask to the input is apparent for head 4, which focuses on late spring for *Spring Cereal* samples and late summer for *Summer Cereal* samples.

## 2.3 Efficient Temporal Attention

An important observation from the previously presented TAE is that increasing the number of heads is beneficial for precision, but incurs a high
<figure><img src="image_50.png" /><figcaption>Figure 5.5: Average attention masks of the TAE heads, obtained from 128 samples of spring (a), and summer (b) cereal parcels. Head 4 adapts to the different classes to detect specific events.</figcaption></figure>

computational and memory cost. We thus propose the Lightweight Temporal Attention Encoder (L-TAE), a stripped-down version of the TAE allowing for the efficient use of many attention heads. We present the main proposed changes below and in Figure 5.6.

**Channel Grouping:** we split the *E* channels of the sequence of spatial embeddings into *H* groups of size *E'* = *E*/*H* with *H* the number of heads. We denote by $e_h^{(t)}$ the groups of input channels for the *h*-th head of the *t*-th element of the input sequence.

**Query-as-Parameter:** We save computations by directly defining the master queries as trainable parameters of the network. While such queries are not adaptive to the input, we can use many heads to maximize the expressivity of the learned features. Only the keys are obtained with a learned linear layer.
<figure><img src="image_51.png" /><figcaption>**Figure 5.6: L-TAE Architecture.** An input sequence **e** of *T* vectors of size *E*, with *H* = 3 heads and keys of size *K*. The channels of the input embeddings are distributed among heads. Each head uses a learned query $\hat{q}_h$, while a linear layer FC<sub>*h*</sub> maps inputs to keys. We concatenate the outputs of all heads into a vector the same size as the input embeddings.</figcaption></figure>

<table><caption>**Table 5.1: Asymptotic Complexity of Sequence Embedding Networks.**
The GRU's memory update complexity is given in the Keys and Mask columns. *X* is the size of the output vector. *M* is the size of the hidden state of the GRU.</caption><thead><tr><th>Method</th><th>Keys</th><th>Mask</th><th>Output</th></tr></thead><tbody><tr><td>L-TAE</td><td>O(TEK)</td><td>O(HTK)</td><td>O(EX)</td></tr><tr><td>TAE</td><td>O(HTEK)</td><td>O(HTK)</td><td>O(HEX)</td></tr><tr><td>Transf.</td><td>O(HTEK)</td><td>O(HT<sup>2</sup>K)</td><td>O(HEX)</td></tr><tr><td>GRU</td><td>O(MT(E + M))</td><td></td><td>O(MX)</td></tr></tbody></table>

**Multi-Head Attention:** The rest of the network operates similarly to the TAE: outputs $o_h$ of each head are defined as the temporal average of inputs weighted by attention scores derived from key-query compatibility. Finally, the outputs of all heads of size $E/H$ are concatenated into a vector of size *E* and processed by a multi-layer perceptron MLP to the desired width.

**Computational Complexity** In Table 5.1, we report the asymptotic complexity of different sequence embedding algorithms. The L-TAE channel grouping strategy removes the influence of *H* in the computation of keys and outputs compared to a TAE or a Transformer. The complexity of the L-TAE is also lower than the GRU's as *M*, the size of the hidden state, is typically larger than *K* (130 vs 8 in the experiments presented in Table 5.2).
## 2.4 Quantitative Results

<table><caption>Table 5.2: **Quantitative Evaluation.** We report the overall Accuracy (OA) and mIoU Classification metrics for different architectures, as well as training time (one epoch) and inference time for the entire dataset, and the size of the dataset on the disk. ¹ disk space required for training and pure inference, ² time for the entire training step, <fnref n="3" /> preprocessing and inference time, ⁴ dataset before and after preprocessing.</caption><thead><tr><th></th><th>OA</th><th>mIoU</th><th>Training (s/epoch)</th><th>Inference (s/dataset)</th><th>Disk Size Gb</th></tr></thead><tbody><tr><td>PSE+TAE (ours)</td><td><strong>94.2</strong> ±0.1</td><td><strong>50.9</strong> ±0.8</td><td>158</td><td><strong>149</strong></td><td><strong>28.6 / 12.3<sup>1</sup></strong></td></tr><tr><td>CNN+GRU [GLGC19]</td><td>93.8 ±0.3</td><td>48.1 ±0.6</td><td>656</td><td>633</td><td>98.1</td></tr><tr><td>CNN+TCNN [PWP19]</td><td>93.3 ±0.2</td><td>47.5 ±1.0</td><td>635</td><td>608</td><td>98.1</td></tr><tr><td>Transformer [RK19]</td><td>92.2 ±0.3</td><td>42.8 ±1.1</td><td>13</td><td>420 + 4<sup>3</sup></td><td><strong>28.6 / 0.22<sup>4</sup></strong></td></tr><tr><td>ConvLSTM [RK18]</td><td>92.5 ±0.5</td><td>42.1 ±1.2</td><td>1283</td><td>666</td><td>98.1</td></tr><tr><td>Random Forest [GBLC20]</td><td>91.6 ±1.7</td><td>32.5 ±1.4</td><td><strong>293<sup>2</sup></strong></td><td>420 + 4<sup>3</sup></td><td><strong>28.6 / 0.44<sup>4</sup></strong></td></tr></tbody></table>

We report in Table 5.2 the performance of our approach and different competing networks on a proposed open-access dataset.

**Sentinel2-Agri.** We propose a dataset comprising 191 703 image time sequences corresponding to agricultural parcels of the T31TFM Tile in southern France. All sequences are composed of 24 observations spanning from January to October 2017, with a spatial resolution of 10m and 10 spectral bands (we discard B01, B09, and B10). We associate each sequence with the label corresponding to the majority culture retrieved from the French Land Parcel Identification System records<fnref n="3" /> with a 20 class nomenclature designed by the subsidy allocation authority of France. The dataset is highly imbalanced: four classes cover 90% of the samples.

We propose two versions of the dataset, depending if images are stored as patches or sets. In the *patch format*, we resize each parcel into a tensor of size $T \times C \times 32 \times 32$ with nearest neighbour spatial interpolation and zero-padding. In the *set format*, the pixels of each parcel are stored in arbitrary order into a tensor of size $T \times C \times N$, with $N$ the total number of pixels in a given parcel. Note that this format neither loses nor creates information, regardless of parcel size. Hence, this setup saves up to 70% disk space compared to the *patch format* (28.6Gb vs. 98.1Gb). This dataset constitutes the first large-scale dataset for object-based agricultural parcel classification and is freely accessible at github.com/VSainteuf/pytorch-psetae.

<footnote n="3">³http://professionnels.ign.fr/rpg</footnote>
<figure><img src="image_52.png" /><figcaption>**Figure 5.7: Precision, Speed, and Model Size.** We report the mIoU averaged over 5 runs for different temporal encoders plotted with respect to the number of FLOPs necessary to process one sequence. The marker size represents the number of trainable parameters in the model. The L-TAE outperforms other models across all model sizes and processing requirements.</figcaption></figure>

**Competing Methods.** We compare our approach to recent algorithms operating on similar datasets. These methods use CNNs or spatial pooling for spatial encoding and recurrent neural networks, temporal convolutions, or transformers for temporal encoding. We parametrize all models to have roughly 150k parameters. We report in Table 5.2 the performance of all methods in terms of precision, speed, and memory usage.

**Efficient Temporal Attention.** We focus on the advantages of the L-TAE by using the same PSE network for spatial encoding and varying the temporal encoders. We report the performance in Figure 5.7.

We propose to adapt the successful attention model from NLP to SITS. Not only does our approach define a new state-of-the-art in terms of precision, but it is four times faster and uses a parsimonious data format only requiring a quarter of disk size. A variation of our temporal model reaches even higher performance with smaller models: a model with only 9k parameters overtakes a GRU with over 3m parameters. Several other studies confirmed the superiority of our models for SITS [SK20, KTR+21].
# 3 Panoptic Segmentation of SITS

We now study crop mapping when both the content and contour of parcels must be retrieved, which is the most common operational setting. This problem can be framed as panoptic segmentation of an image sequence [KGHD19, MV21]. However, since agricultural parcels are geo-referenced and static, we only need to consider the geo-referenced pixels and not all images in the sequences. This removes the complex task of spatially tracking objects [TSA19]. While some approaches propose to perform instance segmentation [Rie17], delineation (border detection) [GPGMLS17, MPT20, WD20], or oversegmentation [PKPG12], there exists no dedicated approach for detecting individual objects consistently across an entire satellite image sequence.

## 3.1 Spatio-Temporal Pixel Encoding

The first step of our pipeline is to learn pixel-wise features encoding the spatio-temporal dynamics of SITS containing multiple parcels. The key insight of our proposed architecture is to use temporal attention to compute meaningful spatio-temporal maps at several resolution levels *simultaneously*. Temporal U-Net networks [SPI+19, RCW+19, PVK21] only process the temporal dimension of the lowest resolution and collapse the skip connections with spatial averaging. This prevents the extraction of spatially adaptive and parcel-specific temporal patterns at higher resolutions. Conversely, convolution-recurrent encoders [RK18] process the temporal dimension at the highest resolution, which results in an increased memory requirement and ignores the low spatial regularity of the data. Instead, we propose to efficiently compute temporal attention masks at all resolutions using an upsampling strategy on the attention masks themselves, see Figure 5.9.

Our model, dubbed U-TAE (U-Net with Temporal Attention Encoder), encodes a sequence $x$ of dimension $T \times C \times H \times W$ in three steps:

(a) Each image in the sequence is embedded simultaneously and independently by a shared multi-level spatial convolutional encoder.

(b) We apply an L-TAE to all the pixels of the lowest resolution feature map independently to obtain pixelwise temporal attention masks. We then use bilinear interpolation to obtain masks at all spatial resolutions and collapse the spatial dimension of the maps at all levels.

(c) The resulting spatial feature maps are processed by a convolutional decoder to produce rich feature maps at all levels.

A subtlety of our spatial encoder is that the image batches mix acquisitions taken at different dates. Consequently, the samples are not identically distributed, which prevents the use of BatchNorms [IS15]. To address this issue, we use Group Normalization [WH18] with 4 groups instead in the
<figure><img src="image_53.png" /><figcaption>Figure 5.8: **Problem Statement.** We consider an image sequence (a) with panoptic annotations (b). Our objective is to learn spatio-temporal pixel features, which can be used to perform panoptic (c) and semantic segmentation (d). We observe some difficult configurations: boundary ambiguities that can not be resolved from a single image (cyan circle ○), and conversely, visually fragmented areas annotated as a single instance (red circle ○).</figcaption></figure>

encoder. As shown later, the impact of this simple change on the performance is drastic for all *temporal U-Net* and not only ours.

**PASTIS Dataset.** We introduce the PASTIS (Panoptic Agricultural Satellite Time Series) dataset, the first large-scale, publicly available SITS dataset with semantic and panoptic annotations. PASTIS is comprised of 2433 sequences of multi-spectral images of size 10 × 128 × 128. Each sequence contains between 38 and 61 observations taken between September 2018 and November 2019, for over 2 billion pixels spanning 4000 km². We estimate that close to 28% of images have at least partial cloud cover.

Each pixel of PASTIS is associated with a semantic label taken from a nomenclature of 18 crop types plus a background class. Each non-
<figure><img src="image_54.png" /><figcaption>Figure 5.9: **Spatio-temporal Encoding.** A sequence of images is processed in parallel by a shared convolutional encoder. At the lowest resolution, an attention-based temporal encoder produces a set of temporal attention masks for each pixel, which are then spatially interpolated at all resolutions. These masks are then used to collapse the temporal dimension of the feature map sequences into a single map per resolution. Finally, a convolutional decoder computes features at all resolution levels.</figcaption></figure>

background pixel has a unique instance label corresponding to its parcel index. In total, 124 422 parcels are available, each with bounding box, pixel-precise mask, and crop type. The French Payment Agency estimates the accuracy of the crop annotations via in situ control to be over 98% and the relative error in surface to be under 0.3%. To allow for fair cross-validation, we split PASTIS into 5 folds with a 1km buffer between images.

**Semantic Segmentation.** By setting the width of the highest resolution layer to the number of classes, the U-TAE can perform semantic segmentation. As seen in Table 5.3 the U-TAE significantly outperforms the other methods from the state-of-the-art. Interestingly, most of the advantage of our approach is attributable to the use of attention maps at all resolution levels. Collapsing the skip connections with simple temporal average results in a performance comparable to the other evaluated methods. The main advantage of temporal attention is that it can be easily used at different resolutions, in contrast to recurrent network-based approaches.

## 3.2 Panoptic Segmentation

We now use the multi-scale feature maps learned with the U-TAE to perform panoptic segmentation. Since the U-TAE operates on many images for each sequence (often over 50), we favour an efficient design. Furthermore, given the relative simplicity of parcel borders, we can avoid com-
<table><caption>Table 5.3: **Semantic Segmentation.** We report for our method and six competing methods the model size in trainable parameters, Overall Accuracy (OA), mean Intersection over Union (mIoU), and Inference Time for one fold of ~ 490 sequences (IT).</caption><thead><tr><th>Model</th><th># param<br>×1000</th><th>OA</th><th>mIoU</th><th>IT (s)</th></tr></thead><tbody><tr><td>U-TAE (ours)</td><td>1 087</td><td><strong>83.2</strong></td><td><strong>63.1</strong></td><td><strong>25.7</strong></td></tr><tr><td>3D-Unet [RCW<sup>+</sup>19]</td><td>1 554</td><td>81.3</td><td>58.4</td><td>29.5</td></tr><tr><td>U-ConvLSTM [RCW<sup>+</sup>19]</td><td>1 508</td><td>82.1</td><td>57.8</td><td>28.3</td></tr><tr><td>FPN-ConvLSTM [MLRF<sup>+</sup>21]</td><td>1 261</td><td>81.6</td><td>57.1</td><td>103.6</td></tr><tr><td>U-BiConvLSTM [MLRF<sup>+</sup>21]</td><td>1 434</td><td>81.8</td><td>55.9</td><td>32.7</td></tr><tr><td>ConvGRU [BYPC16]</td><td>1 040</td><td>79.8</td><td>54.2</td><td>49.0</td></tr><tr><td>ConvLSTM [RK18]</td><td>1 010</td><td>77.9</td><td>49.1</td><td>49.1</td></tr><tr><td>U-TAE w. Skip Mean</td><td>1 074</td><td>82.0</td><td>58.3</td><td>24.5</td></tr></tbody></table>

plex region proposal networks such as Mask-RCNN [HGDG17]. Instead, we adapt the single-stage Object-as-Point instance segmentation network [ZWK19, WXS<sup>+</sup>20] and name our network *Parcels-as-Points* (PaPs) to highlight this inspiration.

**Parcel Detection.** We define a *centerness heatmap* from the ground truth parcel bounding boxes by associating each centroid with a heteroscedastic Gaussian kernel depending on their spatial extent. This allows us to frame the detection problem as the regression of this heatmap, supervised using a modified logistic loss. We predict parcel centerpoints at all spatial local maxima of the predicted centerness heatmap, which can be efficiently retrieved with a single max-pooling operation, see Figure 5.10.

**Size, Class, and Shape Prediction.** We attribute to each detected centerpoint a multi-scale descriptor by concatenating the feature maps at all resolutions in the channel dimension. From this descriptor, we predict a class, a size, and a small shape patch of fixed size *S*, typically 16 pixels. This patch is reshaped using the predicted size to obtain a rough shape prediction. To obtain a pixel-precise instance prediction, we also predict a saliency map at full resolution, shared for all detected centerpoints of the sequence. This saliency map is cropped along the predicted parcel bounding box, and combined with the rough shape prediction in a residual fashion to obtain a soft occupancy mask for each predicted centerpoint, see Figure 5.11

We supervise the class predictions with the cross-entropy and the predicted sizes with a rescaled L1 loss. The shape patches and saliency map are supervised with the binary cross entropy between the predicted and
<figure><img src="image_55.png" /><figcaption>Figure 5.10: **Centerpoint Detection.** The ground truth instance masks (a) is used to construct a target heatmap (b). Our parcel detection module maps the observation (c) to a predicted heatmap (d). The predicted centerpoints (red crosses) are the local maxima of the predicted heatmap (d). The black dots are the true parcel centers.</figcaption></figure>