Fernando De Goes

David Cohen-Steiner

Pierre Alliez

Mathieu Desbrun

An Optimal Transport Approach to Robust Reconstruction and Simplification of 2D Shapes

We propose a robust 2D shape reconstruction and simplification algorithm which takes as input a defect-laden point set with noise and outliers. We introduce an optimal-transport driven approach where the input point set, considered as a sum of Dirac measures, is approximated by a simplicial complex considered as a sum of uniform measures on 0-and 1-simplices. A fine-to-coarse scheme is devised to construct the resulting simplicial complex through greedy decimation of a Delaunay triangulation of the input point set. Our method performs well on a variety of examples ranging from line drawings to grayscale images, with or without noise, features, and boundaries.

Introduction

Shape reconstruction from unorganized point sets is a fundamental problem in geometry processing: despite significant advances, its inherent ill-posedness and the increased heterogeneity of geometric datasets contribute to make current approaches still far from satisfactory. Even the 2D instance of this problem, i.e., the reconstruction of shapes in the plane, remains a challenge in various application areas including computer vision and image processing. Twodimensional point sets are often acquired from sensors or extracted from images, and are thus frequently hampered with aliasing, noise, and outliers. In addition, image-based point sets often exhibit a wide variety of features such as corners, intersections, bifurcations, and boundaries. This combination of noise, outliers, and presence of boundaries and features render most well-known strategies (including Poisson, Delaunay, or MLS-based approaches) deficient.

Shape reconstruction is also intimately linked to shape simplification. While a few authors (in particular in Computational Geometry) have restricted the issue of reconstruction to finding a connectivity of all the input points, the presence of noise and the sheer size of most datasets require the final reconstructed shape to be more concise than the input. Reconstruction and simplification are, however, often performed sequentially rather than in concert.

Instead, we jointly address reconstruction and simplification of 2D shapes through a unified framework rooted in optimal transport of measures. Specific benefits include (i) robust-ness to large amounts of noise and outliers; (ii) preservation of sharp features; (iii) preservation of boundaries; and (iv) guarantee that the output is a (possibly non-manifold) embedded simplicial complex.

Previous Work

To motivate our approach and stress how it fills theoretical and practical needs, we first review previous work in both reconstruction and simplification of 2D point sets.

Reconstruction. For noise free datasets, existing reconstruction methods vary mostly based on sampling assumptions. Uniformly sampled datasets can be dealt with using image thinning [START_REF] Melhi | A novel triangulation procedure for thinning hand-written text[END_REF], alpha shapes [EKS83] or r-regular shapes [Att97]; for non-uniform sampling, most provably correct methods rely on Delaunay filtering [START_REF] Amenta | The crust and the beta-skeleton: Combinatorial curve reconstruction[END_REF], with improvements on computational efficiency [GS01], sampling bounds [DK99, DMR99], and handling of corners and open curves [START_REF] Funke | Reconstructing a collection of curves with corners and endpoints[END_REF]DW02].

Noisy datasets have been tackled by a variety of methods over the last ten years [Lee00, CFG * 03, MD07, MTSM10], with recent successes including the extraction of selfintersecting curves [RVC11]. Most of these reconstruction approaches first perform noise removal through clustering, thinning, or averaging, which often leads to a significant blunting of features. Robustness to outliers has also beento a lesser extent-investigated, with approaches ranging from data clustering [Son10] and robust statistics [FCOS05], to k-order alpha shapes [KP08], spectral methods [KSO04] and 1 -minimization [ASGCO10]-but often at the cost of a Simplification. Polygonal curve simplification has received attention from many fields, including cartography, computer graphics, computer vision and other medical and scientific imaging applications. A common approach to simplification is to proceed in a fine-to-coarse fashion through decimation [AHPMW02,GZ05], whereas coarse-to-fine methods-such as the Douglas-Peucker algorithm (see [START_REF] Hershberger | An O(nlogn) implementation of the Douglas-Peucker algorithm for line simplification[END_REF])-proceed through refinement instead. When topological guarantees are called for, a fine-to-coarse strategy is often preferred; for example, [START_REF] Dyken | Simultaneous curve simplification[END_REF] uses intersection tests at each decimation step to preserve nesting of cartographic contours. Topological guarantees are indeed more difficult to obtain in the coarse-to-fine strategy, especially when sharp angles and singularities are present. Simplification methods can also be classified according to the metric they use to measure the accuracy of the simplified curves. Most methods are based on the Hausdorff distance or the Fréchet distance (see [AH-PMW02] and the references therein). However, these metrics do not handle noise well, and significant processing is required to filter spurious data. A distinct, but related endeavor is to concisely convey semantic aspects of a geometric dataset-this problem is usually referred to as abstraction [MDS09,DGDV11], but it does not fall within the realm of our work even if it is a form of simplification.

While specific data idiosyncrasies such as noise, outliers, features, and boundaries have been successfully dealt with individually, little or none of the previous work can handle them all concurrently.

Contributions

We depart from previous work by leveraging the versatile framework of optimal transport: we view shape reconstruction and simplification as a transport problem between measures (i.e., mass distributions), where the input points are considered as Dirac measures and the reconstructed simplicial complex is seen as the support of a piecewise uniform measure. The use of optimal transport brings forth several benefits, including a unified treatment of noise, outliers, boundaries, and sharp features. Our reconstruction al-gorithm derives a simplicial complex by greedily minimizing the symmetric Wasserstein transport metric between the Dirac masses that the input points represent and a piecewise uniform measure on the 0-and 1-simplices of the reconstructed complex. A fine-to-coarse strategy is proposed for efficiency, starting from the 2D Delaunay triangulation of the input point set and proceeding through repeated edge collapses and vertex relocation. Features thus emerge from our optimal transport procedure rather than from an explicit feature detection scheme. The resulting reconstruction is extracted from the simplified triangulation through edge filtering, guaranteeing an intersection-free output.

Approach

