Julie Digne

David Cohen-Steiner

Pierre Alliez

Mathieu Desbrun

Fernando De Goes

Feature-Preserving Surface Reconstruction and Simplification from Defect-Laden Point Sets

Keywords: Optimal transportation, Wasserstein distance, Linear programming, Surface reconstruction, Shape simplification, Feature recovery

We propose a robust, feature-preserving surface reconstruction algorithm which turns a point set with noise and outliers into a low triangle-count simplicial complex. Our approach starts with a simplicial complex filtered from a 3D Delaunay triangulation of the input points. This initial approximation is iteratively simplified based on the optimal cost to transport the point set to the simplicial complex, both seen as measures (or mass distributions). Our optimal transport formulation allows the recovery of sharp features even in the presence of a large amount of outliers and/or noise in the input set.

Reconstruction et simplification de surfaces à partir de nuages de points imparfaits

Résumé : Nous proposons une méthode robuste de reconstruction de surface qui préserve les bords et les arêtes vives. Cette méthode part d'un nuage de points bruités et contenant des points aberrants pour reconstruire un complexe simplicial parcimonieux. Notre approche débute par la construction d'un complexe simplicial par filtrage d'une triangulation de Delaunay des points initiaux. Cette approximation initiale est ensuite itérativement simplifiée en se basant sur le coût de transport entre le nuage de points et le complexe simplicial, ceux-ci étant vus comme des distributions de masse. Cette formulation basée sur le transport optimal entre les deux distributions permet de retrouver les arêtes vives même en présence de nombreux points aberrants ou de bruit dans le nuage de points initial.

Mots-clés : Transport optimal, Distance de Wasserstein, Programmation linéaire, Reconstruction de surface, Simplification de forme, Recouvrement d'arête et de bord

Introduction

Surface reconstruction is a multi-faceted challenge which precise problem statement depends on the nature and defects of the input data, on the properties of the inferred surface (smooth vs piecewise smooth, with or without boundaries) and on the level of details sought after. Despite a number of major contributions on this topic over the past decade [START_REF] Amenta | The crust algorithm for 3d surface reconstruction[END_REF][START_REF] Kazhdan | Poisson surface reconstruction[END_REF], achieving both feature preservation and robustness to measurement noise and outliers remains a scientific challengeand a pressing requirement for many reverse engineering and geometric modeling applications. Furthermore, low polygon-count reconstructions has only received limited attention despite the increase of point density in 3D scanning technology and the need for efficient further geometry processing.

In this paper we contribute a reconstruction method that simplifies an initial (possibly nonmanifold) triangulation of the input point set based on an error metric that quantifies the distance between the current simplicial complex and the input points through optimal transport of measures. Our reconstruction approach inherits the qualities of the optimal transport derived metric, i.e., resilience to noise, outliers, and imperfect sampling, as well as sensitivity to boundaries and sharp features (creases, corners, tips, etc). We demonstrate these distinguishing properties on a series of examples. The proposed metric can also be used as a post-processing tool for recovering sharp features and boundaries from the output of a smooth reconstruction method.

Previous Work

We first discuss existing surface reconstruction methods, restricting our review to approaches that are both robust to noise and outliers and also feature preserving. We then point out recent, relevant work on geometry processing based on optimal transport.

Surface Reconstruction

A common approach to robust surface reconstruction from defect-laden point sets involves denoising and filtering out outliers, and often requires an interactive adjustment of parameters. Automatic methods such as spectral methods [START_REF] Kolluri | Spectral surface reconstruction from noisy point clouds[END_REF][START_REF] Walder | Implicit surface modelling as an eigenvalue problem[END_REF][START_REF] Alliez | Voronoi-based variational reconstruction of unoriented point sets[END_REF] and graph cut approaches [START_REF] Hornung | Robust reconstruction of watertight 3D models from nonuniformly sampled point clouds without normal information[END_REF][START_REF] Labatut | Robust and efficient surface reconstruction from range data[END_REF], on the other hand, are shown to be extremely robust, but are better suited to the reconstruction of smooth, closed surfaces. More recently, Cohen-Or and co-authors have proposed an excellent series of contributions based on robust norms and sparse recovery [START_REF] Lipman | Parameterization free projection for geometry reconstruction[END_REF][START_REF] Huang | Consolidation of unorganized point clouds for surface reconstruction[END_REF][START_REF] Avron | 1 -sparse reconstruction of sharp point set surfaces[END_REF]]. An interpolating, yet noise robust approach was alternatively proposed by Digne et al. [START_REF] Digne | Scale space meshing of raw data point sets[END_REF] through the construction of a scale space.

Feature preserving methods are typically based on an implicit representation that approximates or interpolates the input points. In [START_REF] Dinh | Reconstructing surfaces using anisotropic basis functions[END_REF], for instance, sharp features are captured through locally adapted anisotropic basis functions. Adamson and Alexa [START_REF] Adamson | Anisotropic point set surfaces[END_REF] proposed instead an anisotropic moving least squares (MLS) method using ellipsoidal mapping functions based on principal curvatures. More recently, Oztireli et al. [START_REF] Oztireli | Feature preserving point set surfaces based on non-linear kernel regression[END_REF] used a kernel regression technique to extend the MLS surface reconstruction, allowing for much sharper features. However, none of these techniques returns truly sharp features: reconstructions are always semi-sharp, with various degrees of smoothness depending on the approach and the sampling density. Moreover, sharpness is always determined locally, often leading to fragmented creases, so the reconstruction quality degrades quickly if defects and outliers are added.

Another way to detect local sharpness in a point set consists in performing a local clustering of estimated normals [START_REF] Ohtake | Multi-level partition of unity implicits[END_REF]: if this process reveals more than one cluster of normals, then the algorithm fits as many quadrics as the number of clusters. Improved robustness is achieved in [START_REF] Fleishman | Robust moving least-squares fitting with sharp features[END_REF] by segmenting neighborhoods through region growing. Lipman et al. [START_REF] Lipman | Data-dependent MLS for faithful surface approximation[END_REF], instead, proposed a systematic enrichment of the MLS projection framework with sharp edges driven by the local error of the MLS approximation. Again, the locality of the feature detection can generate fragmented sharp edges-much like general feature detection approaches (e.g., [START_REF] Gumhold | Feature extraction from point clouds[END_REF][START_REF] Pang | An algorithm for extracting geometric features from point cloud[END_REF]).

To reduce fragmentation a different thread of work aims at extracting long sharp features. Pauly et al. [START_REF] Pauly | Multi-scale feature extraction on point-sampled surfaces[END_REF], for instance, used a multi-scale approach to detect feature points, and construct a minimum-spanning tree to recover the most likely feature graph. Daniels et al. [START_REF] Daniels | Robust smooth feature extraction from point clouds[END_REF] used a robust projection operator onto sharp creases, and grew a set of polylines through projected points. Jenke et al. [START_REF] Jenke | Patch-graph reconstruction for piecewise smooth surfaces. Vision, modeling, and visualization[END_REF] extracted the feature lines by robustly fitting local surface patches and by computing the intersection of close patches with dissimilar normals.

[38] detects features on the point set and uses a Delaunay refinement to get a feature preserving model. This method can also be used on a model built through a noise-robust method such as Poisson [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] for enhancing sharp features of the model. We will show that the metric introduced for reconstruction can also be used in this setting.

Distinctive Approach. In this paper, we adopt a very different methodology: we reconstruct a shape through an iterative, feature-preserving simplification of a simplicial complex constructed from the input point set. To achieve unmatched noise and outlier robustness, the error metric driving the simplification is derived in terms of optimal transport between the input point set and the reconstructed mesh. Next we provide a brief review of optimal transport, and mention its applications to various problems in computer graphics and computer vision.

Optimal Transport

The problem of transporting a measure onto another one as a way to quantify their similarity has a rich scientific history [START_REF] Villani | Topics in Optimal Transportation[END_REF]. For two measures µ and ν defined over R 3 and of equal total mass (i.e., their integrals are the same), the optimal transport from µ to ν consists in finding a transport plan π that realizes the infimum:

inf Ω 2
x -y 2 dπ(x, y) π ∈ Π(µ, ν) , where Π(µ, ν) is the set of all possible transport plans between µ and ν, and L 2 norm is used as the transport cost.

The formulation above is particularly well suited to compare 1D measures such as histograms over the real line or on the circle [START_REF] Rabin | Transportation distances on the circle[END_REF], and it has been used for transferring color and contrast between images [START_REF] Reinhard | Color transfer between images[END_REF][START_REF] Rabin | Regularization of transportation maps for color and contrast transfer[END_REF]. For applications in higher dimensions such as 2D and 3D shape retrieval [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF][START_REF] Rabin | Geodesic shape retrieval via optimal mass transport[END_REF] and segmentation [START_REF] Peyré | Wasserstein active contours[END_REF], the optimal transport formulation is notoriously less tractable. Attempts have been made to design computationally easier surrogates, such as the sliced Wasserstein approach [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF], which consists in projecting the transport problem on a series of 1D problems.

