This PhD tackles two aspects of surface reconstruction from point clouds. Firstly, it addresses the case of large-scale point clouds that are too massive to fit into the memory of one single machine. We present an end-to-end distributed algorithm that can scale up to arbitrarily large point clouds whilst guaranteeing the watertightness of the resulting surface. Secondly, we deal with the assessment of surface reconstruction by proposing two evaluation protocols. The first one requires synthetic data whereas the second one can be set up with data directly obtained from a sensor. These protocols and their respective newly-defined metrics allow the quantification of the quality of surface reconstructions with less bias than previous approaches.

Résumé

Cette thèse de doctorat traite de deux aspects de la reconstruction de surface à partir de nuage de points. Premièrement, elle aborde le cas large échelle où un nuage de points est trop volumineux pour être stocké dans la mémoire d'un seul ordinateur. Nous présentons un algorithme distribué de bout en bout permettant de traiter des nuages de points arbitrairement grands tout en garantissant l'étanchéité de la surface produite. Deuxièmement, cette thèse contribue à l'évaluation de la reconstruction de surface de par la définition de deux protocoles. Le premier nécessite des données synthétiques alors que le deuxième peut être mis en place en ayant uniquement recours à des données provenant de capteurs. Ces protocoles et les nouvelles métriques qui leur sont associées permettent de quantifier la qualité des reconstructions avec un biais moins important que les approches utilisées jusqu'alors.

rent Caraffa. I felt (and I still do!) very lucky to have them as my guides throughout these three years. Laurent helped me settle in at the lab, suggesting activities and giving me tips in order to make things more enjoyable. He was also my cornerman and only supporter at a university boxing competition. I won't forget about that. Bruno helped a lot with establishing project goals, and gave tip top instructions throughout. Our interactions were insightful and I always felt very respected. He is a very kind person who knew when to take more time than he'd planned to provide me with knowledge or support. It was very inspiring to see how he managed my PhD. I learned a lot about human relationships from him and in particular how to be a good supervisor. I feel particularly lucky about having been allowed to give classes on a subject that had nothing to do with my PhD but on which I was keen to gain experience. He also allowed (and encouraged, actually!) me to go to Southampton to work exclusively on a side project for a whole month. Lastly, he's the one who managed to find the right words to cheer me up when I felt under the weather about work or my motivation. Thank you so much Bruno.

A PhD is also about doing some admin work from time to time and as French people, we are quite the specialists in terms of making it hard...! That's why I have to thank Sylvie Cach and Mariam Sidibé, who were very supportive and always provided me with quick and clear answers.

Thank you both for dealing and having dealt with this many people who, for some reason, find computer programming easier than filling in a form and sending it on time...! I cannot forget the tremendous help I received from Alain Sombris, probably one of the three most important employees at the IGN. This guy can solve any kind of problem to the point of making French administration seem easy to follow!! Thanks for your kindness and patience, Alain! Mathieu Brédif contributed to one of the papers that we wrote during my PhD and his pieces of advice were very useful. In addition, I greatly appreciated the discussion we had during my IGN annual interview.

Jean-Pierre Papelard helped me use the stationary LiDAR station at the ENSG and Chapter 4 would not be part of this thesis without him. Thank you very much for the time you spent with me, answering my questions and guiding me throughout both the acquisition and the post-processing work.

I had very few interactions with Emmanuel Clédat (maybe two or three) but wow, what a dynamic guy! Man, you're so passionate and spread so much positive energy that talking to you always gave me a big boost! Thanks for the work you did with Agisoft Metashape.

None of this would have been possible without two teachers I will never forget. Olivier Pichon is a heck of a physics teacher! After he dropped a glue stick and told us that physics was about wondering whether the reason why this object hit the floor was because it was yellow, I was sold! Physics was something I wanted to know about! Hope to see you again soon at "Les poteaux carrés".

Pascale Vérant was my physics teacher in second year of "prépa" and she is the most supportive teacher I've ever had. She always kept an eye on me and kept me on my toes when I was being lazy, dreaming about anything but physics. That second year was really tough for me but she didn't let me drown. Thank you so much. I met Monika Csikos in September 2022 when I signed up to teach "Algorithm Design" at the ESIEE. She was the PhD student in charge of the module and she provided me with all the support I needed to take on this new challenge! I want to give a special thanks to Cédric Demonceaux who introduced me to the world of computer vision. I was initially seeking an internship in image processing "just to get that skill I lacked to be a real robotics engineer". Cédric hired me and there I was! In Le Creusot! What an experience! Talking about 3D convolution in the morning and going for a run at lunchtime was just the best thing! I swear I miss those hills. Cédric also encouraged Bruno to hire me for the PhD although he wanted us to carry on working together... great man! In March 2022, I went to Southampton (UK) for a secondment at the Ordnance Survey, the Bri'ish national mapping agency. I had the best time there and this is in large part thanks to the very warm welcome from David Holland. David introduced me to everyone and made me feel like a real OS employee, leading me to occasionally forget about being French! I enjoyed the discussions at lunchtimes, our walks around the building, the PieFull Tower at Pie Caramba and our bird watching session on a sunny Saturday mornin'! Top bloke! I spent most of my time at the IGN sharing an office with Damien Robert and Romain Loiseau. During covid-related lockdowns and curfews, Damien was one of the very few people I got to see and I am very grateful for those moments we spent together. We were not that many to come to the office every day and I felt very privileged to be able to talk with such a nice guy! Romain started to come more regularly later in the year and I enjoyed insightful discussions about very different topics including religion, living abroad and international conflicts. I then had the privilege to be invited to his wedding! Damien and Romain offered me a coffee pack for my birthday and what I can tell you is that the cafetiere has been used at least as many times as there are days since I turned 25! I can strongly recommend the IGN as a great place to work and this is partly thanks to the research group's team of directors. Sidonie Christophe, Ana-Maria Raimond and Clément Mallet are very inclusive and their status does not prevent them from being modest and funny! Clément did not even show his happiness when AJ Auxerre won their promotion match over Saint-Étienne, sealing our relegation to Ligue 2... very thoughtful. I remember Sidonie was one of the first ones to welcome me and make me feel home at the IGN. She has been (and is still) dealing with an unfair situation, being partially excluded from her working environment for being the victim of assault by a colleague. She showed me what it means to be brave and keep battling. Never back down.

Over the last three years I gave quite a few classes at the ENSG, mainly in image processing, and I want to address a big thanks to Marc Poupée for all the support he gave me, and for trusting me. He is a very dedicated teacher and students are lucky to have him. He gave his best during the covid pandemic dark times, always seeking better ways of delivering classes. Universities need more people like you! Arnaud Le Bris also helped me with a class we gave together, providing me with everything I needed to feel confident in front of the students.

Thanks to Helena for some nice discussions about Brazil and French culture and for keeping the office tidy. Working in a clean environment definitely helps staying motivated! Since working for three years in a research group is also about sharing informal moments with the other people, I'd like to say how much I enjoyed spending time with some of them. Ameur Zaïbi spent three months at the IGN in autumn 2020 and we started sharing very nice coffee breaks together. We became friends pretty quickly and we were soon able to have deep and meaningful conversations, which I enjoyed a lot. Thanks for everything my friend and good luck. Katya joined the team (which was pretty small at that time) in January 2021 and we started having coffee together everyday after lunch. That was the beginning of post-lunch coffee breaks! It's become such a tradition now! We had a lot of fun and I really needed that social interaction fix back in the days when I remember being sad on Friday evenings and happy on Monday mornings... Thank you for every minute you spent with me, for taking me to the music festival in June 2021, that meant a lot to me.

Cheers La Boss!

Anatol Garioud was the first pal from work I started considering as a friend. For six months I enjoyed our Monday lunchtime debriefs of football matches, and going to that Champion's League match (Paris-SG vs. Dortmund) at your bar was so cool!!! My friends still make fun of me though, saying that I turned into some fan of Neymar or something... I feel honoured to have met Paul Chapron (probably at the coffee room) and man, I feel like I could listen to you for years. You have so much to teach me, and it is so cool to talk with you! Once we realised we both adored the Senseï Hugo TSR, I knew we had some good shit to share with each other! I'm not sure we'll ever get to record that lyrics-based interview but I enjoyed the time we spent writing the script, which we ending up losing because the shared file expired! Stay the same, m8! Florent Geniet is definitely the one who made me laugh the most! Imagine a fella who always manages to sneak in unexpected jokes... Well it's impossible, forget about it! But still, imagine!! Always positive, I loved spending time with you, Dora! From our very first interaction, I knew that I would enjoy spending time with Lâmân Lelégard. He's a very likeable person who does not follow the crowd. He is very discrete and modest about what he does best so do not expect him to brag about how well he can draw (I learned that from someone else) or about the last paper he got accepted. Be ready, however, for long explanations about jazz music from the 60's or about how amazing coffee tastes when it is prepared properly. But don't get me wrong, this dude has a lot to teach you, a LOT! There are so many other people in the lab who I've enjoyed spending time with over the last three years. Thanks to Mounia Ahmedi, Luc Beaudoux, Emile Blettery, Charles Villard, Teng Wu, Lulin Zhang, Melvin Hersent, Grégoire Grzeczkowicz, Samuel Mermet, Solenn Tual, Sami Guembour and Mohamed Ali Chebbi to name but a few. I want to whack Helen Mair Rawsthorne in the spotlight for a million reasons. First, as she was my Bri'ish counterpart at the IGN during my secondment at the Ordnance Survey, she proofread my VOLTA project report, correcting a shit-ton of mistakes and explaining them to me. I feel very grateful that she also proofread this thesis. Not only did she help make the manuscript clearer and easier to read but she also taught me a fuchin' lo'! Well, I'm sure you learned some stuff yourself, you Frenchie! Didn't you no'?? The fact that you pulled out all the stops to help me with my English is only the tip of the iceberg. From the moment we met, I felt really excited about talking to you and getting to know you! Having you around has just made my life so much more enjoyable, you have literally brought light into my life. You have made me feel at home. I feel so lucky that our paths crossed, every day with you is a treat. Wow, what a lucky duck I am. Thank you for making every negative aspect of life so much more manageable. Thank you for making me so happy, for being so supportive, for being my number 1 fan, you are the most wonderful person in the whole wide world!!! Caru ti Helen.

We're getting to the end of the list, so it's time for me to speak a bit about those who have been there for me since way before this all started, the ones I can rely on no matter what, the ones I want to be able to spend time with for as long as possible.

Massive thanks to my friends with whom I grew up and who are still there today. I have met a lot of people in my life but I have made very few true friends. Thanks to Jérémy Da Silva, Florent Robert, Maxime Malmezac, you are my brothers. Thanks to Théo Bejannin, with whom I had a lot of fun in Montreuil. I was sad when you left, buddy! I didn't make a lot of friends during my studies but Nabil Mahnane and Guillaume Geoffroy are definitely two who are important to me. Thanks to Alexandre Legrand, you are a truly fantastic and lovable person. Thanks to Benoît Jourjon for all these years [START_REF] Lafarge | Surface reconstruction through point set structuring[END_REF][START_REF] Caraffa | 3D Watertight Mesh Generation with Uncertainties from Ubiquitous Data[END_REF]

3.1

In this situation, the real surface has been partially scanned, resulting in a serious occlusion in the input point cloud. If we uniformly sampled the entirety of both reconstructed meshes R 1 and R 2 and used the Metro tool [1] to compute the distances from them to the ground-truth, the metrics would be dominated by occlusion-based errors (highlighted in bold). . . . . . . .

3.2