We first review the basics of optimal transport before providing details on how we formulate our reconstruction task using this powerful framework.

Optimal Transport Formulation

Optimal transport refers to the problem of optimizing the cost of transportation and allocation of resources [Vil10] (for applications in geometric computing, see [LD10, Mém11, MMDD11]). An intuitive example of optimal transport (used initially by Gaspard Monge in 1781) consists in determining the most effective way to move a pile of sand to a hole of the same volume-"most effective" here meaning that the integral of the distances the sand is moved (one infinitesimal unit of volume at a time) is minimal. This formulation of the problem is referred to as Monge's variational formulation and assumes that the sand is moved through a point-to-point mapping called the transport plan; this restriction was relaxed by Kantorovich who extended the formulation to deal with transport plans between two probability measures µ and ν. In this formulation, a transport plan is a probability measure π on support(µ)×support(ν) whose marginals are µ and ν. If µ and ν have finite support, π(x, y) specifies the amount of mass transferred from x to y. For any such transport plan we can associate a notion of cost. A common cost function is the q-Wasserstein metric wq, defined as

wq(µ, ν) = inf π R d ×R d
xy q dπ(x, y) 1/q . To reuse the example mentioned above, if each distribution is viewed as a unit amount of piled-up sand, this metric is the minimum cost of turning one pile into the other when the Lq distance is used to measure the cost of transport. Note that the measures µ and ν are assumed to be probability distributions, hence both with unit total mass; throughout this paper, we adopt a trivial extension of the w 2 distance in R 2 that incorporates the total amount of mass to be transported. More specifically, for any pair of measures µ and ν with bounded support in R 2 and common total mass M, we define the optimal transport cost W 2 as

W 2 (µ, ν) = √ M w 2 (µ/M, ν/M). (1)
Our general approach to reconstruct a simplicial complex T from a point set S can thus be phrased as follows:

Considering S as a measure µ consisting of Dirac masses, find a coarse simplicial complex T such that µ is well approximated by a linear combination of uniform measures on the edges and vertices of T .

The main advantage of the notion of optimal transport distance in the context of reconstruction is its robustness to noise and outliers [CCSM09]. Moreover, it captures two complementary notions of approximation between S and T . First, it measures a symmetric approximation error between shapes since W 2 (µ, ν) = W 2 (ν, µ), avoiding the notorious shortcomings of asymmetric distances. Second, the optimal transport cost measures the local defect of uniform sampling of the point set along the reconstructed edges; minimizing such a metric favors edges covering uniformly dense regions while preserving boundaries and sharp features (Fig 2).

Transport Plan & Cost from Points to Simplices

Denote by S ={p i } the input point set, for which every point p i is seen as a Dirac measure µ i centered at p i and of mass m i . The point set is thus considered as a measure µ = ∑ i µ i . Let us assume for now that we are given a triangulation T , and a point-to-simplex assignment which maps every input point p i to either an edge e or a vertex v of T -we will explain in Sec. 2.3 how this assignment is automatically de- termined. We regard each 0-and 1-simplex of T as a uniform measure supported over the simplex itself; that is, v is also seen as a Dirac measure, but e is a uniform 1D measure defined over the edge e. Equipped with these measures and knowing the prescribed point-to-simplex assignment, we can now derive a transport plan π from S to vertices and edges of a triangulation T .

Consider an arbitrary vertex v and an arbitrary edge e. Denote by Sv the set of points of S assigned to the vertex v and by Se the set of points of S assigned to the edge e; we assume these sets to be disjoint, with Points to Vertex. For a vertex v, the optimal plan to transport the measure Sv to the Dirac measure centered on v with mass Mv is trivial, and its local W 2 cost easily computed as

W 2 (v, Sv) = ∑ pi∈Sv m i p i -v 2 .
(2)

Points to Edge. For an edge e, the optimal plan (and its associated W 2 cost) is less straightforward to express as it requires transporting a mass Me from Se to a uniform measure of value Me/|e| defined over the edge e, where |e| denotes the edge length-notice that this edge measure integrates to Me, ensuring the existence of a transport plan. Because our transport cost is based on the L 2 distance, we can decompose the transport plan into a normal and a tangential component to e. The normal plan is a simple orthogonal projection, and its transport cost N is expressed as

N(e, Se) = ∑ pi∈Se m i p i -q i 2 , ( 3
)
where q i is the orthogonal projection of p i onto e.

The tangential plan is slightly more involved to derive, but its cost can also be done in closed form. We proceed by first sorting the projected points {q i } along e and partitioning the edge into card(Se) segment bins, where the i-th bin length is l i = (m i /Me)|e|. Each point de Goes et al. / An Optimal Transport Approach to RobustReconstruction and Simplification of 2D Shapes q i is then spread tangentially over the i-th bin so as to result in a uniform measure over e, defining the local optimal transport from Se to e (see inset). Now consider a point p i of mass m i that projects onto q i on edge e; set a 1D coordinate axis along the edge with origin at the center of the i-th bin, and call c i the coordinate of q i in this coordinate axis. The tangential cost t i of p i is then computed as the accumulated distribution resulting from the order-1 moment of all points on the i-th bin with respect to the projection q i , yielding:

t i = Me |e| li/2 -li/2 (x -c i ) 2 dx = m i l 2 i 12 + c 2 i .
We can now sum up the cost t i for every point p i in Se to get the tangential component T of the optimal transport cost for an entire edge:

T (e, Se) = ∑ pi∈Se m i l 2 i 12 + c 2 i . (4)
Note that Eq. ( 4) is exact; it is thus stable under refinement, in the sense that we obtain the same cost if we split each point p i of mass m i into several points that sum up to m i and transport them onto smaller bins whose lengths are function of the new masses.

The total W 2 cost to transport S to T through the transport plan π is thus conveniently written by summing the contributions of every edge and vertex of T :

W 2 (π) = ∑ e∈T N(e,Se) 2 +T (e,Se) 2 + ∑ v∈T W 2 (v,Sv) 2 .
(5)

