The latest generation of aerial-and satellite-based imaging sensors acquires huge volumes of Earth's images with high spatial, spectral and temporal resolution, which open the door to a large range of important applications, such as the monitoring of natural disasters, the planning of urban environments and precision agriculture. In order to fully exploit the potential offered by these sensors, there is a need to develop accurate and time-efficient mathematical models and algorithms for spectral-spatial analysis of the recorded highresolution data.

The main goal of my research is to develop learning approaches, which would help to automatically interpret, or classify, remote sensing images. This manuscript presents several strategies I have explored for this purpose, varying from the use of strong shape priors to detect objects, regularization of classification probabilities on the image graphs, and up to the use of convolutional neural network models capable to learn deep hierarchical contextual features.

The experimental results on diverse benchmarks of images and image time series show the competitiveness of the developed methods when compared to the state-of-the-art approaches. In particular, we have recently created large-scale classification benchmark of aerial images and have demonstrated that the modern deep learning-based methods succeed in generalizing to the dissimilar urban settlements around the Earth. This opens new exciting perspectives towards designing systems which would be able to automatically update world-scale maps from remote sensing data.

Chapter 1

Introduction

Motivations and challenges

The continuous proliferation and improvement of remote (satellite and aerial) data sensors yields a huge volume of Earth's images with high spatial and temporal resolution, as well as with rich spectral information. For instance, the constellation of Pléiades satellites revisits the entire Earth every day, acquiring optical data with a spatial resolution of less than a meter. These data contain a valuable source of information, which opens the door to a large range of important applications, such as the monitoring of natural disasters, the planning of urban environments and precision agriculture. However, petabytes of these massive images are stored into binary files as unstructured raw data. As a result, a large part of them is never used. Thus, it has become crucial to develop methods to automatically analyze these data. ing data is image classification, which consists in the assignment of a class, or thematic label (for example, building, tree, car, etc.), to each pixel in the image. In this work, I focus on supervised classification, in which the learning approaches are used to train the classifier from labeled reference data. In the manuscript, the terms labeling and classification are used interchangeably. I also refer to dense classification or semantic classification to denote the labeling of individual pixels in the image. The term semantic is typically used to emphasize that the classes of interest represent object types (e.g., car, building [START_REF] Volpi | Dense semantic labeling of sub-decimeter resolution images with convolutional neural networks[END_REF]) rather than low-level physical materials (e.g., metal sheets [13]).

While the first classification approaches have been designed more than forty years ago with the appearance of the first satellites [START_REF] Kettig | Classification of multispectral image data by extraction and classification of homogeneous objects[END_REF], recent technological developments raise new challenges. The first challenge consists in adapting classification methods to the characteristics of new sensors, since the state-of-the-art methods may be not directly applicable or may not give satisfactory performances. For instance, first satellite sensors captured data of a low spatial resolution, and significant research was focused on discovering the proportions (or abundances) of physical materials within each single pixel [START_REF] Keshava | A survey of spectral unmixing algorithms[END_REF]. The recently developed sensors acquire images with a high spatial resolution, where pixels become small when compared to the objects. This implies the need of object-based analysis, i.e., the development of the algorithms which learn to analyze the consistency and relationship between groups of pixels in the image, aiming at identifying the whole objects.

The second challenge is the scalability of the systems. With the increasing amount and openness of the data covering a large geographic extent, there is a growing interest in developing scalable methods able to produce classification maps of large areas, in the scale of cities, countries, continents and up to the entire Earth [START_REF] Mattyus | Enhancing road maps by parsing aerial images around the world[END_REF]. While there is an urgent need to design computationally efficient algorithms to process large-scale data, it is also a challenge to design methods that would produce accurate maps for dissimilar areas of the Earth, given the large variability of the appearance of landscapes across the planet, i.e. large intra-class variability (see Fig. 1.2). While the first classification methods in remote sensing tried to distinguish different classes of low-level materials and substances (e.g., water, soil, gravel ), current research studies are in a big part focused on semantic segmentation, where high-level contextual reasoning (e.g., shape and context analysis) must be performed to discriminate more abstract semantic classes.

In this manuscript, I present my research studies on the development of learning methods for remote sensing image classification. In the following, I first summarize the state-of-the-art approaches, and then detail my recent contributions.

Approaches, trends and contributions

The pioneering classification approaches for remote sensing imagery consider each pixel as a pattern and all the data collected by sensors for this pixel are considered as the initial set of features (for instance, spectrum values in the case of optical imagery). Then, the standard pattern recognition scheme [START_REF] Duda | Pattern Classification[END_REF] is applied to classify each pixel:

1. The classification model and its specifications are selected (e.g., neural network, support vector machines (SVM), etc.).

2. Since the initial set of features may be redundant, feature extraction/selection is often performed, aiming at obtaining the features which discriminate at best classes of interest.

3. The next step, called "training of the classifier," consists in partitioning the entire feature space into K exhaustive, non-overlapping regions, so that every point in the feature space is uniquely associated with one of the K classes. This is done by feeding into the classification model examples of samples for each of the classes (so called training samples).

Once the training step is accomplished, every new pixel can be classified according to its feature set. Such classification approach is often referred to as pixelwise, where the term pixelwise has the different meaning than in the previous section: here, by pixelwise we mean that to classify any pixel, no information from its neighborhood is taken into account. Different models have been applied for pixelwise classification in the remote sensing community; among them, probably the most widely used ones are neural networks [7,8,[START_REF] Merényi | Intelligent understanding of hyperspectral images through selforganizing neural maps[END_REF], SVM [START_REF] Vapnik | Statistical Learning Theory[END_REF]21] and random forests [START_REF] Ham | Investigation of the random forest framework for classification of hyperspectral data[END_REF][START_REF] Hansch | When to fuse what? random forest based fusion of low-, mid-, and high-level information for land cover classification from optical and sar images[END_REF].

In a real image, neighboring pixels are related or correlated, both because imaging sensors acquire significant amount of energy from adjacent pixels [START_REF] Richards | Remote Sensing Digital Image Analysis[END_REF] and because objects occurred in the image scene are in most cases large when compared to the size of a pixel. For instance, if a given pixel in the agricultural image scene represents a corn, its adjacent pixels belong with the high probability to the corn field. This spatial contextual information significantly helps to interpret image scene. Therefore, most recent classification approaches are spectral-spatial, i.e. they classify each image pixel based on:

1) its own spectral values (the spectral information) and 2) information extracted from its neighborhood (the spatial information). The use of spectral-spatial classifiers is especially rewarding when processing images with high spatial resolution; that is why these techniques are especially relevant to analyze data captured by modern and near-future sensors.

I refer the reader to my PhD thesis [START_REF] Tarabalka | Classification of Hyperspectral Data Using Spectral-Spatial Approaches[END_REF] for a detailed overview of pixelwise and spectral-spatial classification approaches. There are different ways to group these techniques in the state-of-the-art. In the following, I will present three important families/strategies of spectral-spatial classification, which are successfully used until nowadays and can also be used in combination with each other, and will position my contributions with respect to these strategies:

1. Graph-based methods. In these approaches, an image is represented as a graph, where the nodes typically represent image pixels or regions, while edges either connect adjacent regions [START_REF] Tarabalka | Segmentation and classification of hyperspectral images using minimum spanning forest grown from automatically selected markers[END_REF]10], or keep track of an hierarchy in a multi-scale image representation [START_REF] Valero | Hyperspectral image representation and processing with binary partition trees[END_REF][START_REF] Maggiori | Improved partition trees for multiclass segmentation of remote sensing images[END_REF]. By applying different graph construction and pruning strategies, these techniques typically output a disconnected graph, where each connected subgraph corresponds to the classified object in the resulting classification map. A popular graph-cut technique applies a min-cut algorithm on the image graph to regularize outputs of the pixelwise classifier, by minimizing the Markov Random Field-based energy function associated with the image and the initial classification result [START_REF] Tarabalka | Graph-cut-based model for spectral-spatial classification of hyperspectral images[END_REF]. Another approach consists in building first a hierarchy of image regions at different levels of details, based on standard homogeneity/dissimilarity criteria, and then selecting from this hierarchy the regions at different scales that correspond to the objects in the final classification map [START_REF] Valero | Hyperspectral image representation and processing with binary partition trees[END_REF][START_REF] Maggiori | Optimizing partition trees for multi-object segmentation with shape prior[END_REF].

Such hierarchical data structure can be very flexible, for instance, it can be easily augmented to include textural or shape information for more accurate image anal-ysis. I will present our graph-based methods in Chapters 3 and 4. Chapter 3 will detail a method for shape-aware multi-object multi-class classification, based on a hierarchical data structure known as binary partition trees. Chapter 4 will describe how a graph-cut-based method can be segmentation for the analysis of image time series.

2. Feature engineering. The spatial information can be included into a classifier by analyzing the similarity/correlations between the data in the neighborhood of an image pixel. The initial approaches would typically use a well-defined algorithm, which can be often described by a (set of) pre-defined convolutional mask(s), for this purpose. For instance, Tsai et al. [START_REF] Tsai | Texture analysis for three dimension remote sensing data by 3D GLCM[END_REF] and Huang and Zhang [START_REF] Huang | A comparative study of spatial approaches for urban mapping using hyperspectral rosis images over pavia city, northern italy[END_REF] have investigated the use of texture measures derived from the gray level co-occurrence matrix, such as angular second moment, entropy and homogeneity, for spectralspatial classification. A popular approach in remote sensing explores the principles of mathematical morphology, which by definition aims at analyzing spatial relations between sets of pixels, i.e., extracting information about the size, shape and orientation of structures. For this purpose, morphological profiles [9], extended morphological profiles [START_REF] Fauvel | Spectral and spatial classification of hyperspectral data using svms and morphological profiles[END_REF], attribute profiles [START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF] and extinction profiles [START_REF] Ghamisi | Extinction profiles for the classification of remote sensing data[END_REF] have been successfully applied. In all cases, the core idea is to progressively simplify image based on a pre-defined rule/attribute, which allows to analyze a particular characteristic (for instant, area or standard deviation of regions) within this image. The extracted morphological features typically lay in a very-high-dimensional space, and additional feature extraction must be applied to reduce the dimensionality of the feature space. Another group of techniques introduces a strong shape priors into classification: for example, in remote sensing images buildings can be modelled as rectangles [START_REF] Ortner | Building outline extraction from digital elevation models using marked point processes[END_REF][START_REF] Wang | Building detection from high-resolution polsar data at the rectangle level by combining region and edge information[END_REF] and roads can be modelled by line segments [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF]. Chapter 2 will describe our studies using such approach, where we model curvilinear structures by a set of line segments with variable length and orientations.

3. Deep learned features. Even though feature engineering gives satisfactory results for many applications, when the complexity of the data and the number of classes increases, it is not straightforward to design hand-crafted features that would exhibit high discriminative power for all classes. The recent trend in the state-of-the-art consists in designing classification systems that would be able to automatically learn hundreds of expressive high-level contextual features. In this context, the use of deep learning, in particular convolutional neural networks, is intensively investigated by both image analysis and remote sensing communities.

While some works use deep learning-based features as the input to a conventional classifier, such as SVM or logistic regression classifier [START_REF] Zhao | On combining multiscale deep learning features for the classification of hyperspectral remote sensing imagery[END_REF], most recent studies design an end-to-end system which jointly learns to extract relevant contextual features and conduct the classification [START_REF] Mnih | Machine learning for aerial image labeling[END_REF][START_REF] Maggiori | Convolutional neural networks for large-scale remote-sensing image classification[END_REF][START_REF] Marmanis | Semantic segmentation of aerial images with explicit class-boundary modeling[END_REF]. I will detail our recent contributions on the development of deep learning architectures for remote sensing image classification in Chapter 5. In particular, I will present how we face the challenge to obtain fine-grained high-resolution classification maps by designing appropriate deep learning architectures.

Even though the main focus of my research is remote sensing data classification, I will also show experimental results on other datasets, such as medical or computer vision benchmarks, to illustrate the genericity of the proposed methods.

Selected publications

This section lists my publications after the PhD thesis defense, i.e. after 2010.

Peer-reviewed international journals

1. E. Maggiori, Y. Tarabalka, G. Charpiat, and P. Alliez, "High-resolution aerial image labeling with convolutional neural networks," IEEE Trans. Geosc. Remote Sens., 2017.

2. E. Maggiori, G. Charpiat, Y. Tarabalka, and P. Alliez, "Recurrent neural networks to correct satellite image classification maps," IEEE Trans. Geosc. Remote Sens., 2017.

3. E. Maggiori, Y. Tarabalka, G. Charpiat, and P. Alliez, "Convolutional neural networks for large-scale remote sensing image classification," IEEE Trans. Geosc. Remote Sens., Volume 55 (2), pages 645-657, February 2017.

Reconstruction of Curvilinear Structures

One of the pioneering methods aiming at incorporating information from the spatial neighbourhood to properly classify images used a template-based concept. This implies introducing a strong shape prior into classification. For instance, in moderate-spatialresolution aerial images roads can be modelled by line segments [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF], and buildings can be seen as rectangles [START_REF] Ortner | Building outline extraction from digital elevation models using marked point processes[END_REF].

This chapter covers my research on reconstruction of curvilinear structures based on stochastic modeling and ranking learning system. We assumed that the entire line network can be decomposed into a set of line segments with variable lengths and orientations. This assumption enables us to reconstruct arbitrary shapes of curvilinear structure for different types of datasets.

In the following, I describe two models we have developed. The first model works in an unsupervised way, and is based on Marked point process (MPP) and Monte Carlo sampling. We reconstruct the latent curvilinear structure by sampling disjoint line segments which are associated with with the given image data. Therefore, we aim to maximize a posterior probability density with respect to a set of line segments for given image data.

We constrain local interaction of the line segments to obtain smoothly connected line configuration. The optimization technique consists of two steps to reduce the significance of the parameter selection in our stochastic model. We simulate several Markov chains with different parameter settings to collect line hypotheses on the same configuration space.

Then, we maximize the consensus among line hypotheses to reduce the sampling space and to improve the reliability of the curvilinear structure reconstruction.

The second model is based on supervised learning and graph theory. To extract local curvilinear features, we compute oriented gradient information using steerable filters. We then employ structured support vector machine (SVM) for ordinal regression of the input image patches, where the ordering is determined by shape similarity to latent curvilinear structure. Finally, we progressively reconstruct the curvilinear structure by looking for geodesic paths connecting remote vertices in the graph built on the structured output rankings.

Related work

Curvilinear structure extraction

In early vision, researchers design convolution filters to separate curvilinear structure from the background texture ( [START_REF] Freeman | The design and use of steerable filters[END_REF][START_REF] Perona | Deformable kernels for early vision[END_REF][START_REF] Jacob | Design of steerable filters for feature detection using Canny-like criteria[END_REF]). The main idea behind such filter design is to create simple line shape templates in order to extract features showing high gradient magnitudes with a locally consistent orientation. However, these image filtering responses are unable to discern linear structures from undesirable high-frequency components, e.g., noise and edges. Enhancement filtering (EF) algorithm proposed by [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] analyzes the eigenvalues of the Hessian matrix of the image to evaluate local tubularity score. Optimally Oriented Flux (OOF) measures the amount of gradient flow at the boundaries to find the continuous linear structure ( [START_REF] Law | Three dimensional curvilinear structure detection using optimally oriented flux[END_REF]). Morphological operator also highlights the curvilinear structure by collecting pixels according to the recursive structural similarity ( [START_REF] Talbot | Efficient complete and incomplete path openings and closings[END_REF]).

On the other hand, the graphical models exploited image-based evidence with geometric shape constraints to improve the detection performances. For instance, tree structure can be used to reinterpret complex line networks. The authors of [START_REF] González | Delineating trees in noisy 2D images and 3D image-stacks[END_REF] looked for a path between singular points corresponding to intersections of the latent curvilinear structure. Türetken et al. [START_REF] Türetken | Reconstructing loopy curvilinear structures using integer programming[END_REF] formulated a large linear programming problem to constrain the diverse cases of the local interaction on the line networks.

Furthermore, curves can be approximated by multiple straight line segments. Stochastic models specify the distribution of line objects with pairwise interaction terms ( [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF]).

Reversible Jump Markov chain Monte Carlo sampler proposed by [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF] has been involved to optimize the probability density. However, due to the sensitivity of the parameter setup, the stochastic models are less practical for a large amount of varied datasets.

In recent years, machine learning algorithms have been favoured for designing optimal filter banks to extract curvilinear features ( [START_REF] Rigamonti | Accurate and efficient linear structure segmentation by leveraging ad hoc features with learned filters[END_REF][START_REF] Becker | Supervised feature learning for curvilinear structure segmentation[END_REF]). While a threshold value should be set to reconstruct the centerline of the curvilinear structure from the filtering scores, it is difficult to find the value which would be suitable for different types of applications.

To address this issue, [START_REF] Sironi | Multiscale centerline detection by learning a scale-space distance transform[END_REF] proposed a regression model of distance transform to predict the scale of the curvilinear structure and localize centerlines.

Learning structured information

We want to extract the structured information of the line networks that appear within homogeneous background textures. Since the information embedded in a single pixel is limited to infer the latent structure, Markov Random Fields ( [START_REF] Geman | Stochastic relaxation, Gibbs distribution and the bayesian restoration of images[END_REF][START_REF] Besag | On the statistical analysis of dirty pictures[END_REF]) or Conditional Random Fields (CRF, [START_REF] Lafferty | Conditional random fields: Probabilistic models for segmetaing and labeling sequence data[END_REF][START_REF] Kumar | Discriminative random fields[END_REF]) based models have been developed to enforce the label consistency on the pre-organized output space. e.g., neighborhood pixels are assumed to have the same label with the probability when minimizing the cost function. However, the topology of the curvilinear structure composed of line objects is too intricate to be applicable with such design approaches. It is also impossible to specify the topology of the complex line networks with a few parameters.

Recently, machine learning algorithms have been involved to detect curvilinear structures latent in various types of images. Becker et al. [START_REF] Becker | Supervised feature learning for curvilinear structure segmentation[END_REF] applied a boosting algorithm to obtain an optimal set of convolution filter banks. Sironi et al. proposed a regression model to detect centerlines by learning the scale (width) of the tubular structures with the non-maximum suppression technique [START_REF] Sironi | Multiscale centerline detection by learning a scale-space distance transform[END_REF].