Visualisation of the computation of the best expected surface for a small value of α. Starting from the input point cloud, the ground-truth mesh and two reconstructed meshes R 1 and R 2 , we compute the best expected surface and the relevant parts to assess R α 1 and R α 2 by removing all triangles lying further than α from a point of the input point cloud. . . . . . . . . . A given real laser ray (the thick one) was cast from O 1 and it hit the real surface at Gt. Note that the position of the intersection might be noisy, hence the shift between the real surface and the high-quality point cloud. We compute all intersections between the associated virtual ray (i.e. the extension of the real laser ray) and the reconstructed surface. In this case, it results in six intersections I 1,...,6 . The closest intersection to Gt happens to be I 2 so the "ray distance" metric for that particular ray is the distance (Gt, I 2 ). If (Gt, I 2 ) < d max , Gt is to be counted as a TP, otherwise it will not be taken into account in the evaluation since it is situated after Gt. Note that this piece of surface might still be evaluated thanks to another ray (as one emanating from O 2 , for example

Virtual environment from

List of Tables

Publications International Conferences

Introduction Motivation

Surface reconstruction is the task of producing a continuous digital representation of a real surface about which discrete information has been acquired. This information may already be in the form of point clouds produced by a laser scanner. This includes Timeof-flight [2] and Structured-light [3] devices as well as terrestrial and airborne LiDAR [4] that can scan large environments. Point clouds can also be produced from images using Multi-view stereo [START_REF] Furukawa | Multi-view stereo: A tutorial[END_REF] or Structure from motion [START_REF] Ozyesil | A survey of structure from motion[END_REF].

This task has been extensively studied and a tremendous amount of approaches have been proposed as shown in Chapter 1. Nonetheless, two major aspects are still facing a lack of contributions.

Firstly, very few approaches efficiently handle large point clouds as shown in Chapter 1, Section 1.2. With the improvement in sensor capabilities, it has become much easier to acquire massive amounts of visual data in the last few years. In particular, cameras are able to capture more pixels, LiDAR scanners produce more points, and the platforms onto which these device can be embedded have diversified: static, terrestrial and aerial mobile mapping or drones.

As pointed out in [START_REF] Cignoni | External memory management and simplification of huge meshes[END_REF], random-access memory (RAM) size is often a bottleneck when it comes to dealing with large amounts of data. Most surface reconstruction algorithms build the surface all at once by loading the whole point cloud and then creating additional data structures (Delaunay triangulations, voxel grids, octrees, etc.) to represent the output surface. This implies a strong constraint on the maximum input data size for a given amount of RAM. Processing time inevitably increases with the number of points to process, at best in a linearly fashion. Even without considering real-time applications, it is important to keep the processing time within reasonable bounds when the data size becomes very large.

The second major aspect of surface reconstruction that has not yet been deeply explored is its evaluation in the case of open scenes. While many methods are available to build a 3-dimensional mesh from point clouds, there is a lack of relevant metrics to assess the accuracy of these reconstructions. As surface reconstruction is an ill-posed problem, several outputs are possible for the same input point cloud, and they might not even be equivalent in terms of topology (e.g. number of components). How can we know which model is the best? Can we quantify their accuracy? Although it is possible to assess the quality of the models visually, this approach raises several issues. First, it is a subjective comparison and one might be tempted to favour their own or preferred method over others. Secondly, everyone has a different perception of visual quality making it difficult to come to an agreement, even in the absence of a conflict of interest. Thirdly, while it might sound reasonable to visually evaluate a few different models of a relatively small scene, it is very unlikely that one would be able to carry out a large-scale evaluation involving dozens of models representing large areas. Current metrics are not entirely satisfactory as they often depend heavily on the priors used by the algorithms to reconstruct a surface (closed or not, for example) or the ground truth density. This becomes particularly troublesome in the case of open scenes when the ground truth might not be dense and homogeneous, and when the priors have a massive influence on the output surface. It is thus essential to find relevant metrics.

In the context of remote sensing, this topic was first neglected in favour of 2.5D approaches, where the surface reconstruction problem is merely a question of interpolation of possibly sparse height data sampled on a regular grid. This led to the popular Digital Elevation Models (DEMs) and Digital Surface Models (DSMs) used to represent the geometry of the visible surface of a scene seen from above. This representation is, however, becoming more and more limited: an increasing number of applications require terrestrial data (Mobile Mapping Systems, fixed stations and hand-held cameras) for which the 2.5D setting is completely inappropriate. Moreover, the spread of aerial oblique imagery aiming at acquiring more data on vertical surfaces, and the fact that some aerial LiDAR can scan up to 40 • away from the vertical, call for more generic 3D models in order to have a continuous geometric representation of the underlying scene. Finally, an increasing number of commercial products such as Sure by NFrames or ContextCapture by Bentley already offer full 3D processing pipelines for surface mesh reconstruction from remote sensing data. For all these reasons, 3D surface reconstruction, once a topic mainly studied in the geometry processing community, is becoming more and more widespread in remote sensing.

Why has so little work been carried out on the assessment of surface reconstruction of open scenes? One major reason is the complexity of defining an appropriate ground truth. In real-case scenarios when the goal is to produce a digital model of a real object or scene, there is no ground truth other than the real surface itself. It is thus impossible to directly compute the "distance" or the "difference" between a digital model and the ideal real surface. One option is to use synthetic data as in Chapter 3, where a given mesh is chosen to be the ground truth and is then virtually scanned in a way that simulates the defects of a real acquisition, and then the surface reconstruction algorithms to be evaluated are run on this virtual scan. This makes it more straightforward to compute metrics that assess the difference between the ground truth model and the reconstructed one. Another possibility for working with data representing real scenes is to scan the given scene or object chosen for the evaluation in two different ways such that we can be certain that one scan is significantly "better" than the other one. In particular, this high-quality point cloud has to be denser, more accurate and/or contain fewer occlusions. After running the algorithms to be assessed on the low-quality point cloud, the resulting meshes can be evaluated using the additional information contained in the high-quality one.

Definitions and properties of surfaces

A surface is not defined the same way in all branches of mathematics. In topology [START_REF] Kinsey | Surfaces[END_REF], a surface is a two-dimensional manifold (2-manifold). Definition 0.0.1. An n-manifold is a topological space such that each point has a neighbourhood that is homeomorphic to an n-dimensional open disc. Definition 0.0.2. An n-manifold with boundary is a topological space such that each point has a neighbourhood that is homeomorphic to either an n-dimensional open disc or the half-disc.

We can thus define the boundary of a manifold as the set of points with half-disc neighbourhoods.

In computer graphics, a triangle mesh is the most common way of representing surfaces, so let us bind the two domains by introducing the concepts that lead to the definition of a triangulated surface. Definition 0.0.3. An n-cell is an n-manifold with boundary with the additional property that its boundary must be divided into a finite number of lower-dimensional cells, called the faces of the n-cell.

In particular:

• A 0-dimensional cell is a point • A 1-dimensional cell is a segment • A 2-dimensional
cell is a polygon (we will consider it to be a triangle)

• A 3-dimensional cell is a solid polyhedron (we will consider it to be a tetrahedron) Cells can be bound to each other to form complexes: Definition 0.0.4. A complex K is a finite set of cells K = {σ : σ is a cell} such that:

1. if σ is a cell in K, then all faces of σ are elements of K.

2. if σ and τ are cells in K, then Int(σ) ∩ Int(τ ) = ∅.

The dimension of K is the dimension of its highest-dimensional cell. Definition 0.0.5. A simplicial complex (or triangulation) is a 2-complex structure with only triangular cells, and with the additional condition that any two triangles are either connected along a single edge or at a single vertex, or are disjoint.

In all following chapters, we will always refer to triangle meshes as surfaces. Before extending the meaning of this term, let us properly define what we mean by a triangulated surface. Definition 0.0.6. A triangulated surface (without boundary) is a simplicial 2-complex such that:

1. each edge belongs to exactly 2 triangles.

2. the triangles that meet at a vertex can be labelled T 1 , T 2 , . . . , T n , with adjacent triangles in this sequence connected along an edge and T n connected to T 1 along an edge.

A triangulated surface with boundary can, however, have edges that belong to only one triangle. The boundary is formed by such edges.

In the domain of computer graphics, one may encounter triangle soups (a set of triangles in R 3 ) that are not triangulated surfaces. It is possible to define a list of vertices with connections between them without the resulting triangle soup respecting the criteria for being a surface. We will thus define some of the properties that triangle soups can have.

Edge-manifold if, for each edge, the set of triangles that share this edge form a topological disk.

Vertex-manifold if, for each vertex, the set of triangles that share this vertex form a topological disk.

Manifold if any point on the mesh has a neighbourhood that is homeomorphic to the Euclidean plane R 2 . This is desirable as many mesh-processing algorithms will fail on non-manifold meshes.

Orientable if one can define a consistent continuous orientation of each triangle, meaning that each edge shared by two triangles appears in the opposite order in the index order of each triangle. Once again, mesh processing methods often require orientability, Figure 1: 2D comparison between hard watertightness (left) and soft watertightness (right). The soft-watertight surface has a border (materialised by the blue dots) but only on the boundary of the domain. When trying to model urban environments, the red piece of surface from the hard-watertight surface does not have any significant meaning so soft-watertightness is often better suited to open scenes.

in particular because the widely-used halfedge data structure has this prerequisite.

Intersection-free if all pairs of triangles that do not share an edge or vertex do not intersect. Self-intersection is always seen as a bad property as real surfaces cannot self intersect.

Watertightness: A surface is watertight if it has no border. In the case of a triangle mesh, this means each edge needs to have exactly two incident faces. We will call this property hard watertightness. When trying to reconstruct an open scene, it is often more realistic to authorise the reconstructed surface to intersect the boundary of a domain of interest B (a bounding box, for instance) as illustrated in Figure 1. We thus define soft watertightness as the property that a mesh has when it has no border except on the boundary of the domain ∂B. Only triangle edges lying on the domain boundary can have only one incident face.

Objectives

This thesis aims to address surface reconstruction through the following objectives:

1. Scaling up surface reconstruction by proposing an algorithm that can process arbitrarilylarge point clouds.

2. Defining new metrics to evaluate surface reconstruction.

Main Contributions

The main contributions of this thesis are as follows:

1. An end-to-end distributed out-of-core algorithm guaranteeing watertightness and suited to any kind of large-scale point cloud with an arbitrarily-complex 3D geometry.

2. An evaluation protocol to assess surface reconstruction in the case where a ground truth mesh is available. For open scenes, this is only possible by synthetically scanning a realistic mesh, in our case publicly available 3D meshes of open scenes. We thus developed an aerial LiDAR simulator with different scanning patterns and realistic noise, the code of which is open source and available online1 .

3. An evaluation protocol to assess surface reconstruction using only real data, which we applied to public datasets but also to a dataset specifically designed and acquired for this study. This dataset and the code to evaluate a reconstruction has been publicly released 1 .

Outline

The dissertation is composed of 4 chapters.

Chapter 1 presents a state of the art of the relevant contributions. In particular, we present surface reconstruction algorithms but also general tools to process massive input data. Moreover, we investigate previous work on the evaluation of surface reconstruction.

In chapter 2, we present a distributed and out-of-core algorithm that is able to process large-scale point clouds. We explain how we managed to parallelise the last step of a reference algorithm, which until now was its main bottleneck.

Chapter 3 presents an evaluation protocol using synthetic data as well as some tools we had to develop to carry out such an assessment, in particular an aerial LiDAR simulator and a high-quality mesh dataset.

Another approach to assess surface reconstruction without the need for synthetic data is proposed in Chapter 4.

Finally, we give the conclusion and perspectives of our work.

Chapter 1

State of the art 1.1 Surface Reconstruction

Here we review existing methods to reconstruct a triangle mesh from point cloud and classify them by the paradigm they use. See [START_REF] Berger | A Survey of Surface Reconstruction from Point Clouds, Comput[END_REF] and [START_REF] Khatamian | Survey on 3D Surface Reconstruction[END_REF] for an even deeper analysis of them (even though the most recent methods do not appear). Methods evaluated in Chapters 3 and 4 are typesetted in bold.

Indicator function

Often used to achieve watertight reconstructions, this class of algorithms proceed by computing a space segmentation. The object itself is defined as the region of space where the labelling equals a certain value. The surface is then computed by finding the changes in the segmentation. A popular approach in this field is Poisson reconstruction [START_REF] Kazhdan | Poisson surface reconstruction[END_REF]. Their indicator function χ is defined as 1 inside the object and 0 outside. They show that χ convolved with a smoothing filter has to respect Poisson equation (1.1) where -→ V is a vector field depending on point locations and the associated normals.

∆ χ = ∇. - → V (1.1)
This differential equation is solved numerically and an adaptation of the marching cube algorithm [START_REF] Lorensen | Marching Cubes: A High Resolution 3D Surface Construction Algorithm[END_REF] is used to extract a triangle mesh approximating the χ = γ isosurface, where γ is the average of χ at the sample positions. This approach is screened in [START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF] to incorporate additional constraints on sample locations, which significantly improves the resulting quality. This implementation also supports two boundary conditions:

• Dirichlet specifies the values that χ needs to take along the boundary of the domain ∂M . Hard watertightness is then enforced by imposing χ = 0 along ∂M .

• Neumann specifies the values that ∇χ needs to take along ∂M . While this boundary condition also allows hard watertightness, it is less restrictive because it enables the surface to cross ∂M orthogonally, thereby only imposing soft watertightness.

[14] uses the same input data but their algorithm computes the Fourier coefficients of the indicator function χ, and extracts the surface approximating its isovalue. A differentiable Poisson solver has been introduced in [START_REF] Peng | Shape As Points: A Differentiable Poisson Solver[END_REF] to further increase the robustness to outliers and to serve as a basis for a trainable model to reconstruct surfaces. They predict the effect that shifting the oriented point samples has on the resulting surface in order to find the surface that best approximates the input point cloud. To do so, they use a bi-directional L2 Chamfer Distance L CD between the mesh and the point cloud and iteratively shift points to optimise L CD . They also regularly update the input point cloud by sampling the largest mesh component with the objective of discarding outliers.

Recently, lots of learning-based methods have been proposed. While they often outperform non-learning-based ones on simple geometries, especially closed objects, they have not been proved as able to deal with the complexity of large and open scenes. IM-NET [START_REF] Chen | Learning implicit fields for generative shape modeling[END_REF] is a learning framework which predicts whether any point (x, y, z) is inside or outside the given shape needing to be reconstructed. The input of their network is the 3 coordinates of a point as well as a feature vector that can be computed using PointNET [START_REF] Qi | Pointnet: Deep learning on point sets for 3d classification and segmentation[END_REF].

Occupancy Networks [START_REF] Mescheder | Occupancy networks: Learning 3d reconstruction in function space[END_REF] presents a similar way of computing the so-called occupancy function of the 3D object. However, instead of concatenating a feature vector to the coordinates of points, they use a batch-sampling strategy. An important advantage of the two latter methods is the arbitrary resolution at which the surface can be extracted. [START_REF] Peng | Convolutional occupancy networks[END_REF] introduced a learning-based framework to compute implicit surfaces. Taking as input a point cloud x, their encoder computes features which can be projected in three ways, resulting in a 2D or 3D representation. Given a location p ∈ R 3 , linear interpolation is applied to compute its feature vector ψ(x, p). This vector is known at point cloud locations, allowing for interpolation. Occupancy probability at p is then predicted thanks to a fully-connected network f θ (p, ψ(x, p)) ∈ [0, 1]. A mesh can be extracted by the application of Multiresolution IsoSurface Extraction (MISE) [START_REF] Mescheder | Occupancy networks: Learning 3d reconstruction in function space[END_REF]. A sliding-window implementation enables the algorithm to process large-scale inputs.

Convolutional Occupancy Networks

Recently, [START_REF] Groueix | A Papier-Mâché Approach to Learning 3D Surface Generation[END_REF] proposed a general learning framework dubbed AtlasNet to take as input a 3D point cloud or an RGB image. It proceeds by concatenating this data with a sampling of a patch, namely the unit square, before passing it to multilayer perceptrons (MLPs) with rectified linear unit (ReLU) nonlinearities, producing as output a point cloud of arbitrary resolution. A mesh can be generated either by transferring the connections between vertices of a mesh defined on the patch to their 3D image points or by using Poisson Surface Reconstruction [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] on a sufficiently dense point cloud. A third solution is to sample a 3D sphere instead of patches.

Volumetric Segmentation

This is a sub-discipline of indicator functions as it consists in giving information about whether a region of space is filled by the object or is empty. The data structure can be:

• the Delaunay Triangulation of input samples as in [START_REF] Labatut | Robust and efficient surface reconstruction from range data[END_REF][START_REF] Lafarge | Surface reconstruction through point set structuring[END_REF][START_REF] Caraffa | 3D Watertight Mesh Generation with Uncertainties from Ubiquitous Data[END_REF][START_REF] Kolluri | Spectral Surface Reconstruction from Noisy Point Clouds[END_REF].

• voxels: [START_REF] Holenstein | Watertight surface reconstruction of caves from 3D laser data[END_REF] labels them as free space, occupied or unknown. To achieve this, point locations combined with sensor positions allow computing the ray corresponding to a beam of free space. An interesting feature is that undesirable moving objects such as humans can be erased in the final surface thanks to scans of the same area from different sensor positions.

Robust and Efficient Surface Reconstruction (RESR) [START_REF] Labatut | Robust and efficient surface reconstruction from range data[END_REF] label as inside or outside each tetrahedron of the Delaunay triangulation of the point samples. The triangles separating an empty-labelled tetrahedron from an occupied-labelled one are extracted thanks to a graph-cut optimisation of an energy function defined thanks to the lines of sight (emanating from the vertex and pointing at the laser scanner) and the shape of the triangles. Similarly, Watertight Mesh generation With Uncertainties (WMWU) [START_REF] Caraffa | 3D Watertight Mesh Generation with Uncertainties from Ubiquitous Data[END_REF] label as occupied or empty each tetrahedron t of the Delaunay triangulation T of the point samples. First of all, a set of mass functions m t are computed. Each m t represents the likelihood of tetrahedron t ∈ T being empty or occupied, or the case in which its occupancy is mostly unknown. These are computed based on visibility priors that make use of sensor positions. Binary labels l t are attributed to each tetrahedron t ∈ T minimising an energy function (equation 1.2). Denoting l T = (l t ) t∈T as the labelling of each tetrahedron in the triangulation T and L = {0, 1} the set of possible labels, the problem is formulated as:

l T = arg min l T ∈L T E data (l T ) + λ E prior (l T ) (1.2) E data (l T ) = t∈T ∥l t -m t ∥ 2 (1.3) E prior (l T ) = (t 1 ,t 2 )∈T 2 |t 1 ∩ t 2 |.∥l t 1 -l t 2 ∥ 2 (1.4)
The data term (equation 1.3) enforces a solution close to the overall mass function and the prior term (equation 1.4) minimises the total area of the interfaces. The value of the parameter λ balances the two of them. [START_REF] Sulzer | Scalable Surface Reconstruction with Delaunay-Graph Neural Networks[END_REF] also use this paradigm (equation 1.2) to label Delaunay cells of the input point cloud tetrahedralisation but they claim that the problem with previous methods like RESR [START_REF] Labatut | Robust and efficient surface reconstruction from range data[END_REF] and WMWU [START_REF] Caraffa | 3D Watertight Mesh Generation with Uncertainties from Ubiquitous Data[END_REF] is that the occupancy of the cells is often estimated using visibility models that are not necessarily robust to noise. In this context, the main asset of their method is the estimation of the occupancy of the tetrahedra thanks to a graph neural network. They can thus adapt to different complexities by learning on different point cloud configurations. A graph-cutbased optimisation then reinforces global consistency.

Delaunay-Graph Neural Networks (DGNN)

Signed-distance function

Another way of generating a watertight surface is to compute the signed-distance function f to the surface and to extract its zero-level set. This is the approach chosen in Multi-level Partition of Unity (MPU) [START_REF] Ohtake | Multi-level partition of unity implicits[END_REF] and Smooth Signed Distance (SSD) [START_REF] Calakli | SSD: Smooth signed distance surface reconstruction[END_REF].

MPU [START_REF] Ohtake | Multi-level partition of unity implicits[END_REF] presents an implicit surface reconstruction algorithm which takes as input a normal-equipped point cloud and produces an approximation f : R 3 -→ R of the signeddistance field from the surface. The surface is then approximated by the zero-level set of this function. The space delimiting the region of interest is subdivided using an adaptive octree: if the accuracy of the shape function within a certain cell is not acceptable, the region is subdivided and a new computation is run. This quality-based behaviour offers the possibility of a trade-off between quality preservation (including sharp features) and computation performance (computation time and memory consumption).

Local shape functions:

We denote by Ω the bounded domain of interest included in R 3 which contains the point cloud P = {p 1 , . . . , p N } equipped with the set of normals N = {n 1 , . . . , n N }. As the computation mostly happens within the octree cells, let us analyse the various calculations that will be done for a given cell of centre c and length of the main diagonal d. A support radius (equation 1.5, α being a smoothing parameter typically equal to 0.75) is associated with the cell and only the set of points within the ball B(c, R), which we will denote by P ′ , is considered for further operations.

R = αd (1.5)
For the sake of clarity, we do not detail the conditions under which the cell is supposed to be subdivided or not, let us just say the number of points within B has to be sufficient and as mentioned above, the quality of the fitting procedure needs to respect a user-specified threshold.

A local shape function Q(x) is computed based on P ′ . In order to account for the various conditions of density and distribution of the points, three different types of function are considered:

• a general 3D quadric

• a bivariate quadratic polynomial in local coordinates

• a piecewise quadric surface

Approximation of the global solution: Now that each of these local shape functions Q i have been estimated, the global signed-distance function f is computed by blending all the contributions, weighted by non-negative compactly-supported functions {w i (x)} such that the union of their supports span the whole domain Ω (equation 1.6).

f (x) ≈ i ω i (x) Q i (x) n j=1 ω j (x) (1.6)
SSD [START_REF] Calakli | SSD: Smooth signed distance surface reconstruction[END_REF] consists in using a least-square minimisation of an average of several energy functions E D i (f ), weighted by coefficients λ i (equation 1.7).

E D (f ) = i λ i E D i (f ) (1.7)
The two main energy functions used are those of equations 1.8 and 1.9.

E D 0 (f ) = 1 N N k=1 f (p k ) 2 (1.8)
E D 1 (f ) = 1 N N k=1 ||∇f (p k ) -n k || 2 (1.9)
Equation 1.8 is precisely the one enforcing the condition f (x) = 0 at point locations (near the surface). Equation 1.9 comes from the fact that the gradient of the function f represents the normal field where f (x) = 0. Consequently, near the surface, the condition ∇f (p k ) = n k should be satisfied. A third energy function involving the Hessian matrix of f enforces that the gradient of f should remain almost constant away from the surface. Following this idea, the surface produced is watertight as long as the function f is continuous.

Recently, DeepSDF [START_REF] Park | Deepsdf: Learning continuous signed distance functions for shape representation[END_REF] has shown how to learn the surface distance field. They do so by using the 3D coordinates of a point as well as a latent vector that accounts for the type of shape to which the point belongs.

Unsigned-distance function

[30] presents a method for reconstructing large-scale watertight and manifold surfaces using only point locations. They proceed by estimating a confidence map over a predefined volumetric grid V . This function ϕ : v -→ c ∈ [0, 1] associates to any voxel v its confidence c which can be seen as the pseudo-distance to the nearest point location p. The aim is then to minimise the sum of pseudo-distances over a certain set of voxels (equation 1.10). An algorithm for extracting the corresponding mesh is also proposed.

S opt = arg min S⊂V v∈S ϕ(v) (1.10)

Primitive-based

In this field, PolyFit [START_REF] Nan | Polyfit: Polygonal surface reconstruction from point clouds[END_REF] uses RANSAC [START_REF] Schnabel | Efficient RANSAC for point-cloud shape detection[END_REF] to detect planar segments and refine them.

The surface is extracted by combining the optimisation of an objective function which favours data fitting, point coverage and model complexity and the enforcement of watertightness and manifoldness.

Point Set Structuring (PSS) [START_REF] Lafarge | Surface reconstruction through point set structuring[END_REF] relies on the Delaunay triangulation of input points and the labelling of its tetrahedrons as empty or occupied, but their specificity resides in the extraction of primitives as a pre-processing step, a resampling of the resulting structures and the combination of points from planar regions and unstructured ones in the reconstruction step.

MLS-based

Moving least squares (MLS) was first introduced by Lancaster in [START_REF] Lancaster | Moving Weighted Least-Squares Methods[END_REF], based on the work conducted by, amongst others, Shepard in [START_REF] Shepard | A two-dimensional interpolation function for irregularly-spaced data[END_REF]. Since then, a tremendous amount of extensions have been added as pointed out by a survey conducted in 2008 in [START_REF] Cheng | A survey of methods for moving least squares surfaces[END_REF]. For instance, [START_REF] Levin | The Approximation Power Of Moving Least-Squares[END_REF], [START_REF] Alexa | Point set surfaces[END_REF] and [START_REF] Levin | Mesh-Independent Surface Interpolation[END_REF] significantly contributed to the advances in MLS-based algorithms. As explained in [START_REF] Cheng | A survey of methods for moving least squares surfaces[END_REF], MLS-based algorithms can be roughly classified into two main categories:

Implicit MLS surfaces require the computation of a level set function of which the zero isosurface can be extracted.

Projection MLS surfaces consist of first computing a projection operator that maps any point of the space to a point on the surface. The surface is then made of the set of stationary points. [START_REF] Amenta | Defining point-set surfaces[END_REF] defined a projection operator on points equipped with normals (that they call surfels). That enabled [START_REF] Fiorin | Out-of-core MLS Reconstruction[END_REF] to propose an out-of-core implementation of an MLS-based explicit surface reconstruction method. Their projector, applied to any point p of the space, produces another point p * and a normal n * . n * is defined as the vector that best mimics the normal of the surface patch closest to p. This vector field is computed as a weighted average on a neighbourhood of p as follows:

n P SS (p) = i - → n i ϑ N (p, x i ) (1.11)
where the weighting function is:

ϑ N (p, x i ) = e -∥p-x i ∥ 2 h 2 j e -∥ p-x j ∥ 2 h 2 (1.12)
p * is defined as the point that minimises the energy function defined by equation (1.13).

It is then the nearest point to p that belongs to the surface.

e P SS (p) = i dist(p, x i , n i ) ϑ N (p, x i ) (1.13)
The pair (p * , n * ) enables the definition of a plane P that is, by construction, the best approximation of the nearest surface to p. We are now able to define the distance between the point p and the surface: it is approximated by the distance between p and its projection onto the plane P.

All this being well-defined, the distance between a point and the surface defines the so-called scalar field. An octree is used to implement the out-of-core computing of the surface. This structure is also leveraged at this point: the scalar field is evaluated at the corner of the deepest leaves and if a sign change is detected between two corners of the same leaf, the corresponding voxel is sent to a polygonalisation algorithm to extract the local mesh.

Refinement

These methods are particularly useful when data is too massive or when it needs to be processed on the fly. For example, [START_REF] Allegre | A Dynamic Surface Reconstruction Framework for Large Unstructured Point Sets[END_REF] presents an out-of-core algorithm that interactively processes point clouds that do not fit into memory. Their approach consists in subsampling the initial input point cloud P to produce a new, smaller one: P rep . The Delaunay triangulation (DT) of P rep is computed and the geometric convection algorithm from [START_REF] Chaine | A geometric convection approach of 3-D reconstruction[END_REF] allows to reconstruct a simplified version of the surface implied by P . After dividing P into n regions of equal size such that:

P 1 ∪ P 2 ∪ • • • ∪ P n = P , points of each P i, i=1,.
..,n are inserted in the triangulation and a surface refinement algorithm processes them in order to update the reconstructed surface.

Large-scale reconstruction

As explained in the Introduction, massive point clouds are now available thanks to the rise in sensor capabilities but their processing is challenging mainly because of the increased processing time and the RAM size that is necessary in order to store this data. In this section, we present the work of authors who proposed innovative approaches to limit memory footprint or processing time or both. Work addressing surface reconstruction specifically are summarised and classified in Table 1.1.

Pioneers

Some authors have pioneered Big Data by introducing new paradigms bypassing the typical constraints we discussed. As an example, [START_REF] Hoppe | Surface reconstruction from unorganized points[END_REF] is an algorithm by Hoppe et al. dating back from 1996 that could be performed in streaming given the locality of the information that is necessary to build the surface. The ball-pivoting algorithm [START_REF] Bernardini | The ball-pivoting algorithm for surface reconstruction[END_REF] was proposed in 1999 and has come as a forerunner in terms of out-of-core surface reconstruction. It is an incremental interpolating algorithm which basic idea is the following: a ball pivots around an edge of an initial triangle until it touches another point. If no other point is in the ball, the three points form a triangle. The process is repeated for another edge and then from another seed triangle until the whole point cloud has been explored. The interesting property is that only a small neighbourhood is considered. We can therefore process the data incrementally. To achieve this, two axis-aligned planes π 1 and π 2 define the so-called active region of work for pivoting. When all edges within this space have been tested, the two planes are shifted, data that is not going to be used again is dumped and new points are loaded according to the new active region of work.

In 2005, [START_REF] Pajarola | Stream-processing points[END_REF] introduced a stream-processing model that enables to process point clouds out-of-core. The idea of the so called sweep-plane model is that some operators only need a subset of the whole input data at a given time. Thus, this small number of points is loaded in memory and the function can be performed while the rest of the point cloud stays on disk. This has been the theoretical basis for several algorithms including [START_REF] Bolitho | Multilevel Streaming for Out-of-core Surface Reconstruction[END_REF] which is still considered as a standard in this domain. More precisely, [START_REF] Pajarola | Stream-processing points[END_REF] describes the concept of streaming as a "sliding window" that defines the fraction of the input data that is living ie stored in in-core main memory and being involved in some calculation. This concept of a sweep plane moving forward along a certain axis and finding points on the way implies that the points stored on disk have to be sorted along that axis. This pre-processing step is crucial but obviously cannot be carried out in internal memory. The problem has to be re-thought since classical methods for sorting cannot be used as pointed out by Knuth in [START_REF] Knuth | The art of computer programming[END_REF]. Knuth describes in his book several several solutions to external sorting. [START_REF] Vitter | External memory algorithms and data structures: dealing with massive data[END_REF] reviews it as well as a couple improvements.

Assuming that p 1 , ... , p n ∈ R 3 is the ordered set of points that we want to process, A = {p j-m , ... , p j } is called the active set: these points are the only ones to reside in main memory. The condition that the active set is orders of magnitude smaller than the full data set, |A| = m ≪ n, ensures that we can process arbitrarily large data sets if it can fit on disk.

Two fundamental definitions are given in [START_REF] Pajarola | Stream-processing points[END_REF]. The neighbourhood N i of a point p i is often defined as its k-nearest neighbours or points lying within the sphere centred on p i and of a given radius.

Definition 1.2.1.

A local operator Φ (p i ) performs a function on a point p i that computes or updates a subset of attributes A i associated with p i . As function parameters, Φ (p i ) only accepts p i , A i and a set of points p j ∈ N i within close spatial proximity to p i (and all their associated attributes A j ).

Definition 1.2.2. A local operator Φ k (p i ) is streamable if it is computed in one single invocation on p i and not called recursively on points p j ∈ N i . Additionally, the FIFO semantic of its queue Q k ensures no interference between consecutive operators Φ k±1 .

Indeed, the framework enables to concatenate several operators Φ k each of these acting as a sequential FIFO queue buffer Q k on the point stream. In order to fulfil these requirements, operators need to be organised in a compatible order and semantic rules need to be implemented for each operator buffer based on its class.

Three operators need to be highlighted:

1. The read operator Φ R (p i ) which reads and buffers one new point. It must be the first operator of the stream.

2. The deffered-write operator Φ W (p i ) which removes and writes point p i-m to the output stream. It is the final operator.

3. The neighbourhood operator Φ X (p i ) which computes the set N i . This one can be implemented in a number of ways depending on the definition of the neighbourhood of a point. The k-nearest neighbours is probably the most famous one.

Standardised frameworks

Message Passing Interface (MPI) [START_REF] Pacheco | Parallel programming with MPI[END_REF][START_REF] Gropp | Using MPI: portable parallel programming with the message-passing interface[END_REF] is the standardisation of a set of tools enabling the communication between several processes that own only local memory. The goal is then to share the workload amongst several nodes to save time. In the surface reconstruction field, it has been used by [START_REF] Randrianarivony | Parallel Implementation of Surface Reconstruction From Noisy Samples[END_REF].

MapReduce [START_REF] Dean | MapReduce: simplified data processing on large clusters[END_REF] is a programming model specifically dedicated to processing massive datasets. Data is represented as key/value pairs onto which two types of operation can be carried out: Map and Reduce. The Map operator accepts one pair and produces a set of intermediate pairs that are grouped by key before being sent to the Reduce operator which processes the intermediate values. This paradigm has a master/slave architecture: one particular program, the master, organises the workload across all nodes of the cluster. To achieve this, it first defines the data partitions. It also attributes operations to the nodes. The workers, processes called by the master, then read the necessary data, compute the task they have been assigned and return the result to the master. In case of failure it is the master's responsibility to organise its re-run. As a consequence, a failure of the master is not tolerated.

Spark [START_REF] Zaharia | Spark: Cluster Computing with Working Sets[END_REF] is a cluster computing framework implemented in Scala programming language that uses the same master/slave architecture as MapReduce [START_REF] Dean | MapReduce: simplified data processing on large clusters[END_REF] but it is bettersuited than for two particular types of task:

• iterative applications which necessitates to loop through the same data and apply certain operators in parallel.

• interactive works where data needs to be accessed on demand.

In both these cases, Spark enables to choose a storage strategy that quickens its access and saves an important amount of time. In addition, two advantageous features of MapReduce, scalability and fault-tolerance, are preserved.

Resilient Distributed Datasets A resilient distributed dataset (RDD) [START_REF] Zaharia | Resilient Distributed Datasets: A Fault-Tolerant Abstraction for In-Memory Cluster Computing[END_REF] is an abstraction of a parallel data structure that is fault-tolerant and user-manageable. Essentially, the data is virtually partitioned across several machines (the nodes) and one can create an RDD by defining an operation on data residing in some stable storage device. The result is not necessarily computed but instead, other operations can be defined on this first RDD, resulting in the definition of new RDDs. Such operations are then called transformations and include flatMap which is the equivalent of the Map of the MapReduce paradigm. As opposed to those "lazy" operations, actions actually trigger the computation of necessary pre-defined transformations and either return a result or write data to external memory.

Fault-tolerancy is inherently related to this notion of lazy-computing pattern. Indeed, more than the data itself, its lineage is crucial because it enables any partition that has been lost due to a node failure to be retrieved. The set of transformations guarantees that any RDD can be recomputed based on the original data and the log of all the transformations it went through. Moreover, only the lost partition necessitates such a process.

Operations are run based on data locality: by default, the scheduler favours the node that already contains the necessary input data to execute the task. But on top of that, the user can influence the organisation of the work in two ways:

• persistence can be enforced that is to say one can choose to replicate an RDD. It is particularly useful if it is meant to be transformed several times because in case of failure, the recovery can start from the latest persisted RDD instead of going all the way from the beginning.

• one can choose where to store some data by using partitioning. For instance, an insightful strategy is to maintain in memory on the same nodes partitions that are going to be combined in an operation.

Methods addressing large-scale reconstruction

The refinement-based method [START_REF] Allegre | A Dynamic Surface Reconstruction Framework for Large Unstructured Point Sets[END_REF] already introduced in Section 1.1.7 enables to incrementally refine the surface by treating each of the sub point clouds separately.

The methods presented in Parallel Poisson Surface Reconstruction [START_REF] Bolitho | Parallel Poisson Surface Reconstruction[END_REF] and Streaming MLS [START_REF] Cuccuru | Fast low-memorv streaming MLS reconstruction of point-sampled surfaces[END_REF] are both out-of-core and parallel. The external memory feature of Parallel Poisson [START_REF] Bolitho | Parallel Poisson Surface Reconstruction[END_REF] comes from the fact that it is based on a previous algorithm [START_REF] Bolitho | Multilevel Streaming for Out-of-core Surface Reconstruction[END_REF] while being able to run several computations at the same time is the real contribution of their paper. Their approach consists of formulating the surface reconstruction problem as a Poisson equation as introduced in [START_REF] Kazhdan | Poisson surface reconstruction[END_REF]. However, to enable scalability, they get rid of the necessity to consider the whole point cloud at once and rather partition the space in order to solve the problem locally.

Streaming MLS [START_REF] Cuccuru | Fast low-memorv streaming MLS reconstruction of point-sampled surfaces[END_REF] is based on Moving Least Squares [START_REF] Lancaster | Moving Weighted Least-Squares Methods[END_REF] and has a very straightforward parallel implementation that enables it to be run on a cluster of PCs. Besides, [START_REF] Fiorin | Out-of-core MLS Reconstruction[END_REF] (introduced in Section 1.1.7) also presented an MLS-based out-of-core algorithm.

[57] is also based on the computation of an implicit function and the extraction of a mesh approximating its isovalue. They compute the function by using wavelets. An important result is that their streaming implementation is an order of magnitude faster than Poisson [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] and MPU [START_REF] Ohtake | Multi-level partition of unity implicits[END_REF] still using far less memory.

[51] use Non-Uniform Rational B-Spline (NURBS) approximation to compute a surface. The optimisation problem is formulated such that the knot vector T is part of the unknowns. The solution is parallelised using MPI [START_REF] Pacheco | Parallel programming with MPI[END_REF] on the Chemnitzer Linux Cluster supercomputer [START_REF]Chemnitzer Linux Cluster (CLiC), archive Location: Worldwide Last Modified: 2020-06-02 Library Catalog[END_REF] and the associated speedup is almost linear. Recently, [START_REF] Han | Distributed surface reconstruction from point cloud for city-scale scenes[END_REF] proposed an extension of a segmentation based method to large scale data set based on local Delaunay triangulation and local graph-cut distribution. It improves the tiling by optimising tile borders for them to cross regions of low complexity (i.e. flat areas). However, it does not guarantee watertightness along borders which is an important property for many applications. Moreover, [START_REF] Han | Distributed surface reconstruction from point cloud for city-scale scenes[END_REF] only deals with scenes in which a dimension is highly dominated by the two others (e.g. a city scene enables to "partition the scene on the ground plane" as said in their paper). A particular and interesting way of parallelising algorithms is to make them suitable for a Graphics Processing Unit (GPU) implementation to leverage this specifically-designed hardware which contains a massive set of cores. They are yet often limited to a small portion of algorithms namely Deep-Learning-oriented ones, several of them being introduced in Section 1.1. It is therefore a totally different kind of parallel computation which does not necessarily enable to process arbitrarily large point clouds since they still have to be loaded into the memory of one single machine. As an example, [START_REF] Zhou | Highly Parallel Surface Reconstruction[END_REF] and more recently [START_REF] Williams | Deep Geometric Prior for Surface Reconstruction[END_REF] both offer a GPU implementation of their algorithm.

Method WT OOC SG PL ALS LG AF OS Indicator Function (1 in / 0 out) ∀ x ∈ R 3 Fourier [14] ✓ voxels n i Poisson [11] ✓ octree n i ✓ Screened Poisson [13] ✓ octree n i ✓ Streaming Poisson [46] ✓ ✓ ✓ octree n i Parallel Poisson [55] ✓ ✓ ✓ ✓ octree n i Zhou et al. [56] ✓ ✓(GPU) octree n i IM-NET [16] ✓ voxels ✓ ✓ Occupancy Networks [18] ✓ ✓ octree ✓ ✓ Wavelets [57] ✓ ✓ ✓ octree n i ✓ Volumetric segmentation (1 in / 0 out)
The methods and the properties they offer are summarised in Table 1.1.

Evaluation of surface reconstruction

In order to assess the quality of a reconstruction, there is need for a ground truth, an input point cloud and a means of calculating the difference between a given output surface and the so-called ground truth. Let us present the various possibilities that have so far been considered for these three aspects.

Ground Truth

Ground truth could potentially take any surface form, i.e. implicit field, triangle mesh, volumetric segmentation, point set, deformed model, skeleton curve, primitives. However, only two have so far been considered: triangle mesh [START_REF] Ter Haar | A comparison of systems and tools for 3D scanning, 3D Digital Imaging and Modeling: Applications of Heritage, Industry[END_REF][START_REF] Kazhdan | Reconstruction of solid models from oriented point sets[END_REF] and implicit field [START_REF] Berger | A benchmark for surface reconstruction[END_REF].

Input Point Cloud

Producing point samples from a surface can be carried out in several ways:

• Real scanning: Based on a physical object (or scene), laser-based scanning generates a point cloud directly. Such technologies include Time-of-Flight [2] and Structured-light [3] devices. In addition, terrestrial or airborne LiDAR [4] offer the possibility to deal with large areas.

• Image-based: Multi-view stereo [START_REF] Furukawa | Multi-view stereo: A tutorial[END_REF] and Structure from motion [START_REF] Ozyesil | A survey of structure from motion[END_REF] allow creating a 3D model from images, which can be the starting point for surface reconstruction.

• Model sampling: Based on a continuous digital input model, synthetic sampling has the advantage of making it possible to fully control the data. In particular, one can generate more realistic data by adding noise, outliers, misalignment and occlusions, and by setting the density. In this field, several procedures have been considered: random or uniform sampling [START_REF] Kazhdan | Reconstruction of solid models from oriented point sets[END_REF][START_REF] Manson | Streaming Surface Reconstruction Using Wavelets[END_REF][START_REF] Süßmuth | Surface reconstruction based on hierarchical floating radial basis functions[END_REF], synthetic raytracing [START_REF] Hoppe | Surface reconstruction from unorganized points[END_REF][START_REF] Berger | A benchmark for surface reconstruction[END_REF] or z-buffering [START_REF] Ter Haar | A comparison of systems and tools for 3D scanning, 3D Digital Imaging and Modeling: Applications of Heritage, Industry[END_REF]. Of particular interest is the recent work presented in [START_REF] Winiwarter | Virtual laser scanning with HELIOS++: A novel take on ray tracing-based simulation of topographic full-waveform 3D laser scanning[END_REF] and [START_REF] Manivasagam | LiDARsim: Realistic LiDAR Simulation by Leveraging the Real World[END_REF]. Helios++ [START_REF] Winiwarter | Virtual laser scanning with HELIOS++: A novel take on ray tracing-based simulation of topographic full-waveform 3D laser scanning[END_REF] is an open-source tool for the simulation of airborne, UAV-based, and terrestrial static and mobile laser scanning implemented in C ++ . [START_REF] Manivasagam | LiDARsim: Realistic LiDAR Simulation by Leveraging the Real World[END_REF] developed

LiDARsim: a virtual terrestrial LiDAR platform generating realistic point clouds based on a high-quality mesh, free of moving objects.

Comparison

With regards to comparing an output reconstruction, six main possibilities have been explored:

• Visually: Most of the time, surface reconstruction aims at producing a digital representation as visually similar as possible to a real object. Hence, Poisson [START_REF] Kazhdan | Poisson surface reconstruction[END_REF], MPU [START_REF] Ohtake | Multi-level partition of unity implicits[END_REF] and SSD [START_REF] Calakli | SSD: Smooth signed distance surface reconstruction[END_REF] have simply compared models on a visual basis. This obviously raises the issues of being sensitive to the observer's perception, conflict of interest and the lack of quantitative information.

The following five methods come with the advantage of providing a quantitative quality assessment that is independent of any human bias.

• Point-to-mesh distance computation: When the only ground truth that is available comes in the form of a point cloud, it is relatively straightforward to compute the distance from each of those points to the reconstructed model. [START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF] have evaluated their method by randomly partitioning their point cloud into two equalsized subsets: one of points serving as input for the reconstruction algorithms and one of validation points from which distances to the output meshes are computed.

• Point-to-point distance computation: In the case where the ground truth comes in the form of a point cloud that has been acquired by scanning a real surface using a sensor, or if points have been sampled from a ground-truth surface, it is possible to sample the reconstructed surface and measure the distance between the two resulting point sets. In the computer vision community, the Chamfer distance [START_REF] Lu | Oriented and Directional Chamfer Distance Losses for 3D Object Reconstruction From a Single Image[END_REF] and the Hausdorff distance [START_REF] Huttenlocher | Comparing images using the Hausdorff distance[END_REF] are common ways of measuring the distance between two point sets. Given two point sets P 1 and P 2 , and defining the distance from a point p to a point set P as:

∀p ∈ R 3 , d (p, P) = min q∈P ||p -q|| 2 (1.14)
the Hausdorff distance is defined [START_REF] Huttenlocher | Comparing images using the Hausdorff distance[END_REF] as:

d Hausdorf f (P 1 , P 2 ) = max max p∈P 1 d (p, P 2 ) , max q∈P 2 d (q, P 1 ) (1.15)
the Chamfer distance is defined [START_REF] Lu | Oriented and Directional Chamfer Distance Losses for 3D Object Reconstruction From a Single Image[END_REF] as:

d Chamf er (P 1 , P 2 ) = 1 |P 1 | p∈P 1 d (p, P 2 ) 2 + 1 |P 2 | q∈P 2 d (q, P 1 ) 2 (1.16)
We believe that these metrics are not suited to the specific case of open scene reconstruction because they can be impacted by occlusions. Complex open scenes such as urban environments typically contain many vertical surfaces (building facades, for example) and overhangs (which can be found on bridges). Scanning this type of environment therefore often leads to a tremendous amount of occlusions.

Points associated to correct pieces of surface might be unjustifiably penalised by these metrics if lying in occluded areas, for which no ground truth is available.

• Mesh-to-mesh distance computation: [START_REF] Ter Haar | A comparison of systems and tools for 3D scanning, 3D Digital Imaging and Modeling: Applications of Heritage, Industry[END_REF] and [START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF] use the Metro tool [1] to compute the distance between two meshes. It works as follows: given two meshes (a sampled one M s and a target one M t ), Metro samples M s and measures the shortest distance from each sample to M t . Metro then computes the mean distance, the maximum distance and the Root Mean Square (RMS) distance over all samples.

It is also possible to inverse the roles of M s and M t and perform the same calculations. Instead of evaluating how close the reconstruction is to the ground truth, this will assess how complete the reconstruction is. The Hausdorff distance is defined as the greater value between the maximum distance from M s to M t and the maximum distance from M t to M s . Again, in the case of open scenes, the input point cloud that is fed to the reconstruction method is likely to contain many occlusions. Therefore, care must be taken when sampling the two meshes. Uniformly sampling the meshes is pertinent in the case of a closed object that can be scanned from different viewpoints. However, when dealing with open scenes, sampling the meshes in unobserved regions will lead the metrics to reflect the capacity of algorithms to fill holes rather than their reconstruction accuracy near the input data. This will be illustrated in Chapter 3.

• Mesh-to-implicit distance computation: [START_REF] Berger | A benchmark for surface reconstruction[END_REF] chose to use an implicit field that we will call Ω as the ground truth, and consequently they adapted the Metro methodology in order to compute the distance from a nearly uniform sampling of Ω to the evaluated mesh and vice-versa. The evaluation process answers the question: how well does the reconstructed mesh fit to the implicit surface computed by the MPU [START_REF] Ohtake | Multi-level partition of unity implicits[END_REF] In order to compute these, defining point correspondences between the two surfaces are needed. Let us denote by M the implicit surface and by M the output triangle mesh. As defined in [START_REF] Hildebrandt | On the convergence of metric and geometric properties of polyhedral surfaces[END_REF], the mapping Φ : M -→ M attributes to one point p ∈ M the intersection of the normal line through p and the mesh M . The inverse mapping Φ -1 : M -→ M attributes to Φ(p) its closest neighbour on the implicit surface M . This definition, associated with a sampling P M of M produces a set of nearest neighbour correspondences that we call C M (equation 1.17).

C M = {(x, p)|p ∈ P M , x = Φ(p)} (1.17)
By defining the corresponding operator Ψ : M -→ M and a sampling P M of the reconstructed mesh M , we get C M (equation 1.18).

C M = {(p, x)|x ∈ P M , p = Ψ(x)} (1.18) Denoting |S| = |C M | + |C M |
and with γ(p, x) the angle between the normals N M (p) and N M (x), error measures are the following:

H(M, M ) = max max (x,p)∈C M |x -p|, max (p,x)∈C M |p -x| (1.19) µ(M, M ) = 1 |S|    (x,p)∈C M |x -p| + (p,x)∈C M |p -x|   
(1.20)

H N (M, M ) = max max (x,p)∈C M γ(p, x), max (p,x)∈C M γ(p, x) (1.21) µ(M, M ) = 1 |S|    (x,p)∈C M γ(p, x) + (p,x)∈C M γ(p, x)    (1.22)
• Neural Feature Similarity: Recently, [START_REF] Huang | Surface Reconstruction from Point Clouds: A Survey and a Benchmark[END_REF] proposed a benchmark targeting deeplearning-based surface reconstruction algorithms. They compute the similarity between two surfaces in the deep feature space.

Conclusion

To sum up, whilst surface reconstruction has been extensively studied, few approaches [START_REF] Bolitho | Parallel Poisson Surface Reconstruction[END_REF][START_REF] Han | Distributed surface reconstruction from point cloud for city-scale scenes[END_REF] tackle watertight reconstruction of large open scenes. In addition, while metrics to assess surface reconstruction do exist, they are not adapted to deal with the many occlusions that open scenes typically contain. This calls for more contributions to scalingup and evaluating surface reconstruction from point clouds of open scenes. This PhD Chapter 2

Large scale surface reconstruction

In this chapter, we tackle the task of scaling up surface reconstruction. As highlighted in the Introduction, acquiring massive point clouds has been made possible with the improvement in sensor capabilities but processing them remains a challenge. In Chapter 1, Section 1.2, we surveyed methods contributing to this specific task but very few guarantee that the resulting mesh is watertight and matches the reference one. The major contribution of the work presented in this chapter is the first Delaunay-based watertight surface reconstruction algorithm that can handle arbitrarily large point clouds. Poisson surface reconstruction [START_REF] Bolitho | Parallel Poisson Surface Reconstruction[END_REF] is the only other known approach that scales while ensuring watertightness but it is not suited to urban scenes as it fails to reconstruct sharp edges.

Recently, [START_REF] Han | Distributed surface reconstruction from point cloud for city-scale scenes[END_REF] proposed an extension of a distributed workflow specific to urban scene reconstruction but it does not guarantee watertightness. Moreover, as their objective is to process urban environment only, the tiling strategy presented in their paper is purely 2-dimensional which restricts their algorithm to scenes in which a dimension is highly dominated by the two others. Conversely, the following end-to-end distributed out-ofcore algorithm guarantees watertightness and it is suited to any kind of large-scale point cloud with arbitrarily complex 3-dimensional geometry (overhangs...).

We proceed by first tiling the point cloud with a distributed octree approach which consists in choosing the maximum depth of an octree decomposition and then merging cells for which the total number of points that they contain does not exceed a given value (Figure 2.1a). Then, we compute the exact 3D Delaunay triangulation to discretise space using the algorithm introduced in [START_REF] Caraffa | Tile & Merge: Distributed Delaunay Triangulations for Cloud Computing[END_REF]. Based on a distributed triangulation structure with shared cells between tiles (blue triangles on Figure 2.1b), associated with the adjacency graph of the tiles (red lines on Figure 2.1b), we propose an extension of a distributed graph-cut algorithm [START_REF] Strandmark | Parallel and distributed graph cuts by dual decomposition[END_REF] (initially introduced for images) to carry out the in/out space labelisation that minimises our energy function (equation 2.1). The watertight surface can be extracted on the whole scene (blue line on Figure 2.1c) without any post-processing thanks to the consistent topology ensured by the global Delaunay triangulation.

Surface reconstruction model

In this section, we present the reference (undistributed) surface reconstruction algorithm [START_REF] Caraffa | 3D Watertight Mesh Generation with Uncertainties from Ubiquitous Data[END_REF] of which we present a distributed version in Section 2.2.

Local PCA

Several steps of the pipeline rely on a local PCA (Principal Components Analysis) of the inertia matrix of a carefully chosen neighbourhood of each point p ∈ P. The resulting normalised eigenvectors define a local frame f = (p, ⃗ v 1 , ⃗ v 2 , ⃗ v 3 ) and the corresponding eigenvalues (e 1 , e 2 , e 3 ) indicate how much the point cloud locally spreads in each direction. The eigenvector associated to the smallest eigenvalue provides a local (unoriented) surface normal direction. The position of the optical centre can then be used to consistently orient the normals from occupied to empty space.

Sub-sampling

As planar areas can be represented with large triangles, the input point cloud is subsampled to reduce its size, but without impairing the geometric precision. To this end, we use an adaptive algorithm based on PCA. Points p from the original point cloud P are randomly sampled and added to the set of sub-sampled points P S only if ∀ j ∈ {1, 2, 3}, p j = |(pq).⃗ v j | > σe j where p j are the coordinates of p in the local frame of its closest point q ∈ P S . The parameter σ controls the density of the resulting point cloud: as it increases, the distance threshold becomes larger so less points are added to P S which density decreases.

Delaunay Triangulation (DT)

Space is discretised with the DT [START_REF] Boissonnat | Algorithmic geometry[END_REF] of P S . In 3D, cells of this DT are tetrahedra which circumspheres do not contain any point of P S other than the 4 points lying on them. Additionally, their union exactly covers the convex hull of P S .

Mass computation

Following [START_REF] Caraffa | 3D Watertight Mesh Generation with Uncertainties from Ubiquitous Data[END_REF], the empty/occupied/unknown state of space at any point y ∈ R 3 is defined for each input line of sight associated with point p using Dempster Shafer Theory (DST) masses m p (y) = (e, o, u) ∈ 0, 1 3 such that e + o + u = 1. For example, for a given position y ∈ R 3 :

• e = 1 means certainty that space is empty

• o = 1 means certainty that space is occupied

• u = 1

means total uncertainty

These masses m p (y) are then combined into an overall mass m(y) = (e, o, u) and integrated over each tetrahedron t of the DT with a Monte Carlo sampling, yielding a single occupancy value m t = o/(1 -u) that gives the probability of t to be occupied.

Volumetric segmentation via graph-cut optimisation

The reconstruction method relies on assigning to each tetrahedron t of the DT the label 0 if t is empty or 1 if t is occupied. For simplicity we call T both the DT of P S and the set of tetrahedra of this DT. This labelling problem is formulated as a Boolean energy minimisation (2.1). The energy function is composed of a data term (2.2) encouraging labels to be in accordance with overall masses and a prior term (2.3) smoothing the surface, balanced by a parameter α:

x = arg min x∈{0,1} |T | E data (x) + α E prior (x) (2.1)
where x = (x t ) t∈T ∈ {0, 1} |T | is the set of tetrahedra labels and |T | the number of tetrahedra.

E data (x) = t∈T V t .|x t -m t | (2.
2)

E prior (x) = (t,t ′ )∈T 2 A t∩t ′ .|x t -x t ′ | (2.3)
where V t is the volume of tetrahedron t and A t∩t ′ is the area of the interface between tetrahedra t and t ′ .

This energy (2.1) can be optimised via a graph-cut algorithm. To do so, let us denote by G = (V, E) the dual graph of which the vertices V = {s, t} ∪ T are the source s, the sink t and the tetrahedra of T . The edges E are the triangles making the interface between all pairs of adjacent tetrahedra but also those connecting every tetrahedron to both s and t. Each edge is weighted so that the cost of an s-t cut in the dual graph has the same value as the energy function (2.1) after attributing the label values x. To minimise the energy, we therefore need to minimise the sum of cut edges weights:

(i,j)∈E c i,j x i,j (2.4)
where c i,j is the weight of edge (i, j) and x i,j is a Boolean indicating if the edge (i, j) is cut (in which case x i,j = 1) or not (x i,j = 0). Moreover, x i (the label of tetrahedron i) indicates whether vertex i is linked to s (x i = 0) or to t (x i = 1) after a given cut. We then have:

• x s = 0 • x t = 1 • ∀(i, j) ∈ E, x i,j = |x i -x j |

Surface extraction

The final surface is defined as the set of triangles separating occupied and empty tetrahedra.

Distributed surface reconstruction

In this section, we explain the distribution strategy. Algorithm 1 sums up in detail the different steps and Figure 2.2 presents the complete distributed workflow.

x 1 → 2 ,1 x 2→ 1, 1 x 2 → 3 ,1 x 3→ 2, 1 x 1 → 2 ,2 x 2→ 1, 2 x 2 → 3 ,2 x 3→ 2, 2 P P 1 P I 1 P S 1 T 1 T m 1 T x,0 1 x 0 1 T x,1 1 x 1 1 T x,t 1 x t 1 M 1 P 2 P I 2 P S 2 T 2 T m 2 T x,0 2 x 0 2 T x,1 2 x 1 2 T x,t 2 x t 2 M 2 P 3 P I 3 P S 3 T 3 T m 3 T x,0 3 x 0 3 T x,1 3 x 1 3 T x,t

Algorithm distribution

Tiling: Given a point set P, we note the partition P K = (P k ) k∈K its decomposition into |K| disjoint subsets P k , where K denotes a discrete set of tile indices (Figure 2.1a, Algorithm 1 line 2). We can define the primary tile of a point p ∈ P as the unique tile

P k such that p ∈ P k .
Local PCA: The PCA is distributed by computing it separately on each tile P k . The tile along with the PCA Information is denoted P I k .

Sub-sampling: Each tile P k is sub-sampled separately yielding Sub-sampled tiles P S k from which the surface will be extracted (Algorithm 1 line 6).

Distributed Delaunay triangulation:

We use the algorithm proposed in [START_REF] Caraffa | Tile & Merge: Distributed Delaunay Triangulations for Cloud Computing[END_REF] to compute a distributed DT of the union of the sub-sampled point clouds P S K = k∈K P S k (Algorithm 1 line 8). The distributed DT of a tiled point set P = k∈K P k is defined in [START_REF] Caraffa | Tile & Merge: Distributed Delaunay Triangulations for Cloud Computing[END_REF] as a set of DTs: For a mixed cell t, we call K t the set of indices of tiles containing t. We call a mixed cell main in T k if min(K t ) = k. The algorithm of [START_REF] Caraffa | Tile & Merge: Distributed Delaunay Triangulations for Cloud Computing[END_REF] guarantees that the union of all local and main mixed cells is exactly the DT of P and that these cells are disjoint. We define the adjacency graph of the triangulation G A = (T A , E A ) where T A are the tile-triangulations and E A the edges between them. Two triangulations are connected by an edge if and only if these triangulations share mixed points (see Figure 2.1b).

T K = (T k ) k∈K such that P k is a
T m k ← mass computation(T A k ∪ P I k ) end // Volumetric segmentation t ← 0 for k ∈ K do in parallel T x k,t , x k→l̸ =k,t ← graph cut(T m k ) end t ← 1 while t < max iterations do for k ∈ K do in parallel λ k t ← update lagrange ( l x l→k,t-1 ) T x k,t , x k→l̸ =k,t ← graph cut(T x k,t-1 , λ k t ) end t ← t + 1 end // Surface extraction for k ∈ K do in parallel T ′ k ← aggregate neigbors(T k,t , E A ) M k = extract surf ace(T ′ k ) end return M
Mass computation: Given that the effect of an input point p on a tetrahedron t decreases with the distance, the mass function is only computed locally on each tiletriangulation T k using P I k (Algorithm 1 line 10).

Volumetric segmentation:

The energy function 2.1 is optimised in parallel with independent graph-cuts on each tile (Algorithm 1 line 12) where the labels of mixed cells are encouraged to agree with an iterative scheme as explained in Section 2.2.2.

Surface extraction:

To extract the full surface, each tile is loaded with its respective neighbours regarding the adjacency graph's edges E (Algorithm 1 line 24). We denote by aggregate neigbors(T A k , E A ) the function that aggregates the triangulation and its respective neighbours. Because the iterative scheme does not provide the guarantee to converge on the same label on mixed cells (for which copies exist in multiple tile-triangulations), the label from the main mixed cell is systematically chosen for the surface extraction. Note that the surface is watertight even if the labels do not agree since the optimisation ends up with a unique set of labels which inevitably leads to a watertight mesh.

Distributed graph-cut

The main challenges to distribute a DT-based surface reconstruction method are distributing the DT itself, for which we use the implementation of [START_REF] Caraffa | Tile & Merge: Distributed Delaunay Triangulations for Cloud Computing[END_REF], and distributing the graph-cut optimisation. Graph-cut distribution approaches are highly dependent on the structure of the graph. While many approaches exist to distribute graph cuts as surveyed very recently in [START_REF] Jensen | Review of Serial and Parallel Min-Cut/Max-Flow Algorithms for Computer Vision[END_REF], we chose the dual decomposition paradigm because it offers the best distribution potential. Very little information has to be shared between the different tiles. While dual decomposition has been used for a graph representing the regular 2-dimensional structure of images [START_REF] Strandmark | Parallel and distributed graph cuts by dual decomposition[END_REF], we did not find any work tackling the specific structure of the adjacency graph G of the cells of a DT (see Section 2.1) and of its tilinginduced split into |K| sub-graphs G k = (V k , E k ). Thus, one of the main contributions of this work is the generalisation of the method proposed in [START_REF] Strandmark | Parallel and distributed graph cuts by dual decomposition[END_REF] to our tiled graph G which allows to distribute the graph-cut optimisation across multiple tiles. We use the structure of the tiled DT to split G into |K| subgraphs G k , one for each tiled-triangulation T k . In practice, each G k is defined as the adjacency graph of the set of local and mixed cells of T k (Figure 2.4). The main issue of the distributed problem is to ensure that the mixed cells have the same label relatively to the tiles that share them. For any k ∈ |K|, we denote by x k i,j and x k i the values of x i,j in sub-graph k and x i in sub-graph k respectively. We also define c k i,j = c i,j

|{k:(i,j)∈E k }| as the weight of edge (i, j) normalised by the number of tiles it appears in. Using this decomposition, the optimised function becomes:

f (x) = (i,j)∈E c i,j x i,j = k∈K (i,j)∈E k c k i,j x k i,j (2.5)
under the condition that labels x k i should have the same value amongst all subgraphs they appear into:

∀ k, l ∈ K 2 , ∀ i ∈ V k ∩ V l , x k i -x l i = 0 (2.6)
To split this problem in sub-problems using our tiled graph structure, we solve the graphcut on each sub-graph G k with variables x k i,j . To guarantee the condition expressed in equation 2.6, we add a penalty term to the energy using Lagrangian duality:

L(x, λ) = k∈K (i,j)∈E k c k i,j x k i,j + k∈K l>k i∈V k ∩V l λ k,l i (x k i -x l i ) (2.7) with λ = {λ k,l i : (k, l) ∈ K 2 , i ∈ V k ∩ V l }.
By defining, for any tile index k ∈ K:

L k (x k , λ) = (i,j)∈E k c k i,j x k i,j + l̸ =k i∈V k ∩V l    1 if l > k -1 otherwise    λ k,l i x k i (2.8)
We then have:

L(x, λ) = k∈K L k (x k , λ) (2.9)
Figure 2.3 allows to visualise the dual graph weights update in the case of a single split resulting in two tiles. We define the Lagrange dual function:

g(λ) = min x L(x, λ) (2.10)
Since the Lagrange dual function g(λ) (equation 2.10) is concave [START_REF] Boyd | Convex optimization[END_REF], we can find an optimal distributed solution to our problem (2.5+2.6) by maximising g(λ) via an ascent method. We denote by t the iteration index and by τ = {τ k,l i : (k, l) ∈ K 2 , i ∈ V k ∩V l } the vector of step amplitudes by which we increment λ. We initialise the Lagrange multiplier vector λ(t = 0) = 0 and the vector of steps τ (t = 0) by τ k,l i = τ 0 . As demonstrated in [START_REF] Strandmark | Parallel and distributed graph cuts by dual decomposition[END_REF], (x k -x l ) k,l∈K 2 is a supergradient to g at λ. Consequently, we can solve the distributed graph cut by iterating:

• Solve the graph-cut problems L k (x k , λ) for x k

• Update λ and τ : for k ∈ K and l > k :

λ k,l i,t+1 = λ k,l i,t + τ k,l i,t (x k i,t -x l i,t ) (2.11) τ k,l i,t+1 = τ k,l t 2 if x k i,t -x l i,t ̸ = x k i,t-1 -x l i,t-1
(2.12)

To perform this update at iteration step t, each tile G l sends to its neighbour tiles G k the labels of their shared mixed cells (according to G l ) in a vector x l→k,t .

Implementation details

The workflow is implemented with the Apache Spark framework [START_REF] Zaharia | Spark: Cluster Computing with Working Sets[END_REF]. The input of the algorithm is a point cloud separated in several files stored on the Hadoop distributed filesystem (HDFS) [START_REF] Shvachko | The Hadoop Distributed File System[END_REF]. To ensure a reasonable memory footprint in the tiling, we split large files to ensure that each input file contains less than 1 million points. Each such file is serialised into an RDD [START_REF] Zaharia | Resilient Distributed Datasets: A Fault-Tolerant Abstraction for In-Memory Cluster Computing[END_REF] with a Base64 encoding in a String. The key/value formalism is used: each element of the RDD is represented by a key and a value which is a list of sets (list of point sets, list of tile-triangulations, etc.). A transformation on an RDD (Local graph-cut, mass computation, etc.) is performed with a C ++ call in parallel on each RDD chunk by using the pipe operator. The union operator (∪ in Algorithm 1 and ⋄ in Figure 2.2) is a union followed by a ReduceByKey in Spark.

Experimental results

While memory and computational efficiency are the main goals when trying to distribute an algorithm, the most important aspect of evaluation is to verify that the distributed algorithm reproduces as closely as possible the output of the reference algorithm. For that reason, we will first evaluate the distributed algorithm in comparison with the reference, and then focus on the scalability. Moreover, as PCA and sub-sampling are expected to suffer minor border effects, and the distributed Delaunay triangulation is proved to produce the exact Delaunay Triangulation [START_REF] Caraffa | Tile & Merge: Distributed Delaunay Triangulations for Cloud Computing[END_REF], we mainly focus our evaluation on our major contribution: the distributed graph cut. In all the following experiments, unless otherwise stated, we used the values: α = 0.005, τ 0 = 5, and 30 iterations.

Distributed graph-cut convergence

Graph cut implementations provably find the global minimum cut of the graph [START_REF] Boros | Pseudo-Boolean optimization[END_REF] which equivalently achieves the global minimum of the associated energy (eq 2.1). We will therefore assess our distributed graph cut implementation by comparing the energy obtained with our distributed implementation with this global minimum energy. In order to evaluate our capacity to handle a large proportion of mixed cells (which makes the distributed problem harder), the algorithm is run with different octree depths for a fixed number of points. A larger depth induces more and smaller tiles, resulting in a higher proportion of mixed cells. As our iterative scheme is mainly influenced by the value of τ 0 , we also evaluate its impact on convergence. Figure 2.5 shows the evolution along iterations t of the ratio between the energy (eq 2.1) of the labelling computed with the distributed graph-cut E d and the labelling computed with the reference graph-cut E r . It is important to note that even if two surfaces with the same energy are not necessarily the same, they are considered equally good. For all values of τ 0 and octree depths, the energy ratio increases during the first iterations until it reaches a maximum around iteration 5, then decreases until converging significantly under the initialisation value,which shows the robustness of our distributed algorithm to these two factors. Figure 2.7 shows the resulting mesh along iterations on a LiDAR data set. Shared cells are visually the most impacted after the first graph-cut. In some areas, tiles are not or badly connected, only triangles between local cells at the centre of a tile are well reconstructed. After 3 iterations, the mesh is better regularised but the quality remains bad at the boundary. After 15 iterations, the main errors are removed and the tiles are globally well connected. At iteration 30, small shared areas like the car hood in the example are well reconstructed.

Speedup

The efficiency of the proposed approach is evaluated on a Spark cluster with 28 cores and 100GB of RAM. This configuration is close to the m51 setup of the Amazon EMR service dedicated to general usage that provides 4GB per core. We evaluate the strong scaling speedup which quantifies how the algorithm scales with an increasing number of cores for a fixed number of points. We compute the surface of a 10 million points data set with different cluster configurations: a varying number of executors (7,4,2 and 1) and, for each executor, a varying number of cores (4,3,2,1) for a fixed 12GB of RAM per executor in every case. The total number of cores is the number of executors times the number of cores per executor. The result is shown on Figure 2.6 where the strong scaling factor is defined as t 1 core /(n.t n cores ), where t x cores is the time it takes to compute the solution to the problem using x cores. The strong scaling factor is expressed as a function of the number of cores. As our application is fully distributed, we compare it to the perfect scale factor of 1. Results show that every configuration has a scale factor > 0.6. In other words, the execution time with a single thread is at least divided by 0.6 times the number of cores with the distributed implementation. n.tn cores (right) as a function of the number of cores with 4, 3, 2 and 1 core(s)/executor (12Go RAM/executor) on a 10 million points data set.

Large data set result

We tested our algorithm on a very large data set that combines more than 350 million points acquired with both aerial and terrestrial LiDAR (Figure 2.8), on which our available hardware could not run the reference implementation by lack of RAM. The point cloud was split into 9520 tiles and it took 17 hours to process using 28 cores. The proposed approach produces a watertight mesh of the whole area with highly detailed and complex 3D structures, composed of around 80M vertices. We empirically set the sub-sampling parameter σ = 0.3 to have a good trade-off between accuracy and computing time (a higher σ would reduce computing time at the cost of accuracy).

Conclusion

In this chapter, we presented a fully distributed and out-of-core surface reconstruction method. The extension of the distributed graph-cut method of [START_REF] Strandmark | Parallel and distributed graph cuts by dual decomposition[END_REF] to the structure of the tiled 3D Delaunay Triangulation of [START_REF] Caraffa | Tile & Merge: Distributed Delaunay Triangulations for Cloud Computing[END_REF] leads to a fully distributed method that ensures both watertightness of the resulting mesh and preservation of sharp edges. The evaluation confirms that the proposed method produces meshes of similar quality to the reference (non-distributed) implementation. The fully distributed framework allows an efficient speedup while maintaining a reasonable memory footprint even for very large inputs. The implementation on the state-of-the-art Apache-Spark framework shows very good results on managing very large data sets.

Chapter 3

Evaluating surface reconstruction using synthetic data

This chapter addresses the evaluation of algorithms reconstructing a surface from a point cloud, using synthetic data. The objective is to set a rigorous protocol measuring the quality of the reconstruction. This work is of fundamental interest for two reasons. First, it enables the assessment of different algorithms tackling the same problem to determine the best approach. Secondly, as most algorithms can be tuned via a set of parameters, it makes it possible to quantitatively study the influence of the different settings and adjust them in order to achieve the best performance. This is a new approach, that introduces less bias than previous contributions. We use ground truth surfaces and we proceed by filtering out parts of the meshes that are not relevant to assess, at different thresholds. We compute distances between samplings of these filtered meshes and samplings from the filtered ground truth surfaces.

Evaluation protocol

Input data

In order to assess surface reconstruction, we first assume the availability of a ground truth model taking the form of a triangle mesh M GT . Secondly, a point cloud P representing a realistic scan of M GT must have been generated using a LiDAR simulator like the one we present in section 3.2. Finally, we will denote by M E a reconstructed triangle mesh produced with one of the methods to evaluate.

Intuition

As explained in Chapter 1, one way of assessing the quality of a given reconstruction is to uniformly sample the corresponding mesh M E and to compute the distance from each of these samples to the ground-truth mesh M GT , and then doing the same but inverting the Figure 3.1: In this situation, the real surface has been partially scanned, resulting in a serious occlusion in the input point cloud. If we uniformly sampled the entirety of both reconstructed meshes R 1 and R 2 and used the Metro tool [1] to compute the distances from them to the ground-truth, the metrics would be dominated by occlusion-based errors (highlighted in bold).

roles of M E and M GT . Scanning open scenes such as urban environments often results in severe occlusions because of the presence of many vertical surfaces and overhangs. Sampling points from the reconstructed surface M E in unobserved regions and taking into account the corresponding distances to M GT is very likely to dramatically affect the global performance. This is because the main errors made by the algorithm are very likely to correspond to those occluded areas. This approach would thus give an indication of the capacity of the algorithm to fill holes rather than its reconstruction accuracy near the input data. This is illustrated in Figure 3.1.

Best expected surface

In the general case, as we explained it in Section 3.1.2, only a limited portion of the ground truth surface M GT is covered by the acquisition. If we can expect an algorithm to fill some holes in the data in a reasonable manner, we cannot expect it to recover the shape of the scene far from this covered area. Thus, we need to define the "reconstructible" part M α GT of M GT , which can be non connected and have holes, but is still orientable and manifold. To this end, we propose to define an explicit maximum interpolation (hole filling, super resolution) and extrapolation (uncropping) distance α at which we evaluate the algorithm. We can then compute M α GT by removing all triangles of M GT for which none of its vertices lie closer than α to a point of the sampling P. See algorithm 2 for a description of this mesh processing step. M α GT is then the best surface an algorithm can be expected to reconstruct. Consequently, this is the surface we will compare every reconstructed mesh with. However, watertight reconstruction pipelines interpolate surface even where point samples are absent. While watertightness is desirable for the reconstruction of closed objects, it is less for open scenes that naturally have a border corresponding to the limit of the coverage of the scene by the acquisition. In order to be impartial with every algorithm prior, we will not evaluate those interpolated mesh parts further away than α from the input data. As a consequence, we also need to compute the part of M E that is relevant to assess i.e. the sub-mesh M α E of M E (applying algorithm 2) containing only triangles closer to P than a distance α. This is not limiting as in practice, it is an easy post-processing step that the user can choose to perform if he does not require a fully watertight mesh. For all these reasons, we believe this approach introduces less bias than previous work.

Sampling

Once M α GT and M α E have been computed, we need a way to measure the difference between them. We consider that computing error metrics based on distances between vertices of one mesh and triangles of the other is biased, as different algorithm can produce triangles of very different size. Thus we propose to perform a Poisson-disk sampling of the triangles, which guarantees an even distribution of samples so that we can consider that each sample represents the same amount of surface area. We choose a Poisson-disk radius R significantly smaller than α. We denote by P

(v 0 , v 1 , v 2 ) = vertices(t) for each vertex v t ∈ {v 0 , v 1 , v 2 } ∈ do if v t / ∈ M α then add v t to M α end add t to M α end end

Quality measure

Thanks to notations introduced in Section 3.1.4 and using the point-to-mesh distance d (p, M) defined in equation 3.1, bi-directional distances can be computed (equations 3.2, 3.3).

∀p ∈ R 3 , d (p, M) = min q∈M d(p, q) (3.1) Mean Precision: 1 |P R M α E | p∈P R M α E d (p, M α GT ) (3.2) Mean Recall: 1 |P R M α | p∈P R M α GT d (p, M α E ) (3.3)
Precision measures how close are points from the reconstructed mesh to the ground truth and recall indicates which distance separates points from the ground truth to the reconstruction. While these terms are named in analogy with the machine learning community metrics, it is important to note that they do not measure a ratio of relevant information but a distance and thus the lower, the better. Instead of choosing an arbitrary α, we propose to compute these distances for a range of α values and draw the charts corresponding to each of these. Note that the case where α -→ ∞ corresponds to computing the distances on samplings of the raw meshes M GT and M E , without applying algorithm 2. The resulting curves will both indicate the precision of the algorithm for small α values and the quality of the interpolation/extrapolation/hole filling for larger α values.

Aerial LiDAR simulator

Similar to LiDARsim [START_REF] Manivasagam | LiDARsim: Realistic LiDAR Simulation by Leveraging the Real World[END_REF] for terrestrial scan, we developed our own airborne LiDAR simulation platform in order to generate realistic scans of a given environment.

Virtual environment

We used the open dataset from [START_REF] Et | eurometropole de Strasbourg[END_REF] which was financed by the European Union as part of a FEDER (Fonds Européen de Developpement Régional). It consists of a 3D mesh of a large area covering the metropole of Strasbourg. It was produced by photogrammetry using high resolution (between 4 and 7cm GSD) oblique imagery acquired by helicopter platform such that it presents details at a higher resolution than typical aerial LiDAR acquisitions. Figure 3.3 shows a 250m by 250m tile of this mesh.

Scanning process

In this section, we formalise our aerial LiDAR simulator. We assume the plane trajectory can be modelled as a sequence of straight lines and we implemented two scanning patterns. Among the many possibilities available on the market (zig-zag, elliptical, circular...), we chose to replicate the parallel line and the elliptical scanning pattern produced by a rotating mirror mechanism.

Plane trajectory

We use (O, -→ e x , -→ e y , -→ e z ) as the global coordinate frame, the one in which mesh vertices coordinates are expressed as detailed in Figure 3.4. We model the acquisition by a linear trajectory of the LiDAR optical center M moving from A (x A , y A , z A ) to B(x B , y B , z B ) at constant speed v 0 . We also use a local coordinate frame M,

- → i , - → j , - → k associated to M defined as: - → k = -→ AB -→ AB ; - → j = - → e z ∧ - → k - → e z ∧ - → k ; - → i = - → j ∧ - → k (3.4)
which ensures that -→ j is orthogonal to -→ AB and to the vertical direction -→ e z . Note that -→ k is undefined only when -→ AB and -→ e z are colinear which never happens in practice (a plane does not fly vertically).

-→ i completes the frame so that

- → i , - → j , - →
k is right-handed. M (t) moves in a straight line and thus can be defined as:

∀ t ∈ [0, t B ] , --→ OM (t) = -→ OA + v 0 t - → k (3.5)
Figure 3.4: Global frame and flight trajectory.

Parallel-line scanning pattern

We denote as -→ r the direction of the laser ray. -→ r is rotating around -→ k at constant angular speed ω = θ (see Figure 3.5). Thus:

- → r = cos(θ) - → i + sin(θ) - → j (3.6) ∀ t ∈ [0, t B ] , θ(t) = ω t + θ 0 (3.7)
Such a system is also defined by its field of view i.e. the angle range [θ min , θ max ] to which θ must belong for laser pulses being actually emitted, and by the pulse rate (frequency at which pulses are emitted) f p .

Elliptical scanning pattern

Denoting as -→ r the direction of the laser ray, and using (M, -→ u , -→ v , -→ w ) as the canonical spherical coordinate frame (see Figure 3.6), -→ r is rotating around -→ w at constant angular speed ω = φ with ϕ being the azimuthal angle of -→ r . The polar angle θ is constant.

The laser ray's direction -→ r is:

- → r = sin(θ) cos(ϕ) - → u + sin(θ) sin(ϕ) - → v + cos(θ) - → w - → r = sin(θ) cos(ϕ) -→ -j + sin(θ) sin(ϕ) - → k + cos(θ) -→ -i - → r = -cos(θ) - → i -sin(θ) cos(ϕ) - → j + sin(θ) sin(ϕ) - → k with : • ∀ t, θ(t) = θ 0 (e.g. π -π 9 ) • ∀ t, ϕ(t) = ωt + ϕ 0

Noise model

The quality of LiDAR data has recently been surveyed by [START_REF] Warcho\l | The concept of LiDAR data quality assessment in the context of BIM modeling[END_REF]. We are especially interested in the accuracy of the point coordinates that a real LiDAR can achieve. Most of studies in this area agree to state that it can be split up in an altimetric and a planimetric component. Values are obviously influenced by the modernity of the system but they also vary depending on the type of terrain that is considered (bare soil, low grass, forestry) and the flight parameters (altitude, speed). See [START_REF] Aguilar | Accuracy assessment of LiDAR-derived digital elevation models[END_REF], [START_REF] Aguilar | Modelling vertical error in LiDAR-derived digital elevation models[END_REF] for studies on altimetric error and [START_REF] Toth | Quality assessment of lidar data by using pavement markings[END_REF], [START_REF] Vosselman | Analysis of planimetric accuracy of airborne laser scanning surveys[END_REF] for planimetric error analysis. As suggested by these contributions, we

Implementation

We implemented the method described above in C ++ using CGAL library [START_REF]CGAL User and Reference Manual, 5th Edition[END_REF]. Rays are traced from the virtual aerial station in direction -→ r and their exact intersections with the virtual environment are computed. Gaussian noise following 3.8 is added to these perfect intersections. As oriented normals or sensor positions are required for some reconstruction algorithms, point clouds containing this additional information are generated after any simulation.

Evaluation

The algorithms we assessed have been run on a point cloud we generated using the LiDAR simulator described in section 3.2 and the virtual environment displayed in Figure 3.3. We selected the following ones:

• • SSD: Smooth Signed Distance [START_REF] Calakli | SSD: Smooth signed distance surface reconstruction[END_REF] • PSS: Point Set Structuring [START_REF] Lafarge | Surface reconstruction through point set structuring[END_REF] • RESR: Robust and Efficient Surface Reconstruction [START_REF] Labatut | Robust and efficient surface reconstruction from range data[END_REF] • WMWU: Watertight Mesh generation With Uncertainties [START_REF] Caraffa | 3D Watertight Mesh Generation with Uncertainties from Ubiquitous Data[END_REF] Resulting surfaces are shown on

Experimental parameters and methodology

To generate the point cloud on which we evaluated the various algorithms, we used the parallel-line scanning pattern of the aerial simulator presented in Section 3.2. Values of the experimental parameters we used are to be found in Table 3.1. As regards to the evaluation part, we set to R = 0.3 m the Poisson-disk sampling radius. This parameter is chosen by considering the trade-off between sampling density and computational time.

As every algorithm we assessed can be tuned, we first studied the performance of each of them for a small number of α's, changing some parameter values. We then selected the best version of each algorithm and carried out the full evaluation for which the results are presented in section 3.3.2.

Quantitative results

We plotted the mean precision on Figure 3.8 and mean recall on Figure 3.9 (following equations 3.2 and 3.3) as a function of α. RESR [START_REF] Labatut | Robust and efficient surface reconstruction from range data[END_REF] achieves the best performance both in terms of precision and recall and for any value of α. Additional information provided by sensor positions is certainly helping but we can assume they are making the best use of it, in comparison to WMWU [START_REF] Caraffa | 3D Watertight Mesh Generation with Uncertainties from Ubiquitous Data[END_REF].

Symbol Value Unit

Results of the evaluation of Poisson [START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF] confirm that Neumann (Poisson N) leads to a better reconstruction of open scenes as it enables the surface to extend out to the boundary of the domain. In contrast, Dirichlet (Poisson D) enforces hard watertightness. Consequently, an undesirable closure of the surface from the bottom leads to a poor precision, but only for high α's. We now see how important it is to differentiate between global and local evaluation.

WMWU [START_REF] Caraffa | 3D Watertight Mesh Generation with Uncertainties from Ubiquitous Data[END_REF] evaluation shows a comparable result as those provided in the original article. This method provides good results regarding the recall metrics while precision quickly worsens when α increases. The main explanation is the same as for Poisson with Dirichlet boundary condition: the surface is artificially closed to ensure watertightness, leading a great amount of surface to lie away from the ground truth. SSD [START_REF] Calakli | SSD: Smooth signed distance surface reconstruction[END_REF] curve logically follows the one of Poisson N, as both methods use smooth basis functions to solve their equation. PSS [START_REF] Lafarge | Surface reconstruction through point set structuring[END_REF] could have been expected to perform better since it is based on primitive detection and an urban scene is full of planar regions. However, it is quite sensitive to missing data and facades of the buildings are mostly out of reach for a parallel line scanning pattern LiDAR.

Conclusion

Surface mesh reconstruction from remote sensing point clouds is a challenging task that becomes more and more important as data acquisition starts from a purely vertical perspective. Our evaluation shows important differences between the various state of the art approaches. In this chapter, we evaluate the quality of the interpolation that surface reconstruction algorithms perform between the input points. We argue that this should be measured as a function of the distance at which we expect the algorithm to interpolate The resulting metric take the form of curves indicating two quality criteria as a function of the maximum interpolation distance. Our study shows that [START_REF] Labatut | Efficient Multi-View Reconstruction of Large-Scale Scenes using Interest Points, Delaunay Triangulation and Graph Cuts[END_REF], while quite old, remains a very good choice. The very popular Poisson method [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] is efficient and scales well but does not preserve sharp features that are very present in our urban evaluation scene.

Chapter 4

Evaluating surface reconstruction using real data

Introduction

The evaluation of surface reconstruction from synthetic data, as presented in Chapter 3, offers several advantages including the possibility to directly compute the difference between the reconstructed surface and the ideal one. However, high-quality synthetic meshes representing large scenes are expensive to obtain thus reducing the range of possibilities. This Chapter tackles the evaluation of surface reconstruction in the real-world case where we do not have access to such a synthetic ground truth. We call "real data" the data acquired in the physical world with real sensors. This includes LiDAR scans, images, RGB-D images, etc. As is usually done to address this issue, we assess the reconstruction of real scenes from real data only based on other real data of significantly higher quality. Even though this idea is quite typical, the main contribution of the work presented in this chapter lies in the way that we assess the difference between the reconstructed surface and the high-quality real data as inconsistencies, inspired by recent work on change detection [START_REF] Xiao | Street environment change detection from mobile laser scanning point clouds[END_REF]. This endeavour is hard because of several aspects.

Limits in the quality of the ground truth: the specificity of working with real data is the presence of noise in what we consider as the ground truth. In addition, real data is always sparse and incomplete, which means that we do not know the state of space (occupied by the object or empty) everywhere. This raises the question of how to assess pieces of reconstructed surface in unseen, unobserved regions.

Our contributions are twofold:

1. We propose a setting where the high-quality data used to compute metrics is significantly better than the low-quality data on which surface reconstruction is performed in three separate ways:

• Coverage: we use multiple data sources acquired from multiple points of view to ensure that the high-quality data has a significantly better coverage of the surface to reconstruct than the low-quality data.

• Density: we ensure that the density of the points in the high-quality data is significantly better than that of the low-quality data.

• Noise: we ensure that the noise level is lower in the high-quality data than in the low-quality data.

2. We propose metrics that penalise inconsistencies between the surface to be evaluated and the high-quality data: a piece of surface reconstructed within a volume unseen by the high-quality data will simply not be evaluated as we have no information on it. This does not mean that we do not evaluate the hole-filling capacity of the evaluated methods. As the high-quality data has more coverage, we evaluate hole filling exactly where we have the data to do so.

Assumptions and priors: algorithms make different assumptions about the type of shape that needs to be reconstructed and this leads to very different properties. Consequently, depending on the metrics' definitions, these assumptions can dramatically influence the assessment and sometimes in an unjustified manner. An example of such a situation is shown in Figure 1 where the left model would be attributed a bad mark because of the red piece of surface even though the rest of the model is correct. Does this red piece of surface need to be taken into account when assessing the model? We tackle this issue by defining metrics that only assess hole filling where relevant high-quality data is available.

Our proposed metrics are based on the visibility information contained in the highquality data. We assess the reconstructed surface only where the real one has been observed. For that purpose, we make an extensive use of sensor positions (positions from which points have been acquired). This information is easy to access and it provides us with a full ray along which we know that space is free, instead of just a single position where we know that the real surface lies. We used these newly-defined metrics to assess several open-source surface reconstruction algorithms and a licensed one on different types of scenes. In Section 4.2, we define our metrics and discuss their nature. In Section 4.3, we present the three datasets on which we tested our evaluation protocol and present the experimental setup we used to generate the high-quality data. Results are detailed and analysed in Section 4.4, before we present the conclusions of this work in Section 4.5.

Evaluation protocol 4.2.1 Intuition

Let us give some examples of situations where current metrics are not adequate to evaluate surface reconstruction, and what we suggest would be an improvement. This is going to help understand our metrics' definitions in Section 4.2.3.

First, as presented in Section 1.3, comparing a reconstructed mesh with a ground truth point cloud can be done by computing the distances from those points to the mesh model. While this seems like a good starting point to assess how well holes have been filled around those points, it is inadequate to evaluate the overall accuracy of the surface. Figure 4.1 shows an example of such a situation where a surface would be evaluated as almost perfect even though large portions are clearly incompatible with the ground truth if we take into account the positions from which points have been acquired.

Secondly, it is possible to measure accuracy solely with ground truth points by sampling the reconstructed surface and measuring the distance from these samples to the ground truth points. Nevertheless, large pieces of the reconstructed surface might be judged as being of poor quality (if lying far from the nearest ground truth point) despite being correct, just because of a low ground truth density. We want to assess both accuracy and completeness only in regions where ground truth information is available, and this is possible using sensor positions as shown in Thirdly, as surface reconstruction has often been evaluated visually, we wanted to find metrics that would imitate this human intuition-based assessment. We believe that visibility-based metrics are more appropriate for assessing how the reconstructed surface matches the real one everywhere where we can see and compare them. Figure 4.2 shows four situations for which the piece of reconstructed surface would be marked similarly by a point-to-mesh distance. The distance from the high-quality point to the nearest piece of reconstructed surface is indeed the same in all four situations. However, we are certain that they should correspond to three completely different outcomes and we want our metrics to be able to differentiate between them. In particular:

• in (a), the reconstructed surface lies slightly behind the high-quality point. We consider it as correct at the threshold defined by our tolerance zone.

• in (b), the reconstructed surface lies slightly in front of the high-quality point. As in (a), we consider it as correct even though we can measure a slight error.

• in (c), whilst the piece of real surface corresponding the high-quality point has been recovered just like in (a), the reconstructed surface hides the nearest intersection by crossing the laser ray, resulting in an obvious inaccuracy. What one would see by looking in the direction of the laser ray is not the right piece of reconstructed surface but rather something else far in front of it.

• situation (d) might look similar to (b) but there is actually no intersection between the reconstructed surface and the laser ray. Consequently, what one would see by looking in the direction of the laser ray is not the right piece of reconstructed surface but something else in the background. We thus consider that the surface has not been recovered at all here.

ETH3D benchmark [START_REF] Schops | A Multi-View Stereo Benchmark With High-Resolution Images and Multi-Camera Videos[END_REF] proposed metrics to evaluate how well a two-or multi-view stereo point cloud matches a LiDAR-based one, using sensor positions. While surface reconstruction is a different task to MVS (in particular in terms of the expected output properties), [START_REF] Schops | A Multi-View Stereo Benchmark With High-Resolution Images and Multi-Camera Videos[END_REF] still evaluates the matching between two 3-dimensional structures and we were inspired by their use of sensor positions. However, we chose different paradigms that are more adapted to surfaces, so our evaluation protocol is considerably different.

The ETH3D benchmark [START_REF] Schops | A Multi-View Stereo Benchmark With High-Resolution Images and Multi-Camera Videos[END_REF] defines completeness as the proportion of ground truth points for which the distance to its closest reconstructed point is below a given threshold. We could keep this definition and find the point from the reconstructed surface that minimises the distance to each point of the high-quality point cloud. Nonetheless, given that we know in what direction this point is supposed to be encountered thanks to the sensor position, we find it more relevant to compute the distance between the high-quality point and the nearest intersection between the corresponding laser ray and the mesh model. In other words, while the most natural adaptation of this point-to-point distance would be a point-to-mesh distance, we believe that for each ray that hit the the real surface there should be a piece of reconstructed surface close by and along the ray.

However, we also leverage the information given by each laser ray: the space between each sensor position and its associated high-quality point should be empty. We soften this property by defining a tolerance zone: a piece of reconstructed surface will be considered as correct if its distance along the ray from the high-quality point is smaller than a given threshold d max . Every piece of surface further than d max and situated in front of the high-quality point will affect the accuracy of the model. Contrary to ETH3D benchmark [START_REF] Schops | A Multi-View Stereo Benchmark With High-Resolution Images and Multi-Camera Videos[END_REF], we do not model the shape of a laser beam as a truncated cone. The first reason for this is that we do not need to, since the model we are trying to evaluate (a surface) is continuous instead of discrete (a point cloud). Hence, we are not at risk of missing any part of it. In addition, it makes it simpler to get a point as the intersection between a ray and a piece of surface. This way, every couple (ray -piece of surface) gives the same amount of information.

[87] uses voxels to prevent a "cheating" strategy from achieving both high accuracy and completeness despite raising other issues. For example, regions of low ground truth density contribute as much as high density ones while encapsulating less information. A cheating strategy for surface reconstruction would be to add several parallel layers of surface in regions of high confidence. We do not need to discretise space as in [START_REF] Schops | A Multi-View Stereo Benchmark With High-Resolution Images and Multi-Camera Videos[END_REF] since for each ray we propose to keep only the closest intersection as a potential correct one and penalise all the ones situated in front of it. Note that even layers situated behind the closest intersection might be obstacles to rays pointing at another object in the background.

If the nearest intersection is found behind the high-quality point at a greater distance than d max , or if no intersection is found at all, then we consider that this piece of surface has not been recovered. This therefore affects the completeness of the model.

Definitions and notations

In this chapter, we want to evaluate the quality of surface reconstructions from lowquality data P If , with only having access to high-quality data P GT of the same scene and without having access to the perfect ground truth surface that the algorithms are supposed to produce.

• P If : the low-quality point cloud that will be fed to the evaluated surface reconstruction methods to produce the output surface meshes to be evaluated.

• P GT : the high-quality (ground-truth) point cloud with better coverage, higher density and less noise than P If and for which we know the sensor positions, defining one ray per point.

• M E : the reconstructed mesh to evaluate, produced by an algorithm from P If .

• d max : the maximum distance at which we evaluate the reconstruction. It is a parameter which influences the different metrics as we use it to separate noise (distance We compute all intersections between the associated virtual ray (i.e. the extension of the real laser ray) and the reconstructed surface. In this case, it results in six intersections I 1,...,6 . The closest intersection to Gt happens to be I 2 so the "ray distance" metric for that particular ray is the distance (Gt, I 2 ). If (Gt, I 2 ) < d max , Gt is to be counted as a TP, otherwise it will not be taken into account in the evaluation since it is situated after Gt. Note that this piece of surface might still be evaluated thanks to another ray (as one emanating from O 2 , for example). Besides, we found one intersection I 1 on the way to the closest intersection I 2 which is counted as a FP.

< d max ) from outliers (distance > d max ).

Metrics definitions

Similar to [START_REF] Marchand | Evaluating Surface Mesh Reconstruction of Open Scenes[END_REF], we would like to assess both the precision of each part of the reconstructed surface (each part of the surface should lie near some part of the real surface) and the completeness of the model (there should be as few missing parts of the real surface as possible). As we do not have access to a digital model of the real surface, we cannot compute the actual precision and recall as in [START_REF] Marchand | Evaluating Surface Mesh Reconstruction of Open Scenes[END_REF]. Our knowledge of the ground truth is limited to the high-quality data P GT . However, we also know where the surface is not supposed to lie since we know the sensor position from which every point of P GT has been acquired, thus defining a ray of free space. Therefore, we define a precision metric that penalises inconsistencies between the reconstructed surface and the visibility information contained in P GT . We propose an equivalent of the recall metric: for each [OP ] ray from P GT , we compute the distance from the P to the closest intersection between the [OP ) half line and the reconstructed surface. Here is the formalisation of these metrics in more detail:

• "Ray distance": for each point/ray (p, r) ∈ P GT , we compute the distance from p to the closest intersection (which we will denote as c) between r and M E (among all potential intersections, we choose the one with the smallest distance with respect to p). If this distance is < d max then the piece of reconstructed surface is considered as being correct and we add this distance to an array of distances.

• "Precision": for each point/ray (p, r) ∈ P GT , if the "ray distance" (between p and c) is < d max then we count c as a True Positive (TP) (since there is a piece of surface and it is correct). Otherwise, if the intersected point c lies at a distance greater than d max and it is before the corresponding high-quality point p, we consider it as being false and we count it as a False Positive (FP) (since there is a piece of surface and it is false). We consider that we cannot say anything about the closest point c if it lies at a greater distance than d max and it is situated after the corresponding high-quality point p (neither can we for all intersected points lying after it), so we just ignore them. All intersected points lying before this closest intersection c are also counted as FP as they are inconsistent with the corresponding ray of free space (p, r). This is enough to define the precision ratio (equation 4.1).

• "Recall": it is defined as the ratio between the number of TP and the number of cast rays (equation 4.2. Every high-quality point p ∈ P GT is either mapped to its corresponding TP (in which case the piece of real surface has actually been recovered by the algorithm) or not (meaning a lack of exhaustiveness of the reconstruction and corresponding to a False Negative FN ). The number of rays thus equals the sum of the TP and the FN.

• "Cumulative distances": the cumulative histogram of the ray distances where the x-axis corresponds to the distance and on the y-axis we plot the number of points for which the ray distance is below the x-distance, divided by the total number of rays cast. It contains the information of both mean ray distance and recall. The right-most value will be at a given distance from the y = 1 line, which corresponds to the number of points for which the ray distance is above d max (if "re-multiplied" by the total number of rays).

Precision = T P T P +

F P (4.1) Recall = T P T P + F N = T P number of rays (4.2) F-score = 2 • precision • recall precision + recall (4.3)
As we define precision and recall as ratios, the harmonic mean (F-score, see equation 4.3) allows ranking the methods by taking into account both metrics.

Tuning / Training

The algorithms we evaluate in this chapter are either tunable for the most part or they learn parameters in order to reconstruct surfaces. For example, DGNN [START_REF] Sulzer | Scalable Surface Reconstruction with Delaunay-Graph Neural Networks[END_REF] needs a training dataset in order to learn the parameters of its model and Poisson [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] can be run at different resolutions by changing the depth of the octree that is used. In order to be as fair as possible, we used the same dataset to tune or train the algorithms. We therefore ran the non-learning-based methods with different values for their parameters and carried out an evaluation. A first interesting result is that for some methods, it can be hard to obtain a good performance regarding the precision and the recall metrics at the same time. For PSS, the higher the value of the trade-off parameter λ (therefore the more importance given to the prior term, see equation 1.2), the higher the precision (up to a certain value) but the lower the recall. Maximising the score of one metric (by varying λ) results in minimising the score of the other one (we explain why in Section 4.4). Thus, we selected the two λ values that maximise each of these metrics individually.

The tuning/training dataset that we decided to use is composed of three scenes from STRAS . The reason behind this choice is the availability of ground-truth meshes which are absolutely necessary for the training phase of DGNN. We then also used these three scenes to find the best parameters (in terms of ray distance, precision and recall) for the non-learning-based algorithms.

Input data

We aim to evaluate the algorithms in very different scenarios as the methods' priors might influence the quality of the reconstruction depending on the type of scene or the type of data involved. We therefore compute the metrics introduced above on three significantly varying datasets.

STRAS: Strasbourg dataset and LiDAR simulator

Elliptical aerial LiDAR simulator

In order to control the data itself, we started by using the same synthetic dataset [START_REF] Et | eurometropole de Strasbourg[END_REF] as in Chapter 3. Figure 4.4 shows a 250 metre by 250 metre tile of this mesh. In order to generate the ground-truth and the input point clouds P GT and P If , we used the aerial LiDAR simulator that we introduced in Chapter 3 [START_REF] Marchand | Evaluating Surface Mesh Reconstruction of Open Scenes[END_REF]. This is because such a large urban environment is typically scanned using an airborne LiDAR system. The problem with the parallel line pattern is that facades that are perpendicular to the direction of the plane are not reachable by the laser ray and thus are absent from the resulting point cloud. In order to overcome this issue, we used the elliptical scanning pattern. It is indeed better suited to urban environments, for the laser ray will be able to point at far more facades than with the parallel line pattern, resulting in a better coverage. The values of all the parameters we used can be found in Table 4.1.

Experimental setup

We aim to produce two point clouds ) such that -→ e z represents the ascending vertical direction. In our setup, a trajectory of the LiDAR system for each scene is a single straight pass of the plane along axis -→ e y for x = α x (x max -x min ) , α x ∈ [0, 1]. We generate P If thanks to one pass of the plane with α x = 0.5 and P GT thanks to three passes with α x ∈ {0.25, 0.5, 0.75}. We denote P αx as the point cloud resulting from the flight x = α x (x max -x min ). We then have: P If = P 0.5 and P GT = P 0.25 ∪ P 0.5 ∪ P 0.75 . This way, P If forms part of P GT and contains a lot more occlusions in particular on facades parallel to the direction of the plane.

P

ENSG dataset: indoor and outdoor terrestrial LiDAR scan

As the original goal of our study was to propose an evaluation protocol suited for real data, this dataset is only based on real data that we acquired ourselves. The resulting point clouds intensity channel can be visualised on Figure 4.6.

Stationary LiDAR station

We used the stationary LiDAR station "Leica ScanStation P40 " for which we will give a brief introduction. Once the station is settled, rays are cast 360 • horizontally around its origin and 290 • vertically (the ground area immediately underneath the station remains unobserved during the acquisition). It can acquire up to 1, 000, 000 points per second from 0.4 metre to 270 metres with a 3D position accuracy of 3 millimetres at 50 metres. In order to satisfy the condition of P GT being of higher quality than P If , we decided to acquire points from more viewpoints to generate P GT , using the same LiDAR station. We thus scanned each environment from five positions: the four vertices of an approximately regular 1.5 metres side length tetrahedron (O 1 , O 2 , O 2 , O 4 ) and its centre of mass O 5 as shown in Figure 4.7. In order to evaluate surface reconstruction algorithms in several different scenarios, we repeated this procedure for three scenes:

• "Building": an outdoor scene made of a building, a sloped road, trees and an opening to another scene called "Car park".

• "Car park": an indoor scene with pipes, partially occluded cars and open doors, including one communicating with the outdoor scene "Building".

• "Clutter": a closed indoor scene with a lot of occlusions due to a high density of objects.

Matrix format and sub-sampling

Following the definition of our protocol (see Section 4.2), we need to generate poorerquality point clouds to run reconstructions. Our LiDAR system has a spherical geometry as illustrated in

ETH3D dataset

[87] presents a two-and multi-view stereo benchmark. Their dataset contains several scenes with:

• input images at 24 Megapixel resolution on several scenes

• ground truth 3D laser scan point clouds

Generating the point clouds: Multi-View Stereo and real LiDAR

In order to generate the low-quality data, we used the OpenMVS [START_REF] Cernea | OpenMVS: Multi-View Stereo Reconstruction Library[END_REF] library to generate dense point clouds from images using the provided camera poses of three scenes (Terrace, Courtyard, Pipes) of the ETH3D train dataset. We used the DensifyPointCloud tool of OpenMVS with the standard settings, except for the following parameters: numberviews-fuse = 2, optimize = 0 and resolution-level = 4. We used the provided LiDAR point clouds as our high-quality data. A typical MVS pipeline generates a much sparser and more noisy point cloud than what a laser scan provides. It also contains more outliers. For all these reasons, we consider it relevant to carry out an evaluation on a set of

Results

In this survey, we assessed:

• RESR [START_REF] Labatut | Robust and efficient surface reconstruction from range data[END_REF] • SSD [START_REF] Calakli | SSD: Smooth signed distance surface reconstruction[END_REF] with two combinations of the octree depth and the B-spline degree parameters: (depth = 8, degree = 2) and (depth = 12, degree = 3).

• Poisson [START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF] with two values for the octree depth parameter: 8 and 11. The B-spline degree will always be 2.

• DGNN [START_REF] Sulzer | Scalable Surface Reconstruction with Delaunay-Graph Neural Networks[END_REF] • PSS [START_REF] Lafarge | Surface reconstruction through point set structuring[END_REF] with two values for the trade-off parameter: 0.1 and 0.6

• Agisoft Metashape 1.6.4 (the user manual can be found here) with Extrapolated mode and ultra high resolution.

For the ETH3D dataset only, we assess two other surface reconstruction algorithms as they are part of the OpenMVS [START_REF] Cernea | OpenMVS: Multi-View Stereo Reconstruction Library[END_REF] pipeline. Mesh reconstruction is initiated using Exploiting Visibility Information in Surface Reconstruction to Preserve Weakly Supported Surfaces (WSS) [START_REF] Jancosek | Exploiting visibility information in surface reconstruction to preserve weakly supported surfaces[END_REF]. A refinement step is then carried out using High Accuracy and Visibility-Consistent Dense Multiview Stereo (DMS) [START_REF] Vu | High accuracy and visibility-consistent dense multiview stereo[END_REF]. We thus compute the metrics on the resulting meshes from both these methods.

STRAS PC3E44

STRAS dataset

In accordance with the survey conducted and published in [START_REF] Marchand | Evaluating Surface Mesh Reconstruction of Open Scenes[END_REF] on the same dataset but with different assumptions and metrics, Table 4.3 shows that RESR achieves the best performance again on the urban environment of [START_REF] Et | eurometropole de Strasbourg[END_REF] regarding both precision and recall. As evaluating surface reconstruction from real data only is harder than using synthetic data (we do not have an exhaustive ground truth), it is a very sound validation that the metrics that we defined without access to the ground truth surface show similar tendencies to the metrics that have access to ground truth surfaces.

We also carried out a more detailed evaluation by testing several values for the main parameters of selected methods. We found that the trade-off of PSS [START_REF] Lafarge | Surface reconstruction through point set structuring[END_REF] has a big effect on the metrics. More precisely, the lower we set it (the more confidence we give to the data), the lower the precision but the higher the recall (and vice versa). A high confidence in the data results in a lot more interfaces between occupied tetrahedra and empty ones. Conversely, a higher trade-off λ gives more power to the regularisation term, resulting in fewer couples of adjacent tetrahedra being labelled differently, and so fewer triangles in the output mesh. When more confidence is given to the data, there are a lot more undesired triangles "floating" in regions of free space, which dramatically affects the precision. However, small structures might be erased from the mesh if less confidence is given to the data term, resulting in a poorer recall.

Poisson [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] and SSD [START_REF] Calakli | SSD: Smooth signed distance surface reconstruction[END_REF] are both influenced positively by an increase in the octree depth. This was expected since more points are used to reconstruct the mesh, which results in an increase in the computation time and memory footprint.

DGNN performs a lot better on the STRAS dataset than on the two others, which highlights a problem in its capacity to generalise to scenes that differ from the ones in the training set. However, its poorer F-score performance on ETH3D and ENSG is mostly due to the precision metric. DGNN often succeeds in recovering the scene features but adds too many undesired triangles in the scene. While precision, recall and ray distance provide complete information on the quality of the reconstruction, one might find it more intuitive to start by having a look at the cumulative distances shown in Figure 4.11. The precision at small range can be estimated as the area under the curve. The closer the curve is to the top left-hand corner, the better it is since this means that all the True Positives are actually very close to it. Besides, the highest value of each curve is the recall of the corresponding method so the gap between the cumulative population in the last category and the line y = 1 should be as small as possible.

ENSG dataset

The ENSG dataset, having been generated using a stationary LiDAR system, is the one containing the least amount of noise hence the overall good performance of all of the methods. In particular, we can see that the mean distance is generally a lot smaller than with other datasets even though the scenes themselves have much more complicated geometries and more occlusions. DGNN and PSS have the same problem: whilst having a high recall, denoting their capacity to recover most of the existing pieces of surface (and very accurately, given the very low mean distance metric), they connect too many regions of space with triangles lying in empty space, thus affecting their precision.

The accordance between all of these visual observations and the corresponding quantitative results given by our metrics prove their relevance.

ETH3D dataset

MVS-based point clouds are known to be noisy and contain quite a lot of outliers. This seems to have an effect on the performances of the different methods. achieves a fairly high precision on this dataset. We believe that this is because it is more capable of filtering out the noise with an 8-depth than with an 11-depth octree. That would explain why SSD 8/2 also has a better precision than SSD 12/3. However, their poorer recall indicates that more pieces of real surface have not been recovered.

The method that performs best, however, is DMS, which is not very surprising considering the fact that it is the best version of a real MVS pipeline, fed with images and not with an image-derived point cloud.

The relatively good performance of A. Metashape on ETH3D compared to the other datasets might suggest it copes pretty well with outliers. More generally, considering it is a licensed solution, we might have expected a better overall performance on at least one of our datasets.

General remarks

Overall, RESR is the method that performs best almost everywhere. DGNN shows that learning how to reconstruct large, complex and open scenes is indeed possible, but it faced generalisation problems since the metrics on ENSG and ETH3D datasets are significantly lower than those on STRAS (from which its training set was extracted). However, with RESR and DGNN being the only methods making use of sensor positions, we believe that this is an important reason behind their good results. Sensor positions give an important piece of information that neither the points themselves nor the associated normals provide.

Poisson generally performs structurally better than SSD. Small-scale differences are noticeable when changing the octree depth used by both these algorithms. PSS often reconstructs meshes very close to the real surface but also connects pieces of surface in regions of space that should remain empty. We can assume that we failed to find the right parameter settings because it was definitely the hardest algorithm to tune, but this is the best performance we managed to get.

Conclusion

Surface reconstruction is hard to evaluate since it is impossible to directly compute the difference between the real surface and a reconstructed one. It has often been assessed visually because it seems fairly intuitive to know whether a piece of surface has been accurately recovered. However, human perception can be unfair and a purely visual evaluation lacks quantitative information. In this chapter, we proposed new metrics to assess surface reconstruction. We have leveraged the visual information obtained by combining the acquired points and the associated sensor positions in order to define what we believe are more relevant metrics than the ones currently used. They imitate the process of a human being looking at and comparing the two surfaces (the real one and the reconstructed one). This goal has been achieved since our survey validates behaviours that a human can interpret by just looking at the meshes. In Section 4.4 we drew parallels between the specific visual observations and the quantitative evidence provided by our metrics that confirms them.

Our metrics enable the assessment of the completeness of the reconstructions as well as their precision both locally and globally. One can thus analyse the results from different points of view. As a relevant outcome, our survey also confirms that sensor positions are very relevant when trying to separate occupied from empty space.

As well as all these advantages that make surface reconstruction evaluation more objective, the fact that we only use raw data acquired with basic sensors makes it easy to set up a new experiment. Having access to expensive data is not a requirement. We provide a tool to make surface reconstruction evaluation easier and wish to see it used widely.

General conclusion and perspectives Conclusion

Two major challenges in surface reconstruction have been tackled in this PhD. We first presented an end-to-end distributed surface reconstruction pipeline in order to make it possible to scale up to arbitrarily-large point clouds whilst guaranteeing the watertightness of the model. Secondly, we proposed two protocols to assess the quality of any reconstruction with the aim to make assessments more rigorous and less biased than previous approaches.

Large-scale surface reconstruction

Improvements in sensor capabilities have made it possible to acquire massive point clouds representing large areas. The LiDAR HD project, led by the French Mapping Agency, aims to map all French territory at 10 points per square meter resolution. However, working with such massive data is impossible without re-thinking point cloud processing algorithms and storage. Memory footprint and computation time quickly become prohibitive when working on a standard machine, so we investigated the parallelisation of a surface reconstruction algorithm [START_REF] Caraffa | 3D Watertight Mesh Generation with Uncertainties from Ubiquitous Data[END_REF] and in particular the graph-cut optimisation step.

We generalised the work on the parallelisation of an imaged-based graph-cut optimisation [START_REF] Strandmark | Parallel and distributed graph cuts by dual decomposition[END_REF], which uses the method of Lagrange multipliers [START_REF] Boyd | Convex optimization[END_REF], to the more complex case of a 3-dimensional structure dual graph. We proceeded by splitting the input point cloud into tiles and computing its distributed Delaunay triangulation using the algorithm and data structure presented in [START_REF] Caraffa | Tile & Merge: Distributed Delaunay Triangulations for Cloud Computing[END_REF]. In order to label each tetrahedron as empty or occupied, we needed to minimise an energy function. While it has become standard to use graph-cut optimisation to solve this problem, we wanted to make it possible to split the computation on the nodes of a cluster. We split the Delaunay triangulation dual graph after duplicating the frontier nodes so that the energy function can be divided into as many sub-energies as there are tiles. However, in order to end up with a globally consistent surface, we needed frontier nodes to be labelled identically in all the sub-graphs that they belong to. We showed how to incorporate these constraints into the energy function itself and optimise the global energy function by iterating local graph-cut optimisations and Lagrange multipliers' updates. This leads to a watertight surface that is close to the one that would have been produced by the non-distributed method.

Evaluation of surface reconstruction

Evaluating surface reconstruction is difficult because of the nature of the problem: the goal is to produce a digital model that represents as best as possible a real object or environment, of which discrete information has been acquired using a noisy sensor. For complex open scenes, it is nearly impossible to define the ground truth surface as it should be the real surface itself, so we cannot directly compute the difference between this real surface and the reconstructed one.

We have shown that one possible solution is to work with synthetic data that we consider as being as close to the real surface as possible. We thus considered such a mesh to be the perfect surface that an algorithm is supposed to reproduce, which solves the problem of defining a ground-truth that we can directly compare a reconstruction with. We next had to simulate the acquisition of a point cloud from this ground-truth mesh by a sensor. In order to produce a realistic input point cloud, the simulation includes adding noise to the acquired point positions. We could then run the algorithms that we wish to assess on this simulated data. To achieve this, we developed an aerial LiDAR simulator that replicates the production of a realistic point cloud representing an urban environment. We used it on a publicly available high-quality dataset [START_REF] Et | eurometropole de Strasbourg[END_REF]. However, even after solving this problem, measuring how different the reconstruction is from the ground truth is not straightforward. Point clouds are noisy and contain occlusions so all the features of the scene might not be present in the input point cloud, or they might be scarcely sampled. In this context, cannot expect an algorithm to recover such shapes. To deal with this question, we proposed processing both the ground truth and the reconstructed meshes. We filtered out triangles from both meshes that lie further than a given threshold from the input point cloud, thus defining the "best expected surface" and its reconstructed equivalent. We repeated this process for several threshold values and computed the distance between Poisson-Disk samplings of both these filtered meshes, thus resulting in two curves of distances. This protocol offers the possibility to analyse the precision and the recall of the reconstructions locally and globally. This proved to be insightful since we observed the algorithms we tested had different behaviours depending on the range at which we assessed them.

In order to avoid having to work with synthetic data, we showed that it was possible to set a rigorous protocol to evaluate surface reconstruction exclusively with real data. We only had to make sure that we have access to two different point clouds representing the same scene, with one of them being of significantly better quality than the other: it has to be denser, have fewer occlusions and less noise. In addition, we made extensive use of sensor positions, which are easy-to-access pieces of information in our controlled setting. The combination of high-quality points and the positions from which these points have been acquired allow us to define rays along which space is supposed to be free. Each ray also provides us with a direction along which space is supposed to become occupied, close to the corresponding high-quality point. We therefore expect to find a piece of surface in this region. We proposed computing the intersections between these rays and the reconstructed model, and analysing the ones situated in observed space. These intersections that denote pieces of reconstructed surface can be classified as True Positives or False Positives depending on how close they are to the corresponding highquality point and how compatible they are with the visibility information provided by the laser ray. We defined False Negatives as pieces of surface that were expected but were not recovered, and these entities allow us to calculate precision and recall ratios. We combined these two complementary metrics with the distance from each high-quality point to the nearest intersection, which provides information on the local precision of the reconstruction. We showed that these metrics are relevant when trying to assess surface reconstruction since they allow validating and quantifying human visual assessment.

Perspectives

The end-to-end distributed surface reconstruction pipeline that we presented in Chapter 2 is an improvement on previous Delaunay-based algorithms. We have nevertheless identified a number of perspectives for our work. In particular, there is no theoretical guarantee that the algorithm converges to a unique labelling over all tiles. In other words, there might be residual disagreements between the labelling of a single tetrahedron in the different tiles that it belongs to. To guarantee watertightness, we decided that each of these inconsistently labelled tetrahedra would be labelled according to one particular tile, defined by our distributed Delaunay triangulation structure [START_REF] Caraffa | Tile & Merge: Distributed Delaunay Triangulations for Cloud Computing[END_REF]. In the case where a tetrahedron is split across an odd number of tiles, we could probably change that decision making process to a majority vote between the tiles. However, a shared tetrahedron is much more likely to belong to either two or four tiles, making this source of improvement quite limited. Another possibility to reduce the number of disagreements would be to duplicate more tetrahedra than just the ones defined by the distributed Delaunay triangulation structure [START_REF] Caraffa | Tile & Merge: Distributed Delaunay Triangulations for Cloud Computing[END_REF]. This would add more contextual information to the decision making process, helping to guess the right label of each tetrahedron. The drawback of this approach is the induced increase in the memory footprint and the extra computations carried out for each tile involved.

An ambitious avenue would be to seek a more efficient way to update the Lagrange multipliers. The algorithm, in its current form, needs quite a few iterations to converge and we might be able to save computation time by finding a more efficient strategy.

The contributions to the evaluation of surface reconstruction that we presented in Chapter 3 and 4 could also be strengthened. Our synthetic data evaluation pipeline has been tested on a mesh that represents a real scene (the city of Strasbourg in France). We can theoretically use it on any mesh representing a fictitious environment. An evident continuation of this work would therefore be to use purely handmade synthetic data in order to diversify the type of environment on which we evaluate algorithms.

We could even turn it into an open benchmark by extending to more complex scanning geometries but also to other scanning devices, in particular terrestrial, drone, or satellite platforms. That would also require incorporating the new metrics that we defined on real data, which would result in a multi-scene, multi-data and multi-metric benchmark.

A final perspective that we have identified concerns learning-based algorithms. Such methods are gaining in popularity despite still facing problems when trying to reconstruct large and complex scenes. Our work could provide a way to produce more training data and learn more efficiently thanks to our metrics.

Introduction

La reconstruction de surface à partir de nuages de points consiste à trouver le modèle continu représentant au mieux une surface réelle dont seule une information ponctuelle est disponible. Ce nuage de points peut avoir été produit à partir d'images de la scène grâce aux techniques de Stéréoscopie multi-vues (Multi-view stereo) [START_REF] Furukawa | Multi-view stereo: A tutorial[END_REF] ou de Structure acquise à partir d'un mouvement (Structure from motion) [START_REF] Ozyesil | A survey of structure from motion[END_REF]. Mais nous nous intéresserons surtout au cas où les nuages de points ont été acquis grâce à un dispositif laser de type LiDAR aérien ou terrestre permettant de scanner des zones étendues. D'autres appareils de plus courte portée sont basés sur le temps de vol (time-of-flight) [2] ou la lumière structurée (structured-light) [3].

De nombreuses contributions ont été proposées pour résoudre ce problème. Poisson [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] est une méthode qui consiste à calculer une fonction indicatrice (valant 1 à l'intérieur de l'objet et 0 en dehors) en se basant sur les positions des points du nuage ainsi que sur les normales qui leur sont associées. Un maillage triangulé correspondant à une isovalue est ensuite extrait. [START_REF] Caraffa | 3D Watertight Mesh Generation with Uncertainties from Ubiquitous Data[END_REF] calcule la triangulation de Delaunay du nuage de points et utilise les positions à partir desquelles chaque point a été acquis afin d'attribuer à chaque tétrahèdre un statut (intérieur ou extérieur). La surface est alors définie comme l'ensemble des triangles à l'interface entre un tétraèdre intérieur et un tétraèdre extérieur. Plus précisément, la segmentation des tétraèdres provient de l'optimisation par coupe de graphe d'une fonction énergie composée de deux termes :

• un terme de données qui favorise des statuts (intérieur ou extérieur) en accord avec les hypothèses faites grâce aux données (points et positions du capteur)

• un terme de régularisation qui minimise le nombre de couples de tétraèdres adjacents ayant des statuts différents.

Plus récemment, des algorithmes basés sur l'apprentissage machine ont été développés pour estimer cette segmentation intérieur-extérieur : IM-NET [START_REF] Chen | Learning implicit fields for generative shape modeling[END_REF], Occupancy Networks [START_REF] Mescheder | Occupancy networks: Learning 3d reconstruction in function space[END_REF], Convolutional Occupancy Networks [START_REF] Peng | Convolutional occupancy networks[END_REF].

Les performances des appareils d'acquisition de données cités plus haut ont beaucoup augmenté au cours des dernières années, rendant possible la production massive de données. Les ordinateurs sont quant à eux limités par leur mémoire vive et leur puissance de calcul. Il est par conséquent nécessaire de repenser les algorithmes de traitement de données et certains travaux permettent de répondre à ce besoin. De nouveaux outils et paradigmes ont vu le jour notamment le streaming [START_REF] Pajarola | Stream-processing points[END_REF] mais aussi des procédures standardisées comme le Message Passing Interface (MPI) [START_REF] Pacheco | Parallel programming with MPI[END_REF][START_REF] Gropp | Using MPI: portable parallel programming with the message-passing interface[END_REF] ou encore Spark [START_REF] Zaharia | Spark: Cluster Computing with Working Sets[END_REF], un cadre de programmation en Scala plus adapté que le paradigme MapReduce [START_REF] Dean | MapReduce: simplified data processing on large clusters[END_REF] pour les traitements itératifs et interactifs. D'autre part, des algorithmes de reconstruction de surface dédiés au cas particulier des nuages de points large échelle ont vu le jour : l'algorithme ball-pivoting [START_REF] Bernardini | The ball-pivoting algorithm for surface reconstruction[END_REF] datant de 1999 fait figure de pionnier en termes de méthode out-of-core. La méthode proposée dans [START_REF] Han | Distributed surface reconstruction from point cloud for city-scale scenes[END_REF], basée sur la triangulation de Delaunay, consiste à découper le nuage de points en tuiles et calculer les bouts de surface correspondant aux différentes tuiles en parallèle. Le problème d'une telle approche est que la surface globale résultante n'est pas forcément étanche, propriété importante des surfaces réelles.

Alors que de nombreux algorithmes utilisant un panel d'approches variées ont été développés pour reconstruire un maillage triangulé à partir d'un nuage de points, peu de méthodes d'évaluation sont aujourd'hui disponibles. L'évaluation permet pourtant de déterminer quelle surface est la plus fidèle par rapport à la réalité et idéalement, de quantifier les écarts entre les différentes méthodes. De nombreux auteurs ont simplement évalué la qualité des modèles produits visuellement, ce qui pose évidemment des problèmes d'honnêteté dus aux conflits d'intérêt et de manque de quantification. À condition de disposer d'un maillage que l'on considère idéal, Metro [1] permet de calculer la distance entre un échantillonnage de la surface reconstruite et le maillage vérité terrain (et vice versa).

Passage à l'échelle

L'une des contributions de cette thèse est l'adaptation d'un algorithme de reconstruction de surface [START_REF] Caraffa | 3D Watertight Mesh Generation with Uncertainties from Ubiquitous Data[END_REF] au cas large échelle. En particulier, nous présentons comment paralléliser l'étape d'optimisation de la fonction énergie (équation 4.4). La méthode consiste à attribuer à chaque tétraèdre t de la triangulation de Delaunay T du nuage de points le label 0 si celui-ci est vide et 1 s'il est plein. Pour cela, on minimise la fonction énergie (4.4) par coupe de graphe. On notera G = (V, E) le graphe dual dont les sommets V = {s, t} ∪ T sont la source s, le puits t et les tétraèdres de T . Les arrêtes E sont les triangles à l'interface entre toutes les paires de tétraèdres adjacents mais aussi celles qui connectent chaque tétraèdre à s et t. Chaque arrête est pondérée de manière à ce que le coût d'une coupe s-t dans le graphe dual ait la même valeur que la fonction énergie (4.4) après avoir attribué les valeurs des labels x. Pour minimiser cette énergie, nous devons donc minimiser la somme des poids des arrêtes coupées (équation 4.5).

Algorithme de référence x = arg min

(i,j)∈E c i,j x i,j (4.5) où c i,j est le poids de l'arrête (i, j) et x i,j est un Booléen indiquant si l'arrête est coupée (1 = coupée, 0 = non coupée). De plus, x i (le label du tétraèdre i) indique si le sommet i est relié à s (x i = 0) ou à t (x i = 1) après une coupe. On a alors : x s = 0, x t = 1 et ∀(i, j) ∈ E, x i,j = |x i -x j |.

Algorithme distribué

Pour paralléliser l'étape d'optimisation de la fonction énergie, nous proposons de découper le nuage de points d'entrée P en tuiles. Nous noterons P K = (P k ) k∈K sa décomposition en |K| sous-ensembles disjoints P k , où K représente l'ensemble des indices des tuiles.

La triangulation de Delaunay T du nuage de points peut donc être décomposée elleaussi en |K| sous-graphes G k = (V k , E k ) en utilisant l'algorithme et la structure de [START_REF] Caraffa | Tile & Merge: Distributed Delaunay Triangulations for Cloud Computing[END_REF]. Certains tétraèdres seront partagés par deux tuiles donc le noeud qui les représente apparaîtra dans deux sous-graphes. Le problème principal est alors de garantir que les tétraèdres à l'interface entre deux tuiles obtiendront le même label dans les deux sousgraphes. Pour tout k ∈ |K|, on dénote respectivement par x k i,j et x k i les valeurs de x i,j et x i dans le sous-graphe k. En outre, on définit c k i,j = c i,j |{k:(i,j)∈E k }| : le poids de l'arrête (i, j) normalisé par le nombre de tuiles auxquelles cette arrête appartient. Avec cette décomposition, la fonction optimisée devient :

f (x) = (i,j)∈E c i,j x i,j = k∈K (i,j)∈E k c k i,j x k i,j (4.6)
avec la condition que les labels x k i doivent avoir la même valeur dans tous les sous-graphes dans lesquels ils appartiennent :

∀ k, l ∈ K 2 , ∀ i ∈ V k ∩ V l , x k i -x l i = 0 (4.7)
Pour transformer ce problème en sous-problèmes en utilisant notre structure de graphe tuilé, on optimise par coupe de graphe chaque sous-graphe G k par rapport aux variables

x k i,j . Pour garantir la condition exprimée par l'équation 4.7, on ajoute un terme de pénalité à l'énergie en utilisant la dualité Lagrangienne :

L(x, λ) = k∈K (i,j)∈E k c k i,j x k i,j + k∈K l>k i∈V k ∩V l λ k,l i (x k i -x l i ) (4.8) avec λ = {λ k,l i : (k, l) ∈ K 2 , i ∈ V k ∩ V l }.
En définissant, pour chaque index de tuile k ∈ K : Comme la fonction duale de Lagrange g(λ) (équation 4.11) est concave [START_REF] Boyd | Convex optimization[END_REF], on peut trouver une solution optimale au problème (4.6+4.7) en maximisant g(λ) grâce à une méthode ascendante. On itère alors les phases de :

L k (x k , λ) = (i,j)∈E k c k i,j x k i,j + l̸ =k i∈V k ∩V l    1 if l > k -1 otherwise    λ k,l i x k i (4.
• Résolution des sous-problèmes d'optimisation par coupe de graphe L k (x k , λ) for x k .

• Mise à jour de λ pour aller dans la direction de plus forte montée.

Nous montrons expérimentalement que l'implémentation de cet algorithme avec Spark [START_REF] Zaharia | Spark: Cluster Computing with Working Sets[END_REF] et son exécution sur un cluster d'ordinateurs permet de diminuer le temps de calcul (pas de manière linéaire cependant) en augmentant le nombre de noeuds. Cet algorithme garantit l'étanchéité du modèle grâce à la segmentation unique des tétraèdres de la triangulation. Enfin, le modèle produit par la méthode distribuée converge vers un modèle proche de celui produit par l'algorithme de référence.

Évaluation

Évaluer la reconstruction de surface est un problème difficile par nature étant donné que la seule vérité terrain est la surface réelle elle-même. Comme l'objectif est de produire un modèle numérique représentant cette surface, il est impossible de calculer directement l'"écart" entre les deux. Deux alternatives ont alors été envisagées.

Utilisation de données synthétiques

La ville et eurométropole de Strasbourg a financé la production d'un maillage de son territoire [START_REF] Et | eurometropole de Strasbourg[END_REF] et nous avons utilisé ce jeu de données pour définir des métriques d'évaluation propres à ce cas de figure.

Notre protocole repose premièrement sur la disponibilité d'une surface vérité terrain M GT . Deuxièmement, il faut un nuage de points P représentant une numérisation de M GT réalisée à l'aide d'un simulateur de LiDAR comme celui présenté dans cette thèse. Enfin, nous avons besoin d'un maillage triangulé M E reconstruit grâce à une méthode que l'on souhaite évaluer. En pratique, seule une partie limitée de la vérité terrain est couverte par l'acquisition. Si nous pouvons attendre d'un algorithme qu'il remplisse les trous entre les points, on ne peut cependant pas espérer qu'il devine la forme de la surface loin de tout point d'acquisition. On définit alors la partie "reconstructible" M α GT de M GT au seuil de α. On peut produire M α GT en supprimant tous les triangles de M GT pour lesquels les trois sommets se trouvent à une distance supérieure à α d'un point de P. M α GT est donc la meilleure surface que l'on peut attendre d'un algorithme, au seuil de α. Toute partie reconstruite à une distance supérieure à α d'un point de P sera également volontairement ignorée pour ne pas que les hypothèses faites par les algorithmes impactent l'évaluation. On définira donc M α E comme la partie pertinente à évaluer de M E . La La précision décrit à quel point la surface reconstruite est en accord avec la vérité terrain. Le rappel donne une information d'exhaustivité de la reconstruction c'est-à-dire qu'il permet de savoir à quel point les morceaux de surface de la vérité terrain sont également présents dans la reconstruction.

En calculant ces deux grandeurs pour différentes valeurs de α, on obtient deux courbes qui permettent d'analyser les performances des algorithmes localement (pour les faibles valeurs de α) et plus globalement (pour les hautes valeurs de α). Les différents comportements observés pour les différentes valeurs de α valident la pertinence de ce raisonnement.

Utilisation de données réelles

Une continuation assez naturelle du travail réalisé sur l'évaluation de la reconstruction de surface était de mettre en place un protocole remplissant le même objectif mais sans avoir recours à des données synthétiques. On appelle "données réelles" des données acquises dans le monde réel avec de vrais capteurs. Ceci inclut les images RGB-D, les acquisitions LiDAR, etc. Nous allons évaluer la reconstruction de scènes réelles à partir de données réelles de qualité significativement meilleure que celles qui ont servi à la reconstruction. Même si cette idée est assez classique, la contribution majeure de ce travail réside dans la manière dont nous proposons de calculer la différence entre la surface reconstruite et les données réelles de haute qualité en relevant les incompatibilités et les morceaux de surface manquants. De plus, pour la première fois, nous proposons d'utiliser les positions du capteur pour l'évaluation. Associées aux points de haute qualité, elles permettent de définir un rayon le long duquel l'espace est censé être vide. Tout morceau de surface se trouvant sur le chemin du rayon laser sera donc pénalisé car incompatible avec cette information. Un morceau de surface sera considéré manquant si aucune intersection n'est trouvée entre le rayon laser et la surface reconstruite ou si une telle intersection est trouvée à une distance supérieure à un certain seuil du point de haute qualité correspondant, et située après, le long du rayon. Cela nous permet de classer les intersections trouvées dans des zones observées. Les intersections correctes correspondent à des Vrais positifs, celles incorrectes sont des Faux positifs et les morceaux de surface manquants sont des Faux négatifs. On peut alors calculer la précision et le rappel et combiner ces deux métriques à la distance moyenne des intersections correctes aux points de haute qualité. La Figure 4.15 permet de visualiser ces définitions sur un exemple simple.

Ce protocole d'évaluation, basé sur l'information de visibilité venant des positions du capteur, a été testé sur trois jeux de données différents, contenant chacun trois scènes. Les métriques permettent de confirmer des intuitions visuelles, participant ainsi à valider la pertinence de ce mode d'évaluation.

Conclusion

Cette thèse de doctorat contribue à deux aspects fondamentaux de la reconstruction de surface. Premièrement, un algorithme distribué de bout en bout a été proposé afin de répondre à la demande de traitement des nuages de points large échelle. Il offre notamment la garantie de reconstruire une surface étanche. Deuxièmement, nous avons présenté deux protocoles d'évaluation permettant de quantifier la qualité d'une surface reconstruite. L'un nécessite d'avoir accès à un maillage vérité terrain alors que l'autre fait uniquement usage de données provenant d'un capteur.

Perspectives

L'algorithme de reconstruction de surface distribué de bout en bout constitue une avancée importante par rapport aux approches basées sur la triangulation de Delaunay. Plusieurs pistes d'améliorations sont toutefois envisageables. En particulier, il n'y a pas de garantie théorique sur la convergence de l'algorithme vers une segmentation unique des tétraèdres parmi toutes les tuiles. Pour garantir l'étanchéité du modèle en cas de désaccord, un tétraèdre se verra octroyer le label d'une tuile particulière. Dans le cas où il est partagé par un nombre impair de tuiles, il pourrait être judicieux de prendre cette décision par vote majoritaire entre les tuiles concernées. Une autre possibilité consiste à réduire le nombre de désaccords en dupliquant plus de tétraèdres que ceux définis par la structure de triangulation de Delaunay de [START_REF] Caraffa | Tile & Merge: Distributed Delaunay Triangulations for Cloud Computing[END_REF]. Ceci ajouterait plus de contexte rendant la prise de décision plus fiable. En revanche, l'empreinte mémoire et le temps de calcul à chaque itération se verraient être augmentés à cause de ces tétraèdres dupliqués supplémentaires.

Dans le but d'améliorer la capacité de l'algorithme à converger, réduisant ainsi le temps de calcul, il serait judicieux de trouver une stratégie de mise à jour des multiplicateurs de Lagrange plus efficace.

Les contributions à l'évaluation de la reconstruction de surface peuvent elles aussi êtres approfondies. Le protocole d'évaluation basé sur les données synthétiques a été

2 . 6

 26 ?!) of friendship. I know my grandfather Marcel Gaillard would've loved to attend my PhD defence. I know it would've meant a lot to you. Thank you for all the inspiration, you are forever showing me what tenacity and adaptation could mean! When I moved to Paris to start the PhD, I felt very lucky that my brother Théo was living nearby in Sarcelles. It was the first time were living in the same place as adults and I have some great memories of the weekends we spent together: breakfasts, runs, and watching Rocky ("See ya!")! I miss that stuff. Thank you for all your encouragement and good vibes, brother! My parents Anne Gaillard and Pierre Marchand have given me everything, starting with love. They put way more effort into my education than into any personal aspiration of theirs. They did not worry about their own success, but they cared about mine. I'm not a self-made person, I owe them a lot.4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 STRAS dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2 ENSG dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.3 ETH3D dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.4 General remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . hard watertightness (left) and soft watertightness (right). The soft-watertight surface has a border (materialised by the blue dots) but only on the boundary of the domain. When trying to model urban environments, the red piece of surface from the hard-watertight surface does not have any significant meaning so soft-watertightness is often better suited to open scenes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 2-dimensional example of the proposed approach. . . . . . . . . . . . . . .2.2 The proposed distributed surface reconstruction workflow. P denotes the input point set, P k the tiled point set, P I k and P S k are the tiled point set with PCA information and the simplified one, T k the tile-triangulations, T m k the tile-triangulation with the mass score, T x,t k the labelled tile-triangulation at iteration t, x k the updated mixed cells during the optimisation and M k the final mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Weight update after a single split in the dual graph. We focus on node i. Nodes in bold are the ones that have been duplicated after the split. s and t represent the source and the sink. . . . . . . . . . . . . . . . . . . . . . . 2.4 2D visualisation of the full triangulation (left) being split into 2 tiles (middle and right) with local cells in white and mixed cells in blue. Mixed cells are duplicated in both sub-graphs. . . . . . . . . . . . . . . . . . . . . . . . 2.5 Ratio between the energy of the labelling given by the distributed graphcut E d and the reference graph-cut with one thread E r for different initialisations of the Lagrangian step τ 0 = 1, 2, 3, 5, 10, 20 (left) and different octree depths 1, 2, 3, 4 (right). . . . . . . . . . . . . . . . . . . . . . . . . Execution time (left) and strong scaling factor t 1 core n.tn cores (right) as a function of the number of cores with 4, 3, 2 and 1 core(s)/executor (12Go RAM/executor) on a 10 million points data set. . . . . . . . . . . . . . . . 2.7 Evolution of the resulting mesh along iterations on an urban scene. Each line shows two different parts of the mesh. The first line shows the result with the local graph-cut without iterating, the second line corresponds to iteration n • 3, the third line to iteration n • 15 and the fourth line to iteration n • 30. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Reconstruction result on a 350 million points data set. Tiling is visualised through a random colour per tile. . . . . . . . . . . . . . . . . . . . . . . .

  Strasbourg open dataset . . . . . . . . . . . . . 3.4 Global frame and flight trajectory. . . . . . . . . . . . . . . . . . . . . . . . 3.5 Parallel-line scanning pattern: zoom in on the local frame and laser pulse direction. The red-blured area represents the field of view. . . . . . . . . . 3.6 Elliptical scanning pattern: zoom in on the local frame and laser pulse direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Result of the different algorithms evaluated on a crossroad of the scene . . 3.8 Mean Precision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Mean Recall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 The importance of sensor positions: the dashed part of the reconstructed surface can be identified as wrong by making use of sensor positions when the high-quality point cloud does not provide enough information. . . . . . 4.2 These four cases would be evaluated in the same way by a basic point-tomesh distance. However, they are very different in terms of what a human being would be able to see from the sensor position. . . . . . . . . . . . . .4.3 Toy example to visualise the definitions of the metrics. The real surface has been scanned from two positions: O 1 and O 2 .

4. 1

 1 Values of experimental parameters used. . . . . . . . . . . . . . . . . . . . 4.2 Raw numerical results for d max = 50 cm (MD stands for mean distance, P for precision, R for recall and F1 for F-score) . . . . . . . . . . . . . . . . . 4.3 Average numerical results on the three scenes from STRAS dataset sorted by decreasing F-score for d max = 50 cm . . . . . . . . . . . . . . . . . . . . 4.4 Average numerical results on the three scenes from ENSG dataset sorted by decreasing F-score for d max = 20 cm . . . . . . . . . . . . . . . . . . . . 4.5 Average numerical results on the three scenes from ETH3D dataset sorted by decreasing F-score for d max = 20 cm . . . . . . . . . . . . . . . . . . . .

  (a) Tiled input point cloud with multiple acquisitions. (b) Tiled Delaunay triangulation. The local cells are shown in white and the mixed ones in blue. The adjacency graph is drawn in red. (c) Surface extraction. The cells in red and green are labelled as empty and occupied respectively, the light blue line denotes the extracted surface.

Figure 2 .

 2 Figure 2.1: 2-dimensional example of the proposed approach.

:

  Point set, : DT, : Cell set, : Mesh, : Union operator, : C ++ call.

Figure 2 . 2 :

 22 Figure 2.2: The proposed distributed surface reconstruction workflow. P denotes the input point set, P k the tiled point set, P I k and P S k are the tiled point set with PCA information and the simplified one, T k the tile-triangulations, T m k the tile-triangulation with the mass score, T x,t k the labelled tile-triangulation at iteration t, x k the updated mixed cells during the optimisation and M k the final mesh.

  subset of the vertices of T k . The vertices of T k are said to be local in T k if they belong to P k and foreign in T k if they belong to another tile P l with l ̸ = k. By extension, cells and facets of T k are local or foreign in T k if all their vertices are local in T k or foreign in T k respectively, and mixed otherwise.

  (a) Dual graph before split (b) Dual sub-graphs after split

Figure 2 . 3 :

 23 Figure 2.3: Weight update after a single split in the dual graph. We focus on node i. Nodes in bold are the ones that have been duplicated after the split. s and t represent the source and the sink.

Figure 2 . 4 :

 24 Figure 2.4: 2D visualisation of the full triangulation (left) being split into 2 tiles (middle and right) with local cells in white and mixed cells in blue. Mixed cells are duplicated in both sub-graphs.

Figure 2 . 5 :

 25 Figure 2.5: Ratio between the energy of the labelling given by the distributed graphcut E d and the reference graph-cut with one thread E r for different initialisations of the Lagrangian step τ 0 = 1, 2, 3, 5, 10, 20 (left) and different octree depths 1, 2, 3, 4 (right).

Figure 2 . 6 :

 26 Figure 2.6: Execution time (left) and strong scaling factor t 1 coren.tn cores (right) as a function of the number of cores with 4, 3, 2 and 1 core(s)/executor (12Go RAM/executor) on a 10 million points data set.

Figure 2 . 7 :

 27 Figure 2.7: Evolution of the resulting mesh along iterations on an urban scene. Each line shows two different parts of the mesh. The first line shows the result with the local graph-cut without iterating, the second line corresponds to iteration n • 3, the third line to iteration n • 15 and the fourth line to iteration n • 30.

Figure 2 . 8 :

 28 Figure 2.8: Reconstruction result on a 350 million points data set. Tiling is visualised through a random colour per tile.

Figure 3 .

 3 2 allows to visualise the computation of M α GT and M α E for a given α.

Figure 3 . 2 :

 32 Figure 3.2: Visualisation of the computation of the best expected surface for a small value of α. Starting from the input point cloud, the ground-truth mesh and two reconstructed meshes R 1 and R 2 , we compute the best expected surface and the relevant parts to assess R α 1 and R α 2 by removing all triangles lying further than α from a point of the input point cloud.

Figure 3 . 3 :

 33 Figure 3.3: Virtual environment from Strasbourg open dataset

Figure 3 . 5 :

 35 Figure 3.5: Parallel-line scanning pattern: zoom in on the local frame and laser pulse direction. The red-blured area represents the field of view.

Figure 3 . 6 : 2 z

 362 Figure 3.6: Elliptical scanning pattern: zoom in on the local frame and laser pulse direction.

  Poisson N: Poisson [13] using Neumann boundary condition • Poisson D: Poisson [13] using Dirichlet boundary condition

Figure 3 Figure 3 . 7 :

 337 Figure 3.7: Result of the different algorithms evaluated on a crossroad of the scene

Figure 3 . 8 :

 38 Figure 3.8: Mean Precision

  Figure 4.1.

Figure 4 . 1 :

 41 Figure 4.1: The importance of sensor positions: the dashed part of the reconstructed surface can be identified as wrong by making use of sensor positions when the highquality point cloud does not provide enough information.

Figure 4 . 2 :

 42 Figure 4.2: These four cases would be evaluated in the same way by a basic point-to-mesh distance. However, they are very different in terms of what a human being would be able to see from the sensor position.

Figure 4 . 3 :

 43 Figure 4.3: Toy example to visualise the definitions of the metrics. The real surface has been scanned from two positions: O 1 and O 2 . A given real laser ray (the thick one) was cast from O 1 and it hit the real surface at Gt. Note that the position of the intersection might be noisy, hence the shift between the real surface and the high-quality point cloud.We compute all intersections between the associated virtual ray (i.e. the extension of the real laser ray) and the reconstructed surface. In this case, it results in six intersections I 1,...,6 . The closest intersection to Gt happens to be I 2 so the "ray distance" metric for that particular ray is the distance (Gt, I 2 ). If (Gt, I 2 ) < d max , Gt is to be counted as a TP, otherwise it will not be taken into account in the evaluation since it is situated after Gt. Note that this piece of surface might still be evaluated thanks to another ray (as one emanating from O 2 , for example). Besides, we found one intersection I 1 on the way to the closest intersection I 2 which is counted as a FP.

Figure 4 . 4 :

 44 Figure 4.4: Strasbourg scene (mesh in grey, point cloud in blue)

Figure 4 .

 4 [START_REF] Furukawa | Multi-view stereo: A tutorial[END_REF] gives an example of such a situation.

Figure 4 . 5 :

 45 Figure 4.5: Left: Mesh, Centre: P If , Right: P GT . The facades parallel to the direction of the plane are a lot more occluded in P If .

Figure 4 . 7 :

 47 Figure 4.7: Tetrahedron-like viewpoints. (O 1 , O 2 , O 2 , O 4 ) form an approximately regular 1.5 metres side length tetrahedron and O 5 is positioned at its centre of mass. The black rectangle represents a room in which we installed our stationary LiDAR.

Figure 4 . 8 .

 48 The horizontal resolution and vertical resolution of the scanner define a fixed number w of values for ϕ and a fixed number h of values for θ respectively. Laser rays are thus cast in the directions given by every pair of angles (ϕ, θ). We can thus represent the points as a matrix of height h and width w and then index all points by their (i, j) ∈ [|1, h|] × [|1, w|] coordinates. Figure4.6 shows the intensity of the returns in this matrix-like format. Each raw acquired point cloud did not fit into the memory of our machine so we down-sampled them by keeping odd-indexed points as illustrated in Figure 4.8. Starting from the raw point clouds with origins centred on O 1 , O 2 , O 3 O 4 and O 5 , we down-sampled them all (roughly four times) following this matrix-based scheme and P GT is the union of these five four-time down-sampled point clouds. We repeated this matrix-based down-sampling scheme on the point cloud centred on O 5 to generate P If .

Figure 4 . 8 :

 48 Figure 4.8: Matrix format-based subsampling. Left: spherical representation of a point cloud acquired with the stationary LiDAR station "Leica ScanStation P40 ". Each point acquired is defined by its spherical coordinate angles (ϕ, θ). Right: Matrix format-based sub-sampling: only points represented with a black outskirt are kept in the sub-sampled point cloud.

Figure 4 . 6 :

 46 Figure 4.6: Equirectangular projection of the LiDAR points. Top: Outdoor scene ("Building") from O 3 -Middle: indoor scene ("Car park") from O 2 -Bottom: Indoor Scene ("Clutter") from O 1 . "Building" and "Car park" share some space thanks to what can be seen through the open door in the middle of both images.

  LiDAR Point Cloud With Visibility (d) Image Point Cloud With Visibility

Figure 4 . 9 :

 49 Figure 4.9: Point Clouds with Visibility Information: Terrestrial point clouds (a, b) from the Terrace scene of ETH3D [87]. We visualise some of the sensor positions and lines-of-sight (c, d).

Figure 4 . 10 :

 410 Figure 4.10: Images of the three scenes from the ETH3D dataset we used. Top: Outdoor scene ("Courtyard") -Middle: indoor scene ("Pipes") -Bottom: Outdoor Scene ("Terrace").

Figure 4 . 11 :

 411 Figure 4.11: Cumulative distances over the three scenes from STRAS dataset for d max = 50 cm.

Figure 4 .

 4 Figure 4.12 shows the meshes reconstructed by every assessed algorithm on the Park-

Figure 4 . 12 :

 412 Figure 4.12: Reconstructed meshes from the ENSG Parking Lot scene.

  x∈{0,1} |T | E data (x) + α E prior (x) (4.4)

9 )

 9 On a alors :L(x, λ) = k∈K L k (x k , λ) (4.10)La Figure4.13 permet de visualiser la mise à jour des poids du graphe dual sur un cas simple avec une seule séparation. On définit la fonction duale de Lagrange : g(λ) = min x L(x, λ)(4.11)

  (a) Graphe dual avant séparation (b) Sous-graphes duaux après séparation