Point-to-Simplex Assignment

Given a triangulation T , we still need to define an assignment of the input point set S to 0-and 1-simplices of T which minimizes the total cost of the transport plan π described above. Alas, solving this problem exactly is computationally infeasible since the combinatorial complexity of all possible assignments becomes intractable as the size of the problem increases. We propose a simple heuristic that constructs an assignment from S to edges and vertices of T in two steps. Each point p i is first temporarily assigned to the closest edge of the simplicial complex, resulting in a partition of S into a disjoint union of subsets Se, where each Se contains the input points nearer to e than to any other edge.

We then go over each edge e, and consider the two following assignments: either (i) keep Se assigned to e, or (ii) assign instead each point of Se to its closest endpoint of e. Using Eqs.

(2)-(4) we choose, out of these two possibilities, the local assignment leading to the smallest total cost. Once this simple test has been performed for each edge, the vertices v and edges e of T all have their assigned points (respectively, Se and Sv) from S, and the transport plan (and cost) for this assignment is thus fully determined as derived in Sec. 2.2.

Note that this assignment provides a natural characterization of the edges of T : an edge e is called a ghost edge when Se = ∅, i.e., no input points are transported onto it; conversely, e is called a solid edge if it receives a nonzero amount of mass from the point set. Ghost edges are, in effect, not part of the reconstruction, and are only useful to define the current simplicial tiling of the domain (see Fig. 1, center-right).

Algorithm

With the point-to-simplex assignment and its induced transport plan described in the previous section, we can now devise an algorithm to produce a coarse reconstruction T from the input point set.

Overview. The algorithm proceeds in a fine-to-coarse manner through greedy simplification of a Delaunay triangulation initialized from the input points. Simplification is performed through a series of half-edge collapse operations, ordered so as to minimize the increase of the total transport cost between the input point set and the triangulation. Vertices are also locally optimized after each edge collapse in order to further optimize the local assignment. Finally, the reconstruction is derived through edge filtering of the simplified mesh. We depict the pipeline in Fig. 1 and give pseudocode in Fig. 4.

Initialization

Our algorithm begins by constructing a 2D Delaunay triangulation T 0 of the input point set S. The resulting triangulation is augmented with four vertices placed at the corners of a loose bounding box of the input points; these vertices are not eligible for simplification as they act as pins ensur- ing that the sample points fall within the convex hull of T 0 . With this triangulation is associated a trivial optimal plan π 0 , for which each point is assigned to its corresponding vertex in T 0 for a total transport cost of zero.

Simplification

Simplification of the initial mesh T 0 is performed through half-edge collapses, each one removing one vertex and three edges from the triangulation (Fig. 5). As a collapse turns the current triangulation T k into a new triangulation T k+1 , it induces a local change of the transport plan π k between S and T k and, consequently, a change of the total cost W 2 (π k ) (denoted δ hereafter, with

δ k = W 2 (π k+1 ) -W 2 (π k )).
As our goal is to minimize the increase in total cost, we apply collapses in increasing order of δ. To this end, we initially simulate all feasible collapses to determine their associated δ and fill up a dynamic priority queue sorted in increasing order of δ. Decimation is then achieved by repeatedly popping from the queue the next edge to collapse, performing the collapse, updating the transport plan and cost according to Sec. 2.3, and updating the priority queue. Since our approximation of the transport plan partitions the input points based on their respective nearest edge, updating the transport involves only the edges confined to the one-ring of the removed vertex of the edge collapse (Fig. 6, center). Updating the priority queue, however, is required within a larger stencil around the removed vertex: edges incident to the modified one-ring and emanating from its flap vertices have their respective δ impacted by the new point-to-simplex assignment (Fig. 6, right).

Systematically Collapsing Edges.

A half-edge is said to be collapsible if its collapse creates neither overlaps nor foldovers in the triangulation. In order to preserve the embedding of the triangulation, we must only target collapsible edges during simplification. This is achieved by verifying both topological and geometric conditions: the former corre- sponds to the so-called link condition [DEGN99]; the latter consists of checking whether the target vertex of the halfedge is within the kernel [PS85] of the polygon formed by the one-ring of the source vertex (Fig. 7, left). While simple and commonly used, these validity conditions invalidate around 30% of the candidate collapse operators during the complete course of a simplification. This can severely affect the performance of our greedy approach to optimal transport as it relies on targeting the edge with the least cost first. To overcome this issue and ensure that the greedy decimation is systematically processed in increasing order of cost, we modify the collapse operator through a local edge flip procedure that makes every edge collapsible as we now review.

Suppose we want to collapse a half-edge (x i , x j ). We denote by Px i the (counter-clockwise oriented) polygon formed by the one-ring of x i , and by Kx i its kernel. We say that an edge (a, b) ∈ Px i is blocking x j if the triangle (x j , a, b) has clockwise orientation. We also call an edge flippable if its endpoints and its two opposite vertices form a convex quadrilateral. The idea of the flipping procedure is to elongate Kx i in the direction of (x i , x j ) until it includes x j . For every blocking edge (a, b) ∈ Px i , we call D as the distance from x j to the intersection point between the supporting lines of (x i , x j ) and (a, b). This distance indicates how much Kx i can be elongated by removing (a, b) from Px i . We then sort the blocking edges in a priority queue in decreasing order of distance D, and select the top edge to be removed from Px i . We remove an edge (a, b) from Px i by flipping either the edge (a, x i ) or the edge (b, x i ). Note that, since (a, b) has the greatest distance D within Px i , either (a, x i ) or (b, x i ) is flippable, otherwise (a, b) would not have been at the top of the queue. Between these two choices, we pick the edge with the smallest D when flipped, as it elongates Kx i the most. By repeating edge flips, we iteratively move Kx i towards x j and this procedure is guaranteed to terminate when the modified Px i no longer has blocking edges. It is worth pointing out that x i may not be inside the final Px i ; however, x i is later deleted by the collapse. Fig. 7 illustrates this deterministic flipping procedure.