An alternative and more general strategy to solve the optimal transport problem is based on linear programming. Such numerical tool has been used in applications such as surface comparison [START_REF] Lipman | Surface comparison with mass transportation[END_REF] and displacement interpolation [START_REF] Bonneel | Displacement interpolation using lagrangian mass transport[END_REF].

Distinctive Approach. Our reconstruction method also relies on a linear programming scheme. However, our approach introduces a key distinctive property: our target measure ν is not given, but instead, solved for. More specifically, we search for the simplicial complex of a user-specified size that minimizes the cost of transporting the input pointwise measure (i.e., the initial point set) to the complex simplices.

The closest work to ours was proposed by de Goes et al. [START_REF] De Goes | An optimal transport approach to robust reconstruction and simplification of 2d shapes[END_REF]. Their algorithm reconstructs and simplifies 2D shapes from point sets based also on optimal transport. Nonetheless, their [START_REF] De Goes | An optimal transport approach to robust reconstruction and simplification of 2d shapes[END_REF]. We depict the input point set and simplified triangulation. Red and green line segments depict the transport plan between the input point set and the uniform measure on respectively the vertices and edges of the triangulation. Each input point is simply transported either to its closest edge or to the end vertices of that edge. For each edge the transport plan favors a transport to its end vertices instead of to the whole edge when the corresponding transport cost is lower. The closeup image reveals that near the corner vertex transporting the nearby points to the vertex instead of the nearby edges creates a lower density near the end of that edge, and hence induces a tangential transport cost (red line segments pointing down) and consequently a higher total transport cost. For complex features such as concave corners on surfaces such behavior leads to incorrect reconstructions. Bottom: our transport plan is globally normal to the edges. approach differs from ours in several aspects: 1. their optimal transport involves only points and edges and therefore can be computed in closed form. To our knowledge, no such closed form exists when transporting points to the facets of a simplicial complex. Therefore we use a discretized formulation of the optimal transport problem.

Inria

2. the authors proposed to approximate the optimal transport plan by assigning each input point to its closest edge in the triangulation. Such simplistic scheme can lead to sub-optimal transport plan and cost as illustrated in Figure 1(top); even more so in 3D. Our discretized formulation, combined with a linear programming solver, provides better approximations of both the optimal plan and the optimal cost.

3. their method requires a valid embedding of a 2D triangulation, which they achieved through a recursive edge flip procedure. Such edge flip procedure can not, however, be generalized to 3D triangulations. Alternatively, our method removes the embedding requirement by only employing a (possibly non-manifold) simplicial complex, initially chosen as a subset of a 3D Delaunay triangulation.

Overview

Motivated by the concept of reconstruction introduced in 2D by de Goes et al. [START_REF] De Goes | An optimal transport approach to robust reconstruction and simplification of 2d shapes[END_REF], our 3D method consists in a fine-to-coarse algorithm which reconstructs a surface from a point set through greedy simplification of a 3D simplicial complex. We initialize the complex with a (possibly non-manifold) subset of the 3D Delaunay triangulation of input points, then perform repeated decimations based on half-edge collapse operations. The error metric guiding our simplification is derived from the optimal cost to transport the input point set (seen as Dirac measures) to a constant-per-facet measure defined over the simplicial complex. At each iteration, we collapse the half-edge which minimizes the increase of the total transport cost between the input points and the reconstructed triangulation. Similarly to the formulation presented in [START_REF] De Goes | An optimal transport approach to robust reconstruction and simplification of 2d shapes[END_REF], our optimal transport driven metric brings rare, but desirable properties, such as resilience to noise and outliers, and preservation of boundaries and sharp features.

In the remainder of this paper, we discuss the details of our optimal transport based metric (Sec. 3) and then describe our reconstruction algorithm step by step (Sec. 4). Our method is summarized in Algorithm 1 and its main stages are illustrated in Figure 2.

Optimal Transport as a Linear Program

Our reconstruction method consists in finding a coarse surface simplicial complex C that minimizes the optimal transport cost between the input point set S = {p i } i=1•••N (with associated measures {m i }) and a piecewise uniform measure on the facets of C. In [START_REF] De Goes | An optimal transport approach to robust reconstruction and simplification of 2d shapes[END_REF], the optimal transport cost between S and C is efficiently computed based on closed form expressions for the optimal cost between points and edges. However, to our knowledge, such closed form can not be extended between points and triangles. We thus adopt a discretized formulation of the optimal transport between S and C.

Discretization

We approximate the optimal transport cost between the input point set S and the simplicial complex C using quadrature points, which we call hereafter bins. To this end, we split vertices, edges, and facets of C into a set of bins B and evaluate the optimal cost between S and B as the sum of squared distances between the points in S and the centroid of the bins in B. We discretize the simplices of C as follows: a vertex is considered as its own (single) bin; an edge segment is decomposed into an integer number of equally-spaced bins; and a triangular facet is tiled by Inria a 2D Centroidal Voronoi Tessellation (CVT) with bins centered at the CVT sites (Figure 3). Note that CVT meshes minimize the approximation error given by quadrature points [START_REF] Du | Centroidal voronoi tessellations: Applications and algorithms[END_REF] and therefore it is the optimal point-wise discretization for the simplices in C. The number of bins per edge and facet is set based on a user-defined quadrature parameter (in our experiments, we used 200 bin per unit area). In order to accommodate the nearly uniform distribution of bins, we assign a capacity for each bin in B (i.e., the total amount of mass that a bin can receive).

While we set vertex bins to unit capacity, we assign edge bins proportional to the edge length and facet bins proportional to the area of the Voronoi cells.

Linear Programming Formulation

We now present a linear programming formulation to compute the optimal transport cost between the input point set S and the bin set B. In the following, we denote the simplices of C as , where L and M are the number of simplices and bins respectively. We also define s(j) to be the index of the simplex containing the bin b j (i.e., b j ∈ σ s(j) ). The capacity c j of bin b j is defined as the ratio between the bin's area and the area of its containing simplex. Finally, we denote m ij as the amount of mass transported from point p i ∈ S to the centroid of bin b j (Figure 4). With these definitions, we can now formally refer to a transport plan between S and B as a set of N × M variables m ij such that:

p i m ij4 m ij1 m ij2 m ij3 b j1 b j2 b j4 b j3
∀ij : m ij ≥ 0, (1)
∀i :

j m ij = m i , (2)
∀j 1 , j 2 s.t. s(j 1 ) = s(j 2 ) : 1 c j1 i m ij1 = 1 c j2 i m ij2 , (3)

Inria

where Equation 2 ensures that the entire measure of a point gets transported onto simplices, and Equation 3 ensures a uniform measure over each simplex of C.

An optimal transport plan is then defined as a transport plan π that minimizes the associated transport cost

cost(π) = ij m ij p i -b j 2 .
Finding a transport plan minimizing the transport cost results in a linear program with respect to the m ij , with equality (Eq. 2 and 3) and inequality constraints (Eq. 1). Note that the number of bins, the position of bins, as well as the square distances between input points and bin centroids are all precomputed. In order to enforce the uniformity constraint (Eq. 3), we also introduce L additional variables l i (one per simplex σ i ) indicating the target measure density of the corresponding simplex. The final formulation is thus:

Minimize ij m ij p i -b j 2
w.r.t. the variables m ij and l j , and subject to:

             ∀i : j m ij = m i ∀j : i m ij = c j • l s(j)
∀i, j : m ij ≥ 0, l j ≥ 0

Local Relaxation

The formulation described above is compute-intensive due to the number of variables and constraints involved. Solving this linear programming system requires instantiating a dense matrix (representing the constraints) of size (M × N + L) × (M + N ). For example, computing the optimal transport cost between a simplicial complex containing 782 simplices on which 7, 300 bins are sampled and an input point set of 2, 100 samples involves solving for a linear program of over 15 million variables and 9, 000 constraints. The linear programming solvers at our disposal do not scale to such large numbers.

In order to improve scalability we propose a stencil-based local relaxation strategy (summarized in Algorithm 2) for approximating the global optimization problem. A subset of the global solution space is explored through local solves, until a local minimum of the objective function is reached. Note that we cannot guarantee convergence of this local procedure to the global minimum, but the minima reached in practice have been sufficient in all of our tests to obtain satisfactory results.

Our procedure starts with a trivial transport plan which transports each input point to its nearest vertex of the simplicial complex C. Since no uniformity constraints are imposed on vertices, this transport plan is valid (but obviously suboptimal). Subsequent local optimizations can only decrease the global transport cost or leave it unchanged, as the rest of the plan is kept intact and the local re-assignments have smaller or equal cost after optimization. The transport cost found by local stencil updates is thus an upper bound of the global optimal transport cost. Our experiments showed, unsurprisingly, that convergence depends on the shape of local stencils used: the larger the stencil, the faster the convergence-but with the unfortunate side effect that large stencils increase the size of the corresponding linear program. We found in practice that simply using the 1-ring of a chosen simplex is always a reliable choice. More precisely, the local stencil N of a facet σ includes all the facets incident to σ along with their vertices. Algorithm 2: Solving for a local stencil.