Figure 4 . 13 :

 413 Figure 4.13: Mise à jour des poids du graphe dual après une séparation. On analyse le noeud i en particulier. Les noeuds en gras sont ceux qui ont été dupliqués au cours de la séparation. s et t représentent la source et le puits.

Figure 4 .

 4 14 permet de visualiser le calcul de M α GT et M α E pour un α donné. Nous proposons d'évaluer la qualité de la reconstruction M E en calculant la distance moyenne entre des échantillonnages de type Poisson-disk de M α E et M α GT . En notant P R M α GT and P R M α E les échantillonnages de type Poisson-Disk de rayon R correspondant à M α GT and M α E , et grâce à la définition de la distance point-vers-maillage (équation 4.12), on calcule alors deux métriques (équations 4.13 and 4.14). ∀p ∈ R 3 , d (p, M) = min q∈M d(p, q) (4.12)

Figure 4 . 14 :

 414 Figure 4.14: Visualisation du calcul de la partie "reconstructible" des différents maillages pour un α proche de la valeur de l'espacement moyen entre deux points du nuage. À partir du nuage de points, du maillage vérité terrain et de deux maillages reconstruits R 1 et R 2 , on calcule la partie reconstructible de la vérité terrain, et les parties pertinentes à évaluer R α 1 et R α 2 en enlevant les triangles éloignés de plus de α d'un point du nuage de points.