Vertex Relocation

The simplification scheme presented so far restricts the successive triangulations {T k } obtained through simplification to have vertices on input points. However, noise and miss- ing data usually preempt having the exact location of sharp corners among the input points, making interpolated triangulations poorly adapted to recover features. In order to better preserve features, we perform a vertex relocation (in the same spirit as [START_REF] Garland | Quadric-based simplification in any dimension[END_REF]) after each half-edge collapse (Fig. 8).

The vertex relocation is designed to improve the fitting of vertices and edges of the triangulation to the input data: we move vertices in order to further minimize the normal component of the current W 2 distance (defined in Eq. ( 2) and Eq. ( 3)). Remembering that the square of the normal part of the W 2 cost associated with a vertex v of T can be written as

∑ pi∈Sv m i p i -v 2 + ∑ b∈N1(v) ∑ pi∈S (v,b) m i p i -q i 2 ,
we compute the optimal position v * of v by equating the gradient of the above expression to zero. If we denote by λ i the barycentric coordinate of the projection q i of the input point p i onto an edge (v, b), i.e., q i = (1λ i )v + λ i b, then the optimal position is given by:

v * = ∑ pi∈Sv m i p i + ∑ b∈N1(v) ∑ pi∈S (v,b) m i (1 -λ i )(p i -λ i b) Mv + ∑ b∈N1(v) ∑ pi∈S (v,b) m i (1 -λ i ) 2 .
In practice, we move the vertex v to v * only if the resulting triangulation T k+1 is still an embedding (i.e., v * is inside the kernel of the one-ring of v). Once this relocation is achieved, we proceed as we did for edge collapses: we collect the input points affected by this relocation, assign them to their nearest edge, and determine the new transport plan π k+1 .

Edge Filtering

The triangulation resulting from repeated decimation provides a tessellation of the domain that approximates the inferred shape by a subset of its edges. As described in Sec. 2.3, we already know whether an edge is relevant to the shape: only solid edges carry mass from the input point set. However, the presence of noise, outliers, and scaledependent features may still lead to a few undesirable solid edges (Fig. 1(center-right), in red). While the final edge filtering is most likely quite application-dependent and often unnecessary, we found it satisfactory to offer the option to eliminate solid edges based on a notion of relevance re, with re = Me|e| 2 N(e, Se) 2 + T (e, Se) 2 .

In our experiments, a threshold on relevance was only applied to examples with large amount of outliers (Figs. 1, 10, and 17).

Implementation details

We now describe a few implementation details which improve the efficiency and flexibility of the algorithm.

Initialization. For large datasets, a significant speed-up can be obtained by using only a random subset of the input points for T 0 . This trivial efficiency improvement does not lead to visible artifacts as we optimize the vertices of subsequent triangulations through our relocation procedure.

Multiple-Choice Priority Queue. The reconstruction of large datasets requires the simulation of thousands of collapse operations, and thus, our algorithm initially spends a significant amount of time performing frequent updates of the dynamic priority queue. We can speed up the simplification process through a multiple-choice approach [WK02] without significant degradation of our results: for each atomic simplification step we simulate only a random subset of the half-edges (typically 10) and collapse the half-edge from this subset with the smallest δ. When the mesh has 5 times the targeted number of vertices or less, we switch back to an exhaustive priority queue as the order of collapse for coarse meshes can significantly affect the final results. With this simple strategy, a five-to six-time speedup is typically obtained (see Tab. 1).

Postponing Vertex Relocation. Vertex relocation is effective only when a sufficient number of input points are assigned to each edge of the triangulation. Obviously, this number is very small at the beginning of the algorithm and increases as the triangulation becomes coarser. For this reason, we turn off the vertex relocation step during most of the simplification process, activating it only for the last 100 decimations. This strategy to lower computational cost also has little to no effect on the final results.

Reconstruction Tolerance. It is sometimes desirable to control the maximum error induced by the reconstruction, especially in the case of outlier-free point sets . We add an option which discards from the priority queue each collapse operator which would induce a transport with maximum normal and/or tangential cost above a user-specified tolerance. The algorithm terminates when the priority queue is empty (Fig. 9 is the only figure where this option was activated).

Results

We implemented our reconstruction and simplification algorithm using the CGAL library [CGA]. It takes as input point sets, possibly with mass attributes, and a desired vertex count of the final reconstruction. Timings for typical examples computed on 2.4GHz Intel Core 2 Duo with 2GB RAM are given in the following table.

On closed and noise-free curves, our approach performs The input shape (3K points) has sharp corners subtending small angles as well as boundaries. Our reconstruction is perfect for a noise-free input (left); as noise is added (middle, 2% and 2.5% of bounding box), the output degrades gracefully, still capturing most of the sharp angles; even after adding 4K or 4.5K outliers and 2% of noise (right), the reconstruction remains of high quality, although artifacts start appearing in this regime. comparably with state-of-the-art reconstruction and simplification methods applied in sequence. However, our approach is noticeably better at dealing with cases where the point set contains noise, outliers, and features. For instance, combinatorial methods such as Crust [GS01] fail to properly reconstruct a point set as soon as noise is added (Fig. 11); Poisson-like implicit approaches [START_REF] Kazhdan | Poisson Surface Reconstruction[END_REF] have also been reported to fail in this case unless reliable normals are provided [ACSTD07]. In contrast, our transport-based method is robust to large amounts of noise and outliers as illustrated in Fig. 10 where more outliers than initial points still leads to a good reconstruction-and sharp features are still captured (see also Fig. 12) without the typical blunting that most previous work (based on local denoising) would create. Fig. 13 depicts how the reconstruction gracefully degrades when the point density decreases. Our method also handles complex shapes with branching, intersections, and open curves without additional parameters to tweak (see Fig. 3).