Input : A simplicial complex C and point set S, a threshold ε Output: A locally optimal transport plan π = {m ij } for p i ∈ S do Transport p i to nearest vertex v ∈ C;

δ = ε + 1; while |δ| > ε do for σ j ∈ C do
Build the stencil N of the facet σ j ; Collect sample points and partial measures {p i , mi } that project onto this stencil; Solve the linear program for the optimal transport plan of (p i , mi ) on the bins of N ; Update transport plan π and global cost;

δ = new_cost -old_cost;
Armed with this scalable approximation of the transport cost, we describe next how we put it at work for surface reconstruction through simplification.

Reconstruction through Simplification

Initialization

We begin our reconstruction process by randomly picking a subset of the input points S and computing a 3D Delaunay triangulation. We then construct a simplicial complex C from a subset of facets of this 3D triangulation. To select this subset of facets, we perform two steps: (1) we reuse the local stencil relaxation method (Algorithm 2) to estimate a transport cost from all the input points onto the facets and vertices of the 3D triangulation; (2) we then build C with only the facets containing non-zero transported measure. For the step (1), we use a stencil centered at each facet and containing vertices and edges of the two tetrahedra adjacent to the facet. For an inside facet of the triangulation, for instance, this stencil contains 7 facets and 5 vertices. This stencil-based optimization is repeated in batches until the decrease in transport cost is below a user-specified threshold (set to 10 -5 in all our tests). In practice, the global transport cost decreases rapidly, and only 10 batches are usually needed. For step (2), we convert our data structure to a simplicial complex for two main reasons: first to allow our reconstruction to have long and anisotropic simplices; and more importantly, to remove the difficult issue identified in [START_REF] De Goes | An optimal transport approach to robust reconstruction and simplification of 2d shapes[END_REF] of keeping the triangulation with a valid embedding.

Decimation

From the initial simplicial complex C, we further simplify the reconstruction through a greedy decimation based on half-edge collapse operations. Note, however, that a conventional decimation algorithm (e.g., [START_REF] Garland | Surface simplification using quadric error metrics[END_REF][START_REF] Lindstrom | Evaluation of memoryless simplification[END_REF] and variants) can not be applied in our setup: the presence of outliers and noise renders the typical error metrics inadequate.

Our optimal transport framework provides a solid alternative: we pick the next half-edge to collapse as the one that increases the least the global transport cost. To this end, we simulate the collapse of a candidate half-edge e and evaluate the induced change of transport cost ∆. Since Inria this change cost affects mainly a neighborhood Ω e of e, we can restrict the computation of ∆ only to Ω e . More specifically, setting Ω e to the closure of simplices in the 1-ring of e, we first gather the set of samples p i projecting (partially or entirely) on Ω e (Figure 5), compute the transport cost of this set of samples to Ω e , simulate the collapse of e, and recompute the cost of transporting the set of samples onto the resulting simplices. The change of transport cost ∆ is then set as the difference of transport cost achieved before and after the simulated collapse of e. Once a half-edge is collapsed, we also update the transport plan of the edges with 1-ring intersecting the collapsed edge. Finally, to increase scalability, we employ a multiple choice approach [START_REF] Wu | Fast mesh decimation by multiple-choice techniques[END_REF] to select the next half-edge to be collapsed (as recommended in [START_REF] De Goes | An optimal transport approach to robust reconstruction and simplification of 2d shapes[END_REF]), instead of maintaining an compute-intensive global priority queue. The edge selection process is summarized in Algorithm 3.

Figure 5: Local stencil of an edge. For better depiction we represent a manifold neighborhood of an edge (in red) and do not depict the bins. Simplices in the local stencil are depicted in blue. Point samples are depicted in green. We only solve for the measures transported from and to the stencil (solid green lines) and not for the measures transported outside of the stencil (dash green lines). Note that different transport plans are possible between point samples and bins, i.e., from inside to inside the stencil, inside to outside, and outside to inside.

Algorithm 3: Decimation algorithm.

Input : A simplicial complex C, input point set S, the number of desired vertices V Output: A simplicial complex C f inal with V vertices for each edge e ∈ C do Simulate two half-edge collapse operators; Push the corresponding half-edges to a priority queue P, sorted according to the change of transport cost ∆. repeat Pop half-edge e * out of P;

Collect set E of edges whose neighborhood intersect neighborhood of e * ; Collapse e * and update transport plan on the neighborhood of e * ; Update P by recomputing the change of cost ∆ for all edges in E. until The simplicial complex has V vertices;

Vertex Relocation

So far our method based half-edge collapses limits our reconstruction to be an interpolation of a subset of the input points. This may lead to unsatisfactory results specially in the presence of noise and outliers. We thus couple our decimation with an optimization procedure in order to relocate the position of vertices in the reconstructed simplicial complex C. After the collapse of a half-edge e, the remaining vertex v of e is relocated by iterating two steps: (1) for a given transport plan π, we move v toward the position that improves the optimal cost of π; (2) we update π around v. For the first step, we compute the locally optimal position of v when the transport plan π (i.e., m ij ) is kept fixed. To this end, we express the position of each bin (be it a vertex bin, an edge bin, or a triangle bin) in barycentric coordinates within its containing simplex. For example, finding the optimal position of vertex v of triangle t = (v, v 1 , v 2 ) amounts to minimizing:

min v i j m ij p i -α j v -β j v 1 -γ j v 2 2 ,
where α j , β j , γ j are the barycentric coordinates of bin b j with respect to vertices (v, v 1 , v 2 ). The optimal position with respect to triangle t is thus

v s = i j m ij α j (p i -β j v 1 -γ j v 2 ) i m ij α 2 j .
For the second step, we freeze the vertex locations and update the transport map π by resolving locally the linear program associated to the local transport (Algorithm 2). By alternating these two steps, the vertices move to their optimal position, allowing for the recovery of sharp features and surface boundaries. Figure 6 depicts a simple vertex relocation sequence in 2D.

Facet Filtering

After the decimation terminates, we return as our final reconstructed mesh the subset of facets from C with a non-zero measure. Since some facets may have non-zero measure due the presence of noise and outlier, we instead sort the facets based on the measure density (i.e., the ratio of facet measure and its area) and provide the user with an interactive slider to recover the most dense facets. Figure 7 shows the reconstructed surface obtained with different filtering threshold.

Experimental Results

We implemented our algorithm in C++ using CGAL's 3D Delaunay triangulation library [6] to initialize the reconstruction, and our own data structure for simplicial complexes. We used the Coin-OR Clp library [START_REF]CLP, coin-or linear program solver[END_REF] as our linear programming solver. Our implementation is partially parallelized to accelerate computations, exploiting the fact that all half-edge collapse simulations are independent. The initialization and update of the priority queue are, by far, the most costly operations, as each collapse involves around 120 simulations on average. When using the exhaustive priority queue on a two-core laptop, a point set containing 30, 000 points is reconstructed in around 10 hours (initial and final simplicial complexes containing respectively 3, 000 and 200 vertices). On a 8-core processor, this computation reduces to 2 hours. However, when using a multiple-choice approach with random sets of 40 collapses, the timings are three times faster on average. Memory consumption remains low; e.g., for the experiment mentioned above, the peak memory usage was around 80M o.

Robustness to noise. We tested our method on a point set of a staircase with an increasing amount of uniform synthetic noise (Figure 8). Even in the presence of a large noise the method tends to recover well the creases of the input shape. Only for noise magnitudes larger than 5% of the bounding box size does our method fail (for these high noise levels, spurious facets appear which cannot be discarded by a simple thresholding based on mass density). Our method can thus robustly handle inputs from any current point acquisition devices, which all have much smaller noise magnitudes. For visual clarity we choose a 2D example with a single triangle and only facet bins. We first depict the input point set, here uniformly sampled on a triangle, the initial simplicial complex composed of one facet, and the facet bins and their capacities. For all subsequent images we depict during vertex relocation the transport plan with blue edges connecting the source point samples and their target bin centroids.

Inria

Robustness to outliers. We also tested our method on a point set that samples a cylinder (Figure 9). Results are excellent up to 10% of outliers, but our method can fail when the amount of outliers exceeds 20%-again, current acquisition devices are usually good enough not to reach this amount of outliers.

Feature preservation. Figure 10 depicts the feature preservation property of our approach on the blade model. Our approach performs well even on thin features subtending small angles, where implicit approaches (here, the noise-robust Poisson surface reconstruction method [START_REF] Kazhdan | Poisson surface reconstruction[END_REF]) tend to smooth out features and create spurious topological artifacts in presence of insufficient density in the input point set. On easier features subtending large angles, both approaches produce similar results when comparing our approach to the combination of a smooth implicit reconstruction followed by surface mesh decimation [START_REF] Lindstrom | Evaluation of memoryless simplification[END_REF].