Figure 4 . 15 :

 415 Figure 4.15: Exemple simple pour visualiser les définitions des métriques. La surface réelle a été scannée depuis deux positions : O 1 et O 2 . Un rayon laser réel particulier (le plus épais) a été lancé depuis O 1 et il a intersecté la surface réelle au point Gt. Il est à noter que la position de l'intersection peut être bruitée, d'où le décalage entre la surface réelle et le nuage de points de haute qualité. On calcule toutes les intersections entre le rayon virtuel associé (i.e. l'extension du rayon laser réel) et la surface reconstruite. Dans le cas présent, on trouve six intersections I 1,...,6 . L'intersection la plus proche de Gt se trouve être I 2 donc la métrique de "distance rayon" pour ce rayon particulier est la distance (Gt, I 2 ). Si (Gt, I 2 ) < d max , Gt doit être compté comme un Vrai Positif, sinon il ne sera pas pris en compte dans cette évaluation puisqu'il se situe après Gt. Par ailleurs, on rencontre une intersection I 1 sur le chemin vers l'intersection la plus proche I 2 , ce qui est compté comme un Faux positif.

  .[START_REF] Kazhdan | Screened poisson surface reconstruction[END_REF] Mise à jour des poids du graphe dual après une séparation. On analyse le noeud i en particulier. Les noeuds en gras sont ceux qui ont été dupliqués au cours de la séparation. s et t représentent la source et le puits. . . . . . 4.14 Visualisation du calcul de la partie "reconstructible" des différents maillages pour un α proche de la valeur de l'espacement moyen entre deux points du nuage. À partir du nuage de points, du maillage vérité terrain et de deux maillages reconstruitsR 1 et R 2 ,on calcule la partie reconstructible de la vérité terrain, et les parties pertinentes à évaluer R α 1 et R α 2 en enlevant les triangles éloignés de plus de α d'un point du nuage de points. . . . . . . . 4.15 Exemple simple pour visualiser les définitions des métriques. La surface réelle a été scannée depuis deux positions : O 1 et O 2 . Un rayon laser réel particulier (le plus épais) a été lancé depuis O 1 et il a intersecté la surface réelle au point Gt. Il est à noter que la position de l'intersection

). Besides, we found one intersection I 1 on the way to the closest intersection I 2 which is counted as a FP. . . . . . . 4.4 Strasbourg scene (mesh in grey, point cloud in blue) . . . . . . . . . . . . . 4.5 Left: Mesh, Centre: P If , Right: P GT . The facades parallel to the direction of the plane are a lot more occluded in P If . . . . . . . . . . . . . . . . . .