Our approach performs equally well on unprocessed, nonsynthetic datasets. Images of line drawing, for instance, are particularly challenging for image processing techniques to approximate with polygonal curves. Nevertheless, our method is able to robustly extract such shapes by consid- ering the gray level pixels as points with a mass proportional to the pixel intensity. For instance, Fig. 14 depicts the reconstruction of a complex mechanical illustration (a comparison to a commercial image vectorization product is provided), and Fig. 15 shows the vectorization of a scanned artist drawing. We also show how our method performs on GIS data in Fig. 16, comparing our results with the well known Douglas-Peucker algorithm [START_REF] Hershberger | An O(nlogn) implementation of the Douglas-Peucker algorithm for line simplification[END_REF] and a direct implementation of the Garland-Heckbert simplification approach [START_REF] Garland | Quadric-based simplification in any dimension[END_REF].

Fig. 19(top) also illustrates how our approach infers shape at various scales, from fine to coarse. The original assignment for T 0 implies a vertex discretization of the datasets as all edges are initially classified as ghost edges; then these vertices are regrouped in 1-simplices, first consolidating the reconstruction of highly sampled regions. Finally, as the triangulation is simplified, we are getting a coarse scale approximation of the input data. Note that on inputs that mix anisotropic and isotropic point densities (Fig. 17) we obtain mixed reconstructions composed of line segments and isolated vertices, which also matches the input well.

Finally, Fig. 18 shows the versatility of our approach on a variety of datasets, varying from (kanji) characters, to scanned drawings, and even to gradient magnitudes obtained via finite-differences from an image (which provides a noisy sampling of the main contours of the image).

Conclusion

We have presented a novel optimal transport formulation and derived a practical algorithm to tackle the problem of reconstruction and simplification of 2D shapes from unorganized point sets. We view the point set as a measure that we greedily approximate, in the optimal transport sense, by a sum of uniform measures on the vertices and edges of a simplicial complex.

Strengths. The main value of our approach is the ability to robustly deal with feature preservation, such as sharp corners, cusps, intersections and boundaries. Another interesting property is the consistency between the underlying shape representation model (points and line segments) and the transport distance we considered: as this model inherently approximates features, we do not need any ad-hoc feature detection to preserve them.

Weaknesses. Despite satisfactory robustness to noise and outlier, our method does not deal gracefully with widely variable sampling. We see one such example in Fig. 19, where the reconstruction performs well on areas of high density but may fail to complete the undersampled area. This is quite natural given that our approach makes no assumption that the curve is closed, so undersampling can be interpreted as small features. As a consequence, our approach is also not resilient to missing data as the transport formulation seeks to preserve boundaries; however, very coarse reconstructions do recover the global shape quite well.

Future work. The problem of finding or approximating an optimal transport cost is of interest in many fields. In this work we have proposed a simple approximation strategy that assigns points to simplices in a local fashion, but more global alternatives (possibly through multiresolution) may be desirable. The simplification algorithm could also potentially benefit from a richer set of mesh decimation operators. For the final edge filtering step, we wish to incorporate application-dependent criteria so as to favor or enforce, for example, the reconstruction of 1-manifolds. Finally, the extension of this approach to robust surface reconstruction is probably the most exciting direction for future work; however, new decimation and transport cost evaluation strategies will be needed to offer both robustness and efficiency.

Figure 1 :

 1 Figure 1: Algorithm Pipeline. From left to right: input point set; Delaunay triangulation of input; after simplification, with ghost edges in grey, relevant solid edges in green, discarded solid edges in red; final reconstruction.

Figure 2 :

 2 Figure 2: Boundary and feature preservation. Left: Boundaries are preserved as an edge extending beyond (or falling short of) the point set would incur a large tangential component. Right: When an edge does not fit a sharp corner, the transport cost has both a large tangential component (due to the non-uniform sampling along the edge) and a large normal error. Minimizing transport cost will thus also induce alignment of edges to features.

Figure 3 :

 3 Figure 3: Leaf. Example with sharp corners, boundaries, and branching points.

  We further call Mv the total mass of Sv and Me the total mass of Se-thus, ∑ v∈T Mv + ∑ e∈T Me = M. Finally, we denote by π the plan satisfying the prescribed assignments from points to simplices, and W 2 (π) its transport cost.

Figure 4 :

 4 Figure 4: Pseudocode of the algorithm.

  c 2011 The Author(s) c 2011 The Eurographics Association and Blackwell Publishing Ltd. de Goes et al. / An Optimal Transport Approach to RobustReconstruction and Simplification of 2D Shapes

Figure 5 :

 5 Figure 5: Edge simplification. Half-edge collapse (arrow) triggers a local reassignment of the input points (assignments: red for edges, green for vertices).

Figure 6 :

 6 Figure 6: Updates during collapse. Left: collapse operator. Middle: edges with updated cost. Right: collapse operators pushed to the priority queue.

Figure 7 :

 7 Figure 7: Making edges collapsible. Left: collapsing the edge (in blue) creates fold-overs because of blocking edges (in black). Middle: the flipping procedure makes the edge collapsible. Right: after collapse.

Figure 8 :

 8 Figure 8: Vertex relocation. Noisy skyline point set, before (left) and after (right) vertex relocation.

Figure 10 :

 10 Figure10: Robustness to noise and outliers. The input shape (3K points) has sharp corners subtending small angles as well as boundaries. Our reconstruction is perfect for a noise-free input (left); as noise is added (middle, 2% and 2.5% of bounding box), the output degrades gracefully, still capturing most of the sharp angles; even after adding 4K or 4.5K outliers and 2% of noise (right), the reconstruction remains of high quality, although artifacts start appearing in this regime.

Figure 9 :

 9 Figure 9: Tolerance Control. A normal and tangential tolerance error (tolerance in gray, reconstruction in green) are imposed during simplification on the shark model.

