Noname manuscript No.
(will be inserted by the editor)
Feature-Preserving Surface Reconstruction
and Simplification from Defect-Laden Point Sets
Julie Digne · David Cohen-Steiner · Pierre Alliez ·
Fernando de Goes · Mathieu Desbrun
the date of receipt and acceptance should be inserted later
Abstract We introduce a robust and feature-capturing
surface reconstruction and simplification method that
turns an input point set into a low triangle-count sim-
plicial complex. Our approach starts with a (possibly
non-manifold) simplicial complex filtered from a 3D
Delaunay triangulation of the input points. This ini-
tial approximation is iteratively simplified based on an
error metric that measures, through optimal transport,
the distance between the input points and the current
simplicial complex—both seen as mass distributions.
Our approach is shown to exhibit both robustness to
noise and outliers, as well as preservation of sharp fea-
tures and boundaries. Our new feature-sensitive metric
between point sets and triangle meshes can also be used
as a post-processing tool that, from the smooth output
of a reconstruction method, recovers sharp features and
boundaries present in the initial point set.
Keywords Optimal transportation · Wasserstein dis-
tance · Linear programming · Surface reconstruction ·
Shape simplification · Feature recovery.
Mathematics Subject Classification (2000) 65D17 ·
65D18
1 Introduction
Surface reconstruction is a multi-faceted challenge which
precise problem statement depends on the nature and
defects of the input data, the properties of the inferred
surface (smooth vs piecewise smooth, with or without
boundaries), and the desired level of detail one wishes to
J. Digne · D. Cohen-Steiner · P. Alliez
Inria Sophia Antipolis - M´editerran´ee
F. de Goes · M. Desbrun
California Institute of Technology
capture. Despite a number of major contributions over
the past decade [5,24], achieving both feature preser-
vation and robustness to measurement noise and out-
liers remains a scientific challenge—and a pressing re-
quirement for many reverse engineering and geometric
modeling applications. Furthermore, low polygon-count
reconstructions has only received limited attention de-
spite the increase of point density in 3D scanning tech-
nology and the need for efficient subsequent geometry
processing.
In this paper we contribute a reconstruction method
that simplifies an initial (possibly non-manifold) tri-
angulation of the input point set, based on an error
metric that quantifies through optimal mass transport
the distance between the current simplicial complex and
the input points. Our reconstruction approach inherits
the qualities of our optimal transport based metric: it
is resilient to noise and outliers, can handle uneven
sampling, yet it finely captures boundaries and sharp
features. We demonstrate these distinguishing proper-
ties on a series of examples. Applied to the output of
a feature-lossy reconstruction method, our new metric
can also be used in order to recover sharp features and
boundaries through a vertex relocation process.
2 Previous Work
We first discuss existing surface reconstruction meth-
ods, restricting our review to approaches that are ro-
bust to noise and outliers as well as feature preserving.
We then point out recent, relevant work on geometry
processing based on optimal transport.
2 Julie Digne et al.
2.1 Surface Reconstruction
A common approach to robust surface reconstruction
from defect-laden point sets involves denoising and fil-
tering of outliers, and often requires an interactive ad-
justment of parameters. Automatic methods such as
spectral methods [25,47,4] and graph cut approaches [20,
26] have been proven extremely robust but are better
suited to the reconstruction of smooth, closed surfaces.
More recently, Cohen-Or and co-authors have proposed
a series of contributions based on robust norms and
sparse recovery [29,21,7]. An interpolating, yet noise
robust approach was alternatively proposed by Digne
et al. [13] through the construction of a scale space.
Feature preserving methods are typically based on
an implicit representation that approximates or inter-
polates the input points. In [14], for instance, sharp fea-
tures are captured through locally adapted anisotropic
basis functions. Adamson and Alexa [1] proposed an
anisotropic moving least squares (MLS) method instead,
using ellipsoidal mapping functions based on princi-
pal curvatures. More recently, Oztireli et al. [35] ex-
tended the MLS surface reconstruction through kernel
regression to allow for much sharper features. How-
ever, none of these techniques returns truly sharp fea-
tures: reconstructions are always semi-sharp, that is,
still rounded with various degrees of roundness depend-
ing on the approach and the sampling density. More-
over, the presence of sharpness in the geometry of a
point set is detected only locally, which often leads
to fragmented creases; the reconstruction quality thus
degrades quickly if defects and outliers are present.
Another way to detect local sharpness within a point
set consists in performing a local clustering of estimated
normals [34]: if this process reveals more than one clus-
ter of normals, then the algorithm fits as many quadrics
as the number of clusters. Improved robustness was
achieved in [16] by segmenting neighborhoods through
region growing. Lipman et al. [28], instead, proposed
a systematic enrichment of the MLS projection frame-
work 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., [19,36]).
To reduce crease fragmentation, a different thread
of work aims at extracting long sharp features. Pauly
et al. [37], for instance, used a multi-scale approach
to detect feature points, and constructed a minimum-
spanning tree to recover the most likely feature graph.
Daniels et al. [12] used a robust projection operator
onto sharp creases, and grew a set of polylines through
projected points. Jenke et al. [23] extracted feature lines
by robustly fitting local surface patches and by com-
puting the intersection of close patches with dissimilar
normals.
Shape simplification has also been tackled in [8], but
in a coarse-to-fine manner: a random initial subset of
the input point cloud and a signed distance function
over the set is built. Using this function, points are
added until a significant number of points lie within
an error tolerance. The augmented set is triangulated
and a surface mesh is reconstructed. Thus the method
also interleaves reconstruction with simplification. In
[2] and [3], a surface is reconstructed through point set
simplification and local coarsening or refinement of the
mesh. Salman et al. [45] proposed to detect features
within the point set (as in [33]) and combined Delaunay
refinement over features and Poisson reconstruction on
smooth parts of the inferred surface [24].
Contributions. In this paper, we adopt a very dif-
ferent reconstruction methodology: we reconstruct a sha-
pe through an iterative, feature-preserving simplifica-
tion of a simplicial complex constructed from the input
point set. To achieve noise and outlier robustness, an er-
ror metric driving the simplification is derived in terms
of optimal transport between the input point set and
the reconstructed mesh, both seen as mass distributions
(or equivalently, probability measures) in R
3
.
Next we provide a brief review of optimal transport,
and mention its applications to various problems in
computer graphics and computer vision.
2.2 Optimal Transport
The problem of transporting a measure onto another
one as a way to quantify their similarity has a rich
scientific history. For two measures µ and ν defined over
R
3
and of equal total mass (i.e., their integrals are the
same), the L
2
optimal transport from µ to ν consists
in finding a transport plan π that realizes the following
infimum:
inf
Z
R
3
kx − yk
2
dπ(x, y)
π ∈ Π(µ, ν)
,
where Π(µ, ν) is the set of all possible transport plans
between µ and ν [46]. In a nutshell, a transport plan
π is a displacement that maps every infinitesimal mass
from the input measure µ (here, the set of points) to the
target measure ν (here, the simplicial complex). This
formulation is particularly well suited to comparing 1D
measures such as histograms over the real line or on
the circle [40], and it has been used for transferring
color and contrast between images [43,39]. For applica-
tions in higher dimensions such as 2D or 3D shape re-
trieval [44,41] and segmentation [38], the optimal trans-
port formulation is notoriously less tractable: solving
Feature-Preserving Surface Reconstruction and Simplification from Defect-Laden Point Sets 3
the optimal transport problem requires linear program-
ming (LP). The LP formulation of optimal transport
has been used in applications such as surface compar-
ison [30] and displacement interpolation [9]. Attempts
to design computationally simpler surrogates have also
been made; the sliced Wasserstein approach [42], for
instance, consists in approximating the transport prob-
lem by a series of 1D problems through projection.
Fig. 1 Transport plans. Top: Binary transport plan [18]. 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. A closeup reveals
a spurious tangential component of the transport near corner
vertices (as indicated by red segments pointing towards the
corner), artificially creating a higher total transport cost. For
more complex features such as concave corners on surfaces,
such behavior leads to reconstruction artifacts. Bottom: our
transport plan is, as expected, mostly normal to the edges.
Contributions. Our reconstruction method also re-
lies on a linear programming formulation. However, our
approach introduces a key distinctive property: our tar-
get 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 transport-
ing the input pointwise measure (i.e., the initial point
set) to the complex simplices. This specific setup bears
a resemblance to what is known as the optimal location
problem [31]), where the source measure is given but the
target measure is only partially known. Yet a significant
difference lies in the type of constraints we are enforcing
on the target measure, rendering current computational
methods to solve this problem not appropriate to our
context. Another line of research for finding a trans-
portation plan between an input point set and a set
of discrete sites of various capacities [6,22,32] make
use of power diagrams, and are thus likely to be too
computationally costly for our context.
The closest work to ours was proposed by de Goes et
al. [18]. Their algorithm reconstructs and simplifies 2D
shapes from point sets also based on optimal transport.
Nonetheless, our approach differs from theirs 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 trans-
porting points to the facets of a simplicial complex.
Therefore we use a discretized formulation of the
optimal transport problem.
2. the authors of [18] propose to approximate the op-
timal transport plan by assigning each input point
to its closest edge in the triangulation. Such a sim-
plistic scheme can lead to a sub-optimal transport
plan and cost as illustrated in Figure 1(top)—even
more so in 3D. Our discretized formulation, com-
bined 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 recur-
sive edge flip procedure. Such an edge flip proce-
dure can not, however, be generalized to 3D tri-
angulations. Instead, our method removes the em-
bedding requirement by only employing a (possibly
non-manifold) simplicial complex, initially chosen as
a subset of a 3D Delaunay triangulation.
2.3 Overview
Motivated by the concept of reconstruction introduced
in 2D by de Goes et al. [18], we present a fine-to-coarse
algorithm which reconstructs a surface from a point set
through greedy simplification of a 3D simplicial com-
plex. We initialize the complex with a (possibly non-
manifold) subset of the 3D Delaunay triangulation of
input points, then we perform repeated decimations
based on half-edge collapse operations. The error met-
ric guiding our simplification is derived from the op-
timal cost to transport the input point set (seen as
4 Julie Digne et al.
Dirac measures) to a constant-per-facet measure de-
fined over the simplicial complex. At each iteration,
we collapse the half-edge which minimizes the increase
of total transport cost between input points and re-
constructed triangulation. Just like for the formulation
presented in [18], our optimal transport driven metric
brings desirable properties that are rarely satisfied by
current reconstruction methods, such as resilience to
noise and outliers, and preservation of sharp features
and boundaries.
In the remainder of this paper, we first discuss the
details of our optimal transport based metric (Sec. 3)
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.
Algorithm 1: Algorithm overview.
Input : Point set S, user-specified value V .
Output: Simplicial complex C with V vertices.
Construct 3D Delaunay Triangulation T from S;
Compute transport cost from S to facets of T ;
Construct simplicial complex C from facets of T with
non-zero measure;
Decimate C until desired number of vertices V ;
Filter out facets of C by thresholding mass density.
3 Transport Formulation of Reconstruction
We consider the reconstruction problem of turning an
input point set S into a coarse simplicial complex C. The
point set contains N points at locations {p
i
}
i=1···N
, and
each point is given a mass m
i
that reflects its measure-
ment confidence (all masses are set to a constant if no
confidence is provided). Our reconstruction method is
based on considering both the point set and the com-
plex as mass distributions (or equivalently, probability
measures), where the measure (mass density) of C is
constant per simplex and possibly 0. Our approach then
consists in finding a compact shape C that minimizes
the optimal transport cost between the input point set
S and a uniform measure on each facet and vertex of
C.
In [18], a similar, yet 2D optimal transport cost
between S and C was efficiently approximated based
on closed form expressions for the optimal cost be-
tween points and edges. However, to our knowledge,
such closed form can not be extended between points
and triangles. We present instead, using linear program-
ming (LP), a discretized formulation of the optimal
transport between S and C that we will solve for later
on through local relaxations.
Fig. 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 centroids (green arrows); (e-f) Intermediary decimation
steps; (g-i) Reconstruction with 100, 50, and 22 vertices,
respectively; (j-l) Final transport plan with 100, 50, and 22
vertices, respectively.
3.1 Discretization
We approximate the optimal transport cost between
the input point set S and the simplicial complex C
using quadrature. We start by defining a set B of bins
(small regions of the complex) over C. As we aim at
reconstructing piecewise smooth surfaces from point
sets, facet bins are necessary—edge bins could be used
as well if curves in R
3
were sought after as well; for
simplicity, we do not discuss this extension. However,
vertex bins are useful as well: when outliers are present,
vertices serve as garbage collectors. Vertex and facet
bins are thus used to evaluate the optimal cost between
S and C as a sum of squared distances between the
points in S and the centroids of the bins in B.
Feature-Preserving Surface Reconstruction and Simplification from Defect-Laden Point Sets 5
First, every vertex of C is considered as (the cen-
ter of) its own bin; each triangular facet is, instead,
tiled with bins using a 2D Centroidal Voronoi Tessella-
tion (CVT) (Figure 3); note that our choice of a CVT
tiling stems from the fact that it minimizes the ap-
proximation error given by quadrature points put at
their centroids [15], which will thus provide optimal
approximation of our transport cost. The number of
bins per facet is set based on a user-defined quadrature
parameter. In all our experiments, we used 200 bins
per unit area, the point sets being included in a half-
unit side size box (note that the facet bins of fig 3 are
purposedly generated with a higher density). Finally,
to compensate for a slightly non-uniform distribution
of bins, we assign a capacity for each bin in B (i.e., ratio
of the total amount of mass that a bin can receive over
the total amount of mass transported to the simplex
the bins belongs to): vertex bins are set to unit capacity
(since there is only one vertex per bin), while each facet
bin is given a capacity equal to the ratio between its
area (i.e., the area of the associated centroidal Voronoi
cell) and the area of its containing facet. Finally, the
centroids of the bins in B are computed and stored as
representatives of their bins.
Fig. 3 Bins of a facet. Bins in a facet are defined as cells of
a centroidal Voronoi tessellation. Bin centroids are depicted
as red dots. Capacities of the bins (set proportional to their
areas) are depicted using a thermal color ramp.
3.2 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 {σ
j
}
j=1···L
and the cen-
troids of the bins in B as {b
j
}
j=1···M
, where L and M
are the number of simplices and bins respectively. The
capacity of bin b
j
is denoted c
j
. We also define s(j) to
be the index of the simplex containing the bin b
j
(i.e.,
b
j
∈ σ
s(j)
). Finally, we denote by m
ij
the amount of
mass transported from a given input point p
i
∈ S to
the centroid 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:
∀ij : m
ij
≥ 0, (1)
∀i :
X
j
m
ij
= m
i
, (2)
∀j
1
, j
2
s.t. s(j
1
)=s(j
2
) :
P
i
m
ij
1
c
j
1
=
P
i
m
ij
2
c
j
2
, (3)
where Equation 2 ensures that the entire measure of
an input point gets transported onto the mesh C, and
Equation 3 ensures a uniform measure over each facet
of C.
An optimal transport plan is then defined as a trans-
port plan π that minimizes the associated transport
cost:
cost(π) =
X
ij
m
ij
kp
i
− b
j
k
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, their positions,
as well as the square distances between input points and
bin centroids are all precomputed. In order to enforce
the uniformity constraint (Eq. 3) more sparsely, we also
introduce L additional variables l
i
(one per simplex σ
i
)
indicating the target measure density of the correspond-
ing simplex. The final problem formulation is thus:
Minimize
P
ij
m
ij
kp
i
− b
j
k
2
w.r.t. the variables m
ij
and l
s(j)
, and subject to:
∀i :
X
j
m
ij
= m
i
∀j :
X
i
m
ij
= c
j
· l
s(j)
∀i, j : m
ij
≥ 0, l
s(j)
≥ 0
3.3 Local Relaxation
Solving directly for the formulation described above is
compute-intensive due to the number of variables and
constraints involved: it requires instantiating a dense
matrix (representing the constraints) of size (M × N +
L) × (M + N ). For example, computing the optimal
transport cost between an input point set of 2, 100 sam-
ples and a simplicial complex containing 782 simplices
on which 7, 300 bins are placed involves solving for a
linear program of over 15 million variables and 9, 000