On a cone model (Figure 11) all features (tip, boundaries) are preserved and the simplification is very effective. Similarly, on a cylinder model (Figure 12) the boundaries are preserved and the simplification leads to anisotropic triangles with most edges aligned to minimum curvature directions as expected.

Figure 13 also illustrates boundary and sharp feature preservation. Figure 14 illustrates the behavior of our approach on two intersecting planar polygons. The algorithm behaves well down Figure 7: Facet filtering. We first depict the input point set and all facets with nonzero measure of the simplicial complex returned by our decimation scheme. We then filter out the spurious facets by thresholding facets with low measure density. The final reconstruction is highlighted in a frame. to 10 vertices, and the simplicial complex maintains the initial connectivity during decimation. Going down to 8 vertices (the expected minimum number of vertices) would require a richer set of topological operators in order to disconnect the intersecting edge before pursuing decimation; we did not pursue this particular property.

Weaknesses. Given the current linear programming solvers, results of our approach come at the price of intensive computations, preventing its use on large point sets. Also, there is currently nothing in our formulation that favors 2-manifoldness, as the main data structure is a simplicial complex initialized by the facets of a 3D triangulation; this can lead, in rare occasions, to invalid embedding as well as multiple facets covering the same area (see Figure 15). The latter issue is more complex than just ensuring a 2-manifold reconstruction, as complex features may correspond to non-manifold shapes. One could define a notion of "effectiveness" per facet favoring a high transported mass and a low transport cost, but this notion would lead to a non-linear objective function and would require a richer set of topological operators such as facet deletion.

Feature Recovery

Another application of the proposed metric is to recover sharp features and boundaries from the output of a surface reconstruction specialized to smooth, closed surfaces (e.g., Poisson reconstruction [START_REF] Kazhdan | Poisson surface reconstruction[END_REF]). These approaches are in general noise robust and scalable, but smoothen sharp features and fill holes. For this application we use the vertex relocation and facet filtering steps

Conclusion

We introduced a surface reconstruction and simplification method which exhibits both robustness to noise and outliers, as well as preservation of sharp features and boundaries. Our approach is based on the decimation of a simplicial complex guided by an optimal transportation error metric between the reconstruction and the initial point set. The metric proves also successful for feature recovery after a reconstruction method producing smooth shapes.

The main drawback of our approach is its computational cost: despite our efforts to introduce local relaxation, parallelization, and multiple-choice accelerations, we cannot reconstruct large point sets in reasonable time. The main strength of our approach lies into the simplicity of its formulation, which is expressed directly on the simplicial complex being reconstructed. This departs from common robust operators which require subsequent contouring to obtain the final reconstructed (but not simplified) surface mesh. In addition, our formulation allows the reconstruction of shapes of various dimensions as, e.g., a curve embedded in R 3 would be automatically reconstructed as a set of edges with a large majority of edge bins. Furthermore our formulation provides us with a transport plan, which can be used for further geometry processing of the resulting simplicial complex.

As future work we wish to improve scalability. Reusing the multi-scale approach of Mérigot [START_REF] Mérigot | A multiscale approach to optimal transport[END_REF] is certainly a valid direction but requires significant work to be fully practical.

Figure 1 :

 1 Figure 1: Transport plans. Top: Binary transport plan[START_REF] De Goes | An optimal transport approach to robust reconstruction and simplification of 2d shapes[END_REF]. We depict the input point set and simplified triangulation. Red and green line segments depict the transport plan between the input point set and the uniform measure on respectively the vertices and edges of the triangulation. Each input point is simply transported either to its closest edge or to the end vertices of that edge. For each edge the transport plan favors a transport to its end vertices instead of to the whole edge when the corresponding transport cost is lower. The closeup image reveals that near the corner vertex transporting the nearby points to the vertex instead of the nearby edges creates a lower density near the end of that edge, and hence induces a tangential transport cost (red line segments pointing down) and consequently a higher total transport cost. For complex features such as concave corners on surfaces such behavior leads to incorrect reconstructions. Bottom: our transport plan is globally normal to the edges.

Algorithm 1 :

 1 Algorithm Overview. Input : Point set S, user-specified number of vertices. Output: Simplicial complex C. Construct 3D Delaunay Triangulation T from S; Compute transport cost from S to the facets of T ; Construct simplicial complex C from all facets of T with non-zero measure; Decimate C until desired number of vertices; Filter out facets of C by thresholding mass density.

Figure 2 :

 2 Figure 2: Steps of our algorithm: (a) Initial point set; (b) 3D Delaunay triangulation of a random subset containing 10% of the input points; (c) Initial simplicial complex constructed from facets of the 3D triangulation with non-zero measure; (d) Initial transport plan assigning point samples to bin centroid (green arrows); (e-f) decimation steps; (g-i) Reconstruction with 100, 50, and 14 vertices, respectively; (j-l) Final transport plan with 100, 50, and 14 vertices, respectively.

Figure 3 :

 3 Figure 3: Bins. Left: bins in a facet. Vertices have a single bin each. Edges have equally-spaced bins. Facets have evenly placed bins constructed through a Centroidal Voronoi Tessellation. Inner Voronoi cells are depicted in green. Voronoi cells clipped by the facet boundary are depicted in pink. Bin centroids are depicted as red dots. Right: capacities of the bins are depicted using a thermal color ramp.

Figure 4 :

 4 Figure 4: Transport plan for a single input sample point p i . Variable m ij models the transport of p i to the j th bin of the facet.

Figure 6 :

 6 Figure 6: Vertex relocation.For visual clarity we choose a 2D example with a single triangle and only facet bins. We first depict the input point set, here uniformly sampled on a triangle, the initial simplicial complex composed of one facet, and the facet bins and their capacities. For all subsequent images we depict during vertex relocation the transport plan with blue edges connecting the source point samples and their target bin centroids.

Figure 8 :

 8 Figure 8: Robustness to noise. We increase the amount of synthetic noise from σ = 1% to σ = 2% through σ = 5%, expressed in percentage of the longest edge length of the bounding box. The reconstruction starts failing at σ = 5%.

Inria

Figure 9 :

 9 Figure 9: Robustness to outliers. The reconstruction is effective with 10% outliers (left) but fails from 20%. The outliers are added randomly within a loose bounding box (120%) of the input point set.

Figure 10 :

 10 Figure 10: Reconstruction of the blade model containing 30K sample points. Top: our reconstruction. Bottom: the output of the Poisson reconstruction method (Delaunay refinement being used for contouring the implicit function), and closeup on a sharp crease subtending a small angle, where the implicit approach fails.

Figure 11 :

 11 Figure 11: Reconstruction of a cone. Left: input point set Middle: input point set and final simplicial complex with (nearly uniform) facet densities shown. Right: final simplicial complex.

Figure 12 :

 12 Figure 12: Reconstruction and simplification of a cylinder. Left: 10K noisy sample points and reconstruction with 12 vertices (facet density shown). Middle: transport plan between point samples and bin centroids. Right: simplicial complex and facet density.

Figure 13 :

 13 Figure 13: Reconstruction of a twisted bar. Sharp features are well preserved.

Inria

Figure 16 :

 16 Figure 16: Anchor. Noisy point set (top), Poisson reconstruction (middle left), improved reconstruction (middle right) and associated closeups.

Figure 17 :

 17 Figure 17: Blade. Poisson reconstruction (left) and improved reconstruction (right). The input point set is depicted with black dots. No remeshing nor edge flips are applied: the spurious topological handles shown in Figure 10 are not repaired as the triangles are only pulled toward the point set.

Figure 18 :

 18 Figure 18: Church. Point set (left), Poisson reconstruction (middle) and relocated mesh (right).

Figure 19 :

 19 Figure 19: Church. Point set, Poisson reconstruction and vertex relocation. Filtering combined with vertex relocation allows recovering he surface boundaries.

Figure 20 :

 20 Figure 20: Sharp sphere. Left column: global view, right column: close-up. From top to bottom: Point set, smooth reconstruction and vertex relocation. Features are recovered through vertex relocation.

RR n°7991

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

The input to the algorithm is a surface triangle mesh (the output of a smooth reconstruction algorithm) and the original point set used for reconstruction. Bins are initially sampled on the mesh. The initial assignment is performed through relaxation as described in Section 3.3: each sample is assigned to the nearest mesh vertex, and local reassignments are iterated until a local minimum for the transport cost is reached. Each mesh vertex is then relocated as described in Section 4.3. Through this method sharpness is recovered (Fig. 16).

For open surfaces this method allows for recovering the boundaries of the surface through the last filtering step (section 4.4) as can be seen on the church example (Figure 18 and19). On the latter, the relocation seems incorrect at first glance on the bell tower, but after a more careful examination of the point set the spurious triangles created by the relocation procedure


Wolfgang Bertram
email: wolfgang.bertram@univ-lorraine.fr