Figure 11 :

 11 Figure 11: Bird. Top: For closed, noise-free point sets (left), combinatorial methods such as crust [GS01] (or other crust variants) correctly connects all input points (middle), while we obtain a simplified reconstruction (right). Bottom: With added synthetic noise (0.5% of bounding box), combinatorial approaches fail to recover the shape, while our method (right) returns a result nearly identical to the noise-free case.

Figure 12 :

 12 Figure 12: Noisy star. The lines of the inferred shape near the star center are visually indistinguishable (see closeup). Nevertheless, our method reliably recovers the lines and intersection point by interleaving decimation and relocation.

Figure 13 :

 13 Figure 13: Light house. 20-vertex reconstructions from input points with decreasing sampling density (from left to right: 2K, 400, and 200 points resp.).

Figure 14 :

 14 Figure 14: Mechanical part. Line drawing with complex features (left). Our reconstruction (middle) preserves intersection, corners, and manages to even remove small wiggles and glitches. In comparison, Adobe's Live Trace fails to handle glitches of this sketch (red arrows).

Figure 15 :

 15 Figure 15: Woman. Line drawing with variable thickness and outliers due to scanning artifacts (left). We apply thresholding on the input image (middle top) to reduce the number of input points before applying our reconstruction (right).

Figure 18 :

 18 Figure18: Gallery. Plate of shapes reconstructed through our optimal transport driven approach with input points originating from noiseless dataset (cartoon face contour), scanned drawings (screwdriver, tablecloth), grey-level images (kanji characters, flames), and image gradients computed by finite difference (all others).

Figure 16 :

 16 Figure 16: Australia Map. We compare our approach (top-right) with two curve simplification algorithms: QEM [GZ05] (bottom-left) and Douglas-Peucker [HS94] (bottom-right).

Figure 17 :

 17 Figure 17: Mixed reconstruction. By tuning the edge filtering we can extract a reconstruction made of either edges or of a mixture of edges and vertices.

Figure 19 :

 19 Figure 19: Sampling Density. Top: the transport cost favors the decimation of areas with high point density. As we further simplify the triangulation, regions with lower densities are also recovered. Bottom: areas with missing data are reconstructed as boundaries. ented point sets. In Symposium on Geometry Processing (2007), pp. 39-48. 7

  Simplification (Sec. 3.2) Create T k+1 by merging x i onto x j // Update Transport (Sec. 2.3) π k+1 := π k with local reassignments in N i,1

	// Relocation (Sec. 3.3)
	Optimize position of vertices in N i,1
	// Update Transport (Sec. 2.3)
	π k+1 := π k+1 with local reassignments in N i,2
	k ← k + 1
	until (desired vertex count)
	// Final Extraction (Sec. 3.4)
	Filter edges based on relevance (optional)
	Output: vertices and edges of Tn.

// TRANSPORT-BASED RECONSTRUCTION Input: point set S = {p 1 , . . . , pn}. // Initialization (Sec. 3.1) Construct Delaunay triangulation T 0 of S. Compute initial plan π 0 from S to T 0 k ← 0 repeat Pick best half-edge e = (x i , x j ) to collapse Set N i,m the m-ring of x i //

Table 1 :

 1 Timing. Left to right: Data set, number of points, reconstruction time in seconds using an exhaustive priority queue (EX) vs. a multiple choice (MC) decimation strategy.

	Data	Points	EX (s)	MC (s)
	Shark	440	2.6	0.5
	Bird	800	4.6	0.7
	Horse	1200	6	0.9
	Star	3k	18	3.6
	Australia	5k	39	7
	Falcon	8k	70	15
	Table cloth	20k	248	51

c 2011 The Author(s) c 2011 The Eurographics Association and Blackwell Publishing Ltd.

Acknowledgments. We thank Patrick Mullen for discussions along the way, and Xiaofeng Mi, Doug DeCarlo, and Ravish Mehra for providing data. This work was funded by the European Research Council (ERC Starting Grant "Robust Geometry Processing", Grant agreement 257474), and by NSF grants (CCF-0811373, CMMI-0757106, and CCF-1011944).


Efficient fault monitoring with Collaborative Prediction

Dawei Feng [START_REF] Blodget | Amazon's cloud crash disaster permanently destroyed many customers' data[END_REF], Cécile Germain-Renaud (1), Tristan Glatard [START_REF] Srebro | Maximum-margin matrix factorization[END_REF] (1) Laboratoire de Recherche en Informatique, CNRS Université Paris-Sud [START_REF] Srebro | Maximum-margin matrix factorization[END_REF]

CREATIS, CNRS INSERM Université Lyon 1 INSA Lyon

Overview

Isolating users from the inevitable faults in large distributed systems is critical to Quality of Experience. We formulate the problem of probe selection for fault prediction based on end-to-end probing as a Collaborative Prediction (CP) problem. On an extensive experimental dataset from the EGI grid, the combination of the Maximum Margin Matrix Factorization approach to CP and Active Learning shows excellent performance, reducing the number of probes typically by 80% to 90%.

Problem statement

The recent crash of the Amazon Cloud [START_REF] Blodget | Amazon's cloud crash disaster permanently destroyed many customers' data[END_REF] highlighted the importance of timely discovery of failures in large-scale distributed systems: a local, limited error may result in a global catastrophe. Thus a significant part of the software infrastructure of large scale distributed systems, whether grids or clouds, collects information (monitoring) that will be exploited to discover (knowledge) if, where, and when the system is faulty. This paper addresses the knowledge-building step in the context of end-to-end probing as the class of monitoring techniques. In the end-to-end probing approach, a probe is a program launched from a reliable entry point (probe station), which tests the availability of the components on its path. The only observables are the outcomes of the probes, which are binary: success or failure. For instance, a ping command would test the overall software stacks and network connectivity from the probe station to the target computer. This abstract summarizes our approach to minimize the number of probes for a given discovery performance target. Minimizing the number of probes can be addressed along three avenues: fault prediction, detection and diagnosis. In all cases, the system under consideration is a set of hardware and software components, featuring some dependencies -e.g. a service certainly depends on the hardware it is running on, and possibly of other services. The obvious advantage of detection and diagnosis is that they provide an explanation of the failure, by exhibiting culprits. On the other hand, they strongly rely on a-priori knowledge, namely which components are required for a probe to succeed. For massively distributed systems, where Lamport's famous definition "A distributed system is one in which the failure of a computer you didn't even know existed can render your own computer unusable" applies, assuming such knowledge might be hazardous in principle. Instead, this paper focuses on fault prediction. In this case, the overall infrastructure is a black box, with no a-priori knowledge of its structure.