4.7 Tetrahedron-like viewpoints. (O 1 , O 2 , O 2 , O 4 ) form an approximately regular 1.5 metres side length tetrahedron and O 5 is positioned at its centre of mass. The black rectangle represents a room in which we installed our stationary LiDAR. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Matrix format-based subsampling. Left: spherical representation of a point cloud acquired with the stationary LiDAR station "Leica ScanStation P40 ". Each point acquired is defined by its spherical coordinate angles (ϕ, θ). Right: Matrix format-based sub-sampling: only points represented with a black outskirt are kept in the sub-sampled point cloud. . . . . . . . . . . . 4.6 Equirectangular projection of the LiDAR points. Top: Outdoor scene ("Building") from O 3 -Middle: indoor scene ("Car park") from O 2 -Bottom: Indoor Scene ("Clutter") from O 1 . "Building" and "Car park" share some space thanks to what can be seen through the open door in the middle of both images. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Point Clouds with Visibility Information . . . . . . . . . . . . . . . . . . . 4.10 Images of the three scenes from the ETH3D dataset we used. Top: Outdoor scene ("Courtyard") -Middle: indoor scene ("Pipes") -Bottom: Outdoor Scene ("Terrace"). . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Cumulative distances over the three scenes from STRAS dataset for d max = 50 cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Reconstructed meshes from the ENSG Parking Lot scene. . . . . . . . . . 4peut être bruitée, d'où le décalage entre la surface réelle et le nuage de points de haute qualité. On calcule toutes les intersections entre le rayon virtuel associé (i.e. l'extension du rayon laser réel) et la surface reconstruite. Dans le cas présent, on trouve six intersections I 1,...,6 . L'intersection la plus proche de Gt se trouve être I 2 donc la métrique de "distance rayon" pour ce rayon particulier est la distance (Gt, I 2 ). Si (Gt, I 2 ) < d max , Gt doit être compté comme un Vrai Positif, sinon il ne sera pas pris en compte dans cette évaluation puisqu'il se situe après Gt. Par ailleurs, on rencontre une intersection I 1 sur le chemin vers l'intersection la plus proche I 2 , ce qui est compté comme un Faux positif. . . . . . . . . . . . . . . . . . . . . . . .

  of a volumetric discretisation of the object Kolluri et al.[START_REF] Kolluri | Spectral Surface Reconstruction from Noisy Point Clouds[END_REF]

		✓				tetrahedra
	Labatut et al. [21]	✓				tetrahedra	l.o.s
	Holenstein et al. [25]	✓	✓	✓	✓	voxels + tiles	l.o.s
	Caraffa et al. [23]	✓				tetrahedra	l.o.s
	Han et al. [58]	✓			✓	tetrahedra + tiles	l.o.s
			Signed-distance function
	MPU [27]	✓				octree	n i	✓
	SSD [28]	✓				octree	n i	✓
	DeepSDF [29]	✓				✓		✓
		Unsigned-distance function
	Hornung et al. [30]	✓				octree
	Randrianarivony et al. [51]				✓
			Primitive-based
	PolyFit [31]	✓				primitives		✓
				MLS-based
	Out-of-core MLS [40]		✓			octree	n i
	Streaming MLS [59]		✓	✓	✓	octree	n i
			Refinement (iterative)
	Ball-Pivoting [44]		✓	✓				✓
	Allegre et al. [41]	✓	✓			kd-tree	|n i |
	Table 1.1: Categorisation of surface reconstruction in terms of the properties that are
	guaranteed as regards watertightness (WT), scaling issue i.e. out-of-core (OOC), stream-
	ing (SG) or parallel (PL) implementation as well as the use of an additional localisation
	structure (ALS) ; learning (LG) capabilities, the necessity of additional features (AF) as
	normals: "n i " or lines of sight: "l.o.s" in the computation ; and the open-source (OS)
	availability of the code.

  T A , E A = DDT (P S ) // Mass computation 9 for k ∈ K do in parallel

	Algorithm 1: Distributed surface reconstruction
		Input: Point set P
		// Tiling
	1 for p ∈ P do in parallel
	2	k ← T ileId(p)
	3	P k ← P k {p}
	4 end
		// PCA/Sub-sampling
	5 for k ∈ K do in parallel
	6	P I	k , P S k ← P CA(P k )
	7 end
		// Delaunay Triangulation (DT)