Michael Kinyon
email: mkinyon@du.edu

TORSORS AND TERNARY MOUFANG LOOPS ARISING IN PROJECTIVE GEOMETRY

Keywords: 2010 Mathematics Subject Classification. 06C05, 20N05, 20N10, 51A35 projective plane, Moufang plane, (ternary) Moufang loop, torsor (heap, groud, principal homogeneous space), lattice

We give an interpretation of the construction of torsors from [BeKi10a] in terms of classical projective geometry. For the Desarguesian case, this leads to a reformulation of certain results from [BeKi10a], whereas for the Moufang case the result is new. But even in the Desarguesian case it sheds new light on the relation between the lattice structure and the algebraic structures of a projective space.

1. The geometric construction 1.1. The generic case. In this first chapter, we describe the general construction of torsors and of ternary loops associated to projective spaces; proofs and computational descriptions are given in the two following chapters. We assume that X is a projective space of dimension at least two. For projective subspaces a, b of X , let as usual a ∧ b be the meet (intersection) and a ∨ b be the join (smallest subspace containing a and b).

Definition 1.1. Consider a pair (a, b) of hyperplanes in X and a triple of points (x, y, z), none of them in a or b. Assume that x, y, z are not collinear. Then we define a fourth point w := xz := (xyz) ab as follows: w is the intersection of -the parallel of the line x ∨ y through z in the affine space V a := X \ a, with -the parallel of the line z ∨ y through x in the affine space V b = X \ b; that is:

w = xz = (xyz) ab = (x ∨ y) ∧ a ∨ z ∧ (z ∨ y) ∧ b ∨ x .
Note that this point of intersection exists since all lines belong to the projective plane spanned by x, y, z. For a = b, this is the usual "parallelogram definition" of vector addition in the affine space V a with origin y, that is, xz = x + z in this case. Hence, for a = b, (xyz) ab may be seen as a kind of "deformation of vector addition": we have a sort of "fake parallelogram" with vertices y, x, z, w, as shown in the illustrations below. As for "usual" parallelograms, it is easily seen that, with (xyz) := (xyz) ab for fixed (a, b), the conditions These images admit a spacial interpretation: we may imagine the observer placed in affine space R 3 inside a plane B which is vizualized only by its "horizon", the line b; then we think of the line y ∨ z as lying in a plane B ′ parallel to B, and of the line x ∨ w as lying in another such plane B ′′ ; the other two lines w ∨ z and x ∨ y lie in planes that are parallel to the drawing plane P . This interpretation is not symmetric in x and z: the point z lies "behind" (or "in front of") y, whereas x is considered to be "on the same level" as y.

The product xz is thus in general not commutative, but it is associative: we show that, if X is Desarguesian, then, for any fixed origin y, the binary map (x, z) → xz gives rise to a group law on the intersection of affine parts

(1.3) U ab := X \ (a ∪ b) = V a ∩ V b .
More generally and more conceptually, we show that the ternary law (x, y, z) → (xyz) ab defines a torsor structure on U ab (Theorem 2.1). Naturally, the question arises what we can say for general, non-Desarguesian projective planes, or for still more general lattices.

The most prominent class of non-Desarguesian projective planes are the Moufang planes: we show that in this case we get a kind of "alternative version of a torsor" which we call a ternary Moufang loop (Theorem 3.1). For a = b, these ternary Moufang loops contract to the abelian vector group of an affine plane. For very general projective planes (which need not be "translation planes") it remains an interesting open problem to relate this new algebraic structure to those traditionally considered in the literature: indeed, our definition is closely related to the more traditional ways of coordinatizing projective planes by ternary rings. This is related to the following item. (2) If x = y, then let us choose a point u not belonging to x ∨ y or to a, and we let

w := (xyz) ab := (x ∨ y) ∧ (z ∨ u) ∧ b ∨ (x ∨ y) ∧ a ∨ u ∧ (u ∨ y) ∧ b ∨ x .
(It will be shown below that w does not depend on the choice of u.) (3) If z = y, then we let w := (xyz) ab := (zyx) ba , where the right hand side is defined by the preceding case.

This definition can be interpreted from two different viewpoints:

(A) Algebraic. In the Desarguesian case, the expression in (2) is derived from our first definition by using para-associativity and idempotency:

((xyu)uz) = (xy(uuz)) = (xyz) ,
where now the left hand side can be expressed by using twice Definition 1.1 (see Theorem 2.3 below). This is indeed in keeping with idea explained above of "taking a limit" (imagine u tending towards a point on the line L). The argument still goes through in the Moufang case since one does not need for it full para-associativity, but just a special case which remains valid precisely in the Moufang case (but it breaks down as soon as one wants to go further).

(B) Geometric. The formula in (2) corresponds to classical "constructions of the field associated to a plane". It is known that in the Moufang case the field does not depend on the "off-line" point u. More specifically, we distinguish two cases in (2):

• generic case: if the points L ∧ a and L ∧ b are different, then (xyz) ab is the product zy -1 x on the vector line L with "point at infinity" L ∧ a and "zero point" L ∧ b (in the following illustration, a is the line at infinity and b the horizontal line; in usual textbook drawings, the inverse choice is made. We have marked the points p = (xyu) and q = (puz) = w.)

y x z u p q
• special case: if the line L := x ∨ y intersects a ∧ b, then (xyz) ab is the "ternary sum"

x -y + z in the affine line L (with L ∧ (a ∧ b) as "point at infinity", see the following illustration, which is the limit case of the preceding one, as L becomes parallell to b). The main result of the present work can now be stated as follows:

Definition 1.3. A set G with a map G 3 → G, (x, y, z) → (xyz) is called a torsor if (xxy) = y = (yxx) (T0) (xy(zuv)) = (x(uzy)v) = ((xyz)uv) (T1)
and a ternary Moufang loop if it satisfies (T0) and

(uv(xyx)) = ((uvx)yx) (MT1) (xy(xyz)) = ((xyx)yz) (MT2)
Theorem 1.4. If X is a projective space of dimension bigger than one over a skew-field (i.e., a Desarguesian space), then the preceding constructions define a torsor law on U ab . If X is a Moufang projective plane, then the constructions define a ternary Moufang loop.

1.3. Generalized cross-ratios, and associative geometries. In the Desarguesian case, a very general theory describing torsors of the kind of U ab has been developed in [START_REF] Bertram | Associative Geometries. I: Torsors, Linear Relations and Grassmannians[END_REF].

Comparing with the approach presented here, one may ask for what kinds of lattices there are similar theories -we will, in subsequent work, investigate in more depth the case of Moufang spaces, related to alternative algebras, triple systems and pairs. Returning to the Desarguesian case and to classical projective geometry, the link between the lattice and the structure defined in [START_REF] Bertram | Associative Geometries. I: Torsors, Linear Relations and Grassmannians[END_REF] is surprisingly close; however, one should not forget that for projective lines the lattice structure is completely useless, whereas the structures from [START_REF] Bertram | Associative Geometries. I: Torsors, Linear Relations and Grassmannians[END_REF] are at least as strong as the classical cross-ratio, and hence are much stronger than the lattice structure. Let us briefly explain this. Given a unital ring K and X := X (Ω), the full Grassmannian geometry of some K-module Ω (set of all submodules of Ω), we have associated in [START_REF] Bertram | Associative Geometries. I: Torsors, Linear Relations and Grassmannians[END_REF] to any 5-tuple (x, a, y, b, z) ∈ X 5 another element of X by

(1.4) Γ(x, a, y, b, z) := ω ∈ Ω ∃ξ ∈ x, ∃α ∈ a, ∃η ∈ y, ∃β ∈ b, ∃ζ ∈ z : ω = ζ + α = α + η + β = ξ + β .
In [START_REF] Bertram | Associative Geometries. I: Torsors, Linear Relations and Grassmannians[END_REF], Theorem 2.4, it is shown that the lattice structure is recovered via

(1.5) x ∧ a = Γ(x, a, y, x, a), b ∨ a = Γ(a, a, y, b, b)
for any y ∈ X . On the other hand, in the present work we prove (Theorem 2.3) that, if K is a field, if a, b are hyperplanes and x, y, z one-dimensional subspaces, then Γ(x, a, y, b, z) can be recovered from the lattice structure via

(1.6) Γ(x, a, y, b, z) = (xyz) ab .
Thus, roughly speaking, for Desarguesian projective spaces of dimension bigger than one, Γ and the lattice structure are essentially equivalent data. Summing up, there are two major approaches to our object: the algebraic approach ([BeKi10a]), based on associative algebras and -pairs and on an underlying group structure of the "background" Ω (cf. [START_REF] Bertram | The projective geometry of a group[END_REF]), and the lattice theoretic approach from the present work, keeping close to classical geometric language, and paving the way to incorporate exceptional geometries into the picture.

The Desarguesian case