Methods

A. Motivation for Collaborative Prediction

The motivating application comes from operations management in the European Grid Initiative (EGI), and specifically the Biomed Virtual Organization. Biomed has access to 256 Computing Elements (CEs) and 121 Storage Elements (SEs). CEs are shares of computing resources, implemented as queues of each site manager (e.g. PBS), and SEs are shares of storage resources; the formal definition is part of the Glue Information model [START_REF] Andreozzi | Glue Schema Specification, V.2.0[END_REF]. The probes test the availability of all relations between its endpoints, namely Computing Elements (CEs) and Storage Elements (SEs). Testing the availability of all CE-SE pairs is one of the most challenging issues encountered daily by monitoring operators. Brute force, that is periodically launching a fully distributed all-pairs availability test, for a total of 29512 tests, multiplied by the number of capacities to test at each run, is both intrusive and useless: human operators cannot handle so many results; in practice, only a few issues are reported, with questionable selection criteria. With our method, a massive reduction of the number of tests provides nearly similar availability evaluation performance, creating opportunities for better frequency/intrusiveness trade-off and selection of reported incidents. In this context, fault prediction can be considered as a case for Collaborative Prediction (CP). CP is originally a technique for predicting unknown ratings of products for a particular user, based on observed data from other users and products. The success of CP relies on the hypothesis of a factorial model: hidden and partially shared factors affect the ratings. In our case, two nodes (CE or SE) may share several hidden factors -e.g. location, with the associated network connectivity issues, or use of a particular instance of any middleware service (e.g. brokering, authentication), such that the availability of the CE-SE relation may be affected similarly.

B. Selection and prediction

Minimizing the number of probes encompasses two issues: probe selection, i.e. which subset of the (CE,SE) pairs to actually test; and predicting the availability of all pairs from the outcome of these probes.

Probe selection.

We consider two probe selection strategies.

-Static-Uniform. The probes are selected uniformly at random. In this setting, the prediction step has no influence over the choice of the probes. -Active Probing. With Active Probing (Algorithm 1), the set of probes is constructed dynamically, with an initial set of probes selected for instance by the Static-Uniform method, and run through the system to get basic information; then, additional probes are selected and launched with the goal of maximizing some measure of information; here we used the min-margin heuristic [START_REF] Porteous | Multi-hdp: a non parametric bayesian model for tensor factorization[END_REF], which chooses the probe where the uncertainty of the classification result is maximal, and has been demonstrated to be efficient for CP problems [START_REF] Griffiths | Finding scientific topics[END_REF].

Collaborative Prediction with MMMF

Algorithm 1: Generic active probing algorithm This section sketches Maximum Margin Matrix Factorization as proposed by Srebro et al. in [START_REF] Srebro | Maximum-margin matrix factorization[END_REF]. CP is formalized as a matrix completion problem: if Y is an observed (sparse) matrix, the CP problem is to find a full, real-valued, matrix X of the same size that approximates Y, i.e. that minimizes the "discrepancy" between X and Y, without any external information. Assuming a linear factor model, where k hidden factors define the user preference through a linear combination of them, X is constrained to be of rank k. However, bounding k to small values (low-rank approximation) does not lead to feasible optimization problems for a partially observed matrix and for the binary setting. The key insight in MMMF is to replace the bounded rank constraint with bounding the trace norm of X, under the constraint of no (hard-margin), or small (soft-margin), discrepancy.

Experimental setting

Different capabilities have to be tested; in the following, we consider three of them: probe srm-ls tests the list ability from a CE to a SE, probe lcg-cr tests the read ability from a CE to a SE, and probe lcg-cp tests the write ability alike. Thus, each CE works as a probe station, launching probes to test the functionalities between itself and each SE. For the Biomed grid a whole set of 29512 testing transactions were launched each day for each of the three probe classes. After nearly two months running, information for 51 validated days were collected. The probes themselves are gLite jobs, run by a regular Biomed user. Some of them fail (rejection) in the sense that gLite is not able to complete the job, denoting that some job management services may be down or misconfigured (e.g. authentication, brokering etc.). In the following, we consider only the accepted probes, i.e. those which run to completion, reporting success or failure; this approach amounts to consider that the data access capacities are independent from job management. This is a reasonable hypothesis in a gLite infrastructure because file transfers involved in job management use dedicated storage space independent from the one tested by our probes. Figure 1 illustrates the structure for lcg-cr on day'07-05-2011', where rows represent CEs and columns stand for SEs. Each entry is the probe result between the corresponding CE and SE. Black columns correspond to sustained SE downtimes while black lines are CE failures leading to complete inability to communicate with any SE (e.g. network downtime or configuration issue). These are usually easily detected and reported by human operators with only a few incident reports. The scattered points, however, correspond to local or transient issues, which are very difficult to handle due to the amount of incident reports independently generated. It is thus worth assessing the performance of the methods when sustained systematic errors are eliminated. To do that, we designed a second set of experiments, with curated matrices as the reference fault structure. A curated matrix is a new original matrix, where the lines and columns with only failed entries (black ones in figure 2) have been removed prior to analysis.

Experimental results

