Je tiens à remercier tout particulièrement Sixin Zhang, qui a joué le rôle de co-directeur officieux durant mon odctorat. Il a fait preuve à mon égard d'une infinie patience, d'un soutien sans faille, et de sagesse en toute circonstance.

Je suis également reconnaissant envers mes deux directeurs de thèse, Bartek Błaszczyszyn et Stéphane Mallat, pour leurs encouragements, leur enthousiasme, leur expérience, leur rigueur, et leur curiosité. Ils ont fait de mes années de doctorat une expérience plus riche que je ne pouvais l'imaginer.

Je remercie vivement les membres du jury, et particulièrement les rapporteurs, Frédéric Lavancier et Kai Schneider, pour le temps et l'énergie qu'ils ont consacrés à la lecture de mes travaux, l'intérêt qu'ils y ont porté, et leurs judicieuses remarques. Merci à Brice Ménard pour la confiance qu'il m'a accordé, et son excellent pain au maîs.

Merci à tous les collègues de l'équipe Data: Florentin ma nounou du bureau, et ses claquements de dents dans la cour en hiver pour ne pas me laisser seul, Rudy et ses blagues sur ma calvitie naissante, Gaspar et ses conseils en menuiserie-joaillerie, Simon et son analyse ciselée de mes tenues vestimentaires, John et sa croisade pour la fresque du CSD, et maintenant Etienne, qui m'a montré qu'il y a même des poissons rouges qui jouent mieux au foot que moi. Les anciens également, Tomás, Louis, et Tanguy.

Je remercie Marc Hoffmann, pour m'avoir soutenu et guidé dans des temps difficiles. Par dessus tout, je tiens à remercier ma famille, pour n'avoir jamais perdu espoir; et Charlotte, pour m'avoir montré qu'au fond, le plus important, c'est les pâtes à la sauce tomate. i

Context and motivations

This dissertation discusses the study of different representations of random geometric objects. An elementary class of objects described by their geometry consists of point patterns, i.e. collections of points living in some metric space. In practical applications, these points usually represent the temporal or spatial location of events, or entities of interest. Without formalism, a point pattern is represented by a sum of Dirac atoms, i.e.

φ = N i=1 δ u i ,
where N ∈ N ∪ {∞}, and the u i represent the locations of these entities. When point patterns are considered to be observations of a random phenomenon, this phenomenon is called a point process. Point processes arise in a wide variety of contexts. Time-related processes can represent, for example, the recorded times of financial orders, the spiking times of a neuron, or the times of arrival of customers on a website. For spatial point processes, one can think of the positions of trees in forests, cells in a biological tissue, individuals in a given population, or stars, planets, galaxies in space. One can also study spatio-temporal objects such as the time and position of earthquakes in a given area, or the dynamics of droplets in a cloud. Depending on the phenomenon being observed, the points, or atoms of the process, can tend to repel each other, forming rather regular structures (for instance, trees in a forest tend to grow away from each other, in order to have the most available nutriments in the ground, see e.g. [START_REF] Stoyan | Recent applications of point process methods in forestry statistics[END_REF]). Conversely, they can also regroup together, forming clusters of points of possibly complex geometric structures (for example, galaxies in the universe can tend to group together, and form filamentary structures, see e.g. [START_REF] Tempel | Detecting filamentary pattern in the cosmic web: a catalogue of filaments for the sdss[END_REF]).