Theorem 2.1. Assume X is a Desarguesian projective space of dimension bigger than one, and fix a pair (a, b) of hyperplanes. Then U ab , together with the ternary product (xyz) := (xyz) ab defined above, is a torsor. In particular, if we fix an "origin" y ∈ U ab , then U ab with product xz = (xyz) ab and origin y becomes a group. If a = b, then this is group is not commutative.

Proof. We give three different proofs. In all cases, the main point is to prove "paraassociativity" (T1).

(a) The first proof is rather a "drawing exercise": let us show, just by using Desargues' Theorem, that ((xyz)uv) = (xy(zuv)). We construct first the point ((xyz)uv). This is best visualized by choosing a as line at infinity of our drawing plane, and we may draw the lines y∨x and z∨(xyz) as vertical lines. Then ((xyz)uv) is the point q in the following illustration:

B x y z (xyz) u v q (zuv)
Now, the triangles u, z, (xyz) and v, (zuv), q are in a Desargues configuration, and we conclude that the line q ∨ (zuv) is parallel to z ∨ (xyz), i.e., it is vertical. But then the triangles y, z, (zuv) and x, (xyz), q are also in Desargues configuration, i.e., the intersection points of corresponding sides lie on a common line, which must be b. It follows that (x ∨ q) ∧ b = (y ∨ (zuv)) ∧ b, from which the desired equality follows.

As a next challenge, one may try to establish the remaining equality defining paraassociativity in the same vein.

(b) A computational proof. Let K be the (skew)field of X , and work in the affine space V := V a . If a = b, then (as mentioned above), (xyz) aa = x -y + z is the torsor law of the affine space V a , and the claim is obviously true. If a = b, fix some arbitrary origin o in the affine hyperplane V a ∩ b. There is a linear form β :

V → K such that b ∩ V = ker(β), so that U ab = {x ∈ V a | β(x) ∈ K × }.
Lemma 2.2. For all x, y, z ∈ U ab , in the vector space (V a , o), we have

(xyz) ab = β(z)β(y) -1 (x -y) + z.
Proof. Assume first that x, y, z are not collinear. The parallel of x ∨ y through z in V a is

(x ∨ y) ∧ a ∨ z = {z + s(x -y) | s ∈ K}.
We determine the point (y∨z)∧b. If y∨z is parallel to b, then we get easily from the definition that (xyz) ab = x -y + z is the usual sum, which is in keeping with our claim. Assume that y ∨ z is not parallel to b. Then the intersection point (y ∨ z)∧ is obtaining by solving

β (1 -t)y + tz = 0, whence t = β(y)(β(y) -β(z)) -1 , whence 1 -t = -β(z)(β(y) -β(z)) -1 and (y ∨ z) ∧ b = -β(z)(β(y) -β(z)) -1 y + β(y)(β(y) -β(z)) -1 . The intersection of (z ∨ y) ∧ b ∨ x and (x ∨ y) ∧ a ∨ z = z + K(x -y) is determined by r, s ∈ K such that (1 -r)x + r -β(z)(β(y) -β(z)) -1 y + β(y)(β(y) -β(z)) -1 z = sx -sy + z.
Since both sides are barycentric combinations of x, y, z, we may consider y as new origin. Then, if x and z are linearly independent with respect to this origin, this condition is equivalent to

1 -r = s, r(β(y) β(y) -β(z) -1 ) = 1 whence r = β(y) -β(z) β(y) -1 and s = 1 -r = β(z)β(y) -1
, and finally

(xyz) ab = s(x -y) + z = β(z)β(y) -1 x -y + z ,
proving our claim in the non-collinear case. Now consider the collinear case. As pointed out after Definition 1.2, in this case the definition of (xyz) ab amounts to the geometric definition of the field operations. If the line L spanned by x, y, z is parallel to b, then β(z) = β(y), and the formula from the lemma gives the additive torsor law x -y + x, as required. Else, choose o := L ∧ b as origin, let u ∈ L with β(u) = 1 and write x = ξu, y = ηu, z = ζu with ξ, η, ζ ∈ K × , and then the formula from the lemma gives (xyz) ab = ζη -1 (ξu -ηu) + ζu = ζη -1 ξu, which again corresponds to the definition given in this case. Thus the claim holds in all cases.

Using the lemma, we now prove the torsor laws: first of all, we have (2.1)

β (xyz) = β β(z)β(y) -1 x -y + z = β(z)β(y) -1 β(x) ,
showing that U ab = V a \ker(β) is stable under the ternary law. The idempotent laws follow by an easy computation from the lemma. For para-associativity, using (2.1), a straightforward computation shows that both ((xyz)uv) and (x(uzy)v) are given by

β(v)β(u) -1 β(z)β(y) -1 (x -y) + β(v)β(u) -1 (u -z) + v.
(b') Remark: there is a slightly different version of (b), having the advantage that the cases a = b and a = b can be treated simultaneously, and the drawback that the dependence on y is not visible: choose o := y as origin in V = V a , and a linear form β :

V → K such that b ∩ V = {x ∈ V | β(x) = 1}. The case a = b then corresponds to β = 0. A computation similar as above yields (2.2) xz = (xyz) ab = (1 -β(z))x + z = x -β(z)x + z
from which associativity of the product xz follows easily. Note that Formula (2.2) is a special case of the formulae given in Section 1.4 of [START_REF] Bertram | Associative Geometries. I: Torsors, Linear Relations and Grassmannians[END_REF]. (c) A third proof follows from Theorem 2.3 in [START_REF] Bertram | Associative Geometries. I: Torsors, Linear Relations and Grassmannians[END_REF], combined with the following general result: Theorem 2.3. Let K be a unital ring and X = X (Ω) be the full Grassmannian geometry of some K-module Ω (set of all submodules of Ω), and define, for a 5-tuple (x, a, y, b, z) ∈ X 5 , the submodule Γ(x, a, y, b, z) by Equation (1.4).

(1) Assume that the triple (x, y, z) is in general position, that is,

x ∧ (y ∨ z) = 0, or y ∧ (x ∨ z) = 0, or z ∧ (x ∨ y) = 0.
Then we have the following equality of submodules of Ω:

Γ(x, a, y, b, z) = (x ∨ y) ∧ a ∨ z ∧ (z ∨ y) ∧ b ∨ x .
(2) Assume that z is contained in x ∨ y, i.e., z ∧ (x ∨ y) = z. Then, for any choice of u ∈ U ab satisfying u ∧ (x ∨ y) = 0, we have

Γ(x, a, y, b, z) = ([([ (x ∨ y) ∧ a ∨ u ∧ (u ∨ y) ∧ b ∨ x ] ∨ u) ∧ a] ∨ z) ∧ [((z ∨ u) ∧ b) ∨ (x ∨ y) ∧ a ∨ u ∧ (u ∨ y) ∧ b ∨ x ]
(3) Let a, b be hyperplanes in a vector space and x, y, z lines. Retain assumptions from the preceding item and assume that x = y. Then the expression given there simplifies to

Γ(x, a, y, b, z) = (x ∨ y) ∧ [((z ∨ u) ∧ b) ∨ (x ∨ y) ∧ a ∨ u ∧ (u ∨ y) ∧ b ∨ x ] .
Proof.

(1) We prove first the inclusion "⊂" (which holds in fact for all triples (x, y, z)): on the one hand, (

x ∨ y) ∧ a ∨ z ∧ (z ∨ y) ∧ b ∨ x is the set of all ω ∈ Ω such that we can write ω = α + ζ, ω = β + ξ with α ∈ a, β ∈ b, which in turn can be written α = ξ ′ + η , β = ζ ′ + η ′
with ξ ′ ∈ x, etc. This gives us a system (S) of 4 equations.

On the other hand, by definition, Γ(x, a, y, b, z) is the set of all ω ∈ Ω such that

∃ξ ∈ x, ∃α ∈ a, ∃η ∈ y, ∃β ∈ b, ∃ζ ∈ z : ω = ζ + α = α + η + β = ξ + β
There are several equivalent versions of this system (R) of three equations -see [START_REF] Bertram | The projective geometry of a group[END_REF], Lemma 2.3., from which it is read off that the four conditions from (S) are satisfied for ω ∈ Γ(x, a, y, b, z) if we choose ξ ′ = ξ, η ′ = η, ζ ′ = ζ. Thus the inclusion "⊂" holds always. The other inclusion does not always hold, but the theorem gives a sufficient condition: indeed, if ω belongs to the set on the right hand side, then (S) implies

ω = ξ ′ + η + ζ = ζ ′ + η ′ + ξ , whence ξ -ξ ′ ∈ y ∨ z. If x ∧ (y ∨ z) = 0,
this implies that ξ = ξ ′ , and three of the four equations from (S) are equivalent to (R). If y ∧ (x ∨ z) = 0 or z ∧ (x ∨ y) = 0, then the same argument applies (with respect to another choice of three from the four equations of (S)). In all cases, it follows that ω ∈ Γ(x, a, y, b, z).

(2) From [BeKi10a], Theorem 2.3, we know that Γ is para-associative and satisfies the idempotent law: Γ Γ(x, a, y, b, u), a, u, b, z = Γ(x, a, y, b, Γ(u, a, u, b, z)) = Γ(x, a, y, b, z) .

By assumption, the triple (x, y, u) is in general position, and from this it follows that the triple ((xyu), u, z) is also in general position; therefore the left-hand side may be expressed in terms of the lattice structure by applying twice part (1), which leads to the expression from the claim: in a first step, we get

(xyu) = (x ∨ y) ∧ a ∨ u ∧ (u ∨ y) ∧ b ∨ x ,
and in a second step

((xyu)uz) = (([( (x ∨ y) ∧ a ∨ u ∧ (u ∨ y) ∧ b ∨ x ] ∨ u) ∧ a] ∨ z) ∧ [((z ∨ u) ∧ b) ∨ (x ∨ y) ∧ a ∨ u ∧ (u ∨ y) ∧ b ∨ x ] .
(3) Under the given assumptions, the first term on the right hand side in (2) reduces to the line x ∨ y, and hence the claim follows directly from (2).

Remarks. (a) Not all possible relative positions of (x, y, z) are covered by Theorem 2.3, that is, the lattice theoretic formula for Γ(x, a, y, b, z) does not hold for all triples of submodules of Ω. For instance, if Ω = K 2n and x, y, z are of dimension n, then they cannot be in general position, and in general no u as in (2) exists. This case illustrates the special rôle of "generalized projective lines" (cf. [START_REF] Bertram | Associative Geometries. I: Torsors, Linear Relations and Grassmannians[END_REF]) with respect to lattice approaches.

(b) Both for the definitions given here and in [START_REF] Bertram | Associative Geometries. I: Torsors, Linear Relations and Grassmannians[END_REF], it is not strictly necessary that x, y, z belong to U ab : they may belong to V a , or to V b , or (in [START_REF] Bertram | Associative Geometries. I: Torsors, Linear Relations and Grassmannians[END_REF]) be completely arbitrary. We will not enter here into a discussion of the relation of both definitions if x, y or z does not belong to U ab .

(c) Both approaches lead to their own notions of morphisms. In the situation of Part (3) of Theorem 2.3, both of these notions must lead to the same result: this is precisely the famous "second fundamental theorem of projective geometry".

The Moufang case

Theorem 3.1. Assume X is a Moufang projective plane and (a, b) a pair of lines. Then U ab , together with the ternary product (xyz) ab defined in the first section, is a ternary Moufang loop. In particular, if we fix an element y ∈ U ab as origin, then U ab with product xz = (xyz) ab and origin y becomes a (binary) Moufang loop.

Before proving the theorem, we recall the relevant definitions (cf., e.g., [SB+]): Definition 3.2. A projective plane X is a Moufang plane if it satisfies one of the following equivalent conditions

(1) The group of automorphisms fixing all points of any given line acts transitively on the points not on the line. (2) The group of automorphisms acts transitively on quadrangles.

(3) Any two ternary rings of the plane are isomorphic. (4) Some ternary ring of the plane is an alternative division algebra, i.e., it is a division algebra satisfying the following identities:

x(xy) = (xx)y , (yx)x = y(xx) , (xy)x = x(yx) .

(5) X is isomorphic to the projective plane over an alternative division ring.

(6) The "small Desargues theorem" holds in all affine parts of X .

The set of invertible elements in alternative algebra forms a Moufang loop. A basic reference for loops in general and Moufang loops in particular is [START_REF] Bruck | A Survey of Binary Systems[END_REF]. y) → xy such that for each x, the maps y → xy and y → yx are bijections of Q, and having an element e such that ex = xe = x for all x ∈ Q.