A. Static-Uniform. Figure 2 shows the accuracy (ratio of correctly predicted entries over total number of entries) for five randomly selected days. An excellent performance can be reached with a tiny fraction of the original probes, typically 5%. The baseline, a random guess following the distribution of the sample set, is plotted for comparison purpose, but can be computed a-priori. Figure 3 is the classical visualization of the confusion matrix in the ROC space for all the 51 days at 90% deletion rate (keeping 10% of the probes). Perfect prediction would yield a point in the upper left corner at coordinate (0,1), representing 100% sensitivity (no false negatives) and 100% specificity (no false positives). The srm-ls dataset shows nearly perfect prediction performance, being mostly very close to (0,1); lcg-cp and lcg-cr are vastly better than a random guess, which lies on the diagonal line (note the range of the axes, which cover only the part of the ROC space where the results belong, thus the diagonal line is not visible on the plot). The other classical indicators also show excellent performance: Area Under ROC Curve (AUC) as well as MCC (Matthews Correlation Coefficient) are close to 1. The problem becomes much more difficult when the systematic faults are excluded, thus taking the curated matrices as inputs. The prediction accuracy remains excellent; however, as the number of failed entries left in the curated matrices is much less than in the un-curated ones, e.g. the fraction of failed entries on srm-ls, 04-21-2011 drops from 45.37% to 2.25%, accuracy is not meaningful: predicting all entries as negative would give a similar result. The strategy to tackle this issue is Active Probing. B. Active probing. In this experiment, we compare the Active Probing strategy with the Static one at equal probing cost: first, a Static-Uniform method is applied, in order to get the reference information, then more probes are selected with the min-margin heuristic for Active Probing, while for the Static-Uniform method, the same number of probes are selected uniformly at random. Active Probing does improve accuracy over Static-Uniform. However, as explained in the previous section, the quality of failure prediction is the most important goal in this context. Figure 4 compares two relevant indicators: sensitivity (the probability of detecting an actual failure), and precision (the closer to 1, the smaller the risk of a false alarm). They are detailed for the initial probe fraction equal to 5%, then adding probes by step of 5% fractions. The results as given for a total of 10% and 15% probes. The first important result is that Active Probing always outperforms Static-Uniform. In fact, the performance of the Static-Uniform strategy is bad, except for 07-05-2011 (day 5). The performance bears some relation with the failure rate of the benchmark (table I): larger failure rates in the original curated matrix help uncovering the structure of the faults, even at quite low levels: with 4% failure rate, the 07-05-2011 benchmark exhibits acceptable performance when keeping only 5% of the probes and the Static-Uniform strategy. However, the failure rate does not tell the full story: days 2 and 4 have the same low one (1%), but the performance on day 4 is much better than on day 2. The likely explanation is that faults at day 2 do not present much correlation, while faults day 4 derive from a small number of shared causes. The second result is that Active Probing is quite efficient. We have a good chance of predicting the actual failure, except on day 2, while the false alarm rate remains negligible. A third strategy, omitted here for brevity, harnesses active probing with a differentiated cost policy [START_REF] Morik | Combining statistical learning with a knowledge-based approach -a case study in intensive care monitoring[END_REF][START_REF] Lin | Statistical properties and adaptive tuning of support vector machines[END_REF], where the false positives are more penalized. In all case, except day 2, sensitivity increases by 10-20%, while precision remains almost unharmed. A complete description and more experimental results can be found in [START_REF] Feng | Distributed Monitoring with Collaborative Prediction[END_REF].

Related Work

Collaborative Prediction associated with end-to-end probing, with the components structure considered as a black box, participates in the general Quality of Experience approach [START_REF] Rifai | A brief synthesis of QoS-QoE methodologies[END_REF]. More precisely, an important ingredient separating QoE from QoS is binary (possibly extended to discrete) classification. Most work in this area is devoted to network-based services (e.g. among many others [START_REF] Tran | QoE model driven for network services[END_REF]). Before QoE became a popular keyword, Rish and Tesauro [START_REF] Rish | Estimating end-to-end performance by collaborative prediction with active sampling[END_REF] explored the combination of MMMF and Active probing for the selection of good servers in various distributed and P2P systems. Our work combines the goal of proposing fault-free services to the user exemplified in [START_REF] Rish | Adaptive diagnosis in distributed systems[END_REF], [START_REF] Zhengandm | Collaborative reliability prediction of service-oriented systems[END_REF], and the CP approach of [START_REF] Rish | Estimating end-to-end performance by collaborative prediction with active sampling[END_REF]. Z. Zheng and M.R. Lyu propose explicit users collaboration for estimating failure probabilities [START_REF] Zhengandm | Collaborative reliability prediction of service-oriented systems[END_REF]; while their end-to-end framework is related to ours, the goal is more in line with QoS (estimating the full distribution of probability instead of binary classification), and the correlation amongst users or services is modeled by ad hoc methods. To our knowledge, this work is the first one to exemplify MMMF on fault prediction.

Perspectives

The Achille's heel of large-scale production grids and clouds is reliability. Quality of Experience at reasonable human cost requires extracting the hidden information from monitoring data whose intrusiveness should be limited. Collaborative Prediction is one of the promising and scalable strategies that can address this goal. This paper demonstrates its effectiveness on a large experimental dataset. Further work will consider two avenues. Exploiting the resulting knowledgethe hidden causes uncovered by CP -towards diagnosis might be easier with alternative matrix factorization techniques; specifically Bi-LDA [START_REF] Porteous | Multi-hdp: a non parametric bayesian model for tensor factorization[END_REF] provides a quantified model of latent factors, which, like LDA, could be exploited for interpretation [START_REF] Griffiths | Finding scientific topics[END_REF] more easily than MMMF, but suffers from lower performance. The second avenue targets improving prediction by separating transient failures from more permanent non-systematic ones, with temporal models inspired from the hierarchical Hidden Markov Models approach of text mining.

Figure 2 .Figure 3 .

 23 Figure 2. Accuracy for the Static-Uniform selection (five days)

Figure 4 .

 4 Performance comparison for Static-uniform and Active Probing, curated srm-ls for the five benchmark days.