Being able to understand, and adequately describe these objects is thus of central interest in many areas of science, whether it is to perform prediction of future events, estimation of biological or physical parameters, or modelling. Recent advances in statistical learning, and particularly deep learning, have led to impressive progress in the representation and manipulation of complex objects, yielding great performances on a wide range of tasks, such as classification, regression, segmentation, or generation. These methods are able to handle a vast variety of complex objects (e.g. images, audio signals, languages), including point processes (e.g. [START_REF] Qi | Pointnet: Deep learning on point sets for 3d classification and segmentation[END_REF][START_REF] Xiao | Wasserstein learning of deep generative point process models[END_REF]). However, representations defined by (deep) neural networks are often referred to as 'black boxes', as we do not fully understand what aspects of the geometry these networks learn, which make these representation lack interpretability, and possibly reliability. On the other hand, we find in the literature on point processes many different mathematical representations, called summary characteristics, that each capture a specific aspect of the geometry of point patterns. These characteristics usually rely either on second-order statistics (e.g. the pair correlation function, Ripley's K-function [START_REF] Ripley | The second-order analysis of stationary point processes[END_REF], or Besag's L-function [START_REF] Besag | Contribution to the discussion on dr ripley's paper[END_REF]), or on distances between neighbouring points (such as the spherical contact distribution function or the nearest neighbours distance distribution function). None of them fully capture the entire range of geometric structures that point processes may exhibit, and it can be challenging, depending on the data that one is dealing with, and the task at hand, to find a suitable representation.

Finding concise, easily manipulable representations of complex data is also an active area of research and a central interest in the field of signal processing, where many tools have been developed (see e.g. Mallat [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF]) to extract important information about the data with a relatively 'compact' (or low-dimensional) representation. A useful way to build such a representation is to leverage the possible regularities or invariances present in the data, for which a fundamental tool is the Fourier analysis. The Fourier transform builds a representation of signals based on their frequency, i.e. the repetition of patterns in time or space, with complex sinusoidal functions. However, these oscillating functions are 'delocalized' in space, so they may no be adapted to capture sudden changes, or sharp local transitions in the data, along the time or space axis. To capture such local structures in signals, one can instead compute their correlations with functions that are 'localized' both in space (or time) and frequency, such as wavelets. Wavelet analysis has been extensively studied in the literature, for the processing of many types of signals and stochastic processes, including point processes (see e.g. [START_REF] Abry | Wavelet analysis of long-range-dependent traffic[END_REF][START_REF] Lowen | Estimation and simulation of fractal stochastic point processes[END_REF], in the context of one-dimensional point processes exhibiting fractal behaviour). In this work, we shall present different wavelet-based representations of two-dimensional stationary point processes that, similarly to classical summary characteristics, are designed to capture and compactly represent the geometric structures formed by the data. We study their quality as statistical descriptors through several learning and modelling tasks.

Related works and our contributions

We give an overview of the problems considered in this dissertation, and the current methods that have been developed to address these problems. Relatively to these works, we highlight our contributions, that we shall detail in the following chapters. More specific related works can be found in the introduction of each respective chapter.

Fourier analysis and wavelet-based representations of point processes

A fundamental theory for the analysis of stationary stochastic processes is the Fourier analysis. It is particularly useful, for instance, to study second-order moments, as the covariance matrix of a stationary process is diagonal in the Fourier basis. In the context of point processes, this has been formalized by the covariance measure in the Fourier domain, called the Bartlett spectrum (Bartlett [START_REF] Bartlett | The spectral analysis of point processes[END_REF], [START_REF] Bartlett | The spectral analysis of two-dimensional point processes[END_REF][START_REF] Brémaud | Fourier analysis of stochastic processes[END_REF]), and can be estimated through the Discrete Fourier Transform of observation samples, as we show in Section 2.8.2 of this work. Figure 1.1 illustrates this estimation on a distribution for which the theoretical value of the Bartlett spectrum is known. In order to extract meaningful information about the structure of signals beyond second-order properties, [START_REF] Mallat | Group invariant scattering[END_REF] introduce the scattering moments, which compute a cascade of convolutions of the input signal with wavelets, followed by a complex modulus operator. Informally (see Section 3.2 for a more precise definition), given a zero-mean complex function ψ, a wavelet family can be constructed by dilating and rotating ψ, i.e. ψ j,θ (u) = 2 -2j ψ(2 -j r θ u).

The scattering moments define a representation of an input signal φ, observed in a finite window W , by computing the values

S 1 φ(j, θ) = 1 |W | W |φ ψ j,θ |(u)du S 2 φ(j 1 , θ 1 , j 2 , θ 2 ) = 1 |W | W ||φ ψ j 1 ,θ 1 | ψ j 2 ,θ 2 |(u)du,
where denotes the convolution operator. If Φ is a stationary point process, then the scattering moments of its distribution are defined in expectation, i.e.

E[S 1 Φ(j, θ)]

E[S 2 Φ(j 1 , θ 1 , j 2 , θ 2 )].

In [START_REF] Bruna | Intermittent process analysis with scattering moments[END_REF], the authors study the properties of scattering moments of the one-dimensional homogeneous Poisson point process. We extend some of these results, in the twodimensional case, to a broader class of stationary point processes. In particular, we study the limit behaviour of scattering moments (of first order, see Section 3.2), when the wavelet is asymptotically contracted or dilated. These limits quantify, respectively, the local interactions between points, and their macroscopic growth of variance. In Section 3.3, we show that under some technical conditions:

• When the wavelet ψ j,θ is asymptotically contracted (i.e. j → -∞), then E[S 1 Φ(j, θ)] is proportional to the intensity of the process (i.e. its expected number of atoms in a domain of unit volume). We show that this convergence is controlled by a quantity that depends on the second order moments of Φ.

• When the wavelet ψ j,θ is asymptotically dilated (j → ∞), the limit behaviour of 2 j E[S 1 Φ(j, θ)] is determined by the growth of variance of the number of atoms of Φ in domains of growing volume.

This result, illustrated in Figure 1.2 is also extended to a certain class of marked point process, where to each atom of a point process is attached some numerical value. Furthermore, we give an expansion formula for this class of marked point processes, in the spirit of [START_REF] Błaszczyszyn | Factorial moment expansion for stochastic systems[END_REF][START_REF] Błaszczyszyn | A note on expansion for functionals of spatial marked point processes[END_REF].

Estimating geometric marks

Estimating numerical values associated with point patterns based on their geometry is a well known problem that has been studied, for instance, in quantum chemistry to estimate the energy of molecules from the geometry of its atoms (see e.g. [START_REF] Eickenberg | Solid harmonic wavelet scattering for predictions of molecule properties[END_REF]). In this work, we shall seek to estimate the values of marks of individual points of a given point pattern, with these marks depending on the geometry of the pattern. More precisely, we are interested in measures of the form φ = i m(u i , φ)δ u i , where φ = i δ u i is a point pattern, and m is an unknown marking function. The goal is to estimate m based on a labeled training set. This requires to have a representation not only for the point patterns, but also for the association of the points and their marks. For this problem, we propose the following method: we build a representation Sφ of the unmarked point patterns, and S φ for the marked patterns. From the training set, we approximate the relation between Sφ and S φ:

S φ f (Sφ).
Then, with the knowledge of φ and the approximated S φ, we use an optimization algorithm to reconstruct φ. Our representations S and S are the scattering moments of the point patterns. Figure 1.3 illustrate our proposed method.

Summary characteristics and probabilistic modelling

One of the tasks that we shall study in this dissertation is the modelling of distributions of point processes, using summary characteristics. When dealing with some observation data, one may find itself in need to model the underlying distribution from this data, and sample from this distribution, without using classical parametric model (e.g. if no classical parametric model fit the data satisfyingly, or if such models are hard to sample from). Figure 1.4 (top row) shows examples of observations from distributions that we shall seek to approximate. In this situation, a possible way to model the distribution is to compute summary characteristics of the observation data, and to generate new samples that match the statistics of the observation. In more details, given an realization φ of a point process Φ, these methods consist in computing a representation K( φ), such that

K( φ) E[K(Φ)].
If another point pattern φ satisfies K(φ) K( φ), then K(φ) E[K(Φ)]. If K is a good descriptor of the distribution of Φ, then φ is likely to be a realization of Φ. Such models, which can be formalized as maximum entropy models, also constitute a convenient way to assess wether the summary characteristics capture enough information about the geometry of the data: intuitively, if the generated samples match the statistics that were used to define the model, but do not resemble the observed data, it could indicate that these statistics fail to capture important structural information. The challenges in defining maximum entropy models thus reside in defining a descriptor K that characterizes the distribution of Φ, and finding a sampling method to generate new point patterns which descriptors match the ones of the observation φ.

In the context of point processes, the work of Tscheschel and Stoyan [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF] studies the use of several classical summary characteristics, such as the k-nearest neighbours distance distribution function, to represent the distribution of point processes. In order to generate new samples, this method relies on random search, by sequentially adding and removing points from a Poisson point process sample. While their solution is efficient to generate samples containing a few hundred points, forming simple geometric structures, it is not adapted to model distributions for which samples are formed by a large number of points, exhibiting possibly long-range interactions. In this work, we propose to represent point process distributions with wavelet phase harmonics covariances. These moments were introduced in [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF], and used in [START_REF] Zhang | Maximum entropy models from phase harmonic covariances[END_REF] to model non-Gaussian stationary processes. They Observation Synthesis Spectrum Figure 1.4: Top: samples of several point processes distributions exhibiting interactions at multiple scales. Middle: samples from our model, for each respective distribution. Bottom: Bartlett spectrum density curves for the true (full blue lines) and model (dashed orange lines) distributions. are defined by computing the correlations between phase adjusted wavelet coefficients. The phase harmonics of a complex number z are defined by [z] k = |z|e ikϕ (z) , where |z| denotes the complex modulus of z, and ϕ(z) its argument. Wavelet phase harmonics covariances of a point pattern Φ are defined by the correlations between the phase harmonics of the wavelet coefficients of Φ, i.e. moments of the form

K(φ) = [φ ψ j,θ ] k (u)[φ ψ j ,θ ] k (u -τ )du.
To generate new samples more efficiently than with a random search method, we propose a gradient descent algorithm, that optimizes the positions of the atoms of a sample from a Poisson distribution, in order to minimize the energy

E φ(φ) = 1 2 |K( φ) -K(φ)| 2 .
By noting φ = i δ u i , the positions u i are updated with

u i ← u i -ζ∇ u i E φ(φ),
for some gradient step ζ > 0. Examples of samples from this method are given in Figure 1.4 (middle row), as well as the comparison between the Bartlett densities of the true and model distributions (bottom row). We provide an open-source software implementing the described method.

The overall quality of a model can be defined by its ability to produce patterns similar to the observation, but also with a diversity of samples comparable to that of the true distribution. To evaluate our model with respect to these two aspects, we propose the use of Topological Data Analysis (see e.g. [START_REF] Chazal | An introduction to topological data analysis: fundamental and practical aspects for data scientists[END_REF]). This method is based on the representation of samples from the true distribution, and samples from the model, by their persistence diagrams. Intuitively, persistence diagrams are morphological descriptors, that describe the size of connected components and 'holes' that samples exhibit, when their atoms are replaced by balls of varying size. This representation is then used to compute a notion of 'distance' between the different samples. It allows one to evaluate whether:

• The samples from the model are 'close' to samples from the true distribution.

• The samples from the model are not too close to the single observed sample, and exhibit intra-class distances similar to the ones from the true distribution.

A dimension reduction method is used to approximate this distance matrix with a set of two-dimensional points, in order to visually appreciate the two aforementioned properties of the model.

Modelling texture images

Another type of data defined by complex geometric structures, that can be seen as realizations of a two-dimensional stationary process, are spatially homogeneous images called textures. As for two-dimensional point processes, an approximation of the distribution underlying a given texture can be defined by a maximum entropy model, based on a set of statistical constraints. One of the principal evaluation methods of maximum entropy models for textures is visual inspection, as the human eye is often efficient in distinguishing between 'real' and generated samples. Finding a set of statistics that is able to represent and model texture images is a challenging task, due to the high complexity of geometries that they can exhibit, as illustrated in Figure 1.5 (top row).

In order to assess the ability of a given descriptor to capture geometric information, it is thus interesting to turn to this class of processes. In this framework, the pioneer work of Portilla and Simoncelli [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] on wavelet-based representations has led to impressive results. However, recent models relying on information captured by deep convolutional neural networks (CNN) greatly outperform these representations. In this dissertation, we define a wavelet-based model allowing to generate synthesized images of similar visual quality to such CNN-based representations. This model is defined with a representation that, similarly to wavelet phase harmonics covariances, computes the correlations between modified wavelet coefficients: for a real number α ∈ [0, 2π], we denote the generalized rectifier of a complex number z by ρ α (z) = max(0, Real(ze iα )).

Given a texture image x, we define its representation by the set of coefficients

C(x) = ρ α (x ψ j,θ )(u)ρ α (x ψ j ,θ )(u -τ )du,
for an appropriate set of indices (j, θ, α, j , θ , α , τ ). This set of indices, of coefficients that we choose to correlate, has a crucial impact on the model. We observe that if this set is too small, then the representation may fail to capture important structural information about the image. Conversely, if the index set is too large, the model may memorize parts of the observations, which will be reproduced in the syntheses. We illustrate this trade-off between visual quality and diversity of the syntheses, and define a model that lies between these two extreme cases. Figure 1.5 (bottom row) shows samples of syntheses from our model. We provide an open-source software implementing the described method.

Outline of the dissertation

In order to present our findings, we first need to briefly review the main notions and results from stochastic geometry on which we base our study. Chapter 2 aims to present our principal object of interest: point processes. More specifically, we shall focus on stationary and ergodic point processes on R 2 . As we take interest in representations of such processes, we review classic descriptors, such as second-order moments and other common summary characteristics. We then review maximum entropy models, a particular class of models for point processes that rely on such representations. This leads us to present in Chapter 3 the family of representations that we study in this work, based on the wavelet transform of counting measures. We shall focus on two types of representations: scattering moments, that are constructed by applying a cascade of wavelet convolutions and modulus operators to the input signal, and phase harmonics covariances, which compute correlations between phase-adjusted wavelet transforms of the signal. We derive theoretical results for the behaviour of scattering moments to gain some insights about the information that they capture.

We can then study the adequacy of these representations in the context of statistical learning or modelling tasks. The first experiment, based on [START_REF] Brochard | Statistical learning of geometric characteristics of wireless networks[END_REF], and detailed in Chapter 4, consists in the estimation of geometric marks (i.e. numerical values, associated with each atom of a point process, that depend on the geometry formed by the atoms). The marking function is not known, and shall be estimated from a label training set of realizations of a marked point process (i.e. the data consists of several point patterns, and the label is the associated marks). We propose the use of scattering moments as a representation of both point patterns and marked point patterns. We estimate the scattering moments of marked patterns from the scattering moments of non-marked patterns, and then estimate the marks from the moments of the marked pattern. We compare our method with a more direct approach, using raw geometric information captured by local pair-wise distance matrices, and observe that our method is more adapted to recover geometric marks when the marking function has long-range dependencies.

These representations, similarly to classic summary characteristics, describe the geometry of point patterns. As such, they can also characterize aspects of the distribution of a point process, by considering them as empirical estimators of moments. Based on [START_REF] Brochard | Particle gradient descent model for point process generation[END_REF], we study in Chapter 5 a maximum entropy model built on the description of point process distributions with wavelet phase harmonics covariances. We also present a gradient descent algorithm to efficiently sample from the model. Our findings indicate that wavelet phase harmonics covariances offer a more adequate representation than classic summary characteristics, such as nearest neighbours distances, to represent the distributions of point processes in this context.

Defining a set of statistics that characterize the distribution of a stochastic process is of central importance when defining a maximum entropy model. Following [START_REF] Brochard | Generalized rectifier wavelet covariance models for texture synthesis[END_REF], we study in Chapter 6 a possible improvement of the wavelet phase harmonics covariances. We focus in this chapter on textures, which are stationary random images, that can be seen as realizations of random signed measures, and present complex geometric structures. We show that, using an appropriate set of wavelet-based covariance moments, one can build a maximum entropy model that gives results of similar visual quality to models built on state-of-the-art CNN-based representations, on a wide range of texture images.

Notations

We review the principal notations that shall be frequently used in this work.

• (E, B) shall denote the Borel space associated with some Polish space E, on which point processes are defined in Chapter 2. Our particular focus will be on E = R 2 .

• (Ω, F, P) denotes some probability space.

• M denotes the measurable set of counting measures, and M the measurable set of locally finite measures.

• For a Borel set (E, B) and a measure M on this space, we shall denote by

L p (E, M ) := {f : E → C, E |f | p dM < ∞}. When E = R d , for some integer d ≥ 1,
and M is the Lebesgue measure on E, we shall simply denote this space L p (R d ).

• For an integer N ≥ 1,

X N := {0, • • • , N -1} 2
is the set on which gray-scale images are defined, as elements of R N ×N . Color images are elements of R 3×N ×N .

• For any s > 0,

W s := [-s 2 , s 2 ] 2 .
It shall denote the window of R 2 in which our point processes will be observed.

• For any

d ≥ 1, any two functions f, g ∈ L 1 (R d ) ∪ L 2 (R d ), the convolution of f and g is denoted by f g := R d f (u)g(• -u)du.
• For any function f from R d to C, we note f the function defined, for any u ∈ R d , by f (u) := f (-u).

• For any z ∈ C, or any complex valued function f , z * (respectively f * ) will denote its complex conjugate.

• For any z ∈ C, |z| denotes its complex modulus, and ϕ(z) its argument, i.e. z = |z|e iϕ(z) .

• For

f ∈ L 1 (R d ) ∪ L 2 (R d ), its Fourier transform, noted f , is defined, for ω ∈ R d , by f (ω) := R d f (u)e -iωu du.
Chapter 2

Stochastic geometry

In this first chapter, we recall the basic definitions and results of point process theory that we shall need for the remainder of the work. This is by no means an exhaustive introduction to the overall theory, as it simply highlights the notions that reappear in the chapters that follow. For a more thorough review of the fundamentals about the theory of point processes and random measures, and in particular, proofs of basic results, we refer the reader to [START_REF] Daley | An introduction to the theory of point processes: volume I: elementary theory and methods[END_REF][START_REF] Daley | An introduction to the theory of point processes: volume II[END_REF] or [START_REF] Baccelli | Random measures, point processes, and stochastic geometry[END_REF].

Point process framework

We begin by setting up the mathematical framework that enables us to define and study the objects called point processes and random measures. Informally, a point process can be seen as a random collection of points, living in some topological space E. The aim of this section is to introduce briefly the mathematical framework used to study these objects.

While the topological space E may be more general (one can consider a topological, locally compact, second countable, Hausdorff space), let us assume in this work that it is the d-dimensional Euclidean space R d , for some integer d ≥ 1. More precisely, our applications presented in Chapters 3 to 6 focus on the case E = R 2 . However, most of the definitions and result remain valid for more general spaces, so we shall keep as much as possible the notation E. When the properties of the space R 2 are needed, it will be explicitly stated in the text. E is endowed with its Borel σ-algebra noted B. Let us denote by B c the set of all topologically bounded (i.e. with compact closure) elements of B.

To study a collection of points of E, it is convenient to view it as a sum of Dirac atoms at the locations of the points. In other words, an at most countable, possibly infinite subset {u 1 , . . . , u N } ⊂ E, for some N ∈ N\{0} ∪ ∞, will be identified with the measure

µ = N n=1 δ un .
(2.1)

The collections of points considered in point process theory usually are locally finite, i.e. finite on all bounded Borel sets B ∈ B c . A counting measure on (E, B) is a locally finite measure taking values values in N ∪ ∞. M denotes the set of all counting measures on (E, B). More generally, the set of locally finite measures on (E, B) will be noted M. It is known that any non-null measure µ ∈ M can be expressed as a sum of Dirac atoms located at some points in E, without accumulation points. The link between the collection of points and the counting measure is clear: for any u ∈ E, µ(u) counts the number of points in the collection that are equal to u. One can informally write u ∈ µ to mean µ(u) ≥ 1. In order to consider random variables taking values in M, one needs to endow M with a σ-algebra. Let M be the one generated by the mappings

µ → µ(B), B ∈ B C .
A point process is a measurable mapping from some probability space (Ω, F, P) to the space (M, M). The distribution P Φ of a point process Φ is image of P by Φ, i.e.

∀ Γ ∈ M, P Φ (Γ) = P({Φ ∈ Γ}).

(2.2)

Note that a point process Φ can be viewed as a stochastic process: Φ = {Φ(B)} B∈Bc . Therefore, by the Kolomogorov extension theorem, the distribution of a point process is characterized by its finite dimensional distributions (Φ(B 1 ), Φ(B 2 ), . . . ,

Φ(B k )) k≥1,B 1 ,...,B k ∈Bc .
This notion can be extended to define random measures: let M the space of all locally finite measures on (E, B), endowed with the σ-algebra M, defined similarly to M. A random measure is a measurable mapping Π : (Ω, F, P) → ( M, M).

Basic characteristics of point processes

Many tools and functions have been developed to describe the distributions of point processes. In this section, we review several of them, that shall be useful fro the remainder of this dissertation.

Moment measures

As point processes are random measures, it is natural to look at the expected value of such measures on Borel sets of E, or products of such measures. This can be formalized by the definitions that follow. Definition 2.2.1 (Mean measure). The mean measure M Φ of a point process Φ is defined on (E, B) by

M Φ (B) = E[Φ(B)], B ∈ B. (2.3)
This definition allows us to state one of the fundamental results in the theory of point processes.

Theorem 2.2.2 (Campbell's averaging formula). Let Φ be a point process with mean measure M . For any measurable function f : E → R, non-negative or integrable with respect to M , the integral E f (u)Φ(du) is a.s. well defined, and

E E f (u)Φ(du) = E f (u)M (du).
(2.4)

The mean measure is simply the average number of atoms of Φ present in the set B. This notion can be generalized by looking at the joint number of points present in n-tuples of Borel sets. Definition 2.2.3 (High order moment measures). For any n ≥ 1, the n-th order moment measure of a point process Φ, noted M Φ n , is the measure on (E ×n , B ⊗n ) defined, for all

B 1 , B 2 , • • • , B n ∈ B by: M Φ n (B 1 × B 2 × • • • × B n ) := E[Φ(B 1 )Φ(B 2 ) • • • Φ(B n )].
(2.5)

It is equivalent to the mean measure of the point process Φ n on (E ×n , B ⊗n ) defined by

Φ n (B 1 × B 2 × • • • × B n ) = Φ(B 1 )Φ(B 2 ) • • • Φ(B n ).
For sets B 1 , B 2 , • • • , B n that are not mutually disjoints, atoms of Φ that lie in the intersections of those sets are counted more than once by the n-th order moment measure. To take this into account, one can use instead the following modification. Definition 2.2.4 (Factorial moment measures). Let Φ = i δ u i be a point process and n ∈ N, n ≥ 1. The n-th factorial power of Φ, noted Φ (n) , is the random measure on (E ×n , B ⊗n ) defined by

Φ (n) := ∆ n Φ δ (u i,1 ,u i,2 ,••• ,u i,n ) , ( 2.6)
where

∆ n Φ := {(u i,1 , u i,2 , • • • , u i,n ), u i,j ∈ Φ, ∀j, k, j = k ⇒ u i,j = u i,k }.
The n-the factorial moment measure is defined by:

M Φ (n) (B 1 × B 2 × • • • × B n ) := E[Φ (n) (B 1 × B 2 × • • • × B n )].
(2.7)

Furthermore, if M Φ (n) admits a density w.r.t. the Lebesgue measure, it shall be noted

ρ (n) (x 1 , • • • , x n ).
The joint factorial moment measures characterize the distribution of a point process Φ, provided it has exponential moments, i.e.

∀ B ∈ B C , ∃ : E[e Φ(B) ] < ∞.

Other common characteristics

Other functions can be used to characterize the distribution of a point process. Among them are the two following.

Void probability

The void probability v Φ of a point process Φ is defined by

v Φ (B) = P(Φ(B) = 0), (2.8)
for all B ∈ B. The distribution of a simple point process (i.e. such that P(∀u ∈ E, Φ({u}) ≤ 1) = 1, meaning that, almost surely, Φ has no multiple atoms in E) is fully characterized by the family (v Φ (B)) B∈Bc .

Laplace transform

The Laplace transform L Φ of a point process Φ takes a nonnegative, measurable function f and maps it to R + by

L Φ (f ) := E[e -f dΦ ]. (2.9)
One can show that the Laplace transform of a point process fully characterizes its distribution.

Poisson, clustering, and repulsive distributions

As previously explained, point processes are random collections of points in some metric space. Depending on their distributions, realizations of point processes can exhibit different geometric patterns. The points can tend to regroup together or repel each other, possibly forming complex geometric structures. While there exist many classic distributions in the literature on point processes, we shall, in this section, briefly detail some of the most usual ones, in order to give a panorama of these objects. For a more complete overview of the 'zoology' of point process distributions, see for example [START_REF] Diggle | Statistical analysis of spatial and spatio-temporal point patterns[END_REF].

Poisson point process

The most fundamental example of point process distribution is the Poisson distribution. It describes a process where points are 'placed independently of each other', and identically on the space. It is formally defined as follows.

Definition 2.3.1 (Poisson point process). Let Λ be a deterministic, locally finite measure on (E, B). A point process Φ is a Poisson point process with intensity Λ if it satisfies the following conditions:

1. For any B ∈ B c , Φ(B) is a Poisson random variable with intensity Λ(B).

2. For any k ≥ 1, and for any B 1 , . . . , B k ∈ B c , pairwise disjoint, the random variables Φ(B 1 ), . . . , Φ(B k ) are independent.

Its simple hypotheses (independence, Poissonianity) make the Poisson distribution highly popular, and its properties have been extensively studied (see e.g. [START_REF] Last | Lectures on the Poisson process[END_REF]).

Some properties of the Poisson point process

• Its mean measure M is equal to its intensity measure Λ. We say that the Poisson point process is homogeneous if its intensity measure is proportional to the Lebesgue measure on E.

• Its void probability is µ(B) = e -Λ(B) .

• Its Laplace transform is equal to

L(f ) = e -E 1-e -f (u) Λ(du) .

Other classic distributions

The Poisson point process is sometimes said to have complete independence, as the number of points in any disjoint regions of the space are independent. It is usually considered as a reference process, and other distributions, that exhibit spatial structures, are usually classified in one of the two following categories: clustering or repulsive. The clustering of points refers to their tendency to form groups, well separated from each other, whereas a point process is said to be repulsive if the points tend to repel each other, and are more homogeneously spread in the space. Figure 2.1 displays examples of different point process geometries. In this section, we review several important distributions of point process exhibiting such structures. For more details about comparison methods between point process distributions, see [START_REF] Błaszczyszyn | Clustering comparison of point processes, with applications to random geometric models[END_REF]. Gibbs point processes A Gibbs process is a point process whose probability distribution has a density with respect to the probability distribution of a Poisson point process.

Cox point processes

Informally, one obtains a Gibbs point process by 'filtering' realizations of a Poisson point process, giving more probability to appear to some configuration than others. The clustering or repulsive behaviour of the distribution will depend on the filtering of the configurations. A Gibbs point process can be defined by considering a Poisson process Φ, and some measurable function from

M to R + , such that E[f (Φ)] = 1.
The point process Φ f with distribution defined by P Φ f (dµ) = f (µ)P Φ (dµ) is a Gibbs point process with density f with respect to Φ1 . For a more detailed introduction to Gibbs point processes, see e.g. [START_REF] Dereudre | Introduction to the theory of gibbs point processes[END_REF].

Hard-core point processes Hard-core processes are defined by the property that all atoms of its configurations are further away from each other than some minimal distance. They are called hard-core because the points can represent the center of solid spheres non smaller than some value. Classic examples of hard-core processes are Matérn I, II, and III point processes (see e.g. [START_REF] Matern | Spatial variation. stochastic models and their application to some problems in forest surveys and other sampling investigations[END_REF][START_REF] Matérn | Lecture notes in statistics[END_REF]).

Determinantal point processes

Let Φ be a simple point process on (E, B), Λ a locally finite measure on (E, B), and k : E 2 → C a measurable mapping. Φ is a determinantal point process with background measure Λ if for all n ≥ 1, the n-th order factorial moment measure M Φ (n) of Φ admits a density ρ

(n)
Λ with respect to the product measure Λ n , which equals , for some positive constants α, β, and u, v ∈ R d (see e.g. [START_REF] Lavancier | Determinantal point process models and statistical inference[END_REF]), and the Ginibre point process, with kernel of the form k(u, v) = 1 π e uv * e -(|u| 2 +|v| 2 )/2 , for u, v ∈ C (see e.g. [START_REF] Ginibre | Statistical ensembles of complex, quaternion, and real matrices[END_REF]), both considered w.r.t. the Lebesgue measure Λ(dx) = dx.

ρ (n) Λ (u 1 , • • • , u n ) = det (k(u i , u j ) 1≤i,j≤n , for Λ n -almost all (u 1 , • • • , u n ) ∈ E n . Interpreting Λ as the intensity measure of a Poisson point process Φ 0 , M Φ (n) has density ρ (n) Λ w.r.t. M Φ (n) 0 ([6,

Palm distributions of point processes

When studying point processes, one might consider high order moments, introduced in Section 2.2, that are the expected number of points that would appear simultaneously in several given regions of the space. It may also be of interest to study the conditional number of points in a given region, knowing the number of points in the other region. Going further, on may enquire about the conditional distribution of the point process, knowing that this process has an atom at a given position u ∈ E. However, if u is not a fixed atom of Φ, P(u ∈ Φ) = 0, so one cannot directly consider conditional distribution

P(•|u ∈ Φ).
(E × M, B × M) by C(B × Γ) := E[ E 1 B (u)1 Γ (Φ)Φ(du)], B ∈ B, Γ ∈ M.
We immediately remark that, for any Γ ∈ M, C(B×Γ) ≤ M (B). The Radon-Nikodym theorem tells us that there exists a measurable function u → P u (Γ) such that,

∀B ∈ B, C(B × Γ) = E 1 B (u)P u (Γ)M (du).
Note that for a fixed Γ ∈ M, P u (Γ) is not uniquely defined, and can differ on any set of null M measure. As E is a polish space, P • (•) can be defined as a probability kernel from (E, B) to (M, M), implying that for any u ∈ E, P u (•) is a probability distribution on (M, M), called a Palm distribution of Φ. This definition allows us to state the following results, that extends Campbell's averaging formula to stochastic processes. Theorem 2.4.2 (Campbell-Little-Mecke (CLM)). Let Φ be a point process on E, with Radon mean measure M , and Palm distributions P u , for all u ∈ E. For all non-negative measurable functions f :

E × M → R, E[ E f (u, Φ)Φ(du)] = E×M f (u, µ)P u (dµ)M (du).
(2.10)

The result extends to all functions f for which either of the two sides of the equality eq. (2.10) is finite when f is replaced by |f |.

For any u ∈ E, the quantity M f (u, µ)P u (dµ) shall also be noted E u [f (u, Φ)] (that is, the expectation of f (u, Φ) for some point process Φ having distribution P u on M).

Reduced and high-order Palm distributions

The reduced version of the Palm distributions, noted P ! u can be defined as the Radon-Nikodym derivative of

C (1) (B × Γ) := E[ E 1 B (u)1 Γ (Φ -δ u )Φ(du)], B ∈ B, Γ ∈ M.
(2.11)

The Palm distributions can be extended to higher order Palm distributions P u 1 ,••• ,un by taking the derivative of

C n (B × Γ) := E[ E n 1 B (u)1 Γ (Φ)Φ (n) (du)], B ∈ B ⊗n , Γ ∈ M, (2.12)
and similarly for the reduced high-order versions P ! u 1 ,••• ,un , by taking the derivative of:

C (n) (B × Γ) := E[ E n 1 B (u)1 Γ (Φ - n i=1 δ u i )Φ (n) (du)], B ∈ B ⊗n , Γ ∈ M. (2.13) The quantity M f (u, µ)P u 1 ••• ,un (dµ) shall also be noted E u 1 ••• ,un [f (u, Φ)]
, and the reduced versions

E !u 1 ••• ,un [f (u, Φ)].

Stationary point processes

For this section, we place ourselves in the case where E is some Euclidean space (our particular focus in this work shall be E = R 2 ). A point process Φ is called stationary if its distribution is invariant to translation. That is, for any t ∈ E, any n ≥ 1, and any

B 1 , • • • , B n ∈ B c , the distribution of {Φ(B 1 ), • • • , Φ(B n )} is equal to that of {Φ(B 1 + t), • • • , Φ(B n + t)}, where B + t := {u + t, u ∈ B}. Another way to say this is, noting Φ = i δ x i , Φ is equal in distribution to i δ x i -t .
It can be convenient to consider stationary processes, when one observes a large phenomenon that appear similar in different positions, at smaller scale. For instance, the positions of trees in a forest can be regarded, when observed at a certain scale, as a stationary point process. Another example would be the positions of stars, galaxies, or other cosmological objects, projected on a sphere around the Earth. The theory of stationary point processes gives us tools to better describe their properties. Here, we review the principal definitions and results for such processes.

In order to study the properties of stationary point processes, it is useful do define a stationary framework directly on the probability space. Definition 2.5.1 (Shift operator). For any t ∈ E, and any measure µ on (E, B), the shift operator S t applied to µ is the measure S t µ on (E, B) defined by, for all B ∈ B:

S t µ(B) = µ(B + t), ( 2.14)
where

B + t := {u + t, u ∈ B}. If µ is a counting measure, say µ = i δ u i , then S t µ = i δ u i -t .
The shift operator can be extended to functions f with values in any space by defining S t f (u) := f (u + t).

Definition 2.5.2 (Flow on the probability space). Let (Ω, F) a measurable space. Let us assume that there exists a family {θ t : Ω → Ω} t∈E of measurable mappings satisfying the following conditions:

-For any t ∈ E, θ t is a bijection from Ω to Ω.

-For any t, s ∈ E, θ s • θ t = θ s+t .

-The mapping θ : (E, Ω) → Ω defined by θ(t, ) := θ t ( ) is B ⊗ F-measurable.

Such a family will be called a (measurable) flow on (Ω, F), and the space equipped with this flow will be noted (Ω, F, {θ t }). Definition 2.5.3 (Processes compatible with the flow). A point process Φ : (Ω, F, {θ t }) → (M, M) is said to be compatible with the flow if

∀t ∈ E, Φ • θ t = S t Φ.
(2.15)

Example 2.5.4. Canonical probability space with the flow. The space (Ω, F, {θ t }) = (M, M, {S t }), is the canonical space on which we define Φ(µ) = µ, ∀µ ∈ M, compatible with the flow.

Definition 2.5.5 (Stationary probability). Let (Ω, F, {θ t }) be a measurable space equipped with a flow. Let P be a probability measure on (Ω, F) invariant with respect to any element in the flow, that is

∀ t ∈ E, P • θ -1 t = P. (2.16)
(Ω, F, {θ t }, P) is then called a stationary framework, and P is called the stationary probability on (Ω, F, {θ t }).

It directly follows from the above definition that a point process Φ defined on a stationary framework (Ω, F, {θ t }, P) is stationary.

Example 2.5.6. The distribution of a homogeneous Poisson point process is invariant to any shift, and can equip the canonical probability space.

Palm probabilities in the stationary framework Let Φ be a point process defined on the stationary framework, compatible with the flow. A first observation is that its mean measure is proportional to the Lebesgue measure. The coefficient of proportionality is called the intensity of the process. Definition 2.5.7 (Campbell-Matthes' measure). Let C be a measure on (E × Ω, B × F) defined by, for all B ∈ B and all

A ∈ F C(B × A) := E E 1 B (u)1 A (θ u )Φ(du) (2.17)
Note that Campbell-Matthes' measure is an extension of the mean measure of Φ:

C(B × Ω) = M Φ (B). Moreover, for all A ∈ F, C(B × A) ≤ M Φ (B), which implies that C(•×A) admits a Radon-Nikodym derivative dC(•×A) dM Φ (•) (u) w.r.t. M Φ . Additionally, similarly to M Φ , C(• × A) is invariant to any shift in E,
which gives us that it is proportional to the Lebesgue measure on E (= R d ). Thus, dC(•×A) dM Φ (•) (u) is a constant, depending only on A and λ. Seeing it as a function of A, one can easily check that it constitutes a probability measure on (Ω, F), leading us to the following definition. Definition 2.5.8 (Palm probability). Let Φ be a point process defined on the stationary framework (Ω, F, {θ t }, P), compatible with the flow, with finite non-null intensity λ. The Palm probability of Φ is the unique measure P 0 on (Ω, F) defined by

P 0 (A) := 1 λ|B| E E 1 B (u)1 A (θ u )Φ(du) , ∀A ∈ F, (2.18)
for any B ∈ B on finite, non-null Lebesgue measure |B|. We shall denote by E 0 the expectation under P 0 .

Theorem 2.5.9 (Campbell-Little-Mecke-Matthes). Let Φ be a point process defined on the stationary framework (Ω, F, {θ t }, P), compatible with the flow, with finite non-null intensity λ. Denote P 0 the Palm probability of Φ. For any non-negative measurable function f on E × Ω,

E E f (u, θ u •)Φ(du) = λ E E 0 [f (u, •)]du. ( 2
E 0 [(Φ -δ 0 )(B)] = B κ(u)du.
The function κ is called the reduced second order density of Φ.

Geometric marks for stationary point processes

One can extend the notion of point process, by considering counting measures with values attached to each atom. Such values, possibly random, are usually called the marks of the atoms, and the associated random measure is called a marked point process. In this work, we shall take interest in a particular class of marks, described in what follows. In this section, we consider E = R2 .

Definition 2.6.1 (Geometric marks). Let Φ be a point process defined on the stationary framework (Ω, F, {θ t }, P), compatible with the flow, with finite non-null intensity λ. Let m : R 2 × M → C be a measurable function, satisfying, for any u, v ∈ R 2 , and any φ ∈ M,

m(u -v, S v φ) = m(u, φ).
We say that the function m is translation invariant, and shall be called a geometric mark of Φ. With the convention we consider, the associated marked point process 2 , noted

Φ = u∈Φ m(u, Φ)δ u ,
is stationary. Note that a point process without marks is the special case of m ≡ 1.

In what follows, we briefly recall a few basic examples of geometric marks, that will be considered in this work (more particularly in Sections 3.3 and 3.4 and Chapter 4).

Example 2.6.2 (Shot-noise). For u ∈ R 2 and Φ = i δ u i ∈ M, we can define

m(u, Φ) := i 1(u = u i ) (|u -u i |),
a shot-noise functional, with some non-negative response function , that depends on the distance |u -u i | between u and u i ∈ φ, e.g. (r) = r β , for some β > 2.

Example 2.6.3 (Nearest-neighbour distance). We can define the nearest neighbour distance

m nn (u, Φ) = min{||u -u i ||, u i ∈ Φ}.
Example 2.6.4 (Voronoi cell area and moment of inertia). The Voronoi cell of an atom u ∈ Φ is defined by

V (u) = {v ∈ R 2 : ∀ u ∈ Φ, v -u ≤ v -u }.
One can consider the Voronoi cell volume

m V (u, Φ) = |V (u)|, (2.20)
or its moment of inertia

m I (u, Φ) := V |v -u| 2 dv.
Example 2.6.5 (Voronoi shot-noise). One can consider more sophisticated marks, such as

m Z (u, Φ) = i 1(u i = u) (|u -u i |)m V (u, Φ).

Ergodicity

Ergodic theory studies the conditions under which spatial averages of functions of point processes converge to the expectation of such functions, in order to bridge the gap between the theoretical point of view of probabilistic theory, and real-world measurements. We saw in Section 2.5 that a point process is stationary if it appears similar at different locations in the space. Ergodicity goes even further, saying that one can average a function over all different locations in the space, for a single realization, and by doing so, estimate the expected value of this function. In a sense, it can be viewed as a kind of Law of Large Numbers. This property is useful when we are dealing with a single, large observation, which can happen when collecting the data is difficult.

Recalling the notions of Section 2.5, let (Ω, F, {θ t } t∈E , P) a stationary framework (in this section, E is some Euclidean space). We begin by defining invariant events under this framework. Definition 2.7.1 (Invariant events). An event A ∈ F is said to be ({θ t } t∈E , P)-invariant if, for all t ∈ E, P(A∆θ t A) = 0, where ∆ is the symmetric difference: A∆B := (A ∪ B)\(A ∩ B). The space

I := {A ∈ F : A is ({θ t } t∈E , P)-invariant} is a σ-algebra, called the invariant σ-algebra.
The fundamental theorem in ergodic theory (in the context of point processes), which proof can be found in e.g. [START_REF] Kallenberg | Foundations of modern probability[END_REF], states the following. Theorem 2.7.2 (Birkhoff's individual ergodic theorem). Let (Ω, F, {θ t } t∈E , P) be a stationary framework, I its invariant σ-algebra. Let Φ be a point process on E, compatible withe the flow, with intensity 0 < λ < ∞. Let {B n } n∈N be a convex averaging sequence, i.e. an sequence of bounded, convex Borel sets such that for all n ∈ N, B n ⊂ B n+1 , and sup{r ≥ 0 : B n contains a ball of radius r} ---→ n→∞ ∞. Let f be a measurable, integrable function from M to C. Then

lim n→∞ 1 |B n | Bn f (θ t Φ) dt = E[f (Φ)|I], P-a.s.
We now need the following notions:

Definition 2.7.3 (Metrically transitive, ergodic, or mixing framework). We say that the stationary framework is metrically transitive if I is trivial, i.e. for all A ∈ F,

P(A) ∈ {0, 1}. The framework is ergodic if for all A 1 , A 2 ∈ F, lim a→∞ 1 (2a) d [-a,a] d P(A 1 ∩ θ t A 2 )dt = P(A 1 )P(A 2 ),
where d is the dimension of the Euclidean space E. The framework is said to be mixing if for all A 1 , A 2 ∈ F,

lim |t|→∞ P(A 1 ∩ θ t A 2 ) = P(A 1 )P(A 2 ).
The following result, which proof can be found in e.g. [START_REF] Baccelli | Random measures, point processes, and stochastic geometry[END_REF], states the relation between the three definitions from above. Proposition 2.7.4. Let (Ω, F, {θ t } t∈E , P) be a stationary framework. Then, the following relations are true mixing ⇒ ergodic ⇔ metrically transitive.

The mixing and ergodic properties are useful because they are easier to verify in practice than the metrical transitivity. These notions allow us to state the main result: Corollary 2.7.5. Under the conditions of Theorem 2.7.2, and if the stationary framework is metrically transitive, then

E[f (Φ)|I] = E[f (Φ)].
(2.21)

Linear operators and second-order properties for stationary point processes

Second order characteristics offer a natural way to describe geometric structures in point processes, by studying correlations between pairs of points. In this section, we shall place ourselves in E = R 2 . There have been extensively studied in the literature (see e.g. [START_REF] Ripley | The second-order analysis of stationary point processes[END_REF][START_REF] Bartlett | The spectral analysis of point processes[END_REF][START_REF] Brémaud | Fourier analysis of stochastic processes[END_REF]) and constitute some of the most widely used tools to analyse point processes.

We call a shot noise of a point process Φ any random variable defined by Φ(f ) =

Bartlett spectrum of a point process

The Bartlett spectrum of a stationary point process is the measure that allows one to express the covariance between shot-noises as a function of their Fourier transform. More formally, Proposition 2.8.1. Let Φ be a point process defined on the stationary framework (Ω, F, {θ t }, P), compatible with the flow, with finite non-null intensity λ. Suppose that γ(u)

:= λκ(u) -λ 2 ∈ L 2 (R 2 )
, where κ is the reduced second-order density of Φ (cf. Definition 2.5.11). Let us denote L 2 Φ the set of complex valued, measurable, integrable functions such that

E[Φ(|f |) 2 ] < ∞. For any two functions f, g ∈ L 2 Φ such that f ǧ ∈ L 2 (R 2 ), Cov(Φ(f ), Φ(g)) = 1 (2π) 2 R 2 f (ω)ĝ * (ω)b Φ (ω)(dω), where b Φ (ω) := γ * (ω) + λ. Proof. First remark that L 2 Φ ⊆ L 1 (R 2 , M Φ )∩L 2 (R 2 , M Φ ), so Cov(Φ(f ), Φ(g)) := E[Φ(f )Φ(g * )]- E[Φ(f )]E[Φ(g * )]
is well defined, and both terms in the r.h.s. of the equality are finite. Then, we decompose

Cov(Φ(f ), Φ(g)) = R 2 ×R 2 f (u)g * (v)M Φ (2) (d(u, v)) + λ R 2 f (u)g * (u)du -λ 2 R 2 ×R 2 f (u)g * (v)dudv.
This gives us, with the change of variables z = v -u,

Cov(Φ(f ), Φ(g)) = λ (R 2 ) 2 f (u)ǧ * (-u -z)(κ(z) -λ)dudz + λ (R 2 ) 2 f (u)ǧ * (-u -z)duδ 0 (dz) = λ R 2 f ǧ * (-u)(κ(u) -λ)du + λ R 2 f ǧ * (-u)δ 0 (du) = R 2 f ǧ * (u)Γ Φ (du), where Γ Φ (du) := λκ(u)du -λ 2 du + λδ 0 (du) (2.22) is called the covariance measure of Φ. Recall that γ(u) = λκ(u) -λ 2 ∈ L 2 (R 2 )
, with κ being symmetric, so by Plancherel's equality,

Cov(Φ(f ), Φ(g)) = R 2 f ǧ * (u)γ(u)du + λ R 2 f ǧ * (-u)δ 0 (du) = 1 (2π) 2 R 2 f (ω)ĝ * (ω)γ * (ω)dω + λ (2π) 2 R 2 f (ω)ĝ * (ω)dω = 1 (2π) 2 R 2 f (ω)ĝ * (ω)b Φ (ω)(dω). ( 2

Discrete Fourier Transform

Let Φ be a stationary point process, satisfying the conditions of the previous section.

Let s > 0 be some positive real value, and

W s := [-s/2,
E[U m,s (Φ)] = R 2 G s ( m s -ω)b Φ (ω)dω = G s b Φ ( m s ), with 3 G s (u) := 1 (2π) 2 s 2 sinc 2 ( s 2 u 1 ) sinc 2 ( s 2 u 2 )
, and sinc(u

) := sin(u) u . Furthermore, if b Φ is continuous at ξ ∈ E, then, if m → ∞, s → ∞, such that m s = ξ = 0, E[U m,s (Φ)] → b Φ (ξ).
Proof. We write, using Proposition 2.8.1, and the fact that E[U m,s ] = 0 for m = 0,

E[U m,s (Φ)] = 1 |W s | Var(F * m (Φ)) = 1 (2πs) 2 R 2 | fm | 2 (ω)b Φ (ω)dω,
with f m being the product of a complex exponential and the indicator function of the window W s , so fm (ω) = s 2 sinc( s

2 (ω 1 -m 1 s )) sinc( s 2 (ω 2 -m 2 s ))
, where ω = (ω 1 , ω 2 ), and m = (m 1 , m 2 ). As the sinc function is pair, we obtain

E[U m,s (Φ)] = R 2 G s ( m s -ω)b Φ (ω)dω = G s b Φ ( m s ).
that is, E[U m,s (Φ)] is the convolution of the Bartlett spectral measure with the function G s , evaluated at m s . 3 Gs is the square modulus of the Fourier transform of 1W s .

Suppose now that b

Φ is continuous at ξ ∈ E. Recall that E G s (ξ -ω)dω = 1. Then, let V be a neighbourhood of ξ such that, for any ω ∈ V , |b Φ (ω) -b Φ (ξ)| ≤ . We can write |G s b Φ (ξ) -b Φ (ξ)| ≤ | V G s (ξ -ω)(b Φ (ω) -b Φ (ξ))dω| + | V c G s (ξ -ω)(b Φ (ω) -b Φ (ξ))dω| ≤ V G s (ξ -ω)dω + | V c G s (ξ -ω)(b Φ (ω) -b Φ (ξ))dω|, with V G s (ξ -ω)dω ≤ 1, and | V c G s (ξ -ω)(b Φ (ω) -b Φ (ξ))dω| ≤ | V c G s (ξ -ω)b Φ (ω)dω| + |b Φ (ξ) V c G s (ξ -ω)dω|.
(2.26) Moreover, for any s ≥ 1

G s (ξ -ω)1 V c (ω) ≤ 1 s 2 (ξ 1 -ω 1 ) 2 (ξ 2 -ω 2 ) 2 1 V c (ω) ≤ 1 (ξ 1 -ω 1 ) 2 (ξ 2 -ω 2 ) 2 1 V c (ω) ∈ L 1 (R 2 ) ∩ L 2 (R 2 ),
and G s (ξ -ω) ---→ s→∞ 0. Therefore, the dominated convergence theorem applies.

Example 2.8.4 (Cox-circles point process). The Bartlett spectrum can be calculated for the Cox point process with Poisson points of linear intensity λ on circles in R 2 of radius R and the (total) intensity of points denoted by λ (cf. [START_REF] Chiu | Stochastic geometry and its applications[END_REF]Example 10.6]).

κ BC (u) = λ +      2λ π|u| R √ 4R 2 -|u| 2 for |u| ≤ 2R 0 otherwise, yielding b BC (ω) = λ(1 + 2πRλ J 2 0 (R|ω|)), (2.27)
where J 0 is the 0-th Bessel function of the first kind. Figure 2.2 shows the theoretical Bartlett spectrum, as well as its DFT estimator, averaged over 10 realizations in a window of size 1, for frequencies (averaged in angles) such that |ω| ∈ 2πZ 2 .

Spatial statistics and summary characteristics

In spatial statistics, one is dealing with the observation of one or several collections of points (or more generally signed measures), lying on a bounded window. These collections of points, called point patterns, or configurations, are assumed to be realizations of a point process (or a random signed measure). In order to perform a statistical analysis of the point patterns, it is often necessary to extract important information about this data by representing it in a relatively concise way. In other words, one usually define several functionals of the point patterns, that shall be considered as statistical estimators of moments of the underlying point process, that aim to characterize (fully or partially) its distribution (in Section 2.10, we saw how such estimators could be used to define a model for the underlying distribution). For example, when observing a point pattern φ in a window W , a very simple functional that can be considered is the average number of points per volume unit in φ, i.e. f (φ) = 1 |W | φ(W ), where |W | denotes the volume of W . Here, if φ is assumed to be a realization of a stationary point process Φ, then f (φ) (strictly speaking f (Φ)) can be viewed as an estimator of the intensity of Φ. This value alone obviously doesn't fully characterize the process Φ, and there exist many other classic functionals in the literature on point processes, called summary characteristics, that capture important geometric information about point patterns. In this section, we shall briefly review several of them, which are the most broadly used in the literature. While these functions can be defined for point processes on R d , for any d ≥ 1, we shall focus on the planar case (d = 2).

Second-order statistics

A natural way to extract information beyond the mere intensity of a point process is to look at second-order moments, that capture the pairwise interactions between the atoms of the process. To this end, on can consider Ripley's K-function ( [START_REF] Ripley | The second-order analysis of stationary point processes[END_REF]), defined for a stationary process Φ with intensity 0 < λ < ∞, as

K Φ (r) := 1 λ E 0 [Φ(B 0 (r)) -1], r ≥ 0,
which counts the average number of points present in a ball of radius r around a typical point of Φ. Alternatives to Ripley's K-function are Besag's L-function ( [START_REF] Besag | Contribution to the discussion on dr ripley's paper[END_REF]), defined as L Φ (r) := K Φ (r) π , and the pair correlation function g Φ (r) :=

K Φ (r)
2πr . These three functions offer different representations of the same information about the process. For more details about these functions and their comparisons, see for example [START_REF] Illian | Statistical analysis and modelling of spatial point patterns[END_REF], which also presents many other summary characteristics, among them the ones we shall present in the next section. Note that the Discrete Fourier Transform, defined in Section 2.8.2, also belongs to the family of second-order statistics.

Empty set and nearest neighbour statistics

Other statistics can be defined, relying on the size of sets that do not contain any points, or on the distance between neighbouring points. The two most commonly found statistics are the spherical contact distribution function (d.f.), defined as

H s (r) := P(Φ(B 0 (r)) > 0), r ≥ 0, (2.28)
and the nearest neighbour distance d.f.,

D(r) := P 0 (Φ(B 0 (r)) > 1), r ≥ 0. (2.29)
Both the spherical contact d.f. and the nearest neighbour d.f. measure the probability that at least one point is encountered in the ball of radius r, either centered at a typical point of the plane (for H s ) or at a typical point of Φ (for D). The nearest neighbour distance d.f. can be extended to an arbitrary number of neighbours:

D k (r) := P 0 (Φ(B 0 (r)) > k), r ≥ 0, k ≥ 1.
(2.30)

Remark. Note that these functions are defined by mathematical expectations, and are therefore theoretical. We find in the literature many estimators for these statistics, see e.g. [START_REF] Illian | Statistical analysis and modelling of spatial point patterns[END_REF].

Maximum entropy models

In Section 2.3.2, we presented a few classical distributions of point processes. When dealing with observed data, it is often useful to try to fit a model (usually a well known distribution) to this data. However, one can find itself in the situation where classic parametric models fail to represent the data in a satisfying way. Suggested by statistical physics, maximum entropy models have been introduced in [START_REF] Jaynes | Information theory and statistical mechanics[END_REF] to build models of probability distributions, based on partial knowledge of it, for instance by estimating moments of functionals from the observed data (in the case of point processes, one can use statistics called summary characteristics, that we present in Section 2.9). Intuitively, this means that the model is 'as random as possible' under some prescribed constraints. Two standard models are the macro-canonical model, with expectations constraints, and the micro-canonical model, with path-wise constraints. We review these two models in the following sections.

Macro-canonical model

Let us fix 0 < λ < ∞, and denote L 0 a Poisson distribution on (E, B). We define

M 0 := {L : L L 0 },
where ' ' means 'absolutely continuous with respect to'.

The notion of entropy is well defined for discrete random variables. It can be extended to stochastic processes on Euclidean spaces by the notion of differential entropy, which is the Kullback-Leibler divergence ( [START_REF] Kullback | On information and sufficiency[END_REF]) with respect to the Lebesgue measure. In this work, we shall use the following definition of entropy for point processes. and the inequality becomes an equality if and only if

L 1 = L 2 .
Proof. It follows from Jensen's inequality. By definition,

1 = E L 1 [1] = M dL 1 = M dL 1 dL 2 dL 2 = E L 2 [ dL 1 dL 2 ].
The function

f (u) = u log(u) being convex, 0 = f (E L 2 [ dL 1 dL 2 ]) ≤ E L 2 [f ( dL 1 dL 2 )],
with

E L 2 [f ( dL 1 dL 2 )] = E L 2 [log( dL 1 dL 2 ) dL 1 dL 2 ] = M log( dL 1 dL 2 ) dL 1 dL 2 dL 2 = M log( dL 1 dL 2 )dL 1 = KL(L 1 L 2 ).
Since the function f is strictly convex, we have an equality iff dL 1 dL 2 = 1 L 2 -almost everywhere.

Let K : M → R N , for some N ≥ 1, be a vector of statistics, and a ∈ R N . Note

A := {L ∈ M 0 : E L [K] = a}.
We are interested in solutions, if any exists, to the following problem:

L * ∈ arg max L∈A H(L).
(2.32)

We are going to show the following result:

Proposition 2.10.3. Suppose there exists some γ ∈ R N such that 0 < Z γ := E L 0 [e γ.(K-a) ] < ∞, and define L γ such that

dL γ dL 0 = 1 Z γ e γ.(K-a) . If γ.(K -a) is integrable w.r.t L γ
, and E Lγ [K -a] = 0, then L γ is the unique solution of eq. (2.32).

Proof. Let L ∈ A. First, note that L L γ . Indeed, for any Γ ∈ M, L γ (Γ) =

1

Zγ Γ e γ.(K-a) dL 0 , with e γ.(K(µ)-a) > 0 for all µ ∈ M. Therefore, if L γ (Γ) = 0, then L 0 (Γ) = 0 which implies L(Γ) = 0. Then,

-H(L) = M log( dL dL 0 )dL = M log( dL dL γ dL γ dL 0 )dL = M log( dL dL γ )dL + M log( dL γ dL 0 )dL = KL(L L γ ) + log( 1 Z γ ) M dL + γ. M (K -a)dL = KL(L L γ ) -log(Z γ ) ≥ -log(Z γ ), (2.33)
with the last inequality given by Lemma 2.10.2 and the fact that L L γ . Remark that

KL(L γ , L 0 ) = M log( dL γ dL 0 )dL γ = M log( 1 Z γ e γ.(K-a) )dL γ = -log(Z γ ) + γ.E Lγ [K -a]. ( 2
[K -a] = 0 that H(L) ≤ H(L γ ), with equality iff KL(L L γ ) = 0, that is, L = L γ .
Under some conditions, we can further express the solution of eq. (2.32) as the minimum of Z γ . Proposition 2.10.4. Suppose that, for all γ ∈ R * N , E L 0 [e γ.(K-a) ] < ∞, and

P L 0 (γ.(K - a) = 0) > 0. Then ∃ ! γ * ∈ arg min γ∈R N Z γ .
Moreover, suppose that there exists a neighbourhood V of γ * in which γ → Z γ is differentiable, and

∇Z γ = E L 0 [(K -a)e γ.(K-a) ] < ∞. Then, γ * ∈ arg min γ∈R N Z γ ⇔ L γ * ∈ arg max L∈A H(L). Proof. If for all γ ∈ R N , E L 0 [e γ.(K-a) ] < ∞, then γ → Z γ is defined on R N . It is easy to see that this function is convex. If ∀ γ ∈ R N , P L 0 (γ.(K -a) = 0) > 0 (implied e.g.
by Z γ = 1), the function is strictly convex, and the uniqueness of the minimum follows directly. If the function is differentiable on a neighbourhood of γ * , then γ * is the unique solution of ∇Z γ = 0. Moreover, as

∇Z γ = E L 0 [(K -a)e γ.(K-a) ], ∇Z γ = 0 is equivalent to 0 = E L 0 [(K -a)e γ * .(K-a) ] = Z γ * E L γ * [K -a].
We know that Z γ * > 0, which gives us the result.

Micro-canonical model

The micro-canonical model is defined similarly to the macro-canonical model, but with path-wise constraints rather than constrains in expectations. More formally, let us fix some > 0, and define A := {µ ∈ M : K(µ) -a ≤ }. Suppose that P L 0 (A) > 0, and define

M A := {L L 0 : P L (A) = 1}.
We are interested in solutions, if any exists, to the following problem:

L * ∈ arg max L∈M A H(L). ( 2
dL * dL 0 = 1 A P L 0 (A) . Proof. It is direct to show that L * ∈ M A . Note that KL(L * ||L 0 ) = A log( 1 A P L 0 (A) )dL * = A log( 1 A P L 0 (A) ) 1 A P L 0 (A) dL 0 = -log(P L 0 (A)) 1
P L 0 (A) A dL 0 = -log(P L 0 (A)). Now, let L ∈ A. For any Γ ∈ M, P L * (Γ) = 0 ⇒ P L 0 (Γ ∩ A) = 0 ⇒ P L (Γ ∩ A) = 0 ⇒ P L (Γ) = 0, so L L * .
One can then easily verify that

dL dL 0 = dL dL * 1 A P L 0 (A)
, L-a.s.

Then,

KL(L||L 0 ) = A log( dL dL 0 )dL = A log( dL dL * 1 A P L 0 (A) )dL = KL(L||L * ) + A log(1 A )dL -log(P L 0 (A)) A dL = KL(L||L * ) -log(P L 0 (A)) = KL(L||L * ) + KL(L * ||L 0 ) ≥ KL(L * ||L 0 ).
Maximum entropy models are defined by the description of the distribution, through the function K. This function, or descriptor, can be chosen, for instance, to be the summary characteristics mentioned in Section 2.9. In the next chapter, we shall study another class of descriptors, defined by a multi-scale representation of the geometry, using wavelets.

Chapter 3

Wavelet-based representations

As previously described, the aim of summary characteristics in spatial statistics is to concisely describe the main geometric features of a point pattern, for instance its tendency to exhibit clusters of points, or to form more regular patterns. Applied to realizations of a point process, they allow one to infer information about its distribution, by estimating their expected values. They can, for instance, be used to regress values associated to a point pattern, or to define a maximum entropy model.

Each statistic presented in Section 2.9 describe a particular aspect of a point process, and none of them fully characterizes its distribution. For instance, it is known that second order statistics do not characterize the distribution of a point process. For simple point processes, the void probability of all Borel sets characterize its distribution, but considering only balls may not be enough, especially for anisotropic distributions. Depending on the data and the task at hand, it is therefore a challenge to find the appropriate representation.

In this work, we present a class of representations of stationary point processes, constructed from a family of shot-noise processes, i.e. from a family

{ ψ ξ (u -v)Φ(du)} ξ∈Ξ ,
where the ψ ξ are measurable, integrable functions, and Ξ shall be specified later. This family, that we shall present in the next chapter, is designed to separate different relevant information about the geometric structures of processes. More precisely, these functions are wavelets, which we shall introduce in the following section.

Wavelet transform

Informally speaking, a wavelet ψ is a function averaging to zero, that is localized both is the spatial and the frequency domain (meaning that both ψ and ψ have most of their energy concentrated in compact domains). A wavelet ψ can be constructed in such a way that, by products of dilations and translations of this wavelet, one obtains an orthonormal basis of L 2 (R 2 ) (see e.g. [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF]). Wavelets have been extensively studied in the literature (e.g. [START_REF] Gabor | Theory of communication. part 1: The analysis of information[END_REF][START_REF] Calderón | Intermediate spaces and interpolation, the complex method[END_REF][START_REF] Grossmann | Decomposition of hardy functions into square integrable wavelets of constant shape[END_REF][START_REF] Daubechies | Ten lectures on wavelets[END_REF][START_REF] Daubechies | The wavelet transform, time-frequency localization and signal analysis[END_REF]), as they constitute a popular tool in many areas of science (see e.g. [START_REF] Goupillaud | Cycle-octave and related transforms in seismic signal analysis[END_REF][START_REF] Farge | Wavelet transforms and their applications to turbulence[END_REF][START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF][START_REF] Schneider | Wavelet methods in computational fluid dynamics[END_REF]). In this section, we briefly recall the basic definitions of wavelets, and the wavelet transform, in the context of our work.

While the definition of wavelets differ depending on the context, we shall consider the following:

Definition 3.1.1 (Wavelet family). Let ψ : R 2 → C be a complex valued function such that ψ ∈ L 2 (R 2 ) ∩ L 1 (R 2 )
, such that R 2 ψ(u)du = 0, and ψ 1 = 1. For any j ∈ Z, and any θ ∈ [0, 2π], define

ψ j,θ (u) := 2 -2j ψ(2 -j r θ u). (3.1)
If the function ψ satisfies the following condition

R 2 | ψ| 2 (ω) |ω| 2 dω < ∞, (3.2)
then ψ is called a wavelet. For any discrete subset Θ ⊂ [0, 2π], the family of functions {ψ j,θ } j∈Z,θ∈Θ is called a wavelet family. The condition eq. (3.2) is called the admissibility condition. A wavelet family also includes a low-pass filer (i.e. local averaging function), that we shall denote ψ 0 . It can be chosen, for instance, to be an isotropic Gaussian function.

Usually, one defines ψ such that ψ is centered at a frequency ξ 0 ∈ R 2 , and has fast decay in both the spatial and frequency domains. Informally, we can consider that there exists constants C, C , and such that |ψ|(u

) ≤ if |u| > C, and | ψ|(ω) ≤ if |ω -ξ 0 | > C .
A wavelet family is obtained by dilation and rotation of the wavelet ψ (often called the mother wavelet). Let θ ∈ [0, 2π], r θ the rotation by angle θ and j ∈ Z. The wavelet at scale j and angle θ can also be indexed by its central frequency

ξ := 2 -j r -θ ξ 0 , that is ψ ξ (u) = 2 -2j ψ(2 -j r θ u) ⇒ ψξ (ω) = ψ(2 j r θ ω).
Since ψ(ω) is centered around ξ 0 , it results that ψξ (ω) is centered around the frequency ξ. The wavelet ψ ξ at scale j has negligible amplitude for |u| > 2 j C, and its Fourier transform ψξ has negligible amplitude for ω such that |ω -ξ| > 2 -j C .

By indexing the wavelets by their central frequency, we shall sometimes use the notation {ψ ξ } ξ∈Ξ for the wavelet family, where Ξ denotes the set of frequencies corresponding to (j, θ) ∈ Z × Θ. [START_REF] Goupillaud | Cycle-octave and related transforms in seismic signal analysis[END_REF]), defined as

Example 3.1.2 (Morlet wavelet). A classic example is the Morlet wavelet (see

ψ(u) = 1 2πσ 2 e -|u| 2 2σ 2 (e iξ 0 .u -c ψ ), (3.3)
where c ψ is a constant such that R 2 ψ(u)du = 0. Example 3.1.3 (Bump steerable wavelet). Another example of wavelet that we shall consider in this work is the bump steerable wavelet, introduced in [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF]. It is defined in the frequency domain using polar coordinates. Let L be the even number of angles, and denote its central frequency by ξ 0 = (ω 0 , 0). Then

ψ(ω) = c ψ e - (|ω|-ω 0 ) 2 ω 2 0 -(|ω|-ω 0 ) 2 1 [0,2ω 0 ] (|ω|) cos L 2 -1 (ϕ(ω))1 ϕ(ω)< π 2 , (3.4)
where c ψ is a normalizing constant. Figure 3.3 shows the real and imaginary parts of a bump steerable wavelet, for a given ξ 0 , and Figure 3.4 displays a bump steerable wavelet family in the spatial and frequency domains.

∈ R 2 , ψ(t -•) ∈ L 1 (R 2 , µ).
The wavelet transform of µ is defined by the convolutions

Wµ := {µ ψ 0 , µ ψ j,θ } j∈Z,θ∈[0,2π] , ( 3.5)
with

µ ψ j,θ (t) := R 2 ψ j,θ (t -u)µ(du), ( 3.6)
and similarly for ψ 0 . The family of functions Wµ defined by the wavelet transform of µ are called the wavelet coefficients (at scale j and orientation θ). In this work, for simplicity, we shall only consider wavelets verifying

ψ(u) = ψ * (-u), for u ∈ R 2 .
For a random measure Φ, the wavelet transform is well defined if all ψ ∈ Ψ are integrable with respect to the mean measure M Φ of Φ. If Φ is stationary with finite non-null intensity, then the wavelet transform is always well defined, due to Campbell's averaging formula.

Note that, if Φ is a.s. finite, then for any

ω ∈ R 2 , Φ ψ j,θ (ω) = ψj,θ (ω) Φ(ω),
where

Φ(ω) = 1 2π R 2 e -2iπt
.ω Φ(dt). This highlights the fact that, if ψ j,θ is localized in frequency, then so is Φ ψ j,θ . In other words, the wavelet transform separates the information in Φ into components at different scales and orientations, or equivalently, in separate domains in the frequency plane.

Scattering moments

The wavelet transform is a linear operator that represents a measure by a family of measurable functions on R 2 . When performing statistical analysis or learning tasks, one is usually looking for non-linear representations of the input data. Moreover, if Φ is a stationary process (with non-null finite intensity λ), taking the expectation of the wavelet coefficients would give, for any

u ∈ R 2 , E[Φ ψ j,θ (u)] = λ R 2 ψ j,θ (t)dt = 0.
This shows that, in order to extract meaningful information about the geometry of the process, one needs to apply a non-linear operator to the wavelet transform of Φ. This can be done, for instance, by applying the complex modulus operator to its wavelet coefficients. Then, taking the expectancy of the resulting functions gives us a description of the point process distribution. This defines the (first order) scattering transform of Φ, inspired by [START_REF] Mallat | Group invariant scattering[END_REF]: Definition 3.2.1 (First order scattering moments). Let Φ be a stationary point process (or random measure) on (R 2 , B), and Ψ = {ψ 0 , ψ j,θ } j∈Z,θ∈[0,2π] be a wavelet family such all wavelet coefficients are almost surely integrable stochastic processes. The first order scattering moments of Φ, noted S 1 Φ, are defined, for all j ∈ Z, θ ∈ [0, 2π], by

S 1 Φ(j, θ) := E[|Φ ψ j,θ |], (3.7)
which does not depend on u ∈ R 2 by stationarity of Φ. Intuitively, the first order scattering moments compute 'how much geometric structure at scale j and orientation θ is present in Φ'. For example, if Φ exhibits a certain regularity at large scales (that is, from a certain scale, the process appears uniform, with no irregularities), then the larges scales scattering moments will have small values. Conversely, if Φ is composed of clusters of atoms, well separated from each other, then the values of the scattering moments at this scale will have large value.

If S 1 Φ(j, θ) detects the average regularity or clustering behaviour of Φ at a certain scale and orientation, on can enquire about the regularity or clustering of S 1 Φ(j, θ) at larger scales (i.e. the regularity/clustering of the regularity/clustering). This can be done by computing the scattering transform of S 1 Φ(j, θ). It defines the second order scattering moments of Φ. Definition 3.2.2 (Second order scattering moments). Let Φ be a stationary point process on (R 2 , B). Let us denote

I := {(j 1 , j 2 , θ 1 , θ 2 ) ∈ Z 2 × [0, 2π] 2 : j 2 > j 1 }. If, for all (j 1 , j 2 , θ 1 , θ 2 ) ∈ I, the random process |Φ ψ j 1 ,θ 1 | ψ j 2 ,
θ 2 is integrable, then Φ admits second order scattering moments, defined by

S 2 Φ(j 1 , θ 1 , j 2 , θ 2 ) := E ||Φ ψ j 1 ,θ 1 | ψ j 2 ,θ 2 | , (3.8) for (j 1 , j 2 , θ 1 , θ 2 ) ∈ I.
Note that we only select indices j 2 in I such that j 2 > j 1 . This is heuristically explained by the fact that the modulus operator in the first order scattering transform 'kills' the high frequency modulations of the wavelet coefficients Φ ψ j 1 ,θ 1 . Therefore, the Fourier transform of the |Φ ψ j 1 ,θ 1 | are roughly supported on the low frequencies, which implies that the coefficients |Φ ψ j 1 ,θ 1 | ψ j 2 ,θ 2 are negligible when j 2 ≤ j 1 .

Remark. One can define high order scattering moments, by cascading convolutions with the wavelet family, composed with a modulus operator.

Additionally, the definition of the wavelet transform in Section 3.1 can be extended to stationary marked point processes, by defining, for some wavelet ψ, and

any v ∈ R 2 , Φ ψ(v) : = R 2 ×R ψ(v -u)m(u, Φ)Φ(du) (3.9)
The wavelet transform of Section 3.1 is a particular case, where m ≡ 1. Similarly, under integrability conditions, one can extend the definition of first and second order scattering moments (Section 3.2) to (geometrically) marked point processes, and they do not depend on u ∈ R 2 .

Asymptotic properties of scattering moments

To gain more insights about what information this representation captures, we shall study the limit of first order scattering moments when the scale parameter j goes to ±∞.

For simplicity, in this section we shall assume that the wavelets have compact support, included in some bounded set B ∈ B, and are symmetric, i.e. ψ(-u) = ψ(u) 1 . The results are true for any θ ∈ [0, 2π], and for notation simplification, we omit the subscript θ in the following section.

We consider here random, possibly signed measures of the form (cf. Section 2.6)

Φ = u∈Φ m(u, Φ)δ u ,
where Φ is a stationary point process, with intensity 0 < λ < ∞, and the functions m are stationary, in the sense that, for any

u, v ∈ R 2 , m(u -v, S v Φ) = m(u, Φ).
Typically, m, called a marking function, is a measurable function of the geometry of the points in Φ, such as the ones presented in what follows.

Small scale limit

The small scale limit behaviour (j → -∞) of first order scattering moments can be described by the following result.

m such that B E 0,u [|m(0, Φ)|]|κ(u)|du < ∞. Then, S 1 Φ(j) = λE 0 [|m(0, Φ)|] + O(λ2 2j+1 ψ ∞ B 2 E 0,2 j (v-u) [|m(0, Φ)|]κ(2 j (v -u))dudv).
Proof. We can write the following decomposition:

S 1 Φ(j) := E[| R 2 ψ j (u) Φ(du)|] = E[ R 2 |ψ j (u)|| Φ|(du)] + e j ,
where | Φ| := u∈Φ |m(u, Φ)|δ u , and

e j := S 1 Φ(j) -E[ R 2 |ψ j (u)|| Φ|(du)]. Remark also that m(u, Φ) = m(0, S u Φ) = m(0, Φ • θ u ).
Then, as Φ is stationary, we can apply the CLMM theorem (eq. (2.19)):

E[ R 2 |ψ j (u)|| Φ|(du)] = E[ R 2 |ψ j (u)m(u, Φ)|Φ(du)] = λ R 2 |ψ j (u)|E 0 [|m(0, Φ)|]du. (i) (3.10)
Furthermore, note that ψ j is supported in B j := {2 j u : u ∈ B}. Then,

|e j | = |E[| B j ψ j (u)m(u, Φ)Φ(du)| - B j |ψ j (u)m(u, Φ)|Φ(du)]| =: |E[Z]| ≤ E[|Z|] = E[|Z|1 {Φ(B j )>1} ] ≤ 2 E[ B j |ψ j (u)m(u, Φ)|Φ(du)1 {Φ(B j )>1} ] ≤ 2 ||ψ j || ∞ 2 -2j ψ ∞ E[ B j |m(u, Φ)|Φ(du)1 {Φ(B j )>1} ]. (ii)
Then we write:

E[ B j |m(u, Φ)|Φ(du)1 {Φ(B j )>1} ] ≤ E[ B j |m(u, Φ)|Φ(du)(Φ(B j ) -1)] = E[ B j ×B j |m(u, Φ)|(Φ -δ u )(dv)Φ(du)] = B j ×B j E u,v [|m(u, Φ)|]ρ (2) (u, v)dudv = λ B j ×B j E 0,v-u [|m(u, S u Φ)|κ(v -u)dudv = λ B j ×B j E 0,v-u [|m(0, Φ)|κ(v -u)dudv = λ2 4j B×B E 0,2 j (v-u) [|m(0, Φ)|]κ(2 j (v -u))dudv, (iii)
where the fourth line stems from the stationarity of Φ, and the fifth from the translation invariance property of m. Combining eqs. (i) to (iii) completes the proof.

Remark. In the case of a Poisson point process with f ≡ 1 we find the result from [START_REF] Bruna | Intermittent process analysis with scattering moments[END_REF]:

SΦ(j) = λ(||ψ|| 1 + O j→-∞ (2 2j λ))

Large scale limit

We consider here the case j → +∞. The results of this section rely on some notion of central limit theorem, because the convolution of a wavelet ψ with a marked point process Φ operates a sum over the atoms of Φ. Intuitively, as the scale of the wavelet increases, the sum extends to a wider subset of R 2 . Therefore, in order to apply the CLT, we need some notion of near independence of the positions of the atoms and their marks, as the distance between them increases. Hence, we need to introduce several definitions before stating the main result, based on [START_REF] Błaszczyszyn | Limit theory for geometric statistics of point processes having fast decay of correlations[END_REF]. ), we say that the pair (m, Φ) satisfies the p-moment condition if sup

Definition 3.3.2. Let p ∈ [1, ∞
1≤p ≤ p sup u 1 ,••• ,u p ∈R d E u 1 ,••• ,u p [|m(u 1 , Φ)| p ] < ∞,
where sup signifies here ess sup with respect to M Φ (p) .

Definition 3.3.3. We call the correlation functions of Φ the functions defined for all p ∈ N, and all (k

1 , • • • , k p ) ∈ N p by: m (k 1 ,••• ,kp) (u 1 , • • • , u p ) := E u 1 ,••• ,up [m(u 1 , Φ) k 1 • • • m(u p , Φ) kp ] ρ (p) (u 1 , • • • , u p ) Definition 3.3.4.
The correlation functions of Φ are said to be ω-mixing if there exists a decreasing function ω :

N × R + → R + such that, ∀ n ∈ N, lim u →∞ ω(n, u) = 0 and, ∀ p, q ∈ N s.t. p < q, ∀ (u 1 , • • • , u q ) ∈ R q , ∀ (k 1 , • • • , k q ) ∈ N q , |m (k 1 ,••• ,kq) (u 1 , • • • , u q ) -m (k 1 ,••• ,kp) (u 1 , • • • , u p )m (k p+1 ,••• ,kq) (u p+1 , • • • , u q )| ≤ ω(K, d),
where K := q i=1 k i , and

d := inf i≤p, p<j≤q |u j -u i |.
Definition 3.3.5. Φ is said to have fast decay of correlations if its correlation functions are ω-mixing, with ω(n, u) = Cn f (c n u) for some rapidly decreasing function f : R + → [0, 1] and some constants cn , Cn ∈ R + .

For a stationary point process Φ on R 2 with intensity λ and a marking function m, we define

σ 2 (m, Φ) := λ E 0 [m 2 (0, Φ)] + R 2 (E 0,u [m(0, Φ)m(u, Φ)]κ(u) -E 0 [m(0, Φ)] 2 λ)du
provided the integral exists. We observe that, in the spirit of (2.22), σ 2 (m, Φ) ≥ 0 is the total mass of the covariance measure of Φ. For a wavelet function ψ = Real(ψ) + iIm(ψ), denote by N (ψ) the complex Gaussian random variable with mean 0 and covariance matrix (between real and imaginary parts)

Σ(ψ) := Real(ψ) 2 2 Real(ψ), Im(ψ) Real(ψ), Im(ψ) Im(ψ) 2 2
, where •, • denotes the dot product in L 2 (R 2 ).

Proposition 3.3.6. Let Φ be stationary point process on R 2 with intensity λ, and m a marking function, and ψ a compactly supported wavelet. Suppose that, for all p ∈ (1, ∞), the pair (m, Φ) satisfies the p-moment condition, Φ has fast decay of correlations. Then,

0 ≤ σ 2 (m, Φ) < ∞ and, if furthermore σ 2 (m, Φ) > 0, then lim j→∞ 2 j S 1 Φ(j) = σ(m, Φ)E[|N (ψ)|]. (3.11)
If the (complex) wavelet is such that its real and imaginary parts are orthogonal, with equal L 2 norm (e.g. for Morlet wavelets) then

E[|N (ψ)|] = ψ 2 √ π 2 . Proof. Let us denote Y j := ψ j Φ(0) = ψ j d Φ = 2 -2j u∈Φ∩2 j B ψ(2 -j u)m(u, Φ).
Considering the mother wavelet ψ as a test function (with support contained in a square window), from [22, Theorem 1.12 (ii)] we know that lim j→∞

2 2j Var(Real(Y j )) = σ 2 (m, Φ) Real(ψ) 2 2 (3.12)
and similarly for Var(Im(Y j )). Using CLMM (Theorem 2.5.9),

E[Y j ] = 2 -2j λE 0 [m(0, Φ)] R 2 ψ(2 -j u)du = 0.
Hence, by [START_REF] Błaszczyszyn | Limit theory for geometric statistics of point processes having fast decay of correlations[END_REF]Theorem 1.13 with Remark (vii) Multivariate central limit theorem] (since

ψ is a complex function) 2 j Y j d ---→ j→∞ σ(m, Φ)N (ψ).
We then conclude using the following lemma, from [START_REF] Dasgupta | Asymptotic theory of statistics and probability[END_REF], for which the conditions are satisfied, due to the p-moment condition:

Lemma 3.3.7. Let (X j ) j≥1 be a sequence of random variables, such that X j d ---→ j→∞ X.
If, for some δ > 0, sup

j E[|X j | 1+δ ] < ∞, then E[|X j |] ---→ j→∞ E[|X|].
For a specific wavelet ψ, such that Real(ψ), Im(ψ) = 0 and Real(ψ) 2 2 = Im(ψ) 2 2 , observe that real and imaginary parts of N (ψ) are independent Gaussians

ψ 2 √ 2 N (0, 1). Consequently |N (ψ)| ∼ ψ 2 √ 2 R(1), where R(1) is Rayleigh distribution (i.e. square root of χ 2 (2) -χ-distribution with parameter 2) giving E[|N (ψ)|] = ψ 2 √ 2 E[R(1)] = ψ 2 √ 2 π 2 = ψ 2 √ π 2 .
Remark. In [START_REF] Błaszczyszyn | Limit theory for geometric statistics of point processes having fast decay of correlations[END_REF], the authors give several detailed conditions concerning the point process Φ and marking function m under which Φ satisfies the fast decay of correlations and p-moment conditions, including σ 2 (m, Φ) > 0, (requested in Proposition 3.3.6). The processes satisfying these conditions include unmarked Poisson, Gaussian determinantal, and cluster Poisson processes with fixed range dependence (such as Matérn cluster processes, cf. Section 2.3.2), that we study in Example 3.3.8. These processes admit some local marking functions m, such as statistics of Voronoi cells, or k-nearest neighbors. Other, more regular distributions (such as the Ginibre point process) are shown in [START_REF] Błaszczyszyn | Limit theory for geometric statistics of point processes having fast decay of correlations[END_REF] to satisfy the fast decay of correlations and p-moment conditions, but with σ 2 (m, Φ) = 0. In this case, it is possible to obtain the Central Limit Theorem under some technical conditions, and we refer the reader to [START_REF] Błaszczyszyn | Limit theory for geometric statistics of point processes having fast decay of correlations[END_REF] for more details.

There also exist Central Limit Theorem results for Y i involving more clustering point process, not satisfying conditions of Proposition 3.3.6, such as intersections of Poisson line process (see [START_REF] Paroux | Quelques théoremes centraux limites pour les processus poissoniens de droites dans le plan[END_REF][START_REF] Heinrich | Central limit theorems for poisson hyperplane tessellations[END_REF]) .

More generally, observe that the asymptotic of the variance of Y j = ψ j d Φ depends on the wavelet ψ and the marked point Φ. In the case of m ≡ 1, using Proposition 2.8.1, we can express the variance of Y j as

2 2j Var(Y j ) = 1 (2π) 2 R 2 | ψ(ω)| 2 b Φ (2 -j ω)dω, (3.13)
where b Φ is the density of the Bartlett spectrum of Φ. If Φ is a Poisson point process, then b Φ (ω) = 1, confirming that 2 2j Var(Y j ) → ψ 2 2 b Φ (0), as shown in [START_REF] Bruna | Intermittent process analysis with scattering moments[END_REF]. If b Φ is continuous at 0, with 0 < b Φ (0) < ∞ (e.g. Gaussian determinantal point processes and Matérn cluster processes) then we obtain the same limit of the variance 2 2j Var(Y j ). If lim |ω|→0 b Φ (ω) = 0 or lim |ω|→0 b Φ (ω) = ∞, the asymptotic of S 1 Φ(j) may be inferior or superior to 2 -j respectively, which the numerical experiments of

scattering moments (averaged in angle) for Poisson point processes of different intensities, a

Gaussian determinantal point process, and a Matérn cluster process. We can observe that when j → -∞ (j = 0 is the smallest scale in our implementation), the scattering moments converge towards log 2 (λ) (since log 2 ( ψ 1 ) = 0), and when j → ∞, the slope goes to -1, independently of the distribution of Φ. However the value at the origin of this affine limit changes from one distribution to another. Let us recall that this value is:

a := 1 2 log 2 λ π 4 ψ 2 2 (1 + R 2 (κ(u) -λ)du) . Recall from eq. (2.22) that the covariance measure of Φ is defined as Γ Φ (du) = λδ 0 + λ(κ(u) -λ). Therefore a = 1 2 log 2 π 4 ψ 2 2 Γ Φ (R 2 )
measures the total mass of the covariance measure of Φ, indicating the general repulsion of points (Γ Φ (R 2 ) < λ) or their attraction (Γ Φ (R 2 ) > λ). In Figure 3.5, we observe that while the Gaussian determinantal, Poisson, and Matérn cluster point processes have the same intensity λ = 256, and the same slope of scattering moments when j goes to ∞, they exhibit different values of a, that are increasing according to the repulsion, Poisson-like independence, and clustering behaviour. For Poisson distributions,

R 2 (κ(u) -λ)du = 0, so a = 1 2 log 2 λ π 4 ψ 2 2 .
In our simulations, using a Morlet mother wavelet with central frequency (3π/4, 0), ψ 2 169, giving for a Poisson point process with intensity λ = 256, a = log 2 (256 π 4 169) 11.2. For j = 7, we find log 2 (S 1 Φ(j)) 4.2, so j + log 2 (S 1 Φ(j)) 11.2. Example 3.3.9 (Hyperuniform and hyper clustering processes). In Figure 3.6, we compare a Poisson distribution with different distributions that may not satisfy the assumptions of Proposition 3.3.6: a hyperuniform matching process (see [START_REF] Klatt | Hyperuniform and rigid stable matchings[END_REF]), a Ginibre point process (which we approximate using eigenvalues of Hermitian matrices with complex Gaussian entries, see e.g. [START_REF] Ginibre | Statistical ensembles of complex, quaternion, and real matrices[END_REF]), a Poisson lines crossing process (the points of the process are the intersections of random lines uniformly and independently drawn in the window), and a Poisson point process geometrically marked with the areas of its Voronoi cells (cf. Example 2.6.4).

Factorial moment expansion for scattering moments

In this section, we shall extend the result from [START_REF] Błaszczyszyn | Factorial moment expansion for stochastic systems[END_REF], expressing functionals of a point process Φ through its factorial moments measures. In the context of first order scattering moments of marked point process, we consider functions of the form

F (Φ) = | R 2 ψ(u)m(u, Φ)Φ(du)|.
The idea of the factorial moment expansion is to express the function F (Φ) as a telescoping sum depending on the number of atoms in Φ. In analogy to the Taylor expansion for functions of Euclidean spaces, on can think of this telescoping sum as a sum of derivatives, expressed at the origin. However in our case, the origin is the point measure with no atoms, noted ∅. For some marking functions, it is not well defined to apply them to null, or even finite measures (think of the Voronoi cell of an atom with no points around it). To remedy this problem, we shall express such functions in the form F (Φ) = f (Φ, Φ), where f is a function from M × M to R. This shall allow us to circumvent the telescoping for the value of the marking function.

For any u ∈ E, and µ ∈ M, we note µ u ∈ M the counting measure defined by:

∀B ∈ B, µ| u (B) = µ(B ∩ B 0 (u)),
where

B 0 (u) = {v ∈ R 2 : v < u }. Let f : M×M → R. We say that f is continuous at ∞ if ∀ν, µ ∈ M, and ∀(u n ) n ∈ (R 2 ) N s.t. u n ---→ n→∞ ∞, f (ν, µ| un ) ---→ n→∞ f (ν, µ).
We define the first order differential operator of f at u by: f (1) u :

M × M → R (ν, µ) → f (ν, µ| u + δ u ) -f (ν, µ| u ).
We recursively define the nth order differential operator by:

f (n) u 1 ,••• ,un (ν, µ) = f (n-1) u 1 ,••• ,u n-1 (1) un (ν, µ).
We remark that if f is continuous at ∞, then so is f [START_REF] Błaszczyszyn | Factorial moment expansion for stochastic systems[END_REF][START_REF] Błaszczyszyn | A note on expansion for functionals of spatial marked point processes[END_REF] for a proof of this fact).

(n) u 1 ,••• ,un , for all n ≥ 1 and (u 1 , • • • , u n ) ∈ (R 2 ) n (see e.g
Observe that, if

µ = n i=1 δ u i , with u 1 ≤ • • • ≤ u n , then, for any ν ∈ M, f (ν, µ) = f (ν, ∅) + n i=1 f (ν, µ| u i + δ u i ) -f (ν, µ| u i ) = f (ν, ∅) + R 2 (f (ν, µ| u + δ u ) -f (ν, µ| u )µ(du) = f (ν, ∅) + R 2 f (1) u (ν, µ)µ(du).
Now, let µ ∈ M, and f continuous at ∞ and such that E |f

(1)

u (µ, µ)|µ(du) < ∞. Then, for any (v n ) n ∈ E n such that v n ---→ n→∞ ∞, we can write f (ν, µ) = lim n→∞ f (ν, µ| vn ) = f (ν, ∅) + lim n→∞ R 2 f (1) u (ν, µ| vn )µ| vn (du) = f (ν, ∅) + lim n→∞ R 2 f (1) u (ν, µ)µ| vn (du) = f (ν, ∅) + R 2 f (1) u (ν, µ)µ(du),
where the last equality follows by dominated convergence.

Remark. As we can see, the measure ν is unchanged in the r.h.s. of the above equality. This shall allow one to compute the value of the marking function on the entire measure Φ.

Let Φ be a point process such that

R 2 E !u [|f (1)
u (Φ + δ u , Φ)|]M Φ (du) < ∞,
where M Φ denotes the mean measure of Φ. Then, using the Campbell-Little-Mecke theorem, we obtain

E[f (Φ, Φ)] = E[f (Φ, ∅)] + R 2 E u [f (1) u (Φ, Φ)]M Φ (du) = E[f (Φ, ∅)] + R 2 E !u [f (1) u (Φ + δ u , Φ)]M Φ (du), (3.14)
This leads us to the following result: Proposition 3.4.1. Let Φ be a simple point process, and let f be a measurable function, continuous at ∞, such that, for any

j ∈ {1, • • • , n + 1}, (R 2 ) j E !u [|f (j) u (Φ + δ u , Φ)|]M Φ (j) (du) < ∞. (3.15)
Then,

E[f (Φ, Φ)] =E[f (Φ, ∅)] + n j=1 (R 2 ) j E !u [f (j) u (Φ + δ u , ∅)]M Φ (j) (du) + (R 2 ) n+1 E !u [f (n+1) u (Φ + δ u , Φ)]M Φ (n+1) (du), (3.16)
with, for simplicity,

δ u denoting i δ u i . If furthermore, (R 2 ) n E !u [f (n) u (Φ + δ u , Φ)]M Φ (n) (du) ---→ n→∞ 0, (3.17) then E[f (Φ, Φ)] = E[f (Φ, ∅)] + ∞ j=1 (R 2 ) j E !u [f (j) u (Φ + δ u , ∅)]M Φ (j) (du). (3.18)
Proof. Equation (3.14) tells us that the result holds for n = 0. Suppose that it is true for some n ∈ N, and eq. (3.15) holds for n

+ 1. Let u ∈ (R 2 ) n , v ∈ R 2
, and let Φ u a point process with distribution

P ! u . u∈(R 2 ) n v∈R 2 E !v [|(f (n) u ) (1) v (Φ u + δ u + δ v , Φ u )|]M Φu (dv)M Φ (n) (du) = (u,v)∈(R 2 ) n+1 E !u,v [|f (n+1) (u,v) (Φ + δ u + δ v , Φ)|]M Φ (n+1) (d(u, v)) < ∞ Thus, E !v [|(f (n) u ) (1) v (Φ u + δ u + δ v , Φ u )|]M Φu (dv) < ∞, M Φ (n) -a.
s. We can therefore apply eq. (3.14) to f

(n) u (Φ u + δ u , Φ u )
, and obtain the result for n + 1.

We can now apply this result to first order scattering moments of marked point processes.

Corollary 3.4.2. Les Φ be a simple stationary point process, ψ j a wavelet, m(u, Φ) a marking function, and

f (ν, µ) = | R 2 ψ j (u)m(u, ν)µ(du)|,
such that the conditions of Proposition 3.4.1 are satisfied, including the condition given by eq. (3.17). Then,

S 1 Φ(j) = λE 0 [|m(0, Φ)|] + v < u E u,v [|ψ j (u)m(u, Φ) + ψ j (v)m(v, Φ)|]ρ (2) (u, v)dudv - v < u E u,v [|ψ j (u)m(u, Φ)| + |ψ j (v)m(v, Φ)|]ρ (2) (u, v)dudv + R, with R = ∞ n=3 1 n! J∈{2 n } (-1) n-|J| (R 2 ) n E u J [| i∈J ψ j (u i )m(u i , Φ)|]ρ (n) (u)du, noting u = (u i ) i=1,••• ,n , and u J = (u i ) i∈J .
Proof. The proof follows directly from Proposition 3.4.1, and the following Lemma ( [START_REF] Błaszczyszyn | Factorial moment expansion for stochastic systems[END_REF]).

Lemma 3.4.3. Let f be a measurable function, continuous at ∞. Let n ≥ 1, and

u = (u 1 , • • • , u n ) ∈ E n . Then, noting {2 n } the set of all subsets of {1, • • • , n}, for any ν, µ ∈ M, if u n < • • • < u 1 , f (n) u (ν, µ) = J∈{2 n } (-1) n-|J| f (ν, µ| un + i∈J δ u i ), (3.19) otherwise f (n) u (ν, µ) = 0.
u n < • • • < u 1 , f (n) u (ν, µ) = J∈{2 n }\{1} (-1) n-1-|J| f (ν, µ| un + δ u 1 + i∈J δ u i ) -(-1) n-1-|J| f (ν, µ| un + i∈J δ u i ), (3.20)
where {2 n }\{1} denotes all subsets of {2, • • • , n}. Therefore,

|f (n) u (ν, µ)| ≤ J∈{2 n }\{1} |f (ν, µ| un + δ u 1 + i∈J δ u i ) -f (ν, µ| un + i∈J δ u i )|,
which gives us,

|f (n) u (ν, µ)| ≤ J∈{2 n }\{1} || R 2 ψ j (u)|V (u)|(µ J + δ u 1 )(du)| -| R 2 ψ j (u)|V (u)|(µ J )(du)|| ≤ J∈{2 n }\{1} | R 2 ψ j (u)|V (u)|(µ J + δ u 1 -µ J )(du)| ≤ J∈{2 n }\{1} |ψ j (u 1 )||V (u 1 )| ≤ 2 n-1 |ψ j (u 1 )||V (u 1 )|, (3.21)
where we noted µ J = µ| un + i∈J δ u i . Applying eq. (3.21), we obtain

(R 2 ) p E !u [|f (p) u (Φ + δ u , Φ)|]M Φ (p) (du) ≤ 2 p-1 (p -1)! (R 2 ) p 1 |ψ j (u 1 )|E u [|V (u 1 )|]M Φ (p) (du), (3.22) where we note (R 2 ) p 1 = {(u 1 , • • • , u p ) : ∀ i = 1, u i < u 1 }.

Now, in order to bound the expected value of the Voronoi cell surface

E u [|V (u 1 )|],
we need to explicitly write its dependency on Φ:

E u [|V (u 1 )|] = E !u [|V Φ+δu (u 1 )|] (3.23) ≤ E !u [|V Φ+δu 1 (u 1 )|] (3.24) = R + P ! u (|V Φ+δu 1 (u 1 )| > t)dt = R + P ! u (|V Φ+δu 1 (u 1 )| > πr 2 )2πrdr ≤ R + P ! u (∃i ∈ {1, • • • , 6} : Φ(∠ i (r)) = 0)2πrdr (*) ≤ 6 i=1 R + P ! u (Φ(∠ i (r)) = 0))2πrdr, (3.25)
where

{∠ i (r)} i=1,••• ,6 forms a partition of B u 1 (r)
, with every {∠ i (r)} having the same volume. The explanation for eq. (*) is geometric: if every subset

{∠ i (r)} of B u 1 (r) contains an atom, then the Voronoi cell of u 1 is contained in B u 1 (r) 2 .
Let Φ be a Poisson point process on (R 2 , B), with intensity 0 < λ < ∞. We know that, for all B ∈ B, )] . Furthermore, we know that, for any N ≥ 1, and any u ∈ (R 2 ) N , the reduced Palm version of Φ at u, noted Φ ! u , is equal in distribution to Φ. Thus, for all B ∈ B,

P(Φ(B) = 0) = e -E[Φ(B
P(Φ ! u (B) = 0) = e -E[Φ(B)] .
Now, let Φ K be a stationary determinantal point process, with Hermitian kernel K, and intensity λ (that is, for all

B ∈ B, E[Φ K (B)] = E[Φ(B)]).
Then, for all B ∈ B, we know that ( [START_REF] Błaszczyszyn | On comparison of clustering properties of point processes[END_REF]),

P(Φ K (B) = 0) ≤ P(Φ(B) = 0) = e -E[Φ(B)] = e -E[Φ K (B)] . Let u = (u 1 , • • • , u p ) ∈ R p . One can show that the reduced Palm version of Φ K at u, noted Φ K! u
, is also a determinantal point process, which kernel shall be noted K u . Denote also Φ Ku! P the Poisson point process with the same mean measure as Φ K! u . Then, for all B ∈ B ( [START_REF] Błaszczyszyn | On comparison of clustering properties of point processes[END_REF]), where the last inequality comes from the fact that [22, eq.(1.6)

P ! u (Φ K (B) = 0) = P(Φ K! u (B) = 0) ≤ P(Φ Ku! P (B) = 0) = e -E[Φ Ku! P (B)] = e -E[Φ K! u (B)] ≤ e -(E[Φ K (B)]-p) = e -(E[Φ(B)]-p) ,
E[Φ K! u (B)] ≥ E[Φ K (B)] -p (see

in the Supplemental content], for a proof).

We can then write, for any r ≥ 0, and any i ∈ {1, • • • , 6},

P ! u (Φ K (∠ i (r)) = 0) ≤ e -(E[Φ(∠ i (r))]-p)
≤ e p-λ πr 2 6 .

Injecting this inequality in eq. (3.25), it gives us

E !u [V Φ+δu (u 1 )] ≤ 6 R + e p-λ πr 2 6 2πrdr = c p ,
where c p = 6e p R + e -λ πr 2 6 2πrdr = ce p < ∞. Going back to eq. (3.22), we obtain

E p E !u [|f (p) u (Φ + δ u , Φ)|]M Φ (p) (du) ≤ c 2 p-1 (p -1)! e p (R 2 ) p 1 |ψ j (u 1 )|M Φ (p) (du) ≤ c 2 p-1 (p -1)! e p λ p (R 2 ) p 1 |ψ j (u 1 )|du = c 2 p-1 (p -1)! e p λ p π p-1 R 2 |ψ j (u)||u| 2(p-1) du, (3.26)
where the second inequality comes from the fact that, for a determinantal point process, for any u ∈ E p , ρ (p) (u) ≤ λ p . Now, as ψ j is a Morlet wavelet (cf. eq. (3.3)), we can give an upper bound for |ψ j (u)|:

|ψ j (u)| ≤ 1 2πσ 2 2 -2j e -2 -2j |u| 2 /(2σ 2 ) |e i2 -j ξ 0 .u -c ψ | ≤ 2 -2j c ψ e -2 -2j |u| 2 /(2σ 2 ) ,
where c ψ is a positive constant. Therefore,

R 2 |ψ j (u)||u| 2(p-1) du ≤ 2 -2j c ψ R 2 e -2 -2j |u| 2 /(2σ 2 ) |u| 2(p-1) du ≤ 2 -2j c ψ [0,2π] dθ R + e -2 -2j r 2 /(2σ 2 ) r 2p-1 dr ≤ 2 -2j c ψ 2π(2 2j 2σ 2 ) p R + e -r 2 r 2p-1 dr ≤ 2 -2j c ψ 2π(2 2j 2σ 2 ) p 1 2 R + e -r r p-1 2 r -1 2 dr ≤ 2 -2j c ψ π(2 2j 2σ 2 ) p R + e -r r p-1 dr ≤ πc ψ 2 2j(p-1) σ 2p 2 p (p -1)!
Injecting this into eq. (3.26), we obtain

E p E !u [|f (p) u (Φ + δ u , Φ)|]M Φ (p) (du) ≤ c 1 (4πσ 2 eλ2 2j ) p , (3.27)
where c 1 is a positive constant, independent of p. This shows that, for all p,

(R 2 ) p E !u [|f (p) u (Φ + δ u , Φ)|]M Φ (p) (du) < ∞, so the condition Equation (3.15) of proposition 3.4.1 is satisfied. If furthermore 4πσ 2 eλ2 2j < 1, then (R 2 ) p E !u [|f (p) u (Φ + δ u , Φ)|]M Φ (p) (du) ---→ p→∞ 0,
and eq. (3.18) holds. In conclusion, factorial moment expansions of any order exist for all scattering moments S 1 Φ(j) of any scale j, for a determinantal process marked by the volumes of the Voronoi cells. Also, the infinite expansion exists for all sufficiently small scales with respect to the intensity of the determinantal process. Remark that the fact that Φ is a determinantal point process introduces the factor e p in the bound eq. (3.27) w.r.t. a Poisson point process.

Remark. Note that the wavelet transform (eq. (3.9)) is a linear function (shot-noise functional) of the (marked) point process. In other words, it is a first-order U -statistics of Φ m ; cf [START_REF] Last | Lectures on the Poisson process[END_REF]Section 12.3]. The fact that scattering moments S 1 Φ(j, θ) are defined via the modulus of the wavelet transform makes them dependent on all higher-order correlation functions of Φ: the factorial moment expansions of S 1 Φ(j, θ) involve all moment measures (in contrast to the square norm of the wavelet transforms, which can be represented using the first and second-order U -statistics and thus their expansion involves only the first two correlation functions). Consequently, the first order scattering moments are supposed to capture more information regarding the intrinsic dependence of the points than just the pairwise correlations. We shall see in our numerical study that this information allows one to efficiently recover some geometric marks which do depend on higher order correlations, e.g, the Voronoi cell characteristics.

Wavelet phase harmonics covariance

We shall now take interest in another class of representations based on the wavelet transform of point processes. Scattering moments decompose the information of a point process into different features at different scales, and those features in different subfeatures at different scales. However, this representation does not take into account the dependencies between the different features at different scales. Furthermore, the modulus operator introduced in the definition of scattering moments remove the phase information of the wavelet coefficients, which can contain important information about the geometry of the observed pattern, as illustrated in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF][START_REF] Zhang | Maximum entropy models from phase harmonic covariances[END_REF]. For this reason, in [START_REF] Zhang | Maximum entropy models from phase harmonic covariances[END_REF], the authors introduce covariance moments between wavelet coefficients, which phases are adjusted by a non-linear operator. In this section, we shall review this representation, for stationary random measures defined on (R 2 , B).

As stated in Section 3.1, for a wavelet ψ with central frequency ξ 0 , one builds a wavelet family by dilating and rotating ψ. The wavelet at scale j and angle θ has central frequency ξ := 2 -j r -θ ξ 0 . In this section, we shall, for notations simplicity, identify the central frequency ξ with the indices of the wavelets, (j, θ). We are interested, in this section, in capturing the dependencies between wavelet coefficients of a stationary point process Φ, at possibly different scales, orientations, and spatial positions. To study such dependencies, one can naturally look at the covariance moments between wavelet coefficients, that is, moments of the form

CΦ(ξ, ξ , τ ) = E[Φ ψ ξ (u) Φ ψ ξ (u -τ ) * ], (3.28) for (ξ, ξ ) ∈ (Z × [0, 2π]) 2 , τ ∈ R 2 .
Note that in this case, we do not need to subtract the expected averages, since the wavelets ψ ξ have zero means. Furthermore, as Φ is stationary, the value in the right hand side of eq. (3.28) does not depend on u, that can be set to 0. Suppose that Φ satisfies the conditions of Section 2.8, that is, it is stationary, has a non-null finite intensity, and with a σ-finite second-order factorial moment measure M Φ (2) absolutely continuous w.r.t. the Lebesgue measure. Then, using eq. (2.23), we obtain

CΦ(ξ, ξ , τ ) = R 2 ψξ (ω) ψ * ξ (ω)e -iω.τ b Φ (ω)dω. (3.29)
Therefore, computing CΦ(ξ, ξ , τ ), for τ ∈ R 2 , is equivalent to computing the Fourier transform of ψξ ψ * ξ b Φ . Recall also that the wavelets considered in this work are localized in frequency, i.e. there exists a constant C such that ψ(ω) 0 if |ω -ξ 0 | > C . This means that, the (approximate) frequency support of ψ ξ is contained in a ball centered around 2 -j r -θ ξ 0 and of radius 2 -j C . If C is small enough, then for ξ = ξ , ψ ξ and ψ ξ approximately have disjoint frequency supports, meaning that

ψξ (ω) ψ * ξ (ω) 0, ∀ ω ∈ R 2 .
This implies that the covariance between wavelet coefficients will be close to zero, even though they may not be independent.

To remedy this problem, one can apply a non-linear operator to the wavelet coefficients, in order to modify their frequency support, so that they may overlap. The phase harmonics operator, that we shall now define, is particularly suited for that purpose. Definition 3.5.1 (Phase harmonics operator). Let z ∈ C, and denote |z| its modulus, and ϕ(z) its phase, i.e. z = |z|e iϕ (z) . The phase harmonics operator is defined, for every

k ∈ Z, by [z] k := |z|e ikϕ(z) . (3.30) Note that, for any z ∈ C, [z] 0 = |z|, [z] 1 = z, and [z] -1 = z * . More generally, for any k ∈ Z, ([z] k ) * = [z] -k and |[z] k | = |z| for k ∈ Z.
For k ∈ Z and z ∈ C, the phase harmonics operator has the same effect on the phase of z than taking z to the power k. However, using estimators of high order moments can prove problematic, due to their possibly large variance. This operator, however, does not change the modulus of z.

We can apply the phase harmonics operator to the wavelet coefficients of µ ∈ M to adjust their phase. Definition 3.5.2 (Wavelet phase harmonics). Let µ be a measure on (R 2 , B), and ψ ∈ L 1 (R 2 , µ). We define the wavelet phase harmonics (WPH) of µ by

[µ ψ ξ ] k = |µ ψ ξ |e ikϕ(µ ψ ξ ) , (3.31) for ξ ∈ Λ = Z × [0, 2π],
and k ∈ Z.

The phase harmonics operator was introduced in [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF], and used in [START_REF] Zhang | Maximum entropy models from phase harmonic covariances[END_REF] to build a class of representations, and define maximum entropy models for non-Gaussian stationary processes. The idea of multiplying the phase of complex wavelet coefficients to compute their dependencies was already present in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF], although it was not conceptualized in an operator, that can be computed for any k ∈ Z.

Frequency transposition

The phase harmonics operator modifies the phase of complex wavelet coefficients, without modifying its modulus. In [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF], the authors explain that such an operator, called a frequency transposition, performs a dilation in the Fourier domain. They show that, in the frequency domain, if ψj 1 ,θ 2 is approximately supported on the ball B ξ (C ) of radius C , centered at ξ ∈ C, then [ ψj 1 ,θ 2 ] k is approximately supported on B kξ (kC ). This property, illustrated in [START_REF] Zhang | Maximum entropy models from phase harmonic covariances[END_REF] in the case of realizations of stationary processes, allows us to capture dependencies of wavelet coefficients across scales and angles, as we shall detail next. Definition 3.5.3 (WPH covariance moments). Let Φ be a stationary point process on (R 2 , B), and Ψ = {ψ 0 , ψ ξ } ξ∈Ξ=Z×[0,2π] be a wavelet family such that all wavelet coefficients have finite second-order moments. For any ξ, ξ ∈ Ξ, k, k ∈ Z, τ ∈ R 2 , let us note γ = (ξ, ξ , k, k , τ ). The WPH covariance of Φ is defined by

K γ Φ := Cov [Φ ψ ξ ] k (u), [Φ ψ ξ ] k (u -τ ) . (3.32)
Note that, as Φ is stationary, the r.h.s. of eq. (3.32) does not depend on u ∈ R 2 .

Note that the WPH covariance representation does not need to include the covariance between coefficients for all indices ξ, ξ ∈ Ξ, k, k ∈ Z, τ ∈ R 2 , because some of these covariance moments are nearly zero (as explained by the frequency transposition property). A way to choose the coefficients in the representation shall be detailed in Section 5.4, in the case where this representation is used to define a maximum entropy model.

Estimation and coefficient selection

Estimation The wavelet-based representations presented in Sections 3.2 and 3.5 capture information about the distribution L Φ of Φ by computing mathematical expectancies of non-linear functions applied to the the wavelet transform of Φ (cf. Section 3.1). In other words, we build representations of the form

R(Φ) = E[f (Φ)],
for non-linear functions f .

If Φ is a stationary, ergodic point process, these expectations can be estimated from a realization φ of Φ, observed on a finite window W , by spatial averaging over the window. That is, we build our estimators as

RW (Φ) := 1 |W | W f (Φ| W )(u)du,
where Φ| W = Φ ∩ W . Note that, if we were given Φ over the whole space R 2 instead of the bounded window W , we could define

RW,R 2 (Φ) := 1 |W | W f (Φ)(u)du,
and this estimator would be unbiased. Furthermore, Birkhoff's ergodic theorem (cf. eq. (2.21)) would tell us that

lim |W |→∞ RW,R 2 (Φ) = R(Φ), P-a.s.
As we only observe Φ| W , that is the points of Φ that lie inside W , the wavelet coefficients Φ| W ψ j,θ (u) at points u ∈ W near the border of W cannot take into account points outside W . One can therefore consider that the observation sample has no point outside W , or consider other types of border corrections. For instance, in this work, we shall usually perform a periodic boundary correction, i.e; consider that the sample repeats itself outside W , when W is a square window. This can be done by using periodic wavelets {ψ s ξ }, defined as

ψ s ξ (u) := n∈Z 2 ψ ξ (u + sn), (3.33)
where s > 0 is the size of the square window W . We shall however often make the approximation that |W | is large enough, so that border effects are negligible, and the ergodic theorem holds.

Finite representation

In Section 3.1, the wavelet transform has been defined with an infinite wavelet family {ψ j,θ }, with j ∈ Z and θ ∈ [0, 2π]. For practical purposes, we need to select a finite number of such coefficients.

For the angular parameter θ, this can be done by simply selecting discrete, evenly spaced angles in [0, 2π]. Thus θ shall be parametrized by θ = 2π l L , with l ∈ {0, • • • , L-1}, for some L ≥ 1. In our simulations, we shall usually take L = 4 or L = 8.

The scale parameter j will be taken in {0, • • • , J -1}, for some J ≥ 1. The highest frequency captured by the wavelet transform is therefore determined by the central frequency ξ 0 of the mother wavelet. Equivalently, the wavelet family will not capture structures smaller than C, the size of the approximate spatial support of ψ. The parameter ξ 0 must therefore be chosen according to the task at hand. In practical applications, one deals with samples observed inside a bounded window W . The maximal scale parameter J is typically chosen so that the approximate support of the wavelets ψ J-1,θ (those with maximal support) is contained in the observation window, i.e. 2 J C ≤ s, for a square window W of size s. Scales equal or larger than J are carried by the low-pass filter ψ 0 , whose frequency support is centered at ξ = 0.

Properties The estimators of scattering moments and wavelet phase harmonics covariance moments enjoy properties that make them suitable for statistical learning of non-linear functions of a (possibly marked) point pattern, a point process, or its distribution. First, being defined with a global spatial average, these representations are invariant to translations, which can be a desirable property (e.g. if the function to approximate is itself also translation invariant). Another property that is sometimes required of such descriptors is some notion of stability. In other words, if a point pattern undergoes some kind of small deformation, then the associated value of the function to approximate might not change too much. Thus, the descriptor of this point pattern should not be too affected by this deformation either. In the case of estimators of first order scattering moments, such stability property can be easily proven.

Proposition 3.6.1. Let φ ∈ M be a finite counting measure, i.e. φ = N i=1 δ u i . Let m ∈ R N , m = (m 1 , • • • , m N ). We define φ = N i=1 m i δ u i .
Let us also define the marks and positions perturbations: u ∈ (R 2 ) N and m ∈ R N . We note φ = N i=1 δ u i +u ,i , and

φ = N i=1 (m i + m ,i )δ u i +u ,i . Finally, note Ŝ1 φ = R 2 | φ ψ|(u)du,
for some wavelet ψ, and similarly for φ . Suppose that ψ is integrable, and such that ∀ u, t ∈ R 2 , |ψ(u + t) -ψ(u)| ≤ t ς(u), for some integrable function ς. Then,

| Ŝ1 φ -Ŝ1 φ | ≤ ς 1 m 2 |u| 2 + ψ 1 m 1 where |u| = ( u ,i ) i=1,••• ,N . Proof. Let us denote ξ = N i=1 m i δ u i +u ,i . Using the fact that, for all z, z ∈ C, ||z|-|z || ≤ |z -z |, we write | Ŝ1 φ -Ŝ1 φ | = | Ŝ1 φ -Ŝ1 ξ + Ŝ1 ξ -Ŝ1 φ | ≤ | Ŝ1 φ -Ŝ1 ξ| + | Ŝ1 ξ -Ŝ1 φ | ≤ i |m i | R 2 |ψ(u -u i ) -ψ(u -u i -u ,i )|du + i |m ,i | R 2 |ψ(u -u i -u ,i )|du ≤ ς 1 i |m i | u ,i E + ψ 1 i |m ,i | ≤ ς 1 m 2 |u| 2 + ψ 1 m 1 .
In this chapter, we have presented two classes of representations of stationary point processes, based on a multi-scale decomposition of the geometric information with a wavelet transform.

In the remainder of this work, we shall study the adequacy of these representations to perform statistical learning and probabilistic modelling tasks.

Chapter 4

Scattering moments for point process marks regression

In this chapter, we propose to address the problem of learning marks of point processes: an unknown function, called a marking function, takes as argument a point pattern, and produces numerical values, called marks, associated with each point of the pattern. Our goal is to learn this function, from examples of point patterns with observed marks, in order to predict unknown marks for new point patterns.

Our particular motivation comes from the problem of learning how the cell loads in wireless networks depend on the geometry of the base stations (and possibly of the traffic demand). One may want to learn this dependence directly from real data, collected in the existing operational networks, to predict the loads of base stations for different configurations of base stations, and/or different traffic demand (see more details in the Related works paragraph). In this work, we consider some simple, generic marks, often appearing in this context. They are produced by the standard shot-noise interference model, the nearest neighbour distance, and some characteristics of the Voronoi cells (cf. Section 2.6 Figure 4.1). Our goal is to to understand the efficiency of our proposed method to represent the dependence of the marks on the geometry of point patterns, rather than to solve the original cell load problem for some particular real network data.

Learning geometric marks via local or global representation

The problem of learning the marking function, as posed above, does not seem amenable to any direct statistical approach, due to the structure of the space of marked point measures. To overcome this difficulty, we need some suitable representation of the geometry of point patterns, and we propose two approaches in this regard.

As a baseline, we propose to estimate the mark of each point using a statistical regression model based on the local distance matrix of a suitable chosen vicinity of this point. In this relatively simple approach, the training data set consists of an ensemble of central marked points, surrounded by some number, say K, of their (non-marked) neighbours. Observe that the dimension of the representation of the local geometry of the point (local distance matrix) increases as K 2 . For highly non-local marking functions (e.g. related to power-law shot-noise function), one might need to chose a large number of neighbours K, proportional to the total number of points, possibly making this baseline approach inefficient.

When the marks have non-local dependence, or when we have no prior knowledge of the dependence range, we further propose to use the scattering moments representation (cf. Section 3.2) to capture the geometry of marked point patterns. Recall that the scattering moments are defined as the spatial average of a cascade of convolutions with a wavelet family, followed by a complex modulus operator. The property of the modulus operator makes this representation locally invariant to small translations ( [START_REF] Bruna | Invariant scattering convolution networks[END_REF]), and the global averaging implies a global translation invariance. Furthermore, they have been shown ( [START_REF] Mallat | Group invariant scattering[END_REF][START_REF] Bruna | Invariant scattering convolution networks[END_REF], Proposition 3.6.1) to be stable to small deformations of the input signal. These invariance and stability properties make them useful in signal processing, in particular in relation to statistical learning. Indeed, if the information content of a signal is typically not (strongly) affected by translations, rotations, and small deformations, similar properties of the signal representation allows one to capture this content in a more concise way, and hopefully learn its intrinsic structure from a smaller number of signal samples. The pertinence of this approach has already been demonstrated in various contexts (see e.g. [START_REF] Bruna | Invariant scattering convolution networks[END_REF][START_REF] Sifre | Rotation, scaling and deformation invariant scattering for texture discrimination[END_REF][START_REF] Hirn | Wavelet scattering regression of quantum chemical energies[END_REF]).

Using scattering moments to represent marked point processes, our learning problem can be addressed in the following two steps (also depicted on Figure 4.2):

• We build a statistical regression model on scattering moments of the marked point patterns, with the explanatory variables being the scattering moments of the nonmarked point patterns. This model is computed on the training data consisting of point patterns with observable marks, and is meant to be used to estimate the (marked) scattering moments of new point patterns, for which marks are not observed.

• We estimate (reconstruct) the marks of new point patterns, for which marks are not observed, from their estimated scattering moments. It is a non-convex optimization problem, that we numerically solve using the L-BFGS-B algorithm [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF][START_REF] Zhu | Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization[END_REF].

Compared to our baseline approach, for which the size of the representation of the local geometry of a point grows quadratically with its number of neighbours, the scattering representation grows logarithmically with the neighbouring size. However, it may suffer from possible errors, introduced in the reconstruction phase (absent in the benchmark approach). The overall benefits of the scattering approach become significant for highly non-local marking functions. In general, the quality of the regression of the scattering moments depends on the sensitivity of the marking function to small point pattern deformations. On the other side, in the reconstruction phase (recovering of marks from the true or estimated scattering moments) significant errors consist in swapping a large and a small mark of two neighbouring points (not leading to a significant modification of the scattering moments). The quality of this phase of the approach depends thus on the existence of clusters of points in point patterns -a problem already studied in the literature; cf Section 4. The Poisson point process, considered in this paper, exhibits a baseline type of clustering.

Related works

Stochastic-geometric study of cellular networks expanded rapidly in recent years, primarily through analytic results regarding Poisson network models ( [START_REF] Błaszczyszyn | Stochastic geometry analysis of cellular networks[END_REF]). Performance evaluation of operational wireless cellular networks, in particular the quality of service perceived by users in function of the traffic demand, is a complex problem involving stochastic and geometric modelling of several network layers. A key element of this problem is the analysis of the cell loads, which non trivially depend on the geometric configuration of serving base stations and their traffic demands, and capture in a concise way the quality of service offered by the individual cells. A detailed physical model of cell loads was proposed in [START_REF] Siomina | Analysis of cell load coupling for lte network planning and optimization[END_REF] and revisited in [START_REF] Błaszczyszyn | How user throughput depends on the traffic demand in large cellular networks[END_REF][START_REF] Błaszczyszyn | Spatial disparity of qos metrics between base stations in wireless cellular networks[END_REF][START_REF] Błaszczyszyn | Performance analysis of cellular networks with opportunistic scheduling using queueing theory and stochastic geometry[END_REF], including the validation with respect to some data collected in operational networks.

Recently, the prediction of the cell loads has been the subject of a machine learning study published in [START_REF] Awan | A robust machine learning method for cell-load approximation in wireless networks[END_REF], though not as the geometric problem posed in our paper. Indeed, in this earlier paper, the number and the locations of base stations are fixed, and the problem consists in learning the loads of this given configuration of stations in function of their traffic demands. This (non-geometric) problem is cast and solved in the framework of monotone interpolation of Lipschitz functions.

The problem formulated and studied in our paper allows one to address the orthogonal question of learning the cell loads, given constant traffic demand, in function of the geometry of the network. Combining the two approaches, left for future work, may allow one to predict the performance of new geometric configurations of base stations and traffic demand, by collecting the cell load data in existing operational networks, thus offering an alternative to building and solving complex physical models.

Observe that, even for a fixed number of base stations, the two dimensional location of each base station on the plane makes the dimension of our input data two times higher than the one considered in [START_REF] Awan | A robust machine learning method for cell-load approximation in wireless networks[END_REF]. The use of wavelet scattering transforms allows one to leverage geometric invariants, such as translation and rotation, and thus significantly reduce the dimension of the representation. Scattering moments have been proven useful in this context in a number of tasks such as quantum molecular energy prediction [START_REF] Hirn | Wavelet scattering regression of quantum chemical energies[END_REF][START_REF] Eickenberg | Solid harmonic wavelet scattering: Predicting quantum molecular energy from invariant descriptors of 3d electronic densities[END_REF] or texture classification [START_REF] Sifre | Rotation, scaling and deformation invariant scattering for texture discrimination[END_REF].

Geometric marks (also called score functions) have received a lot of attention in stochastic geometry and spatial statistics, where they represent some interaction of a given point with the whole point pattern. If this interaction is local in some sense, and the underlying point process exhibits some decay of correlations, then it is possible to establish asymptotic results (including central limit theorems) regarding the sums of the geometric marks in increasingly large windows (cf e.g. [START_REF] Penrose | Weak laws of large numbers in geometric probability[END_REF][START_REF] Baryshnikov | Gaussian limits for random measures in geometric probability[END_REF][START_REF] Błaszczyszyn | Limit theory for geometric statistics of point processes having fast decay of correlations[END_REF]. These results can be used to study large-scale asymptotics of the scattering moments, as in [START_REF] Bruna | Intermittent process analysis with scattering moments[END_REF].

We mentioned earlier in this introduction that the quality of reconstruction of the marks depends on whether point process is regular, or exhibits clusters of points. These notions are formalised in [START_REF] Błaszczyszyn | On comparison of clustering properties of point processes[END_REF][START_REF] Błaszczyszyn | Clustering comparison of point processes, with applications to random geometric models[END_REF], and it is argued that Poisson point process can be considered as a reference type of point clustering, to which more regular point processes (e.g. determinantal ones) and more clustering ones (e.g. permanental and Cox processes) can be compared.

The remaining part of this chapter, based on [START_REF] Brochard | Statistical learning of geometric characteristics of wireless networks[END_REF], is organized as follows: in Section 4.1, we describe several models of geometric marks and formulate the problem of their statistical learning. We then present our main approach to solving this problem, based on scattering transforms of marked point processes. We also describe the benchmark approach, to which our method shall be compared. The numerical results of both approaches applied to the considered mark models are presented in Section 4.2.

Statistical learning of geometric marks

In this section, we formally present the problem of learning geometric marks of point processes, and formulate possible solutions to address this problem. Our main approach, described in Section 4.1.4, uses the scattering moments of marked point patterns, which are presented in Section 4.1.3. An alternative method is proposed in Section 4.1.5, based on a distance matrix representation of point neighbourhoods. This alternative approach is used as a comparative benchmark for our numerical study in Section 4.2.

Geometric characteristics of wireless networks

In stochastic-geometric modeling of networks (e.g. wireless networks), one usually represents locations of objects (e.g. transmitters/receivers) as points of a point process, and their characteristics as marks. These characteristics can depend not only on the given object, but also on the geometry of the network (at least locally, e.g. the surface of the cell served by the base station, the extra-cell interference). These marks, called in what follows geometric marks (cf.Section 2.6), are the object of interest in this study.

In this work, we shall consider consider stationary marks. Recall from Section 2.6, that they are defined as mappings m : R 2 × M → R such that, for any u, v ∈ R 2 , and

any φ ∈ M, m(u -v, S v φ) = m(u, φ).
We denote by M i := m(u i , Φ) the geometric mark of the point u i of Φ. Note that the process Φ := i δ (u i ,M i ) ∼ i M i δ u i is a stationary marked point process.

In the remainder of this chapter, we shall consider the basic examples presented in Section 2.6. For a better understanding of the context, we interpret the points u i as locations of base stations of some cellular network. The shot-noise mark (cf. Example 2.6.2) can be considered as the total power (usually interpreted as interference) received at a given base station in the network, transmitting with unit power. The Voronoi cell (cf. Example 2.6.4) is a fundamental model for cellular networks. Its area and moment of inertia can be interpreted as very simple proxys for the traffic demand and the cell load of a base station (see [START_REF] Błaszczyszyn | Stochastic geometry analysis of cellular networks[END_REF]Section 4.1.8]). The Voronoi shot-noise mark (cf. Example 2.6.5), can be also interpreted as the interference, however with the stations transmitting the signals with a power proportional to their Voronoi cells. (Similar kinds of dependence, with more complex expressions, can be recognized in the cell load model of [START_REF] Siomina | Analysis of cell load coupling for lte network planning and optimization[END_REF].)

Problem formulation

Let us now formulate the main problem studied in this chapter, that is the problem of the learning of geometric marks. Suppose the marking function m is not known explicitly. One observes only some realizations of the marked point process Φ with points restricted to some finite observation window W . Denote these realizations by

φ k = i m i (k)δ u i (k) , with u i (k) ∈ W , k = 1, . . . , n,

Scattering moments of marked point processes

Recall from Section 3.2, the wavelet transform of a marked point process is defined by

Φ ψ : = R 2 ψ(v -u)m(Φ, u) Φ(du) = u i ∈Φ M i ψ(u -u i ),
and that first and second order scattering moments (Section 3.2) can similarly be extended to marked point processes. Scattering moments have not yet been extensively theoretically studied in the context of the theory of point processes. Some asymptotic properties at small and large scale (j → -∞ and j → ∞ respectively) have been presented in section 3.3, extending results of [START_REF] Bruna | Invariant scattering convolution networks[END_REF].

As explained in Section 3.6, for a given realization φ of Φ observed in a finite window W , empirical scattering moments Ŝ1 Φ(j, θ) and Ŝ2 Φ(j 1 , θ 1 , j 2 , θ 2 ) are obtained by replacing the expectation by a spatial averaging of the fields | φ ψ (j,θ) (u)| and

|| φ ψ (j 1 ,θ 1 ) | ψ (j 2 ,θ 2 ) (u))|,
respectively, over u ∈ W . When W increases, these (empirical) moments become asymptotically non-biased estimators of S 1 Φ(j, θ) and S 2 Φ(j 1 , θ 1 , j 2 , θ 2 ), respectively, provided Φ is ergodic (under some technical conditions about the wavelet ψ), (cf ergodic theorem for point processes [START_REF] Daley | An introduction to the theory of point processes: volume II[END_REF]Theorem 13.4.III], or Section 2.7). For some other properties of the scattering moment estimators see [START_REF] Bruna | Intermittent process analysis with scattering moments[END_REF]Section 5.1].

The empirical scattering moments are calculated over a finite square window W by periodizing the input signal. As a consequence, the value of the mark is calculated as if the point process were periodic.

As explained in Section 3.6, in order to define a finite wavelet family, we need to choose a finite number of scales j ∈ {0, • • • , J -1}. Assume that there is also a wavelet at scale j min small enough so that it separates all the points in every pattern (theoretically corresponding to j min ∼ -∞). At this scale, Ŝ Φ(j min , θ) is proportional to the empirical mean 1 |W | W 2 mφ(du), and does not depend on θ. The first and second empirical scattering moments calculated on Φ form finitedimensional vectors. We denote by Ŝ Φ, and Ŝ2 Φ, respectively, the vector of the first moments and the joint vector of the first and second moments. If there is no ambiguity, for simplicity, in what follows we shall simply call them scattering moments of Φ (without distinction between first and second order).

Learning of marks via scattering moments

Recall from Section 4.1.2 that our goal is to learn the marking function m from the training set of data, which consists of examples of realizations of a marked point process Φ in a finite window. Note that this problem consists in the interpolation of the function m on the space R 2 × M, and, due to the complexity of this space, it is not amenable to any direct statistical approach. To overcome this difficulty, we map this original problem to some finite dimensional regression problem, and solve it using classical tools.

More specifically, in our main approach, we shall capture the function m through the relation between the vector of the first order scattering moments Ŝ φ of the marked point pattern, and the two order moments Ŝ2φ of the non-marked one. This relation is established using a regression model described in Section 4.1.4. In order to be able to use it to estimate the marks m i = m(u i , φ), we need next to solve an inverse problem described in Section 4.1.4. It consists in reconstructing marks from the regressed scattering moments, knowing the positions of the points in φ.

Linear regression

Let X k := Ŝ2φ k and Y k := Ŝ φ k , k = 1, . . . , n, be the vectors of the (empirical) scattering moments, calculated for the training data set consisting of n realizations of the point process Φ, where

φ k = i δ u i (k) and φ k = i m i (k)δ u i (k) (cf.
[β(λ)X k + β 0 (λ) -Y k (λ)] 2 + γ(λ)||β(λ)|| 2 ,
for some regularization parameter γ(λ) ≥ 0, chosen by cross-validation (to minimize this squared residuals on the validation set: a subset of the training set), where • is the Euclidean norm. This model admits a well known explicit solution in the form

[ β0 (λ), β(λ)] t = (X t X + γ(λ)I) -1 X t Y(λ), (4.2)
where, for all k = 1,

• • • , n, X(k, •) = [1, X k (•)], Y(λ)(k) = Y k (λ)
= i m i (k)δ u i observing only its points φ = i δ u i Ŝ φ(λ) := β(λ) Ŝ2φ + β0 (λ),
where Ŝ2φ is the vector of the first and second order scattering moments calculated on φ (without marks). Remember, expression (4.2) requires the tuning of the regularization parameters γ(λ) ≥ 0, usually needed in high dimensional regression problems when the matrix X t X is not invertible. The ordinary least square (OLS) estimator corresponding to γ(λ) = 0 is usually not performing well in this case.

Reconstruction

Using linear ridge, we calculate approximations Ŝ φ(j, θ) of the scattering moments of a new marked configuration φ, observing only φ. Denote the whole vector of Ŝ φ(j, θ) by Ŝ φ. From Ŝ φ we estimate (reconstruct) unknown marks m i of φ by looking for a solution to the following minimization problem arg min

φ :φ =φ || Ŝ φ -Ŝ φ|| 2 , ( 4.3)
where we minimize over all arbitrarily marked configurations φ having the same atoms as φ, and Ŝ φ denotes the scattering moment calculated for φ . It should be noted that (4.3) is a non convex optimization problem. To solve it, we use the L-BFGS-B algorithm, which is a limited-memory algorithm for solving large nonlinear optimization problems subject to simple bounds on the variables (cf. [START_REF] Liu | On the limited memory bfgs method for large scale optimization[END_REF][START_REF] Zhu | Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization[END_REF]).

Learning via local distance representation

In the following approach, taken as a benchmark, we consider each marked point (u i (k), m i (k)) of each realization φ k , k = 1, . . . , n, of the training set, along with some neighborhood, as one element of the new training data set. Therefore, in this case, the training set consists of point patterns having a marked point in their center, surrounded by some neighboring points. Using the linear ridge regression method described in Section 4.1.4, we regress the marks of those central points with respect to the vectors containing all inter-point distances in the considered point neighborhood, ordering the points according to the distance to the central point. Note that, in this approach, there is no reconstruction phase, as the marks are directly approximated. The main parameter of this approach is the number of neighbouring points K of each central point.

Numerical results

In this section, we provide the details of our numerical study of the problem considered in this chapter. We begin by describing some general assumptions and procedures.

General numerical framework

Scattering moments approach

For our numerical experiments, for each specific geometric model, we create a data set (denote it by X ) of Poisson point patterns φ k (with constant intensity to be specified) using the R software and its package Spatstat for point process analysis [START_REF] Baddeley | Spatstat: an r package for analyzing spatial point patterns[END_REF]. The points are considered in the unit square window, and their marks are analytically computed (according to the given model) using Spatstat with the window mapped to the torus, thus leading to the marked point patterns φ k , k = 1, . . . , n. The size of the data set X is n = 10k marked point patterns. The (empirical) scattering moments are computed on these point patterns (with and without marks) using ScatNet software [START_REF] Andén | Scatnet[END_REF] developed in Matlab (no Spatstat implementation is available yet). It uses the Morlet wavelet ( [START_REF] Goupillaud | Cycle-octave and related transforms in seismic signal analysis[END_REF], cf. eq. (3.3)). Because this software uses raster images, we convert each marked point pattern of X into images of size 2 7 × 2 7 pixels (removing images with points corresponding to the same pixel).

The following family of scattering moments are computed. Recall from Section 4.1.3 that there is a minimal scale j min such that, at this scale, the first order scattering moments correspond to the empirical mean measure, and do not depend on the angle θ (we can therefore take only one). We further choose J = 7, and L = 8 (i.e. θ ∈ {0, • • • , 7π/8}). Thus, there are 1 + 8 × 7 = 57 first order scattering moments (the dimension of the regressed vectors Ŝ φk is 57), and 8 × 8 × 7×6 2 = 1344 second order scattering moments (thus making the dimension of the explanatory vectors Ŝ2φ k equal to 57 + 1344 = 1401).

We use the linear ridge regression, described in Section 4.1.4, on the data set X . To optimize the regression parameters, we use a 5-fold cross-validation [97, Section 1.4.8] on X .

Having calculated the estimators Ŝ φ of the first order scattering moments Ŝ φ for the point patterns in the test set, we use the L-BFGS-B algorithm to solve the inverse problem (4.3), i.e. to reconstruct the marks. This is a steepest descent algorithm for which it is important to optimize (ideally via cross-validation) the number of iterations; to be explained in Section 4.2.3. Using a rule of thumb (no formal cross-validation), we fix the number of steps so as to minimize the mean square error (MSE) on the test data set.

Benchmark

As explained in Section 4.1.5, for every image in the data set X , we consider each point of the image, along with its K neighbours as an element of a new data set X . More precisely, we take the first 20k points of the images of X , and create X with K neighbours of the chosen point. The value of K depends on the range of dependence of a given mark model; it is experimentally discovered for each mark, as will be explained in Section 4.2.3. We use X to regress the marks of the central points directly, using linear ridge regression with respect to the vector of dimension K(K -1) of the local distance matrix. Note that there is no reconstruction phase in the benchmark method.

Validation methodology

To test our main approach, we produce an independent data set of 100 marked point pattern realizations, and we use the first 100 points of this set for the benchmark approach.

For all points in the respective test sets, we compute Q-Q plots, the root mean square error (RMSE), the normalized RMSE with normalization by the range (maxmin) of the marks (NRMSE1), and the normalized RMSE with normalization by the mean of the marks (NRMSE2).

In order to study the quality of the representation of the mark point patterns by their first order scattering moments (independently of the regression), we also perform the reconstruction of marks 4.1.4 directly from the exact (and not regressed) scattering moments Ŝ φ. Note, these latter are not exact marks but marks reconstructed from exact empirical scattering moments of the given point pattern. The gap between them and the marks reconstructed from the regressed moments allows one to apprehend the error introduced by the regression and the one introduced by the reconstruction.

Results

We now present our numerical study of different mark models, presented in Section 2.6. The observed results are discussed in Section 4.2.3.

Shot-noise

We consider shot-noise marks (cf. Example 2.6.2), with the response function (r) := max(10r, 0.6) -3 . For this example the data set X consists of 10k realizations of Poisson point process with intensity 40 (recall, we consider the unit square window). In order to observe how the performance of our approach depends on the size of the data set, we first use only 5k elements of X as the training set for the linear ridge regression method. Next, we use the whole training set, with linear ridge regression. For the benchmark method, we choose a neighbourhood of K = 15 points. The number of iterations for the reconstruction from the exact scattering moments and estimated ones is taken, respectively, 30 and 4; see Section 4.2.3 for the explanation. 4.1 gives the reconstruction errors. Recall that the Q-Q plot (c) and the last column of the table represents the error of the reconstruction of the marks from the exact (and not regressed) scattering moments.

Nearest neighbour distance

The data set for this example is also made of 10k realizations of a Poisson point process with intensity 40. In this case, in Figures 4. and in Table 4.3. The neighborhood for the benchmark is very large K = 35 (more than the average number of points) because the sum of the areas of the Voronoi cells in the finite window is constant equal to the total window surface area, which introduces a strong global dependence for this mark. The number of iterations for the reconstruction from the exact scattering moments and estimated ones is taken, respectively, 30 and 6.

Voronoi cells moment of inertia

For this example, similarly to the previous example, the training set X consists of 10k realizations of a Poisson point process with intensity 30. The results are presented on Figures 4.9 and 4.10 and in Table 4.4. The neighborhood for the benchmark is K = 15 (note that the global dependence specific for the Voronoi surface area does exist here).

The number of iterations for the reconstruction, from the exact scattering moments and the estimated ones, are chosen as 150 and 8, respectively. The number of iterations for the reconstruction, from the exact scattering moments and the estimated ones, are chosen as 50 and 5, respectively.

Discussion

The following remarks can be formulated regarding the observed results.

Reconstruction from the exact and estimated moments

Observing the Q-Q plots and (N)RMSE's of the marks reconstructed from the exact first order scattering moments, we see 'how much information' they effectively carry regarding the marking function. While all marks are relatively well reproduced in this way, the quality of the reconstruction depends on the type of dependencies represented by a given mark. For example, the shot-noise and the surface areas of the Voronoi cell are more easy to represent than the nearest neighbour. This can be explained by different sensitivity (stability) to small deformations of the point pattern, with a precise formulation yet to be theoretically studied on the ground of point processes.

A typical significant error in the signal reconstruction, both from the exact and the estimated scattering moments, consists in the swap of a large and a small mark of two neighbouring points (e.g. the points number 36 and 37 on Figure 4.3), not leading to a significant modification of the considered scattering moments. This effect can be seen also on the Q-Q plots where many points significantly far from the diagonal appear in symmetric pairs. We believe these swapping errors should occur less often for more regular (less clustering) point processes than the considered Poisson one. Examples of such point processes are determinantal point processes; see [START_REF] Błaszczyszyn | On comparison of clustering properties of point processes[END_REF][START_REF] Błaszczyszyn | Clustering comparison of point processes, with applications to random geometric models[END_REF] for a clustering comparison theory. Future studies should investigate this issue.

As already mentioned, it is important to properly tune the number of iterations of the steepest descent algorithm used in the reconstruction phase, preventing it from going too deeply into potential local minima. We observe the following local-global reconstruction quality trade-off: while some number of initial iterations makes all the marks approach their right values, further iterations improve the quality of approximation of some subset of marks at the price of degrading this quality for the remaining ones, as illustrated in Figure 4.13. This trade-off can be observed on average by watching the RMSE on the test set, which first decreases and then increases. We use this observation to choose an optimal number of iterations. It is larger for the reconstruction from the exact scattering moments than for the reconstruction from the estimated ones, where this effect is less visible.

Choice of neighbourhood

A crucial benchmark parameter is the number of neighbors to be considered in the local representation. They should be selected in function of the type of mark dependence. If no a priori information is available, this can be done observing the RMSE on the test set as shown in Table 4

Quality of the regression

The regression relative errors

m (j, φ) := | Ŝ φ(j, θ) -Ŝ φ(j, θ)| Ŝ φ(j, θ)
calculated on the training and test set, presented on Figure 4.14, show that there is no overfitting in the regression of the scattering moments.

Scattering moments versus benchmark

With respect to the RMSE, the scattering moments approach outperforms the benchmark for all marks except for the Voronoi surface area with the neighborhood taking almost all points and, for obvious reasons, the nearest neighbor distance. A significant difference can be observed on the Q-Q plots of the shot-noise example in Figure 4.4, and the Voronoi shot-noise example in Figure 4.12, where the benchmark essentially fails to capture the marking function. This shows that the benchmark approach might not be appropriate in the case of long-range dependent marks.

Discussion

Motivated by the stochastic-geometry problems related to wireless networks, in this study we have discussed how to learn a point-marking function dependent on the configurations of points. We propose two different approaches to address the problem, using tools from statistical learning. The baseline approach extracts the geometric information for each point based on the matrix distance of its nearby points. It is then solved using the ridge linear regression method, and the difficulty lies in the choice of the number of nearby points. The other approach uses multi-scale wavelet scattering moments to define a global feature vector for all the points in the domain, and another feature vector for the points with marks. The relation between the two feature vectors is also learned using linear regression. We then solve a non-convex optimization problem to reconstruct the marks, from the predicted feature vector with marks. These feature vectors are translation-invariant and stable to deformation of the domain.

Depending on the nature of the clustering of the points and the regularity of the marking function, we find that the scattering moments predict as good as the baseline approach, showing that they capture well the geometric property. In case where marks depend on points in a non-local way, the scattering moments method outperforms the benchmark approach.

Regarding practical applications, it seems promising to combine geometry with local demand, in order to obtain an operational method for cell load prediction in cellular networks.

Future directions to be explored may include a combination of the two approaches, e.g. by representing the neighbourhood of each point not by its distance matrix, but by its scattering moments. Replacing the linear regression by others such as kernel regression or neural network approaches are also of potential interest. Additionally, the impact of regularity or clustering of the underlying point patterns on the quality of scattering representation ought to be better understood.

Limitations of scattering moments

This work demonstrated the ability of wavelet scattering moments to adequately capture important geometric structures in point patterns, and their amenability to be used in a learning task involving point processes. However, they do not fully characterize these processes. Indeed, scattering moments do not capture the dependencies between wavelet coefficients at different scales and orientations. Interpreting the first order scattering moments as measuring 'how much structure at each scale and orientation is present in the data', second order coefficients measure the partition of these structures between larger scales, but not the joint presence of structures at different scales. Additionally, due to the wavelet symmetries and the use of the modulus operator, this representation is invariant to changes of sign in the marking function, and to space reversal (u → -u), which probably explains the error due to the swapping of marks between nearby points. Other wavelet-based representations, such as the WPH covariance presented in Section 3.5, may help capture this kind of missing information. To study the quality of WPH covariance moments as descriptors of the geometry in point processes, we shall use them to define a model for stationary point process distributions, as detailed in the next chapter.

Chapter 5

Point process synthesis using wavelet phase harmonics and gradient descent

In order to generate new realizations of a stochastic process of which we have only one realization, one needs to build a probabilistic model which approximates the distribution of this process, and from which we can efficiently sample from. In this chapter, we are interested in generative models for stationary, ergodic point processes. Such models are of interest in a wide range of applications [64, Section 6], for instance biology [START_REF] Baddeley | Multitype point process analysis of spines on the dendrite network of a neuron[END_REF][START_REF] Diggle | Modelling the bivariate spatial distribution of amacrine cells[END_REF], ecology [START_REF] Wiegand | Handbook of spatial point-pattern analysis in ecology[END_REF], turbulent flows in atmosphere science [START_REF] Ducasse | Inertial particle collisions in turbulent synthetic flows: quantifying the sling effect[END_REF][START_REF] Oujia | Divergence and convergence of inertial particles in high-reynolds-number turbulence[END_REF][START_REF] Matsuda | Turbulent enhancement of radar reflectivity factor for polydisperse cloud droplets[END_REF][START_REF] Ducasse | Intermittent particle distribution in synthetic free-surface turbulent flows[END_REF], or cosmology [START_REF] Tempel | Bisous model-detecting filamentary patterns in point processes[END_REF][START_REF] Stoica | Detection of cosmic filaments using the candy model[END_REF]. In many of these domains, the observed patterns exhibit complex structures, with a large number of particles (such as filaments in cosmology, or vortices in turbulent flows).

More particularly, we seek to generate realizations formed by a large number of particles, with both short and long range interactions. Figure 5.1 shows some examples of distributions that we shall consider. Currently, for such diverse geometries exhibiting long-range interactions, which naturally appear e.g. in cosmology or turbulent flows in physics and atmosphere science, no model has been proposed in the literature on point processes. To address this problem, we shall introduce a statistical model inspired from the maximum-entropy principle ( [START_REF] Jaynes | Information theory and statistical mechanics[END_REF], cf. Section 2.10), to approximate such point process distributions and simulate new realizations. Recall from Section 2.10, maximum entropy models are based on the description of the distribution with a set of moments, or their estimators. Intuitively, this means that the model is 'as random as possible' under certain constraints, based on the information captured by the moments.

Generation of point processes with maximum entropy models

When defining a generative model based on the maximum entropy principle, there are three underlying problems:

1. Choosing the moments that will describe the distribution. They should be informative enough to capture the geometric structures characterizing the distribution. On the other hand, they should be accurately estimated from a single observation, so the number of moments should not be too large.

2. Specifying a model deriving from these moments. This can be done by defining a maximum entropy model such as the macro-canonical model (maximizing the entropy under expectation constraints), or the micro-canonical model (maximizing the entropy under path-wise constraints).

3. Generating new samples from the model. In the micro-canonical setup, this can be done by minimizing an energy, that defines the set of admissible realizations. The minimization method must make it possible to generate diverse low energy samples without being too costly in terms of calculation.

Recall from Section 2.10, under some technical assumptions (in particular, having density with respect to the reference Poisson distribution), the solution of the macrocanonical model is given by the Gibbs point process. Sampling from the macro-canonical model is usually computationally very expensive (c.f. [START_REF] Bruna | Multiscale sparse microcanonical models[END_REF]). Therefore, we shall place our model in the micro-canonical setup, as detailed in Section 5.2. The main challenges then reside in the problems 1 and 3.

In Section 5.3, we present our method to address the problem of generating new samples: we minimize the energy of a new sample by moving the particles of an initial random configuration using the gradient of its energy with respect to the particles positions. In the point process literature, a classical method (see [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF]) consists in updating an initial random configuration by successively replacing the particles one by one, with new particles located at random positions (we shall call this method random search). The major drawback of this method is its computational cost, as the optimization, which does not use gradient information to minimize the energy, requires a large number of energy evaluations. In fact, this method has been applied to generate point processes formed by a few hundred particles. On the other hand, advanced methods in the modelling of textures and non-Gaussian stationary processes allow for fast sampling by first drawing from an initial Gaussian distribution, and minimizing an energy by gradient descent on the amplitudes of the pixels of the image (see e.g. [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF][START_REF] Gatys | Texture synthesis using convolutional neural networks[END_REF][START_REF] Bruna | Multiscale sparse microcanonical models[END_REF][START_REF] Zhang | Maximum entropy models from phase harmonic covariances[END_REF]). Our approach leverages the efficiency of this sampling method, while ensuring that the resulting samples are atomic measures. The idea of moving the points according to their gradient is often used in molecular dynamics (e.g. [START_REF] Zhang | Ground states of stealthy hyperuniform potentials: I. entropically favored configurations[END_REF][START_REF] Zhang | Ground states of stealthy hyperuniform potentials. ii. stacked-slider phases[END_REF]), however it requires knowledge of the physical mechanisms behind the underlying process. Our statistical modeling approach has a potential to simulate new, complex particle configurations directly from one observation, when the underlying physical phenomena are very complicated to model. This brings us to the other challenge that we address in this work: choosing the moments that we shall use to characterize the distribution. We shall use the wavelet phase harmonics (WPH) covariance moments, presented in Section 3.5, and briefly reviewed in Section 5.4. Recall that, similarly to scattering moments (cf. Section 3.2, and used in Section 4.1.3), they are spatial statistics based on coefficients computed from a wavelet transform of atomic measures. It is known that the covariance between the wavelet coefficients capture only second-order correlations, which are equivalent to the Bartlett spectrum (cf. eq. (3.29) and Section 2.8). To capture information beyond second-order correlations, we apply a non-linear phase harmonic operator to the wavelet coefficients. This operator acts on the complex phase of the wavelet coefficients, without changing their amplitude. The covariance between the resulting coefficients allows one to capture particles interactions across different scales.

Compared to high-order correlation functions, (see [START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF]Section 12.4.2]), our moments have the potential to define a sufficient set of statistics, while maintaining a small estimation error, which is similar to the second-order statistics. Other statistics often used in the point process literature (cf. Section 2.9), such as the k nearest neighbour distribution function suggested in [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF], have a number of elements that grows with the intensity. Since there is only one observation, the number of moments should be limited, in order to control their estimation variance. The wavelet transform allows for direct control over the scales of the structures that we wish to capture, regardless of the intensity of the process. This property allows one to model point processes formed by a large number of particles with a limited number of moments.

The wavelet phase harmonic covariance descriptors are defined as spatial averages evaluated over a point process realization. In practice, the calculation of such descriptors can be done by discretization of the observation window in the form of a grid of pixels. However, making these descriptors differentiable with respect to the positions of the particles remains a challenge. In this regard, we describe in Section 5.5 a complete numerical scheme allowing one to solve this problem. It is based on a differentiable discretization of atomic measures. We further present a multi-scale optimization in the gradient descent, intended to avoid unwanted shallow minima of the energy.

In Section 5.6, we evaluate our model on some distributions exhibiting various geometric structures, like Cox point processes on the edges of Poisson-Voronoi tessellations and on the Boolean model with circular grains. Other processes we consider are Matern hard-core and cluster processes driven by Poisson processes with turbulent intensities. Their intensities are sampled from a turbulent flow simulated from Navier-Stokes equations [START_REF] Schneider | Coherent vortex extraction and simulation of 2d isotropic turbulence[END_REF]. Besides the visual inspection of the samples from our generative model, we evaluate second order correlations and compare the persistent homology diagrams, that has been proven useful for topological data analysis (see e.g. [START_REF] Chazal | An introduction to topological data analysis: fundamental and practical aspects for data scientists[END_REF]).

To evaluate the possible benefits of our methods relatively to the well established model of [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF], Section 5.7 presents a comparison between the two methods, in terms of speed of convergence of the optimization methods, and of quality of the models.

For simplicity, we present our model with strong assumptions about the data (stationarity, periodicity, square window). In practical cases, these assumptions would probably not be fully met, so we discuss in Section 5.9 how to relax these assumptions.

Finally, in Section 5.8, we propose a modified version of our model, that consists in adding a regularization term intended to improve the diversity of the generated samples. Moreover, we give a proof of a result stating the translation invariance of our model (see Section 5.3.2).

Micro-canonical framework for point processes

This section details the specific context of this study, and the definitions that we shall need.

We consider here stationary, ergodic point processes on E = R 2 (cf. Sections 2.5 and 2.7). Recall that a point process Φ is stationary if its distribution is invariant to translation, i.e.

L StΦ = L Φ ,
where L Φ denotes the distribution of Φ, and S t Φ(B) := Φ(B + t), for all B ∈ B, and all t ∈ R 2 . Note that S t is a simple example of a push-forward operator on Φ, which is generally noted F # µ for a point measure µ and a measurable function

F : R 2 -→ R 2 ,
and is simply the displacement of its atoms by the function F

F # µ = i δ F (u i ) .
Recall also that Φ is said to be ergodic if the empirical averages of integrable functions f , over windows of increasing sizes

W s = [-s, s[×[-s, s[, converge almost surely to the expectation of f , i.e. lim s→∞ 1 |W s | Ws f (S u Φ) du = E[f (Φ)] = M f (µ) L Φ (dµ), P-a.s. (5.1)
For a given s > 0, we denote by M s the set of counting measures on W s , and M s its induced σ-algebra. We will consider W s with addition and scalar multiplication modulo W s . Also we shall denote by Su the corresponding shift operator on M s with torus correction on the window W s .

For a point process Φ on R 2 , one can only observe realizations of Φ on bounded subsets of R 2 . For the remainder of this chapter, we shall consider realizations of point processes observed on a finite square window W s = [-s, s] 2 , for some s > 0. We denote by Φ the restriction of Φ to W s , that is Φ is a point process on W s such that, ∀ n ∈ N, ∀ (B 1 , ..., B n ) ∈ B(W s ) n , (Φ(B 1 ), ..., Φ(B n )) = ( Φ(B 1 ), ..., Φ(B n )) in distribution (where B(W s ) stands for the Borel σ-algebra on W s ). A realization of Φ observed on W s is therefore a realization of Φ, and will be noted φ.

As mentioned in Section 5.1, we shall place ourselves in the context of maximum entropy models (cf. Section 2.10). Such models are based on the following intuitive idea: given an observation pattern, we aim at finding new patterns that are similar to, but different from the observation. To this end, we define a notion of similarity by choosing a set of statistics that will be computed on the observation and on the new patterns. The two will be considered similar if their statistics match. Furthermore, if the chosen statistics describe sufficiently well the point process behind our observation, we do not want to add any more constraints, that is, we want to find new patterns 'as random as possible', under the constraints defined by the statistics. This can be formalized by maximizing the entropy of the model. They are used in large classes of stochastic models [START_REF] Geman | Stochastic relaxation, gibbs distributions, and the bayesian restoration of images[END_REF], and will inspire our particle gradient descent model. In this work, we place ourselves in the micro-canonical setting.

Consider a mapping

K : M s → C d ,
for some d ≥ 1, that shall define our set of statistics (e.g. statistics presented in Section 2.9, or wavelet-based statistics in Chapter 3).

Micro-canonical model

The micro-canonical model is defined by replacing the expectation constraints E(K(Ξ)) = a with pathwise constraints. Let φ ∈ M s be our observation sample, of unknown distribution. For all µ ∈ M s , we define the energy of µ as:

E φ(µ) := 1 2 |K(µ) -K( φ)| 2 .
(5.

2)

The micro-canonical set of level , for some > 0, is defined as

Ω := {µ ∈ M s : E φ(µ) ≤ }. (5.
3)

The micro-canonical model is defined as the distribution L that minimizes the KL divergence with respect to the reference distribution L 0 under pathwise constraints requiring L to be supported on Ω :

arg min L KL(L, L 0 ) (5.4)
given

M s 1(µ ∈ Ω )L(dµ) = 1, (5.5)
where 1(•) is the indicator function. As stated in Proposition 2.10.5, if L 0 (Ω ) > 0, the solution to this problem (5.4) (5.5), is the measure L having a uniform density on M s given by dL dL 0

(µ) = 1 L 0 (Ω ) 1(µ ∈ Ω ), L 0 -a.s.
In order to consider the micro-canonical model as a good approximation of the observation distribution, one usually aims at finding K satisfying the following properties: A natural framework allowing one to address (P1) and (P2) consists in defining the descriptors K = (K 1 , . . . , K d ) as a vector of empirical averages

K i (µ) = 1 |W s | Ws f i ( Su µ) du µ ∈ M s , ( 5.6)
for a sufficiently rich class of functions f i on M s , and relying on the ergodic assumption (5.1) regarding Φ. These properties are needed in order to have a model that reproduces typical geometric structures in Φ, and generates diverse samples.

In this study, we shall consider that the number of points of our model in W s is fixed. In such case it is customary to take homogeneous Poisson point process distribution conditioned to have exactly n points in W s as the reference measure, which is equivalent to n points sampled uniformly, independently in W s . We will note this distribution L n 0 . Sampling from the uniform density on Ω efficiently remains an open problem. In the literature of stochastic process modelling, most sampling algorithms rely on the following method: one first samples from an initial, high-entropy measure, and iteratively minimize the energy (cf. (5.2)) of this sample until it reaches Ω . By choosing a high entropy initial measure, one hopes that the resulting model also has a high entropy. Recall that the micro-canonical model has the highest entropy supported on Ω . Contrary to the classical method in the point process literature [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF][START_REF] Koňasová | Stochastic reconstruction for inhomogeneous point patterns[END_REF], which relies on random search, popular methods in image modelling use gradient descent to perform fast sampling in the micro-canonical set. However, optimizing the values of the image pixels does not guarantee that the resulting sample is an atomic measure. For these reasons, in what follows we propose a model based on the transport of a Poisson point process via a gradient descent algorithm.

Particle gradient descent model

In Section 5.3.1, we introduce the particle gradient descent model, that uses gradient descent on the positions of the particles of the sample. This model consists in using the gradient of a prescribed energy to move the points of an initial random configuration, until we obtain a pattern similar (in an informal sense) to the observation. We then present in Section 5.3.2 a theorem stating that this model preserves some basic invariances of the original distribution. This result, allowing us to gain some understanding about the entropy of our model, extends the results of [START_REF] Bruna | Multiscale sparse microcanonical models[END_REF]. While we focus here on planar point processes, our approach can readily be extended to any dimensions.

Particle gradient descent model

As in Section 5.2, let φ ∈ M s be our observation sample of unknown distribution, and K : M s -→ C d , for some d ≥ 1, a mapping defining our descriptors. We note the resulting energy E φ (cf. eq. ( 5.2)).

Let φ0 sample from an initial distribution that we choose as L φ(Ws) 0 (i.e. the same number of particles as φ, drawn uniformly and i.i.d.). We minimize the energy of φ0 through its gradient with respect to the particles positions. More precisely, we define the mapping F :

M s -→ M s µ = i δ u i -→ i δ u i -ζ∇u i E φ(µ) (5.7)
for some gradient step ζ > 0. The measure F (µ) can be seen as the push-forward F µ# µ of the measure µ by the mapping F µ (u) := u -ζ∇ u E φ(µ). Note that the function F µ depends on the measure µ which is pushed forward. For any initial point measure φ0 ∈ M s we define the successive point measures:

φn := F φn-1 # φn-1 , n ≥ 1.
(5.8)

Pushforward of the point process distributions

The pushforward operation F #µ on M s induces the corresponding pushforward operation on the probability measures on M s , which are distributions of point processes. We denote this latter by F # : For a probability law L on M s F # L(Γ) := L({µ ∈ M s : F µ# µ ∈ Γ}), for any Γ ∈ M s . Then, for an initial probability law L φ0 on M s we define the successive probability laws

L Φn := F # L Φn-1 , n ≥ 1.
(5.9)

Note that L Φn is the distribution of the point process Φn obtained by n iterations of (5.8) starting from Φ0 having law L Φ0 = L Φ(Ws) 0

. Our model is defined by setting a fixed number of iterations as a stopping rule. Figure 5.2 illustrates the particle gradient descent method to generate samples in Ω . Observe that our model takes inspiration from the micro-canonical model, however there is no guarantee that our model reaches Ω (defined in (5.3)), for any > 0. By setting a fixed number of iterations and not rejecting any configuration, we make the implicit assumption that our model reaches a low energy level.

Leveraging invariances

One can leverage some a priori known invariance properties of Φ (for instance stationarity or isotropy), by building a model that satisfies the same invariance properties as Φ. By using the descriptor K with the same invariance, the particle gradient model respects these invariance properties. In particular, we obtain a stationary point process model when K is defined by the empirical averaging (5.6).

This requires some explanation, since invariance properties of the distribution of Φ do not, in general, imply any natural invariance of its restriction Φ to W s . Indeed, while some invariances can be observed on the torus for the distribution of Φ on R 2 (the most popular being translation invariance), it does not imply the same for Φ with respect to the translation on W s . The latter, that we call circular stationarity, requires also Φ to be periodic. However, circular stationarity of the generated point process on large window W s (as a distributional approximation of Φ) can be considered as a desirable ersatz of the stationarity of Φ. Indeed, in what follows we shall formulate a result saying that, when K and the distribution of Φ0 are invariant with respect to some subset of rigid circular transformations on W s , then the resulting model satisfies this property as well.

More specifically, a rigid circluar transformation on W s is an invertible operator T on W s of the form T u := Au + u 0 for some orthogonal matrix A with entries in {-1, 0, 1} and u 0 ∈ W s . Note that the matrix A is restricted in integer entries for T to be a well defined invertible operator. It encapsulates translations, flips, and orthogonal rotations.

We say that:

• The initial probability law L Φ0 of the model is invariant to the action of

T if ∀ Γ ∈ M s , L Φ0 (T -1 # (Γ)) = L Φ0 (Γ). • The descriptor K is invariant to the action of T if ∀ µ ∈ M s , K(T # µ) = K(µ).
Theorem 5.3.1. Let T be a rigid circular transformation. Let Φ0 be a point process on W s such that its distribution L Φ0 is invariant to the action of T and let K be a descriptor invariant to the action of T . Then, for all n ∈ N, L Φn defined as the push-forward of L Φ0 by (5.9) is invariant to the action of T . This result also holds true for the distributions of our model, defined by section 5.3.1.

The proof is given in Section 5.8, for a more general model. The result itself is inspired from [START_REF] Bruna | Multiscale sparse microcanonical models[END_REF], where the preservation of invariance is proven for the gradient descent model in the pixel domain. Observe, the invariance of the distribution of the point process Φn increases the diversity of the generative model samples. Any descriptor taking the form of eq. (5.6), using the periodic shift operator Su , is invariant to all circular translations.

A drawback is that such a boundary condition introduces a statistical bias to the spatial average (5.6) as an estimator of E[K( Φ)] in the case of a non periodic Φ over W s . One can expect, however, that when the window size is large enough and spatial correlations of the patterns are not too large, this border effect becomes negligible.

Wavelet phase harmonic descriptors

Classical descriptors for spatial point process usually include statistics more or less directly related to the pair correlation function, such as Ripley's K-function, Besag's L-function, or the radial distribution function (see [START_REF] Chiu | Stochastic geometry and its applications[END_REF]Section 4.5]). All of these functions only capture second order correlations of the process. Other usual functions are the empty space function or the k-nearest neighbors function; (see Section 2.9, and [37, Section 2.3.4 and 4.1.7]). In [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF], the authors advocate the use of the k-nearest neighbors distribution function, with a k significantly greater than 1. In addition to being non-differentiable, these moments suffer from another drawback. If one wants to capture geometric structures formed by the particles, up to a fixed scale, the number of moments (i.e. the k nearest neighbours) will grow linearly with the number of particles forming such structures. This can become a problem if the intensity of the process is large, both computationally, and from a statistical point of view, as the variance of the moments may become large when estimated from a single observation.

For this reason, we choose to use descriptors for which the spatial range of structure captured is independent of the intensity of the process, and the computational time is linear in the number of points. As a result, this method would become much faster for large samples, as the number of statistics would remain constant. In this work, we shall use the wavelet phase harmonics (WPH) covariance moments (cf. eq. (3.32)) to represent the distributions of our point processes. Recall from Section 3.5, these moments are based on the wavelet transform of a point process (eq. (3.5)), and in this study we shall used periodic wavelets for the observation window W s (cf. eq. (3.33)). The WPH covariance moments use a non-linear operator to adjust the phase of the wavelet coefficients, in order to capture dependencies between coefficients at different scales and orientations.

More precisely, to capture edge-like structures in a planar point process, we shall use bump steerable wavelets. Recall from Example 3.1.3, they are defined by the translations, dilations and rotations of a complex analytic function ψ(u) ∈ C with compact frequency support. The finite wavelet family is defined as {ψ ξ } ξ∈Ξ , for ξ = (j, θ), with 0 < j < J, and J such that all wavelets have (approximate) support contained in the observation window W s . Scales equal or larger than J are carried by a low-pass filter ψ 0 whose frequency support is centered at 0. The angles are chosen such that θ = 2π L , with ∈ {0, • • • , L -1}. For any ξ, ξ ∈ Ξ, k, k ∈ Z, τ ∈ W s , we note γ = (ξ, ξ , k, k , τ ), and estimate

K γ Φ = Cov [ Φ ψ ξ ] k (u), [ Φ ψ ξ ] k (u -τ ) ,
(5.10) which does not depend on u ∈ W s , as the process is (circular) stationary. Recall that the phase harmonics operator [•] k is defined, for all z ∈ C, as

[z] k = |z|e ikϕ(z) .
As explained in Section 3.5, the overlap between the frequency support of [ Φ ψ ξ ] k and that of [ Φ ψ ξ ] k is necessary for the wavelet phase harmonic covariance to be large. Due to the frequency transposition property of the wavelet phase harmonics, it is empirically verified that the covariance at kξ k ξ is often non-negligible when the process is non-Gaussian (i.e. has structures beyond second order correlations). We shall follow this empirical rule to define a set of WPH coefficients to correlate.

Let

v ξ,k = E([ Φ ψ ξ (u)] k ).
We define the descriptors K(µ) using (5.6), as empirical estimators of moments. Additionally, let us denote µ ξ,k (u) := [µ ψ ξ (u)] k -v ξ,k . Taking the spatial average (5.6) gives the descriptor of the form:

K(µ) = 1 |W s | Ws µ ξ,k (u) µ ξ ,k (u -τ ) * du (ξ,k,ξ ,k ,τ )∈Υ . (5.11)
As Φ is circular-stationary, it suffices to choose τ = 0 to measure spatial shifts (i.e., there is no need to consider u -τ , for some non null τ ). Note also that K(µ) is invariant with respect to any circular translation Su of µ ∈ M s on u ∈ W s .

In the numerical computation, we shall replace v ξ,k in (5.11) by vξ,k = 1 |Ws| Ws [ φ ψ ξ (u)] k du as a plug-in estimator for the first-order moment v ξ,k . The K( φ) modified in this way becomes an empirical estimator of the covariances in eq. (5.10). This is a good approximation of K(µ) as the estimation variance of the covariance moments is typically much larger than that of the first-order moments.

Note that the wavelet phase harmonics do not increase the amplitude of the wavelet coefficients with k > 1. This approach may thus significantly reduce the variance of the descriptor K (to satisfy (P1)) compared to the higher order correlations, while still capturing information beyond second-order correlations (to satisfy (P2)). (5.11) Rather than detailing the full list of elements in Υ, we provide an intuitive way to choose the set Υ. For the full list, see [START_REF] Brochard | Particle gradient descent model for point process generation[END_REF]Section 4]. Overall, the total number of elements in Υ is in the order of O(L 2 J 2 ). Note that the smallest structures that the descriptors can capture depend on the spatial support C of the wavelet ψ. Information about structures smaller than C can be added in a post-processing step that will be explained in Section 5.

Choice of the covariance set Υ in

• Choice of J: The covariance set Υ depends on the parameter J, which is the maximal scale of the wavelet transform. A suitable choice for this parameter J would be one allowing for a good trade-off between satisfying the sufficiency of K, while maintaining the concentration property (cf. properties (P1) and (P2) from Section 2.10).

• The parameter τ is chosen so that each wavelet is translated in a particular direction in order to capture correlations along nearby edges in the observation.

• Choice of (ξ, ξ , k, k ): These parameters are chosen in order to capture 2nd-order correlations, as well as dependencies between wavelet coefficients at different scales and orientations, both with and without phase information, based on a rule of thumb that kξ ≈ k ξ due to the frequency transposition property of the wavelet phase harmonics.

Numerical scheme for particle gradient descent

Calculating the wavelet phase harmonic covariances can be computationally demanding (due to the calculation of two integrals). In order to gain some efficiency, we can perform the computations in a discrete domain. However, the energy needs to remain differentiable with respect to the positions of the points in the pattern. We propose a method, consisting of a Gaussian smoothing of the configuration of points, to address this problem. Building on that method, we then present two technical aspects of the sampling method.

In this section, we discuss a complete numerical scheme to generate samples from the particle gradient descent model, defined with the wavelet phase harmonic descriptors presented in Section 5.4. It is composed of the following ideas:

• Discretization for an approximate calculation of the covariance of the wavelet phase harmonics: necessary to accelerate the calculation of the descriptor and the gradients.

• Multiscale optimization: allowing one to avoid shallow local minima in the gradient descent model. At each scale, we use a quasi-Newton gradient-descent method for greater efficiency.

• Final blurring (optional): to add a priori information on structures whose size is smaller than C into the model samples. It helps to get rid of some clusterisation (clumping) artifact caused by the initial discretization.

Discretization

Differentiable discretization of atomic measures

To compute the descriptor K in (5.11) for a point measure µ, we need to integrate functions over the observation window W s (first for the convolution operators, then for the averages). Computationally efficient integration requires discretization of the atomic measure. The main difficulty is to do it in such a way that the (periodic) convolutions of the discretized atomic measures with wavelets, as in (6.2), remain differentiable with respect to the positions of the original atoms in µ, so that we can still perform gradient descent. Classical finite element methods may not achieve this goal efficiently.

We are going to approximate our atomic measures on W s by matrices (images) of given size N × N (the image resolution), and then use the automatic differentiation software Pytorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] to perform the following operations. It allows one to compute the derivative of a modified energy w.r.t. any point u i in µ. The following paragraph details this discretization:

We first map a given point measure µ on W s to a continuous function µ σ by the convolution

µ σ (u) := µ g σ (u) = u i ∈µ e -|u-u i | 2 2σ , x ∈ W s (5.12)
with a (periodized) Gaussian function g σ of given standard deviation σ. Then we evaluate µ σ on the N × N regular grid inside W s and denote the resulting matrix µ N σ , with entries called (values of) pixels. The convolution with a Gaussian function makes each entry of µ N σ smoothly depend on the atom positions of µ. We then compute K(µ N σ ) instead of K(µ), where K is this discrete analogy of the descriptor (5.11) (cf. [120]1 ). Note that, because the value of a pixel continuously depends on the positions of the atoms, this discretization makes our descriptor only invariant to discrete translations (multiple of the pixel size 2 s N ), for which Theorem 5.3.1 applies. The gradient of the energy

| K(µ N σ ) -K( φN σ )| 2
with respect to each atom position of µ can therefore be computed using automatic differentiation (with the Pytorch software [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF]). Indeed, we know that K is differentiable w.r.t. each entry of µ N σ , as a combination of linear and non-linear operators. Moreover, for any i, j ∈ {1, N } 2 , noting ĩ = -s + 2si/N, j = -s + 2sj/N , (5.12) gives us that

µ N σ (i, j) = µ g σ ( ĩ, j) = u i ∈µ e -|( ĩ, j)-u i | 2 2σ , ( 5.13)
which is differentiable w.r.t. any u i in µ. This discretization step is illustrated in the Figure 5. In signal processing, the Gaussian function acts as a low-pass filter. It is needed to cut-off high frequency information of µ so that µ σ can be discretized into an image with negligible aliasing effect. This means that µ σ carries the information on the positions of µ up to some precision which depends on σ. The subsequent evaluation of µ σ on the grid N × N in µ N σ implies that σ cannot be taken too small. Indeed, we take σ min = s N as the lowest value of σ.

Wavelet discretization and choice of scales

As stated in Section 5.4, the family of wavelets used in our descriptor is constructed by dilating the mother wavelet ψ in the range of the scales 0 ≤ j < J. Based on the choice of N , we set C = 2s N . In this way, the spatial support of ψ has a radius C of one pixel of the image. As a consequence, this smallest-scale wavelet ψ can also be discretized (without significant aliasing) in order to compute the discretized descriptor K.

The choice of the largest scale J can be decided based on the visual structures in the observation. For example, if we want to model structures whose spatial size is close to the size of the window [0, 1/8] 2 ⊂ W s , we shall set 2 J C = 1/8, i.e. J = log 2 (N ) -3 -log 2 (2s).

Multiscale optimization

Phase harmonic covariance moments of point process images (i.e. point patterns converted into regular pixel grids, as described above) may have large values at high frequencies (large values of |ξ|, |ξ | in (5.11)), due to the fact that the point-images are composed of local spikes when σ is small. This implies that these high frequency statistics have an important impact on the gradient of K, which in turn can lead to the gradient descent model being trapped at shallow local minima, where only the high frequencies are well optimized to match the observation.

This optimization issue can be overcome by matching the descriptors from low frequency to high frequency in a sequential order, through an appropriate modulation of the parameter σ ∈ σ 0 , σ 1 , ..., σ J-1 of the Gaussian functions used to discretize µ, introduced in Section 5.5.1.

Indeed, since Gaussian functions are low-pass filters, we can interpret the convolution in (5.12) as a blurring, limiting the space localization of Dirac measures. When such smoothing of the point pattern is done by a Gaussian function that has a large σ, the high frequencies of the signal function are close to 0 and the same holds true for the phase harmonics, because wavelets are localized in frequency. Therefore the wavelet phase harmonics are dominated by the low frequencies. Thus, by smoothing the observed sample and generating the optimal one with high variance Gaussian function, we create a new objective leading, in the gradient descent optimization, to a point configuration for which only low frequencies moments (small values of |ξ|, |ξ | in (5.11)) are matched with the ones of our observed sample. Thus, we propose a multiscale gradient descent procedure that consists in choosing first a high value for precision parameter σ, run the optimization algorithm, and then reduce the value of σ to run the optimization again, starting from the result of the previous run (and repeat this operation until σ = σ J-1 ). We choose σ j := s N 2 J-j-2 . Note that σ J-1 is equal to σ min . For numerical efficiency, we perform the gradient descent procedure using the L-BFGS optimization algorithm [START_REF] Liu | On the limited memory bfgs method for large scale optimization[END_REF].

Empirical evidence in Section 5.6.2 shows that the above multi-scale optimization procedure allows one to reconstruct (modulo translation) the observed sample when using K defined with J = log 2 (N ) -2, which is not the case when simultaneously optimizing all frequencies. In order to preserve the ability to reproduce geometric structures at all scales, we shall also apply this multiscale optimization method to our model defined with

Final blurring

We observed that the contrast between the continuous nature of our objects and the discrete approximation described in Section 5.5.1 creates undesired artificial structures at frequencies higher than the image resolution: when the number of particles in a configuration is large with respect to the number of pixels in the image, or if the configuration exhibits strong clustering behaviour, several pixels may contain more than one particle. In such cases, our algorithm produces samples having an artificial clustering structure inside each of these pixels (see [START_REF] Brochard | Particle gradient descent model for point process generation[END_REF]Section 5] for an illustration of this phenomenon).

To remove this artificial clustering, we chose to force these high frequencies to be "as random as possible", i.e. to have Poisson-like structure. To this end, we introduce a uniform i.i.d. perturbation of the positions of points after the last optimization run. It can be viewed as an additional, this time stochastic, measure transport, following the deterministic one from the particle gradient descent. This final randomization can be viewed as enforcing a-priori information on high frequency structures of the process: Poisson-like structure.

Numerical experiments

In this section we present numerical experiments involving our generative model. We begin by presenting in Section 5.6.1 our numerical settings, in particular the distributions of point processes whose samples are used as original point patterns. We next evaluate how well our generative model with the phase harmonic covariance descriptor can generate samples similar to those given by the original point processes. In Section 5.6.2, we evaluate these models by comparing samples from the original distributions to samples from our models, visually as well as by estimating their power spectrum (defined as the average over angles of the squared modulus of the Discrete Fourier Transform, cf. Section 2.8). The power spectrum gives information equivalent to the second order correlation function of the process (cf. [START_REF] Brémaud | Mathematical Principles of Signal Processing: Fourier and Wavelet Analysis[END_REF]), which captures clustering or repulsive behaviour between atoms of a realization. Such information cannot always be detected visually. In order to further quantify how well our model captures visual geometric structures, and to gain some insight into the ability of our model to produce diverse samples, we shall use the topological data analysis (TDA), derived from the theory of persistent homology. The comparison will be done in Section 5.6.3.

Numerical settings

We first describe the original point processes on which we shall evaluate the particle gradient-descent model, then specify the parameters of the model in the numerical experiments.

Original point process distributions

For our experiments, we choose point process distributions that show complex geometric structures, for which we can visually recognize geometric structures. We begin by presenting results for Cox (double-stochastic Poisson) processes with Poisson points living on one dimensional structures generated by two famous stochastic geometric models, namely edges of the Voronoi tessellation, and the Boolean model with circular grains of fixed radius considered in Example 2.8.4. Both underlying geometric models are generated by a Poisson parent process within the observation window W s , and we construct these models in a periodic way to avoid border effects. We call the respective Cox processes Voronoi and Circle processes. Note that, for these two processes, Poisson points live on different geometric shapes: polygons for the Voronoi and possibly overlapping circles for other one. Additionally, we consider two different radii of circles.

Then, we take interest in distributions having turbulent intensity (derived from the simulations of a decaying isotropic turbulent vorticity field driven by 2d Navier-Stokes equations, see e.g. [START_REF] Schneider | Coherent vortex extraction and simulation of 2d isotropic turbulence[END_REF]). Such fields exhibit complex mulsticale structures, and are representations of physical phenomena, known to be difficult to model faithfully. Furthermore, the distributions we consider have much greater intensities that the previous Cox models. From the turbulent intensity, we sample three different processes, exhibiting distinct microscopic structures (repulsive, independent or clustering): a Matern cluster process, a Poisson point process and a Matern II hard-core process, see [37, Example 5.5 and Section 5.4, respectively]. We study the ability of our model to reproduce simultaneously the macroscopic (i.e. the turbulent intensity) and microscopic structures (i.e. at small scales) of the process.

The number of points in the Cox Voronoi, Small circles, Big circles, and the Turbulent Hardcore, Poisson, ad Cluster processes are around, respectively, 1 900, 2 500, 2 000, 1 700, 3 800 and 13 000. Note that, for comparison, point patterns considered in [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF] have around 400 points.

Image resolution, number of iterations and computation time

As discussed in Section 5.5.1, point configurations are convoluted with Gaussian densities and evaluated on N × N grid (as images) in order to efficiently compute our descriptors, and move the particles with gradient descent. For simplicity, we fix s = 1/2 for all the examples that we shall consider. The ultimate Gaussian variance (precision) of this mapping is thus σ min = 1 2N . The larger N is, the more information we are able to keep (in high frequencies), but the larger the computation time. We chose for our experiments a resolution of N = 128. We show one example where a higher resolution, N = 256, is used to capture most of the high frequency information. The number of iterations of the L-BFGS optimization is chosen to be 100 for each scale σ j (a total of 400 iterations for N = 128, and 500 iterations for N = 256). The number of angles in the steerable wavelets is L = 8. The average computation time on 4 GPU (Nvidia Tesla P100) for a sample for a turbulent process having roughly 13 000 points with resolution N = 256 is between 5 and 10 minutes while the same task at the resolution N = 128 takes between 1 and 2 minutes.

Visual evaluation and spectrum comparison

We evaluate the ability of our model to capture and reproduce geometric structures exhibited by realizations of the point processes described in Section 5.6.1. A natural first method to assess the sufficiency of a generative model (property (P2)) is visual evaluation, which is widely used in image analysis but subjective. We then compare the power spectra of our models and the original distributions. To estimate the power spectra, we generate (for each original distribution) 10 i.i.d. samples from the same model (i.e. from the same observation sample φ, but with different initial configurations φ0 ). We average the power spectra of the 10 syntheses, and compare it to the average of 10 i.i.d. samples from the original distribution. All these samples will also serve in Section 5.6.3 to compare their geometric similarities. Figure 5.4 shows a study of our three Cox distributions. The first line presents samples from the original distributions. The second line presents samples from the model using our descriptor with J = log(N ) -2. In this setup, the concentration property is not satisfied (see property (P1) in Section 5.2), and the result is the memorization of transform, which is included in the descriptor K(µ) (c.f. (5.11)). Moreover, the variance of the empirical information at small k can create extra error since it can be far away from its expectation. Similarly, because the wavelet convolutions average the spectral information over different frequency bands when using a reduced number of τ in K(µ), our descriptor does not capture fast oscillations in the power spectra in the range of k ≤ N/2 = 64 (see [START_REF] Zhang | Maximum entropy models from phase harmonic covariances[END_REF] for more details about how to capture these oscillations). This is observed in the cases of Small and Big Circles Cox processes. See example 2.8.4 for a theoretical formula of the power spectrum in the case of (Small or Big) circles. When k > N/2 = 64, the descriptor K (cf. 5.5.1) does not contain accurate spectral information due to a finite image resolution N . In this regime, we observe a smooth decay of the (log) power spectrum towards 0. We observed that if we apply the final blurring (cf. Section 5.5.3), then the error of the model spectrum becomes larger. Therefore, for these three processes, no final blurring has been applied.

All models discussed up to now are Cox processes, with Poisson (hence independent) points sitting on some random macroscopic structures. Figure 5.5 presents our analysis of three turbulent point processes having different microscopic structures: a hardcore, a noncorrelated (Poisson) and a clustering one. We see that our generated samples capture to some extent this microscopic structure. For the clustering model, the presented synthesis is done with N = 256. We see that our model (using J = log(N ) -3) can generate samples with similar macroscopic and microscopic structure. The power spectrum at small k has larger errors, as we have observed in the Cox models. However, since the power spectra are mostly smooth in these Turbulent processes, we observe a relatively small spectrum error over a wide range of k. This is also due to the use of the final blurring which helps to remove some artificial spectrum errors for k > N/2. For the clustering model, we also compare the power spectrum of two models with different resolutions: N = 128 and N = 256. We see that setting a higher resolution reduces significantly the error, allowing to match the spectrum up to k 80. we still observe some small error when k ≥ 80, probably because of the final blurring (cf. Section 5.5.3), which may also impact the high frequencies that we optimize. Overall, both the visual and the spectral analysis suggest that our model can generate well various Turbulent points processes.

Persistent homology and topology analysis

As previously mentioned, power spectrum evaluation corresponds to the comparison of second order moments, which only partially capture geometric structures. Visual evaluation can be more discriminate, but is subjective. To evaluate more precisely the ability of our model to capture the geometric structures of the given distributions, we shall use a representation of objects derived from persistent homology theory, which is a powerful algebraic tool for studying the topological structure of shapes, functions, or in our case point clouds. We shall perform this evaluation by comparing the persistence diagrams of the generated samples to those of the original ones. Furthermore, this representation allows us to evaluate in a simple way the ability of our model to produce diverse samples.

We begin by a brief, intuitive presentation of persistence diagrams, and the whole comparison method that will be simply referred to as topology data analysis (TDA). For more details we refer the reader to [START_REF] Boissonnat | Geometric and topological inference[END_REF]Section 11.5]. We then present the TDA of our point process distributions and models. TDA can be seen as a complementary tool with respect to the spectrum analysis, being more consistent with visual perception (see [29, Appendix C] for more details about this link).

Persistence diagram Persistent homology theory describes a way to encode the topological structure of a point cloud through a representation called persistent diagram (PD). It is constructed, for a given point configuration φ ∈ M s , from the family (G r ) r≥0 of Gilbert graphs, where the vertices are the positions of atoms of φ, and the edges are pairs of points closer to each other than r. (In our case we use the periodic metric.) Then, we fill-in the triangles (triplets of points joined by edges) of the graph. Points, edges and filled-in triangles constitute the so-called 2-skeleton of the Vietoris-Rips (VR) complex. For any r ≥ 0, we study two characteristics of the skeleton: its connected components, and its holes (this latter notion is well formalized in the algebraic topology, in our case they correspond to the natural idea of a hole). Each connected component "is born" at time (radius) r = 0 and it "dies" at some time r > 0 when it is merged with another connected component. Similarly, each hole has a birth time (r > 0) corresponding to the minimal radius at which it appears, and a (larger) death time corresponding to the minimal radius for which the hole is completely filled-in by the triangles. The persistence diagram of φ is the collection of pairs of birth and death times of the connected components and holes. It is hence a point process in the positive orthant of the plane, offering a multiscale (as our wavelet-base descriptor) description of the topology of φ. As our descriptor, it is also stable to small deformations of φ. It is hence interesting to use this alternative tool to evaluate our generative model.

Topological data analysis

Our approach in this matter is inspired by [START_REF] Chazal | An introduction to topological data analysis: fundamental and practical aspects for data scientists[END_REF], and we refer the reader to this paper for a more detailed description. We use the 'holes' birthdeath process, as it appears more relevant to capture information in the Cox distributions, such as the polygons and the circles.

In order to compare the distributions of our models to the original distributions, we compute the PDs of our samples from each distribution (cf. Section 5.6.2 for a description of these samples). Recall, these PDs can be viewed again as point clouds in two dimensions. Therefore, a distance between two PDs can be computed, and we use in this regard a periodic version of the Wasserstein distance between two point clouds on the plane (we found that the bottleneck distance also suggested in [START_REF] Chazal | An introduction to topological data analysis: fundamental and practical aspects for data scientists[END_REF] is not sufficiently discriminating for our point patterns). We obtain in this way a distance matrix between different PDs (reflecting topological similarities or differences of the point processees realizations for which PDs were calculated). We then apply a standard dimension reduction algorithm (namely Multi Dimensional Scaling) to this distance matrix, to represent every PD (and hence the corresponding sample) as one point on the plane, and we visualize the representation of all samples.

TDA of our experiments

In the plots on the first line of Figure 5.6, we study separately the Cox Voronoi and Big circles processes, and the turbulent hardcore process. In each plot, we observe 20 dots (having different shapes), each representing one configuration of points in W s (the term "dot" is used to avoid confusion with points in W s ). For each model there are 10 dots representing i.i.d. realizations of the original distribution and 10 representing realizations from the generative model. For each plot, the sample additionally marked with a black dot represents the observation used in our model to distribution as it is to the original Big circles distribution. Nevertheless, the original distribution of the Big circles and the Small circles are well separated, suggesting that there is a topological distinction that is not well respected in the small circles. However, this kind of error is hard to perceive visually. On the other hand, the Voronoi case is well separated from the other two, suggesting that the model is better than the ones of the circles distributions. For the turbulence case, we observe that the three distributions (both original and model) are well separated, and each respective model is closer to its original distribution than to the other distributions. This agrees with our earlier visual and spectral analysis, suggesting that our model is able to capture complex geometric structures formed by a large amount of points.

Remark: We used the R packages TDAstats [START_REF] Wadhwa | TDAstats: R pipeline for computing persistent homology in topological data analysis[END_REF] to calculate the PDs of our point patterns and TDA [START_REF] Fasy | TDA: statistical tools for topological data analysis[END_REF] to calculate their Wasserstein distances. Due to memory constraints, for second line, the analysis was done using a random thinning to reduce the number of points of each sample to 2 000, which could artificially impact the results. The experiments were repeated several times and the variability in the random thinning did not impact our conclusions.

Numerical comparisons

In this section, in order to illustrate the advantages of the method presented here, we present a brief comparison between the method presented in [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF] and ours.

Differences between the two methods

The method developed in [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF] and ours are based on the same idea: to produce similar but different point patterns to a given observation, one first defines what should be 'similar' between the observation and the synthesis, by choosing a set of statistical constraints, computed on the observation. Then, starting from an initial random configuration of points, one iteratively modifies this configuration in order to match the set of prescribed statistics (by minimizing an energy E, related to the square difference between the statistics of the original and the synthesised point patterns). If the set of statistics does not describe the observation itself, but rather its underlying distribution, then the output of the optimization procedure should be a new point pattern, similar but different to the observation.

However, the two methods differ on two points:

• First, the optimization method to match the set of statistics. In [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF], the optimization steps can be described as follows: given the point configuration being

synthesized φ k = N i=1
δ u i,k at some step k, a point in the configuration is chosen uniformly at random, say u j,k , for j ∈ {1, • • • , N }. A candidate for a new point v ∈ W is chosen uniformly at random in the observation window. Then, if the energy of φk := φ k -δ u j,k + δ v is lower than the energy of φ k , we define φ k+1 = φk . Otherwise, φ k+1 = φ k . We call this optimization random search (RS).

In our method,

φ k+1 = N i=1 δ u i,k -ζ∇u i,k E(φ k )
, where ∇ u i,k E(φ k ) is the gradient of the energy with respect to the point u i,k of φ k , see (5.7). This optimization method will be noted (GD).

• Second, the set of statistical constraints used to describe the geometry of the point patterns. In [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF], the authors use the k nearest neighbour distance distribution functions (d.f.) D k (r), for k ∈ {1, • • • , k max }, k max ≥ 1, and evaluated at a sequence of radii r ∈ {r 0 , • • • , r max }, r max > 0, (see [106, p. 267]). We call this statistical descriptor nearest neighbour distances (NND).

Our statistics are based on the covariance between phase harmonics of the wavelet phase harmonics coefficients (WPH) of the point patterns (see 5.4).

Preliminary discussion

Before presenting a numerical comparison between the two methods, we briefly explain what the limitations of the method in [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF] are, and why our method could overcome such limitations. The optimization method in [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF] is based on random search. This implies that for some configuration φ k at some step k, there may be a lot of failing new candidates to replace some point in φ k before finding one that reduces the energy of φ k . This means that there may be a lot of energy evaluations before updating the current configuration. Furthermore, each iteration, requiring one energy evaluation, only moves one point in the configuration. Conversely, at each iteration, our algorithm computes the energy as well as the gradient of the energy, and all the points are moves according to the gradient. This implies that, for some energy level e > 0, the gradient descent method may reach this level in less iterations than [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF].

Moreover, the statistical constraints used in [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF] are based on 3 parameters: the number of neighbours for the points in the configuration, the maximal radius at which to evaluate whether or not there is a neighbour, and the number of radii between 0 and this maximal radius. While this latter relates to the precision of the d.f.'s, and can be pretty much fixed for any observation, this is not the case for the first two parameters. Indeed, k max and r max have to be chosen so as to describe the geometry formed by the points, up to some scale. The maximal radius r max can be seen as the maximal scale up to which the statistical constraints describe the geometric structure. This can be fixed depending on the observation, but should not be too large, in order to satisfy the ergodic averaging property. However, for a fixed sequence of radii, the parameter k max can change significantly depending on the observation. Even if two configuration exhibit structures up to similar scales, the number of nearest neighbours inside some ball may differ depending on the intensity of the process. For instance, consider the Cox Circles distribution, where points are located on circles of fixed radius r 0 , with the center of those circles forming a Poisson point process. If the observed pattern has around 10 points per circle, one would probably need to fix r max = 2r 0 , and k max = 10. However, if the circles contain an average of 100 points, then one would have to increase k max up to 100, even if the circles have the same size as before. This would increase significantly the number statistics to compute at every step. Conversely, our descriptors only depend on a number of scales at which we compute the wavelet coefficients, which does not depend on the intensity of the process. The only parameters to fix are the maximal scale J, and the minimal scale, set by the resolution N , which relates to the precision set by the number of radii in the k nearest neighbours case. In the above Cox Circles example, if the points have Poisson distribution on the circles, then the size of the descriptor will not change between the two setups (10 or 100 points per circles).

Numerical comparison

In this section we present a numerical comparison between the method of [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF] and ours. This comparison aims at illustrating the two following points:

1. For the same energy, the gradient descent optimization method reaches low energy levels in less time (i.e. less evaluations of energy value) than the optimization method used in [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF]. To highlight this point, we shall consider the Cox Voronoi example, and use the WPH descriptors (c.f. (5.11)) to define the energy, and compare the RS and GD optimization methods.

2. The amount of information captured by the kth nearest neighbours d.f.'s depends on the intensity of the process, regardless of the scales of the structures. Considering the Turbulent Poisson example with different intensities, we shall see that, for a fixed descriptor (i.e. fixed k max , r max , number of r), the quality of the syntheses decreases with the intensity of the process.

These two points are illustrated in the two following sections, respectively. A last section will present an overall comparison of the two methods, on the Cox Voronoi and Turbulence Poisson examples.

Comparison between RS and GD (Cox Voronoi example)

For this experiment, we study the Cox Voronoi example, and define the energy from the wavelet phase harmonics covariances, presented in Section 5.4. Let K be our descriptor (defined in (5.11)), φ our observation sample, and E φ(•) = 1 2 |K(•) -K( φ)| 2 | the corresponding energy. We define the relative energy by

e(•) = 2 E Φ(•) |K( φ)| 2 = |K(•) -K( φ)| 2 |K( φ)| 2 .
(5.14)

We ran the optimization of the energy with the random search method from [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF], and observe the relative energy of the syntheses (for 10 syntheses), after n = 19870 and n = 29805 iterations, i.e. respectively 10 and 15 iterations per point. After n = 19870 iterations, the algorithm reaches a relative energy of e = 9, 00.10 -4 (with a std of 9, 23.10 -5 ), and after n = 29805 iterations, we found e = 4, 76.10 -4 (std= 2, 47.10 -5 ). It took an average of 1h04min and 1h36min respectively. We observed the respective relative energies, and ran our gradient descent optimization algorithm (without the multi-scale procedure) until the relative energy reaches the levels from the random search method. The results and comparisons with our method are summarized in Table 5.1.

Dependence of WPH and NND on the intensity of points (Turbulence Poisson example)

To illustrate our second point, for NND descriptor we fix the parameters of the k th nearest neighbours d.f.'s to k max = 16, r max = .125 (on a window of size 1), and discretize D k (r), r ∈ (0, r max ], regularly by 250 values of radii r. With this fixed descriptor, we perform a synthesis using RS optimization for three different observations: the Turbulence Poisson observation randomly thinned to have 500 points, the same observation thinned Random search Gradient descent e = 9, 00.10 -4 19870 (1h04m) 52 (0m35s) e = 4, 76.10 -4 29805 (1h36m) 69 (0m45s) The time per iteration in the gradient descent method is larger, due to the possible several energy (and gradient) evaluations for the line search. However, the total amount of time is much lower. The computations have been run on a single GPU Nvidia Tesla P100.

to 2000 points, and the raw observation, which contains 3784 points. We set the number of iterations to 400 iterations per point. Figure 5.7 shows examples of syntheses for the 3 different patterns, as well as syntheses from our method using WPH descriptor with GD multiscale optimization. We observe that the method using the RS+NND fails to reproduce the geometric structures in the example with the largest number of points. Figure 5.7 also presents the estimations of the spherical contact distribution function (SCDF), defined for a point process Ξ as H s (r) := 1 -P(Ξ ∩ B(0, r) = ∅), where B(0, r) denotes the ball of radius r, centered at 0. For the observation with the largest number of points, the model from [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF] presents a larger error than ours.

Direct comparison

On the two distributions considered above, we shall illustrate the overall performance of both methods. We ran 10 simulations of syntheses from the method of [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF], with k max = 64 for the Cox Voronoi example, and k max = 128 for the Turbulence Poisson example, with 400 iterations per point. For both models we perform also the syntheses from our multiscale gradient descent method (with 400 iterations).

The two families of examples of syntheses share the same corresponding, original pattern of points coming from the original distributions. We also simulate other 9 patterns for the original distributions to compare the averaged statistics and the diversity among the original distribution and the models.

We compare the different distributions with visual evaluation (Figure 5.8), TDA (Figure 5.9), and by estimation of the SCDF (Figure 5.10). For the evaluation with TDA, in addition to the visualization of a 2-dimensional representation of the distance matrix between (the PD of) all original and synthesized point patterns (cf. Section 5.6.3), we also computed the average Wasserstein distance between all pairs of point patterns belonging to different distributions. In more details, for a given distribution (Cox Voronoi or Turbulence Poisson), let M := M orig/GD+W P H be the 10×10 distance matrix between the 10 realizations of the original distribution and the 10 realizations of our model. We compute d orig/GD+W P H = 1

Regularized model and proof of Theorem 5.3.1

Note that eq. (5.2) admits at least one global minimum, which is the observation φ itself. By building the model defined in Section 5.3, one hopes to generate new configurations φ such that E φ( φ ) 0, and φ = φ. If the descriptor K is well chosen, then local minima of E φ correspond to configurations that are likely to have been drawn from L Φ (i.e. the model is a good approximation of the true distribution). The landscape of E φ is defined by the descriptor K. One might find itself in the situation where a certain descriptor K gives poor local minima, (i.e. K does not satisfy the sufficiency property (P2), so samples from the model defined with K are not similar to φ), but incorporating additional constrains to K result in a model for which the only local minimum is φ itself. This is the case, for instance, when adding new constraints to K leads to a violation of the concentration property (P1), due to the high variance of the new descriptor, as it will be illustrated in Figure 5.42 .

To address this issue, we propose an alternative model allowing for the use of a descriptor K that may not satisfy (P1). It consists in adding a regularization term to the energy E φ(µ) in (5.2), so as to prevent the gradient descent algorithm converge to the global minimum which is φ. This can be realized, for example, using a Wassertein-type distance between the current and the initial configuration of points, as described in what follows.

The regularization that we introduce in this section aims at preventing the successive iterations of the gradient descent model considered in Section 5.3.1 to move too far away from its initial configuration. It is achieved by adding the distance between the initial configuration and optimized configuration as a penalization term in the objective. It thus forces the model to explore local minima of E φ(•) around the initial configuration. As all the configurations are measures with equal mass, it is natural to consider the Wasserstein distance between them. However, as this distance is computationally expensive to compute, we choose to replace it with the sliced Wesserstein distance [START_REF] Bonneel | Sliced and radon wasserstein barycenters of measures[END_REF]. For counting measures with same mass, the 2-Wasserstein distance is defined as:

W 2 (µ, ν) = min σ∈S N N i=1 ||u i -v σ(i) || 2 1 2
, where the minimum is taken over the set S N of permutations of (1, . . . , N ).

The square SW 2 distance between two planar counting measures µ and ν is then defined as follows:

SW 2 2 (µ, ν) := θ∈S 1 W 2 2 (θ # µ, θ # ν)dθ,
where θ # designates the pushforward by the orthogonal projection operator in the direction θ. The use of SW 2 instead of W 2 is motivated by the fact that applying the projection operators leads us to compute W 2 on R instead of R 2 , which is much faster, as we simply need to sort the Dirac measures to obtain the optimal matching. It has been shown [START_REF] Bonnotte | Unidimensional and Evolution Methods for Optimal Transportation[END_REF] that SW 2 is indeed a distance, and that it induces the same topology as W 2 for compact domains. One can approximate SW 2 by choosing a certain number of fixed directions [START_REF] Bonneel | Sliced and radon wasserstein barycenters of measures[END_REF].

We use this distance with respect to the initial configuration φ0 ∈ M s as a regularization term in the minimization of E φ(•) defined in Eq. (5.2). More precisely, we consider the following optimization problem arg min

f ∈F (R 2 ) 1 2 |K(f # φ0 ) -K( φ)| + η 1 2 SW 2 2 (f # φ0 , φ0 ), (5.15)
where the minimization is done over the set F(R 2 ) of measurable functions from R 2 into itself, and η ≥ 0 is a regularization parameter. This new optimization problem implies that we do not try to sample from the microcanonical set here, as we suppose that it does not contain enough realizations of Φ.

Note, in (5.15) we optimize the transport f # of a given initial configuration φ0 towards φ, by minimizing the value of the function

E φ0 φ (f ) = E φ(f # φ 0 ) + η 1 2 SW 2 2 (f # φ0 , φ0 ) = 1 2 |K(f # φ0 ) -K( φ)| + η 1 2 SW 2 2 (f # φ0 , φ0 ).
In a similar way as eq. (5.7), for any initial measure φ0 and any measurable function f , we define a function from

R 2 to R 2 x -→ ∇ u E φ0 φ (f ) := η∇ u SW 2 (f # φ0 , φ0 ) + ∇ u E φ(f # φ0 ).
The sequence of point-configuration transports induced by gradient descent algorithm is defined recursively by taking an initial configuration φ0 , f 0 := I d (identity) and for n ≥ 0

   f n+1 : x -→ u -ζ∇ u E φ0 φ (f n • . . . • f 0 ) φn+1 = f n+1 # φn = (f n • . . . • f 0 ) # φ0 .
, (5.16) where ζ > 0 is some gradient step size. In analogy to (5.8) and (5.9), for all n ≥ 1, φn transports the initial distribution L 0 (of φ0 ) to some L n . Note L k = L Φn is the distribution of the point process Φn obtained by n iterations of (5.16) started of Φ0 having law L 0 = L Φ0 .

Proof of Theorem 5.3.1

Proof. In order to prove Theorem 5.3.1, we need to formally define Eq. (5.7). Recall that in this section and in what follows, W s is interpreted as endowed with the addition and scalar multiplication modulo W s .

For µ ∈ M s and any u ∈ Supp(µ) 3 , we define the following functions:

h µ u : R 2 -→ M s v -→ µ -δ u + δ u+v , K µ u : R 2 -→ C d ≡ R 2d v -→ K • h µ u (v),
Using (5.22) and (5.23) in (5.21) we then obtain

∇ T Fn(u i ) E φ(T F n ) = ∇ Fn(u i ) K(F n )A t t K(F n ) -K( φ) = A ∇ Fn(u i ) K(F n ) t K(F n ) -K( φ) = A∇ Fn(u i ) E φ(F n ).
The computations for the regularization term follow the same lines, as the gradient of the regularization function is equivariant to T and, the regularization itself is invariant. Finally, we obtain

f n+1 (T u, T # φ0 ) = T u -A∇ u SW 2 (F n (•, φ0 ) # φ0 , φ0 ) -A∇ u E φ(F n (•, φ0 ) # φ0 ) = T (u -∇ u SW 2 (F n (•, φ0 ) # φ0 , φ0 ) -∇ u E φ(F n (•, φ0 ) # φ0 )) = T f n+1 (u, φ0 ),
which concludes the proof.

Relaxing the assumptions on the data

In order to present our model in a simple setting, strong theoretical assumptions have been made on the data. However, in real world applications, the data will most likely not satisfy these assumptions. This section presents ideas on how to adapt our model in such cases.

Non-periodic boundaries Recall that our descriptor, defined in (5.11), applies periodic boundary correction to point patterns in a square window. If the structure of the observed pattern is not periodic, one can modify the descriptor by applying non-periodic integrals in (5.11) over some smaller window. In particular, we suggest a scale-dependent reduction of the integration window, pertinent when the wavelet ψ has a compact (or approximately compact) spatial support. Specifically, we consider a new descriptor K by considering the integrals in (5.11) with i = (ξ, k, ξ , k , τ ) ∈ Υ over smaller windows W s i ⊂ W s , such that boundary effects are negligible. Our current software can also handle such non-periodic boundary conditions.

More general observation windows Additionally, we considered that the observed pattern lies in a square observation window. If this is not the case, one could use a similar idea to the non-periodic case: embed the observation window in a square window and considering integrals in (5.11) over the observation window.

Non stationary process

In [START_REF] Koňasová | Stochastic reconstruction for inhomogeneous point patterns[END_REF], the authors focus on building a model for non stationary point processes inspired by [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF]. Similarly, one might adapt our method to model non stationary processes. This could be done by modifying two aspects of the method. First, the initial distribution Φ 0 (cf. Section 5.3.1) could be chosen as a non stationary Poisson point process, estimating the intensity with a kernel estimator, such as in [START_REF] Koňasová | Stochastic reconstruction for inhomogeneous point patterns[END_REF]. In addition, as pointed out in [START_REF] Koňasová | Stochastic reconstruction for inhomogeneous point patterns[END_REF], the descriptor should be adapted not to be translation invariant. This could be done, for example, by applying local integrals over patches of the observation window in (5.11) (this requires some notion of "local stationarity" of the process). Another method that may be useful in this scenario is the regularization proposed in Section 5.8, where a regularization term is added to the energy. This term consists of the (Sliced Wasserstein) distance [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF] between the initial configuration and the current configuration (the one being optimized). By adding this regularization term to the energy, the points of the configuration are forced not to move too far away from the initial configuration, which could help preserve the non stationarity of the initial distribution in the distribution of the model.

Discussion

In this chapter, we presented a particle gradient descent model to simulate stationary and ergodic point processes, based on a single observation in a square window. This model is able to synthesize processes formed by a large number of points, exhibiting interactions at multiple scales. Our method is built upon recent works on gradient descent methods to approximate the micro-canonical model. To characterize complex geometric point patterns, we use the wavelet phase harmonic descriptors that allow to explicitly control the scales of the structures to model. Compared to the classical approaches developed in [START_REF] Torquato | Random Heterogeneous Materials: Microstructure and Macroscopic Properties[END_REF][START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF], our approach brought a new perspective to the modeling of point processes, through the lens of wavelet analysis and image modeling.

While numerical results on Cox and Turbulent distributions illustrate the ability of the model to capture various geometric structures in the observation, these distributions are not perfectly well approximated by our model, as shown in Figure 5.6. In order to study what geometric information is missing from our wavelet covariance descriptor, we shall take interest in the modelling of texture images, as they can exhibit complex structures, often more easily perceptible.

Chapter 6

Generalized rectifier wavelet covariance models for texture synthesis

In Chapter 5, we used wavelet phase harmonics covariance statistics (cf. Section 3.5) to represent the distribution of point processes, and built a maximum entropy model based on this representation. We shall now study what possible improvements can be made to this set of descriptors, by focusing on texture images, which can be view as stationary random fields, and can exhibit more complex geometric structures than point processes. The model developed in this chapter can, however, be readily applied to random measures and stationary point processes.

Modelling texture images

Textures ares spatially homogeneous images, consisting of similar patterns forming a coherent ensemble. In texture modeling, one of the standard approaches to synthesize textures relies on defining a maximum entropy model [START_REF] Jaynes | Information theory and statistical mechanics[END_REF] using a single observed image [START_REF] Raad | A survey of exemplar-based texture synthesis[END_REF]. Recall from 2.10, maximum entropy models consist of computing a set of prescribed statistics from the observed texture image, and then generating synthetic textures producing the same statistics as the observation. If the statistics correctly describe the structures present in the observation, then any new image with the same statistics should appear similar to the observation. A major challenge of such methods resides in finding a suitable set of statistics, that can generate both high-quality and diverse synthetic samples. This problem is fundamental as it is at the heart of many texture related problems. For example, in patch re-arrangement methods for texture modeling, these statistics are used to compute high-level similarities of image patches [START_REF] Li | Combining markov random fields and convolutional neural networks for image synthesis[END_REF][START_REF] Raad | A survey of exemplar-based texture synthesis[END_REF]. Such statistics are also used in texture interpolation for probing visual perception [START_REF] Vacher | Texture interpolation for probing visual perception[END_REF], style transfer and image inpainting [START_REF] Gatys | Image style transfer using convolutional neural networks[END_REF][START_REF] Laube | Image inpainting for high-resolution textures using cnn texture synthesis[END_REF].

A key question along this line of research is to find what it takes to generate natural textures. This problem was originally posed in [START_REF] Julesz | Visual pattern discrimination[END_REF], in which the author looks for a statistical characterization of textures. In the classical work of [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF], the authors presented a model whose statistics are built on the wavelet transform of an input texture image. These statistics were carefully chosen, by showing that each of them captured a specific aspect of the structure of the image. This model produces satisfying results for a wide range of textures, but fails to reproduce complex geometric structures present in some natural texture images. Figure 6.1 presents a typical example composed of radishes, and synthetic images from three state-of-the-art models developed over the last few decades. To address this problem, the work of [START_REF] Gatys | Texture synthesis using convolutional neural networks[END_REF] proposes to use statistics built on the correlations between the feature maps of a deep convolutional neural network (CNN), pre-trained on the ImageNet classification problem [START_REF] Deng | Imagenet: A large-scale hierarchical image database[END_REF][START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF]. While this model produces visually appealing images, these statistics are hard to interpret. The work of [START_REF] Ustyuzhaninov | What does it take to generate natural textures?[END_REF] made a significant simplification of such statistics, by using the feature maps of a one-layer rectifier CNN with random filters (without learning). A crucial aspect of this simplification relies on using multi-scale filters, which are naturally connected to the wavelet transform. In this chapter, we propose a wavelet-based model, more interpretable than CNN-based models (with learned or random filters), to synthesize textures with complex geometric structures. It allows one to bridge the gap between the classical work of [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF], and state-of-the-art models. This model is built on the recent development of the phase harmonics for image representations and non-Gaussian stationary process modeling, introduced in [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF][START_REF] Zhang | Maximum entropy models from phase harmonic covariances[END_REF], and presented in Section 3.5 in the stochastic-geometric framework. Recall that the phase harmonics are non-linear transformations that adjust the phase of a complex number. In [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF][START_REF] Zhang | Maximum entropy models from phase harmonic covariances[END_REF], the authors illustrate that the phase dependencies between wavelet coefficients across scales contain important information about the geometric structures in textures and turbulent flows, and that they can be captured by applying the phase harmonics to complex wavelet coefficients. Remarkably, [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF] show that the phase harmonics admit a dual representation, closely related to the rectifier non-linearity in CNNs.

In this work, we present the following results:

• We develop a family of texture models based on the wavelet transform and a generalized rectifier non-linearity, that significantly improves the visual quality of the classical wavelet-based model of [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] on a wide range of textures. It relies on introducing spatial shift statistics across scales to capture geometric structures in textures.

• By changing the number of statistics in our models, we show explicitly the trade-off on the quality and diversity of the synthesis. When there are too many statistics, our model tends to memorize image patches. We further investigate such memorization effects on non-stationary images and find that it sometimes relies on what statistics are chosen, rather than on how many.

• Through the modeling of geometric structures in gray-scale textures, our model indicates the possibility of reducing significantly the number of statistics in the works of [START_REF] Gatys | Texture synthesis using convolutional neural networks[END_REF] and [START_REF] Ustyuzhaninov | What does it take to generate natural textures?[END_REF], to achieve a similar visual quality.

The rest of this chapter is organized as follows: Section 6.2 reviews the framework of microcanonical maximum-entropy models, build upon a general family of covariance statistics. We then present our model for both gray-scale and color textures in Section 6.3. Section 6.4 shows synthesis results of our model, compared with state-of-the-art models. Finally, in Section 6.5, we discuss possible improvements of our model.

Notations

For the remainder of this chapter, N denotes a positive integer. A gray-scale image x is an element of R N ×N , i.e. x = x(u), u ∈ X N , with

X N := {0, • • • , N -1} 2 .
Recall, a color image x = {x c } c=1,2,3 is an element of R 3×N ×N , or equivalently, each x c ∈ R N ×N . We shall denote x the observed texture (observation), which is assumed to be a realisation a random vector X. For any complex number z ∈ C, z * is the complex conjugate of z, Real(z) its real part, |z| its modulus, and ϕ(z) its phase. In this work, for the sake of clarity, the model proposed in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] shall be abbreviated PS, the model presented in [START_REF] Ustyuzhaninov | What does it take to generate natural textures?[END_REF] shall be noted RF (as they use random filters), and the VGG-based model of [START_REF] Gatys | Texture synthesis using convolutional neural networks[END_REF] shall be referred to as VGG.

Microcanonical covariance models

We briefly review the standard framework of micro-canonical maximum-entropy models for textures. To reliably estimate the statistics in these models, we assume that a texture is a realization of a stationary and ergodic process (restricted to X N ). We then review a special family of statistics that are used in the state-of-the-art texture models (mentioned in Figure 6.1), based on covariance statistics of an image representation.

Framework

Given a observation texture x, we aim at generating new texture images, similar but different from x. To that end, a classical method is to define a set of statistics C x, computed on the observation, and try to sample from the microcanonical set {x : Cx -C x ≤ }, where • denotes the L 2 norm.

Under the stationary and ergodic assumption of X, one can construct Cx as a statistical estimator of E(CX), from a complex-valued representation Rx(γ, u) ∈ C as a function of (γ, u) in an index set Γ × X N . The set of covariance statistics Cx of a model can then be constructed by computing an averaging over the spatial variable u, i.e.

Cx(γ, γ , τ

) := 1 |X N | u∈X N Rx(γ, u)Rx(γ , u -τ ) * , ( 6.1)
for (γ, γ , τ ) ∈ Υ ⊆ Γ×Γ×X N . The statistics Cx(γ, γ , τ ) can be interpreted as estimating the covariance (resp. correlations) between RX(γ, u) and RX(γ , u -τ ) for zero-mean RX (resp. non-zero mean RX). The ergodicity assumption ensures that when N is large enough, the approximation C x E(CX) over Υ should hold with high probability. Under these conditions, it makes sense to sample the microcanonical set in order to generate new texture samples. To draw the samples, we follow gradient-based sampling algorithms, suitable in high-dimensions [START_REF] Bruna | Multiscale sparse microcanonical models[END_REF], to minimize the objective function Cx -C x 2 , starting from an initial sample from a normal distribution. This sampling method is similar to the one presented in Section 5.3 to sample from point processes (recall that the method proposed in Section 5.3 use a differentiable mapping from point patterns to images, cf. Section 5.5.1). This framework encompasses a wide range of state-of-the-art texture models, such as [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF][START_REF] Gatys | Texture synthesis using convolutional neural networks[END_REF][START_REF] Ustyuzhaninov | What does it take to generate natural textures?[END_REF][START_REF] Zhang | Maximum entropy models from phase harmonic covariances[END_REF]. In particular, the PS model developed in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] takes inspiration from the human early visual system to define a multi-scale representation based on the wavelet transform of the image [START_REF] Heeger | Pyramid-based texture analysis/synthesis[END_REF]. We next review a family of covariance model which generalizes the statistics in the PS model. We write C M the statistics for a specific model M that uses the representation R M .

Wavelet phase harmonic covariance models

We review a family of microcanonical covariance models defined by a representation built upon the wavelet transform and phase harmonics (Sections 3.1 and 3.5). The wavelet phase harmonics covariance representation is defined for texture images in the same way as for point processes. It defines a class of covariance statistics that capture dependencies between wavelet coefficients across scales.

Wavelet transform

For texture modeling, we consider oriented wavelets to model geometric structures in images at multiple scales. They include the Morlet wavelets (cf. Example 3.1.2) and steerable wavelets (such as Example 3.1.3), proposed in [START_REF] Goupillaud | Cycle-octave and related transforms in seismic signal analysis[END_REF][START_REF] Simoncelli | The steerable pyramid: a flexible architecture for multi-scale derivative computation[END_REF][START_REF] Unser | A unifying parametric framework for 2d steerable wavelet transforms[END_REF]. In particular, the Simoncelli steerable wavelets have been used to model a diverse variety of textures in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF].

Oriented wavelets are defined by the dilation and rotation of a complex function ψ : R 2 → C on a plane. Let r θ denote the rotation by angle θ in R 2 . They are derived from ψ with dilations by factors 2 j , for j ∈ {0, 1, • • • , J -1}, and rotations r θ over angles θ = π/L for 0 ≤ < L, where L is the number of angles in [0, π). The wavelet at scale j and angle θ is defined by

ψ j,θ (u) = 2 -2j ψ(2 -j r θ u), u ∈ R 2
Scales equal or larger than J are carried by a low-pass filter φ J . These continuous wavelets are then discretized (with periodic boundary condition) to define a discrete wavelet transform in the spatial domain X N .

The wavelet transform of an image x ∈ R N ×N is a family of functions obtained by the convolution of x with the discrete wavelets. Let Λ

:= {0, • • • , J -1} × π L {0, • • • , L -1} be
the index set of the wavelet coefficients. For notation convenience, we denote the resulting coefficients by

x ψ j,θ (u) = v∈X N x(u -v)ψ j,θ (v), u ∈ X N , (j, θ) ∈ Λ. (6.2)
The low-pass coefficients x φ J are defined similarly.

Wavelet phase harmonics and the PS model

To model natural textures, it has been shown ( [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF][START_REF] Zhang | Maximum entropy models from phase harmonic covariances[END_REF]) that it is crucial to capture statistical dependencies between wavelet coefficients across scales. This can be achieved by using a wavelet phase harmonic representation, which is defined by the composition of a linear wavelet transform of x, and a non-linear phase harmonic transform.

In [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF], the authors introduce the phase harmonics to adjust the phase of a complex number z ∈ C (cf. Section 3.5). Recall that the phase harmonics {[z] k } k∈Z of a complex number z ∈ C are defined by multiplying its phase ϕ(z) of z by integers k, while keeping the modulus constant, i.e. ∀ k ∈ Z, [z] k := |z|e ikϕ(z) ..

The wavelet phase harmonic representation (WPH) is then defined by

R WPH x(γ, u) = [x ψ j,θ (u)] k -µ γ , γ = (j, θ, k) ∈ Γ = Λ × Z, (6.3)
where µ γ is defined as the spatial average of [x ψ j,θ ] k . It is shown in [START_REF] Zhang | Maximum entropy models from phase harmonic covariances[END_REF] that the PS model can be regarded as a low-order wavelet phase harmonics covariance model, which considers only a restricted number of pairs (k, k ). In the next section, we shall use a dual representation of the phase harmonic operator to define a covariance model to capture high-order phase harmonics.

Generalized rectifier wavelet covariance model

In the previous section, we presented a class of models, built from the wavelet phase harmonic representation. A dual representation of the phase harmonic operator [•] k can be defined via a generalized rectified linear unit, that we review in Section 6.3.1. We then discuss in Section 6.3.2 how to define an appropriate index set of Γ for gray-scale textures. Section 6.3.3 extends the model to color textures.

From phase harmonics to the generalized rectifier

The generalized rectified linear unit of a complex number z, with a phase shifted by α ∈ [0, 2π], is defined by ρ α (z) = ρ(Real(e iα z)), (6.4) where ρ is a rectified linear unit, i.e. for any t ∈ R, ρ(t) := max(0, t). In [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF], it is shown that applying a Fourier transform on ρ α (z) along the variable α results in the phase harmonics of z (up to some normalization constant). This suggests an alternative model, defined by coefficients of the form

R ALPHA x(γ, u) = ρ α (x ψ j,θ (u)) -µ γ , γ = (j, θ, α), ( 6.5)
for γ ∈ Γ = Λ×[0, 2π], and, as before, µ γ is defined as the spatial average of ρ α (x ψ j,θ (u)) over u ∈ X N .

Relation with high-order phase harmonics Based on the duality between the phase harmonics k ∈ Z and the phase shift variable α ∈ [0, 2π], we now present the relation between C ALPHA and the high-order phase harmonics in C WPH . The following result was first proved in [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF], Proposition 6.3.1. There exists a complex-valued sequence {c k } k∈Z such that for all (j, θ, α) ∈ Γ, (j , θ , α ) ∈ Γ, and all τ ∈ X N ,

C ALPHA x((j, θ, α), (j , θ , α ), τ ) = (k,k )∈Z 2 c k c * k C WPH x((j, θ, k), (j , θ , k ), τ )e i(kα-k α ) .
Proof. Using the fact that

ρ(Real(ze iα )) = ρ(Real(|z|e i(ϕ(z)+α) )) = |z|ρ(cos(α + ϕ(z)),
and taking the Fourier transform of ρ α (z) in the variable α, we obtain

F(ρ α (z))(k) := 1 2π [0,2π] ρ α (z)e -ikα dα = |z| 1 2π [0,2π] ρ(cos(α + ϕ(z))e -ikα dα = |z|e ikϕ(z) c k = [z] k c k ,
where c k is the Fourier transform of h(.) := ρ(cos(.)) at the frequency k. The function α → ρ α (z) being periodic in α, we have its decomposition in Fourier series

ρ α (z) = k∈Z F(ρ α (z))(k)e ikα = k∈Z c k [z] k e ikα .
We can then write, for any z, z ∈ C, and α, α ∈ [0, 2π],

ρ α (z)ρ α (z ) * = k,k ∈Z 2 c k c * k [z] k [z ] -k e i(kα-k α ) .
Replacing z and z by any two wavelet coefficients x ψ j,θ (u) and x ψ j ,θ (u -τ ), we thus obtain the relation in Proposition 6.3.1.

We remark that the sequence {c k } k∈Z is uniquely determined by the rectifier nonlinearity ρ, and they are non-zero if k is even ( [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF]). This result shows that for a suitable choice of (α, α ), the covariance statistics C ALPHA x can implicitly capture C WPH x with a wide range of k and k .

Relation with second order statistics Using a simple decomposition of wavelet coefficients into their positive, negative, real and imaginary parts, we can further show that the covariance statistics C ALPHA x capture the classical second order statistics of wavelet coefficients, also used in the PS model (with phase harmonic coefficients k = k = 1). Proposition 6.3.2. Let I = {0, π 2 , π, 3π 2 }. There exists a finite complex-valued sequence {w α,α } (α,α )∈I 2 such that for all (j, θ) ∈ Λ, and all τ ∈ X N ,

(α,α )∈I 2 w α,α C ALPHA x((j, θ, α), (j , θ , α ), τ ) = 1 X N u∈X N x ψ j,θ (u) x ψ j ,θ (u-τ ) * . (6.6)
the size of Υ without redundancies, in Section 6.4.1 we provide an upper bound on the non-redundant statistics in our models. This upper bound is used to count the number of statistics in our models. To keep this number from being too large, instead of taking all shifts in a square box, such as in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF], we choose to select only shifts of dyadic moduli, and with the same orientations as the wavelets. , and the spatial shift correlations are only considered for (j, θ) = (j , θ ). There are two notable differences in the statistics included in the ALPHA S and ALPHA I models. The first one is the range of scales being correlated. It has been shown in [START_REF] Zhang | Maximum entropy models from phase harmonic covariances[END_REF] that constraining correlation between a wider range of scales induces a better model for non-Gaussian stationary processes, and a better estimation of cosmological parameters from observed data [START_REF] Allys | New interpretable statistics for large-scale structure analysis and generation[END_REF]. The second difference, which has a significant impact on the number of statistics (it increases the model size by a factor ∼10), is the number of spatial shifts in the correlations. In the ALPHA I model, spatially shifted correlations are computed for all pairs of coefficient (γ, γ ). For both stationary textures and non-stationary images in gray-scale, shape and contours of salient structures and objects are better reproduced with ALPHA I , as illustrated in Figure 6.2 and Figure 6.3.

(θ, θ ) ∈ Θ 2 (α, α ) ∈ A 4 ×A 1 τ ∈T if (j,θ)=(j ,θ ) τ =0 otherwise. ( J|Θ| 2 + J|Θ||T|)|A 4 | ALPHA I (j, j ) ∈ J 2 (θ, θ ) ∈ Θ 2 (α, α ) ∈ A 4 ×A 1 τ ∈ T J 2 |Θ| 2 |A 4 ||T| ALPHA L (j, j ) ∈ J 2 (θ, θ ) ∈ Θ 2 (α, α ) ∈ A 2 4 τ ∈ T J 2 |Θ| 2 |A 4 | 2 |T|
ALPHA I vs. ALPHA L As we observe in Figure 6.2, the ALPHA I model, containing 4 times less coefficients than the ALPHA L , suffers less from memorization effects, while still capturing most of the geometric information in the images. This small loss of information can be partially explained by the frequency transposition property of the phase harmonics operator [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF], for compactly supported wavelets in the frequency domain, as detailed in the following. First, remark that the simple linear relation

α∈A 4 ρ α (z)e ikα ρ 0 (z) = α∈A 4 ρ α (z)ρ 0 (z)e ikα (6.7)
tells us that computing all correlations for α ∈ A 4 and α = 0 gives us at least the information contained in the r.h.s. of eq. (6.7). Furthermore, if α ∈ A 4 = {0, • • • , 3π 4 }, we make the following approximation

α∈A 4 ρ α (z)e ikα [0,2π] ρ α (z)e -ikα dα = 2πc k [z] k .
Recall also from the proof of Proposition 6.3.2, that F(ρ α (z))(k) = [z] k c k , where the Fourier transform is taken along the variable α. Therefore,

ρ 0 (z) = k∈Z c k [z] k . Therefore, α∈A 4 ρ α (z)e ikα ρ 0 (z) [0,2π] ρ α (z)e -ikα dα ρ 0 (z) = 2πc k [z] k k ∈Z c k [z ] k = 2πc k k ∈Z c k [z] k [z ] k .
Then, replacing z and z with wavelet coefficients, we get

α∈A A e -ikα C ALPHA x((j, θ, α), (j , θ , 0), τ ) 2πc k k ∈Z c k C WPH x((j, θ, k), (j , θ , k ), τ ).
Using Plancherel's theorem, we can write that

C WPH x((j, θ, k), (j , θ , k ), τ ) = 1 |X N | ω∈ 2π N X N F([x ψ ξ ] k )(ω)F([x t τ ψ ξ ] -k )(ω),
where the Fourier transform of an image x is defined by F(x)(ω) := u∈X N x(u)e -iωu , and t τ denotes the translation by τ , i.e. t τ f (•) = f (• -τ ). Now, suppose that the wavelets ψ α have disjoint compact frequency support, in balls B ξ (2 -j C ), where ξ = 2 -j r -θ ξ 0 , and ξ 0 is the central frequency of the mother wavelet ψ (cf. Section 3.1). Suppose also that frequency transposition property of the phase harmonics operator (cf. [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF]) is such that [x ψ ξ ] k has (approximately) frequency support in B kξ (k2 -j C ). Then, for all ξ, ξ , and all k ∈ Z, there exists only one k * such that the frequency supports of [x ψ ξ ] k and [x ψ ξ ] -k * are not disjoint, i.e. such that C WPH x((j, θ, k), (j , θ , k ), τ ) = 0. This tells us that

α∈A A e -ikα C ALPHA x((j, θ, α), (j , θ , 0), τ ) c k c k * C WPH x((j, θ, k), (j , θ , k * ), τ ).
Thus, computing all correlations for α ∈ A 4 , and α = 0 gives us (approximately) all the information contained in WPH coefficients for any pair k, k .

This result lies on several approximations, and strong assumptions about the wavelets, which are not fully met in practice. For this reason, setting α = 0 instead of α ∈ A 4 effectively reduced the amount of information captured by the statistics, and therefore increases the diversity of the model. However, as we observe in Figure 6.2, there is not too much information lost, and the resulting model still captures most of the important geometric structures in texture images. In order to avoid this memorization effect, we shall, in the rest of the chapter, consider only the intermediate model.

examples, and discuss their quality through visual inspection. A quantitative evaluation of the quality of the syntheses, based on the synthesis loss of the VGG model, and proposed in [START_REF] Ustyuzhaninov | What does it take to generate natural textures?[END_REF], is also discussed here.

Experimental Setup

For our experiments, we choose gray-scale and color textures with various visual structures. Our natural texture examples were obtained from the following three sources: CNS NYU1 , Textures.com2 , Describable Textures Dataset model3 and the Github page of [START_REF] Berger | Incorporating long-range consistency in cnn-based texture generation[END_REF] 4 . In the gray-scale examples, we also include a turbulent vorticity flow, which is simulated from Navier-Stokes equations in two dimensions ( [START_REF] Schneider | Coherent vortex extraction and simulation of 2d isotropic turbulence[END_REF]). These examples all contain complex geometric structures that are hard to model by the classical PS model.

The texture images presented in this work have a size of N = 256, giving a number of pixels of ∼65k. For all the ALPHA models, we use Morlet wavelets. This choice differs from the PS model, which uses Simoncelli steerable wavelets. As illustrated in Figure 6.5, this choice can have a visible impact on the quality of the textures. We observe that, while on the first example, the coherence of the structures appear similar for the three wavelet families, the second example shows that the wavelets used in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] The number of orientations is set to L = 4, and the maximal scale is J = 5. The maximal scale parameter also plays an important role in the definition of the wavelet transform. It determines the scales of the structures being captured by the transform. If this parameter is too small, large structures in the observation image might not be captured and reproduced in the model syntheses. Conversely, if J is too large, then the large scale statistics may have a high variance, inducing a memorization effect in the syntheses. Figure 6.6 illustrates this point on two examples from Section 6.4.2. By setting J = 4 (i.e. the maximal range of structures captured by the wavelets is of size 2 4 = 16), we observe on the first example that the larger structures (bubbles) are not well reproduced. When J is set to 6, the observation is almost identically reproduced by the synthesis. Similarly on the second examples, several parts of the synthesis appear very similar to ones in the observation. We found that a suitable trade-off consists in setting J = 5 for images of size N = 256. Similarly to [START_REF] Gatys | Texture synthesis using convolutional neural networks[END_REF] and [START_REF] Ustyuzhaninov | What does it take to generate natural textures?[END_REF], we use the L-BFGS algorithm [START_REF]Updating quasi-newton matrices with limited storage[END_REF] for the optimization of the objective. As in the VGG model ( [START_REF] Gatys | Texture synthesis using convolutional neural networks[END_REF]), we further apply a histogram matching procedure as post-processing after optimization.

Number of statistics

Here, we detail the number of statistics in the different ALPHA models. We begin by giving the formula for each model (note that we do not include the low-pass statistics, which numbers are negligible). To (partially) avoid redundancy in the coefficients, for all models, we compute only the correlations for indices in Υ such that j 2 ≥ j 1 . This gives us the following formulas for the number of statistics:

• #(ALPHA S ) = (2J -1)|Θ 4 | 2 |A 4 | + J|Θ 4 ||A 4 ||T| • #(ALPHA I ) = 1 2 J(J + 1)|Θ 4 | 2 |A 4 ||T| • #(ALPHA L ) = 1 2 J(J + 1)|Θ 4 | 2 |A 4 | 2 |T| • #(ALPHA C ) = 9 2 J(J + 1)|Θ 4 | 2 |A 4 ||T|
Note that, for all ALPHA models, we also compute first order statistics µ γ , i.e. the spatial averages of R ALPHA x(γ, u). There are J|Θ 4 | 2 statistics of this sort in every model, which is negligible with respect to the total number of second order statistics.

Note also that there are still some redundancies in these statistics, as for j = j , all correlations for θ, θ , α, α , τ are counted twice. The number of such statistics is 5 , for all model, superior to the number of first order statistics, which shows that our formula is in fact an upper bound for the exact number of statistics.

This leads us to the following descriptor sizes: for gray-scale textures, the number of the statistics used in the ALPHA I model is ∼35k. The small model ALPHA S has ∼3.5k statistics, while the large model ALPHA L has ∼142k statistics. For comparison, the PS, RF and VGG models have ∼7k, ∼525k , ∼177k statistics respectively. For color textures, the number of statistics of the ALPHA C model is ∼320k, while for the PS, RF and VGG models, it is ∼36k, ∼525k, ∼177k respectively.

Non-periodic boundaries in natural images

The convolution operation in the wavelet transform ((6.2)) is performed using the Fast Fourier Transform. Additionally, recall from Section 6.2.1 that spatial shifts are defined with periodic boundary conditions. This implies periodicity of the input image x. However, natural texture images are not periodic, so one needs to adapt the computation of coefficients to take into account possible border effects. To that end, instead of averaging over all u ∈ X N as in eq. (6.1), each correlation coefficient is averaged over a sub-window inside X N (also proposed for the point process synthesis model, cf. Section 5.9), which size depends on the scales of the coefficients being correlated. More precisely, let γ = (j, θ, α) and γ = (j , θ , α ). Note j m := max(j, j ). We define X jm := {u = (u 1 , u 2 ) ∈ X N : 2 jm ≤ u i < N -2 jm , i = 1, 2}. Then, for non periodic images, we compute (6.9) where the spatial shifts are defined periodically. Note that the spatial averages µ γ and µ γ are also performed on X jm .

C ALPHA x(γ, γ , τ ) = 1 |X jm | u∈X N 1 X jm (u)1 X jm (u -τ )R ALPHA x(γ, u)R ALPHA x(γ , u -τ ),

Results

In Figure 6.7, we present examples of syntheses from the ALPHA I (or ALPHA C ), PS, RF and VGG models, for both gray-scale and color textures, as well as for non-stationary images. We observe that our model ALPHA I produces texture syntheses of similar visual quality to the RF and VGG models. It also significantly outperforms the PS model in terms of the visual quality, without introducing visible memorization effects. As the model PS uses the statistics closer to ALPHA S compared to ALPHA I , the performance of PS is somehow expected. More synthesis examples can be found in Figures 6.9 and 6.10.

Note that for the tiles example (the fifth row) in Figure 6.7, the VGG model produces less convincing textures, because the long-range correlations present in the image (aligned tiles) are not reproduced. To remedy this issue, it has been proposed in [START_REF] Berger | Incorporating long-range consistency in cnn-based texture generation[END_REF] to add spatial shifts to the correlations of the network feature maps. These shift statistics are similar to the parameter τ in our model. We also observe that, in the case of the sixth row example (flowers), all models fail to reproduce complex structures at object-level. Possible improvements of such models is further discussed in Section 6.5.

For non-stationary images, we find that certain image patches can be more or less memorized by the RF, VGG and ALPHA I models, as illustrated in the seventh row example. Understanding such memorization effect of non-stationary images is a subtle topic, as we find that in some binary images (x(u) ∈ {0, 1}), only the PS model can reproduce the observation, even though it has a much smaller number of statistics (the last row example). We find that this is related to the spatial correlation statistics in PS (non-zero τ ). By removing this constraint, x is no longer always reproduced. This simple example suggests that sometimes it is very important to choose the right statistics to capture specific geometric structures. The non-stationary nature also appear in some logo near the boundary of some textures (e.g. bottom left in the observation of the first and fourth rows). Although this logo is reproduced by RF, VGG, ALPHA I and ALPHA C , it is a very local phenomenon, as we do not find visible copies of the textures when there is no logo, and it is likely related to the way one addresses the boundary effect (cf. Section 6.4.1).

VGG score

In [START_REF] Ustyuzhaninov | What does it take to generate natural textures?[END_REF], the authors proposed to use the synthesis loss of the VGG model to evaluate the quality of syntheses from any model. The goal is to define a quantitative, and more objective evaluation method than mere visual inspection. Since the VGG model produces syntheses almost indistinguishable from real textures, it is natural to consider its loss to asses the quality of a synthesis. We computed this loss for the examples presented in Figure 6.7, and report it in Table 6.2. Note however that this loss is not exactly the same as the one used in [START_REF] Ustyuzhaninov | What does it take to generate natural textures?[END_REF], as the layers selected to compute the loss are different. In this work, we chose to use the layers suggested in [START_REF] Gatys | Texture synthesis using convolutional neural networks[END_REF], i.e. 'conv1_1', 'pool1', 'pool2', 'pool3', and 'pool4' of the VGG-19 network ( [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF]), and compute the relative VGG loss 6 .

We notice that this score is not always consistent with visual inspection, as there are texture examples and models for which the syntheses do not look much like the observation image, yet produce a small VGG loss (see e.g. the first and last rows of Figure 6.7, the RF model syntheses have the smallest loss). It should also be noted that the VGG loss reported on the VGG syntheses is not the synthesis loss after optimization, as a histogram matching (HM) procedure is performed as post-processing after optimization. We observed that the VGG loss of the syntheses from the VGG model after HM was considerably higher than the one for syntheses before it, while being visually very similar as illustrated in Figure 6.8. These observations suggest that the VGG score suffers from instabilities after reaching a certain level (that is, if the VGG loss is small enough, small perturbations of the values of the image pixels might have a strong impact on the loss).

Discussion

In this work, we presented a new generative model for texture images built on a waveletbased representation of geometric structures. Inspired from the classical work of [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF], this model incorporates a wider range of statistics, by computing covariances between rectified wavelet coefficients, at different scales, phases and positions. We showed that this model is able to capture and reproduce complex geometric structures, present in natural textures or physical processes, producing syntheses of similar quality to state-of-the-art models that use CNN-based representations.

Being defined with a wavelet family instead of multi-scale random filters, the proposed model uses less statistics than the RF model for color textures. For the gray-scale textures, our model has about 15 times less statistics, as it focuses on capturing the geometric structures present in images. However, further reduction of the number of statistics remains an open problem. For instance, the size of the model could be reduced with a more sophisticated selection of the spatial shift parameter τ , depending on the scales of the coefficients being correlated. Additionally, our color model is defined by computing covariances between coefficients of all different color channels, and is therefore nine times larger than the gray-scale model. In model of [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] extended to color images performs a Principal Component Analysis on the color channels ( [START_REF] Vacher | The portilla-simoncelli texture model: towards understanding the early visual cortex[END_REF]). A similar idea could also be applied to our model for color images. Another line of research to reduce the size of the model could be to incorporate a second layer of wavelet transform ( [START_REF] Leonarduzzi | Maximum-entropy scattering models for financial time series[END_REF][START_REF] Sifre | Rotation, scaling and deformation invariant scattering for texture discrimination[END_REF]). Such statistics could also capture additional information to improve the quality of the synthesized textures.

Finding a minimal set of statistics to define a texture model is important because a large number of statistics can result in a high variance of the estimators, and the associated model can suffer from memorization effects. That is, parts of the observation image can be reproduced in the syntheses. This is a problem because the aim of the model is to approximate the underlying distribution of the observation, and therefore produce diverse textures. In this regard, the mere visual evaluation of the synthetic textures can fail to take this aspect of the model into account. A more general evaluation of texture models is an interesting line of research, left for future works.

Chapter 7 Conclusion

This dissertation studied different wavelet-based representations of random geometric processes, and in particular point processes. Similarly to classical summary characteristics, these representations capture geometric structures by applying non-linear transformations to the data to extract relevant information about the inter-point distances. However, the wavelet transform separates the structures present in the data at different scales, which allows one to decompose the information into different parts, making such representations suitable for statistical tasks such as regression or probabilistic modelling.

Summary of findings

We first studied the scattering moments representation of point processes (cf. Section 3.2), which decompose the geometric structures into different frequency bands, and those bands into different frequency sub-bands. We observed that this relatively concise representation ( 1k statistics) is able to efficiently represent the geometry of point processes, and is suitable to perform the task we considered, which was to estimate different geometric marks associated with point patterns.

The scattering moments, however, do not capture the dependencies between the different frequency parts of the patterns. Additionally, the complex modulus operator discards the phase of the wavelet coefficients, and induces symmetries which may lead to the loss of important information about the data. For these reasons, we took interest in the wavelet phase harmonics covariance moments, which compute correlations between phase adjusted wavelet coefficients at different frequency bands. We found that this representation is able to capture the geometry of various point process distributions, and allows to efficiently model these distributions, compared to the summary characteristics usually considered in the literature.

To further improve the representation of complex geometric objects, we turned to the modelling of texture images, where state-of-the-art methods are based on the representation of images with convolutional neural networks. By noticing the important amount of information carried by spatially shifted wavelet coefficients, we built a wavelet covariance model that performs similarly to these CNN-based state-of-the-art models, illustrating that such wavelet-based model can capture most of the geometric information present in random processes such as natural texture images.

Future perspectives

While this work presents empirical evidence that wavelet-based moments can adequately describe point patterns and point processes, little is known about their mathematical properties, such as their discriminating power, their limit behaviour, their stability to deformations, their dependency on different geometric marks, or their ability to fully characterize a certain class of point processes. Recent advances have been made in this area (see e.g. [START_REF] Perlmutter | Scattering statistics of generalized spatial poisson point processes[END_REF]), and further efforts in this direction is an interesting line of research.

The optimization method presented in Chapter 5 differs from classic gradient-descent models, such as in Chapter 6, as it moves mass around in space, rather than changing the values of each pixel. Preliminary results indicate the possibility of synthesizing physical fields (such as a turbulent flows, or an observation of cosmic gas) with this method, by approximating this field by a sum of Gaussian functions. More precisely, we generate x(u) = i exp(-(u -u i ) 2 /2σ 2 i )/(σ i √ 2π), such that Cx ∼ C x, by optimizing the parameters x i and σ i , that is, the positions and variances of the Gaussian functions. It could be interesting to study if this method can improve synthesis models, e.g. by avoiding shallow local minima of the objective function.

Another promising line of research is the use of Topological Data Analysis (cf. Section 5.6.3) to evaluate the error between a distribution and a model. This method seems to be able to distinguish between various types of point process distributions, but also between different samples from a same distribution, enabling to evaluate the diversity of a model. TDA has been proposed for analysing images in learning tasks (see e.g. [START_REF] Kim | Pllay: Efficient topological layer based on persistent landscapes[END_REF], so it could be promising to study this tool as an evaluation of texture synthesis models.

The scattering moments and phase harmonics covariances offer representations of somewhat different nature: the former captures rich information of processes at each scale, with an architecture similar to these of CNNs by using a cascade of linear and non-linear operators, while the other efficiently captures correlations between the different scales. It is easy to see that spatially shifted correlations, as in eq. (6.1), have a direct link with a second layer of convolution, by noting that Rx(γ, •) ψ(u) = Rx(γ, u -τ )ψ(τ )dτ . A unifying way to encode the information contained in both types of moments (scattering and phase harmonics covariance) could be an interesting way to further improve the representation of geometric objects.

MOTS CLÉS

processus ponctuels, traitement du signal, apprentissage statistique, synthèse RÉSUMÉ Cette dissertation présente une classe de représentations de processus ponctuels. Inspirés par le succès des méthodes d'ondelettes en traitement du signal, ces descripteurs sont basés sur la convolution d'un processus ponctuel avec un famille d'ondelettes. À partir de ces convolutions sont construits des ensembles de descripteurs statistiques de processus ponctuels stationnaires, en appliquant des opérateurs non linéaires, suivis d'un moyennage spatial. Tout comme les caractéristiques classiques pour les processus ponctuels, ces statistiques sont conçues pour extraire les informations contenues dans le processus en un nombre relativement faible de valeurs numériques, en décrivant sa géométrie. Leur but est de décrire la façon dont les atomes du processus ont tendance à se repousser, ou bien se regrouper, et de la sorte former des formes géométriques complexes. De par leur construction, ces descripteurs bénéficient de plusieurs propriétés qui les rendent adaptés à des tâches d'apprentissage et d'analyse statistique. Afin d'illustrer la qualité de ces représentations en tant que descripteurs statistiques, nous étudions plusieurs problèmes impliquant l'analyse statistique de processus ponctuels. Dans un première expérience, nous cherchons à estimer une fonction inconnue qui prend en argument une configuration de points, et renvoie une version marquée de cette configuration, c'est-à-dire pour laquelle une valeur numérique est associée à chaque atome de la configuration. Nous utilisons une représentation en ondelettes de ces configurations pour estimer la relation entre leur versions non marquées et marquées. Dans un second temps, nous étudions la capacité de ces descripteurs à modéliser certaines distributions de processus ponctuels, en définissant un modèle de maximum d'entropie défini par des statistiques d'ondelettes, calculées sur une unique observation. Pour ces deux problèmes, nous observons que les représentations que l'on propose amènent de meilleures performances que les statistiques classiques couramment utilisées dans la littérature sur les processus ponctuels. Enfin, pour étudier à quel point de telles représentations peuvent capturer les structures géométriques présentes dans les textures, nous définissons un model de maximum d'entropie qui s'appuie sur des statistiques d'ondelettes similaires, produisant des synthèses de qualité comparable à celles de l'état de l'art, dont les modèles sont basés sur des representations à partir de réseaux convolutionels profonds.

ABSTRACT

This dissertation presents a class of representations of spatial point processes. Inspired from the success of wavelet methods in signal processing, these descriptors rely on the convolution of a point process with a family of wavelet filters. From these convolutions are built sets of statistical descriptors of stationary point process, by applying non-linear operators, followed by a spatial averaging. Much like classical summary characteristics for point processes, these statistics are designed to extract information about the process with a relatively small number of numerical values, by describing its geometry. Their goal is to describe whether the atoms of the process tend to repel each other, or cluster together, and by doing so, form possibly complex geometric shapes. By construction, these descriptors enjoy several properties that make them suitable for statistical analysis and learning tasks. To illustrate the quality of these representations as statistical descriptors, we study several problems involving statistical analysis of point processes. In a first experiment, we seek to estimate an unknown function that takes as input a point pattern, and returns a marked version of this pattern, where a numerical value is associated to each atom of the pattern. We use a wavelet-based representation of point patterns to estimate the relation between their non-marked and marked version. We then study, in a second experiment, the ability of such representations to model the distribution of a point process, by defining a maximum entropy model defined by a set of wavelet-based statistics, computed on a single observation. For these two problems, we observe that our representations lead to better performance than summary statistics commonly used in the literature on point processes.

Figure 1 . 1 :

 11 Figure 1.1: Left: sample from a Cox-cirles point process (see Example 2.8.4). Right: plot of the (theoretical) Bartlett spectrum density, and its estimation using the Discrete Fourier Transform.

Figure 1 . 2 :

 12 Figure 1.2: Illustration of the asymptotic behaviour of first order scattering moments. Left: Poisson processes (PPP) with different intensities, a Gaussian determinantal (DPP), and a Matérn cluster (MCP) process exhibiting the same type of macroscopic variance. Right: a Poisson process, compared to hyperuniform processes (more regular), see Section 3.3.

Figure 1 . 3 :

 13 Figure 1.3: Summary of the learning of geometric marks. Left: overview of the method using the scattering moments representation. Right: an example of estimated marks, compared to the exact ones, which depend on the Voronoi cell areas of each atom (see Example 2.6.4).

Figure 1 . 5 :

 15 Figure 1.5: Top: examples of gray-scale and color texture images. Bottom: examples of syntheses from our texture model.

  Chapter 5]). Classical examples of determinantal point processes are the Gaussian determinantal point process, with a kernel of the form k(u, v) = αe -u-v 2 β

Figure 2 . 2 :

 22 Figure 2.2: Theoretical (black full line) and estimated (blue dotted line) Bartlett spectrum for the Cox-circles point process with parameters R = 10/256 and 2πRλ = 25. The density of the Bartlett spectrum is normalized by the total intensity λ, thus, by (2.27), the curve is invariant to the density of circle centers (a particular numerical scenario considered here is presented in Figure 1.1 left). Observe in particular, log 10 b BC (0)/λ = log 10 (1 + 2πRλ ) = log 10 (26) ≈ 1.415.

Figure 3 .

 3 [START_REF] Abry | Wavelet analysis of long-range-dependent traffic[END_REF] shows the real and imaginary parts of a Morlet wavelet, for a given ξ 0 , and Figure3.2 displays a Morlet wavelet family in the spatial and frequency domains.

Figure 3 . 1 :

 31 Figure 3.1: A Morlet wavelet. Left: Real part, right: imaginary part.

Figure 3 . 2 :Figure 3 . 3 :Definition 3 . 1 . 4 (

 3233314 Figure 3.2: A Morlet wavelet family. Left: spatial domain (real part), right: frequency domain, centered around 0 (real part).

Figure 3 . 4 :

 34 Figure 3.4: A bump steerable wavelet family. Left: spatial domain (real part), right: frequency domain, centered around 0 (real part).

Figure 3 .

 3 6 indicate. Example 3.3.8 (Poisson, Gaussian determinantal and Matérn cluster processes). Figure 3.5 illustrates asymptotics of scattering moments for three different distributions of Φ verifying σ 2 (1, Φ) > 0 with no marking (m = 1). It shows the log 2 of

Figure 3 . 5 :

 35 Figure 3.5: Asymptotics of scattering moments for Poisson, Gaussian determinantal (DPP), and Matérn cluster (MCP) processes, estimated on 500 realizations.

Figure 3 . 6 :

 36 Figure 3.6: Asymptotics of scattering moments for various processes, estimated on 500 realizations. Left, a Poisson process is compared with a matching process, middle: with a marked Poisson process, with Voronoi area marks, right: with a Ginibre process and a crossing process.

2

  Indeed, the partition {∠i(r)}i=1,••• ,6 can be defined by 6 points v1, • • • , v6, such that for all i, the angle viu1vi+1 equals π 3 . Let v be the atom of Φ lying inside the subset of Bu 1 (r) inside {∠i(r)}. Since viu1vi+1 = π 3 , and v -u1 ≤ r the bisector of [u1, v] intersects [u1, vi] and [u1, vi+1].

Figure 4 . 1 :

 41 Figure 4.1: Exact (blue) and reconstructed (orange) marks being the surface areas of the Voronoi cells.

Figure 4 . 2 :

 42 Figure 4.2: An overview of our method for the reconstruction of geometric marks via scattering moments.

  forming a training set of data. The problem consists in learning the function m from the training set, so as to be able to calculate approximations of the unobserved marks m i = m(φ, u i ) for a new realization φ = i δ u i of Φ.

  Sections 4.1.2 and 4.1.3). Our goal is to find a common relation between X k and Y k for all samples k, and the simplest possible one is a linear relation, represented by some matrix B and vector β 0 such thatBX k + β 0 ≈ Y k for all k = 1, . . . , n. (4.1)If the linear function does not allow one to capture the dependence, we can use kernel regression, or more advanced machine learning methods. For simplicity, in this study, we focus on the linear ridge regression, detailed in what follows.To find the relation (4.1), we will use linear ridge model; cf [97, Section 7.5]. For λ = (j, θ), denote by β(λ) the line of the matrix B (4.1), and similarly the component β 0 (λ) for the vector β 0 . Let Y k (λ) := Ŝ φk (λ) be the (j, θ)-component of Ŝ φk . The ridge model consists in minimizing the regularized sum of the squared residuals n k=1

Figure 4 .

 4 3 presents an example of the reconstruction of marks of one given image, with our method, and Figure4.4 shows Q-Q plots for various reconstruction variants. Table

5

 5

Figure 4 . 3 :

 43 Figure 4.3: Reconstructed image example (mapped to 2D; the peaks of the lines correspond to the values of marks of points numbered in lexicographic order; orange curve -reconstruction, blue -exact values).

Figure 4 . 4 :

 44 Figure 4.4: Shot-noise reconstruction. Q-Q plots (reconstructed mark in function of its true value) for 100 test images (a) estimated scattering moments, (b) exact scattering moments, (c) benchmark.

Figure 4 . 5 :Figure 4 . 6 :

 4546 Figure 4.5: Reconstructed image example for the nearest neighbour model.

Figure 4 . 7 :

 47 Figure 4.7: Reconstructed image example for the Voronoi cell surface model.

Figure 4 . 8 :

 48 Figure 4.8: Voronoi cell surface area reconstruction; example and Q-Q plots using: (a) estimated scattering moments, (b) exact scattering moments, (c) benchmark.

Figure 4 . 9 :

 49 Figure 4.9: Reconstructed image example for the Voronoi cell moment of inertia model.

Figure 4 . 10 :

 410 Figure 4.10: Voronoi cell moment of inertia reconstruction; example and Q-Q points using (a) estimated scattering moments, (b) exact scattering moments, (c) benchmark.

Figure 4 . 11 :

 411 Figure 4.11: Reconstructed image example for the Voronoi shot noise model.

Figure 4 . 12 :

 412 Figure 4.12: Voronoi shot-noise reconstruction; example and Q-Q points using (a) estimated scattering moments, (b) exact scattering moments, (c) benchmark.

Figure 4 . 13 :

 413 Figure 4.13: Reconstruction of a sample of the Voronoi shot-noise marks from the estimated scattering moments after 1, 5 and 50 iterations. 5 iterations give the smallest RMSE of 4.17e-2, while 50 iterations used for the reconstruction from the exact scattering moments improve the reconstruction of some marks but give worse RMSE of 5.26e-2.

Figure 4

 4 Figure 4.14: Relative errors m (j, θ) for (j, θ) in order (0, 0), (1, π/8), . . . , (2, 0), . . . , (7, 7π/8); calculated on training and test set.

Figure 5 . 1 :

 51 Figure 5.1: Samples of point processes of various geometries. The number of points ranges from 1000-13000.

( P1 )

 P1 Concentration property: The value of K( Φ) should concentrate around its mean, i.e. K( Φ) E[K( Φ)] with high probability. A natural assumption is that the variance of K( Φ) is small. (P2) Sufficiency property: The moments E[K( Φ)], should characterize the unknown distribution as completely as possible. It requires that K has a strong (distributional) discriminate power.

Figure 5 . 2 :

 52 Figure 5.2: Overview of the particle gradient descent method: starting from φ0 ∈ M s , we iteratively modify the configuration until reaching Ω .

  3.

Figure 5 . 3 :

 53 Figure 5.3: The differentiable discretization of a courting measure µ into a point-image µ N σ on a N × N regular grid.

Figure 5 . 4 :

 54 Figure 5.4: Three Cox models; original sample, reconstruction, and synthesis. For the power spectrum plots, the full lines correspond to the original distributions, and the dashed lines correspond to the models. The y-axis is presented in log scale.

100 iFigure 5 . 7 :

 10057 Figure 5.7: Syntheses from the model of[START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF] and our model, from observations containing 500, 2000 and 3784 points. The last line shows estimations of the spherical contact distribution function H s (r), for the corresponding models and methods, averaged over 10 realizations. For the three plots, the full line (blue) corresponds to the observation distribution, the (green) dashed line to the method (RS+NND), and the (orange) dashdotted line to (GD+WPH).

Figure 5 . 8 : 61 Figure 5 . 9 :

 586159 Figure 5.8: Visualization of syntheses for the Cox Voronoi and Turbulence Poisson examples. Left: observations, middle: RS+NND, right: GD+WPH.Cox VoronoiTurbulence Poisson

Figure 5 . 10 :

 510 Figure 5.10: Estimation of the spherical contact distribution functions (SCDF), for the Cox Voronoi example (left), and the Turbulence Poisson example (right). For the two plots, the full (blue) line corresponds to the observation distribution, the dashed-dotted (green) line to the method (RS+NND), and the dash-dotted (orange) line to (WPH+GD).

Figure 6 . 1 :

 61 Figure 6.1: Example of syntheses from three texture models in chronological order. From left to right: the observed texture in gray-scale, synthesis from PS ([92]), from VGG ([54]), and from RF ([112]).

Figure 6 . 5 :

 65 Figure 6.5: Comparison between different wavelets families used in the ALPHA I model. Central zooms of syntheses using the same covariance model, with three different wavelet families. From left to right: observation, Simoncelli steerable wavelets, bump steerable wavelets (cf. Example 3.1.3), and Morlet wavelets (Example 3.1.2).

Observation J = 4 J = 5 J = 6 Figure 6 . 6 :

 45666 Figure 6.6: Syntheses from the ALPHA C model defined with three different maximal scale parameter J ∈ {4, 5, 6}.

  Data / Model ALPHA I VGG PS RF Radishes 5.02e-05 1.87e-05 2.37e-04 1.13e-05 Cherries 4.86e-05 1.47e-06 6.68e-04 9.65e-06 Gravel 5.97e-05 3.08e-06 7.25e-04 1.29e-05 Turbulence 5.59e-05 5.97e-05 2.42e-04 3.95e-05 Table 6.2: Relative VGG loss of each model on examples in Figure 6.7. Each row corresponds to the same row in Figure 6.7.

  Finally, to study to what extent such representations can capture geometric structures of texture images, we define a maximum entropy model relying on similar wavelet statistics, yielding syntheses of similar visual quality to state-of-the-art models based on deep convolutional neural networks representations. KEYWORDS point processes, signal processing, statistical learning, synthesis

  Figure 2.1: Examples of different point process distributions. Left: Matern II point process (repulsive), middle: homogeneous Poisson point process, right: modified Thomas point process (clustering).is a Gaussian function, then the resulting process is called a modified Thomas process. If the density is uniform in a ball, it is called a Matérn cluster process.

Also called doubly stochastic point processes, Cox point processes constitute a rich class of processes, often used to model patterns exhibiting more clustering than the Poisson point process. Let Π be a random, locally finite, non null measure on (E, B). A Cox point process Cox Π on (E, B) generated by (or directed by) Π is the mixture of Poisson point processes with respect to Π, that is, conditioned on Π = Λ, Cox Π is the Poisson point process with mean (deterministic) measure Λ. A rich class of Cox processes is constituted by Neyman-Scott processes (see e.g.

[START_REF] Neyman | Statistical approach to problems of cosmology[END_REF]

, or [37, Section 5.3]), where 'parent points' are drawn from a homogeneous Poisson distribution, and from each parent point p i descend N i children, with N i being a Poisson r.v., and the location of the children a determined by a given density f . For instance, if this density f

  To overcome this issue, Palm distributions were introduced, based on the Radon-Nikodym theorem. (Campbell measure and Palm distributions). Let Φ be a point process on E, with Radon mean measure M . The Campbell measure of Φ is defined on

	Definition 2.4.1

Corollary 2.5.10. Under

  the assumptions of Theorem 2.5.9, 0 ∈ Φ P 0 -almost surely.

	Definition 2.5.11 (Reduced second-order density). Let Φ be a point process defined on
	the stationary framework (Ω, F, {θ t }, P), compatible with the flow, with finite non-null
	intensity λ. Suppose that M Φ (2) is σ-finite, and has a density w.r.t. the Lebesgue measure.
	As Φ is stationary, we can write M Φ (2) (d(u, v)) = λκ(v -u)dudv, and for any B ∈ B,

.

[START_REF] Błaszczyszyn | How user throughput depends on the traffic demand in large cellular networks[END_REF]

The results extends for any function f for which either one of the two sides of the above equation is finite when f is replaced by |f |.

Definition 2.10.1 (

  Kullback-Leibler divergence and entropy of a point process). Let L 1 , L 2 be two probability distributions on M, such that L 1

						L 2 . The Kullback-Leibler
	divergence (or KL divergence) of L 1 w.r.t. L 2 is defined by
		KL(L 1 ||L 2 ) :=	M	log(	dL 1 dL 2	)dL 1 .
	It is well defined, as shown in Lemma 2.10.2. If L 0 is the Poisson distribution on (E, B),
	and L	L 0 , then
		H(L) := -KL(L||L 0 )
	is called the entropy of L.
	We shall need the following result:
	Lemma 2.10.2. Let L 1 , L 2 be two probability distributions on M, such that L 1	L 2 .
	Then,
		KL(L 1 L 2 ) ≥ 0,	(2.31)

.36) Proposition 2.10.5. The

  solution to eq. (2.36) is given by L * such that

  , and I is the identity matrix of dimension n + 1. Using (4.2), one can calculate approximations Ŝ φ(λ) of the scattering moments of a new marked configuration φ

Table 4 .

 4

	method	scattering scattering exact 5 000 10 000 scatt.	bench-mark
	RMSE	1.99	1.98	1.64	1.98
	NRMSE1 9.21e-2	9.15e-2	7.56e-2 9.14e-2
	NRMSE2 3.32e-1	3.29e-1	2.75e-1 3.29e-1
	method	estimated scatt. exact scatt. benchmark
	RMSE	3.14e-2	2.06e-2	-
	NRMSE1 1.21e-1	7.94e-2	-
	NRMSE2 3.94e-1	2.58e-1	-

1: Shot-noise reconstruction errors using different methods. because the nearest neighbour distance is equal to the first element of the distance matrix. Thus, the scattering moment approach cannot do better. The number of iterations for the reconstruction from the exact scattering moments and estimated ones is taken, respectively, 250 and 8.

Table 4 . 2

 42

: Nearest neighbour distance; reconstruction error Voronoi cell surface area For this example, the training set X consists of 10 000 realizations of Poisson point process with intensity 30. The reconstruction results are presented on Figures 4.

7

and 4.8

Table 4 .

 4

	method	estimated scatt. exact scatt. benchmark
	RMSE	9.71e-3	7.33e-3	9.60e-3
	NRMSE1 6.60e-2	4.98e-2	6.52e-2
	NRMSE2 2.87e-1	2.16e-1	2.85e-1

3: Voronoi cell surface reconstruction errors.

Table 4 .

 4 4: Voronoi cell moment of inertia reconstruction errors.For this example, similarly to the previous example, the training set X consists of 10k realizations of a Poisson point process with intensity 30. The results are presented in Figures 4.11 and 4.12 and in Table4.5. The neighborhood for the benchmark is K = 15 .

	method	estimated scatt. exact scatt. benchmark
	RMSE	2.42e-4	1.37e-4	2.65e-4
	NRMSE1 5.98-2	3.39e-2	6.53e-2
	NRMSE2 7.12e-1	4.04e-1	7.74e-1
	Voronoi shot-noise

Table 4 . 5

 45

	method	estimated scatt. exact scatt. benchmark
	RMSE	4.18e-2	3.02e-2	4.29e-2
	NRMSE1 9.24-2	6.69e-2	9.50e-2
	NRMSE2 3.76e-1	2.72e-1	3.86e-1

: Voronoi shot-noise reconstruction errors.

Table 4 . 6 :

 46 .6. Observe that the constant sum of the Voronoi cell areas makes them globally dependent, unlike the Voronoi moments of inertia. Choice of the number of neighbours K for the benchmark approach and the corresponding RMSE.

	model	K	10	15	20	35
	Shot-noise		2.00	1.98	1.98	-
	Voronoi area		1.04e-2 1.00e-2 9.80e-3 9.60e-3
	Voronoi inertia		2.76e-4 2.65e-4 2.80e-4 -
	Voronoi shot-noise 4.32e-2 4.30e-2 4.29e-2 -

Table 5 .

 5 1: Speed comparison between random search and gradient descent, in number of iterations (computation time in parenthesis) for the synthesis of Poisson Voronoi patterns.

Table 6 .

 6 1: List of indices in Υ for different ALPHA models. ALPHA S vs. ALPHA I The small model ALPHA S is inspired from the PS model, as it only takes into account of the interactions between nearby scales (i.e. |j -j| ≤ 1)

In principle, one can consider Φ to be a homogeneous Poisson point process on E. However, one cannot obtain a stationary Φ f = Φ with this construction. Indeed, in this case f has to be translation invariant. Therefore, by ergodicity, f ≡ 1. For this reason, such a construction is often considered in a finite window.

Rigorously, a marked point process is a sum of atoms Φ = i δ (u i ,m i ) , with (mi)i being a sequence of random variables. When the mi ∈ C, we adopt the convention Φ = u∈Φ m(u, Φ)δu, which is formally not a marked point process.

R 2 f (u)Φ(du), for any measurable, integrable function f : R 2 → C. In this section, we shall study second order moments based on shot noises for stationary processes, define the Bartlett spectrum, and show how we can estimate it using the Discrete Fourier Transform on a finite window.

More precisely, the wavelets we shall consider, such as Morlet wavelets, verify ψ = ψ * . The results remain valid for these wavelets, and we make the symmetric assumption for mere notation simplicity.

Computations are produced using a software available at https://github.com/kymatio/ phaseharmonics

Note that the violation of (P1) does not imply memorizing φ.

For any counting measure µ, Supp(µ) is defined as the union of all singletons where µ has an atom, i.e. Supp(µ) = x : µ(x) > 0.

http://www.cns.nyu.edu/~lcv/texture/

https://textures.com/

https://www.robots.ox.ac.uk/~vgg/data/dtd/index.html

https://github.com/guillaumebrg/texture_generation

The number of such moments is of the order of J|Θ| 3 for the small model, J|Θ|[START_REF] Klatt | Hyperuniform and rigid stable matchings[END_REF] for the intermediate and large models, and 9J|Θ| 4 for the color model.

Using the code from https://github.com/ivust/random-texture-synthesis/blob/master/vgg_ loss.py (function style_loss_relative).

Hard-core

the observation sample φ. Indeed, this line shows quite faithful reconstructions of the original samples subjected to a periodic translation, up to some precision error due to a finite image resolution N . This is however not a good model because it essentially only contains the observation φ. It suggests that we need to improve the concentration property (P1) of the descriptors in order to enlarge the ensemble Ω . (Note that, in the work of [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF], the authors use the term 'reconstruction' to refer to random sampling method, which we call in this work 'synthesis'.)

In order to improve (P1), we shall reduce the parameter J in the wavelet transform. The third line of Figure 5.4 shows realizations sampled using J = log(N ) -3 for different original distributions. Our analysis in Section 5.5.1 suggests that this range of J can model structures whose spatial size is at most 1/8 of the window W s . Observe that most polygons and circles are well reproduced in the synthesis of Voronoi and Small circles. The Big Circles are harder to model since the size of each circle is slightly larger than 1/8 of W s .

In the last line of produce the 10 syntheses. We see in the first two Cox examples a clear separation between the original process and the model, implying a lack of sufficiency (P2) in our model. This is probably because, in order to satisfy (P1) and produce diverse samples, we have chosen to reduce J, and therefore lose some information about large scale structures of the process. On the other hand, we observe that this separation is smaller for the Turbulence hard-core case, where there is a better balance between (P1) and (P2). The error in the Cox models is probably due to the difficulty to reproduce highly constrained structures (perfect circles or convex polygons). These observations agree with our visual evaluation of the syntheses: the Voronoi and (more particularly so) the Big circles models are easily discriminated from their original distributions (their highly constrained structures are not perfectly reproduced in the syntheses), but this discrimination is harder for the Turbulent hardcore case. Moreover, these figures indicate, by the spread of the dots representing the syntheses, that our model reproduces, to some extent, the diversity in the samples of the original distributions (suggesting a certain entropy in our model).

To further reduce the distances between the model samples and the original samples of a process, while maintaining a similar diversity (and hence a similar entropy between the model distribution and the original distribution) remains an interesting problem for future works.

The two plots in the second line present the TDA of the three Cox distributions together, and the three turbulent distributions together. We observe that, for the Cox distributions, the Small circles model is about as close to the original Small circles

The function K µ u can be complex valued. However, as our energy function is the square Euclidean norm, it is equivalent to consider that K µ u has values in R 2d . Moreover, we assume in what follows that the function K is such that for all µ ∈ M s and all u ∈ Supp(µ), K µ u is differentiable. We can then define from chain rule, for any µ ∈ M s and any u ∈ W s

and

where Jac[f ] denotes the Jacobian matrix of the function f . When u ∈ Supp(µ), the chain-rule gives Jac

). Additionally, for φ 0 ∈ M, u ∈ φ 0 , and any f ∈ F(R 2 ), we define the functions

and

For any u ∈ φ0 , and any

The result follows from the fact that, for each n ≥ 0, the function u -→ f n (u) = f n (u, φ0 ) defined in (5.16) satisfies

and, by consequence, the same holds for the composition

Then for any Γ ∈ M s , as Φ0 is invariant to the action of T :

where the second to last equality comes from the invariance of the distribution of Φ0 to the action of T . It remains to prove (5.19), which can be done by induction. It is trivially true for n = 0, as f 0 = Id. Assume now that (5. [START_REF] Błaszczyszyn | How user throughput depends on the traffic demand in large cellular networks[END_REF]) is satisfied for all f k , k = 1, . . . , n and, consequently, for

Let us fix, without loss of generality, ζ = η = 1, then

(5.20)

, then the gradients are null, so the result is trivial.

If u ∈ F n (•, φ0 ) # φ0 , then there exists u i ∈ φ0 such that u = F n (u i , φ0 ). We study the second and third term of the r.h.s. of eq. (5.20) separately. Let us note A the linear component of T , i.e. T u

In what follows, we shall simplify the notation replacing F n (•, φ0 ) # φ0 by F n and F n (u i , φ0 ) by F n (u i ).

We can write, as K is invariant to T ,

with

where we used the fact that A -1 = A t . Furthermore, using the fact that K is invariant to T ,

(5.23)

Proof. Let z ∈ C, and recall from eq. (6.4), that ρ α (z) = ρ(Real(ze iα )). Note that we have the following relation

We can then write

Replacing z with x ψ j,θ (u), z with x ψ j ,θ (u -τ ), and injecting this relation in eq. ( 6.1) gives us the desired result. This shows that using only four α uniformly chosen between [0, 2π] is sufficient to capture second order statistics. Because the wavelet transform is an invertible linear operator (on its range space), computing the r.h.s of eq. ( 6.6) for all (j, θ, τ ), as well as the low-pass coefficients carried out by Φ J , is equivalent to computing the correlation matrix of x.

Relation with the RF model Setting aside the subtraction by the spatial mean µ γ , the RF model can be viewed as a particular case of models defined by eq. (6.5). Indeed, the statistics of the RF model take the form of eq. ( 6.1), with

where {ψ f } being a family of multi-scale random filters (samples from a uniform distribution, see [START_REF] Ustyuzhaninov | What does it take to generate natural textures?[END_REF] for more details). By writing ρ α (x ψ j,θ (u -τ )) = ρ(x Real(ψ τ j,θ e iα )(u)), with ψ τ j,θ denoting the translation of ψ j,θ by τ , we see that the models are similar, the difference being that our models use wavelet-based filters instead of random ones.

Defining an appropriate Υ

The choice of the covariance set Υ is of central importance in the definition of the model. Intuitively, a too small set of indices will induce a model that could miss important structural information about the texture that we want to synthesize. Conversely, if Υ contains too many indices, the syntheses can have good visual quality, but the statistics of the model may have a large variance, leading to the memorization of some patterns of the observation. There is a trade-off between these two aspects: one must capture enough information to get syntheses of good visual quality, but not much, so as not to reproduce parts of the original image. To illustrate this point, we shall study the model ALPHA defined with three different sets Υ : A smaller model ALPHA S with a limited amount of elements in Υ, an intermediate model ALPHA I , and a larger model ALPHA L .

To precisely define these models, let us note

, from which the spatial shift shall be selected. Table 6.1 summarizes the conditions that all parameters have to satisfy to be contained in these sets. Additionally, these models include large scale information through the covariance of a low-pass filter, i.e. the spatial average of x φ J (•)x φ J (• -τ ), for τ ∈ T. To count Observation ALPHA S (3.5k) ALPHA I (35k) ALPHA L (142k) Figure 6.2: Examples of syntheses from the ALPHA models defined in Table 6.1. Visually similar image patches in textures are highlighted by red circles. The number of statistics is given in brackets next to each model name. From top to bottom: a Julesz counterexample, a stationary texture image, a stationary turbulent field and a non-stationary image. We see that the model ALPHA I achieves a balance between the visual quality and diversity on these examples.

Modelling color interactions

In order to generate color images, the covariance model ALPHA I defined in Section 6.3.2 could be directly applied to each R, G and B color channel independently. However, it would not take into account the color coherence in the structures of the observation.

To capture color interactions in the observation image, we shall impose the covariance between the coefficients of (6.5) for all indices in Υ and all color channels. More precisely, let x = {x c } c=1,2,3 be a color image, with the parameter c representing the color channel. The ALPHA C color model is defined by correlations between coefficients of the form:

The set of indices is defined as

Reduced ALPHA C The model ALPHA C has a large size as it computes correlations between all coefficients for all color channels. This size can be significantly reduced by

Numerical experiments

In this section, we compare our intermediate model to the state-of-art models (PS, RF and VGG) on both gray-scale and color textures. We first specify the experimental setup and detail the size of each model. We then present the synthesis results of various