Structured learning systems have been employed in image segmentation models based on random fields [11,[START_REF] Lucchi | Structured image segmentation using kernelized features[END_REF][START_REF] Kim | Image segmentation using higherorder correlation clustering[END_REF]. More specifically, Structured Support Vector Machine (SSVM) [START_REF] Tsochantaridis | Large margin methods for structured and interdependent output variables[END_REF] has been used to predict model parameters for inference of the structured information between input image space and output label space. Exploring all possible combinations of the labels in the output space is computationally intractable. Thus, the random field models define the pairwise relationship in a neighborhood system to enforce the labeling consistency. While such prior models based on random fields are successful to describe convex shaped objects, they should be adapted to detect thin and elongated shaped curvilinear objects.

Stochastic modeling of curvilinear structure reconstruction

In the following, I describe the proposed stochastic model based on MPP framework for curvilinear feature extraction in a fully automatic way, where the performance is not biased by the hyperparameter selection.

Marked point process modeling

We briefly review the definition of MPP [START_REF] Van Lieshout | Markov point processes and their application[END_REF][START_REF] Stoyan | Stochastic geometry and its applications[END_REF] to provide a mathematical description of the proposed model. Definition 1 (Spatial point process). A realization of point process consists of an unordered set of points in a compact set F ⊂ R d . A point process on F maps from a measurable probability space (F, B, µ) onto the configuration space Ω = ∪ ∞ n=0 Ω n , where B denotes σ-algebra of subset of F, and µ is the Lebesgue measure. In other words, for all bounded Borel sets B ⊆ B, the number of points falling in B is a finite random variable.

Definition 2 (Marked point process). In the MPP framework, each point is associated with additional information which describe a shape of the object. Specifically, we reconstruct curvilinear structures as smoothly connected line segments. Let s i = (x i , m i ) be a line segment specifying its center point x i = (x i , y i ) in the image sites F with a label of the length and the orientation m i = ( i , θ i ), where the label is sampled from the mark space M with a probability measure µ M . We now define a marked point process on F × M as a finite random configuration s = {s 1 , . . . , s n } ∈ Ψ.

The probability distribution of the MPP is defined based on an image I and spatial interactions between line segments. Given an image, we look for an optimal configuration ŝ which maximizes the unnormalized probability density f (s) as follows:

ŝ = argmax s∈Ψ f (s) = argmin s∈Ψ #(s) i=1 U d (s i ) + i∼j U p (s i , s j ), (2.1)
where #(s) is the cardinality of the configuration, and i ∼ j represents the symmetry relationship between interacting line segments s i and s j . U d and U p denote the data energy (which we also call data likelihood ) and the prior energy, respectively. In general, Monte

Carlo samplers [START_REF] Gamal-Eldin | Multiple birth and cut algorithm for point process optimization[END_REF][START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF][START_REF] Green | Delayed rejection in reversible jump Metropolis-Hastings[END_REF][START_REF] Verdié | Detecting parametric objects in large scenes by Monte Carlo sampling[END_REF] are employed in MPP models to maximize the proposed density function f (s). Each state of a discrete Markov chain (X t ) t∈N corresponds to a random configuration on the Ψ. The chain is locally perturbed by transition kernels, and is evolved to converge to the stationary distribution which is identical to the proposed probability density.

Data energy

We define the data energy of the line segment s i as a weighted sum of the rotated gradient magnitudes U m d and the intensity variance U v d along the line:

U d (s i ) = ω m d U m d (s i ) + ω v d U v d (s i ), (2.2)
where ω m d and ω v d are weighting coefficients corresponding to U m d and U v d , respectively. We obtain the rotated gradient information by convolving the input image with steerable filters [START_REF] Freeman | The design and use of steerable filters[END_REF][START_REF] Jacob | Design of steerable filters for feature detection using Canny-like criteria[END_REF]. Steerable filters are generated from a linear combination of basis filters. In this work, we use second-order derivatives of an isotropic Gaussian function as the basis filters. Let g θ i x; σ 2 be a steerable filter associated with an orientation θ i and ∇I θ i = g θ i * I be its filtering response, which adaptively accentuates gradient magnitudes corresponding to the angle θ i . Then, the gradient magnitude energy U m d is defined as

U m d (s i ) = 1 0 ∇I θ i (p i (t)) dt, (2.3)
where p i (t) represents points on the line segment s i . Note that

p i (t) = (1 -t)u i + tv i is a function of the endpoints u i and v i with parameter t ∈ [0, 1[.
When the input image is heavily corrupted by noise or composed of uneven textures, observing gradient distribution often fails to detect linear structures. To ease this problem, we also measure the intensity variance along the line segment. This is because intensities are likely to be homogeneous, if pixels are laid on the same line. We can write:

U v d (s i ) = 1 i 1 0 I (p i (t)) -E[I(si)] 2 dt, (2.4)
where E[I(si)] denotes the intensity mean of the line segment s i , and i is the line length.

Prior energy

In this section, we propose the prior energy to define spatial interactions on a local configuration. We want to obtain smoothly connected lines with a small curvature as a final solution. We compute the overlapping area Υ(s i , s j ) to reject congestion of lines and the coupling energy states c ij to evaluate attraction between line segments (see Fig. 2.1).

The prior energy U p (s i , s j ) is defined as where w p denotes a vector of weighting factors which control relative importance of each element in c ij . We assume that a line segment only correlates with the other ones within a certain distance. Thus, a neighborhood system consists of pairs of line segments, such that their center distance is smaller than half the sum of their lengths. In other words,

U p (s i , s j ) = Υ(si, s j ) + w ⊺ p c ij , ∀ i ∼ j, (2.5) (a) (b) (c) (d) (e) (f) (g)
i ∼ j = (s i , s j ) ∈ Ψ 2 ∶ 0 < x i -x j 2 ≤ i + j 2 + , (2.6)
where denotes the marginal distance to be connected with each other.

In order to evaluate an overlapping area between line segments, we dilate the line segments with a three pixel-radius disk, and then count up the number of pixels falling in the same image site. Suppose that we have a set of points A(s i ) which is a dilated version of the line segment s i , and n(A(s i )) denotes the number of pixels in A(s i ).

As shown in Fig. 2.1 (e)-(g), we penalize a configuration {s i , s j }, when a portion of the overlapping area is greater than 10% of min{n(A(s i )), n(A(s j ))}. However, almost perpendicular line segments are excluded from this penalty. The criteria for rejection are then given as

Υ(si, s j ) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 0 if θ ⊥ ij < τ , 0 if n(A(s i )∩A(s j )) min{n(A(s i )),n(A(s j ))} < 0.1, ∞ otherwise, (2.7)
where θ ⊥ ij = π 2 -θ ij represents an angle difference between s i and the perpendicular line of s j , τ is the maximum angle difference for segments to be aligned.

The coupling energy states c ij of the lines are composed of the singularity, connec-tivity, curvature, and perpendicularity:

c ij = [1, ϕ(d ij , ), ϕ(θ ij , τ ), ϕ(θ ⊥ ij , τ )] ⊺ , ϕ(u, v) = min{0, ( u v) 2 -1}, (2.8)
where d ij denotes the minimum distance from endpoints of s i to a point on the line s j , and θ ij is the angle difference between line segments. The function ϕ(u, v) tests a firmness of the coupling state u by comparing with the given tolerance value v.

The weighting factors w p = [ω s p , ω c p , ω a p , ω r p ] ⊺ can be derived from their role in the prior energy. Specifically, ω s p penalizes birth of a single line segment in the final configuration; hence its value is affected by the average gradient magnitude and the noise level of the input. ω c p encourages adjacent segments within to become connected. ω a p promotes segments being aligned with a small curvature in the final configuration. ω r p supports perpendicularly approaching line segments. Although the selection of w p values is critical for the performances of the MPP model, it is hard to estimate the coefficients because of hidden dependencies among them.

Monte Carlo sampler with delayed rejection

We employ the Reversible jump Markov chain Monte Carlo (RJMCMC) sampler [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF] to obtain an optimal line configuration which maximizes the probability density function.

The RJMCMC sampler is an iterative method that locally perturbs a current configuration s with a transition kernel. The transition kernel consists of multiple sub-transition kernels, namely, birth-and-death (BD) and linear transform (LT). A new configuration s ′ is proposed according to the transition kernel, given by

ξ s, s ′ = m p m ξ m s, s ′ , (2.9)
where p m denotes a probability to choose m-th type of sub-transition kernel ξ m (s, s ′ ). For each sub-transition kernel, the detailed balance condition [START_REF] Green | Reversible jump Markov chain Monte Carlo computation and Bayesian model determination[END_REF] is required to ensure the reversibility of the Markov chain. Acceptance ratio α m (s, s ′ ) is compared with a stochastic value rand[0, 1] to take a new configuration into account. The RJMCMC sampler is coupled with the simulated annealing (SA) algorithm [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] to secure the convergence of the Markov chain via relaxation parameter T (temperature); the temperature gradually decreases as the iteration goes on. To compute an acceptance ratio of the transition kernel, we use a density f (s) 1 T instead of f (s). The acceptance ratio is Note that we refuse a birth of the line lying on singular points, which have zero gradient magnitudes. On the other hand, the death kernel removes a line segment which is randomly picked from the current configuration. Thus, a new configuration s ′ = s {s} is proposed by the death kernel. We compute the acceptance ratio of the birth kernel α B and the death kernel α D in the same way as proposed in [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF], given by

α m (s, s ′ ) = min 1, ξ m (s ′ , s) ξ m (s, s ′ ) f (s ′ ) 1 T f (s) 1 T . (2.10) (a) s = {s 1 , s 2 , s 3 } (b) s {s 2 } ∪ {s ′ } (c) Searching (d) s {s 2 } ∪ {s}
α B (s, s ′ ) = min ⎛ ⎝ 1, p D p B µ(F) #(s) + 1 f (s ′ ) 1 T f (s) 1 T ⎞ ⎠ , (2.11)
α D (s, s ′ ) = min ⎛ ⎝ 1, p B p D #(s) µ(F) f (s ′ ) 1 T f (s) 1 T ⎞ ⎠ . (2.
12)

The LT kernel chooses a line segment s randomly, and then modifies its model pa-

rameters: s = (x, ( , θ)) → s ′ = (x ± dx, ( ± d , θ ± dθ))
, where dx, d , and dθ denote changes of center position, length, and orientation, respectively. The LT kernel draws a new configuration s ′ = s {s} ∪ {s ′ }. The acceptance ratio of the LT kernel is defined by

α LT (s, s ′ ) = min 1, f (s ′ ) 1 T f (s) 1 T . (2.13)
Algorithm 1 RJMCMC sampler with delayed rejection The LT kernel can be extended by the delayed rejection scheme [START_REF] Green | Delayed rejection in reversible jump Metropolis-Hastings[END_REF]. The main idea of the delayed rejection scheme is to give a second chance to a rejected sample point.

The acceptance ratio of delayed rejection is defined by

α 2 LT (s, s ′ , s) = min 1, ξ LT (s, s ′ ) ξ LT (s, s ′ ) ξ 2 LT (s ′ , s, s) ξ 2 LT (s, s ′ , s) [1 -α LT (s, s ′ )] [1 -α LT (s, s ′ )] f (s) 1 T f (s) 1 T , ≃ min 1, f (s) 1 T -f (s ′ ) 1 T f (s) 1 T -f (s ′ ) 1 T . (2.14)
where s ′ = s {s} ∪ {s ′ }, s = s {s} ∪ {s}, and ξ 2 LT (s, s ′ , s) is the transition kernel for the delayed rejection. In order to reduce the burn-in time, we add heuristics to design the delayed rejection kernel. When we propose an alternative line segment s, we look for the closest endpoints from both ends of s ′ , which is rejected from the first trial. The line segment s is generated by interpolation of the retrieved points; we force the connectivity of the neighboring segments, so that a probability of being accepted increases in terms of prior energy. Fig. 2.2 summarizes the process of the delayed rejection kernel, and Algorithm 1 provides the pseudo-code of the RJMCMC sampler with delayed rejection.

Curvilinear structure extraction via integration of line hypotheses

While the MPP allows to design complex prior knowledge of the object distribution, its performance is very sensitive to the selection of modeling parameters and hyperpa-rameters. For clarity, we note that the modeling parameters are related to the physical characteristics of the line segments (e.g., range of length and orientation). The hyperparameters denote the weighting coefficients of energy terms (i.e., w m d , w v d , and w p ). The modeling parameters can be chosen empirically since the values are related to the image resolution (see Sec. 2.2.4); however, it is hard to estimate the hyperparameters via trial-and-error for different types of dataset. Our goal is to maximize the probability density without estimating hyperparameters.

Generation of K line hypotheses

Let w = [ω m d , ω v d , ω s p , ω c
p , ω a p , ω r p ] ⊺ be a hyperparameter vector which consists of the weighting coefficients of the proposed probability density. Suppose that we have K different hyperparameter vectors, w 1 , . . . , w K . For each hyperparameter vector, we substitute kth hyperparameter vector w k into the proposed probability density f (s; w k ). Then, we look for its optimal configuration ŝk via Monte Carlo sampler proposed in Sec. 2.2.1.3.

Combination of line hypotheses into a probability map

We now have a family of line hypotheses Ŝ = {ŝ 1 , . . . , ŝK } obtained from K different hyperparameter vectors. We jointly use the image data and the line hypotheses. More specifically, the final solution s * maximizes not only the probability density but also the consensus among line hypotheses. For each optimal configuration ŝk , we compute a probability map P k of being a line in the image site. Then, we integrate K probability maps into a mixture density P Ŝ :

P k (x) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1 if ∃s k i ∈ ŝk , x ∈ s k i , 1 2 if ∃s k i ∈ ŝk , x ∈ A(s k i ), 0 otherwise, P Ŝ (x) = 1 K K k=1 P k (x).
(2.15) line hypotheses in the sense that the pixels corresponding to line structures are more highlighted when compared to [6,[START_REF] Talbot | Efficient complete and incomplete path openings and closings[END_REF].

We assume that the most promising hyperparameter vector draws a configuration which is more akin to the mixture density. We compute the correlation-coefficient (CC)

between P Ŝ and P k 's to analyze coherence of line detection results. That is

k * = argmax k={1,...,K} CC(P Ŝ , P k ), (2.16
) where k * represents the index of the most reliable hyperparameter vector.

CC(P Ŝ , P k ) = ∑ x P Ŝ (x) -E[PŜ] (P k (x) -E[Pk]) ∑ x P Ŝ (x) -E[PŜ] 2 ∑ x (P k (x) -E[Pk]) 2 , (2.17

Curvilinear structure extraction from reduced sampling space

The line hypotheses span a configuration space S ⊂ Ψ which will be considered as a new sample space. Since the size of S is significantly reduced compared to the original sample space Ψ, the optimization process becomes more tractable in terms of convergence time and detection accuracy.

We redefine the data energy by adding a new energy term as follows:

U ′ d (s i ) = U d (s i ) + U h d (s i ), U h d (s i ) = 1 0 -log P Ŝ (s i (t)) dt, (2.18)
where U h d (s i ) quantifies the consensus among line hypotheses with respect to the line segment s i . We stimulate the modified probability density over the reduced sample space S with the most promising hyperparameter vector w k * :

s * = argmin s∈S #(s) i=1 U ′ d (s i ) + i∼j U p (s i , s j ; w k * ). (2.19)

Experiments

We test the proposed method on a wide range of datasets: road cracks, facial wrinkles, DNA filaments1 , and retinas. Test images of the defects on the road pavements and ground-truth are courtesy of Chambon et al. [START_REF] Chambon | Road crack extraction with adapted filtering and Markov model-based segmentation[END_REF]. The facial wrinkle images are collected on the Internet, and forehead areas are manually selected for the experiments. We use the DRIVE dataset [START_REF] Staal | Ridge based vessel segmentation in color images of the retina[END_REF] to test the proposed algorithm on retina images.

For all test sequences, we fix the modeling parameters as follows: min is set to 5 pixels and max = 20 pixels. The orientation θ is varying from -90 ○ to 90 ○ with increments of 2 ○ . The marginal distance of connected segments is fixed to 2 pixels, and the maximum angular difference of aligned segments τ is 30 ○ . For the SA, the initial temperature T 0 is set to 100, and it follows the logarithm cooling schedule T t = T 0 log(1+t), where t denotes the number of the current iteration. We start the sampling process with the empty configuration. However, careful choice of initial segments can speed up the convergence of the algorithm (see Fig. 2.4). The computational time depends on the image resolutions; it takes less than a minute for the experimental images having 300 × 400 pixels, approximately. We use a PC with a 2.9 GHz CPU (4 cores) and 8 GB RAM.

To compare the performances of the proposed method with the state-of-the-art techniques, we apply the path opening operator [START_REF] Talbot | Efficient complete and incomplete path openings and closings[END_REF] on the gradient magnitude images by controlling the length parameters. For the supervised feature learning algorithm [6], we train 15 images for each dataset. In our experiments, we use the original implementations of path opening operator2 and supervised feature learning algorithm 3 . Fig. 2.5 shows the precision-and-recall curves for four test images. To obtain the curve of the comparison methods [6,[START_REF] Talbot | Efficient complete and incomplete path openings and closings[END_REF], we tune thresholds on line detection results.

The baseline MPP is selected from the line hypotheses among which it shows the best performance. The performances of the supervised learning algorithm are controlled by the quality of the training set; hence, it shows low performances on wrinkle and dna datasets, which are composed of noisy images with various sizes. In particular, the ground truth set of the wrinkle dataset is based on subjective perception. While the morphology operator enhances linear structures on gradient magnitude images, it is required to specify the length of the linear structures according to the target applications.

Since the pixelwise comparison fails to incorporate the geometry similarity with the ground-truth, the proposed algorithm shows lower scores on the crack and retina

Recall

Precision

PathOpening [START_REF] Chen | Semantic image segmentation with deep convolutional nets and fully connected CRFs[END_REF], L=10 PathOpening [START_REF] Chen | Semantic image segmentation with deep convolutional nets and fully connected CRFs[END_REF], L=90 Supervised learning [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for robust semantic pixel-wise labelling[END_REF] Baseline MPP Proposed method Human expert (d) retina Figure 2.5: Precision-and-recall curves for pixelwise segmentation of curvilinear structures using path opening operator [START_REF] Talbot | Efficient complete and incomplete path openings and closings[END_REF] with different setups of length, supervised feature learning [6], baseline MPP, and the proposed method. datasets. More specifically, the proposed algorithm detects slightly shifted lines for the crack image. Fig. 2.12 compares the detection results of the proposed MPP model with the manually labeled image by human expert, morphology operator [START_REF] Talbot | Efficient complete and incomplete path openings and closings[END_REF], supervised feature learning algorithm [6], and baseline MPP. For a fair comparison, we set the threshold values of the competing algorithms [6,[START_REF] Talbot | Efficient complete and incomplete path openings and closings[END_REF] to obtain the closest recall scores to the proposed algorithm. Blue pixels denote perfectly matching regions as compared with the ground-truth. Green and red pixels show over-detected and under-detected results, respectively. The main strength of the proposed algorithm is that it ensures stable performances for all datasets without any parameter estimation procedure. The proposed algorithm extracts the most salient line structures in the input image. On the other hand, the proposed algorithm suffers from under-detection when the width of the line structure is varying, for example, see the result for the retina. Such drawback can be overcome if we introduce an additional parameter for width of the line segment in our MPP model. sample space, so that the final solution interprets underlying curvilinear structures more faithfully. We have shown line detection results on a wide rage of datasets, and compared the performances of the proposed method with morphological filtering [START_REF] Talbot | Efficient complete and incomplete path openings and closings[END_REF], supervised learning [6], and baseline MPP method. The whole optimization process is friendly designed to the parallel implementation; thus, the computational time can be further reduced by applying the parallel Monte Carlo sampler [START_REF] Verdié | Detecting parametric objects in large scenes by Monte Carlo sampling[END_REF].

Concluding remarks

The main limitation of the method concerns its accuracy and suitability for real-life applications. Even though we showed that the proposed unsupervised approach performs fairly well for different kinds of images, the obtained accuracy may be not sufficient to be applied in a real-world scenario. Therefore, the current trend consists in applying supervised learning for applications dealing with curvilinear structure extraction.

Inference of curvilinear structure

In the following, I present our model based on supervised learning. More specifically, we employ structured SVM to learn the ranking function that predicts the correspondence between the given line segments and the latent curvilinear structures. To reconstruct the curvilinear structure, we build a graph using the output rankings of the line segments, and then look for the longest geodesic paths within this graph.

Overview of the proposed method

We first define notations and provide an overview of the proposed method (see Figure 2.7).

Assume that an image I contains a curvilinear structure. We denote the latent curvilinear structure Ω ∶ I ↦ {0, 1} for any pixel x of the image I:

Ω(x) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1 if x is on the curvilinear structure, 0 otherwise. (2.20)
This function is also related to a ground truth map, which is manually labeled, for the machine learning framework and for the performance evaluation. We compute a curvilinear feature map φ ∶ I ↦ R that represents oriented gradient information (see Section 2.3.2.1). Since information embedded in a single pixel is limited to infer the latent spatial patterns, we exploit image patches to compute input feature vectors. Let P x be a patch of the feature map values within

√ M × √ M size of square window centered at x, i.e., P x = {φ(x ′ ) x -x ′ ∞ ≤ √ M
2 }. Using a rotation matrix R θ defined on Euclidean image space, we can rotate a patch with respect to the given orientation θ such as Figure 2.7: Overview of the proposed inference method.

P x,θ = {φ(x ′ ) R ⊺ θ (x -x ′ ) ∞ ≤ √ M
2 }. We use graph theory for shape simplification of the complex curvilinear structure.

For the pixels showing higher rankings, we build an undirected and weighted graph G = (V, E). Then, we look for the longest geodesic path which corresponds to the coarset curvilinear structure in the image. We iteratively reconstruct the curvilinear structure by collecting paths that connect the remotest vertices in the graph. Consequently, the proposed algorithm can represent the different levels of detail in the latent curvilinear structure using the minimum number of pixels.

Curvilinear Feature

This section is devoted to compute the curvilinear feature descriptor that is used for the inputs of the learning system. We perceive the latent curvilinear structure based on inconsistency of background textures and its geometric characteristics. In other words, a sequence of pixels corresponding to the curvilinear structure has different intensity values compared to its surroundings, and shows thin and elongated shape. Thus, we compute multi-direction and multi-scale image gradients to detect locally oriented image features.

Curvilinear feature extraction

We obtain oriented gradient maps using a set of steerable filters [START_REF] Freeman | The design and use of steerable filters[END_REF][START_REF] Jacob | Design of steerable filters for feature detection using Canny-like criteria[END_REF][START_REF] Perona | Deformable kernels for early vision[END_REF]. Before applying the convolution operations, we normalize the training and test images to remove the effects of various illumination factors:

Ĩ = 1 1 + e - I-E[I] max(I)-min(I) , (2.21)
where E[I] is the sample mean of the image.

The steerable filters are created by the second order derivative of isotropic 2D Gaussian kernels. Let f θ,σ be a steerable filter that accentuates image gradient magnitude for direction θ at scale σ. To take into account varying orientations and widths of the curvilinear structure, a feature map φ is combined by multiple filtering responses:

φ = 1 Θ θ max σ {f θ,σ * Ĩ}, (2.22)
where Θ denote the total number of orientations. We first find the maximum filtering responses of the scale spaces, and then average for all directional responses. In this work we consider 8 orientations and 3 scales.

Local orientation estimation

In the previous subsection, we evaluate the presence of curvilinear structure from the amount of gradient magnitudes in image patches. Complex shaped curvilinear structure also consists of various local orientations.

Assume that B ∈ {0, 1} M be a basis spatial pattern of the curvilinear structure for the baseline orientation ( θ = 0 ○ ). We manually design B as a simple line shape template which highlights the middle of image patch (see Figure 2.8). If an image patch P x contains a part of curvilinear structure, pixels on the curvilinear structure are sequentially accentuated towards a particular direction θ. Recall that the geometric properties of the curvilinear structure, it is rotatable and symmetric. We can rotate image patch to be aligned with the basis spatial pattern. More specfically, we estimate the local orientation shown in image patches using statistical difference of the areas determined by B. For a rotated image patch P x,θ , we compute normalized histograms of the pixels which are labeled as curvilinear structure p x,θ = hist P x,θ B and its counterpart q x,θ = hist P x,θ ¬B , where hist (⋅) denotes histgram of the given distribution with 32 bins. Similar to [2], we employ χ 2 test to compute statistical difference of two distributions:

θ = argmax θ χ 2 p x,θ , q x,θ , (2.23)
χ 2 (p, q) = 1 2 i (p(i) -g(i)) 2 p(i) + g(i) . ( 2

Learning

We regulate the orientation and the shape of the input patches with respect to the basis spatial pattern B. In this section, we aim to learn a function that predicts structured output rankings of the input image patches. Thus, the output rankings indirectly infer the spatial patterns of the curvilinear structure for image patch.

Let vec (⋅) be an operator to convert an image patch into a column vector. For a pixel

x i , the input feature vector is defined as z i = vec P x i ,θ i B ∈ R N and y i ∈ R + denote the corresponding rank value. For a feature vector, we exploit subsampled pixels on the linear template to avoid data imbalance. Thus, the dimension of feature vector N is small when compared to the total number of pixels within the patch, i.e., N ≤ M .

For the setup of machine learning, a training dataset D = {(z i , y i )} K i=1 consists of the K input-and-output pairs. Let {z 1 , . . . , z K } ∈ Z be an unordered list of the input feature vectors and {y 1 , . . . , y K } ∈ Y be the corresponding classes. The input feature vector z is built from the curvilinear feature map φ. Our goal is to learn a ranking function h(z) which assigns a global ordering (ranking) of feature vectors:

h(z i ) > h(z j ) ⇔ y i > y j .
The inner product of the model parameter and a feature vector w ⊺ z is used to predict ranking score of the given image patch.

In our work Structured SVM framework [START_REF] Joachims | Training linear SVMs in linear time[END_REF][START_REF] Mittal | Taxonomic multi-class prediction and person layout using efficient structured ranking[END_REF] is employed to exploit the structure and dependencies within the output space Y. To encode the structure of the output space, a loss value ∆ of each input feature vector is defined:

h(z i ) > h(z j ) ⇔ ∆ i < ∆ j .
The loss value evaluates the quality of the learning system. Intuitively, the input image patches containing curvilinear structure are comparable to the basis spatial pattern.

Hence, as an image patch is similar to the basis spatial pattern, we should put it on a higher ranking. Specifically, we compute the loss value as the overlapping ratio between image patch from the groundtruth map Ω and the basis B:

∆ i = 1 - Ω x i ,θ i ∩ B Ω x i ,θ i ∪ B , (2.25)
where Ω x i ,θ i denotes a rotated image patch referring to groundtruth map Ω. The objective function of Structured SVM with a single slack variable ξ is given by:

min w,ξ≥0 1
2 w ⊺ w + Cξ s.t. 1 N ∑ (i,j)∈N c ij w ⊺ (z i -z j ) ≥ 1 N ∑ (i,j)∈N c ij -ξ ∆ j -∆ i , ∀(i, j) ∈ N , ∀c ij ∈ {0, 1} (2.26)
where c ij is the indicator variable to reduce the number constraint in linear complexity O(K):

c ij = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1 if ∆ i < ∆ j and w ⊺ z i -w ⊺ z j < 1, 0 otherwise. (2.27)
Cutting plain algorithm [START_REF] Tsochantaridis | Large margin methods for structured and interdependent output variables[END_REF] is employed to solve (2.26). For the details about the Structured SVM optimization, we refer the reader to [START_REF] Joachims | Training linear SVMs in linear time[END_REF][START_REF] Mittal | Taxonomic multi-class prediction and person layout using efficient structured ranking[END_REF][START_REF] Tsochantaridis | Large margin methods for structured and interdependent output variables[END_REF]. In the following section, we plot the initial segmentation map from the structured rankings scores. We also exploit these ranking scores for the shape simplification based on graph traversal algorithm.

Curvilinear Structure Reconstruction

To detect and represent the latent curvilinear structure, various local measures have been proposed. The score of the local measure is regarded as a likelihood probability whether a pixel is on the latent curvilinear structure. Most of the previous works represent the curvilinear structure as a binary map based on the linearity scores of the pixels. This approach easily misinterprets the topological features. In this section, we propose a novel structured score map based on the output ranking scores and the basis spatial pattern.

We also develop a graph-based model which is able to organize the topological features of the structure in different levels of detail. Figure 2.9 shows step-by-step processing results

to reconstruct the latent curvilinear structure. In other words, we have the information how each pixel is relatively important to represent the underlying curvilinear structure.

Structured score map

We have obtained a score function h(z; w) that evaluates the compatibility between the input features and the underlying curvilinear structure based on the structured output rankings. From the training dataset, we pre-compute the proportion ρ of pixels being part of the curvilinear structure. Note that the value ρ maximizes F 1 test scores for the groundtruth maps in the training dataset. During the test phase, we retain the output ranking according to the output ranking scores from the highest to the ρ I -th rank. As shown in Figure 2.9 (c), a binary map can be disconnected.

To avoid unwilling breakpoints, we composite the structured score map Π ∶ I ↦ R using the binary map, induced by ranking scores w ⊺ z, and the basis spatial pattern B.

Recall that we initially estimate the ranking score using the shape similarity between oriented image patch of ground truth and the basis spatial pattern in (2.25). Thus, it can be seen as the inverse mapping from output ranking score to the latent curvilinear structure. Let b = vec (B) be the column vector version of the binary mask B and Q x i ,θ i be the patch of structured score map Π centered at x i with θ i . We minimize the following cost function with respect to π i = vec Q x i ,θ i :

J = min x i ∈I ′ b -(w ⊺ z i )π i 2 2 .
(2.28) This is a least square problem, so that the solution is found at points which satisfy ∂J ∂π = 0. To reconstruct the structured segmentation map, we synthesize obtained patches Q x i ,θ i = vec -1 ( πi ) at the subsampled grid points on image x i ∈ I ′ , where πi = (w ⊺ z i ) -1 b.

We refer the readers to [START_REF] Kwatra | Texture optimization for examplebased synthesis[END_REF] for the implementation details of the related texture synthesis Algorithm 2 Progressive curvilinear path reconstruction Compute the longest geodesic path T in G ′ using [START_REF] Corneil | Diameter determination on restricted graph families[END_REF] 10:

1: Inputs: 2: G ′ = (V ′ , E ′ ) ∼
longest_path_length ← T w(u, v) ← 0, ∀{u, v} ∈ T , Update all edge weights on the path T as 0 for the given subgraph G ′ 16: end while technique.

Progressive curvilinear path reconstruction

We consider the graph G = (V, E) where V is the set of pixels. Two pixels u and v are connected by the edge if and only if uv ≤ √ 2. Moreover, we assign a weight for each edge {u, v} ∈ E as w(u, v) = uv Π(u)+Π(v) 2

. A path P in the graph G is a sequence of distinct vertices such that consecutive vertices are adjacent. The length of P is the sum of the weights of its edges, and the distance dist (u, v) between two vertices u and v is the minimum length of a path from u to v. The eccentricity ecc (v) denotes the maximum distance from the vertex v to a vertex u ∈ V , i.e., ecc (v) = max u∈V dist (u, v). The diameter diam (G) of G equals max v∈V ecc (v), i.e., it is the maximum distance between two vertices in G.

To simplify the latent curvilinear structure, we look for long geodesic paths in the subgraph G ′ of G induced by the pixels with structured score map Π. More precisely, our algorithm computes a diameter of G ′ , i.e., a shortest path T with length diam (G ′ ). This path T is added in the simplified curvilinear structure, then the path T is contracted into a single (virtual) vertex. For the implementation, the weight of all edges of T become 0.

We repeat this process till the diameter of the subgraph is larger than pre-defined path length ˆ . The entire procedure is summarized in Algorithm 2.

The time-consuming part of the proposed algorithm is the computation of a diameter of G ′ . Rather than computing all pair distances (which requires a linear number of applications of Dijkstra's algorithm), we use an efficient heuristic algorithm called 2sweep algorithm [START_REF] Corneil | Diameter determination on restricted graph families[END_REF]. The 2-sweep algorithm randomly picks a vertex t in G ′ , then performs Dijkstra's algorithm from t to find a vertex u at the maximum distance from t, i.e., dist (u, t) = ecc (t). Then, it computes (using Dijkstra's algorithm) a path from u to a vertex v at the maximum distance from u. The length of the second path (from u to v) is a good estimation of diam (G ′ ). Note that this algorithm is able to compute the exact diameter if the considered graph has a tree structure [16]. Topologically speaking, most of latent curvilinear structures are very close to trees [START_REF] González | Delineating trees in noisy 2D images and 3D image-stacks[END_REF][START_REF] Türetken | Automated reconstruction of dendritic and axonal tress by global optimization with geometric priors[END_REF]. Therefore, the 2-sweep algorithm is well adapted to reconstruct the tree-like curvilinear structures. For a better understanding, we schematically explain the intermediate steps of the proposed curvilinear structure simplification algorithm in Figure 2.10.

(a) (b) (c) (d) (e) (f)

Experimental results

In this section, we first discuss the parameters of the proposed algorithm and datasets.

We then compare the quantitative and qualitative results of the proposed algorithm and those of competing models proposed by [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF], [START_REF] Law | Three dimensional curvilinear structure detection using optimally oriented flux[END_REF], [START_REF] Becker | Supervised feature learning for curvilinear structure segmentation[END_REF], and [START_REF] Sironi | Multiscale centerline detection by learning a scale-space distance transform[END_REF].

Parameters and Datasets

The proposed algorithm requires few parameters to compute the curvilinear feature descriptor φ. We use 8-different orientations in this work: Θ = {0 All patch sizes are fixed to M = 33 2 . C which controls the relative importance of slack variables is set to 0.1 for all datasets. We choose the set of parameters for the SSVM training via 3-fold cross validation [START_REF] Kohavi | A study of cross-validation and bootstrap for accuracy estimation and model selection[END_REF] which maximizes the average F 1 score [START_REF] Martin | Learning to detect natural image boundaries using local brightness, color, and texture cues[END_REF] of the training set.

We test our curvilinear structure model on the following public datasets:

• Aerial [START_REF] Sironi | Multiscale centerline detection by learning a scale-space distance transform[END_REF]: The dataset contains 14 remote sensing images of road networks. We select 7 images for the training and 7 images for the test, respectively. We use ˆ = 80, σ 2 = {4, 8, 12}, and τ = 9.

• Cracks [START_REF] Chambon | Road crack extraction with adapted filtering and Markov model-based segmentation[END_REF]: Images of the dataset correspond to road cracks on the asphalt surfaces. We use 6 images to train and test the algorithms on different 6 images. For this dataset, we set to ˆ = 30, σ 2 = {2, 4, 8}, and τ = 3, respectively.

• DRIVE [START_REF] Staal | Ridge based vessel segmentation in color images of the retina[END_REF]: The dataset consists of 40 retina scan images with manual segmentation by ophthalmologists to evaluate the blood vessel segmentation algorithms.

We use 20 images for the training and 20 images for the test, respectively. The path length ˆ is set to 40. The scale factors and thickness are set to σ 2 = {2, 4, 8} and τ = 5, respectively.

• RecA [START_REF] Jeong | Marked point process model for curvilinear structures extraction[END_REF]: We collect electron microscopic images of RecA proteins on DNA which contain filament structure. We use 4 training images and 4 test images. We use ˆ = 30, σ 2 = {4, 8, 12}, and τ = 5.

Evaluations

The proposed algorithm progressively reconstructs the curvilinear structure by adding a long path on the subgraph. Figure 2.11 shows the intermediate steps of the proposed curvilinear structure simplification algorithm for DRIVE dataset. Unlike the previous models, the proposed algorithm is able to show different levels of detail for the latent curvilinear structure. Such information to visualize shape complexity of the curvilinear structure cannot be retrieved by setting a threshold. In practice, a few number of iterations is required to converge the algorithm and each step to find a long path takes less than milliseconds for the computation. For the experiments, we used a PC with a 2.9 GHz CPU (4 cores) and 8 GB RAM. Moreover, we visually compare the performance of the proposed algorithm with the competing algorithms. Figure 2.12 shows the results for Aerial, Cracks, DRIVE, and RecA datasets, respectively. The proposed algorithm is the most suitable to show the topological information of the latent curvilinear structures.

For the quantitative evaluation, we provide F 1 scores of the proposed algorithm and the state-of-the-art models in Table . 2.1. The measure of true positive is sensitive for the misalignment; therefore, we consider surrounding pixels of the detection results as the true positive if a predicted point is falling into the ground truth with a small radius (equivalent to its thickness parameter τ ) similarly to [START_REF] Sironi | Multiscale centerline detection by learning a scale-space distance transform[END_REF]. We also provide the average proportion of pixels to represent the curvilinear structures. It shows that the proposed algorithm efficiently draws the curvilinear structures using smaller number of pixels than the other algorithms.

The proposed algorithm achieved the best F 1 scores for all except the DRIVE dataset.

It is because the proposed algorithm use the fixed thickness τ to describe linear structure in the binary mask B; however, the images consisting of DRIVE dataset exhibit varied thicknesses of blood vessels and many junction points. In other words, we obtain a model parameter regarding to the manually designed binary pattern B. To overcome this drawback, we plan to exploit generative binary patterns based on the training images in the future, which is possible in our framework.

It is worth mentioning the risk of over-fitting in curvilinear segmentation task based on machine learning. Manual segmentation (ground truth) contains many errors around boundaries and minutiae components. Also, the number of training data employed in this work is relatively small due to the difficulties of making the accurate annotations.

Specifically, since RecA dataset has only four training images, there is a high risk of over-fitting to the training set. It is remarkable that the proposed algorithm shows good performance for all datasets without over-fitting problem when compared to other learning based algorithms [START_REF] Becker | Supervised feature learning for curvilinear structure segmentation[END_REF][START_REF] Sironi | Multiscale centerline detection by learning a scale-space distance transform[END_REF]. The second approach is based on the ranking learning system and graph theory. We learned a SSVM-based ranking function corresponding to the plausibility of the latent curvilinear structure. Furthermore, we proposed a novel graphical model that infers the curvilinear structure according to the topological importance. Across the various types of datasets, our model showed good performances to reconstruct the latent curvilinear structure with a smaller number of pixels when compared to the state-of-the-art algorithms.

Concluding remarks

I
The main limitation of the proposed methods lies in the lack of genericity, which is a well-known limitation of the template-based approaches. The approaches in this chapter have been tuned for applications dealing with curvilinear structures. In the next chapter, I will present a graph-based model which allows to introduce different soft shape priors into a classification process.

Chapter 3

Hierarchical Model for Image Classification

In this chapter, I present our study on the use of soft shape priors to improve the quality of classification. While classification is usually formulated as the minimization of an objective function, or energy, it is difficult to optimize such an energy if it involves shape priors because of their non-local nature [START_REF] Lempitsky | Image segmentation by branch-and-mincut[END_REF]. The state-of-the art methods require either the design of an optimizer specific to the particular shape prior [START_REF] Gorelick | Convexity shape prior for segmentation[END_REF][START_REF] Kemal | Uninitialized, globally optimal, graph-based rectilinear shape segmentation-the opposing metrics method[END_REF], or a complex way of incorporating it, when possible, into energies minimizable by standard techniques [START_REF] Das | Semiautomatic segmentation with compact shape prior[END_REF][START_REF] Veksler | Star shape prior for graph-cut image segmentation[END_REF]. Moreover, in the context of remote sensing we require to treat classification as a multi-object multi-class problem, since we typically find multiple instances of the different classes in a single image.

We have proposed a method for shape-aware classification, which takes into account the typical shape of object classes in terms of shape descriptors, such as rectangularity, compactness and solidity. For this, we use a hierarchical data structure, the binary partition tree (BPT) [START_REF] Salembier | Binary partition tree as an efficient representation for image processing, segmentation, and information retrieval[END_REF], which can be easily augmented to include shape information.

However, their traditional greedy construction approach does not yield a good utilization of the shape constraints. We thus proposed a method to iteratively optimize the structure of BPTs to produce better partitions with shape constraints. This enables us to perform multi-object multi-class classification with shape prior, and to enhance BPTs so that they better represent the underlying scenes.

Related work

Shape and optimization

In the literature regarding shape features for segmentation, we can distinguish contributions that focus on explicit shape models in a template matching manner and others geared at incorporating discriminative features, or soft shape priors (e.g., convexity, compactness).

In the context of template matching, a number of approaches have proposed to use the active contours framework [START_REF] Cremers | Kernel density estimation and intrinsic alignment for shape priors in level set segmentation[END_REF][START_REF] Michael E Leventon | Statistical shape influence in geodesic active contours[END_REF]. The evolution of the curves is usually slow and prone to get trapped in poor local minima. A second family of approaches finds global optimal solutions, but on high-dimensional specially constructed graphs. For example, Schoenemann and Cremers [START_REF] Schoenemann | Globally optimal image segmentation with an elastic shape prior[END_REF] proposed to look for minimal ratio cycles in the product graph of the input image and the shape template. These methods find a globally optimal segmentation in polynomial time on the high-dimensional graph. However, besides their complexity, these contributions are in general not suitable for fitting multiple templates in the same image [START_REF] Karantzalos | Recognition-driven two-dimensional competing priors toward automatic and accurate building detection[END_REF]. Regarding graph-based optimization, a template fitting method to fit was proposed by Fredman and Zhang [START_REF] Freedman | Interactive graph cut based segmentation with shape priors[END_REF], which must be run repeatedly to account for non-rigid deformations. Other iterative models have been also proposed by coupling graphical models with shape cues [START_REF] Huang | A graphical model framework for coupling MRFs and deformable models[END_REF][START_REF] Kumar | Obj cut[END_REF]. These approaches are computationally intensive, since every loop contains expensive operations.

To include discriminative features (e.g., compactness, ellipticity), Slaubaugh and Unal [START_REF] Slabaugh | Graph cuts segmentation using an elliptical shape prior[END_REF] proposed to repeatedly fit ellipses on minimal cuts for an elliptical prior. More recent contributions have managed to express certain shape priors (star-shape [START_REF] Veksler | Star shape prior for graph-cut image segmentation[END_REF], compactness [START_REF] Das | Semiautomatic segmentation with compact shape prior[END_REF][START_REF] Funka-Lea | Automatic heart isolation for ct coronary visualization using graph-cuts[END_REF]) in energies minimizable by traditional s-t cuts. These approaches rely on the ability to express the shape prior in the pairwise interaction term of a Markov random field which might be, when feasible, algorithmically complex (e.g., analyzing intersections with a rotating discrete line around a used-defined point [START_REF] Veksler | Star shape prior for graph-cut image segmentation[END_REF]). The trust regions framework [START_REF] Gorelick | Fast trust region for segmentation[END_REF] can minimize high-order functionals, and has been adapted for certain shape priors (volume, shape moments [START_REF] Gorelick | Fast trust region for segmentation[END_REF]). It requires to provide a linear approximation of the energy around the current solution. Gorelick et al. [START_REF] Gorelick | Convexity shape prior for segmentation[END_REF] expressed the convexity prior as the count of '1-0-1' configurations along a discrete set of lines. They proposed a linear approximation and a dynamic programming algorithm for its efficient computation. It is not clear however how to directly adapt the framework to further priors (e.g., rectangularity index) or to combine different priors under the same scheme.

Moreover, most of these techniques do not contemplate the occurrence of multiple object instances. For example, multiple convex segments are indeed penalized in [START_REF] Gorelick | Convexity shape prior for segmentation[END_REF] when counting the '1-0-1' sequences in the scene.

Multi-object segmentation

To cope with the simultaneous detection of multiple objects, the marked point process (MPP) framework has been used to fit an unknown number of parametric shape models.

Rectangles [START_REF] Ortner | Building outline extraction from digital elevation models using marked point processes[END_REF] and ellipses [START_REF] Descamps | Automatic flamingo detection using a multiple birth and death process[END_REF] are some of the geometric elements that have been considered, and I have described our MPP model for line elements in the previous chapter.

The optimal solution is sought by stochastic optimization. Karantzalos and Paragios [START_REF] Karantzalos | Recognition-driven two-dimensional competing priors toward automatic and accurate building detection[END_REF] have extended the active contours framework to fit multiple templates in a single image.

To our best knowledge, most recent contributions cited above (e.g., [START_REF] Freedman | Interactive graph cut based segmentation with shape priors[END_REF][START_REF] Gorelick | Convexity shape prior for segmentation[END_REF][START_REF] Veksler | Star shape prior for graph-cut image segmentation[END_REF]) require to isolate every object (with prior knowledge on their location) and segment it individually.

Hierarchical trees

The notion of hierarchy has been particularly exploited in several image analysis applications, such as the detection of objects by its parts [START_REF] Felzenszwalb | Object detection with discriminatively trained part-based models[END_REF] and depth ordering [START_REF] Palou | Occlusion-based depth ordering on monocular images with binary partition tree[END_REF].

A number of data structures have been designed to represent images as a hierarchy of regions, such as min/max-trees [START_REF] Salembier | Antiextensive connected operators for image and sequence processing[END_REF], α-trees [START_REF] Lefevre | Hyperspectral image classification from multiscale description with constrained connectivity and metric learning[END_REF] and binary partition trees [START_REF] Salembier | Binary partition tree as an efficient representation for image processing, segmentation, and information retrieval[END_REF]. In our work, we have focused on binary partition trees (BPTs), which have been particularly useful in various domains including remote sensing. Multi-label classifications can be extracted efficiently from a BPT by performing a cut on the tree that covers all pixels at arbitrary scales. Recent works have explored the use of shape descriptors during the cuts [START_REF] Valero | Object recognition in urban hyperspectral images using binary partition tree representation[END_REF][START_REF] Vilaplana | Binary partition trees for object detection[END_REF].

BPTs are constructed by successively merging regions with similar colors. Even though BPTs can constitute a good hierarchical approximation to the underlying structure of an image, the bottom-up approach propagates and amplifies the errors produced at the lower scales. As a result, it is very likely that nodes in the BPT will not represent complete significant objects [START_REF] Lu | Binary partition tree for semantic object extraction and image segmentation[END_REF][START_REF] Vilaplana | Binary partition trees for object detection[END_REF]. Previous works mitigated this problem by adding penalties, such as the growth of perimeters [START_REF] Vilaplana | Binary partition trees for object detection[END_REF], the elongation of the regions [START_REF] Kurtz | Extraction of complex patterns from multiresolution remote sensing images: A hierarchical top-down methodology[END_REF] or edge information [START_REF] Tilton | Incorporating edge information into best merge region-growing segmentation[END_REF]. The authors of [START_REF] Vilaplana | Binary partition trees for object detection[END_REF] fitted templates on partially detected objects. These approaches can only alleviate the effect.

As another consequence of the bottom-up approach, shape information cannot be used during construction: the ultimate shape of an object in a branch cannot be predicted by a portion of it. As a result, the criteria used at construction and processing of the trees are different (as in [START_REF] Valero | Object recognition in urban hyperspectral images using binary partition tree representation[END_REF][START_REF] Vilaplana | Binary partition trees for object detection[END_REF]), which limits the feasibility of the tool as it is to perform Figure 3.1: A binary partition tree (BPT) is a hierarchical subdivision of an image. An exhaustive partitioning can be extracted by "cutting" branches at different scales.

segmentation with shape priors.

Optimization of hierarchical trees

Optimization of hierarchical trees has been carried out in the area of computational phylogenetics, on the construction of phylogenetic trees. These trees represent the evolutionary relationships among species [START_REF] Li | Random tree optimization for the construction of the most parsimonious phylogenetic trees[END_REF]. Several common optimization algorithms (hill climbing, simulated annealing, genetic algorithms) have been applied on the tree structures [START_REF] Giribet | Efficient tree searches with available algorithms[END_REF]. The standard moves (known as branch-swapping or swappers) are nearest neighbor interchange, subtree pruning and regrafting (SPR) and tree bisection and reconnection (TBR). SPR is the pruning/paste of a subtree into another location, while TBR rearranges the subtree before pasting. The most common objective is to maximize the parsimony of the tree, i.e., to explain the observed data with the least evolutionary change. We have incorporated the idea of regrafting tree branches. Our optimization objective is however different, and our context requires to define moves that preserve the parent, child and spatial adjacency relations of BPTs.

Background on binary partition trees

Binary partition trees (BPTs) were pioneered by Salembier and Garrido [START_REF] Salembier | Binary partition tree as an efficient representation for image processing, segmentation, and information retrieval[END_REF] as a means to represent a set of meaningful image regions in a compact and structured manner. The root node corresponds to the entire image, the following level represents the subdivision of the entire image into two disjoint regions, and so on. It represents then a hierarchical abstraction of an image, which can be navigated to extract meaningful regions at different scales. The typical workflow involves an initial tree construction stage, followed by a second stage of information extraction from the tree. For example, once a tree is constructed, an exhaustive segmentation of the image can be obtained by performing a horizontal "cut" on the structure (see Figure 3.1). In this procedure, commonly referred to as pruning, branches can be selected at different scales, an inherent advantage of such hierarchical structure. Visual browsing [START_REF] Salembier | Binary partition tree as an efficient representation for image processing, segmentation, and information retrieval[END_REF], object localization [START_REF] Vilaplana | Binary partition trees for object detection[END_REF] and depth ordering [START_REF] Palou | Occlusion-based depth ordering on monocular images with binary partition tree[END_REF] are some of the application domains where BPTs have been used, beyond remote sensing.

The construction of a BPT is performed in a bottom-up fashion, by iteratively clustering pairs of similar regions. The starting point is an initial subdivision of the image represented by a region adjacency graph (RAG), where every node conveys a region and the edges link spatial neighbors (i.e., candidates for merging). The typical initial RAG is the pixel grid, though nothing prevents the approach from being used with other inputs too (e.g., a RAG of small regions containing similar pixels, known as superpixel segmentation). Every edge in the RAG is labeled with a dissimilarity value that compares the two associated regions.

BPTs are constructed by following a global mutual best fitting region merging approach [START_REF] Lassalle | A scalable tile-based framework for region-merging segmentation[END_REF]: at each iteration, the two most similar regions in the current subdivision are merged together (i.e., the least weighted edge out of all edges in the RAG). When a merge occurs, a new region is added to the BPT, connected to its two corresponding children, as illustrated in Figure 3.2. The process finishes when there are no more edges left in the RAG. A BPT thus records the history of merges that occurred during the execution of a region merging algorithm.

The overall process can be implemented efficiently by using an updatable priority queue structure on top of the RAG edges to keep track of the highest priority element.

Such a structure is first constructed in linear time and every subsequent update incurs in a logarithmic time cost. When two regions R 1 and R 2 are merged into a new region R 12 , one must update the RAG (and the associated priority queue). The edge connecting R 1 and R 2 must be removed, but let us also remark that all the edges adjacent to R 1 and R 2 must also be eliminated from the RAG, since none of both regions exists anymore.

We must then add the adjacency relations of the new region R 12 . The computation is straightforward: the neighbors of the new region R 12 are nothing but the union of the neighbors of the old R 1 and R 2 (with special care to remove any duplicates that may arise). The dissimilarity value associated to each of these edges must be computed and pushed to the priority queue. The complexity of the overall BPT construction process is O(n log(n)M ), n being the initial number of nodes and M the maximum number of neighbors of a merged region during the construction. Given that typically M ≪ n, the algorithm is quasilinear in practice. Recent advances have further accelerated the algorithm [1]. Among other improvements, the priority queue only stores the edge with the lowest dissimilarity value for each node. The other edges are brought back into consideration if the surrounding nodes are affected by a merge.

The final tree contains exactly 2n -1 nodes, which is a very space-efficient representation. Let us remark though that only a subset of all possible planar subdivisions is represented by the tree, which highlights the importance of research efforts to construct a good initial tree that conveys meaningful objects of the underlying image.

The key elements to define the behavior of a BPT are the region model, i.e. how regions are represented, and the dissimilarity function, i.e. the function to compare the region models, used to define the priority of the merges during tree construction. We review next the contributions related to these two elements.

Region model

The object-based nature of BPTs yields rich representations of the region that go beyond pixel spectra. Every BPT node can convey regional information, describing the region as a whole and not as a set of individual pixels. Examples of the regional data that can be associated to every node are the standard variation of the spectral signatures in the region, or shape features such as compactness.

To represent the spectrum of a region (and then compare it to the spectra of other regions) there are essentially two alternatives: parametric and non-parametric models.

A parametric model makes assumptions about the homogeneity or Gaussian distribution inside the regions. A typical parametric model is to represent the spectrum of a region as the mean spectrum of its pixels. Non-parametric models, on the contrary, consist of perband histograms of the pixel values, hence they represent the real observed distributions.

In hyperspectral imagery, non-parametric models have a better performance since they can describe the internal variability of a region [START_REF] Valero | Hyperspectral image representation and processing with binary partition trees[END_REF]. For example, a texture might correspond to several peaks in the histogram. When averaging spectra in regions with high variability, one might end up representing the region with a "false" spectrum that is not present in any of the individual pixels. In addition to spectral data, the model usually stores the area of the region, since it is commonly used in the dissimilarity function.

Dissimilarity function

To establish a priority for merging during BPT construction, it is required to provide a means to compare models of two regions. A dissimilarity function O(R 1 , R 2 ) typically used for this purpose comprises two factors as follows:

O(R 1 , R 2 ) = min( R 1 , R 2 ) β D(R 1 , R 2 ), (3.1)
where R i denotes the area of region R i . The first part of (3.1), min( R 1 , R 2 ) β , is the so-called area-weighting factor. This is an agglomerative force intended to cluster regions that are very small compared to the rest of the elements in the RAG. When no area-weighting is used (i.e., β = 0), the resulting BPT might isolate small noisy areas and connect them to the rest only near the root of the tree. With moderate values of β, small regions are merged at some point, forcing the tree to better look like a hierarchical subdivision. When β is too large, the trees may become too biased toward being balanced, hampering their representation capabilities. Even though this parameter is barely discussed in the literature, being mostly set to β = 0.5 or β = 1, we must point out that it is an arbitrary parameter that has to be selected. In our experience, no area-weighting leads to poor representations (e.g., the root containing two children: one noisy pixel and all the rest of the image), while low values of β solve this issue without biasing the trees too much. Alternatively, Calderero and Marques [20] proposed to keep track of the out-of-scale regions and force their merging at some point, while Valero et al. [START_REF] Valero | Hyperspectral image representation and processing with binary partition trees[END_REF] used a weighted sum of pixel values in a window to initialize the histograms, as a way of smoothing out outliers.

The second factor, D(R 1 , R 2 ), compares both regions based on their spectra. Kullback-Leiber divergence and Bhattacharyya distance are popular choices both in hyperspectral imagery and other types of images [20,[START_REF] Valero | Hyperspectral image representation and processing with binary partition trees[END_REF]. Spectra are seen as probability distributions and compared using standard information theory concepts. Every bin of one histogram is compared against the corresponding bin of the other histogram. However, using crossbin measures, which go beyond individual bins, has proven to be more robust [START_REF] Valero | Hyperspectral image representation and processing with binary partition trees[END_REF]. The average of Earth Mover's Distances [START_REF] Rubner | A metric for distributions with applications to image databases[END_REF] among histograms of all bands has also been used as a robust and efficient cross-bin dissimilarity function [START_REF] Palou | Occlusion-based depth ordering on monocular images with binary partition tree[END_REF]. Every distribution is seen as a pile of dirt, and the difference between two distributions is seen as the amount of work required to turn one pile into the other one.

Proposed method

In the following, we first pose our problem as the minimization of an energy, based on probability distributions of spectrum and shape features. We then describe how we build an initial BPT and process it to find the best partition represented by this tree. Finally, we describe an optimization algorithm to modify the initial BPT, in order to extract a classification that minimizes our energy.

Energy formulation

Let I = (I j ) 1⩽j⩽n be an input image containing n pixels. We assume we are given a set of possible object classes, as well as priors for each class. Multi-label segmentation consists in an exhaustive partitioning of the pixels into a non-overlapping set of regions R = (R i ), together with associated class labels L = (L i ), where labels L i belong to a set Ω of available labels. It can be stated as an optimization problem: minimize

E(R, L) = E C (I, R, L) + R i ∈R E S (R i , L i ), (3.2)
where E C expresses the color prior (quantifying how the segmentation fits the image spectrum), and E S , the shape prior.

Color prior

For each object class, we suppose we are given training examples, from which the color distribution can be estimated and used as a prior. Given a candidate segmentation (R, L), the color prior is defined as follows:

E C (I, R, L) = R i ∈R I j ∈R i -log P (L i I j ). (3.3)
One way of obtaining the posterior P (L i I j ) is to train classifiers based on the samples' colors, using support vector machines (SVM), and to extend them to output probability estimates as usual in classification problems [START_REF] Ting-Fan | Probability estimates for multiclass classification by pairwise coupling[END_REF].

Low-complexity shape features in BPTs

Similarly, the shape prior term is defined as follows:

E S (R i , L i ) = -R i log P (L i S i ), (3.4)
R i being the area of region R i , and P (L i S i ) being the probability of assigning the label L i to the region R i , given a vector S i of shape features of that region. Common regularization (such as boundary length [START_REF] Nieuwenhuis | A survey and comparison of discrete and continuous multi-label optimization approaches for the Potts model[END_REF]) can be incorporated as part of this term.

The weight on the area makes the per-pixel contribution of the color prior and the perregion contribution of the shape prior equally important.

We wish to enrich the nodes of BPTs by including shape information of the corresponding regions. Given that the optimization of the trees involves recomputing region descriptors, we must design a pool of features that can be computed efficiently from children nodes. In addition, the errors in estimating the shape descriptors in the finer levels should not be amplified in the upper ones.

Area can be efficiently computed by adding the areas of the children. Rectangularity and elongatedness shape descriptors have been used in the context of hierarchical methods [START_REF] Kurtz | Extraction of complex patterns from multiresolution remote sensing images: A hierarchical top-down methodology[END_REF][START_REF] Tarabalka | Improved hierarchical optimization-based classification of hyperspectral images using shape analysis[END_REF][START_REF] Valero | Object recognition in urban hyperspectral images using binary partition tree representation[END_REF], though no details on how to efficiently implement these features are provided. An oriented rectangle B of height h and width w is said to be the minimum area enclosing rectangle (MAR) of a region R if it is the rectangle of minimal area that entirely contains R. Rectangularity measures the resemblance to a rectangle, and is computed as R (h ⋅ w). Elongatedness measures the resemblance to a line, and is defined as w h.

We propose to store the convex hull of the region at every node. When two regions are merged, the convex hull of the new region can be computed by merging the convex hulls of its children. This can be done in linear time in the size of the input polygons by using rotating calipers [START_REF] Godfried | The rotating calipers: An efficient, multipurpose, computational tool[END_REF]. The MAR can be efficiently computed after it [START_REF] Godfried | The rotating calipers: An efficient, multipurpose, computational tool[END_REF].

Rectangularity, elongatedness and other descriptors like solidity [START_REF] Yang | A survey of shape feature extraction techniques[END_REF] (filled fraction of the convex hull) are then directly derived. In a balanced tree, which is enforced by our construction function, convex hulls incur in an O(n log(n)) increase of the storage required. Their computation does not increase the complexity of tree construction (proofs

in Appendix A.2).
In a discrete environment, the errors of computing the region areas converge as the areas grow, so do the convex hulls. As a consequence, the aforementioned features tend to be more precise in the upper levels of the trees.

Another useful shape descriptor is compactness [START_REF] Raul | State of the art of compactness and circularity measures[END_REF], related to the resemblance to a circle. It is typically defined as δR 2 (4π R ), where δR is the perimeter of R. However, this formulation imposes difficulties in a discrete environment [START_REF] Raul | State of the art of compactness and circularity measures[END_REF] due to the fact that the error in the estimation of δR does not converge. Li et al. [START_REF] Li | An efficient measure of compactness for two-dimensional shapes and its application in regionalization problems[END_REF] proved the robustness of computing compactness as R 2 (2πI g ), I g being the moment of inertia of the shape with respect to its centroid. Given that the centroid and moment of inertia of a region can be computed in constant time from the children, we propose this method to measure compactness in BPTs.

Let us assume that a probability density function p(s L) is available for every feature and class. These densities can be obtained by smoothing histograms of training samples [START_REF] Bernard W Silverman | Density estimation for statistics and data analysis[END_REF]. Let us call S = s 1 , ..., s m a vector of shape features. Assuming features' conditional independence, we have:

P (L S) ∝ m k=1 P (L s k ) = m k=1 p(s k L) ∑ L j ∈L p(s k L j ) .
(3.5)

Total energy

Combining Eq. (3.2), (3.3) and (3.4), the energy criterion to minimize is formulated as:

E(R, L) = - R R i ∈R j∈R i log P (L i I j ) + R i log P (L i S i ) (3.6)

Tree construction and processing

In this section we describe our proposed BPT construction approach, as well as the algorithm to extract the optimal classification from a BPT.

Supervised BPT construction

In Section 3.2, it was mentioned the region merging algorithm used to construct BPTs requires to define a model to represent the regions and a dissimilarity function to compare the regions. We choose a non-parametric model (i.e., spectral histograms) to represent the region, due to their well-known advantages (see Section 3.2). We have proposed, however, a modified dissimilarity function to construct the BPT [START_REF] Maggiori | Improved partition trees for multiclass segmentation of remote sensing images[END_REF].

In BPT-based classification methods (e.g., [START_REF] Valero | Hyperspectral image representation and processing with binary partition trees[END_REF]), the dissimilarity measure used to construct the trees is purely based on the comparison of spectral histograms. In fact, commonly used dissimilarity functions (Eq. 3.1) always penalize the merging of dissimilar regions. Let us recall that non-parametric models were introduced to represent and compare inhomogeneous regions, useful for textures and light gradients. However, internal class variability (e.g., an object composed by areas of different spectra) is not at all considered. In an unsupervised context, where there is no notion of object class, there is little hope to deal with this, since there is no reason to cluster dissimilar regions. However, when class probabilities are available we propose to include an additional force that clusters regions belonging to the same class, despite being spectrally dissimilar. The new function is as follows:

O(R 1 , R 2 ) = min( R 1 , R 2 ) β (1 -α)D(R 1 , R 2 ) -α log P (L R 1 =L R 2 ) . (3.7)
As in the original dissimilarity function (3.1), there is an area-weighting factor and an

unsupervised term D(R 1 , R 2 )
, which is computed by comparing spectral histograms of regions without any preliminary training. We here use the Earth's Mover Distance, due to its robustness to changes in illumination and its efficient computation [START_REF] Rubner | A metric for distributions with applications to image databases[END_REF][START_REF] Palou | Occlusion-based depth ordering on monocular images with binary partition tree[END_REF].

Equation 3.7 adds a supervised term

P (L R 1 =L R 2 R 1 , R 2 )
, the probability of assigning the same label to both regions. This way, while the unsupervised term penalizes spectral dissimilarity, the supervised term will encourage merging regions that are likely to belong to the same class. The trade-off between both terms is controlled by parameter α.

The term P (L R 1 =L R 2 R 1 , R 2
) is computed by marginalizing over the classes as follows:

P (L R 1 = L R 2 R 1 , R 2 ) = K j=1 P (L j R 1 )P (L j R 2 ), (3.8)
where K denotes the number of classes and P (L j R k ), with k ∈ {1, 2}, represents the probability of assigning a certain label L j to segment R k . We must now define a way to compute P (L j R k ) based on the posteriors of the individual pixels contained in the region. One way to do this is to compute the probability of assigning the label to all pixels, conditioned by the fact that all labels are known to be equal inside the region:

P (L j R k ) = x i ∈R k P (L j x i ) ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ωm∈Ω x i ∈R k P (ω m x i ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ . (3.9)
Alternatively, one can estimate P (L j R k ) by averaging the individual pixel probabilities:

P (L j R k ) = 1 R k x i ∈R k P (L j x i ). (3.10)
While the first expression is closer to a strict Bayesian interpretation, we found the second one to be a simple yet useful approximation.

By introducing (3.7) we expect to better cluster semantically significant objects, according to the classifier's output.

Best segmentation represented by a given tree

The common processing on BPTs consists in selecting the highest or lowest branches satisfying a given condition [START_REF] Lu | Binary partition tree for semantic object extraction and image segmentation[END_REF][START_REF] Valero | Hyperspectral image representation and processing with binary partition trees[END_REF]. However, some contributions have formulated the problem in terms of energy minimization [START_REF] Salembier | Low-level processing of PolSAR images with binary partition trees[END_REF][START_REF] Salembier | Binary partition tree as an efficient representation for image processing, segmentation, and information retrieval[END_REF]. In particular, Salembier et al. [START_REF] Salembier | Low-level processing of PolSAR images with binary partition trees[END_REF] have interpreted segmentation as a horizontal s-t cut on the tree (see Fig. 3.1), i.e., with a source at every leaf and a sink at the root. Let us denote τ a tree and C(τ ) the energy of the cut on τ with minimal (3.6) among all possible cuts. Our task is to find such a minimal cut.

Considering that the branches in the tree are independent, the globally optimal cut can be found by a dynamic programming algorithm. Let us denote by

E(R) = min L∈Ω E({R}, {L})
the lowest possible energy of a region R (by assigning the label that incurs the lowest cost). The tree is traversed in a bottom-up manner. Whenever a region R is visited, the following property is evaluated:

E(R) ⩽ C(R lef t ) + C(R right ), (3.11)
where R lef t and R right are the children of R. If the property does not stand, we set

C(R) = C(R lef t ) + C(R right )
and keep the best cuts of both children. Otherwise, we set

C(R) = E(R)
and replace the cuts by R with label L. This process is executed recursively until reaching the root of the tree. The overall algorithm is linear in the image size, since only one BPT traversal is required, and it guarantees the optimal cut in the space of solutions represented by the BPT.

Optimizing the trees

Even though the globally optimal cut on a BPT can be found efficiently, the organization of the nodes in the tree structure restricts the possible cuts that can be done on them.

In Fig. much better results, unless the classifier properly labels all regions. If we wish to include shape priors to enhance the results, let us point out that at the moment a and b were merged, it was impossible to know that b would eventually form a more significant object under a different sequence of merges, with the typical shape of a building. In other words, following this greedy construction approach we cannot predict the final shape of the objects in the output classification map. The resulting tree does not allow to perform any cut that would include the whole building into the same object, not even using strong shape priors, given that it is split through different branches. This is why we now propose to optimize the tree itself.

Our method consists in constructing an initial BPT with the usual region merging greedy algorithm, and then optimize it to extract a classification map that minimizes Eq. 3.6, thus incorporating the shape prior. To optimize BPTs we follow a local search approach, in which a solution is iteratively modified by performing local transformations on the trees, named moves.

Moves and associated updates

We propose a prune-and-paste move that prunes a branch of the tree and inserts it into another part of the tree. The pruned node must be pasted in a spatially adjacent location. In a balanced tree there are at most O(log(n)) ancestors of α and β, and, as stated before, for some of these ancestors the model of the regions must be updated. If we denote by K the complexity of updating a model (i.e., merging two children), the computation

of the new costs C is O(K log(n)). Usually K ≪ log(n), therefore the time is O(log(n))
in practice.

The adjacency relations in a BPT can be derived from the children by observing that the ancestors of two adjacent nodes are also adjacent, until the LCA. At the finest level the adjacency is known (e.g., a 4-connected grid).

Properties of the moves

We now explore some properties of the prune-and-paste move. In particular, we will

show that finding all possible energy decreasing moves does not require to exhaustively evaluate the energy gain of every possible move. Proofs are available in Appendix A.

Let us consider the situation of Fig. 3.4(b). In the second tree, a node was pasted at the position τ i . We wish to compare the effect of pasting higher instead (as in the third tree). Let us denote by τ i < τ j the is-a-descendant-of relation.

Proposition 1. Given a tree τ , suppose a node R m is pasted at τ i < τ 1 leading to a new tree ϕ. Let us consider an alternative move that pastes R m at τ j , with τ i < τ j < τ 1 , producing a tree ψ. In the cases where either

C(ϕ 1 ) -C(τ 1 ) ⩽ 0 or C(R m ) ⩾ C(ϕ 1 ) - C(τ 1 ), then C(ψ 1 ) ⩾ C(ϕ 1 ).
The proposition states that, if a move reduces the energy in the branch, a higher paste place will not do better. Under certain assumptions, if the move increases the energy, pasting higher will also increase it. Intuitively, pasting lower is more general.

Proposition 2. Let us consider a case where Prop. 1 hypotheses do not apply. There might then exist a higher paste place τ α so that C(ψ 1 ) < C(ϕ 1 ). Let us suppose that instead of pasting at τ α we paste at τ β , with τ α < τ β < τ 1 , leading to a tree ρ. Then C(ρ 1 )

would monotonously decrease as the paste place τ β is located higher.

When the hypotheses of Prop. 1 do not apply, there might exist a favorable paste place higher in the branch. However, we know that the higher it is, the most beneficial it can be. As a result, we can just consider the highest possible paste place (right below the LCA). However, any paste place between the location of the original cut and the LCA would lead to the same energy. These cases happen when it is preferable to cut the pruned node apart. The higher we paste R m , the less we condition the way the rest of the tree must be cut.

Following these properties, an exhaustive search of energy decreasing moves can be achieved as follows: for every possible pruning place we check the energy gain of the moves for only the lowest paste places. If the move is beneficial, we keep that move as a candidate. If the move is not beneficial we can discard it as well as all the paste places in the branch (but we do this only after performing one additional check with the paste place at some point between the cut and the LCA).

Optimization algorithm

We propose the following optimization scheme, which must be iterated:

1. Construct a heap of all moves according to the branch cost variation ∆C.

(a)

Either apply the best move, or (b) apply the best k moves (when still appropriate).

The first step involves exploring the whole search space, featuring an O(n log 2 (n))

complexity. Moves are simulated to measure the energy gain but are not applied. Propositions (1-2) can be used to reduce the execution time of this step, and the energy gain can be evaluated in parallel.

In the second step, energy decreasing moves are applied. As soon as a move is applied, the tree is restructured and the effect of some other moves may be affected. New energy decreasing moves may also arise. 2a) just applies the best move and reiterates. As a shortcut, we can just update the entries in the heap of the moves that may have been affected. This option is still costly because ∆C must be recomputed for any move that could have been possibly affected, even though in practice this might be the case for just a few of them. This approach guarantees to apply the best move each time. Considering the fact that there might be many unrelated energy decreasing moves in the tree, 2b)

proposes to apply a number k of best moves, but verifying for each move that ∆C did not increase as a result of the previous transformations done on the tree. This second approach will apply a number of independent moves first, ignoring the fact that some new energy decreasing moves might arise, which will be dealt with in the next iteration.

The loop stops when there are no more moves whose ∆C is negative.

The lowest scale both on the pruning and paste side can be set to the pixel level. A coarser scale on the pruning side can be used to adjust the precision of the moves. If a coarser scale is set on the paste side, we limit the minimum object size of the partitions, since we ignore moves that lower the cut below a certain level. This can be adjusted by observing the density functions of the area feature.

Experiments

We first validate our supervised BPT construction approach. We do this on hyperspectral images, since they contain a rich and discriminant spectrum, and evaluate whether our function (3.7) indeed clusters regions together based on the classifier's probabilities.

However, such classification power based solely on pixel's values is notoriously degraded in other types of images. We therefore introduce shape features in order to enhance classification on different types of inputs, including urban aerial remote sensing images.

Supervised BPT construction

We perform experiments on the Pavia Center hyperspectral dataset, acquired with the Reflective Optics System Imaging Spectrometer (ROSIS-03). The image has spatial dimensions 400 × 300 and contains 102 bands, covering a range from 0.43 to 0.86 µm and 1.3 m spatial resolution. A color composition of the image is shown in Figure 3.

5(a).

A reference image that labels entire objects was built, including four classes (see

E(R, L) = λ R - R i ∈R I j ∈R i log P (L i I j ), (3.12)
where λ controls the coarseness of the output. We experimentally set λ = 40 to optimize the accuracy. The resulting classification map with supervised tree construction (α = 0.5) is shown in Figure 3.5(d). We also executed a graph cut with α-expansion [19] on the probabilities derived from the SVM, which proved to be effective in the past for hyperspectral image classification [START_REF] Tarabalka | Graph-cut-based model for spectral-spatial classification of hyperspectral images[END_REF]. Its regularity parameter was also set empirically to optimize the accuracy. From these illustrations we can clearly appreciate that including class probabilities during BPT construction yields an improved tree that better clusters the objects together.

To validate this numerically we compute the overlap between every building (belonging either to tiles or bitumen classes) in the reference data and the most overlapping building region in the BPT output. The overlap is measured with Dice's coefficient [START_REF] Dice | Measure of the amount of ecological association between species[END_REF] defined

as: 2 R 1 ∩ R 2 ( R 1 + R 2 ).
The resulting overlap coefficients are averaged over all reference buildings to estimate how the BPT output matches the reference data from an object-based perspective. The numerical results, together with the overall accuracy, are summarized in Table 3.1, which also includes the values for SVM and graph cut. A first observation we can make is that BP T α=0.5 performs better than BP T α=0 , corroborating the visual result from Figs. 3.5(e-f). Secondly, while graph cut is known to improve the SVM classification, we can verify that while this is true from a pixelwise perspective (in terms of overall accuracy), it is not from an object-based perspective (in terms of building overlap). Finally, the use of BP T α=0.5 outperforms the other methods in terms of object overlap. This validates the idea of including class probabilities during tree construction.

Classification with shape priors

In a first series of experiments we constrain our method to binary segmentation with convexity shape prior, in order to compare it with a recent state-of-the-art technique designed for this purpose. In a second series of experiments we move to multi-object multi-class segmentation in remote sensing imagery and compare the behavior of our algorithm against the common approaches in the field.

An algorithm for a soft convexity shape prior in image segmentation was introduced by Gorelick et al. [START_REF] Gorelick | Convexity shape prior for segmentation[END_REF]. We will first show that we can perform similarly in a single convex object, while our technique is not specifically designed for this prior. A natural image extracted from [START_REF] Gorelick | Convexity shape prior for segmentation[END_REF] and the markers used are shown in Fig. 3.7(a). We set ω = 1, a contrast sensitive Potts model (weight 10 -4 ) and an 8 × 8 orientation stencil, the parameters required by [START_REF] Gorelick | Convexity shape prior for segmentation[END_REF]. We found this parameter setting to yield the most visually pleasant result. We ran our BPT optimization approach with a data term learned from markers' histograms and with density functions favoring solidity and compactness for the foreground class, which convey the notion of convexity. For BPT construction we set α = 0.5 throughout all the experiments, to make the terms in (3.7) equally relevant.

We sample moves that involve regions of at least ten pixels. In this case we apply all good moves in the queue at every iteration, which are in practice less than 10. We did not impose hard constraints on the marker locations. The resulting segmentations are wish to identify cell nuclei (green markers). If we apply the method by Gorelick et al. [START_REF] Gorelick | Convexity shape prior for segmentation[END_REF] to the whole image with all the markers, the result is inaccurate because the technique is not designed to segment more than one object at once (Fig. 3.8-b). In Fig. 3.8(c) we overlap the result of applying [START_REF] Gorelick | Convexity shape prior for segmentation[END_REF] to different fragments of the image that include each individual object (parameter ω in [START_REF] Gorelick | Convexity shape prior for segmentation[END_REF] is set to 0.01). The technique individually outlines each of the nuclei in the absence of the rest, but prior knowledge about their location is required.

We constructed a BPT on this image (α in Eq. 3.7 is set to 0.5 in all the experiments).

As expected, the best-cut segmentation without optimizing the tree (Fig. 3.8-d) is not competitive, because many tree nodes do not satisfactorily represent objects. We ran our BPT optimization approach with a data term learned from markers' histograms and with density functions favoring solidity and compactness for the foreground class, which convey the notion of convexity. We sample moves that involve regions of at least ten pixels. In this case we apply all energy decreasing moves in the queue at every iteration, which are in practice less than 10. We did not impose hard constraints on the marker locations. After optimizing the BPT with the method proposed in Sec. 3.3.3, each object is segmented individually and adjacent objects are delineated separately (Fig. 3.8-e).

Our approach also produces accurate boundaries despite the low foreground/background contrast. The method by Gorelick et al. tends to either oversmooth the boundary to enforce convexity, or produce a very non-convex object, depending on the parameter ω (see the amplified nucleus in Fig. 3.8-f).

In the context of multi-object multi-class segmentation, we tested our method on images of urban scenes extracted from Google Maps screenshots. Figs. [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for robust semantic pixel-wise labelling[END_REF].9(a) and 3.10(a)

show two color images acquired over New York City and the area of Brest, respectively.

For both images, the manual segmentation was based on visual inspection combined with cadastral records available online. The list of the considered object classes is given in Fig. 3.9(b) (no instance of the internal road class is present in the Brest image). In the particular case of buildings, cadastral information was used to delineate every object independently even when they are spatially adjacent.

To evaluate the performance of the proposed method, we use two criteria:

1) Overall accuracy A, defined as the proportion of correctly classified pixels.

2) Building's overlap D. For every building in the manually segmented image, we search for the most overlapping building region in the segmentation map in terms of Dice's coefficient [START_REF] Dice | Measure of the amount of ecological association between species[END_REF]. Criterion D is estimated by averaging the computed coefficients.

As in Section 3.4.1, an SVM with a Gaussian radial basis function kernel was used for the data term, tuned by tenfold cross-validation. The criterion (3.4) involved area, rectangularity and elongatedness shape descriptors. The distributions were trained by kernel density estimation on a set of sample objects from an adjacent image. In the area covered by these objects, 100 random pixels per class were selected to train the SVM.

We compared the performance of the proposed approach with the following methods:

1) SVM; 2) graph cut with α-expansion [19] (GC); 3) cut on the BPT, regularized by the number of regions without using shape priors (TC) [START_REF] Salembier | Low-level processing of PolSAR images with binary partition trees[END_REF]; 4) cut on the same BPT with our shape formulation (3.6), but without tree optimization (TSC). The SVM classification exhibits many issues, notably the assignment of some roof parts to the wrong class. GC and TC smooth the results though do not correct the main mistakes in the classification. In the initial cut with shape priors on the unoptimized tree (TSC), some regions are enhanced but some others are significantly deteriorated with respect to the previous methods. This is due to the faulty tree construction that does not represent the entire objects in unique nodes. Figs. 3.9(g-h) show that the optimization of the tree copes with these issues, not only enhancing the initial cut on the tree but also outperforming the other techniques.

The evolution of the energy (3.6), the accuracy A and the building's overlap D with respect to the number of iterations are depicted in Figs. [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for robust semantic pixel-wise labelling[END_REF].9(i-k). As expected, the energy curve becomes less smooth as k increases. For small enough values of k, the segmentation maps are almost identical. This validates the fact that many branch moves are independent and can be applied prior to reconstructing the heap.

The BPT for this image is constructed in 1.25 seconds on an 8-CPU 2.7 GHz processor.

The optimization time, summarized in Fig. 3.9(l), is considerably faster when Propos. 1-2 are used. We highlight in another portion of image over Brest (Fig. 3.11) that our method is able to separate an entire building blob into rectangles. In addition, the use of shape features permits to "switch" a small building to the correct class, even though its color was assigned to road by the classifier.

Concluding remarks

I have presented a hierarchical model for the multi-class multi-object classification of images with shape priors. The problem is formulated as an energy minimization task, using the learned probability distributions of spectrum and shape features. We construct a BPT of the image and search for the horizontal cut on the tree which minimizes the proposed energy function. One of the key innovations of the proposed method is that instead of using the standard BPTs, we optimize them by pruning and regrafting their branches at different scales in order to minimize the best classification that can be extracted from the tree. BPTs are a good departure point for optimization, since they allow to perform moves of variable region sizes. Conversely, in an MRF approach defined on the pixel grid, each time we perform a move we would have to search for a relevant region to perform a label switch (e.g., by using minimum spanning forests). A BPT can in fact be seen as a way of dynamically storing a set of candidate adjacent regions at different scales. Our theoretical study permits to reduce the space of branch moves.

We also proposed an improved BPT construction function that adds a force to cluster together regions that are similar in terms of the classifier's output.

The experiments show that the method effectively incorporates shape information into classification, even when multiple classes and objects are present in the image, thus gaining a competitive edge with respect to recent literature in the field. The main limitation of the method concerns its scalability, both in terms of potential for parallelisation and of the choice of features. We will address how to handle these limitations later in this manuscript.

In the next chapter, I will present another graph-based model, which aims at incorporating a shape prior for segmentation of time series of noisy images. This time, we focus on segmenting shapes which only grow or shrink in time, and our model is designed to suit this growth or shrinkage constraint, respectively.

Chapter 4

Video Segmentation with Shape Growth or Shrinkage Constraint

Automatic segmentation of objects in videos is a difficult endeavor in computer vision [START_REF] Lézoray | Image Processing and Analysis with Graphs: Theory and Practice[END_REF]. This task becomes even more challenging for image sequences with low signal-to-noise ratio or low contrast between intensities of spatially adjacent objects in the image scene. Such challenging data are recorded frequently, for instance, in satellite remote sensing or medical imaging.

In this chapter, we aim at segmenting objects with shapes which can only grow or

Related work

Segmentation of time series

Image segmentation methods applied independently to each frame [START_REF] Winnemöller | Real-time video abstraction[END_REF][START_REF] Chen | Real-time edge-aware image processing with the bilateral grid[END_REF] produce unstable results, while temporal coherence in video sequences yields a lot of information not available for a single image. There are two main categories of approaches for the spatiotemporal segmentation of image sequences. Causal, or feedforward, techniques consider only past data for segmenting each next frame [START_REF] Paris | Edge-preserving smoothing and mean-shift segmentation of video streams[END_REF][START_REF] Tarabalka | Shape-constrained segmentation approach for arctic multiyear sea ice floe analysis[END_REF]. Omniscient approaches take advantage of both past and future data by analyzing the video as a 2D+T = 3D spatiotemporal pixel volume [START_REF] Shi | Motion segmentation and tracking using normalized cuts[END_REF][START_REF] Dementhon | Spatio-temporal segmentation of video by hierarchical mean shift analysis[END_REF][START_REF] Grundmann | Efficient hierarchical graph based video segmentation[END_REF]. In this case, the segmentation of the entire image set supports each of the individual segmentations.

Graph-based methods gained popularity among omniscient approaches, in particular those using hierarchical model [START_REF] Grundmann | Efficient hierarchical graph based video segmentation[END_REF], normalized cuts [START_REF] Shi | Motion segmentation and tracking using normalized cuts[END_REF] or graph cuts [START_REF] Wolz | Measurement of hippocampal atrophy using 4d graph-cut segmentation: Application to ADNI[END_REF]. These techniques do not impose any shape prior knowledge. It was proven that introducing shape priors into image segmentation, i.e. favoring segmentations similar in some sense to a given shape, allows to drastically improve segmentation of objects in the presence of strong noise and occlusions [START_REF] Cremers | Diffusion snakes: introducing statistical shape knowledge into the Mumford-Shah functional[END_REF][START_REF] Rousson | Shape priors for level set representations[END_REF]. However, imposing shape priors increases significantly both algorithmic and computational complexity of segmentation algorithms [START_REF] Felzenszwalb | Representation and Detection of Shapes in Images[END_REF].

Schoenemann and Cremers proposed to compute minimal ratio cycles in a large product graph spanned by the image and the shape template for finding globally optimal segmentations, which are consistent both with edge information and with a shape prior [START_REF] Schoenemann | Globally optimal image segmentation with an elastic shape prior[END_REF][START_REF] Schoenemann | Globally optimal shape-based tracking in realtime[END_REF].

In order to speed up the algorithm, they implemented it in parallel on graphics hardware.

This algorithm can be applied for segmenting shapes in time series in a causal way, so that the template determined for the last frame is matched to the next. However, it does not guarantee globally optimal solution over the whole temporal sequence.

Globally optimal segmentation

Extensions of graph cuts to multi-class segmentation have been proposed but generally do no guarantee optimal solutions, except e.g. in the case of Ishikawa's construction [START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF],

which requires labels to be ordered and the interaction term to be a convex function of their differences. This graph construction makes intensive use of infinite links to constrain the min-cut solutions to satisfy desired properties required to interpret them as image segmentation solutions. This was the source of inspiration for our work.

A study related to shape constraints can be found in [START_REF] Delong | Globally optimal segmentation of multi-region objects[END_REF], where one image has to be segmented in several possibly-overlapping objects. Infinite links are used for imposing common boundaries, inclusion or exclusion conditions between objects in a same single image. A similar approach [START_REF] Li | Optimal surface segmentation in volumetric images -a graph-theoretic approach[END_REF] segments jointly two surfaces in a same volumic image, under the constraint that they should be separated by a given minimal margin. There is however no work related to shape growth or shrinkage in time series.

Wolz et al. [START_REF] Wolz | Measurement of hippocampal atrophy using 4d graph-cut segmentation: Application to ADNI[END_REF] applied graph cuts for simultaneous segmentation of serially acquired MRI volumes. They defined temporal edge weights as the intensity differences of voxels at the same spatial locations. The same smoothness constraint was applied both in space and time, and the segmentations at different timepoints were forced to be consistent in areas where a small intensity difference between the images exist. This type of temporal constraint is suboptimal in image series where intensity distributions of foreground and background vary significantly over time. To the best of our knowledge, our work is the first to use infinite links to enforce a temporal growth constraint, and we illustrate in Sec. 4.3 the advantage of the new method over previous approaches such as [START_REF] Wolz | Measurement of hippocampal atrophy using 4d graph-cut segmentation: Application to ADNI[END_REF].

Proposed method

In the following, I first recall a formulation of the graph cut, and then present the proposed method. Graph cut is an optimization tool coming from graph theory, based on the rewriting of image segmentation problems as (s,t)-min-cuts in graphs, on the , between all pixels of coordinates (x, y, t, s 1 ) and (x, y, t, s 2 ). equivalence of (s,t)-min-cut and max-flow problems, and on the existence of efficient algorithms to solve the latter ones [18,[START_REF] Dantzig | On the max-flow min-cut theorem of networks[END_REF][START_REF] Goldberg | A new approach to the maximum-flow problem[END_REF]. In practice in computer vision, it can be used to find the globally optimal binary segmentation of images where the segmentation criterion is related to a Markov Random Field with submodular interaction terms, i.e. a criterion E of the form:

E(L) = pixels i V i (L i ) + i∼j W i,j (L i , L j ), (4.1)
where L is the binary labelling function to be found (L i is the label of pixel i), individual potentials V i are any binary real-valued functions measuring the disagreement between a prior probabilistic model and the observed data, i ∼ j denotes a pair of neighboring pixels (any neighborhood system can be used), and W i,j are any real-valued interaction terms between neighboring pixels expressing spatial coherency of labels, satisfying

W i,j (0, 0) + W i,j (1, 1) ⩽ W i,j (0, 1) + W i,j (1, 0). (4.2)
A directed infinite link between two pixels expresses precisely the constraint that this pair of pixels cannot have the pair of labels (0,1), by assigning an infinite cost to such an interaction. The remaining possible pairs of labels are thus (0,0), (1,1) and (1,0), which means that either both pixels have the same label, or the order of labels is predefined

(1 for first pixel and 0 for the second one). In the case of image binary segmentation, if 0 stands for the background and 1 for the foreground object, then this means that the second pixel may belong to the foreground only if the first one already does.

Growth/shrinkage constraint

Given a sequence of images I(t) preliminarily aligned, shape growth can be easily expressed as the property that the foreground object cannot lose any pixel when time advances. Otherwise said, if a pixel belongs to the foreground object at time t 1 , then it belongs also to the foreground object for all times t 2 > t 1 . Equivalently, and simpler: a pair of pixels ((x, y, t), (x, y, t + 1)), sharing the same location and immediately successive in time, cannot have the pair of labels (1,0), with the same binary segmentation notations as above. This can be enforced by setting monodirectional infinite links from all pixels to their immediate predecessor in time.

Given T images I(t), with t ∈ [1, T ], and as many associated submodular segmentation criteria E t , we transform the problem of segmenting independently each image I(t) according to its criterion E t , into a joint segmentation of all images together, by enforcing the shape growth constraint with directed infinite links (see Fig. 4.1(a)). Thus, instead of applying graph cut T times independently to planar grids of the size of the images W × H, we apply graph cut once to a 3D grid W × H × T , consisting of the same nodes and edges, but with additional monodirectional infinite links in time. The criterion to be minimized is then E = ∑ t E t under the constraint of shape growth:

E(L) = pixels i V i (L i ) + i∼j W i,j (L i , L j ) + ∞ t δ L t i >L t+1 i . (4.3)
Since the problem is binary and submodular, the solution found by graph cut is globally optimal.

Manifestly, one can enforce shape shrinkage instead of shape growth, by reversing the direction of the infinite links. Another straightforward extension, needed in Sec. 4.3.3, consists in applying this approach to the case of sequences of 3D images. The directed infinite links are then set for all pairs of voxels of the form ((x, y, z, t), (x, y, z, t -1)) to enforce 3D shape growth.

In some applications, it may happen that growth (or shrinkage) is only very probable, but not with probability 1, i.e. growth should be considered as a probable hint but should not be enforced strictly at all locations at all times. In that case, one may replace directed infinite links by directed finite links: the weights of these links will encourage growth (more or less strongly depending on the weight), but sufficiently disagreeing potentials V i may make the shape locally shrink instead. Thus, shrinkage would be discouraged but not forbidden.

Inter-sequences inclusion constraint

It is also possible to segment jointly several image sequences I(s)(t) with the constraint that the foreground object in some sequences should be included in the foreground object of some other sequences. This can be done similarly by considering together the graphs associated to all sequences, and, for each inclusion constraint, by adding directed infinite links between pixels of the desired sequences s 1 and s 2 , sharing same location and time:

such links from (s 1 , x, y, t) to (s 2 , x, y, t) for all x, y, t will force the foreground object in sequence s 1 to contain the one of sequence s 2 (see Fig. Naturally, instead of imposing an inclusion constraint over the whole time span and the whole image space, it is possible to specify spatio-temporal domains of constraints, for instance to express that the inclusion property between two sequences has to be satisfied inside a pre-defined region and/or during a pre-defined time span [t 1 , t 2 ] only, by adding directed infinite links in these sets only.

Weighting frames by reliability

As presented earlier, enforcing shape growth explicitly is particularly important when facing noisy sequences of images, as a joint segmentation incorporates information from different frames. With our approach, a fully noisy frame in the middle of several goodquality frames will be automatically ignored in practice, because statistically the (bad) potentials V i at any pixel i in that frame will be neglectable with respect to all other (correct) potentials at the same pixel in other frames. Hence, the solution we bring will be strongly robust to any kind of noise, provided that the noise is not coherent in time during long time spans.

The quantity of noise that the method can handle depends on the number of neighboring frames with good-quality information at the same spatial location. One way to increase even more the trade-off between admissible level of noise and quantity of good information consists in estimating the reliability of the information given, in each image or even at each pixel, and in weighting accordingly the associated energies :

E = ∑ t w t E t .
For instance, if a strong level of noise is detected globally in one image I(t), one may multiply the corresponding energy E t by a small reliability factor w t < 1, in order to make it less influential than other frames. Such a reliability factor could be computed for instance from preliminary image processing steps (for example, correlation between an image I(t) and its neighboring ones {I(t ′ ) ; 0 < tt ′ ⩽ δt} ).

Complexity

The precise theoretical worst case complexity depends on the graph cut algorithm used, and is the same as for usual single-image graph cuts, but with T times more nodes However, the computational time observed in practice is known to be much faster on typical image segmentation problems. We applied the binary graph-cut algorithm of Boykov and Kolmogorov [18], and we report a linear observed complexity with the total number of pixels T ×W ×H (see for instance Fig. 4.5 for experiments of section 4.3.1). As a consequence, enforcing shape growth on a long sequence, or on complementary shorter bits of the same sequence, will take approximately the same time.

In the case of long sequences of big images, the memory space required may exceed the capacities of a computer. This is however not an issue, as there exist graph cut implementations for massive grids [START_REF] Delong | A scalable graph-cut algorithm for n-d grids[END_REF] meant for such cases, where all information is not stored in the memory at all times. This was not required for the experiments presented in this paper though.

Rewriting as a multi-label problem

We show now another point of view on sequence segmentation with growth constraint.

The successive labels L i (t) of a given pixel i over time might change only once, and only from 0 (background) to 1 (foreground object). Hence, this vector of labels L i (t) is of the form (0, 0, . . . , 0, 1, . . . , 1, 1) and can be represented by just the time index τ i of the first 1, i.e. the earliest time at which the pixel starts belonging to the object. This time τ i is in [1, T +1], with T +1 meaning "never". We have thus transformed a binary optimization problem on a sequence of images with shape growth constraint into a multi-label problem defined on one single image, without any constraint. This new problem can be expressed in the Markov Random Field form (4.1) with V i (τ i ) ∶=

t<τ i V t i (0) + t⩾τ i V t i (1)
and

W i,j (τ i , τ j ) ∶= t<min(τ i ,τ j ) W t i,j (0, 0) + τ i ⩽t<τ j W t i,j (1, 0) + τ j ⩽t<τ i W t i,j (0, 1) + t⩾max(τ i ,τ j ) W t i,j (1, 1)
where (V t i ) and (W t i,j ) define the energy E t at time t, and where sets of summation may be empty (note that any time t appears in exactly one summation set only). The submodularity of the binary interaction terms W t in each frame implies the submodularity of the multilabel interaction term W , i.e.

W i,j (τ 1 , τ 2 ) + W i,j (τ ′ 1 , τ ′ 2 ) ⩽ W i,j (τ 1 , τ ′ 2 ) + W i,j (τ ′ 1 , τ 2 )
for all labels satisfying τ 1 ⩽ τ ′ 1 and τ 2 ⩽ τ ′ 2 (proof in Appendix). Thus, this energy can be minimized globally efficiently with now standard techniques (e.g. [START_REF] Darbon | Global optimization for first order Markov random fields with submodular priors[END_REF]). Note that in the particular case where interaction terms W t do not depend on t, the interaction term W of this multi-label energy above can be rewritten as a convex function g of (τ i -τ j ), and then Ishikawa's construction [START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF] can be applied. It turns out that the graph built this way is precisely the graph that we built in our initial binary multi-frame problem.

Our initial formulation is however more flexible, in that interaction terms can depend on t, and more natural, in that inclusion constraints can easily be enforced in spatial or/and time subregions only, while this would not be expressible with the multi-label formulation.

Experimental results

We applied the proposed method to three different applications, in order to validate its use to enforce shape shrinkage or growth, for 2D or 3D image sequences:

(Application 1) We impose shrinkage constraint for a 2D sequence, to segment a melting multiyear ice from a time series of satellite measurements.

(Application 2) We enforce shape growth for a 2D sequence, to segment burned areas from optical satellite images.

(Application 3) We use growth constraint for 3D image sequences, to segment growing brain tumors from multimodal MRI volumes.

The performance of our framework with monodirectional links, [Mono=const], is compared with other graph-cut-based methods:

• [w/o] : Graph cut with no temporal links, i.e. independent segmentation of each frame.

• [Feedforward] : After segmenting the first frame with graph cut approach, foreground/background pixels are marked as seeds with infinite unary costs in the next frame for enforcing shape growth/shrinkage. both directions (from t to t + 1 and from t + 1 to t), with a constant weight (finite or infinite).

• [Bi=variable] : As proposed in [START_REF] Wolz | Measurement of hippocampal atrophy using 4d graph-cut segmentation: Application to ADNI[END_REF], bidirectional temporal links are computed based on intensity differences between pixels in successive image frames, i.e. in the same way as spatial links.

For comparison between the methods we used the Dice score [START_REF] Dice | Measure of the amount of ecological association between species[END_REF], D = (2 M ∩ M ) ( M + M ), where M and M are manually and automatically segmented foreground regions, respectively.

Application 1: melting sea ice in satellite images

The melting of sea ice is correlated to the increases in sea surface temperature and to associated climatic changes. Thus it is very important to monitor sea ice evolution and to develop methods for automated analysis of satellite measurements. Previous works on ice floe segmentation attempted to extract temporal information by using ice percentages, area and shape parameters of the ice floes at the previous time moments as prior information to segment the floe at the next time moment [START_REF] Tarabalka | Shape-constrained segmentation approach for arctic multiyear sea ice floe analysis[END_REF][START_REF] Haverkamp | Using temporal information in an automated classification of summer, marginal ice zone imagery[END_REF]. These feedforward approaches were unable to accurately estimate melting ice profiles because of the low signal-to-noise ratio and of the lack of contrast in satellite data, producing area estimations that are noisy in time, while a multiyear ice floe can actually only melt in shows MODIS and AMSR-E images at two time moments from the considered data set.

The floe was tracked from the AMSR-E data, where multiyear ice has a low microwave emissivity (dark area in Fig. 4.2), and is in this way distinguishable from clouds and younger ice which has a higher emissivity (white area in Fig. 4.2). However, the low spatial resolution of these data does not allow to quantify the ice floes areas accurately. In accordance with the tracking measurements, a time series of T = 75 MODIS and upscaled AMSR-E images with the ice floe was built, with spatial dimensions of 800 × 800 pixels.

We denote each MODIS image by I t , and each upscaled AMSR-E image smoothed by Gaussian filter by A t , t = 1, ..., T .

The objective is to compute T segmentation maps L t , where each pixel (x, y) has label L t (x,y) = 1 if it belongs to the floe at time t, and 0 otherwise. In order to apply the proposed method with a shrinkage constraint to the selected time series, the images must be aligned, so that the property that the floe in the image I t+1 is included in the floe of the image at the previous time moment I t can be expressed directly in terms of pixel locations.

∀t, x, y, L t+1 (x,y) = 1 ⇒ L t (x,y) = 1.
For this purpose, we estimated a reliable region of the foreground (i.e., a region which can be considered with high probability as a part of the floe), R F , and a reliable region of the background, R B , from the AMSR-E images, where a multiyear ice floe is darker than water, young ice and clouds. From the AMSR-E tracking measurements, we derived a foreground seed point and an approximate area of the floe for each time moment t. We then grew a reliable region of the floe R F in each A t from the foreground seed point, until it reached approximately half of the size of the floe. Another region, R B , was grown in A t from the same foreground seed point, but until an area half larger than the approximate We then computed the histograms of the intensities I t of the floe, p t (I F ), and of the background, p t (I B), respectively, and a map of floe probabilities as

p t (F I) = p t (I F )P t (F ) p t (I F )P t (F ) + p t (I B)P t (B) , (4.4)
P t (B) = A t -min x,y A t (x,y) max x,y A t (x,y) -min x,y A t (x,y) , P t (F ) = 1 -P t (B).
The quantities above depend on the pixel location. The images I t were aligned by exhaustive searching over rigid motions (rotations and translations) to maximize the cor-CHAPTER 4. SEGMENTATION WITH SHAPE GROWTH/SHRINKAGE CONSTRAINT 78 relation between maps of foreground probabilities at the current and previous moments.

This search is largely affordable given the low frame rate. We computed potentials and interaction terms between neighboring pixels as:

V t i (1) = -ln[p t (F I)], V t i (0) = -ln[p t (B I)], (4.5)
W t i,j = δ L i ≠L j β exp ⎡ ⎢ ⎢ ⎢ ⎣ - (I t i -I t j ) 2 2σ 2 ⎤ ⎥ ⎥ ⎥ ⎦ , (4.6)
where σ is a standard deviation of I t , β is a parameter that controls the importance of the spatial interaction energy term. We found experimentally that setting β = 2 yields robust results. The proposed method was applied with monodirectional temporal links to enforce floe shrinkage, as described in Sec. constant shape for all time moments, and the Dice score decreases.

When the proposed shrinkage constraint is used instead, the segmentation accuracy increases with w, and w = ∞ yields results with monotonous shrinkage of the shape area (see Fig. 4.4(a)). Moreover, the proposed method copes well with rapid shrinkage events, without underestimating preceding images, or overestimating the event itself at onset.

Another advantage of using monodirectional infinite links is that there are no additional parameters to quantify temporal coherency.

Application 2: growing burned areas in satellite observations

Biomass burning has a significant impact on the Earth's climate system. Satellite remote sensors acquire data for the continuous monitoring of burned areas at both regional and global scales. Thus, there is a need to develop methods for automated fire mapping.

While most of the existing techniques for mapping burned areas analyze temporal evolution of each pixel in an image scene [START_REF] Roy | Burned area mapping using multitemporal moderate spatial resolution data -a bi-directional reflectance model-based expectation approach[END_REF], recent studies have proved the advantage of considering spatial contextual classification for accurate fire classification [START_REF] Lewis | Fire mapping for managers in North Australian savanna; an object orientated classification model for MODIS imagery[END_REF][START_REF] Giglio | An active-fire based burned area mapping algorithm for the modis sensor[END_REF].

Both works [START_REF] Lewis | Fire mapping for managers in North Australian savanna; an object orientated classification model for MODIS imagery[END_REF][START_REF] Giglio | An active-fire based burned area mapping algorithm for the modis sensor[END_REF] map fires from MODIS data by detecting and classifying persistent changes in a daily vegetation-index time series. Giglio et al. [START_REF] Giglio | An active-fire based burned area mapping algorithm for the modis sensor[END_REF] exploit the closest fixed pixel's neighborhood to refine fire classification. Lewis [START_REF] Lewis | Fire mapping for managers in North Australian savanna; an object orientated classification model for MODIS imagery[END_REF] segments and analyzes change detection maps between two consecutive time moments. Manual post-processing is needed to correct classification errors, which are a consequence of either a cloud cover, or low contrast between burned and unburned areas.

In our study we analyzed two time series of Terra MODIS atmospherically-corrected Level 2G daily surface reflectance measurements over the tropical savannas in the Northern Australia ("MOD09GA" product, tile h31v10), each of the data sets being acquired during forty days of the dry season in September-October (days 244-283) of 2011 and 2013, respectively. Wildfires in this region of Australia are frequent and extensive. We used MODIS band 5 (1.240 µm) 500-m land surface reflectance data as they provide the highest burned-unburned separability and are largely insensitive to smoke aerosols [START_REF] Roy | Burned area mapping using multitemporal moderate spatial resolution data -a bi-directional reflectance model-based expectation approach[END_REF].

Each time series comprised a set of T = 40 images with spatial dimensions W × H = 400 × 400 pixels. computing an initial histogram of burned areas and comparing fire maps, respectively.

The MCD64A1 product contains fire classification maps, where each pixel is associated with either an estimated day of burn, or an unburned flag, or an unmapped flag due to insufficient data. These maps are computed by applying the approach from [START_REF] Giglio | An active-fire based burned area mapping algorithm for the modis sensor[END_REF] on two 500-m MODIS channels coupled with 1-km MODIS active fire observations.

We segmented each of the considered image time series [t 1 , ..., t T ] by applying the following iterative procedure (the flowchart is shown in Fig. 4.6):

Initialization: k ∶= 0. The initial training mask of burned areas R B k is built using MCD64A1 product, by selecting the pixels burned during the days [t 1 -D, t 1 -1]. This mask can also be created based on ground observations of the considered area on the day (t 1 -1).

1. The training mask of unburned areas R U k is constructed by dilating R B k with a disk of radius r [START_REF] Soille | Morphological Image Analysis[END_REF] and selecting the complementary of the resulting image. the MODIS band 5 for burned p t (I B) and unburned p t (I U ) areas are computed, using the masks R B k and R U k , respectively. If the data for some pixels is missing, these pixels are not considered when computing histograms.

For a subset of T

′ images t = [t 1 + kT ′ , t 1 + (k + 1)T ′ -1],
3. For the images [t 1 + kT ′ , t 1 + (k + 1)T ′ -1], potentials are computed, assuming equal priors p t (B) = p t (U ) = 1 2: Interaction terms are computed using Eq. 4.6.

V t i (1) = -ln[p t (B I t i )] = -ln p t (I t i B) p t (I t i B) + p t (I t i U ) , (4.7)
V t i (0) = -ln[p t (U I t i )] = -ln p t (I t i U ) p t (I t i B) + p t (I t i U ) . ( 4
4. The graph-cut optimization is applied on a joint graph of the images [t 1 , t 1 + (k + 1)T ′ -1], yielding (k + 1)T ′ segmentation maps.

5. If the whole set of T images is segmented, exit the algorithm. Otherwise: k ∶= k +1.

The segmentation map L t 1 +kT ′ -3 is used as the new training mask of burned areas R B k . We do not use the segmentation result of the last frame (t 1 + kT ′ ), because extreme frames benefit from less information from neighboring frames, and are therefore more subject to segmentation errors. Go to step 1. sulting dice scores for both data sets. We have found that a weak temporal regularization (both bi-/monodirectional) outperforms a segmentation without temporal constraint. Increasing the bidirectional temporal regularization towards high values, however, decreases the performance. On the opposite, introducing monodirectional infinite links to impose shape growth yields the most accurate results. As can be seen from Fig. 4.7 and 4.8(b, d), the proposed approach achieves comparable results with the method [START_REF] Giglio | An active-fire based burned area mapping algorithm for the modis sensor[END_REF] by using only one MODIS channel and no post-processing, while [START_REF] Giglio | An active-fire based burned area mapping algorithm for the modis sensor[END_REF] used two MODIS bands coupled with 1-km MODIS active fire observations and post-processing. Furthermore, the new method copes better with the missing or noisy data thanks to the introduced spatiotemporal graph.

Application 3: growing tumor in 3D medical scans

Glioma is the most frequent primary tumor of the brain. The tumor is known to grow steadily, and lesions are evaluated with respect to volume change in different magnetic resonance image (MRI) modalities. In our experiment we evaluated a set of 760 multimodal image volumes -each comprising T1 MRI, contrast-enhanced T1 MRI (T1c), T2 MRI, and T2 FLAIR MRI -acquired from ten patients initially diagnosed with low grade glioma. The time series have 3-14 time points, with 3-6 months in between any two acquisitions. All image volumes were rigidly registered and three 2D slices intersecting with the tumor center were manually annotated through an expert in every volume, representing an approximate truth. Full 3D segmentations for images of each individual time point were obtained using a generative model for multimodal brain tumor segmentation [START_REF] Menze | A generative model for brain tumor segmentation in multi-modal images[END_REF]. This algorithm models the lesion with a latent atlas class [START_REF] Riklin-Raviv | Segmentation of image ensembles via latent atlases[END_REF] amending the tissue atlas of the standard EM segmenter [START_REF] Williams | Adaptive segmentation of MRI data[END_REF][START_REF] Kapur | Segmentation of brain tissue from magnetic resonance images[END_REF]. We applied it to each multimodal data set of each time point in an independent fashion. The segmentation model delineates the lesion individually in each modality. It assumes that changes of the core (visible in T1c) will occur within the larger edema regions (visible in T2 or FLAIR) and, hence, to only have class transitions from healthy to edema and from edema to core. As the tumor grows steadily, we can assume that negative volume changes stem from imaging artifacts, such as local intensity changes, a common problem in MRI. To this end we model the tumor to be either stable, in this case regularizing the segmentation along time and suppressing noise, or to expand in volume between any two time points.

We identified the foreground label F with tumor (edema and core) and background B with healthy tissue. Then the potential V s,t i (L s,t i ) of label L s,t i at voxel i, time point t, and imaging sequence s is equal to V s,t i (0) = p s,t (F I s,t ) and V s,t i (1) = p s,t (B I s,t ) = 1 -p s,t (F I s,t ). The tumor probabilities were calculated from image volumes I using the generative model [START_REF] Menze | A generative model for brain tumor segmentation in multi-modal images[END_REF], and we identified tumor subclasses with p(F I s=T 1,t ) for core, and p(F I s=T 2,t ) with edema. We modeled the 3D spatial constraints through a 26 neighbor-Figure 4.9: Two time series of T2 and FLAIR MR image volumes, each of which acquired about 3-6 months apart. The left case is rapidly growing between the second and fourth scene, the right case displays intensity modifications in the last scene, leading to a suboptimal performance of the initial multimodal segmentation (yellow) [START_REF] Menze | A generative model for brain tumor segmentation in multi-modal images[END_REF]. The proposed multi-temporal segmentation with growth constraints (green) delineates areas similar to the manual evaluation (magenta), being more robust against intensity variations of the MR images. It does not smooth out outlines of rapidly growing tumors as conventional bi-directional temporal constraints would do. hood (N ) linking the central voxel with all its immediate neighbors. Interaction terms W s,t i,j (L s,t i , L s,t j ) between neighboring voxels in each sequence s ∈ [T 1, T 2] are computed from the channel-specific intensity differences as

δ L i ≠L j β α(i, j) α tot exp ⎛ ⎝ - I s,t (i) -I s,t (j) A 2 ⎞ ⎠
with β = 0.5, α(p, q) = 1 distance(p,q) , α tot = ∑ q∈N (pixel p) α(p, q) and A = 1 3 max I s,tmin I s,t . We impose growth constraint in 3D+t as explained in Sec. 4.2.1, and inclusion constraints as in Sec. 4.2.2 : the foregrounds in T1 and T1c modalities are required to be included in the one of T2, which is included in the one of FLAIR.

In our test we first segmented the images from each time point independently, and calculated Dice scores averaged over all images of a time series as an estimate of the baseline performance. Then we tested different regularizations in time, i.e., [w/o], [Mono] and [Bi], calculate Dice scores, and compared them against the baseline results from the unregularized segmentations. Fig. 4.9 shows results for two exemplary time series, and Fig. 4.10(a) reports differences between Dice scores of regularized segmentations and baseline segmentations. In this experiment "weak" regularizations refer to small regularization parameters (w ≪ 1), while "strong" regularizations refer to w ≫ 1, representing the infinite mono-and bi-directional links. Fig. 4.10(a) shows results for all ten time series.

As for the previous two applications, both bi-/monodirectional temporal regularization with low values of w yields better accuracies when compared to the results without tem- right) shows the exponential growth of the tumor core (dotted lines) that is associated with this disease [15], for this patient indicating a rapid tumor progression starting at about day 500. Longitudinal image segmentation, as obtained for the ten time series in our test, can be further analyzed, for diagnosis and treatment monitoring, e.g., through algorithms estimating the speed of the tumor outlines under anatomical constraints [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF].

Extension of the current 5D segmentation could integrate this speed estimation, or extend the multimodal tumor segmentation [START_REF] Menze | A generative model for brain tumor segmentation in multi-modal images[END_REF] for longitudinal data sets.

Conclusion and discussion

In this chapter, we addressed the problem of shape segmentation in 2D and 3D sequences of very noisy/low-contrast images, where shapes monotonously grow or shrink in time. In order to enforce shape growth or shrinkage, we proposed a new graph-cut-based method for computing the globally-optimal spatio-temporal segmentation satisfying that constraint. The main idea was to introduce monodirectional infinite links between pixels at the same spatial locations in successive image frames, which prohibit a shape to shrink or grow over time, and then to perform a graph cut optimization on the constructed graph.

The limitation of the proposed method is that it can be applied only to a time series of images on the same scale and perfectly aligned with respect to the foreground object, so that in the case of shape growth the foreground at the moment (t + 1) contains all the foreground pixels at the previous moment t, sharing the same spatial locations. Thus, if the foreground object moves over time, images must be aligned before applying the new graph-cut-based technique. We also demonstrated the possibility to impose intersequences object inclusion constraints by adding directed infinite links to the joint graph associated to all sequences.

We validated the performance of the proposed approach for the segmentation of growing (burned areas) and shrinking (ice floe) shapes from 2D time series of satellite images, and for the segmentation of growing 3D tumor volumes from MR sequences. The method proved to be robust to important noise and low contrast, and to cope well with missing data. Moreover, it showed linear complexity in practice, so that globally optimal shape-consistent segmentations of image time series are obtained in a matter of seconds.

Another limitation of the presented method is that it is designed and tuned for the specific shape prior, and is thus not generic. Furthermore, even though the complexity of the method is linear, it is not easy to parallelise it, therefore its scalability is limited.

During the last decade, the use of deep convolutional neural networks (CNNs) for machine learning tasks has been growing in the vision community. CNNs are able to automatically learn hundreds of task-specific features, which may allow us to easily scale the pool of contextual features considered. Moreover, they are able to incorporate high-level features in terms of simple operations that are highly parallelized in GPUs. With the growing availability of high volumes of satellite images, it is clear that highly parallel approaches will gain attention in the near future. In the next chapters I will describe our research outcomes on the use of the deep neural network-based methods for remote sensing image classification.

Chapter 5 Deep Learning for Large-Scale Image Classification

There is a growing interest in analyzing remote sensing data at a large scale, i.e. from the scale of cities to the world-scale [START_REF] Lassalle | A scalable tile-based framework for region-merging segmentation[END_REF][START_REF] Mattyus | Enhancing road maps by parsing aerial images around the world[END_REF][START_REF] Michel | Stable mean-shift algorithm and its application to the segmentation of arbitrarily large remote sensing images[END_REF]. This implies several challenges, such as the computational complexity of the methods and intra-class variability. Due to the large spatial extent of the large-scale datasets, the objects belonging to the same semantic class may have considerably different appearance. To address this challenge and conduct an accurate classification, the high-level contextual information must be taken into consideration, such as the shape of objects and their location with respect to other objects.

There is an incipient trend both in the vision and remote sensing literature to design classifiers that automatically learn expressive high-level contextual features. In this context, deep learning and, in particular, convolutional neural networks (CNNs), have gained significant attention in the image analysis community over the last years. These have been originally devised for the image categorization problem, i.e., the assignment of one label to an entire image. For example, they have been used to categorize objects in natural scenes (e.g., airplane, bird, person) or land use in the case of aerial images (e.g., forest, beach, tennis court). CNNs jointly learn to extract relevant contextual features and conduct the categorization. In addition to the suppression of the feature design process, which is an interesting advantage itself, this technique has consistently beaten alternative methods in a wide range of problems [START_REF] Matthew | Visualizing and understanding convolutional networks[END_REF]. Nowadays, one can reasonably expect to find CNN-based techniques scoring the best positions in the leaderboards of online image-related contests.

One of the challenges to apply CNN-based methods for remote sensing image classification consists in obtaining fine-grained high-resolution classification maps. This is due to two main reasons: a) There is a structural limitation of CNNs to carry out fine-grained classification.

If we wish to keep a low number of learnable parameters, the ability to learn long-range contextual features comes at the cost of losing spatial accuracy, i.e., a trade-off between detection and localization. This is a well-known issue and still a scientific challenge [START_REF] Chen | Semantic image segmentation with deep convolutional nets and fully connected CRFs[END_REF][START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF]. b) In the specific context of remote sensing imagery, there is a significant lack of spatially accurate reference data for training. For example, the OpenStreetMap collaborative database provides large amounts of free-access maps over the Earth, but irregular misregistrations and omissions are frequent all over the dataset. In such circumstances, CNNs cannot do better than learning rough estimates of the objects' locations, given that the boundaries are hardly located on real edges in the training set.

In this chapter, I present two approaches to address this challenge. One approach consists in designing new architectures specifically intended to provide high-resolution outputs. We conduct an in-depth analysis of the recent dense semantic labelling CNNs to establish the desired properties of an ideal semantic labeling CNN, and assess how the existing methods stand with regards to these properties. Out of these observations, we then derive a CNN framework specifically adapted to the high-resolution dense classification problem. In addition to learning features at different resolutions, it learns how to combine these features. By integrating local and global information in an efficient and flexible manner, it outperforms previous techniques when evaluated on public benchmarks of high-resolution aerial image labeling. I also present an aerial image labeling dataset we have developed to assess the capability of CNNs to generalize to different geographic regions.

Another approach is to use first the base CNN as a rough classifier of the objects' locations, and then process this classification using the original image as guidance, so that the output objects better align to real image edges. Different iterative enhancement algorithms have been presented in the literature to progressively improve the coarse CNN outputs, seeking to sharpen object boundaries around real image edges. However, one must carefully design, choose and tune such algorithms. Instead, we have proposed to directly learn the iterative process itself. For this, we formulate a generic iterative enhancement process inspired from partial differential equations, and observe that it can be expressed as a recurrent neural network (RNN). Consequently, we train such a network from manually labeled data with the purpose to correct image classification maps. The designed RNN automatically discovers relevant data-dependent features to enhance the classification and the specifics of an iterative process to do so.

High-resolution semantic labeling with convolutional neural networks

Background on convolutional neural networks

An artificial neural network is a system of interconnected neurons that pass messages to each other. When the messages are passed from one neuron to the next one without ever going back (i.e., the graph of message passing is acyclic), the network is referred to as feed-forward [14], which is the most common type of network in image categorization. An individual neuron takes a vector of inputs x = x 1 . . . x n and performs a simple operation to produce an output a. The most common neuron is defined as follows:

a = σ(wx + b), (5.1)
where w denotes a weight vector, b a scalar known as bias and σ an activation function.

The weights w and biases b are the parameters of the neurons that define the function.

The goal of training is to find the optimal values for these parameters, so that the function computed by the neural network performs the best on the task assigned.

The most common activation functions σ are sigmoids, hyperbolic tangents and rectified linear units (ReLU). For image analysis, ReLUs have become the most popular choice due to some practical advantages at training time, but novel activation functions have been recently proposed as well [START_REF] Clevert | Fast and accurate deep network learning by exponential linear units (elus)[END_REF].

Instead of directly connecting a huge set of neurons to the input, it is common to organize them in groups of stacked layers that transform the outputs of the previous layer and feed it to the next layer. This enforces the networks to learn hierarchical features, performing low-level reasoning in the first layers (such as edge detection) and higher-level tasks in the last layers (e.g. [START_REF] Felzenszwalb | Object detection with discriminatively trained part-based models[END_REF], assembling object parts). For this reason, the first and last layers are often referred to as lower and upper layers, respectively.

In an image categorization problem, the input of the network is an image (or a set of features derived from an image), and the goal is to predict the correct category of the entire image. We can view the pixelwise semantic labeling problem as taking an image patch and categorizing its central pixel. Finding the optimal neural network classifier reduces to finding the weights and biases that minimize a loss L between the predicted labels and the target labels in a training set. Let Ω be the set of possible semantic classes.

Labels are typically encoded as a vector of length Ω with values '1' at the position of the correct label and '0' elsewhere. The network contains thus as many output neurons as possible labels. A softmax normalization is performed on top of the last layer to guarantee that the output is a probability distribution, i.e. the label values are between zero and one and sum to one. The multi-label problem is then seen as a regression on the desired output label vectors.

The loss function L quantifies the misclassification by comparing the target label vectors y (i) and the predicted label vectors ŷ(i) , for p training samples i = 1 . . . p. In this work we use the common cross-entropy loss, defined as:

L = - 1 p p i=1 Ω k=1 y (i) k log ŷ(i) k . (5.2)
Training neural networks by optimizing this criterion converges faster, compared with, for instance, the Euclidean distance between y and ŷ. In addition, it is numerically stable when coupled with softmax normalization [14].

Note that in the special case of binary labeling we can produce only one output (with targets '1' for positive and '0' for negative). In this case a sigmoid normalization and cross-entropy loss are analogously used, although a multi-class framework can also be used for two classes.

Once the loss function is defined, the parameters (weights and biases) that minimize the loss are found via gradient descent, by computing the derivative ∂L ∂w i of the loss function with respect to every parameter w i , and updating the parameters with a learning rate λ as follows:

w i ← w i -λ ∂L ∂w i . (5.
3)

The derivatives ∂L ∂w i are obtained by backpropagation, which consists in explicitly computing the derivatives of the loss with respect to the last layer's parameters and using the chain rule to recursively compute the derivatives of each layer's outputs with respect to its weights and inputs (the inputs being the previous layer's outputs). In practice, instead of averaging over the full dataset, the loss (5.2) is estimated from a random small subset of the training set, referred to as mini-batch. This learning technique is known as stochastic gradient descent.

Convolutional neural networks (CNNs) are a particular type of neural network. The overall success of CNNs lies mostly in the fact the the networks are forced by construction to learn hierarchical contextual translation-invariant features, which is particularly useful in the context of image analysis.

CNNs contain so-called convolutional layers, a specific type of layer that imposes a number of restrictions compared to a more general fully connected layer, where every neuron is connected to all outputs of the previous layer [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF]. Every neuron is associated to a spatial location (i, j) in the input image (see Fig. 5.1). The output a ij associated with location (i, j) in a convolutional layer is computed as:

a ij = σ((W * X) ij + b), (5.4)
where W denotes a kernel with learned weights, X the input to the layer and ' * ' the convolution operation. Notice that this is a special case of the neuron in Eq. 5.1 with the following constraints:

• The connections only extend to a limited spatial neighborhood determined by the kernel size;

• The same filter is applied to each location, guaranteeing translation invariance.

Multiple convolution kernels are usually learned in every layer, interpreted as a set of spatial feature detectors. The responses to every learned filter are thus referred to as feature maps. Note that the convolution kernels are actually three-dimensional: in addition to their spatial extent (2D), they span along all the feature maps in the previous layer (or eventually through all the bands in the input image). As this third dimension can be inferred from the previous layer it is rarely mentioned in the architecture descriptions.

Compared to the fully connected layer, a convolutional layer highly reduces the number of parameters by enforcing the aforementioned constraints. This results in an easier optimization problem, without losing much generality. This opened the door to using the image itself as an input without any feature design and selection process, as CNNs discover the relevant spatial features to conduct classification.

Increasing the receptive field

In addition to convolutional layers, state-of-the-art networks such as Imagenet [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] involve some degree of downsampling, i.e., a reduction in the resolution of the feature maps with the purpose to increase the so-called receptive field of the neurons. The receptive field denotes the spatial extent of the input image connected to a certain neuron, possibly indirectly through other neurons in previous layers: it is the set of pixels on which a neuron depends. In other words, it quantifies how far a neuron can "see" in the image. In most applications, a large amount of spatial context must be taken into account in order to successfully label the images. For example, to deduce that a certain pixel belongs to a rooftop, it might not be enough to just consider its individual spectrum: we might need to observe a large patch around this pixel, taking into account geometry and structure of the objects, to infer its correct class.

Neural networks for image analysis should thus be designed to accumulate, through their layers, a large enough receptive field. While a straightforward way to do it is to use large convolution kernels, this is not a common practice mostly due to its computational complexity. Besides, this would aim at learning large filters all at once, with millions of parameters. It is preferable to learn a hierarchy of small filters instead, reducing the number of parameters while remaining expressive, and thus making the optimization problem easier.

The most common approach to reduce the number of parameters for a given receptive field size is to downsample the feature maps throughout the network. This is commonly achieved by interleaving downsampling layers with convolutional layers. This way, the resolution of the feature maps gets lower and lower as we traverse the layers from input to output.

To downsample the feature maps, the most popular approach is to use the so-called max pooling layer [17]. A max pooling layer takes a group of neighbors in the feature map and condenses them into a single output by computing the maximum of all incoming activations in the window. The pooling windows in general do not overlap, hence the output map is downsampled (see Fig. 5.2). For instance, if pooling is performed in a 2 × 2 window, the feature map is reduced to half of its resolution.

Computing the maximum value is inspired by the idea of detecting objects from their parts. For example, in a face detector it is important to identify the constituents of a face, such as hair or nose, while the exact locations of these components should not be such a determinant factor. The max pooling layer conveys then to which extent there is evidence of the existence of a feature in a vicinity. Other less popular forms of downsampling include average pooling and applying convolutions with a stride, i.e., "skipping" some of them (e.g., applying every other convolution).

Pooling operations (and downsampling in general) hard-code robustness to spatial deformations, a virtue that boosted the success of CNNs for image categorization. However, spatial precision is lost when downsampling. The increased receptive field (and thus recognition capability) comes at the price of losing localization capability. This well-reported trade-off [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF][START_REF] Chen | Semantic image segmentation with deep convolutional nets and fully connected CRFs[END_REF] is a major concern for dense labeling.

We could still imagine a downsampling network that preserves localization: it would learn features of the type "a corner at the center of the receptive field", "a corner one pixel left of the center of the receptive field", "a corner two pixels left of the center of the receptive field", and so on, multiplying the number of features to be learned. This would however discredit the use of downsampling to gain robustness to spatial variation in the first place. The recognition/localization trade-off must thus be properly addressed to design a high-resolution semantic labeling network.

Overview of high-resolution labeling CNNs

Fully convolutional networks (FCNs) [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF] have become the standard in semantic labeling, including in remote sensing image analysis [START_REF] Volpi | Dense semantic labeling of sub-decimeter resolution images with convolutional neural networks[END_REF][START_REF] Kampffmeyer | Semantic segmentation of small objects and modeling of uncertainty in urban remote sensing images using deep convolutional neural networks[END_REF]. They contain only convolutional layers, i.e., no fully connected layers. Therefore, they can be applied to images with various sizes: inputting a larger image patch produces a larger output, the convolutions being performed on more locations. When an FCN has downsampling layers, the output contains fewer elements than the input, since the resolution has been decreased. This gave birth to the so-called deconvolutional (or upconvolutional) layer, which upsamples a feature map by interpolating neighboring elements (as the last layer in Fig. 5.2). However, simply adding a deconvolutional layer to upsample the output on top of a network provides dense outputs but imprecise labeling results, because the upsampling is performed in a naive way from the coarse classification. This is dissatisfying in many applications, such as high-resolution aerial image classification, where the goal is to precisely identify and outline tiny objects such as cars. The open question is then which type of FCN would be able to conduct fine predictions that provide detailed high-resolution outputs, while still taking large amounts of context into account and without exploding the number of trainable parameters.

We now describe what we consider to be the elementary principle from which to derive efficient dense classification architectures. Let us then first observe that while our goal is to take large amounts of context into account, we do not need this context at the same spatial resolution everywhere. Suppose we want to classify the central pixel of the patch in Fig. 5.3(a). Such a gray pixel, taken out of context, could be easily confused with an asphalt road. Considering the whole patch at once helps to infer that the pixel belongs indeed to a gray rooftop. However, two significant issues arise if we take a fullresolution large patch for context: a) it requires many computational resources that are actually not needed for an effective labeling, and b) it does not provide robustness to spatial variation (we might actually not care about the exact location of certain features to determine the class). Conducting predictions from low-resolution patches instead is not a solution as it produces inaccurate coarse classification maps. Nevertheless, it is actually not necessary to observe all surrounding pixels at full resolution: the farther we go from the pixel we want to label, the lower the requirement to know the exact location of the objects. For example, in the patch of Fig. 5.3(b) it is still possible to classify the central pixel, despite the outer pixels being blurry. Therefore, we argue that a combination of reasoning at different resolutions is necessary to conduct fine labeling, if we wish to take a large context into account in an efficient manner.

In the following, we analyze the main families of high-resolution classification networks that have been proposed in the past two years. For each of them we discuss the following aspects:

• How a solution to the fine-grained labeling problem is provided;

• Where this solution stands with respect to the principle of Fig.

• General advantages and disadvantages, and computational efficiency.

Dilation Networks

Dilation networks are based on the shift-and-stitch approach or à trous algorithm [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF].

This consists in conducting a prediction at different offsets to produce multiple lowresolution outputs, which are then interleaved to compose the final high-resolution result.

For example, if the downsampling factor of a network is S, one should produce S 2 classification maps by shifting the input horizontally and vertically. Such an interleaving can also be implemented directly in the architecture, by using "dilated" operations [START_REF] Yu | Multi-scale context aggregation by dilated convolutions[END_REF],

i.e., performing them on non-contiguous elements of the previous feature maps. This principle is illustrated in Fig. 5.4.

Dilations have been used with two purposes:

1. As an alternative to upsampling for generating full-resolution outputs [START_REF] Dubrovina | Computational mammography using deep neural networks[END_REF][START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF].

2. As a means to increase the receptive field [START_REF] Chen | Semantic image segmentation with deep convolutional nets and fully connected CRFs[END_REF][START_REF] Yu | Multi-scale context aggregation by dilated convolutions[END_REF], by enlarging the area covered by a convolution kernel, without increasing the number of trainable parameters.

Regarding the first point, we must mention that there is no theoretical improvement compared to an FCN with naive upsampling, because the presence of pooling layers still reduces spatial precision. Executing the prediction multiple times at small offsets still keeps predictions spatially imprecise.

Regarding the second point, we must remark that while dilated convolutions increase the receptive field, this does not introduce robustness to spatial variation per se. For example, a network with only dilated convolution layers would have a large receptive field but would only be able to learn filters of the type "a building in the center, with a car exactly five pixels to the left". This robustness would have to be thus learned, hopefully, by using a larger number of filters.

The use of an interleaved architecture at training time, implemented with dilations, has been however reported to be beneficial. In the context of aerial image labeling, Sherrah [START_REF] Sherrah | Fully convolutional networks for dense semantic labelling of highresolution aerial imagery[END_REF] recently showed that it outperforms its FCN/upsampling counterpart1 .

The major improvement compared to the FCN/upsampling network was measured in the labeling capabilities of the car class, which is a minority class with tiny objects, difficult to recognize [START_REF] Volpi | Dense semantic labeling of sub-decimeter resolution images with convolutional neural networks[END_REF]. While the dilation strategy is not substantially different from an architectural point of view compared to naive upsampling, some advantages in training might explain the better results: In the upsampling case the network is encouraged to provide a coarse classification that, once upsampled, is close to the ground truth. In the dilation network, on the contrary, the interleaved outputs are directly compared to individual pixels in the ground truth, one by one. The latter seems to better avoid suboptimal solutions that absorb minority classes or tiny objects.

The computational time and memory required by dilation networks are significant, to the point that using GPUs might become impractical even with moderately large architectures. This is because the whole network rationale is applied to many contiguous locations.

Overall, while dilation networks have been reported to exhibit certain advantages, they are computationally demanding and do not particularly address the principle of

Deconvolution Networks (unpooling)

Instead of naively upsampling the classification score maps with one deconvolutional layer, a more advanced approach is to attach a multi-layer network to learn a complex upsampling function. This idea was simultaneously presented by different research groups [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for robust semantic pixel-wise labelling[END_REF][START_REF] Noh | Learning deconvolution network for semantic segmentation[END_REF] and later extended to different problems (e.g., [START_REF] Yang | Object contour detection with a fully convolutional encoder-decoder network[END_REF]). A simple way to implement this idea is to "reflect" an existent FCN, with the same number of layers and The other elements in the unpooling window are set to zero, leading to sparse feature maps, as illustrated in Fig. 5.6. Unpooling was first introduced as part of a framework to analyze and visualize CNN features [START_REF] Matthew | Visualizing and understanding convolutional networks[END_REF]. The arrows in Fig. 5.5 represent the communication of the pooling indices from the pooling layer to the unpooling layer. In the case of average pooling, the corresponding unpooling layer simply outputs at every location the input activation value divided by the number of elements in the target unpooling window(see Fig. 5.6). In this case, there is no need to transmit a location from pooling to unpooling.

This concept can be thought of as an "encoder-decoder", where the middle layer is seen as a common representation to images and classification maps, while the "encoder"

and "decoder" ensure the translation between this representation and the two modalities.

When converting an FCN to a deconvolution network, the final classification layer of the FCN is usually dropped before reflecting the architecture. This way the interface between the encoder and the decoder is a rich representation with multiple features. The first layer of the encoder takes as input as many channels as the ones in the input image, and the last layer of the decoder produces as many feature maps as the number of classes required. In [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for robust semantic pixel-wise labelling[END_REF][START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for image segmentation[END_REF], alternatively, the network outputs a larger set of features that are then classified with additional layers.

While pooling is used to add robustness to spatial deformation, the fact of "remembering" the location of the max activation helps to precisely locate objects in the deconvolution steps. For example, the exact location of a road might be irrelevant to do any higher-level reasoning later on, but once the network decides to label the road as a semantic object we need to recover the location information to outline it with high precision. This illustrates how deconvolution networks balance the localization/recognition trade-off.

Note however that if one happens not to choose max pooling for downsampling, then the unpooling scheme is not able to recover per se the lost spatial resolution. There is no memory about the location of the higher-resolution feature. Even though max pooling is very common, it has been shown that average or other types of pooling might be more effective in certain applications [17]. In fact, recent results [START_REF] Tobias Springenberg | Striving for simplicity: The all convolutional net[END_REF] suggested that max pooling can be emulated with a strided convolution and achieve similar performance.

The deconvolution network idea is however leveraged when max pooling is chosen for downsampling.

This certainly does not mean that a deconvolutional network is incapable of learning without max pooling layers. Convolution/deconvolution architectures without max pooling have been successfully used in different domains [START_REF] Volpi | Dense semantic labeling of sub-decimeter resolution images with convolutional neural networks[END_REF][START_REF] Simo-Serra | Learning to simplify: fully convolutional networks for rough sketch cleanup[END_REF]. For example, a recent submission to the ISPRS Semantic Labeling Challenge [START_REF] Volpi | Dense semantic labeling of sub-decimeter resolution images with convolutional neural networks[END_REF] is such type of network. The recognition/localization trade-off is not really alleviated in this case: the encoder should encode features of the type "an object boundary 5 (or 7, 10...) pixels away to the left", so that the decoder can really leverage this information and reconstruct a high-resolution classification map.

The depth of deconvolution networks is significantly larger, roughly twice the one of the associated FCN. This often implies a slower and more difficult optimization, due to the increase in the number of trainable parameters introduced by deconvolutional layers.

While the decoding part of the network can be simplified [START_REF] Paszke | Enet: A deep neural network architecture for real-time semantic segmentation[END_REF], this adds arbitrariness to the design, since instead of just reflecting the encoder we must also design the decoder.

To conclude, the deconvolution scheme does address the recognition/localization trade-off, but only in the case where max pooling is used for downsampling. The in-creased network depth can be a concern for an effective training.

Skip Networks

In the original paper about fully convolutional networks, Long et al. [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF] proposed the so-called "skip" architecture to generate high-resolution classification outputs. The idea is to build the final classification map by combining multiple classification maps, obtained from intermediate features of the network at different resolutions (and not just the last one).

The last layer of an FCN outputs as many feature maps as classes, which are interpreted as score or "heat" maps for every class. Intermediate layers, however, tend to have many more features than the number of classes. Therefore, skip networks add extra layers that convert the arbitrarily large number of features of intermediate layers into the desired number of heat maps. This approach allows us to extract multiple score maps for each class from a single network, at different resolutions. The lower-level score maps are fine but have a small receptive field, while the higher-level ones can see farther but with less detail. The score maps are then combined pairwise, from the lower scales to the higher scales. At every step, the lower-resolution score maps are upsampled to match the higher-resolution ones. They are then added elementwise. This is repeated until all intermediate maps are processed. A skip network is illustrated in Fig. 5.7.

Skip networks address the trade-off between localization and precision quite explicitly: the information at different resolutions is extracted and combined. The original paper introduces this methodology as "combining what and where". This approach is closer to the principle described in Fig. 5.3 than the previous approaches reviewed above. The skip network mixes observations at different resolutions, without unnecessarily increasing the depth or width of the architecture (as in deconvolution and dilation networks respectively) and it does not impose a particular type of downsampling (as in deconvolution networks).

While the idea of extracting different resolutions is certainly very relevant, the skip model seems to be inflexible and arbitrary in how to combine them. First of all, it combines classification verdicts, instead of a rich set of features, coming from each of the resolutions. For example, it combines how a layer evaluates that an object is a building by using low-level information, with how another layer evaluates whether the same object is a building by using higher-level information. Let us recall that we use deep multi-layer schemes with downsampling because we actually consider that certain objects can only be detected at the upper layers of the network, when a large amount of context has been taken into account and at a high level of abstraction. It seems thus contradictory to try to refine the boundaries of an object detected at a high level, by using a classification conducted at a lower level, where the object might not be detected at all. Moreover, the element-wise addition restricts the combination of resolutions to be simply a linear combination. The skip links to combine resolutions are in fact parameterless (besides the addition of the scoring layers). We could certainly imagine classes that require a more complex nonlinear combination of high-and low-level information to be effectively classified.

It is worth noting that the creation of the intermediate score maps has also been referred to as a dimensionality reduction step [START_REF] Badrinarayanan | Segnet: A deep convolutional encoder-decoder architecture for image segmentation[END_REF]. It is however not by chance that the amount of reduced features coincides with the amount of classes: even though it is technically a dimensionality reduction, its spirit is to create a partial classification, not just to reduce the number of features. This is confirmed by the name of these layers in the public implementation by Long et al. [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF]: "score" layers. Moreover, if this operation were indeed intended to be just a reduction of dimensionality, we could imagine outputting different amounts of feature maps from different resolutions. However, in that case there would be no way of adding them element by element as suggested. Hariharan et al. [START_REF] Hariharan | Hypercolumns for object segmentation and fine-grained localization[END_REF] also extract intermediate features ("skip" links) from a CNN. Additional convolutional layers are used to derive a fixed number of features from each of these intermediate layers,

and then added to give the final output.

To conclude, the skip network architecture provides an efficient solution to address the localization/recognition trade-off, yet this could be done in a more flexible way that enables a more complex combination of the features.

Learning to combine resolutions

In this section I describe our proposed alternative scheme for high-resolution labeling, derived as a natural consequence of the observations described hereabove. In particular, this architecture leverages the benefits of the skip network while addressing its potential limitations.

Taking multiple intermediate features at different resolutions and combining them seems to be a sensible approach to specifically address the localization/recognition tradeoff, as done with skip networks. In such a scheme, the high-resolution features have a small receptive field, while the low-resolution ones have a wider receptive field. Combining them constitutes indeed an efficient use of resources, since we do not actually need the high-resolution filters to have a wide receptive field, following the principle of Fig. The skip network combines predictions derived from the different resolutions, i.e., score maps for each of the classes. For example, we try to refine the "blobby" building outputted by the coarse classifier, via a higher-resolution classification. However, it is unclear how effectively the higher-resolution classifier detects such building, considering its reduced receptive field and shallow reasoning.

Concatenate Learn to combine features

We thus argue that a more relevant way of performing fine semantic labeling is to combine features, not classification maps. For example, to refine the boundaries of a coarse building, we would use high-resolution edge detectors and not high-resolution building detectors.

In our proposed scheme, intermediate features are extracted from the network and treated equally, creating a pool of features that emanate from different resolutions. A neural network then learns how to combine these features to give the final classification verdict. This adds flexibility to learn more complex relations between the different resolutions and generalizes the element-wise addition of the skip architecture.

The overall process is depicted in Fig. 5.8. First, a subset of intermediate features are extracted from the network. These are naively upsampled to match the resolution of the highest-resolution features. They are then concatenated to create the pool of features.

Notice that while the spatial dimensions of the feature maps are all the same, they originally come from different resolutions. This way, the variation of the feature responses across space will be smoother in certain maps and sharper in others. Note that while it is practical to store in memory the upsampled responses, this is not intrinsically necessary.

For example, we could imagine a system that answers to a high-resolution query by outputting the nearest neighbor in the coarser map or by interpolating neighboring values on the fly.

From the pool of features, a neural network predicts the final classification map (we could certainly use other classifiers, but this lets us train the system end to end). We assume that all the spatial reasoning has been conveyed in the features computed by the initial CNN. This is why we operate on a pixel-by-pixel basis to combine the features.

Any need to look at neighbors should be expressed in the spatial filters of the CNN. This way we conceptually and architecturally separate the extraction of spatial features from their combination.

We can think of the multi-layer perceptron (MLP) with one hidden layer and a non-linear activation function as a minimal system to learn how to combine the pool of features. Such MLPs can learn to approximate any function and, since we do not have any particular constraints, it seems an appropriate choice. In practice, this is implemented as a succession of convolutional layers with 1 × 1 kernels, since we want the same MLP to be applied at every location. By introducing the MLP and executing it at a fine resolution, we must expect an overhead in processing time compared to the skip network.

The proposed technique is intended to learn how to combine information at different resolutions, not how to upsample a low-resolution classification. An example of the type of relation conveyed by this scheme is as follows: "label a pixel as building if it is red and belongs to a larger red rectangular structure, which is surrounded by areas of green vegetation and near a road".

Finally, let us discuss the CNN from which features are extracted (the topmost part of This assumes that the higher-level features can be derived from the lower-level ones. It is basically a part-based model [START_REF] Felzenszwalb | Object detection with discriminatively trained part-based models[END_REF], where we consider that an object can be detected by its parts, and we are using those same parts as the higher-resolution features inputted to the MLP. This seems to be a sensible assumption, yet we must mention that we could eventually think of separate networks to detect features at different resolutions instead of extracting intermediate representations of a single network (as, e.g., in [START_REF] Farabet | Learning hierarchical features for scene labeling[END_REF]). While we adopt the model of Fig. 5.8 in this work, the alternative could be also considered.

It would be certainly interesting to study to which extent it is redundant to learn the features in separate networks and, conversely, how results could be eventually improved by doing it.

Experiments on Potsdam and Vaihingen datasets

We evaluate the aforementioned architectures on two benchmarks of aerial image labeling:

Vaihingen and Potsdam, provided by Commission III of the ISPRS [START_REF]ISPRS. 2D semantic labeling contest[END_REF]. The Vaihingen dataset is composed of 33 image tiles (of average size 2494 × 2064), out of which 16 are fully annotated with class labels. The spatial resolution is 9 cm. Near infrared (NIR), red (R) and green (G) bands are provided, as well as a digital surface model (DSM), normalized and distributed by [START_REF] Gerke | Use of the stair vision library within the ISPRS 2d semantic labeling benchmark (vaihingen)[END_REF]. We select 5 images for validation (IDs: 11, 15, 28, 30, 34) and the remaining 11 images for training, following [START_REF] Sherrah | Fully convolutional networks for dense semantic labelling of highresolution aerial imagery[END_REF][START_REF] Volpi | Dense semantic labeling of sub-decimeter resolution images with convolutional neural networks[END_REF][START_REF] Paisitkriangkrai | Effective semantic pixel labelling with convolutional networks and conditional random fields[END_REF].

Potsdam dataset consists of 38 tiles of size 6000 × 6000 at a spatial resolution of 5 cm, out of which 24 are annotated. It provides an additional blue channel (which we here exclude for simplicity) and the normalized DSM. We select the same 7 validation tiles as in [START_REF] Sherrah | Fully convolutional networks for dense semantic labelling of highresolution aerial imagery[END_REF] (IDs: 2_11, 2_12, 4_10, 5_11, 6_7, 7_8 7_10) and the remaining 17 tiles for training. Both datasets are labeled into the following six classes: impervious surface, building, low vegetation, tree, car and clutter/background.

In order to account for labeling mistakes, another version of the ground truth with eroded boundaries is provided, on which accuracy is measured. To evaluate the overall performance, overall accuracy is used, i.e., the percentage of correctly classified pixels.

To evaluate class-specific performance, the F1-score is used, computed as the harmonic mean between precision and recall [START_REF] Csurka | What is a good evaluation measure for semantic segmentation?[END_REF]. We also include the mean F1 measure among classes, since overall accuracy tends to be too insensitive to minority classes in imbalanced datasets.

Network architectures

To conduct our experiments we depart from a base fully convolutional network (FCN)

and derive other architectures from it. Table 5.1 summarizes our base FCN for the Vaihingen dataset. The architecture is borrowed from [START_REF] Kampffmeyer | Semantic segmentation of small objects and modeling of uncertainty in urban remote sensing images using deep convolutional neural networks[END_REF], except for the fact that we increased the size of the filters from 3 to 5 in the first layer, since it is a common practice to use larger filters if there is a stride. Every convolutional layer (except the last one) is followed by a batch normalization layer [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] and a ReLU activation. We did not optimize the architecture of the base FCN. Padding refers to the amount of zero-valued pixels added around the input to a convolutional layer, so as to preserve the dimension of the feature maps (otherwise the convolution is not applied near the border and the maps become slightly smaller).

The total downsampling factor is 16, out of which 8 is the result of the max pooling layers and 2 of the stride in the first layer. The conversion of the last set of features to classification maps (the "score" layer) is performed by a 1 × 1 convolution. To produce a dense pixel labeling we must add a deconvolutional layer to upsample the predictions by a factor of 16, thus bringing them back to the original resolution.

To implement a skip network, we extract the features of layers Conv-*_2, i.e., produced by the last convolution in each resolution and before max pooling. Additional scoring layers are added to produce classification maps from the intermediate features.

The resulting score maps are then combined as explained in Section 5. We also created a deconvolution network that exactly reflects the base FCN (as in [START_REF] Noh | Learning deconvolution network for semantic segmentation[END_REF]). This is straightforward, with deconvolutional and unpooling layers associated to every convolutional and pooling layer. The only particularity is that the last layer outputs as many maps as required classes and not as input channels. We here call it unpooling network, to differentiate it in the experiments from another method that uses a stack of deconvolutions but without unpooling [START_REF] Volpi | Dense semantic labeling of sub-decimeter resolution images with convolutional neural networks[END_REF], which we simply refer to as deconvolution network. To cover the last family of architectures of Sec. 5.1.2, the dilation network, we incorporate the results recently presented by Sherrah [START_REF] Sherrah | Fully convolutional networks for dense semantic labelling of highresolution aerial imagery[END_REF].

In both datasets we use the same four input channels: DSM, NIR, R and G. Notice Considering the difference in resolution in both datasets, in the case of Potsdam we downsample the input and linearly upsample the output by a factor of 2 (following [START_REF] Sherrah | Fully convolutional networks for dense semantic labelling of highresolution aerial imagery[END_REF]).

We use the same architecture as for Vaihingen (besides the different number of output classes) between the downsampling and upsampling layers. This is to cover a roughly similar receptive field in terms of meters (and not pixels) for both datasets.

Training

The networks are trained by stochastic gradient descent [14]. In every iteration a group of patches is fed to the network for backpropagation. We sample random patches from the images, performing random flips (vertically, horizontally or both) and transpositions, augmenting the data 8 times. At every iteration we group five patches in the mini-batch, of size 256 × 256 for Vaihingen dataset and 512 × 512 for Potsdam (to roughly cover the same geographical area, considering the difference in resolution). In all cases, gradient descent is run with a momentum of 0.9, and an L2 penalty on the network's parameters of 0.0005. Weights are initialized following [START_REF] He | Delving deep into rectifiers: Surpassing human-level performance on imagenet classification[END_REF] and, since we use batch normalization layers before ReLUs, there is no need to normalize the input channels.

We start from a base learning rate of 0.1 and anneal it with an exponential decay.

The decay rate is set so that the learning rate is divided by ten every 10,000 iterations in the case of Vaihingen and every 20,000 iterations in Potsdam. We decrease the learning rate more slowly in the case of Potsdam because the total surface covered by the dataset is larger, thus we assume it must take longer to explore. Training is stopped after 45,000 iterations in the first dataset and 90,000 in the second one, when the error stagnates on the validation set.

To train the unpooling, skip and MLP networks we initialize the weights with the pretrained base FCN, and jointly retrain the entire architecture. We start this second training phase with a learning rate of 0.01, and stop after 30,000 and 65,000 iterations for Vaihingen and Potsdam datasets respectively. We verified that the initialization with the pretrained weights is indeed beneficial compared to training from scratch.

Numerical results

In this section I first present how the base FCN network compares to its derived architectures: unpooling, skip and MLP. I then position MLP with respect to other results reported in the literature, including a dilation network, thus completing the evaluation over all four families of techniques. I finally discuss our submission to the ISPRS contest.

Comparison of a base FCN to its derived unpooling, skip and MLP networks

The classification performances on the validation sets are included in Tables 5.2 and 5.3, for Vaihingen and Potsdam datasets, respectively. The MLP network exhibits the best performance in almost every case. The skip network effectively enhances the results compared with the base network, yet it does not outperform MLP. Let us remark that the unpooling strategy does not necessarily improve the base FCN. This might be a result of the increased training difficulty due to the depth of the network and the sparsity of the unpooled maps. Our attempts to modify the training scheme did not conduce to improving its performance.

Overall, the numerical results show that the incorporation of lower-resolution features significantly improves the classification accuracy. MLP is the most competitive method, boosting the performance by learning how to combine these features.

Comparison with other methods Tables 5.4 and 5.5 (for Vaihingen and Potsdam datasets respectively) incorporate the numerical results reported by other authors using the same training and validation sets. Since not every method was applied to both datasets, the tables do not display exactly the same techniques. The MLP approach also outperforms the dilation strategy, in both datasets, thus positioning it as the most competitive category among those presented in Sections 5.1.2, 5.1.3 (dilation, unpooling, skip, MLP).

In the case of Vaihingen dataset, we also report the results of the deconvolution network [START_REF] Volpi | Dense semantic labeling of sub-decimeter resolution images with convolutional neural networks[END_REF], commented in Sec. 5.1.2, which performs upsampling by using a series of deconvolutional layers. Contrary to the unpooling network, the decoder does not exactly reflect the encoder and no unpooling operations are used. Additionally, we include the performance of other methods recently presented in the literature: the CNN+RF approach [START_REF] Paisitkriangkrai | Effective semantic pixel labelling with convolutional networks and conditional random fields[END_REF], which combines a CNN with a random forest classifier; the CNN+RF+CRF approch, which adds CRF post-processing to CNN+RF; and Dilation+CRF [START_REF] Sherrah | Fully convolutional networks for dense semantic labelling of highresolution aerial imagery[END_REF], which adds CRF post-processing to the dilation network. As depicted in the table, the MLP approach outperforms these other methods too.

For Potsdam dataset, Table 5.5 reports the performance of two other methods, presented in [START_REF] Sherrah | Fully convolutional networks for dense semantic labelling of highresolution aerial imagery[END_REF]. In both cases, a pretrained network based on VGG [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF] is applied to the IR-R-G channels of the image, and another FCN is applied to the DSM, resulting in a huge hybrid architecture. An ordinary version (with upsampling at the end) and a dilation version are considered ('VGG pretr.' and 'VGG+Dilation' in Table 5.5, respectively). In the latter version, the dilation strategy could only be applied partially as it is too memory intensive. While MLP outperforms the non-pretrained simpler dilation network, the VGG+Dilation variants exhibits the best overall performance (though not on all of the individual classes). This suggests that the VGG component might be adding a competitive edge, though the authors stated that this is not the case on the Vaihingen dataset.

Overall, MLP provides better accuracies than most techniques presented in the literature, including dilation networks, ensemble approaches and CRF post-processing.

Submission to the ISPRS challenge

We submitted the result of executing MLP on the Vaihingen test set to the ISPRS server (ID: 'INR'), which can be accessed online [START_REF]ISPRS. 2D semantic labeling contest[END_REF].

Our method scored second out of 29 methods, with an overall accuracy of 89.5%. Note that our MLP technique is very simple compared to other methods in the leaderboard, yet it scored better than them. For example, an ensemble of two skip CNNs was pretrained on large natural image databases [START_REF] Marmanis | Semantic segmentation of aerial images with an ensemble of cnns[END_REF], with over 20 convolutional layers and separate paths for the image and the DSM. Despite being simpler, our MLP network outperforms it in the benchmark.

Visual results

Fig. 5.9 shows visual comparisons on closeups of classified images of both datasets. As expected, the base FCN tends to output "blobby" objects, while the other methods provide sharper results. This is particularly noticeable for the cars of Rows 2, 5 and 6, and for the thin road at the lower left corner of Row 4. We also observe that the incorporation of reasoning at lower resolutions allows the derived networks to discover small objects that are otherwise lost. This is particularly noticeable in the 4th row, where there is a set of small round/cross-shaped objects of the clutter class (in red) that are omitted or grouped together by the base FCN.

The unpooling technique seems to be prone to outputting artifacts. These are often very small in size, even isolated pixels. This is well observed for example for the car of Row 3. This effect could be a natural consequence of the max unpooling mechanism, as depicted in Fig. 5.6, which upsamples into sparse matrices and delegates the task of reconstructing a smoother output to the deconvolutional layers.

At first sight it is more challenging to visually assess why MLP outperforms the skip network in almost every case in the numerical evaluation. Taking a closer look we can however observe that boundaries tend to be more accurate at a fine level in the case of MLP. For example, the "staircase" shape of one of the buildings in Row 1 is noticeably better outlined by the MLP network.

We can also observe that the ground truth itself is often not very precise. For example, the car in Row 3 does not seem to be labeled accurately, hence it is difficult to imagine that a network would learn to finely label that class. In Row 5, an entire lightwell between buildings has apparently been omitted in the ground truth (labeled as part of the building), yet recognized as an impervious surface by the CNNs. The general recognition capabilities of CNNs can also be well appreciated in these fragments. For example, in Row 4, while there are tiny round objects both on the roof of the building and outside the building, CNNs correctly label as building the ones on the roof and as clutter the other ones.

In Fig. 5.10 the classification of entire tiles of the Vaihingen set are shown, obtained from the test set submissions. We include the deconvolution [START_REF] Volpi | Dense semantic labeling of sub-decimeter resolution images with convolutional neural networks[END_REF] and dilation [START_REF] Sherrah | Fully convolutional networks for dense semantic labelling of highresolution aerial imagery[END_REF] network results, together with our MLP. We can see, for instance, that a large white building in the first image is recognized by MLP but misclassified or only partially recovered by the other methods. In the second tile, the Dilation method outputs some holes in the buildings which are not present in the MLP results. A better combination of the information coming from different resolutions might explain why MLP successfully recognizes that these entire surfaces do belong to the same object.

Running times

Table 5.6 reports the running times for training and testing on both datasets. The training time of the architectures derived from the base FCN comprises the time to pretrain the base FCN first and the time to then train the whole system altogether (see details in Sec. 5.1.4). The architectures were implemented using Caffe [START_REF] Jia | Caffe: Convolutional architecture for fast feature embedding[END_REF] and run on an Intel I7 CPU @ 2.7Ghz with a Quadro K3100M GPU (4 GB RAM). We also add for comparison the results reported by the author of the Dilation network [START_REF] Sherrah | Fully convolutional networks for dense semantic labelling of highresolution aerial imagery[END_REF], run on a larger 12 GB RAM GPU. To classify large images we crop them into tiles with as much overlap as the amount of padding in the network, to avoid tile border effects.

As reported in the table, the unpooling, skip and MLP networks introduce an overhead to the base FCN. MLP is the slowest of the derived networks, followed by the unpooling and skip networks. MLP, which provides the highest accuracy, classifies the entire Vaihingen validation set in about 30 seconds and the Postdam validation set in 2 minutes. This is substantially faster than the dilation network. Incorporating the principle of Fig. 5.3 allows us to better allocate computational resources, not spending too much time and space in conducting a high-resolution analysis where it is not needed, boosting accuracy and performance.

Can classification methods generalize to any city?

Over the last few years, there has been a growing interest in processing remote sensing imagery at a global scale, often the entire Earth at once [START_REF] Mattyus | Enhancing road maps by parsing aerial images around the world[END_REF]. New perspectives in remote sensing have particularly highlighted this interest, such as the use of aerial imagery for autonomous driving [START_REF] Mattyus | Enhancing road maps by parsing aerial images around the world[END_REF]. The improvements in the algorithms, and the use of clusters and GPUs have made the processing time less of a constraint. One of the current challenges is to design methods that generalize to different areas of the earth, considering the important intra-class variability encountered over large geographic extents.

The standard way of evaluating and comparing classification methods is to split the labeled data into two sets: one used for training and the other one for testing. For example, in the hyperspectral literature it is particularly common to randomly extract certain pixels from the labeled data and use them for training (ranging from as little as 50 pixels [START_REF] Fauvel | Advances in spectral-spatial classification of hyperspectral images[END_REF] to as much as 20% of all the labeled data [START_REF] Valero | Hyperspectral image representation and processing with binary partition trees[END_REF]), while the rest is used for testing. The Pavia and Indian pines datasets [START_REF] Fauvel | Advances in spectral-spatial classification of hyperspectral images[END_REF] have become the standard benchmarks in the hyperspectral literature. They are mostly geared at distinguishing materials (e.g., bitumen building and bricks), thus leveraging the properties of hyperspectral imagery.

However, those images cover limited geographic areas and the evaluation procedure does not assess how the methods generalize to different contexts or more abstract semantic classes.

With the goal of comparing classification methods over large areas, Mnih [START_REF] Mnih | Machine learning for aerial image labeling[END_REF] created building and road classification datasets over Massachusetts, covering 340 km 2 and 2600 km 2 respectively. For testing, several randomly selected tiles were removed from the reference data. The training set thus covers a geographic surface with "holes", which are used for testing. This situation is analogous to the procedure used for the aforementioned hyperspectral datasets, though taken to a larger scale. While the Massachusetts datasets indeed cover a large surface with significant intra-class variability, the image tiles tend to be self-similar and with uniform color histograms. As shown in [START_REF] Mnih | Machine learning for aerial image labeling[END_REF], a CNN trained on the Massachusetts dataset generalizes poorly to images over Buffalo, and a fine-tuning of the CNN to the new dataset is required.

In the context of high-resolution image classification, the Vaihingen and Potsdam datasets [START_REF] Volpi | Dense semantic labeling of sub-decimeter resolution images with convolutional neural networks[END_REF] have gained increasing attention over the last year, as discussed previously in this chapter. While they provide exhaustive reference data with multiple object classes, the area covered is limited (roughly 1.5 km2 and 3.5 km 2 respectively). The Bavaria and Aerial KITTI datasets [START_REF] Mattyus | Enhancing road maps by parsing aerial images around the world[END_REF], used for road labeling, also cover small surfaces (5 km 2 and 6 km 2 , respectively).

In our experience, and in accordance to [START_REF] Mnih | Machine learning for aerial image labeling[END_REF], training a classifier with images over a particular region and illumination conditions tends to generalize poorly to other images.

For example, Fig. 5.11 depicts a classification map over Zurich into the building/not building classes, created by using the FCN trained over Forez, France. We can observe Lake Zurich being mostly classified as building. Even though there were buildings and body waters in the French imagery, the CNN seems to have learned what a building looks like in that particular images and not simply what a building looks like.

Our goal is to provide a common framework to evaluate classification techniques and, in particular, their generalization capabilities. We created a benchmark database of labeled imagery that covers varied urban landscapes, ranging from highly dense metropolitan financial districts to alpine resorts. The data, referred to as the Inria Aerial Image Labeling Dataset 2 , includes urban settlements over the United States and Austria, and is labeled into building and not building classes. Contrary to all previous datasets, the training and test sets are split by city instead of excluding random pixels or tiles. This way, a system trained, for example, on Chicago, is expected to classify imagery over San

Francisco (with a significantly different appearance). The test set reference data is not publicly released, and a contest has been launched for researchers to submit their results.

In the following, we describe the dataset and then assess the performance of the networks presented previously in this chapter.

The dataset

One of the first key points to decide when creating the dataset was which geographic areas to include and which semantic classes to consider. The criteria were as follows:

• Recent orthorectified imagery available;

• Recent official cadastral records available;

• Precise registration between the cadastral records and the orthorectified imagery;

• Open-access data, both for the images and the cadaster (free to access and distribute);

• Cover varied urban landscapes and illumination.

Let us first highlight the fact that we can only focus on regions where both the images and the reference data are available. In addition, we require the data to be open access in order to freely share our derived dataset with the community. After extensive research, we found that certain US and Austrian areas satisfy those requirements. In the case of the US, public domain orthoimages have been released by USGS through the National Map service (nationalmap.gov) in most urban areas of the country. Vectorial cadastral records have been released through certain local or statewide geographic information system (GIS) websites. We must focus on the zones where such reference data are available in addition to the images.

In the case of Austria, the different provinces have shared images through their respective GIS agencies. We focus, in particular, on Tyrol and Vienna provinces, since open vectorial cadastral data are also on hand. We obtained the images through the WMS services provided by the GIS departments3 as well as the associated reference shapefiles.

The original US imagery is provided at either 15 or 30 cm resolution with three or four spectral bands (RGB/RGB-Infrared), depending on the area, and Vienna imagery contains three bands (RGB) at a resolution of 10 or 20 cm. We took out the common factor and built our dataset with 30 cm images (average resampling if needed) and using the three color bands.

We consider two semantic classes: building and not building. For this we must extract the so-called building footprints from the cadaster. While there are other classes present in some areas (e.g., trees and roads), the building class is the only one that is consistent across different areas. Roads, for example, are often represented with a line, but it is very often not located at the center of the road and its width is usually not specified. This makes it difficult to derive a pixelwise semantic labeling for roads, and is an active research problem itself [START_REF] Mattyus | Enhancing road maps by parsing aerial images around the world[END_REF].

Once we selected a number of candidate areas for the dataset, we visually inspected them to assess whether the cadaster is properly aligned with the images. In some regions, there are irregular shifts that led us to exclude them (e.g., Seattle and Spokane cities).

This may be the result of errors or imprecision in the terrain model used to orthorectify the images, or in the digitization of the cadaster. Note that we have only considered official image and cadaster data sources, ignoring, e.g., OpenStreetMap (OSM) data.

Curiously enough, while we find the official San Francisco building footprints to be perfectly aligned with the USGS imagery, a team of OSM collaborators manually modified 150,000 buildings from these footprints prior to their inclusion in OSM 4 , arguing that they were inaccurate. We now observe them to be misaligned with the imagery. A possible explanation for this is that, at the time of the edit, the Bing images used as the base layer of the OSM editor may have not been geographically precise.

The regions included in the dataset and their distribution into training and test subsets are depicted in Table 5.7. Note first that the amount of data in each of the subsets is the same. This stresses our goal of properly assessing classification methods that generalize to different areas and images. The regions were split in such a way that each of the subsets contains both European and American landscapes, as well as high-density (e.g., Chicago/San Francisco and Vienna/Innsbruck) and low-density (e.g., Kistap/Bloomington, West/East Tyrol) urban settlements. While aerial images over Tyrol are present in both subsets, they have been obtained at different flights over the country, thus exhibiting different illumination characteristics. We have also selected dissimilar images inside some of the groups (e.g., Kitsap County contains tiles from two different flights with very dissimilar characteristics). The reference data was created by rasterizing the shapefiles with GDAL. Fig. 5.12 shows closeups of the images in the dataset.

.12: Close-ups of the dataset images and their corresponding reference data.

We consider two evaluation measures to assess the performance of different methods on the dataset: the accuracy and the intersection over union (IoU) of the positive (building) class. We compute accuracy and IoU on the overall dataset and for every region independently (e.g., San Francisco).

Experiments

We experimented with convolutional neural networks on the dataset. We created a validation set by excluding the first five tiles of each area from the training set (e.g., Austin{1-5}). We first trained the base fully convolutional network (FCN) proposed from Potsdam dataset in Section 5.1.4, for 120,000 iterations on randomly sampled patches of our dataset (momentum is set to 0.9, the L2 penalty to 0.0005 and the learning rate to 0.001). To provide a finer classification, we derived an MLP network on top of the base FCN, as explained in Section 5.1.3 and illustrated in Fig. 5.8. The pretrained FCN was used to initialize the corresponding parameters in the MLP network, and then the overall system was trained for an extra 250,000 iterations, which took 50 hours on a single GPU.

We started with a learning rate of 0.0001, multiplying it by 0.1 every 50k iterations.

The numerical results are summarized in Tables 5.8 and 5.9, for the validation and test sets, respectively. We also include the performance of a skip network as an alternative way of combining features to refine the predictions of the coarse base FCN (see Section 5.1.2). The MLP network reaches about 60% IoU on the entire test set. This means that the output objects overlap the real ones by 60%, as assessed over a significant amount of test data. While there is certainly room for improvement, these values suggest that the current network does generalize well to different cities.

Concluding remarks

In this section, I have presented an overview of different families of dense classification convolutional neural network (CNN) prototypes. We have studied which relevant constraints can be imposed in the architecture by construction, reducing the number of parameters and improving the optimization. We observed that existing networks often spend efforts in learning invariances that could be otherwise guaranteed, and reason at a high resolution even when it is not needed. While previous methods are already competitive, we can devise more optimal approaches. We derived a model in which spatial features are learned at multiple resolutions accuracy with low computational requirements, leading to a win-win situation. Some of the outperformed methods are in fact significantly more complex than our approach, proving once again that striving for simplicity is often the way to go when using CNN architectures.

We also created a dataset for the classification of aerial images. This dataset highlights the need for methods that generalize to the dissimilar appearance of urban settlements around the Earth. Contrary to previous work, the testing is not performed over excluded areas of the training surface, but over entirely different cities instead. We cover a wide range of urban densities, on both European and American cities.

Our experiments with deep neural networks show their satisfactory generalization capability. However, there is still a large room for improvement, as observed in the numerical results. In the following, I present a method, designed with the goal to improve classification results, when there is a significant lack of spatially accurate reference data to train the classification system.

Recurrent neural networks to correct classification maps

In this section, we explore incorporating image information a posteriori in an enhancement module that sharpens the coarse classification maps around objects. The pioneering approach in this direction, called Deeplab [START_REF] Chen | Semantic image segmentation with deep convolutional nets and fully connected CRFs[END_REF], uses a fully connected conditional random field (CRF) to perform the enhancement. Zheng et al. [START_REF] Zheng | Conditional random fields as recurrent neural networks[END_REF] recently reformulated the fully connected CRF of Deeplab as an RNN, and Chen et al. [START_REF] Chen | Semantic image segmentation with task-specific edge detection using CNNs and a discriminatively trained domain transform[END_REF] designed an RNN that emulates the domain transform filter [START_REF] Eduardo | Domain transform for edge-aware image and video processing[END_REF]. Such a filter is used to sharpen the classification maps around image edges, which are themselves detected with a CNN. In these methods the refinement algorithm is designed beforehand and only few parameters that rule the algorithm are learned as part of the network's parameters. The innovating aspect of these approaches is that both steps (coarse classification and enhancement) can be seen as a single end-to-end network and optimized simultaneously.

Instead of predefining the algorithmic details as in previous works, we formulate a general iterative refinement algorithm through an RNN and let the network learn the specific algorithm. To our knowledge, little work has explored the idea of learning an iterative algorithm. In the context of image restoration, the preliminary work by Liu et al. [START_REF] Liu | Learning PDEs for image restoration via optimal control[END_REF][START_REF] Liu | Toward designing intelligent pdes for computer vision: An optimal control approach[END_REF] proposed to optimize the coefficients of a linear combination of predefined terms. Chen et al. [START_REF] Chen | On learning optimized reaction diffusion processes for effective image restoration[END_REF] later modeled this problem as a diffusion process and used an RNN to learn the linear filters involved as well as the coefficients of a parametrized nonlinear function. Our problem is however different, in that we use the image as guidance to update a classification map, and not to restore the image itself. Besides, while we drew inspiration on diffusion processes, we are also interested in imitating other iterative processes like active contours, thus we do not restrict our system to diffusions but consider all PDEs.

Learning iterative processes to enhance classification

Let us assume we are given a set of score (or "heat") maps u k , one for each possible class k ∈ Ω, in a pixelwise labeling problem. The score of a pixel reflects the likelihood of belonging to a class, according to the classifier's predictions. The final class assigned to every pixel is the one with maximal value u k . Alternatively, a softmax function can be used to interpret the results as probability scores: P (k) = e u k ∑ j∈Ω e u j , as usually done in CNN literature. Fig. 5.14 shows a sample of the type of fuzzy heat map outputted by a CNN for the class 'building'.

Our goal is to combine the score maps u k with information derived from the input image (e.g., edge features) to sharpen the scores near the real objects in order to enhance the classification. One way to perform such a task is to progressively enhance the score maps by using partial differential equations (PDEs). We first describe different types of PDEs we could certainly imagine to design in order to solve our problem. Instead of discussing which one is the best, we later propose a generic iterative process to enhance the classification maps without specific constraints on the algorithm rationale. Finally, we show how this equation can be expressed and trained as a recurrent neural network (RNN).

Partial differential equations (PDEs)

We can formulate a variety of diffusion processes applied to the maps u k as partial differential equations. For example, the heat flow is described as:

∂u k (x) ∂t = div(∇u k (x)), (5.5)
where div(⋅) denotes the divergence operator in the spatial domain of x. Applying such a diffusion process in our context would smooth out the heat maps. Instead, our goal is to design an image-dependent smoothing process that aligns the heat maps to the image features. A natural way of achieving this is to modulate the gradient in Eq. 5.5 by a scalar function g(x, I) that depends on the input image I:

∂u k (x) ∂t = div(g(I, x)∇u k (x)). (5.6)
Eq. 5.6 is similar to the Perona-Malik diffusion [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] with the exception that Perona-Malik uses the smoothed function itself to guide the diffusion. g(I, x) denotes an edge-stopping function that takes low values near borders of I(x) in order to slow down the smoothing process there.

Another possibility would be to consider a more general variant in which g(I, x) is replaced by a matrix D(I, x), acting as a diffusion tensor that redirects the flow based on image properties instead of just slowing it down near edges:

∂u k (x) ∂t = div(D(I, x)∇u k (x)). (5.7)
This formulation is related to the so-called anisotropic diffusion process [START_REF] Weickert | Anisotropic diffusion in image processing[END_REF].

Alternatively, one can draw inspiration from the level set framework. For example, the geodesic active contours technique formulated as level sets translates into:

∂u k (x) ∂t = ∇u k (x) div g(I, x) ∇u k (x) ∇u k (x) . (5.8)
Such a formulation favors the zero level set to align with minima of g(I, x) [START_REF] Caselles | Geodesic active contours[END_REF]. Schemes based on Eq. 5.8 could then be used to improve heat maps u k , provided they are scaled so that segmentation boundaries match zero levels.

As shown above, many different PDE approaches can be devised to enhance classification maps. However, several choices must be made to select the appropriate PDE and tailor it to our problem. For example, one must choose the edge-stopping function g(I, x) in Eqs. 5.6, 5.8. Common choices are exponential or rational functions on the image gradient [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF], which in turn requires to set an edge-sensitivity parameter. Extensions to the original Perona-Malik approach could also be considered, such as a popular regularized variant that computes the gradient on a Gaussian-smoothed version of the input image [START_REF] Weickert | Anisotropic diffusion in image processing[END_REF]. In the case of opting for anisotropic diffusion, one must design D(I, x).

Instead of using trial and error to perform such design, our goal is to let a machine learning system discover by itself a useful iterative process for our task.

Generic classification enhancement

PDEs are usually discretized in space by using finite differences, which represent derivatives as discrete convolution filters. We build upon this scheme to write a generic discrete formulation of an iterative enhancement process.

Let us consider that we take as input a score map u k (for class k) and, in the most general case, an arbitrary number of feature maps {g 1 , ..., g p } derived from image I. In order to perform differential operations, of the type { ∂ ∂x , ∂ ∂y , ∂ 2 ∂x∂y , ∂ 2 ∂x 2 , ...}, we consider convolution kernels {M 1 , M 2 , ...} and {N j 1 , N j 2 , ...} to be applied to the heat map u k and to the features g j derived from image I, respectively. While we could certainly directly provide a bank of filters M i and N j i in the form of Sobel operators, Laplacian operators, etc., we may simply let the system learn the required filters. We group all the feature maps that result from applying these convolutions, in a single set:

Φ(u k , I) = M i * u k , N j l * g j (I) ; ∀i, j, l . (5.9)
Let us now define a generic discretized scheme as:

∂u k (x) ∂t = f k Φ(u k , I)(x) , (5.10)
where f k is a function that takes as input the values of all the features in Φ(u k , I) at an image point x, and combines them. While convolutions M i and N j i convey the "spatial" reasoning, e.g., gradients, f k captures the combination of these elements, such as the products in Eqs. 5.6 and 5.8.

Instead of deriving an arbitrary number of possibly complex features N j i * g j (I) from image I, we can think of a simplified scheme in which we directly operate on I, by considering only convolutions: N i * I. The list of functionals considered in Eq. 5.10 is then

Φ(u k , I) = M i * u k , N j * I ; ∀i, j (5.11)
and consists only of convolutional kernels directly applied to the heat maps u k and to the image I. From now on, we here stick to this simpler formulation, yet we acknowledge that it may be eventually useful to work on a higher-level representation rather than on the input image itself. Note that if one restricts functions f k in Eq. 5.10 to be linear, we still obtain the set of all linear PDEs. We consider any function f k , thus introducing non-linearities.

PDEs are usually discretized in time, taking the form:

u k,t+1 (x) = u k,t (x) + δu k,t (x), (5.12)
where δu k,t denotes the overall update of u k,t at time t.

Note that the convolution filters in Eqs. 5.9 and 5.11 are class-agnostic: M i , N j and N j l do not depend on k, while f k may be a different function for each class k. Function f k thus determines the contribution of each feature to the equation, contemplating the case in which a different evolution might be optimal for each of the classes, even if just in terms of a time-step factor. In the following, we detail a way to learn the update functions δu k,t from training data.

Iterative processes as RNNs

We now show that the generic iterative process can be implemented as an RNN, and thus trained from labeled data. This stage requires to provide the system with a piece of accurately labeled ground truth (see e.g., Fig. 5

.15).

Let us first show that one iteration, as defined in Eqs. 5.10-5.12, can be expressed in terms of common neural network layers. Let us focus on a single pixel for a specific class, simplifying the notation from u k,t (x) to u t . Fig. 5.16 illustrates the proposed network architecture. Each iteration takes as input the image I and a given heat map u t to enhance at time t. In the first iteration, u t is the initial coarse heat map to be improved, outputted by another pre-trained neural network in our case. From the heat map u t we derive a series of filter responses, which correspond to M i * u t in Eq. 5.11. These These filters are then "concatenated", forming the pool of features Φ coming from both the input image and the heat map, as in Eq. 5.11, and inputted to f k in Eq. 5.10.

We must now learn the function δu t that describes how the heat map u t is updated at iteration t (see Eq. 5.12), based on these features.

Eq. 5.10 does not introduce specifics about function f k . In (5.5)-(5.8), for example, it includes products between different terms, but we could certainly imagine other functions. We therefore model δu t through a multilayer perceptron (MLP), because it can approximate any function within a bounded error [14]. We include one hidden layer with nonlinear activation functions followed by an output neuron with a linear activation (a typical configuration for regression problems), although other MLP architectures could be used. Applying this MLP densely is equivalent to performing convolutions with 1 × 1 kernels at every layer. The implementation to densely label entire images is then straightforward.

The value of δu t is then added to u t in order to generate the updated map u t+1 (see Fig. 5.16). This addition is performed pixel by pixel in the case of a dense input. Note that although we could have removed this addition and let the MLP directly output the updated map u t+1 , we opted for this architecture since it is more closely related to the equations and better conveys the intention of a progressive refinement of the classification map. Moreover, learning δu t instead of u t+1 has a significant advantage at The overall iterative process is implemented by unrolling a finite number of iterations, as illustrated in Fig. 5.17, under the constraint that the parameters are shared among all iterations. Such a sharing is enforced at training time by a simple modification to the back-propagation training algorithm where the derivatives of every instance of a weight at different iterations are averaged [START_REF] Werbos | Backpropagation through time: what it does and how to do it[END_REF]. Note that the spatial features are shared across the classes, while a different MLP is learned for each of them, following Eq. 5.10. As depicted by Fig. 5.17 and conveyed in the equations, the features extracted from the input image are independent of the iteration.

The RNN of Fig. 5.17 represents a dynamical system that iteratively improves the class heat maps. Training such an RNN amounts to finding the optimal dynamical system for our enhancement task.

Experimental results

Implementation details

We first describe the CNN used to produce the coarse predictions, then detail the proposed RNN. The employed CNN is based on a previous remote sensing network presented by Mnih [START_REF] Mnih | Machine learning for aerial image labeling[END_REF]. We have designed a fully convolutional [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF] version of Mnih's network, since recent remote sensing work has shown the theoretical and practical advantages of this type of architecture [START_REF] Kampffmeyer | Semantic segmentation of small objects and modeling of uncertainty in urban remote sensing images using deep convolutional neural networks[END_REF][START_REF] Maggiori | Fully convolutional neural networks for remote sensing image classification[END_REF]. The CNN takes 3-band color image patches at 1m 2 resolution and produces as many heat maps as classes considered. The resulting fourlayer architecture is as follows: 64 conv. filters (12 × 12, stride 4) → 128 conv. filters (3 × 3) → 128 conv. filters (3 × 3) → 3 conv. filters (9 × 9). Since the first convolution is performed with a stride of 4, the resulting feature maps have a quarter of the input resolution. Therefore, a deconvolutional layer [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF] is added on top to upsample the classification maps to the original resolution. The activation functions used in the hidden layers are rectified linear units. This network is trained on patches randomly selected from the training dataset. We group 64 patches with classification maps of size 64 × 64 into mini-batches (following [START_REF] Mnih | Machine learning for aerial image labeling[END_REF]) to estimate the gradient of the network's parameters and back-propagate them. Our loss function is the cross-entropy between the target and predicted class probabilities. Stochastic gradient descent is used for optimization, with learning rate 0.01, momentum 0.9 and an L2 weight regularization of 0.0002. We did not however optimized these parameters nor the networks' architectures.

We now detail the implementation of the RNN described in Sec. 5.2.1. Let us remark that at this stage we fix the weights of the initial coarse CNN and the manually labeled tile is used to train the RNN only. Our RNN learns 32 M i and 32 N j filters, both of spatial dimensions 5×5. An independent MLP is learned for every class, using 32 hidden neurons each and with rectified linear activations, while M i and N j filters are shared across the different classes (in accordance to Equations 5.10 and 5.11). This highlights the fact that M i and N j capture low-level features while the MLPs convey class-specific behavior. We unroll five RNN iterations, which enables us to significantly improve the classification maps without exhausting our GPU's memory. Training is performed on random patches and with the cross-entropy loss function, as done with the coarse CNN. The employed gradient descent algorithm is AdaGrad [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF], which exhibits a faster convergence in our case, using a base learning rate of 0.01 (higher values make the loss diverge). All weights are initialized randomly by sampling from a distribution that depends on the number of neuron inputs [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF]. We trained the RNN for 50,000 iterations, until observing convergence of the training loss, which took around four hours on a single GPU.

Dataset and results

We perform our experiments on images acquired by a Pléiades satellite over the area of Forez, France. An RGB color image is used, obtained by pansharpening [START_REF] Wang | A comparative analysis of image fusion methods[END_REF] the satellite data, which provides a spatial resolution of 0.5 m 2 . Since the networks described hereabove are designed for 1 m 2 resolution images, we downsample the Pléiades images before feeding them to our networks and bilinearly upsample the outputs.

From this image we selected an area with OpenStreetMap (OSM) [START_REF] Haklay | Openstreetmap: User-generated street maps[END_REF] misaligned or omitted, while many roads in the ground truth are not visible in the image (or the other way around). Moreover, since OSM's roads are represented by polylines, we set a fixed road width of 7 m to rasterize this class (following [START_REF] Mnih | Machine learning for aerial image labeling[END_REF]), which makes their classification particularly challenging. This dataset is used to train the initial coarse CNNs.

We manually labeled two 2.25 km 2 tiles to train and test the RNN at enhancing the predictions of the coarse network. We denote them by enhancement and test sets, respectively. Note that our RNN system must discover an algorithm to refine an existing classification map, and not to conduct the classification itself, hence a smaller training set should be sufficient for this stage. The enhancement set is depicted in Fig. 5.18 while the test set is shown in Figs. [START_REF] Becker | Supervised feature learning for curvilinear structure segmentation[END_REF].22(a)/(f).

In the following, we report the results obtained by using the proposed method on the Pléiades dataset. The initial classification of roads has an overlap of less than 10% with the roads in the ground truth, as shown by its individual IoU. The RNN makes them emerge from the background class, now overlapping the ground truth roads by over 50%. Buildings also become better aligned to the real boundaries, going from less than 40% to over 70% overlap with the ground truth buildings. This constitutes a multiplication of the IoU by a factor of 5 for roads and 2 for buildings, which indicates a significant improvement at outlining and not just detecting objects.

Additional visual fragments before and after the RNN refinement are shown in Fig. 5.21.

We can observe in the last row how the iterative process learned by the RNN both thickens and narrows the roads depending on the location. To validate the importance of using a recurrent architecture, and following Zheng et al. [START_REF] Zheng | Conditional random fields as recurrent neural networks[END_REF], we retrained our system considering every iteration of the RNN as an independent step with its own parameters. After training for the same number of iterations, it yields a lower performance on the test set compared to the RNN and a higher performance on the training set. If we keep on training, the non-recurrent network still enhances its training accuracy while performing poorly on the test set, implying a significant degree of overfitting with this variant of the architecture. This provides evidence that constraining our network to learn an iterative enhancement process is crucial for its success.

Feature visualization

Though it is difficult to interpret the overall function learned by the RNN, especially the part of the multi-layer perceptron, there are some things we expect to find if we analyze the spatial filters M i and N j learned by the RNN (see Eq. 5.11). Carrying out this analysis is a way of validating the behavior of the network.

The iterative process learned by the RNN should combine information from both the heat maps and the image at every iteration (since the heat maps constitute the prior on where the objects are located, and the image guides the enhancement of these heat maps). A logical way of enhancing the classification is to align the high-gradient areas of the heat maps with the object boundaries. We expect then to find derivative operators among the filters N j applied to the heat maps. Concerning the image filters N j , we expect to find data-dependent filters (e.g., image edge detectors) that help identify the location of object boundaries.

To interpret the meaning of the filters learned by the RNN we plot the map of responses of a sample input to the different filters. Color image Feature responses

Concluding remarks

I have described a recurrent neural network (RNN) that learns how to refine the coarse output of another neural network, in the context of pixelwise image labeling. The inputs are both the coarse classification maps to be corrected and the original color image. The output at every RNN iteration is an update to the classification map of the previous iteration, using the color image for guidance.

Little human intervention is required, since the specifics of the refinement algorithm are not provided by the user but learned by the network itself. For this, we analyzed different iterative alternatives and devised a general formulation that can be interpreted as a stack of common neuron layers. At training time, the RNN discovers the relevant features to be taken both from the classification map and from the input image, as well as the function that combines them.

The experiments on satellite imagery show that the classification maps are improved significantly, increasing the overlap of the foreground classes with the ground truth.

To conclude, we demonstrated that RNNs succeed in learning iterative processes for classification enhancement tasks. I can conclude that for specific applications, it is often advantageous to incorporate and enforce the known shape priors (e.g., as we have seen in Chapter 4, if we know that a shape of the object will only grow in time series, enforcing this shape constraint within the classification system naturally boosts its performance). However, if we aim at designing methods to automatically classify remote sensing images on a world-scale, the developed learning approaches must be generic and highly scalable. Within this context, deep learning techniques, in particular convolutional neural networks have gained significant attention, since they have demonstrated their ability to learn expressive multi-scale contextual features and have shown a remarkable computational performance. We have shown in Chapter 5 that the modern learning methods succeed in classifying new areas on the Earth belonging to the cities the classifier has never seen before. However, there is still significant work to do for designing systems which would be able to automatically create and update maps (e.g., geographic information system (GIS) maps) on a world-scale, using remote sensing images.

Future Work

To explore these ideas further, I have proposed the project EPITOME (Efficient reP-resentatIon TO structure large-scale satellite iMagEs), which was funded for the period of 2017-2021 by the French National Research Agency ANR. The goal of the EPITOME project is to devise novel effective representations for large-scale satellite images, which would represent and structure the essential content (which we call epitome) of input images, thus enabling novel possibilities for the exploitation of remote sensing datasets, including novel efficient navigation modes through space, scale, and perhaps time, or more powerful indexation mechanisms. Such representations must be highly generic and suited to massive image data, to be applicable for images from all over the world and for a wide range of applications. At the same time, they should be structure-preserving, i.e. best represent the meaningful objects in the image scenes, with the possibility to be enriched by semantic or other kinds of information.

The vector-based representation is well-known to provide several important advantages over raster images, the three most salient ones being compactness, scalability and easiness of updating [START_REF] Liao | A subdivision-based representation for vector image editing[END_REF]. These are the main reasons the GIS maps are stored as the database of vector primitives. Thus, I opt to devise a multi-resolution vector-based representation, together with the methods for its efficient generation, manipulation and continuous navigation between scales. Each level (scale) of this representation would convey different level of details from the image. Here I describe a non-exhaustive list of milestones I have identified to address the presented goal:

1. I plan to investigate the use of a large source of free-access online map data, such as OpenStreetMap [START_REF] Haklay | Openstreetmap: User-generated street maps[END_REF], to automatically derive weakly labeled training data for learning about both geometric and semantic structures, their relations in the image scene and their accurate scale alignment. In particular, my current work deals with the development of a data alignment framework, by exploring CNN architectures to learn how to align maps with satellite images.

2. I envision to study how to adapt and extend the state-of-the-art image vectorization methods, for instance the algorithm based on iterative simplification of triangular mesh [START_REF] Xia | Patch-based image vectorization with automatic curvilinear feature alignment[END_REF], to my task of representing satellite data. This implies enforcing the preservation of geometric structures and semantic labels discovered by applying a machine learning algorithm. In our initial work, we have proposed a novel polygonization algorithm, which vectorizes the input classification maps and is based on the approximation of a triangular mesh to the considered map [START_REF] Maggiori | Polygonization of remote sensing classification maps by mesh approximation[END_REF].

3. It would certainly be interesting to design neural networks that output directly geometric primitives instead of raster classification maps. This objective is inspired by a few recent works on using neural networks to output geometric objects. For example, a recurrent neural network has been used to generate the convex hull of a point cloud [START_REF] Vinyals | Pointer networks[END_REF] and in [START_REF] Erhan | Scalable object detection using deep neural networks[END_REF] the coordinates of object bounding boxes are regressed by a CNN. I envision to further bridge the gap between deep machine learning and geometric modeling tools, by devising novel methods which allow learning directly in a space of vector primitives. A possible approach consists in defining a discrete finite set of vector primitives within admitted discretization error bounds, and then designing an appropriate neural network architecture capable to learn how to output the desired representations in this space of primitives.

R m . The cost c ψ 1 associated with this move will then be

c ψ 1 = Σ i j=2 c τ j + c τ i + c R m , (A.8)
considering the location of the new cut. Analogously, the cost when pasting R m k units up of i is:

c ρ 1 = Σ i-k j=2 c τ j + c τ i-k + c R m . (A.9)
Notice in the previous expression that we consider that the cut is still as low as R m . The cut could not be higher, because if it were the case, then it would have already been cut there before.

Let us now see if the cost c ρ 1 could increase as k advances:

Σ i-k j=2 c τ j + c τ i-k + c R m -(Σ i-k+1 j=2 c τ j + c τ i-k+1 + c R m ) = -c τ i-k+1 + c τ i-k -c τ i-k+1 > 0 ⇔c τ i-k > c τ i-k+1 + c τ i-k+1 , (A.10)
contradicting the algorithm to compute the cuts, hence the cost at R 1 must monotonously decrease.

Proposition 3. Let us suppose we paste R m at or over the initial cut of tree τ , leading to tree ϕ. Let us consider we paste higher instead, producing tree ψ. It must then be c ψ 1 ⩾ c ϕ 1 . Proof. Let us resume the proof of Proposition 2. It was shown that c ψ 1 < c ϕ 1 if and only if III was false. At that point we could not contradict III but show that under certain conditions it would be contradicted. Now we will show that the fact that we know the first cut was at or below τ i will contradict III.

After III and γ we have:

e τ i-1 + c ϕ 1 -Σ i j=2 c τ j -c τ i ⩽ e τ i-1 + c R m < c ϕ 1 -Σ i-1 j=2 c τ j ⇔e i-1 τ -c τ i -c τ i ⩽ e τ i-1 + c R m < 0.
(A.11)

If we now add the knowledge about the cut being below τ i-1 , it must be e τ i-1 > c τ i + c τ i , then

c τ i + c τ i -c τ i -c τ i < e τ i-1 -c i τ -c τ i ⩽ e τ i-1 + c R m < 0 ⇔ 0 < 0. (A.12)
Therefore, III cannot be true, which proves the proposition.

instead and reverse names i and j. Three cases are possible:

1. τ i ⩽ τ ′ i ⩽ τ j ⩽ τ ′ j 2. τ i ⩽ τ j ⩽ τ ′ i ⩽ τ ′ j 3. τ i ⩽ τ j ⩽ τ ′ j ⩽ τ ′ i
In case (1), the terms of the submodularity condition write:

-t τ i τ ′ i τ j τ ′ j W i,j (τ i , τ j ) ∶ AAABBBBBBBBDDDDDDD W i,j (τ ′ i , τ ′ j ) ∶ AAAAAAABBBBB B BDDDD W i,j (τ ′ i , τ j ) ∶ AAAAAAABBBBDDDDDDD W i,j (τ i , τ ′ j ) ∶ AAABBBBBBBBB B BDDDD
and the submodularity condition is satisfied if the sum of the two first rows is less than or equal to the sum of the two last rows. It turns out that in the particular case above, the sums are equal for each instant t, and consequently the submodularity condition is satisfied as being an equality.

In case (2), the terms of the submodularity condition write:

-t τ i τ j τ ′ i τ ′ j W i,j (τ i , τ j ) ∶ AAABBBBDDDDDDDDDDD W i,j (τ ′ i , τ ′ j ) ∶ AAAAAAAAAAAB B BDDDD W i,j (τ ′ i , τ j ) ∶ AAAAAAACCCCDDDDDDD W i,j (τ i , τ ′ j ) ∶ AAABBBB B B B B B B BDDDD This time, the two sums, at any instant t, are equal if t < τ j or t ⩾ τ ′ i . For other values of t, that is to say when τ j ⩽ t < τ ′ i , the sum A + D is to be compared with B + C, i.e. the question is to compare W t i,j (0, 0) + W t i,j (1, 1) with W t i,j (0, 1) + W t i,j (1, 0). As E t is binary submodular by hypothesis, and by definition of binary submodularity in equation (4.2), we have A + D ⩽ B + C, and thus the sum of the two first rows is less than or equal to the sum of the two last ones, which means precisely W i,j (τ i , τ j ) + W i,j (τ ′ i , τ ′ j ) ⩽ W i,j (τ i , τ ′ j ) + W i,j (τ ′ i , τ j ), and thus the submodularity condition for E is checked.

Finally, in case (3), the terms of the submodularity condition write:

-t τ i τ j τ ′ j τ ′ i W i,j (τ i , τ j ) ∶ AAABBBBDDDDDDDDDDD W i,j (τ ′ i , τ ′ j ) ∶ AAAAAAAAAAACCCDDDD W i,j (τ ′ i , τ j ) ∶ AAAAAAACCCCCCCDDDD W i,j (τ i , τ ′ j ) ∶ AAABBBB B B B BDDDDDDD
The two sums are equal for all instants t < τ j or t ⩾ τ ′ j . When τ j ⩽ t < τ ′ j , the sum A + D is to be compared with B + C, as previously, and consequently the submodularity condition for E is checked again. In all cases (1), ( 2) and ( 3), E is proven to be submodular, and this concludes the proof.

Output◻

  One of the most important problems for the automatic interpretation of remote sens-Input ⇒ Impervious surf. ∎ Building ∎ Low veget. ∎ Tree ∎ Car ∎ Clutter

Figure 1 . 1 :

 11 Figure 1.1: Example of pixelwise classification of an aerial image.

  Fig. 1.1 shows an example of such pixelwise classification of an aerial image into six classes of interest.

Figure 1 . 2 :

 12 Figure 1.2: Examples of aerial images around the Earth, where large intra-class variability can be appreciated for 'building' class.

  4. H. Aghighi, J. Trinder, S. Lim and Y. Tarabalka, "Fully spatially adaptive smoothing parameter estimation for Markov random field super-resolution mapping of remotely sensed images," International Journal of Remote Sensing, vol. 36, no. 11, pp. 2851-2879, 2015. 5. B. B. Damodaran, R. R. Nidamanuri and Y. Tarabalka, "Dynamic ensemble selection approach for hyperspectral image classification with joint spectral and spatial International conference papers 1. E. Maggiori, Y. Tarabalka, G. Charpiat and P. Alliez. "Polygonization of Remote Sensing Classification Maps by Mesh Approximation", IEEE ICIP'17. Beijing, China, 2017. 2. E. Maggiori, Y. Tarabalka, G. Charpiat and P. Alliez, "Can semantic labeling methods generalize to any city? The Inria aerial image labeling benchmark," IGARSS'17, Fort Worth, Texas, USA, 2017. 3. E. Maggiori, Y. Tarabalka, G. Charpiat and P. Alliez, "High-resolution image classification with convolutional networks," IGARSS'17, Fort Worth, Texas, USA, 2017. 4. E. Maggiori, Y. Tarabalka, G. Charpiat and P. Alliez, "Fully convolutional neural networks for remote sensing image classification," IGARSS'16, Beijing, China, 2016. 5. E. Maggiori, Y. Tarabalka and G. Charpiat, "Optimizing partition trees for multiobject segmentation with shape prior," BMVC'2015, Swansea, UK, 2015. 6. S-G. Jeong, Y. Tarabalka and J. Zerubia, "Stochastic model for curvilinear structure reconstruction using morphological profiles," IEEE ICIP'15, Quebec City, Canada, 2015. 7. E. Maggiori, Y. Tarabalka and G. Charpiat, "Improved partition trees for multiclass segmentation of remote sensing images," IGARSS'15, Milan, Italy, 2015. 8. S-G. Jeong, Y. Tarabalka and J. Zerubia, "Marked point process model for curvilinear structures extraction," EMMCVPR'15, Hong Kong, 2015. 9. S-G. Jeong, Y. Tarabalka and J. Zerubia, "Marked point process model for facial wrinkle detection," IEEE ICIP'14, Paris, France, 2014. 10. Y. Tarabalka and A. Rana, "Graph-cut-based model for spectral-spatial classification of hyperspectral images," IGARSS'14, Quebec City, Canada, 2014. 11. H. Aghighi, J. Trinder, K. Wang, Y. Tarabalka and S. Lim, "Smoothing parameter estimation for Markov random field classification of non-Gaussian distribution image," ISPRS TC VII Symposium, Istanbul, Turkey, 2014. 12. Y. Tarabalka, G. Charpiat, L. Brucker and B. H. Menze, "Enforcing monotonous shape growth or shrinkage in video segmentation," BMVC'2013, Bristol, UK, 2013.Chapter 2

Figure 2 . 1 :

 21 Figure 2.1: Examples of the line configurations with different prior energies: (a)-(c) show preferable line configurations composed of aligned lines (a), almost perpendicular lines (b), and adjacent lines (c). (d)-(g) depict unfavourable line configurations which are penalized because of a singular segment (d), acute corner (e), overlap (f), and parallel (g), respectively.

Figure 2 . 2 :

 22 Figure 2.2: Given configuration (a), if a line segment s ′ proposed by LT kernel is rejected (b), the delayed rejection kernel searches for the nearest extremes in the rest of line segments (c). An alternative line segment s, which enforces connectivity, will be proposed by interpolation of the retrieved points (d). The BD kernel changes the dimensionality of the current configuration s by adding a new line segment or removing an existing line segment. When the birth kernel proposes a new configuration s ′ = s ∪ {s}, the length and the orientation of the new line segment are uniformly sampled from the mark space M = [ min , max ] × [θ min , θ max ], where min and max are the minimum and maximum length of the line segment, respectively. θ min and θ max denote the minimum and maximum orientation of the line segment, respectively.

1 .

 1 For a practical reason related to the implementation, we bound the values of w. Specifically, we sweep the weighting coefficients of the prior energy w p according to the gradient magnitude and noise level of the input image. Let χ =min ×E[∇I]+Var[I σ 2 ] be a baseline to accept a new line segment into the current configuration without considering spatial interaction, where I σ 2 denotes a smoothed image using a Gaussian kernel with σ 2 = {1.5, 2.25, 3.5}. To reduce computation overhead, we fix the weighting factors of data likelihood energy as ω m d = -1 and ω v d = We set w 1 = [-1, 1, χ , 0.1 χ , 0.01 χ , 0.01 χ ] ⊺ , and gradually change χ by 10% of increments, i.e., w 2 = [-1, 1, χ 2 , 0.1 χ 2 , 0.01 χ 2 , 0.01 χ 2 ] ⊺ , where χ 2 = 1.1 χ . In our experiments, we set K = 15 to create line hypotheses.

Figure 2 . 3 :Fig. 2 .

 232 Figure 2.3: Given the input image (a), we compute the gradient magnitude (b).Mathematical morphology operator, path opening[START_REF] Talbot | Efficient complete and incomplete path openings and closings[END_REF], is applied on such gradient magnitude image (c). Linearity score of each pixel is drawn by the supervised feature learning algorithm[6] (d). We provide line hypotheses (e)-(g) associated with different hyperparameter vectors. Composition result (h) is equivalent to mixture probability density, and it highlights pixels corresponding to linear structures.

  )

  (a) T 0 = 100 (b) T 10 2 = 21.67 (c) T 10 3 = 14.47 (d) T 10 4 = 10.86

Figure 2 . 4 :

 24 Figure 2.4: We provide intermediate sampling processes when the temperature parameter T t is decreasing. The results shown in first row are obtained without specifying seed segment. For the second row, we randomly set 20 seed segments and run the algorithm. For the third row, we initialize 20 line segments which are highly corresponding to underlying curvilinear structures. The algorithm converges toward almost the same solution regardless of the initial state.

IFigure 2 . 6 :

 26 Figure 2.6: We visualize the localization of the curvilinear structures on input images (a). We compare with the results of a manually labeled image by a human expert (b), morphological filtering [177] (c), supervised feature learning [6] (d), baseline MPP (e), and the proposed algorithm (f). Threshold values of (c) and (d) are chosen to achieve the closest recall scores to the proposed method. We use blue pixels to indicate areas which are completely corresponding to (b). Green and red pixels denote over-detected and under-detected areas, respectively, as compared with ground-truth. The name of the test images is from top to bottom: crack, wrinkle, dna, and retina.

  45 ○ θ = 0 ○ θ = 45 ○ θ = 90 ○ χ 2

Figure 2 . 8 :

 28 Figure 2.8: In this figure, we compare the statistics of pixel values in the rotated image patches with respect to the given binary mask B ∈ {0, 1} M . Assume that B encodes a basis spatial pattern of the curvilinear structure toward the baseline orientation ( θ = 0 ○ ), where darker gray color refers to 1. To estimate local orientation of the input image patch, we rotate it and find the orientation which maximizes the statistical difference of pair distributions. If the input image patch does not contain curvilinear structure (upper row), there is no meaningful statistical difference for any orientation.

Figure 2 . 9 :

 29 Figure 2.9: For the given image (a), we compute (b) the structured output ranking with Structured SVM. We retain binary map (c) from the highest to ρ I -th rankings. However, it could be broken. We exploit the basis spatial pattern B to reconstruct (d) the structured score map. Graph-based representation provides a tool to analyze the topological features of the curvilinear structures. (e) The coarsest curvilinear structure evolves to define (f) branches.

Figure 2 . 10 :

 210 Figure 2.10: Toy example of the proposed curvilinear structure inference algorithm: (a) input image contains a curvilinear structure which is denoted by gray color; (b) subgraph G ′ is induced from the structured segmentation map; (c) and (d) show the intermediate processes of the 2-sweep algorithm starting from vertex t to find a diameter of the subgraph; (e) we assign 0 weight for all edges on the path; and (f) we repeat the process and add branches if the path length is larger than pre-defined length ˆ .

Figure 2 . 11 :

 211 Figure 2.11: Intermediate steps of the curvilinear structure reconstruction for a retina image. We iteratively reconstruct the curvilinear structure according to topological importance orders. As the iteration goes on, detail structures (layer) appear.

Figure 2 . 12 :

 212 Figure 2.12: Illustration of the curvilinear structure reconstruction results on Road, Crack, DRIVE, and RecA datasets (top to bottom). We compare the results of Frangi et al.'s [52], Law & Chung's [110], Becker et al.'s [5], Sironi et al.'s [169], and the proposed algorithm (left to right).

Figure 3 . 2 :

 32 Figure 3.2: A BPT is constructed by iteratively removing edges in a region adjacency graph (RAG). The resulting BPT encodes the history of the merges.

Figure 3 . 3 :

 33 Figure 3.3: Faulty BPT: the object (bcde) is not represented in a single node, since a part of it (b) merged first to something else.

Fig. 3 .

 3 Fig.3.4(a) illustrates such a move: α is the paste place and β is the pruning place. We denote by LCA(α, β) their lowest common ancestor in the tree. The move creates a new node αβ in the paste side that comprises α and β. In the pruned side, the tree is collapsed after β is removed. In a balanced tree, which is encouraged by setting α > 0 in (3.7), the number of possible moves is bounded by O(n log(n)) (see proof in Appendix A). The neighborhood system is much richer than, for instance, Markov random fields (MRFs) on the pixel grid, considering that it comprises pairs of adjacent regions at several scales. We store at each node R the branch cost C(R) of the best possible cut within its branch. When applying a move as depicted in Fig.3.4(a), it is necessary to recompute the branch cost C till the ancestry of α and β only. The rest of the branches are unaffected,

Figure 3 . 4 :

 34 Figure 3.4: Family of moves to optimize BPTs.

Figure 3 .Figure 3 . 5 : 5 Figure 3 . 6 :

 335536 Figure 3.5-b). This reference was constructed by combining the labeling of isolated pixels provided with the original image, visual inspection and publicly available official Italian records of building boundaries. Since the boundaries of buildings are well defined,

Figures 3 .

 3 Figures 3.5(e-f) and the close-ups of Figure 3.6 compare the results obtained by applying the unsupervised and supervised approaches for BPT construction. These figures isolate the tile objects from the rest and assign a random color to every individual object.

  (a) Input (400 × 267) (b) Gorelick et al. (c) BPT Opt

Figure 3 . 7 :

 37 Figure 3.7: Convex object. Method [69] and optimized BPTs.

  (a) Input (457 × 454) (b) Gorelick et al. [69] (c) Gorelick et al. [69], on usersupplied windows (d) BPT (e) Optimized BPT Input Gorelick et al. (ω = 0.01) Gorelick et al. (ω = 0.001) BPT opt.(f) Amplified fragment

Figure 3 . 8 :

 38 Figure 3.8: Multiple convex objects (cell nuclei). Gorelick et al. [69] and BPT optimization.

Figs. 3

 3 .9(c-f) illustrate the output of these techniques for the image of New York. Figs. 3.9(g-h) depict the results obtained by applying our optimization method with different values of k (see Sec. 3.3.3.3).

Figure 3 . 9 :

 39 Figure 3.9: Experimental results for the image over New York City.

Fig. 3 .

 3 Fig. 3.10 illustrates experimental results for the image of Brest. Our method was executed with k = 600. Fig. 3.10(c) shows the obtained segmentation map. Two fragments of the map (boxed in Fig. 3.10(b)) are amplified for comparison in Fig. 3.10(d). These results validate the previous observations. The BPT is built in 13 seconds and the optimization takes 84 seconds using Propositions. 1-2 (see Section 3.3.3.2), against 214 s.

Figure 3 . 10 :

 310 Figure 3.10: Experimental results for the image over Brest.

Figure 3 . 11 :

 311 Figure 3.11: BPT optimization detects individual buildings.

  shrink in time, from sequences of extremely noisy images. Examples of growing shapes are forest fires or expanding cities in satellite images and organ development in medical imaging. In the image sequences we consider, both foreground and background intensity distributions can vary significantly over time: foreground can be heavily occluded or undistinguishable from a part of the background, and data for some pixels can be missing (see Fig.4.2 and 4.7 for examples of such sequences). Most of previously-proposed spatiotemporal methods rely on coherence of foreground/background intensity distributions in successive image frames, and are therefore not suited for segmenting such noisy data sets. Few approaches have been specifically designed for spatio-temporal segmentation of magnetic resonance image (MRI) sequences with low signal-to-noise ratio[START_REF] Riklin-Raviv | Segmentation of image ensembles via latent atlases[END_REF][START_REF] Wolz | Measurement of hippocampal atrophy using 4d graph-cut segmentation: Application to ADNI[END_REF]. Applied to multi-temporal time series that show a monotonously growing or shrinking structure, however, these smoothing methods bias results towards the mean shape obtained from averaging consecutive segmentations and, hence, underestimate rapid growth or shrinkage events.To address this issue, we have proposed a new omniscient segmentation framework based on graph cuts for the joint segmentation of a multi-temporal image sequence. It introduces growth or shrinkage constraint in graph cuts by using directed infinite links, which connect pixels at the same spatial locations in successive image frames. By minimizing a submodular energy computed on the resulting spatio-temporal graph of the image sequence, the proposed method yields a globally optimal solution. Differently from the state-of-the-art omniscient techniques, it does not rely on the coherence of the intensity in time, but only on the coherence of the shape.We have validated the performance of the proposed framework on three applications with very noisy image sequences. The first one deals with the segmentation of multiyear sea ice floes in a set of satellite images acquired through different satellite sensors. The new method returns accurate melting profiles of sea ice, which is important for building climate models. The second application segments growing burned areas from time series of optical satellite images with missing data. The third application addresses the segmentation of brain tumors from longitudinal sets of multimodal MRI volumes, where we impose additional inter-modal inclusion constraints for the joint segmentation of different image structures (brain tissues).

CHAPTER 4 . 70 (Figure 4 . 1 :

 47041 Figure 4.1: (a) Enforcing shape growth in an image sequence. (b) Segmenting jointly two sequences S1 and S2, by enforcing the foreground of S1 to contain the foreground of S2, with directed infinite links from S1 towards S2, between all pixels of coordinates (x, y, t, s 1 ) and (x, y, t, s 2 ).

  4.1(b)). An example of such application is given in Sec. 4.3.3, where image sequences correspond to four different MRI modalities, aligned both in space and time.

  and edges. Denoting by N the number of nodes and M the number of edges, this worst case complexity is O(N M 2 ) for the Edmonds-Karp algorithm, O(N 2 M ) for the Dinitz blocking flow algorithm, which goes down to O(N M log(N )) using dynamic trees.

•Figure 4 . 2 :

 42 Figure 4.2: (a, b) Reprojected MODIS images captured on the 244th and 252th days, respectively. (c, d) Corresponding AMSR-E images.

CHAPTER 4 . 77 DayFigure 4 . 3 :

 47743 Figure 4.3: (From top to bottom) Aligned images and segmentation contours (red) for four time moments (days 230, 233, 235 and 267, respectively) computed by the graph-cut methods: (From left to right) Aligned image, [Feedforward], [w/o], [Bi=variable], [Bi=4], [Bi=16], proposed method with monodirectional infinite links. Manual segmentation is shown in green. The rightmost part of the white area in the third row is not part of the object, but another ice floe who temporarily collided.

4 . 2 . 1 .Figure 4 . 4 :

 42144 Figure 4.4: Results for Application 1. (a) Mean and standard deviation for the dice score as a function of the temporal link's weight, when using mono-(green) and bidirectional (red) temporal links. (b) Area of a multiyear ice floe as a function of time, computed by using mono-and bidirectional links with different weights.

Fig. 4 .

 4 Fig. 4.5 depicts the computational time for the proposed graph-cut-based optimization as a function of the number of frames. The total computational time grows linearly with the number of frames, and is approximately twice the time that would be taken by the independent segmentation of each frame.

Figure 4 . 5 :

 45 Figure 4.5: Computational time for the proposed segmentation of the ice floe set as a function of the number of frames.

Fig. 4 .Figure 4 . 6 :

 446 Figure 4.6: Flowchart of the segmentation algorithm applied for Application 2.

Figure 4 . 7 :

 47 Figure 4.7: Two time series of MODIS images acquired in 2011 and 2013 (left and right, respectively). (a) MODIS band 5 images for three time moments: days 251, 265 and 279. Black pixels denote missing data. (b) Maps from the MCD64A1 product of areas burned during the days 213-251, 213-265 and 213-279, respectively (white pixels = burned areas). (c-d) Segmentation maps for images (a) computed by: (c) the proposed method with growth constraint; (d) [w/o] method with no temporal constraints.

Figure 4 . 8 :

 48 Figure 4.8: Results for Application 2, for two MODIS time series acquired in 2011 (a-b) and 2013 (c-d), respectively. (a, c) Mean and standard deviation for the dice score as a function of the temporal link's weight, when using mono-(green) and bidirectional (red) temporal links. (b, d) Burned area as a function of time, when using no temporal links (blue), monodirectional infinite links (green) and MCD64A1 product (red).

  We applied the described algorithm with the parameters empirically set as D = 31, r = 20, β = 2 and T ′ = 20, and with different temporal regularizations, i.e., [w/o], [Mono=0.25...∞] and [Bi=0.25...∞]. Neither [Feedforward] nor [Bi=variable] methods are suited for segmenting image sequences with missing data. Fig. 4.8(a, c) gives the re-

CHAPTER 4 . 86 -Figure 4 . 10 :

 486410 Figure 4.10: Results for Application 3. (a) Boxplots report changes in the average segmentation performance of the ten image sequences when testing different regularization approaches. (The box indicated quartiles, the whiskers outliers.). Using a strong monodirectional regularization acting as a growth prior yields the best results. (b) Volume-time plot for a patient with 14 observations (cmp. video in supplement). Solid lines indicate edema, dashed indicate tumor core that starts growing with constant rate at around day 500. The segmentation with growth constraint (red) returns results similar to the manual segmentation (green). Segmentations obtained by evaluating image volumes of each time point individually (blue) show significant variation, even obscuring the overall trend.

Figure 5 . 1 :

 51 Figure 5.1: Convolutional layer. Connections are limited to a local spatial neighborhood and the weights are shared across the different neurons.

  3x3 conv. 2x2 pooling 3x3 conv. 2x deconv.

Figure 5 . 2 :

 52 Figure 5.2: Lateral view of a fully convolutional network (dashed lines indicate inputs that have been padded in convolutional layers and cropped in the deconvolutional layer to preserve spatial dimensions).

Figure 5 . 3 :

 53 Figure 5.3: To classify the central gray pixel of this patch (and not to confuse it, e.g., with an asphalt road), we need to take into account a spatial context (a). However, we do not need a high resolution everywhere in the patch. It can be lower as we go away from the central pixel and still identify the class (b).

Figure 5 . 4 :

 54 Figure 5.4: A dilated convolution (i.e., on non-adjacent inputs) with a dilation factor S = 4.

Fig. 5

 5 Fig. 5.3.

Figure 5 . 5 :Figure 5 . 6 :

 5556 Figure 5.5: Deconvolution network. The CNN is "mirrored" to learn the deconvolution.

Figure 5 . 7 :

 57 Figure 5.7: Skip network: multiple classification scores are obtained from intermediate CNN features at different resolutions, and are combined by element-wise adding and upsampling.

Figure 5 . 8 :

 58 Figure 5.8: MLP network: intermediate CNN features are concatenated, to create a pool of features. Another network learns how to combine them to produce the final classification.

Fig. 5 .

 5 Fig. 5.8). The different features are extracted from intermediate layers of a single CNN.

  that we simply add the DSM as an extra band. In the case of Vaihingen we predict five classes, ignoring the clutter class, due to the lack of training data for that class. In the case of Potsdam we predict all six classes.

Figure 5 . 9 :Figure 5 . 10 :

 59510 Figure 5.9: Classification of closeups of Vahingen (1-3) and Potsdam (4-6) validation sets. Classes: Impervious surface (white), Building (blue), Low veget. (cyan), Tree (green), Car (yellow), Clutter (red).

Figure 5 .

 5 Figure 5.11: A CNN trained on a different dataset misclassifies most of Lake Zurich as a building.

Fig. 5 .

 5 Fig. 5.13 depicts close-ups of the classification on the test set, i.e., on regions never "seen" by the neural network at training time. While the FCN produces fuzzy results, it successfully identifies buildings in varied images. The MLP network provides finer outputs, as confirmed both numerically and visually.

(Figure 5 . 13 :

 513 Figure 5.13: Visual close-ups on test set. (a) Color input. (b) Reference data. (c) FCN results. (d) MLP results.

Figure 5 . 14 :

 514 Figure 5.14: Sample classification of buildings with a CNN.

Figure 5 . 15 :

 515 Figure 5.15: Samples of reference data for the building class. Imprecise OpenStreetMap data vs manually labeled data.

Figure 5 . 17 :

 517 Figure 5.17: Modules of Fig. 5.16 are stacked (while sharing parameters) to implement an RNN.

coverage to create a 22 Figure 5 . 18 :Figure 5 . 19 :

 22518519 Figure 5.18: Manually labeled tile used to train the RNN for the classification enhancement task.

Fig. 5 .

 5 19 provides closeups of results on different fragments of the test dataset. The initial and final maps (before and after the RNN enhancement) are depicted, as well as the intermediate results through the RNN iterations. We show both a set of final classification maps and some single-class fuzzy probability maps. We can observe that as the RNN iterations go by, the classification maps are refined and the objects better align to image edges. The fuzzy probabilities become more confident, sharpening object boundaries. To quantitatively assess this improvement we compute two measures on the test set: the overall accuracy (proportion of correctly classified pixels) and the intersection over union (IoU)[START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF]. Mean IoU has become the standard in semantic segmentation since it is more reliable in the presence of imbalanced classes (such as background class, which is included to compute the mean)[START_REF] Csurka | What is a good evaluation measure for semantic segmentation?[END_REF]. As summarized in the table of Fig.5.20(a), the performance of the original coarse CNN (denoted by cnn) is significantly improved by attaching our RNN (cnn+rnn). Both measures increase monotonously along the intermediate RNN iterations, as depicted in Fig.5.20(b).

  Evolution through RNN iterations

Figure 5 . 20 :Figure 5 . 21 :Figure 5 . 22 :

 520521522 Figure 5.20: Quantitative evaluation on Pléaiades images test set over Forez, France.

Fig. 5 .

 5 24(a) illustrates fragments of heat maps of the building and road classes, and the responses to two of the filters M i learned by the RNN. When analyzing these responses we can observe that they act as gradients in different directions, confirming the expected behavior. Fig.5.24(b) illustrates a fragment of the color image and its response to two filters N j . One of them acts as a gradient operator an the other one highlights green vegetation, suggesting that this information is used to enhance the classification maps. Filter M 1 acts like a gradient operator in the South-East direction and M 2 in the North direction (top: building, bottom: road).

N 1 N 2 (

 2 b) N 1 acts like a gradient operator in the North direction and N 2 highlights green vegetation.

Figure 5 . 24 :

 524 Figure 5.24: Feature responses (red: high, blue: low) to selected M i and N j filters, applied to the heat maps and input image respectively (see Eq. 5.11).

Chapter 6

 6 Conclusion and Future WorkI have presented my research on spectral-spatial image classification, mainly applied for automatic interpretation of remote sensing data. I have explored different strategies for this purpose, varying from the use of strong shape priors for object discrimination, regularization of classification probabilities on the image graphs, and up to the use of deep learning approaches which are capable of learning complex shape features from training data.

  ○ , 22.5 ○ , 45 ○ , 57.5 ○ , 90 ○ , 112.5 ○ , 135 ○ , 157.5 ○ }. The size of steerable filters is fixed to 21 × 21 pixels. The scale factors σ 2 and the minimum path length ˆ are adaptively selected for each dataset to obtain the best performances. The binary spatial pattern B is based on the physical attribute of the dataset in that thickness τ of curvilinear structure shows various depending on the dataset. For Structured SVM training, we sample 2000 image patches on each dataset.

  have described two classification approaches which are based on a template concept and aim to reconstruct curvilinear structures from image data. The structures are modeled as a set of line segments with different orientations. The first approach is unsupervised and based on Marked point process and Monte Carlo sampling. One of key innovations

			DRIVE RecA Aerial Cracks
		Frangi et al. [52]	0.33	0.33	0.32	0.056
		Law & Chung [110]	0.43	0.21	0.25	0.085
	(a) F 1	Becker et al. [5]	0.50	0.45	0.53	0.23
		Sironi et al. [169]	0.55	0.50	0.55	0.27
		Proposed (Graph)	0.36	0.59	0.59	0.38
		Frangi et al. [52]	37.09	16.14 22.19	28.64
		Law & Chung [110]	19.60	29.82 29.51	33.00
	(b) ρ(%) Becker et al. [5]	12.67	9.39	9.76	32.13
		Sironi et al. [169]	5.41	5.47	0.83	17.12
		Proposed (Graph)	0.17	0.57	0.35	1.07

Table 2 .

 2 1: Comparison of the quantitative performances over the datasets. We provide results with various performance measures: (a) F 1 scores, and (b) average proportion ρ of the pixels being a part of curvilinear structure. Boldfaced numbers are used to show the best score in each test.

of the proposed method is that instead of manually fixing hyperparameters, we proposed a new optimisation algorithm which seeks to maximize the consensus between various line hypotheses, obtained by simulating several Markov chains with different parameter settings.

Table 3 .

 3

	Building overlap 0.51	0.51	0.54	0.56
	Overall accuracy 0.88	0.94	0.91	0.94

1: Numerical evaluation on the Pavia Center dataset. SVM Graph cut BPT α=0 BPT α=0.5

regularize the cut by the total number of regions:

Table 4 .

 4 1: Results for Application 1. Mean and standard deviations of the Dice scores for the proposed method [Mono = ∞] and graph-cut-based approaches used for comparison.

	Method	Feedforward	w/o	Bi = variable	Bi = 16	Mono = ∞
	Dice score	.554 ± .128	.933 ± .099	.958 ± .048	.978 ± .007 .980 ± .007
	the summer period.
	We aimed at segmenting a multiyear ice floe from a 45-day sequence (summer period
	from the 227th to 271th day of 2008) of measurements over the polar regions by two
	Aqua satellite sensors: Advanced Microwave Scanning Radiometer -Earth Observing
	System (AMSR-E, 89 GHz, 6.25 km spatial resolution) and Moderate-Resolution Imaging
	Spectroradiometer (MODIS, band 1, 0.620-0.670 µm, 250 m spatial resolution). Fig. 4.2

Table 5 .

 5 1: Architecture of our base FCN.

	Layer	Filter size Number of filters Stride Padding
	Conv-1_1	5	32	2	2
	Conv-1_2	3	32	1	1
	Pool-1	2		2
	Conv-2_1	3	64	1	1
	Conv-2_2	3	64	1	1
	Pool-2	2		2
	Conv-3_1	3	96	1	1
	Conv-3_2	3	96	1	1
	Pool-3	2		2
	Conv-4_1	3	128	1	1
	Conv-4_2	3	128	1	1
	Pool_4	2		2
	Conv-Score	1	5	1

Table 5 .

 5 2: Numerical evaluation of architectures derived from our base FCN on the Vaihingen validation set.

		Imp. surf. Building Low veg.	Tree	Car	Mean F1 Overall acc.
	Base FCN	91.46	94.88	79.19	87.89	72.25	85.14	88.61
	Unpooling	91.17	95.16	79.06	87.78	69.49	84.54	88.55
	Skip	91.66	95.02	79.13	88.11	77.96	86.38	88.80
	MLP	91.69	95.24	79.44	88.12 78.42	86.58	88.92

Table 5 .

 5 3: Numerical evaluation of architectures derived from our base FCN on the Potsdam validation set.

		Imp. surf. Building Low veg. Tree	Car	Clutter Mean F1	Acc.
	Base FCN	88.33	93.97	84.11	80.30 86.13	75.35	84.70	86.20
	Unpooling	87.00	92.86	82.93	78.04 84.85	72.47	83.03	84.67
	Skip	89.27	94.21	84.73	81.23 93.47	75.18	86.35	86.89
	MLP	89.31	94.37	84.83	81.10 93.56 76.54	86.62	87.02

Table 5 .

 5 4: Comparison of MLP with other methods on the Vaihingen validation set.

		Imp. surf. Build. Low veg.	Tree	Car	F1	Acc.
	CNN+RF [146]	88.58	94.23	76.58	86.29	67.58	82.65	86.52
	CNN+RF+CRF [146]	89.10	94.30	77.36	86.25	71.91	83.78	86.89
	Deconvolution [198]						83.58	87.83
	Dilation [163]	90.19	94.49	77.69	87.24	76.77	85.28	87.70
	Dilation + CRF [163]	90.41	94.73	78.25	87.25	75.57	85.24	87.90
	MLP	91.69	95.24	79.44	88.12 78.42 86.58 88.92

Table 5 .

 5 5: Comparison of MLP with other methods on the Potsdam validation set.

		Imp. surf. Build. Low veg. Tree	Car Clutter	F1	Acc.
	Dilation [163]	86.52	90.78	83.01	78.41 90.42	68.67	82.94 84.14
	VGG pretr. [163]	89.84	93.80	85.43	83.61 88.00	74.48	85.86 87.42
	VGG+Dilation [163]	89.95	93.73	85.91	83.86 94.31 74.62 87.06 87.69
	MLP	89.31	94.37	84.83	81.10 93.56 76.54 86.62 87.02

Table 5 .

 5 6: Execution times.

		Train [s]	Test [s/ha]
		Vaih. Pots. Vaih. Pots.
	Base FCN	3.9	9.8	0.81	1.44
	Unpooling	8.4	21.0	1.38	1.84
	Skip	6.6	16.9	0.81	1.48
	MLP	10.0	24.5	1.70	2.0
	Dilation*	62	400	4.81	17.2
	*As reported in [163] (see details in Sec. 5.1.4.)

Table 5 .

 5 7: Dataset statistics. *Tile size: 1500 2 px. (0.3 m resolution).

	Train	Tiles* Total area	Test	Tiles* Total area
	Austin, TX	36	81 km 2	Bellingham, WA	36	81 km 2
	Chicago, IL	36	81 km 2	San Francisco, CA	36	81 km 2
	Kitsap County, WA	36	81 km 2	Bloomington, IN	36	81 km 2
	Vienna, Austria	36	81 km 2	Innsbruck, Austria	36	81 km 2
	West Tyrol, Austria	36	81 km 2	East Tyrol, Austria	36	81 km 2
	Total	180	405 km 2	Total	180	405 km 2

Table 5 .

 5 8: Numerical evaluation on small validation set.

			Austin Chicago Kitsap Co. West Tyrol Vienna Overall
	FCN	IoU	47.66	53.62	33.70	46.86	60.60	53.82
		Acc.	92.22	88.59	98.58	95.83	88.72	92.79
	Skip	IoU	57.87	61.13	46.43	54.91	70.51	62.97
		Acc.	93.85	90.54	98.84	96.47	91.48	94.24
	MLP IoU	61.20	61.30	51.50	57.95	72.13	64.67
		Acc.	94.20	90.43	98.92	96.66	91.87	94.42

Table 5 .

 5 9: Numerical evaluation on test set.

		Bellingham Bloomington Innsbruck S. Francisco East Tyrol Overall
	FCN IoU	44.83	35.38	36.50	44.92	43.69	42.19
	Acc.	94.48	94.07	92.97	82.60	95.14	91.85
	Skip IoU	52.91	46.08	58.12	57.84	59.03	55.82
	Acc.	95.14	94.95	95.16	86.05	96.40	93.54
	MLP IoU	56.11	50.40	61.03	61.38	62.51	59.31
	Acc.	95.37	95.27	95.37	87.00	96.61	93.93

  Figure 5.16: One enhancement iteration represented as common neural network layers.responses are found by computing the dot product between a set of filters M i and the values of u k,t (⋅) in a spatial neighborhood of a given point. Analogously, a set of filter responses are computed at the same spatial location on the input image, corresponding to the different N j * I of Eq. 5.11. These operations are convolutions when performed densely in space, N j * I and M i * u t being feature maps of the filter responses.

	u t	Conv.	. . . M i * u t	Concat.	. . .	...	δ u t	+	u t +1
						MLP
		Conv.	. . .
			N j * I
	Image I

=

  97.78 65.30 59.12 39.03 97.74 CNN+RNN 98.24 72.90 69.16 51.32 98.20

		Overall Mean	Class-specific IoU
	Method	accuracy IoU Build. Road Backg.
	CNN	96.72	48.32 38.92 9.34	96.69
	CNN+CRF	96.96	44.15 29.05 6.62	96.78
	CNN+RNN

https://www.biochem.wisc.edu/faculty/inman/empics/dna-prot.htm

http://hugues.zahlt.info/91.html

http://cvlab.epfl.ch/page-108936-en.html

While such architecture is named a "no-downsampling" network in[START_REF] Sherrah | Fully convolutional networks for dense semantic labelling of highresolution aerial imagery[END_REF], a more appropriate name would probably be "no-upsampling", because there is indeed downsampling due to the max pooling layers.

project.inria.fr/aerialimagelabeling

https://gis.tirol.gv.at/arcgis/services/Service_Public/orthofoto/MapServer/WMSServer; http://maps.wien.gv.at/wmts/1.0.0/WMTSCapabilities.xml

https://www.mapbox.com/blog/status-san-francisco-complete/

Color image

Coarse CNN CNN+CRF CNN+RNN = CNN+RNN Ground truth We also compare our RNN to the approach in [START_REF] Chen | Semantic image segmentation with deep convolutional nets and fully connected CRFs[END_REF] (here denoted by cnn+crf),

where a fully-connected CRF is coupled both to the input image and the coarse CNN output, in order to refine the predictions. This is the idea behind the so-called Deeplab network, which constitutes one of the most important current baselines in the semantic segmentation community. While the CRF itself could also be implemented as an RNN [START_REF] Zheng | Conditional random fields as recurrent neural networks[END_REF], we here stick to the original formulation because the CRF as RNN idea is only interesting if we want to train the system end to end (i.e., together with the coarse prediction network). In our case we wish to leave the coarse network as is, otherwise we risk overfitting it to this much smaller set. We thus simply use the CRF as in [START_REF] Chen | Semantic image segmentation with deep convolutional nets and fully connected CRFs[END_REF] and tune the energy parameters by performing a grid search using the enhancement set as a reference. Five iterations of inference on the fully-connected CRF were preformed in every case.

To further analyze our method, we also consider an alternative enhancement RNN in which the weights of the MLP are shared across the different classes (which we denote by "class-agnostic cnn+rnn"). This forces the system to learn the same function to update all the classes, instead of a class-specific function. Proofs related to BPT optimization I first elaborate the proofs that supplement the paragraphs that describe the space of moves in the optimization approach. I then add the proofs of the space and computational complexities of including convex hulls in Binary Partition Trees (BPTs).

Numerical results are included in the table of

A.1 Properties of prune-and-paste moves Proposition 1. Given a tree τ , suppose a node R m is pasted at τ i < τ 1 leading to a new tree ϕ. Let us consider an alternative move that pastes R m at τ j , with τ i < τ j < τ 1 , producing a tree ψ. In the cases where either

Proof. Let us abbreviate E(τ i ) as e τ i and C(τ i ) as c τ i . Following (3.11) in Chapter 3,

where c τ i denotes the sibling of c τ i (see Fig. 3.4-b). This implies:

Let us now suppose that a node R m is pasted at τ i . Then

Let us now paste R m one position upper than before. We wish to check if it is possible that this move will be better that the previous one (c ψ 1 < c ϕ 1 ):

This can be true if and only if:

In this expression, I contradicts β. Considering that e ψ y = e ϕ i-1 (see Fig. 3.4-b), the term II also contradicts β. The term IV contradicts γ. We must now analyze III. By combining III and α:

Proposition 2. Let us consider a case where Prop. 1 hypotheses do not apply. There might then exist a higher paste place τ α so that C(ψ 1 ) < C(ϕ 1 ). Let us suppose that instead of pasting at τ α we paste at τ β , with τ α < τ β < τ 1 , leading to a tree ρ. Then C(ρ 1 )

would monotonously decrease as the paste place τ β is located higher.

Proof. If Prop. 1 hypotheses do not apply, then the term III in its proof must be true.

This term implies that when pasting at τ j , the cut on the tree will be located at or below Corollary: Combining Proposition 2, where C monotonously decreases (⩽) as the paste location gets higher, and Proposition 3, where C cannot decrease (⩾), it becomes evident that pasting a node anywhere between the initial cut and the lowest common ancestor produces the same effect on the energy. Proof. Let us call N R the number of neighbors of the region R. At the lowest scale and in a discrete environment we can suppose that the number of neighbors is equal to its boundary length (δR). We are interested in knowing the number of neighbors at all scales. In a balanced tree it can be assumed that the number of neighbors at every scale is half the number at the following one. As a result:

In a discrete implementation (assuming 4-connectivity):

The summation of the neighbors of all regions in a tree T is then

Following (A.18), the total number of neighbors (the possible cut/paste moves) is a factor of n log(n).

A.2 Complexity of incorporating convex hulls

Proposition 5. The storage space required to add the convex hull to every node of a balanced BPT in a discrete environment is bounded by O(n log(n)).

Proof. Let us call CH(R) the convex hull of a region R. In the extreme case (the most compact region), CH(R) can be as large as the perimeter δR of R which, in a discrete implementation (assuming 4-connectivity) does not contain more points than four times the area of the region:

As a consequence, the points of the convex hull of all regions in a tree T must be

If we observe that in a balanced tree

then (A.17) is bounded by a factor of n log(n).

Proposition 6. The complexity of computing the convex hull of every region represented in a balanced BPT in a discrete environment, is bounded by O(n log(n)).

Proof. Let us call CH(R) the convex hull of a region R. In the extreme case (the most compact region), CH(R) can be as large as the perimeter δR of R which, in a discrete implementation (assuming 4-connectivity) does not contain more points than four times the area of the region:

The time to compute CH(R i ) is linear on the number of points in the polygons of the children:

The time to compute the convex hull of every node in the tree is then bounded by a factor of:

Following (A.18), the execution time is then a factor of n log(n).

Appendix B

Submodularity proof

W i,j can be written as a sum of terms over time t, which is, if τ i ⩽ τ j :

and

W t i,j (0, 0) if t < min(τ i , τ j ) (A) W t i,j (0, 1) if τ j ⩽ t < τ i (C) W t i,j (1, 1) if t ⩾ max(τ i , τ j ) (D) otherwise.

We can thus represent W i,j (τ i , τ j ) by the sequence of cases A, B, C or D chosen for each t:

-t τ i τ j W i,j (τ i , τ j ) = t AAAAA BBBBBBB DDDDDD when τ i ⩽ τ j , and otherwise :

-t τ j τ i W i,j (τ i , τ j ) = t AAAAA CCCCCCC DDDDDD Now, to represent the submodularity condition, let τ i , τ j , τ ′ i , τ ′ j be any times, satisfying τ i ⩽ τ ′ i and τ j ⩽ τ ′ j . We will suppose moreover that τ i ⩽ τ j : otherwise consider W j,i (τ j , τ i )