Definition 3.3. A loop (Q, •) is a set Q with a binary operation Q 2 → Q; (x,
A Moufang loop is a loop Q that satisfies any, and hence all of the following equivalent identities (the Moufang identities):

z(x(zy)) = ((zx)z)y (M1) x(z(yz)) = ((xz)y)z (M2) (zx)(yz) = (z(xy))z (N1) (zx)(yz) = z((xy)z) (N2)
The left and right multiplication maps (sometimes called translations) in a loop are defined, respectively by L x y := xy =: R y x. The Moufang identities can be written in terms of the left and right multiplication maps. For instance, the first two identities state that

L z L x L z = L zxz and R z R y R z = R zyz .
Moufang's Theorem implies that Moufang loops are diassociative, that is, for any a, b, the subloop a, b generated by a, b is a group. This can be seen as a loop theoretic analog of Artin's Theorem for alternative algebras. Two particular instances of diassociativity are the left and right inverse properties

x -1 (xy) = y (LIP) (xy)y -1 = x , (RIP)
where x -1 is the unique element satisfying xx -1 = x -1 x = e. The following lemma gives the Moufang analog of the well-known relation between torsors and groups: Lemma 3.4. Let Q be a Moufang loop, and define a ternary operation (• • •) : Q 3 → Q by (xyz) := (xy -1 )z. Then the following three identities hold:

(xxy) = y = (yxx) (MT0) (uv(xyx)) = ((uvx)yx) (MT1) (xy(xyz)) = ((xyx)yz) (MT2)
Conversely, if M is a set with a ternary operation (•••) : M 3 → M satisfying (MT0), (MT1) and (MT2), then, for every choice of "origin" e ∈ M, the binary operation x • y := (xey) and the unary operation x -1 := (exe) define the structure of a Moufang loop on M with neutral element e.

Proof. Firstly assume Q is a Moufang loop. The leftmost identity in (MT0) is trivial while the rightmost follows immediately from (RIP). For (MT1), we compute

(uv(xyx)) = (uv -1 )((xy -1 )x) = (uv -1 )(x(y -1 x)) ( x, y is a group) = (((uv -1 )x)y -1 )x (M2) = ((uvx)yx) . For (MT2), (xy(xyz)) = (xy -1 )((xy -1 )z) = ((xy -1 )(xy -1 ))z ( xy -1 , z is a group) = (((xy -1 )x)y -1 )z ( x, y is a group) = (((xyx)y)z) .
Conversely, suppose M is a set with a ternary operation (•••) : M 3 → M satisfying (MT0), (MT1) and (MT2). Fix e ∈ M and define x • y := (xey) and x -1 := (exe) for all x, y ∈ M. By (MT0), we see that e is neutral element for the binary operation.

First we establish the following identities:

x • y -1 = (xye) , (3.1) (x • y -1 ) • z -1 = (xyz -1 ) , (3.2) (x -1 ) -1 • x = e , (3.3) ((x -1 ) -1 xy -1 ) = y . (3.4)
For (3.1) we compute x • y -1 = (xe(eye)) = ((xee)ye) = (xye) using (MT1) in the second equality and (MT0) in the third. For (3.2), we have (x • y -1 ) • z -1 = ((xye)ze) = (xy(eze)) = (xyz -1 ), using (3.1) (twice) and (MT1). For (3.3), (x -1 ) -1 •x = ((x -1 ) -1 xe) = ((ex -1 e)xe) = (ex -1 (exe)) = (ex -1 x -1 ) = e, using (3.1) in the first equality, (MT1) in the third and (MT0) in the fourth. Finally, for (3.4), ((x -1 ) -1 xy -1 ) = ((x -1 ) -1 • x -1 ) • y -1 = e • y -1 = y -1 , using (3.2) followed by (3.3).

Next we prove (3.5) (xy(y -1 ) -1 ) = x .

Indeed, (xy(y -1 ) -1 ) = ((x(y -1 ) -1 (y -1 ) -1 )y(y -1 ) -1 ) = (x(y -1 ) -1 ((y -1 ) -1 y(y -1 ) -1 ))

= (x(y -1 ) -1 (y -1 )y -1 ) = x , where we have used (MT0), (MT1), (3.4) and (MT0). Taking y = x in (3.5) and applying (MT0), we obtain

(3.6) (x -1 ) -1 = x .
From this it follows that e -1 = e, since e -1 = e • e -1 = (e -1 ) -1 • e -1 = e. Now in (3.2), replace z with z -1 and use (3.6) to obtain

(3.7) (x • y -1 ) • z = (xyz) .
Replacing y with y -1 and then setting z = y -1 in (3.7), we obtain (x•y)•y -1 = (xy -1 y -1 ) = x using (MT0). Thus the right inverse property (RIP) holds.

Next we almost obtain the Moufang identity (M2) as follows:

((x • y) • z) • y = (((x • e) • y) • z) • y = ((xey)z -1 y) = (xe(yz -1 y)) = (x • e) • ((y • z) • y) ,
using (3.7), (MT1) and (3.7) again. In loop theory, this is known as the right Bol identity.

We also have the left alternative law: (2) For an associative torsor (• • •) : M 3 → M, the groups determined by different choices of "origin," that is, fixed middle slot, are all isomorphic. The analog of this does not hold for ternary Moufang loops. Instead the different Moufang loops are isotopic [START_REF] Bruck | A Survey of Binary Systems[END_REF]. In fact, it is straightforward to show that for a Moufang loop Q, each isotope of Q is isomorphic to an isotope with multiplication given by x• z = (xy -1 )z for some y ∈ Q. Thus just as alternative triple systems encode all homotopes of an alternative algebra into a single structure, so do ternary Moufang loops encode all isotopes of Moufang loops.

(x • x) • y = (((x • e) • x) • e) • y = ((xex)ez) = (xe(xez)) = (x • e) • ((x • e) • y) = x • (x • y) , using (3
(3) Though we did not bother to state this in the lemma, it is clear from the proof that if we start with a Moufang loop Q with neutral element e, construct the corresponding ternary operation (• • •) and then construct the binary and unary operations induced by (• • •) with origin e, we recover the original loop operations. Similarly, if we start with a ternary Moufang loop M, construct the binary and unary operations with origin e and then construct the corresponding ternary operation induced by the loop structure, we recover the original ternary Moufang loop.

Proof of Theorem 3.1. In principle, the first two strategies of proof of Theorem 2.1 carry over:

(a) A proof in the framework of axiomatic geometry. Instead of the full Desargues theorem we now can only use the Little Desargues theorem. The drawings will become more complicated than above since one has to introduce auxiliary points. We will not pursue this proof here.

(b) A computational proof. Let K be the alternative division ring belonging to the plane. Then the affine space V := V a is isomorphic to K 2 , and affine lines can be described as in the Desarguesian case, eg. x ∨ y = {(1 -t)x + ty | t ∈ K}, where multiplication by "scalars" in K 2 is componentwise. If a = b, then (xyz) aa = x -y + z is the torsor law of the abelian group V a ∼ = (K 2 , +), and the claim is true. If a = b, fix some arbitrary origin o in the affine hyperplane V a ∩ b. There is a linear form β : V → K such that b ∩ V = ker(β), so that U ab = {x ∈ V a | β(x) ∈ K × }. (To fix things, one may choose coordinates such that β = pr 1 is the projection onto the first coordinate of K 2 , so b is the vertical axis.) Lemma 3.6. Let notation be as above. Then, for all x, y, z ∈ U ab , we have

(xyz) ab = (β(z)β(y) -1 ) • (x -y) + z.
Proof. The proof of Lemma 2.2 carries over without any changes -associativity of the ring has not been used there, only some elementary properties of inverses which are direct consequences of the left and right inverse properties (LIP) and (RIP).

From the lemma we get, as before, the formula

(3.8) β (xyz) = [β(z)β(y) -1 ]β(x) = β(z)β(y)β(x)
which means that β induces a homomorphism from U ab to the ternary Moufang loop K × . This formula is crucial in the proof of the alternative laws of U ab : essentially, it implies that identities holding in K × will carry over to U ab ; but the unit loop of K is a ternary Moufang loop, and hence so will be U ab . For instance, for the proof of (MT1), (uv(xyx)) = ((uvx)yx), write both sides, using the lemma: one sees that equality holds iff, for the vector w := u -v ∈ K 2 and for all x, y, v ∈ K 2 \ ker(β), we have (((β(x)β(y) -1 )β(x))β(v) -1 )w = (β(x)β(y) -1 )((β(x)β(v) -1 )w) But this amounts to an identity in K (or, if one prefers, two identities, one for each component of w), of the same form as the one we want to prove; this identity holds since K is an alternative algebra.

Prospects

In subsequent work, we will investigate more thoroughly the geometry corresponding to alternative algebras and alternative pairs (cf. [START_REF] Loos | Jordan Pairs[END_REF]): "alternative geometries" correspond to such algebras in a similar way as the associative geometries from [START_REF] Bertram | Associative Geometries. I: Torsors, Linear Relations and Grassmannians[END_REF] correspond to associative algebras and associative pairs. They play a key rôle in the construction of exceptional spaces corresponding to Jordan algebas and Jordan pairs. In the following, we briefly mention some topics to be discussed in this context. 4.1. Structure of the torsors and ternary Moufang loops. First of all, it is easy to understand the structure of the groups U ab in the Desarguesian case: for a = b, U ab = V a is a vector group (this is true even in the Moufang case), and for a = b, U ab is isomorphic to the dilation or ax + b-group

(4.1) Dil(E) := {f : E → E | f (x) = ax + b, b ∈ E, a ∈ K × }
of the affine space E = a b = a \ b (where a ∩ b is considered as hyperplane at infinity of a). This dilation group, in turn, is a semidirect product of K × with the translation group of E. The resulting homomorphism U ab → K × can be described in a purely geometric way (cf. [START_REF] Bertram | The projective geometry of a group[END_REF], Theorem 7.4 for the case of very general Grassmannians). For Moufang planes, partial analogs of this hold: there is a split exact sequence of ternary Moufang loops

a b → U ab → K × ,
where any line L in X which intersects a ∪ b in exactly two different points provides a splitting. But, if the plane is not Desarguesian, the set Dil(E) defined by (4.1) is then no longer a group, nor is it contained in the automorphism group of the plane. However, it remains true in the Moufang case that one obtains a symmetric plane (defined in [START_REF] Löwen | Symmetric planes[END_REF]; in the rough classification of symmetric planes by H. Löwe [START_REF] Löwe | A rough classification of symmetric planes[END_REF], our spaces appear among the split symmetric planes.)

4.2. Duality. Carrying out our geometric construction from Chapter 1 in the dual projective space, by general duality principles of projective geometry, we get again torsors, respectively ternary Moufang loops. Remarkably, the description of the torsors in the Desarguesian case by equation (1.4) does not change, except for a switch in a and b. In other words, up to this switch, the map Γ is "self-dual", which is in keeping with results on anti-automorphisms from [START_REF] Bertram | Associative Geometries. II: Involutions, the classical torsors, and their homotopes[END_REF]. For the moment, it is an open problem whether a similar "self-dual description" exists also in the Moufang case.

4.3. General projective planes. Our definition of (xyz) ab in the generic case (Definition 1.1) makes sense for any projective plane (and even for any lattice if we admit 0 as possible result). What, then, are its properties? In particular, what is its relation with the "ternary field" associated to a quadruple of points in the plane? Put differently, how do we have to modify the definition in the collinear case (Definition 1.2)? Does the "split exact sequence" a b → U ab → K × survive in some suitable algebraic category? Can one re-interprete the classical Lenz-Barlotti types of projective planes (cf., e.g., [SB+], p. 142) in terms of (xyz) ab ? 4.4. Perspective drawing. Our construction also has aspects that should be interesting for applied sciences: as already pointed out, our two-dimensional drawings have a "spacial interpretation". This can be explained by observing that the torsors U AB living in a three-dimensional space KP 3 can be mapped homomorphically onto torsors U ab living in a projective plane P (by choosing P ⊂ KP 3 intersecting A ∧ B in a single point and projecting from a point q ∈ A ∧ B, q / ∈ P , onto P ; then let a := P ∧ A and b := P ∧ B). A careful look shows that the torsor structure thus represented on P is quite often implicitly used in two-dimensional "perspective representations" of three-dimensional space; however, to our knowledge, the underlying algebraic structure has so far not yet been clearly recognized.

y

  w = (xyz) , y = (zwx) , z = (yxw) , x = (wzy) are all equivalent. Note also that it is obvious from the definition that(1.2) (xyz) ba = (zyx) ab .If we represent a and b by affine lines, intersecting in the affine drawing plane, the construction is visualized like this:If we choose a as "line at infinity" of our drawing plane, and if we choose to draw b horizontally, then we get the following image:

  .7), (MT2) and (3.7) again. The rest of the argument is standard. A magma satisfying the right Bol identity and (RIP) is a loop, called a right Bol loop (see, e.g., [Kie02], Theorem 3.11, suitably dualized). A right Bol loop satisfying the left alternative law is a Moufang loop [Rob66]. Definition 3.5. A set M with a map (• • •) : M 3 → M satisfying the three identities from Lemma 3.4 will be called a ternary Moufang loop.Remarks. (1) The axioms (MT1) and (MT2) for ternary Moufang loops are precisely the identities (AP2) and (AP3) in Loos' axiomatization of an alternative pair[START_REF] Loos | Jordan Pairs[END_REF].

  1.2. The collinear case. We have not yet defined what (xyz) ab should mean if x, y, z are collinear. If X is a topological projective plane, then one would like to complete our definition simply "by continuity", e.g., by taking the limit of (xyz) ab as y, not lying on the line x ∨ z, converges to a point on x ∨ z. This is indeed what happens in the classical planes over the division algebras R, C, H, O. Since we do not know whether in very general cases such a "limit" exists, we restrict ourselves here to the Moufang case, and leave the general case for later work.

Definition 1.2. Assume that X is a Moufang plane or a projective space of dimension bigger than 2. Consider a pair (a, b) of hyperplanes and a collinear triple (x, y, z) of points, none of them in a or b.

(1) If x = y = z, let (xyz) ab := x.