[START_REF] Kinsey | Surfaces[END_REF]

Table 3 .

 3 1: Values of experimental parameters we used (parallel-line pattern).

	Description

Table 4 .

 4 3 (d max = 50 cm) 2: Raw numerical results for d max = 50 cm (MD stands for mean distance, P for precision, R for recall and F1 for F-score)

	Method	TP	FP	TP + FN MD (cm) P (%) R (%) F1 (%)
	RESR	1278734 57089	1365120	6.58	95.73 93.67	94.69
	DGNN	1246228 72899	1365120	6.79	94.47 91.29	92.85
	Poisson 11/2 1267301 128888 1365120	8.91	90.77 92.83	91.79
	PSS 0.6	1217138 81549	1365120	7.39	93.72 89.16	91.38
	SSD 12/3	1271969 172285 1365120	9.56	88.07 93.18	90.55
	SSD 8/2	1203150 167204 1365120	13.04	87.80 88.14	87.97
	Poisson 8/2	1174069 151134 1365120	12.77	88.60 86.00	87.28
	A. Metashape 1033179 195725 1365120	19.41	84.07 75.68	79.66
	PSS 0.1	1280230 686825 1365120	7.47	65.08 93.78	76.84

Table 4 .

 4 3: Average numerical results on the three scenes from STRAS dataset sorted by decreasing F-score for d max = 50 cm

		STRAS (d max = 50 cm)
	Method	mean distance (cm) precision (%) recall(%) F-score(%)
	RESR	5.98	96.68	94.88	95.77
	DGNN	6.13	96.08	92.56	94.29
	Poisson 11/2	7.99	93.19	94.31	93.75
	PSS 0.6	6.67	95.08	91.35	93.17
	SSD 12/3	8.53	90.72	94.63	92.63
	SSD 8/2	11.57	90.11	90.48	90.30
	Poisson 8/2	11.49	91.11	88.39	89.73
	A. Metashape	16.21	87.92	80.63	84.11
	PSS 0.1	6.76	72.72	95.00	82.28
	mean methods	9.04	90.40	91.36	90.67

Table 4 .

 4 4: Average numerical results on the three scenes from ENSG dataset sorted by decreasing F-score for d max = 20 cm ing Lot scene (part of the ENSG dataset). One can fairly easily interpret the performance achieved by these methods by analysing the type of mistake they made on the corresponding scene. RESR succeeds at reconstructing most of the visible parts, and very few undesired triangles lie in free space (most of them are connecting the pipes to the wall and the ceiling). At first glance, Poisson 11/2 reconstruction seems to be a lot more accurate than Poisson 8/2 so it might not be obvious why they have the same precision. This situation actually shows the interest of the mean distance metric. While the two reconstructed models are structurally the same, the difference between them is visible at close range: under the threshold d max . Consequently, Poisson 8/2 has a much higher mean distance than Poisson 11/2 but achieves a similar precision. The same kind of argument holds for explaining the relatively poor performance of SSD 12/3 and A. Metashape: whilst being locally more accurate than Poisson 8/2, the meshes are structurally not in accordance with the visibility information provided by the high-quality point cloud. The precision metric is dramatically affected by large portions of surface lying in free space.

		ENSG (d max = 20 cm)
	Method	mean distance (cm) precision (%) recall(%) F-score(%)
	RESR	0.45	93.10	95.99	94.51
	Poisson 11/2	0.90	78.38	96.96	86.22
	Poisson 8/2	2.72	78.66	88.16	83.05
	SSD 12/3	1.27	72.41	96.04	82.23
	A. Metashape	1.63	79.00	83.95	81.37
	SSD 8/2	3.04	71.95	87.77	78.96
	DGNN	0.49	52.99	96.22	68.28
	PSS 0.6	0.54	28.44	97.92	43.94
	PSS 0.1	0.54	22.30	98.19	36.18
	mean methods	1.29	64.14	93.47	72.75

  The ones performing best on ENSG and STRAS seem to struggle more, and surprisingly, Poisson 8[START_REF] Kazhdan | Poisson surface reconstruction[END_REF]

Table 4 .

 4 5: Average numerical results on the three scenes from ETH3D dataset sorted by decreasing F-score for d max = 20 cm

		ETH3D (d max = 20 cm)
	Method	mean distance (cm) precision (%) recall(%) F-score(%)
	DMS	2.04	95.33	93.72	94.47
	Poisson 11/2	2.56	93.91	94.04	93.96
	RESR	2.62	94.12	92.65	93.29
	A. Metashape	2.57	95.09	91.18	93.00
	WSS	2.66	91.31	93.79	92.51
	SSD 12/3	2.66	90.65	94.44	92.48
	SSD 8/2	4.21	93.38	89.03	91.13
	Poisson 8/2	3.99	94.66	86.94	90.59
	DGNN	2.61	79.47	93.63	85.95
	PSS 0.6	2.52	67.41	94.63	78.52
	PSS 0.1	2.56	53.79	95.16	66.99
	mean methods	2.82	86.28	92.66	88.44

Source code is available at: GitHub/SurfaceReconEval

https://aws.amazon.com/ec2/instance-types/m5/

Acknowledgements

Three years have passed since I arrived in Paris and stepped into the IGN building for the first time. I would like to express my gratitude to all the people who gave me support and enriched my experience.

This list has to start with my supervisors Bruno Vallet and Lau-testé sur un maillage représentant une vraie scène (la ville de Strasbourg, en France). Il est théoriquement possible de l'utiliser avec un maillage représentant une scène imaginaire. Nous pourrions alors envisager d'utiliser des scènes créées à la main afin de diversifier le type d'environnement sur lesquels nous évaluons les algorithmes.

Nous pourrions même en faire un "benchmark" ouvert en incluant des géométries d'acquisition plus complexes mais aussi en proposant d'autres dispositifs comme des Li-DAR terrestres ou des drones. Cela serait aussi l'occasion d'incorporer les nouvelles métriques basées sur les données réelles que nous avons définies, permettant ainsi de créer une plateforme d'évaluation multi-scène, multi-données et multi-métriques.

Enfin, les algorithmes d'apprentissage machine connaissent une popularité grandissante même s'ils sont souvent limités à la reconstruction d'objets fermés simples. Dans ce contexte, ces travaux sur l'évaluation pourraient devenir un moyen de produire davantage de données d'entraînement et d'apprendre plus efficacement en utilisant les métriques présentées.

thesis proposes new approaches to tackle both these challenges.

Résumé de la thèse en français

Passage à l'échelle et évaluation de la reconstruction de surface à partir de nuages de points de scènes ouvertes