Surface Reconstruction from Unorganized Points

Hugues Hoppe Tony DeRose Tom Duchampy

John McDonald z Werner Stuetzlez

University of Washington
Seattle, WA 98195

Abstract complex objects, including those as simple as a coffee cup

with a handle (a surface of genus 1), or the object depicted
We describe and demonstrate an algorithm that takes as input an in Figure 1a (a surface of genus 3), cannot be accomplished
f
unorganized set of points x1 : : : xn g IR3 on or near an un- by either of these methods. To adequately scan these objects,
known manifold M, and produces as output a simplicial surface that multiple view points must be used. Merging the data generated
approximates M. Neither the topology, the presence of boundaries, from multiple view points to reconstruct a polyhedral surface
nor the geometry of M are assumed to be known in advance — all representation is a non-trivial task [11].
are inferred automatically from the data. This problem naturally
arises in a variety of practical situations such as range scanning
an object from multiple view points, recovery of biological shapes
Surfaces from contours: In many medical studies it is com-
from two-dimensional slices, and interactive surface sketching. mon to slice biological specimens into thin layers with a mi-
crotome. The outlines of the structures of interest are then
digitized to create a stack of contours. The problem is to
CR Categories and Subject Descriptors: I.3.5 [Computer reconstruct the three-dimensional structures from the stacks
Graphics]: Computational Geometry and Object Modeling. of two-dimensional contours. Although this problem has re-
Additional Keywords: Geometric Modeling, Surface Fitting, ceived a good deal of attention, there remain severe limitations
Three-Dimensional Shape Recovery, Range Data Analysis. with current methods. Perhaps foremost among these is the
difficulty of automatically dealing with branching structures
[3, 12].
1 Introduction
Broadly speaking, the class of problems we are interested in can
Interactive surface sketching: A number of researchers, in-
be stated as follows: Given partial information of an unknown cluding Schneider [21] and Eisenman [6], have investigated
surface, construct, to the extent possible, a compact representation the creation of curves in IR2 by tracing the path of a stylus or
of the surface. Reconstruction problems of this sort occur in diverse mouse as the user sketches the desired shape. Sachs et al. [19]
scientific and engineering application domains, including: describe a system, called 3-Draw, that permits the creation of
free-form curves in IR3 by recording the motion of a stylus fitted

Surfaces from range data: The data produced by laser range with a Polhemus sensor. This can be extended to the design of
scanning systems is typically a rectangular grid of distances free-form surfaces by ignoring the order in which positions are
from the sensor to the object being scanned. If the sensor recorded, allowing the user to move the stylus arbitrarily back
and object are fixed, only objects that are “point viewable” and forth over the surface. The problem is then to construct
can be fully digitized. More sophisticated systems, such as a surface representation faithful to the unordered collection of
those produced by Cyberware Laboratory, Inc., are capable points.
of digitizing cylindrical objects by rotating either the sensor
or the object. However, the scanning of topologically more
Department of Computer Science and Engineering, FR-35 Reconstruction algorithms addressing these problems have typi-
y Department of Mathematics, GN-50 cally been crafted on a case by case basis to exploit partial structure
z Department of Statistics, GN-22 in the data. For instance, algorithms solving the surface from con-
This work was supported in part by Bellcore, the Xerox Corporation, tours problem make heavy use of the fact that data are organized into
IBM, Hewlett-Packard, the Digital Equipment Corporation, the Depart- contours (i.e., closed polygons), and that the contours lie in paral-
ment of Energy under grant DE-FG06-85-ER25006, the National Library of lel planes. Similarly, specialized algorithms to reconstruct surfaces
Medicine under grant NIH LM-04174, and the National Science Foundation from multiple view point range data might exploit the adjacency
under grants CCR-8957323 and DMS-9103002. relationship of the data points within each view.
In contrast, our approach is to pose a unifying general problem
that does not assume any structure on the data points. This approach
has both theoretical and practical merit. On the theoretical side,
abstracting to a general problem often sheds light on the truly critical
aspects of the problem. On the practical side, a single algorithm
that solves the general problem can be used to solve any specific
problem instance.
1.1 Terminology In contrast to implicit reconstruction techniques, parametric re-
construction techniques represent the reconstructed surface as a
By a surface we mean a “compact, connected, orientable two- topological embedding f ( ) of a 2-dimensional parameter domain
dimensional manifold, possibly with boundary, embedded in IR3 ” into IR3 . Previous work has concentrated on domain spaces with
(cf. O’Neill [17]). A surface without boundary will be called a simple topology, i.e. the plane and the sphere. Hastie and Stuet-
closed surface. If we want to emphasize that a surface possesses a zle [9] and Vemuri [26, 27] discuss reconstruction of surfaces by a
non-empty boundary, we will call it a bordered surface. A piecewise topological embedding f ( ) of a planar region into IR3 . Schudy
linear surface with triangular faces will be referred to as a simplicial
kk
surface. We use x to denote the Euclidean length of a vector x,
and Ballard [22, 23] and Brinkley [4] consider the reconstruction
of surfaces that are slightly deformed spheres, and thus choose
and we use d(X Y) to denote the Hausdorff distance between the to be a sphere. Sclaroff and Pentland [24] describe a hybrid im-
sets of points X and Y (the Hausdorff distance is simply the distance plicit/parametric method for fitting a deformed sphere to a set of
between the two closest points of X and Y). points using deformations of a superquadric.
f g
Let X = x1 : : : xn be sampled data points on or near an Compared to the techniques mentioned above, our method has
unknown surface M (see Figure 1b). To capture the error in most
sampling processes, we assume that each of the points xi 2 X is
several advantages:
2
of the form xi = yi + ei , where yi M is a point on the unknown
It requires only an unorganized collection of points on or near
2
surface and ei IR3 is an error vector. We call such a sample X the surface. No additional information is needed (such as
k k
 -noisy if ei  for all i. A value for  can be estimated in most normal information used by Muraki’s method).
applications (e.g., the accuracy of the laser scanner). Features of M
Unlike the parametric methods mentioned above, it can recon-
that are small compared to  will obviously not be recoverable. struct surfaces of arbitrary topology.
It is also impossible to recover features of M in regions where
insufficient sampling has occurred. In particular, if M is a bordered

Unlike previously suggested implicit methods, it deals with
boundaries in a natural way, and it does not generate spurious
surface, such as a sphere with a disc removed, it is impossible to surface components not supported by the data.
distinguish holes in the sample from holes in the surface. To capture
the intuitive notion of sampling density we need to make another
f
definition: Let Y = y1 : : : yn g M be a (noiseless) sample of a 2.2 Surface Reconstruction vs Function Recon-
surface M. The sample Y is said to be -dense if any sphere with struction
radius  and center in M contains at least one sample point in Y. A
f
 -noisy sample x1 : : : xn g IR3 of a surface M is said to be - Terms like “surface fitting” appear in reference to two distinct
f
dense if there exists a noiseless -dense sample y1 : : : yn g M classes of problems: surface reconstruction and function recon-
k k
such that xi = yi + ei , ei  , i = 1 : : : n.
struction. The goal of surface reconstruction was stated earlier. The
goal of function reconstruction may be stated as follows: Given a
f 2 g f 2 g
surface M, a set xi M , and a set yi IR , determine a function
1.2 Problem Statement f :M ! 
IR, such that f (xi ) yi .
The goal of surface reconstruction is to determine a surface M0 (see The domain surface M is most commonly a plane embedded in
Figure 2f) that approximates an unknown surface M (Figure 1a), IR3 , in which case the problem is a standard one considered in
using a sample X (Figure 1b) and information about the sampling approximation theory. The case where M is a sphere has also been
process, for example, bounds on the noise magnitude  and the extensively treated (cf. [7]). Some recent work under the title
sampling density . surfaces on surfaces addresses the case when M is a general curved
surface such as the skin of an airplane [16].
We are currently working to develop conditions on the original
surface M and the sample X that are sufficient to allow M to be Function reconstruction methods can be used for surface recon-
reliably reconstructed. As that work is still preliminary, we are un- struction in simple, special cases, where the surface to be recon-
able to give guarantees for the algorithm presented here. However, structed is, roughly speaking, the graph of a function over a known
the algorithm has worked well in practice where the results can be surface M. It is important to recognize just how limited these spe-
compared to the original surface (see Section 4). cial cases are — for example, not every surface homeomorphic to a
sphere is the graph of a function over the sphere. The point we want
to make is that function reconstruction must not be misconstrued to
2 Related Work solve the general surface reconstruction problem.

2.1 Surface Reconstruction 3 A Description of the Algorithm

Surface reconstruction methods can be classified according to the
way in which they represent the reconstructed surface. 3.1 Overview
Implicit reconstruction methods attempt to find a smooth func-
tion f : IR3 ! f g
IR such that x1 : : : xn is close to the zero set
Our surface reconstruction algorithm consists of two stages. In the
first stage we define a function f : D! IR, where D IR3 is a region
Z(f ). They differ with respect to the form of f and the measure of near the data, such that f estimates the signed geometric distance to
closeness. Pratt [18] and Taubin [25] minimize the sum of squared the unknown surface M. The zero set Z(f ) is our estimate for M.
Hausdorff distances from the data points to the zero set of a poly- In the second stage we use a contouring algorithm to approximate
nomial in three variables. Muraki [15] takes f to be a linear combi- Z(f ) by a simplicial surface.
nation of three-dimensional Gaussian kernels with different means
jj
Although the unsigned distance function f would be easier to
and spreads. His goodness-of-fit function measures how close the
values of f at the data points are to zero, and how well the unit jj
estimate, zero is not a regular value of f . Zero is, however, a regular
normals to the zero set of f match the normals estimated from the value of f , and the implicit function theorem thus guarantees that
data. Moore and Warren [13] fit a piecewise polynomial recursively our approximation Z(f ) is a manifold.
and then enforce continuity using a technique they call free form The key ingredient to defining the signed distance function is to
blending. associate an oriented plane with each of the data points. These
tangent planes serve as local linear approximations to the surface. the surface is assumed to consist of a single connected component,
Although the construction of the tangent planes is relatively simple, the graph should be connected. A simple connected graph for a set
the selection of their orientations so as to define a globally consistent of points that tends to connect neighbors is the Euclidean Minimum
orientation for the surface is one of the major obstacles facing the Spanning Tree (EMST). However, the EMST over the tangent plane
algorithm. As indicated in Figure 2b, the tangent planes do not f g
centers o1 : : : on (Figure 1c) is not sufficiently dense in edges
directly define the surface, since their union may have a complicated to serve our purposes. We therefore enrich it by adding a number
non-manifold structure. Rather, we use the tangent planes to define of edges to it. Specifically, we add the edge (i j) if either oi is
the signed distance function to the surface. An example of the in the k-neighborhood of oj , or oj is in the k-neighborhood of oi
simplicial surface obtained by contouring the zero set of the signed f
(where k-neighborhood is defined over o1 : : : on as it was for g
distance function is shown in Figure 2e. The next several sections X). The resulting graph (Figure 1d), called the Riemannian Graph,
develop in more detail the successive steps of the algorithm. is thus constructed to be a connected graph that encodes geometric
proximity of the tangent plane centers.
3.2 Tangent Plane Estimation A relatively simple-minded algorithm to orient the planes would
be to arbitrarily choose an orientation for some plane, then “propa-
The first step toward defining a signed distance function is to com- gate” the orientation to neighboring planes in the Riemannian Graph.
pute an oriented tangent plane for each data point. The tangent plane
Tp (xi) associated with the data point xi is represented as a point oi , In practice, we found that the order in which the orientation is prop-
called the center, together with a unit normal vector n ^ i . The signed agated is important. Figure 3b shows what may result when prop-
distance of an arbitrary point p 2
IR3 to Tp (xi ) is defined to be
agating orientation solely on the basis of geometric proximity; a
; 
disti (p) = (p oi ) n^ i . The center and normal for Tp (xi) are deter- correct reconstruction is shown in Figure 3c. Intuitively, we would
like to choose an order of propagation that favors propagation from
mined by gathering together the k points of X nearest to xi ; this set Tp (xi) to Tp (xj) if the unoriented planes are nearly parallel. This
is denoted by Nbhd (xi ) and is called the k-neighborhood of xi . (We can be accomplished by assigning to each edge (i j) in the Rieman-
currently assume k to be a user-specified parameter, although in Sec-
tion 5 we propose a method for determining k automatically.) The
nian Graph the cost 1 n ;j  j
^ i n^ j . In addition to being non-negative,
f
center and unit normal are computed so that the plane disti (p) = 0 g this assignment has the property that a cost is small if the unoriented
is the least squares best fitting plane to Nbhd (xi ). That is, the cen-
tangent planes are nearly parallel. A favorable propagation order
ter oi is taken to be the centroid of Nbhd (xi ), and the normal n ^i can therefore be achieved by traversing the minimal spanning tree
is determined using principal component analysis. To compute n ^i, (MST) of the resulting graph. This order is advantageous because it
the covariance matrix of Nbhd (xi ) is formed. This is the symmetric
tends to propagate orientation along directions of low curvature in

X
the data, thereby largely avoiding ambiguous situations encountered
3 3 positive semi-definite matrix when trying to propagate orientation across sharp edges (as at the tip
of the cat’s ears in Figure 3b). In the MST shown in Figure 2a, the
CV = (y ; oi ) (y ; oi ) edges are colored according to their cost, with the brightly colored
y2Nbhd (xi ) edges corresponding to regions of high variation (where n 
^ i n^ j is
somewhat less than 1).
where denotes the outer product vector operator1 . If 1i 2
i To assign orientation to an initial plane, the unit normal of the
3
i denote the eigenvalues of CV associated with unit eigenvectors plane whose center has the largest z coordinate is forced to point
v^ i1 v^ i2 v^ i3 , respectively, we choose n ;
^ i to be either v^ i3 or v^ i3. The toward the +z axis. Then, rooting the tree at this initial node,
selection determines the orientation of the tangent plane, and it must we traverse the tree in depth-first order, assigning each plane an
be done so that nearby planes are “consistently oriented”. orientation that is consistent with that of its parent. That is, if
during traversal, the current plane Tp (xi ) has been assigned the
3.3 Consistent Tangent Plane Orientation orientation n^ i and Tp (xj) is the next plane to be visited, then n^ j is
replaced with n ; 
^ j if n^ i n^ j < 0.
2
Suppose two data points xi xj X are geometrically close. Ideally,
This orientation algorithm has been used in all our examples
when the data is dense and the surface is smooth, the corresponding
tangent planes Tp (xi ) = (oi n^ i ) and Tp (xj) = (oj n^ j ) are nearly and has produced correct orientations in all the cases we have run.
parallel, i.e. n  
^ i n^ j 1. If the planes are consistently oriented,
The resulting oriented tangent planes are represented as shaded
then n  
^ i n^ j +1; otherwise, either n^ i or n^ j should be flipped. rectangles in Figure 2b.
The difficulty in finding a consistent global orientation is that this
condition should hold between all pairs of “sufficiently close” data 3.4 Signed Distance Function
points. 2
The signed distance f (p) from an arbitrary point p IR3 to a known
We can model the problem as graph optimization. The graph surface M is the distance between p and the closest point z M, 2
contains one node Ni per tangent plane Tp (xi ), with an edge (i j) multiplied by 1, depending on which side of the surface p lies.
between Ni and Nj if the tangent plane centers oi and oj are suf- In reality M is not known, but we can mimic this procedure using
ficiently close (we will be more precise about what we mean by the oriented tangent planes as follows. First, we find the tangent
sufficiently close shortly). The cost on edge (i j) encodes the de- plane Tp (xi ) whose center oi is closest to p. This tangent plane is
gree to which Ni and Nj are consistently oriented and is taken to be a local linear approximation to M, so we take the signed distance
n 
^ i n^ j . The problem is then to select orientations for the tangent f (p) to M to be the signed distance between p and its projection z
planes so as to maximize the total cost of the graph. Unfortunately, onto Tp (xi ); that is,
this problem can be shown to be NP-hard via a reduction to MAX-
CUT [8]. To efficiently solve the orientation problem we must f (p) = disti (p) = (p ; oi )  n^ i :

therefore resort to an approximation algorithm.

Before describing the approximation algorithm we use, we must If M is known not to have boundaries, this simple rule works
decide when a pair of nodes are to be connected in the graph. Since well. However, the rule must be extended to accommodate surfaces
that might have boundaries. Recall that the set X = x1 : : : xn f g
1 If a
and b have components ai and bj respectively, then the matrix is assumed to be a -dense,  -noisy sample of M. If there was no
a
b has ai bj as its ij-th entry. noise, we could deduce that a point z with d(z X) >  cannot be
a point of M since that would violate X being -dense. Intuitively, 4 Results
the sample points do not leave holes of radius larger than . If
the sample is  -noisy, the radius of the holes may increase, but by We have experimented with the reconstruction method on data sets
no more than  . We therefore conclude that a point z cannot be obtained from several different sources. In all cases, any structure
a point of M if d(z X) >  +  . If the projection z of p onto (including ordering) that might have been present in the point sets
the closest tangent plane has d(z X) >  +  , we take f (p) to be was discarded.
undefined. Undefined values are used by the contouring algorithm Meshes : Points were randomly sampled from a number of existing
of Section 3.5 to identify boundaries. simplicial surfaces3 . For instance, the mesh of Figure 3a was
Stated procedurally, our signed distance function is defined as: randomly sampled to yield 1000 unorganized points, and these
in turn were used to reconstruct the surface in Figure 3c. This
i index of tangent plane whose center is closest to p particular case illustrates the behavior of the method on a bor-
f Compute z as the projection of p onto Tp(xi ) g
dered surface (the cat has no base and is thus homeomorphic
z oi ; ((p ; oi )  n
to a disc). The reconstructed knot (original mesh from Rob
^ i ) n^ i Scharein) of Figure 3d is an example of a surface with simple
topology yet complex geometrical embedding.
if d(z X) <  +  then
f (p) ;
(p oi ) n ^i  f= kp ; zkg Ray Traced Points : To simulate laser range imaging from mul-
tiple view points, CSG models were ray traced from multiple
else
f (p) undefined eye points. The ray tracer recorded the point of first intersec-
endif tion along each ray. Eight eye points (the vertices of a large
cube centered at the object) were used to generate the point set
The simple approach outlined above creates a zero set Z(f ) that of Figure 1b from the CSG object shown in Figure 1a. This
is piecewise linear but contains discontinuities. The discontinuities is the point set used in Section 3 to illustrate the steps of the
result from the implicit partitioning of space into regions within algorithm (Figures 1a-2f).
which a single tangent plane is used to define the signed distance Range Images : The bust of Spock (Figure 3e) was reconstructed
function. (These regions are in fact the Voronoi regions associated from points taken from an actual cylindrical range image (gen-
with the centers oi .) Fortunately, the discontinuities do not adversely erated by Cyberware Laboratory, Inc.). Only 25% of the orig-
affect our algorithm. The contouring algorithm discussed in the inal points were used.
next section will discretely sample the function f over a portion
Contours : Points from 39 planar (horizontal) slices of the CT
of a 3-dimensional grid near the data and reconstruct a continuous
scan of a femur were combined together to obtain the surface
piecewise linear approximation to Z(f ).
of Figure 3f.
The algorithm’s parameters are shown in the next table for each
3.5 Contour Tracing of the examples. The execution times were obtained on a 20 MIPS
workstation. The parameter  +  and the marching cube cell size
Contour tracing, the extraction of an isosurface from a scalar func- are both expressed as a fraction of the object’s size. The parameter
tion, is a well-studied problem [1, 5, 28]. We chose to implement  +  is set to infinity for those surfaces that are known to be closed.

a variation of the marching cubes algorithm (cf. [1]) that samples

the function at the vertices of a cubical lattice and finds the contour Object n k  + cell size time
intersections within tetrahedral decompositions of the cubical cells. (seconds)
cat 1000 15 .06 1=30 19
To accurately estimate boundaries, the cube size should be set so knot 10000 20 1 1=50 137
that edges are of length less than  +  . In practice we have often
found it convenient to set the cube size somewhat larger than this
mechpart 4102 12 1 1=40 54
spock 21760 8 .08 1=80 514
value, simply to increase the speed of execution and to reduce the femur 18224 40 .06 1=50 2135
number of triangular facets generated.
The algorithm only visits cubes that intersect the zero set by push-
ing onto a queue only the appropriate neighboring cubes (Figure 2c).
5 Discussion
In this way, the signed distance function f is evaluated only at points
close to the data. Figure 2d illustrates the signed distance function 5.1 Tangent Plane Approximation
by showing line segments between the query points p (at the cube The neighborhood Nbhd (xi ) of a data point xi is defined to consist
vertices) and their associated projected points z. As suggested in of its k nearest neighbors, where k is currently assumed to be an in-
Section 3.4, no intersection is reported within a cube if the signed put parameter. In the case where the data contains little or no noise,
distance function is undefined at any vertex of the cube, thereby k is not a critical parameter since the output has been empirically
giving rise to boundaries in the simplicial surface. observed to be stable over a wide range of settings. However, it
The resulting simplicial surface can contain triangles with ar- would be best if k could be selected automatically. Furthermore, al-
bitrarily poor aspect ratio (Figure 2e). We alleviate this problem lowing k to adapt locally would make less stringent the requirement
using a post-processing procedure that collapses edges in the sur- that the data be uniformly distributed over the surface. To select
face using an aspect ratio criterion.2 The final result is shown in and adapt k, the algorithm could incrementally gather points while
Figure 2f. Alternatively, other contouring methods exist that can monitoring the changing eigenvalues of the covariance matrix (see
guarantee bounds on the triangle aspect ratio [14]. Section 3.2). For small values of k, data noise tends to dominate,
the eigenvalues are similar, and the eigenvectors do not reveal the
surface’s true tangent plane. At the other extreme, as k becomes
2 The edges are kept in a priority queue; the criterion to minimize is
the product of the edge length times the minimum inscribed radius of its 3 Discrete inverse transform sampling [10, page 469] on triangle area was
two adjacent faces. Tests are also performed to ensure that edge collapses used to select face indices from the mesh, and uniform sampling was used
preserve the topological type of the surface. within the faces.
large, the k-neighborhoods become less localized and the surface on the sample and the original surface. To further improve the
curvature tends to increase the “thickness” 3i of the neighborhood. geometric accuracy of the fit, and to reduce the space required
Another possible criterion is to compare 3i to some local or global to store the reconstruction, we envision using the output of our
estimate of data noise. Although we have done some initial exper- algorithm as the starting point for a subsequent spline surface fitting
imentation in this direction, we have not yet fully examined these procedure. We are currently investigating such a method based on
options. a nonlinear least squares approach using triangular Bézier surfaces.
If the data is obtained from range images, there exists some
knowledge of surface orientation at each data point. Indeed, each
data point is known to be visible from a particular viewing direction,
References
so that, unless the surface incident angle is large, the point’s tangent [1] E. L. Allgower and P. H. Schmidt. An algorithm for piecewise
plane orientation can be inferred from that viewing direction. Our linear approximation of an implicitly defined manifold. SIAM
method could exploit this additional information in the tangent plane Journal of Numerical Analysis, 22:322–346, April 1985.
orientation step (Section 3.3) by augmenting the Riemannian Graph
with an additional pseudo-node and n additional edges. [2] J. L. Bentley. Multidimensional divide and conquer. Comm.
ACM, 23(4):214–229, 1980.
5.2 Algorithm Complexity [3] Y. Breseler, J. A. Fessler, and A. Macovski. A Bayesian
approach to reconstruction from incomplete projections of a
A spatial partitioning Abstract Data Type greatly improves per- multiple object 3D domain. IEEE Trans. Pat. Anal. Mach.
formance of many of the subproblems discussed previously. The Intell., 11(8):840–858, August 1989.
critical subproblems are (with their standard time complexity):
[4] James F. Brinkley. Knowledge-driven ultrasonic three-

EMST graph (O(n2 )) dimensional organ modeling. IEEE Trans. Pat. Anal. Mach.

k-nearest neighbors to a given point (O(n + k log n)) Intell., 7(4):431–441, July 1985.

nearest tangent plane origin to a given point (O(n)) [5] David P. Dobkin, Silvio V. F. Levy, William P. Thurston, and
Allan R. Wilks. Contour tracing by piecewise linear approxi-
Hierarchical spatial partitioning schemes such as octrees [20] mations. ACM TOG, 9(4):389–423, October 1990.
and k-D trees [2] can be used to solve these problems more ef-
ficiently. However, the uniform sampling density assumed in our [6] John A. Eisenman. Graphical editing of composite bezier
data allows simple spatial cubic partitioning to work efficiently. The curves. Master’s thesis, Department of Electrical Engineering
axis-aligned bounding box of the points is partitioned by a cubical and Computer Science, M.I.T., 1988.
grid. Points are entered into sets corresponding to the cube to which [7] T.A. Foley. Interpolation to scattered data on a spherical do-
they belong, and these sets are accessed through a hash table in- main. In M. Cox and J. Mason, editors, Algorithms for Ap-
dexed by the cube indices. It is difficult to analyze the resulting proximation II, pages 303–310. Chapman and Hall, London,
improvements analytically, but, empirically, the time complexity of 1990.
the above problems is effectively reduced by a factor of n, except
[8] Michael R. Garey and David S. Johnson. Computers and
for the k-nearest neighbors problem which becomes O(k).
Intractability. W. H. Freeman and Company, 1979.
As a result of the spatial partitioning, the Riemannian Graph can
be constructed in O(nk) time. Because the Riemannian Graph has [9] T. Hastie and W. Stuetzle. Principal curves. JASA, 84:502–
O(n) edges (at most n+nk), the MST computation used in finding the 516, 1989.
best path on which to propagate orientation requires only O(n log n) [10] Averill M. Law and W. David Kelton. Simulation Modeling
time. Traversal of the MST is of course O(n). and Analysis. McGraw-Hill, Inc., second edition, 1991.
The time complexity of the contouring algorithm depends only [11] Marshal L. Merriam. Experience with the cyberware 3D dig-
on the number of cubes visited, since the evaluation of the signed itizer. In NCGA Proceedings, pages 125–133, March 1992.
distance function f at a point p can be done in constant time (the
closest tangent plane origin oi to p and the closest data point xj to [12] David Meyers, Shelly Skinner, and Kenneth Sloan. Surfaces
the projected point z can both be found in constant time with spatial from contours: The correspondence and branching problems.
partitioning). In Proceedings of Graphics Interface ’91, pages 246–254,
June 1991.
[13] Doug Moore and Joe Warren. Approximation of dense scat-
6 Conclusions and Future Work tered data using algebraic surfaces. TR 90-135, Rice Univer-
sity, October 1990.
We have developed an algorithm to reconstruct a surface in three-
dimensional space with or without boundary from a set of unorga- [14] Doug Moore and Joe Warren. Adaptive mesh generation ii:
nized points scattered on or near the surface. The algorithm, based Packing solids. TR 90-139, Rice University, March 1991.
on the idea of determining the zero set of an estimated signed dis- [15] Shigeru Muraki. Volumetric shape description of range data
tance function, was demonstrated on data gathered from a variety using “blobby model”. Computer Graphics (SIGGRAPH ’91
of sources. It is capable of automatically inferring the topological Proceedings), 25(4):227–235, July 1991.
type of the surface, including the presence of boundary curves.
[16] Gregory M. Nielson, Thomas A. Foley, Bernd Hamann, and
The algorithm can, in principle, be extended to reconstruct mani- David Lane. Visualizing and modeling scattered multivariate
folds of co-dimension one in spaces of arbitrary dimension; that is,
;
to reconstruct d 1 dimensional manifolds in d dimensional space.
data. IEEE CG&A, 11(3):47–55, May 1991.
Thus, essentially the same algorithm can be used to reconstruct [17] Barrett O’Neill. Elementary Differential Geometry. Academic
curves in the plane or volumes in four-dimensional space. Press, Orlando, Florida, 1966.
The output of our reconstruction method produced the correct [18] Vaughan Pratt. Direct least-squares fitting of algebraic sur-
topology in all the examples. We are trying to develop formal faces. Computer Graphics (SIGGRAPH ’87 Proceedings),
guarantees on the correctness of the reconstruction, given constraints 21(4):145–152, July 1987.
 (a) Original CSG object (b) Sampled points (xi ) (n = 4102)

(c) EMST of tangent plane centers oi (d) Riemannian Graph over oi

Figure 1: Reconstruction of ray-traced CSG object (simulated multi-view range data).

[19] Emanuel Sachs, Andrew Roberts, and David Stoops. 3-Draw: [24] Stan Sclaroff and Alex Pentland. Generalized implicit func-
A tool for designing 3D shapes. IEEE Computer Graphics tions for computer graphics. Computer Graphics (SIGGRAPH
and Applications, 11(6):18–26, November 1991. ’91 Proceedings), 25(4):247–250, July 1991.
[20] Hanan Samet. Applications of Spatial Data Structures. [25] G. Taubin. Estimation of planar curves, surfaces and nonpla-
Addison-Wesley, 1990. nar space curves defined by implicit equations, with applica-
tions to edge and range image segmentation. Technical Report
[21] Philip J. Schneider. Phoenix: An interactive curve design sys- LEMS-66, Division of Engineering, Brown University, 1990.
tem based on the automatic fitting of hand-sketched curves.
Master’s thesis, Department of Computer Science, U. of Wash- [26] B. C. Vemuri. Representation and Recognition of Objects
ington, 1988. From Dense Range Maps. PhD thesis, Department of Electri-
cal and Computer Engineering, University of Texas at Austin,
[22] R. B. Schudy and D. H. Ballard. Model detection of cardiac 1987.
chambers in ultrasound images. Technical Report 12, Com- [27] B. C. Vemuri, A. Mitiche, and J. K. Aggarwal. Curvature-
puter Science Department, University of Rochester, 1978. based representation of objects from range data. Image and
[23] R. B. Schudy and D. H. Ballard. Towards an anatomical model Vision Computing, 4(2):107–114, 1986.
of heart motion as seen in 4-d cardiac ultrasound data. In [28] G. Wyvill, C. McPheeters, and B. Wyvill. Data structures
Proceedings of the 6th Conference on Computer Applications for soft objects. The Visual Computer, 2(4):227–234, August
in Radiology and Computer-Aided Analysis of Radiological 1986.
Images, 1979.
(a) Traversal order of orientation propagation (b) Oriented tangent planes (Tp (xi ))

(c) Cubes visited during contouring (d) Estimated signed distance (shown as p ; z)

(e) Output of modified marching cubes (f) Final surface after edge collapses
Figure 2: Reconstruction of ray-traced CSG object (continued).
(a) Original mesh (b) Result of naive orientation propagation

(c) Reconstructed bordered surface (d) Reconstructed surface with complex geometry

(e) Reconstruction from cylindrical range data (f) Reconstruction from contour data
Figure 3: Reconstruction examples.
Eurographics Symposium on Geometry Processing (2006)
Konrad Polthier, Alla Sheffer (Editors)

Poisson Surface Reconstruction

Michael Kazhdan1 , Matthew Bolitho1 and Hugues Hoppe2

1 Johns Hopkins University, Baltimore MD, USA

2 Microsoft Research, Redmond WA, USA

Abstract

We show that surface reconstruction from oriented points can be cast as a spatial Poisson problem. This Poisson
formulation considers all the points at once, without resorting to heuristic spatial partitioning or blending, and
is therefore highly resilient to data noise. Unlike radial basis function schemes, our Poisson approach allows a
hierarchy of locally supported basis functions, and therefore the solution reduces to a well conditioned sparse
linear system. We describe a spatially adaptive multiscale algorithm whose time and space complexities are pro-
portional to the size of the reconstructed model. Experimenting with publicly available scan data, we demonstrate
reconstruction of surfaces with greater detail than previously achievable.

0
1. Introduction 0 0 0
0
0 1
Reconstructing 3D surfaces from point samples is a well 0
0 1 1
studied problem in computer graphics. It allows fitting of 1
0
0 1
scanned data, filling of surface holes, and remeshing of ex- 0 0
1
isting models. We provide a novel approach that expresses 0 0
surface reconstruction as the solution to a Poisson equation. Oriented
G points Indicator gradient Indicator function Surface

Like much previous work (Section 2), we approach the V ’FM FM wM

problem of surface reconstruction using an implicit function Figure 1: Intuitive illustration of Poisson reconstruction in 2D.
framework. Specifically, like [Kaz05] we compute a 3D in-
dicator function χ (defined as 1 at points inside the model, tion χ whose Laplacian (divergence of gradient) equals the
and 0 at points outside), and then obtain the reconstructed divergence of the vector field
V ,
surface by extracting an appropriate isosurface.
∆χ ≡ ∇ · ∇χ = ∇ ·
V .
Our key insight is that there is an integral relationship be-
tween oriented points sampled from the surface of a model We will make these definitions precise in Sections 3 and 4.
and the indicator function of the model. Specifically, the gra- Formulating surface reconstruction as a Poisson problem
dient of the indicator function is a vector field that is zero offers a number of advantages. Many implicit surface fitting
almost everywhere (since the indicator function is constant methods segment the data into regions for local fitting, and
almost everywhere), except at points near the surface, where further combine these local approximations using blending
it is equal to the inward surface normal. Thus, the oriented functions. In contrast, Poisson reconstruction is a global so-
point samples can be viewed as samples of the gradient of lution that considers all the data at once, without resorting
the model’s indicator function (Figure 1). to heuristic partitioning or blending. Thus, like radial basis
The problem of computing the indicator function thus re- function (RBF) approaches, Poisson reconstruction creates
duces to inverting the gradient operator, i.e. finding the scalar very smooth surfaces that robustly approximate noisy data.
function χ whose gradient best approximates a vector field But, whereas ideal RBFs are globally supported and non-

V defined by the samples, i.e. minχ
∇χ −
V
. If we apply decaying, the Poisson problem admits a hierarchy of locally
the divergence operator, this variational problem transforms supported functions, and therefore its solution reduces to a
into a standard Poisson problem: compute the scalar func- well-conditioned sparse linear system.

c The Eurographics Association 2006.
 Kazhdan et al. / Poisson Surface Reconstruction

Moreover, in many implicit fitting schemes, the value Local fitting methods consider subsets of nearby points at
of the implicit function is constrained only near the sam- a time. A simple scheme is to estimate tangent planes and
ple points, and consequently the reconstruction may con- define the implicit function as the signed distance to the tan-
tain spurious surface sheets away from these samples. Typ- gent plane of the closest point [HDD∗ 92]. Signed distance
ically this problem is attenuated by introducing auxiliary can also be accumulated into a volumetric grid [CL96]. For
“off-surface” points (e.g. [CBC∗ 01, OBA∗ 03]). With Pois- function continuity, the influence of several nearby points
son surface reconstruction, such surface sheets seldom arise can be blended together, for instance using moving least
because the gradient of the implicit function is constrained at squares [ABCO∗ 01,SOS04]. A different approach is to form
all spatial points. In particular it is constrained to zero away point neighborhoods by adaptively subdividing space, for
from the samples. example with an adaptive octree. Blending is possible over
an octree structure using a multilevel partition of unity, and
Poisson systems are well known for their resilience in the
the type of local implicit patch within each octree node can
presence of imperfect data. For instance, “gradient domain”
be selected heuristically [OBA∗ 03].
manipulation algorithms (e.g. [FLW02]) intentionally mod-
ify the gradient data such that it no longer corresponds to any Our Poisson reconstruction combines benefits of both
real potential field, and rely on a Poisson system to recover global and local fitting schemes. It is global and therefore
the globally best-fitting model. does not involve heuristic decisions for forming local neigh-
There has been broad interdisciplinary research on solv- borhoods, selecting surface patch types, and choosing blend
ing Poisson problems and many efficient and robust methods weights. Yet, the basis functions are associated with the am-
have been developed. One particular aspect of our problem bient space rather than the data points, are locally supported,
instance is that an accurate solution to the Poisson equation and have a simple hierarchical structure that results in a
is only necessary near the reconstructed surface. This allows sparse, well-conditioned system.
us to leverage adaptive Poisson solvers to develop a recon- Our approach of solving an indicator function is sim-
struction algorithm whose spatial and temporal complexities ilar to the Fourier-based reconstruction scheme of Kazh-
are proportional to the size of the reconstructed surface. dan [Kaz05]. In fact, we show in Appendix A that our basic
Poisson formulation is mathematically equivalent. Indeed,
2. Related Work the Fast Fourier Transform (FFT) is a common technique
for solving dense, periodic Poisson systems. However, the
Surface reconstruction The reconstruction of surfaces FFT requires O(r3 log r) time and O(r3 ) space where r is
from oriented points has a number of difficulties in prac- the 3D grid resolution, quickly becoming prohibitive for fine
tice. The point sampling is often nonuniform. The positions resolutions. In contrast, the Poisson system allows adaptive
and normals are generally noisy due to sampling inaccuracy discretization, and thus yields a scalable solution.
and scan misregistration. And, accessibility constraints dur-
ing scanning may leave some surface regions devoid of data. Poisson problems The Poisson equation arises in numer-
Given these challenges, reconstruction methods attempt to ous applications areas. For instance, in computer graph-
infer the topology of the unknown surface, accurately fit (but ics it is used for tone mapping of high dynamic range im-
not overfit) the noisy data, and fill holes reasonably. ages [FLW02], seamless editing of image regions [PGB03],
Several approaches are based on combinatorial structures, fluid mechanics [LGF04], and mesh editing [YZX∗ 04].
such as Delaunay triangulations (e.g. [Boi84, KSO04]), al- Multigrid Poisson solutions have even been adapted for effi-
pha shapes [EM94, BBX95, BMR∗ 99]), or Voronoi dia- cient GPU computation [BFGS03, GWL∗ 03].
grams [ABK98, ACK01]. These schemes typically create a The Poisson equation is also used in heat transfer and
triangle mesh that interpolates all or a most of the points. diffusion problems. Interestingly, Davis et al [DMGL02]
In the presence of noisy data, the resulting surface is often use diffusion to fill holes in reconstructed surfaces. Given
jagged, and is therefore smoothed (e.g. [KSO04]) or refit to boundary conditions in the form of a clamped signed dis-
the points (e.g. [BBX95]) in subsequent processing. tance function d, their diffusion approach essentially solves
Other schemes directly reconstruct an approximating sur- the homogeneous Poisson equation ∆d = 0 to create an im-
face, typically represented in implicit form. We can broadly plicit surface spanning the boundaries. They use a local iter-
classify these as either global or local approaches. ative solution rather than a global multiscale Poisson system.
Global fitting methods commonly define the implicit Nehab et al [NRDR05] use a Poisson system to fit a 2.5D
function as the sum of radial basis functions (RBFs) centered height field to sampled positions and normals. Their ap-
at the points (e.g. [Mur91, CBC∗ 01, TO02]). However, the proach fits a given parametric surface and is well-suited to
ideal RBFs (polyharmonics) are globally supported and non- the reconstruction of surfaces from individual scans. How-
decaying, so the solution matrix is dense and ill-conditioned. ever, in the case that the samples are obtained from the union
Practical solutions on large datasets involve adaptive RBF of multiple scans, their approach cannot be directly applied
reduction and the fast multipole method [CBC∗ 01]. to obtain a connected, watertight surface.

c The Eurographics Association 2006.
 Kazhdan et al. / Poisson Surface Reconstruction

3. Our Poisson reconstruction approach

The input data S is a set of samples s ∈ S, each consisting of a
point s.p and an inward-facing normal s.
N, assumed to lie on
or near the surface ∂ M of an unknown model M. Our goal is
to reconstruct a watertight, triangulated approximation to the
surface by approximating the indicator function of the model
and extracting the isosurface, as illustrated in Figure 2.
The key challenge is to accurately compute the indicator Figure 2: Points from scans of the “Armadillo Man” model (left),
function from the samples. In this section, we derive a rela- our Poisson surface reconstruction (right), and a visualization of the
tionship between the gradient of the indicator function and indicator function (middle) along a plane through the 3D volume.
an integral of the surface normal field. We then approximate surface geometry. However, the input set of oriented points
this surface integral by a summation over the given oriented provides precisely enough information to approximate the
point samples. Finally, we reconstruct the indicator function integral with a discrete summation. Specifically, using the
from this gradient field as a Poisson problem. point set S to partition ∂ M into distinct patches Ps ⊂ ∂ M,
we can approximate the integral over a patch Ps by the value
Defining the gradient field Because the indicator function at point sample s.p, scaled by the area of the patch:
is a piecewise constant function, explicit computation of its 
gradient field would result in a vector field with unbounded ∇(χM ∗ F̃)(q) = ∑ F̃p (q)
N∂ M (p)d p
values at the surface boundary. To avoid this, we convolve s∈S Ps
(2)
the indicator function with a smoothing filter and consider
the gradient field of the smoothed function. The following
≈ ∑ |Ps | F̃s.p (q) s.
N ≡
V (q).
s∈S
lemma formalizes the relationship between the gradient of
It should be noted that though Equation 1 is true for any
the smoothed indicator function and the surface normal field.
smoothing filter F̃, in practice, care must be taken in choos-
Lemma: Given a solid M with boundary ∂ M, let χM de- ing the filter. In particular, we would like the filter to satisfy
note the indicator function of M,
N∂ M (p) be the inward two conditions. On the one hand, it should be sufficiently
surface normal at p ∈ ∂ M, F̃(q) be a smoothing filter, and narrow so that we do not over-smooth the data. And on the
F̃p (q) = F̃(q− p) its translation to the point p. The gradient other hand, it should be wide enough so that the integral over
of the smoothed indicator function is equal to the vector field Ps is well approximated by the value at s.p scaled by the
obtained by smoothing the surface normal field: patch area. A good choice of filter that balances these two

  requirements is a Gaussian whose variance is on the order of
∇ χM ∗ F̃ (q0 ) = F̃p (q0 )
N∂ M (p)d p. (1) the sampling resolution.
∂M

Proof: To prove this, we show equality for each of the com- Solving the Poisson problem Having formed a vector field
ponents of the vector field. Computing the partial derivative
V , we want to solve for the function χ̃ such that ∇χ̃ =
V .
of the smoothed indicator function with respect to x, we get: However,
V is generally not integrable (i.e. it is not curl-
   free), so an exact solution does not generally exist. To find
∂ 
 ∂  the best least-squares approximate solution, we apply the di-
χ ∗ F̃ = F̃(q − p)d p
∂ x q0 ∂ x q=q0 M
M
vergence operator to form the Poisson equation
  
∂ ∆χ̃ = ∇ ·
V .
= − F̃(q0 − p) d p
M ∂x

 In the next section, we describe our implementation of
= − ∇ · F̃(q0 − p), 0, 0 d p
M these steps in more detail.
 
 
= F̃p (q0 ), 0, 0 ,
N∂ M (p) d p.
∂M
4. Implementation
(The first equality follows from the fact that χM is equal to
zero outside of M and one inside. The second follows from We first present our reconstruction algorithm under the as-
the fact that (∂ /∂ q)F̃(q − p) = −(∂ /∂ p)F̃(q − p). The last sumption that the point samples are uniformly distributed
follows from the Divergence Theorem.) over the model surface. We define a space of functions with
high resolution near the surface of the model and coarser
A similar argument shows that the y-, and z-components resolution away from it, express the vector field
V as a linear
of the two sides are equal, thereby completing the proof.
sum of functions in this space, set up and solve the Poisson
equation, and extract an isosurface of the resulting indicator
Approximating the gradient field Of course, we cannot function. We then extend our algorithm to address the case
evaluate the surface integral since we do not yet know the of non-uniformly sampled points.

c The Eurographics Association 2006.
 Kazhdan et al. / Poisson Surface Reconstruction

4.1. Problem Discretization node function. Since the sampling width is 2−D and the sam-
ples all fall into leaf nodes of depth D, the error arising from
First, we must choose the space of functions in which to dis-
the clamping can never be too big (at most, on the order of
cretize the problem. The most straightforward approach is
half the sampling width). In the next section, we show how
to start with a regular 3D grid [Kaz05], but such a uniform
the error can be further reduced by using trilinear interpola-
structure becomes impractical for fine-detail reconstruction,
tion to allow for sub-node precision.
since the dimension of the space is cubic in the resolution
while the number of surface triangles grows quadratically. Finally, since a maximum tree depth of D corresponds to a
Fortunately, an accurate representation of the implicit sampling width of 2−D , the smoothing filter should approxi-
function is only necessary near the reconstructed surface. mate a Gaussian with variance on the order of 2−D . Thus, F
This motivates the use of an adaptive octree both to repre- should approximate a Gaussian with unit-variance.
sent the implicit function and to solve the Poisson system For efficiency, we approximate the unit-variance Gaussian
(e.g. [GKS02, LGF04]). Specifically, we use the positions of by a compactly supported function so that (1) the resulting
the sample points to define an octree O and associate a func- Divergence and Laplacian operators are sparse and (2) the
tion Fo to each node o ∈ O of the tree, choosing the tree and evaluation of a function expressed as the linear sum of Fo at
the functions so that the following conditions are satisfied: some point q only requires summing over the nodes o ∈ O
1. The vector field
V can be precisely and efficiently repre- that are close to q. Thus, we set F to be the n-th convolution
sented as the linear sum of the Fo . of a box filter with itself resulting in the base function F:
2. The matrix representation of the Poisson equation, ex- 1 |t| < 0.5
pressed in terms of the Fo can be solved efficiently. F(x, y, z) ≡ (B(x)B(y)B(z))∗n with B(t) =
0 otherwise
3. A representation of the indicator function as the sum of
the Fo can be precisely and efficiently evaluated near the Note that as n is increased, F more closely approximates
surface of the model. a Gaussian and its support grows larger; in our implemen-
tation we use a piecewise quadratic approximation with
Defining the function space Given a set of point samples n = 3. Therefore, the function F is supported on the domain
S and a maximum tree depth D, we define the octree O to be [-1.5, 1.5]3 and, for the basis function of any octree node,
the minimal octree with the property that every point sample there are at most 53 -1 = 124 other nodes at the same depth
falls into a leaf node at depth D. whose functions overlap with it.

Next, we define a space of functions obtained as the span

of translates and scales of a fixed, unit-integral, base func- 4.2. Vector Field Definition
tion F : R3 → R. For every node o ∈ O, we set Fo to be the
To allow for sub-node precision, we avoid clamping a sam-
unit-integral “node function” centered about the node o and
ple’s position to the center of the containing leaf node and
stretched by the size of o:
  instead use trilinear interpolation to distribute the sample
q − o.c 1 across the eight nearest nodes. Thus, we define our approxi-
Fo (q) ≡ F .
o.w o.w3 mation to the gradient field of the indicator function as:
where o.c and o.w are the center and width of node o.
V (q) ≡ αo,s Fo (q)s.
N
∑ ∑ (3)
This space of functions FO,F ≡ Span{Fo } has a multires- s∈S o∈Ngbr (s)
D

olution structure similar to that of traditional wavelet repre- where NgbrD (s) are the eight depth-D nodes closest to s.p
sentations. Finer nodes are associated with higher-frequency and {αo,s } are the trilinear interpolation weights. (If the
functions, and the function representation becomes more neighbors are not in the tree, we refine it to include them.)
precise as we near the surface.
Since the samples are uniform, we can assume that the
Selecting a base function In selecting a base function F, area of a patch Ps is constant and
V is a good approxima-
our goal is to choose a function so that the vector field
V , tion, up to a multiplicative constant, of the gradient of the
defined in Equation 2, can be precisely and efficiently repre- smoothed indicator function. We will show that the choice
sented as the linear sum of the node functions {Fo }. of multiplicative constant does not affect the reconstruction.
If we were to replace the position of each sample with the
center of the leaf node containing it, the vector field
V could 4.3. Poisson Solution
be efficiently expressed as the linear sum of {Fo } by setting:
Having defined the vector field
V , we would like to solve for
q
F(q) = F̃ D . the function χ̃ ∈ FO,F such that the gradient of χ̃ is closest
2
to
V , i.e. a solution to the Poisson equation ∆χ̃ = ∇ ·
V .
This way, each sample would contribute a single term (the
normal vector) to the coefficient corresponding to its leaf’s One challenge of solving for χ̃ is that though χ̃ and the

c The Eurographics Association 2006.
 Kazhdan et al. / Poisson Surface Reconstruction

coordinate functions of
V are in the space FO,F it is not 4.4. Isosurface Extraction
necessarily the case that the functions ∆χ̃ and ∇ ·
V are. In order to obtain a reconstructed surface ∂ M̃, it is necessary
To address this issue, we need to solve for the function χ̃ to first select an isovalue and then extract the corresponding
such that the projection of ∆χ̃ onto the space FO,F is closest isosurface from the computed indicator function.
to the projection of ∇ ·
V . Since, in general, the functions We choose the isovalue so that the extracted surface
Fo do not form an orthonormal basis, solving this problem closely approximates the positions of the input samples. We
directly is expensive. However, we can simplify the problem do this by evaluating χ̃ at the sample positions and use the
by solving for the function χ̃ minimizing: average of the values for isosurface extraction:
2 2  1
∑ ∆χ̃ − ∇ ·
V , Fo = ∑ ∆χ̃ , Fo − ∇ ·
V , Fo . ∂ M̃ ≡ {q ∈ R3  χ̃ (q) = γ } with γ = ∑ χ̃ (s.p).
o∈O o∈O
|S| s∈S

Thus given the |O|-dimensional vector v whose o-th coordi- This choice of isovalue has the property that scaling χ̃ does
nate is vo = ∇ ·
V , Fo , the goal is to solve for the function not change the isosurface. Thus, knowing the vector field
V
χ̃ such that the vector obtained by projecting the Laplacian up to a multiplicative constant provides sufficient informa-
of χ̃ onto each of the Fo is as close to v as possible. tion for reconstructing the surface.

To express this in matrix form, let χ̃ = ∑o xo Fo , so that To extract the isosurface from the indicator function, we
we are solving for the vector x ∈ R|O| . Then, let us define the use a method similar to previous adaptations of the March-
|O| × |O| matrix L such that Lx returns the dot product of the ing Cubes [LC87] to octree representations (e.g. [WG92,
Laplacian with each of the Fo . Specifically, for all o, o ∈ O, SFYC96, WKE99]). However, due to the nonconforming
the (o, o )-th entry of L is set to: properties of our tree, we modify the reconstruction ap-
proach slightly, defining the positions of zero-crossings
  
∂ 2 Fo ∂ 2 Fo ∂ 2 Fo along an edge in terms of the zero-crossings computed by
Lo,o ≡ , Fo + , Fo + , Fo . the finest level nodes adjacent to the edge. In the case that an
∂ x2 ∂ y2 ∂ z2
edge of a leaf node has more than one zero-crossing associ-
Thus, solving for χ̃ amounts to finding ated to it, the node is subdivided. As in previous approaches,
we avoid cracks arising when coarser nodes share a face with
min
L x − v
2 .
x∈R|O| finer ones by projecting the isocurve segments from the faces
of finer nodes onto the face of the coarser one.
Note that the matrix L is sparse and symmetric. (Sparse
because
the Fo are compactly supported, and symmetric be- 4.5. Non-uniform Samples
cause f g = − f g .) Furthermore, there is an inherent
multiresolution structure on FO,F , so we use an approach We now extend our method to the case of non-uniformly dis-
similar to the multigrid approach in [GKS02], solving the tributed point samples. As in [Kaz05], our approach is to es-
restriction Ld of L to the space spanned by the depth d func- timate the local sampling density, and scale the contribution
tions (using a conjugate gradient solver) and projecting the of each point accordingly. However, rather than simply scal-
fixed-depth solution back onto FO,F to update the residual. ing the magnitude of a fixed-width kernel associated with
each point, we additionally adapt the kernel width. This re-
sults in a reconstruction that maintains sharp features in ar-
Addressing memory concerns In practice, as the depth in-
eas of dense sampling and provides a smooth fit in sparsely
creases, the matrix Ld becomes larger and it may not be prac-
sampled regions.
tical to store it in memory. Although the number of entries in
a column of Ld is bounded by a constant, the constant value
can be large. For example, even using a piecewise quadratic Estimating local sampling density Following the ap-
base function F, we end up with as many as 125 non-zero proach of [Kaz05], we implement the density computation
entries in a column, resulting in a memory requirement that using a kernel density estimator [Par62]. The approach is to
is 125 times larger than the size of the octree. estimate the number of points in a neighborhood of a sam-
ple by “splatting” the samples into a 3D grid, convolving the
To address this issue, we augment our solver with a block “splatting” function with a smoothing filter, and evaluating
Gauss-Seidel solver. That is, we decompose the d-th dimen- the convolution at each of the sample points.
sional space into overlapping regions and solve the restric-
We implement the convolution in a manner similar to
tion of Ld to these different regions, projecting the local so-
Equation 3. Given a depth D̂ ≤ D we set the density esti-
lutions back into the d-dimensional space and updating the
mator to be the sum of node functions at depth D̂:
residuals. By choosing the number of regions to be a func-
tion of the depth d, we ensure that the size of the matrix used WD̂ (q) ≡ ∑ ∑ αo,s Fo (q).
by the solver never exceeds a desired memory threshold. s∈S o∈Ngbr (s)
D̂

c The Eurographics Association 2006.
 Kazhdan et al. / Poisson Surface Reconstruction

Since octree nodes at lower resolution are associated with 5.1. Resolution
functions that approximate Gaussians of larger width, the We first consider the effects of the maximum octree depth
parameter D̂ provides away for specifying the locality of the on the reconstructed surface.
density estimation, with smaller values of D̂ giving sampling
density estimates over larger regions. Figure 3 shows our reconstruction results for the “dragon”
model at octree depths 6, 8, and 10. (In the context of recon-
Computing the vector field Using the density estimator, struction on a regular grid, this would correspond to reso-
we modify the summation in Equation 3 so that each sam- lutions of 643 , 2563 , and 10243 , respectively.) As the tree
ple’s contribution is proportional to its associated area on the depth is increased, higher-resolution functions are used to fit
surface. Specifically, using the fact that the area is inversely the indicator function, and consequently the reconstructions
proportional to sampling density, we set: capture finer detail. For example, the scales of the dragon,
which are too fine to be captured at the coarsest resolution

V (q) ≡ 1
∑ ∑ αo,s Fo (q). begin appearing and become more sharply pronounced as
s∈S WD̂ (s.p) o∈Ngbr (s) the octree depth is increased.
D

However, adapting only the magnitudes of the sample

contributions results in poor noise filtering in sparsely sam-
pled regions as demonstrated later in Figure 7. Therefore,
we additionally adapt the width of the smoothing filter F̃ to
the local sampling density. Adapting the filter width lets us
retain fine detail in regions of dense sampling, while smooth-
ing out noise in regions of sparse sampling.
Using the fact that node functions at smaller depths corre-
spond to wider smoothing filters, we define

V (q) ≡ 1
∑W (s.p) ∑ αo,s Fo (q).
s∈S D̂ o∈NgbrDepth(s.p) (s)

In this definition, Depth(s.p) represents the desired depth of

a sample point s ∈ S. It is defined by computing the average
sampling density W over all of the samples and setting:


Depth(s.p) ≡ min D, D + log4 (WD̂ (s.p)/W )
so that the width of the smoothing filter with which s con-
tributes to
V is proportional to the radius of its associated
surface patch Ps .

Selecting an isovalue Finally, we modify the surface ex-

traction step by selecting an isovalue which is the weighted
average of the values of χ̃ at the sample positions:
1
 ∑ W (s.p) χ̃ (s.p)
∂ M̃ ≡ {q ∈ R3  χ̃ (q) = γ } with γ = D̂
1
.
∑ W (s.p)
D̂
Figure 3: Reconstructions of the dragon model at octree depths 6
(top), 8 (middle), and 10 (bottom).

5. Results
To evaluate our method we conducted a series of experi- 5.2. Comparison to Previous Work
ments. Our goal was to address three separate questions:
We compare the results of our reconstruction algorithm
How well does the algorithm reconstruct surfaces? How
to the results obtained using Power Crust [ACK01], Ro-
does it compare to other reconstruction methods? And, what
bust Cocone [DG04], Fast Radial Basis Functions (Fas-
are its performance characteristics?
tRBF) [CBC∗ 01], Multi-Level Partition of Unity Implicits
Much practical motivation for surface reconstruction de- (MPU) [OBA∗ 03], Surface Reconstruction from Unorga-
rives from 3D scanning, so we have focused our experiments nized Points [HDD∗ 92], Volumetric Range Image Process-
on the reconstruction of 3D models from real-world data. ing (VRIP) [CL96], and the FFT-based method of [Kaz05].

c The Eurographics Association 2006.
 Kazhdan et al. / Poisson Surface Reconstruction

the sample points, result in reconstructions with spurious

surface sheets. Non-interpolatory methods, such as the ap-
proach of [HDD∗ 92] (e), can smooth out the noise, although
often at the cost of model detail. VRIP (f), the FFT-based
approach (g), and the Poisson approach (h) all accurately re-
construct the surface of the bunny, even in the presence of
noise, and we compare these three methods in more detail.

Figure 5: Reconstructions of a fragment of the Forma Urbis Ro-

mae tablet using VRIP (left) and the Poisson solution (right).

Comparison to VRIP A challenge in surface reconstruc-

tion is the recovery of sharp features. We compared our
method to VRIP by evaluating the reconstruction of sam-
ple points obtained from fragment 661a of the Forma Ur-
bis Romae (30 scans, 2, 470, 000 points) and the “Happy
Buddha” model (48 scans, 2, 468, 000 points), shown in Fig-
ures 5 and 6. In both cases, we find that VRIP exhibits a
“lipping” phenomenon at sharp creases. This is due to the
fact that VRIP’s distance function is grown perpendicular to
the view direction, not the surface normal. In contrast, our
Poisson reconstruction, which is independent of view direc-
tion, accurately reconstructs the corner of the fragment and
the sharp creases in the Buddha’s cloak.

Comparison to the FFT-based approach As Fig-

ure 4 demonstrates, our Poisson reconstruction (h) closely
Figure 4: Reconstructions of the Stanford bunny using Power
Crust (a), Robust Cocone (b), Fast RBF (c), MPU (d), Hoppe et al.’s
matches the one obtained with the FFT-based method (g).
reconstruction (e), VRIP (f), FFT-based reconstruction (g), and our Since our method provides an adaptive solution to the same
Poisson reconstruction (h). problem, the similarity is a confirmation that in adapting
the octree to the data, our method does not discard salient,
Our initial test case is the Stanford “bunny” raw dataset of high-frequency information. We have also confirmed that
362,000 points assembled from ten range images. The data our Poisson method maintains the high noise resilience al-
was processed to fit the input format of each algorithm. For ready demonstrated in the results of [Kaz05].
example, when running our method, we estimated a sample’s
Though theoretically equivalent in the context of uni-
normal from the positions of the neighbors; Running VRIP,
formly sampled data, our use of adaptive-width filters (Sec-
we used the registered scans as input, maintaining the regu-
tion 4.5) gives better reconstructions than the FFT-based
larity of the sampling, and providing the confidence values.
method on the non-uniform data commonly encountered in
Figure 4 compares the different reconstructions. Since the 3D scanning. For example, let us consider the region around
scanned data contains noise, interpolatory methods such as the left eye of the “David” model, shown in Figure 7(a). The
Power Crust (a) and Robust Cocone (b) generate surfaces area above the eyelid (highlighted in red) is sparsely sam-
that are themselves noisy. Methods such as Fast RBF (c) and pled due to the fact that it is in a concave region and is seen
MPU (d), which only constrain the implicit function near only by a few scans. Furthermore, the scans that do sample

c The Eurographics Association 2006.
 Kazhdan et al. / Poisson Surface Reconstruction

Figure 7: Reconstruction of samples from the region around the

left eye of the David model (a), using the fixed-resolution FFT ap-
proach (b), and Poisson reconstruction (c).

memory and time requirements of our algorithm are roughly

quadratic in the resolution. Thus, as we increase the oc-
tree depth by one, we find that the running time, the mem-
ory overhead, and the number of output triangles increases
roughly by a factor of four.
Figure 6: Reconstructions of the “Happy Buddha” model using
VRIP (left) and Poisson reconstruction (right).
Tree Depth Time Peak Memory # of Tris.
7 6 19 21,000
the region tend to sample at near-grazing angles resulting 8 26 75 90,244
in noisy position and normal estimates. Consequently, fixed- 9 126 155 374,868
resolution reconstruction schemes such as the FFT-based ap- 10 633 699 1,516,806
proach (b) introduce high-frequency noise in these regions.
In contrast, our method (c), which adapts both the scale and Table 1: The running time (in seconds), the peak memory usage (in
the variance of the samples’ contributions, fits a smoother re- megabytes), and the number of triangles in the reconstructed model
for the different depth reconstructions of the dragon model. A kernel
construction to these regions, without sacrificing fidelity in
depth of 6 was used for density estimation.
areas of dense sampling (e.g. the region highlighted in blue).

Limitation of our approach A limitation of our method The running time and memory performance of our method
is that it does not incorporate information associated with in reconstructing the Stanford Bunny at a depth of 9 is com-
the acquisition modality. Figure 6 shows an example of this pared to the performance of related methods in Table 2. Al-
in the reconstruction at the base of the Buddha. Since there though in this experiment, our method is neither fastest nor
are no samples between the two feet, our method (right) most memory efficient, its quadratic nature makes it scalable
connects the two regions. In contrast, the ability to use sec- to higher resolution reconstructions. As an example, Fig-
ondary information such as line of sight allows VRIP (left) ure 8 shows a reconstruction of the head of Michelangelo’s
to perform the space carving necessary to disconnect the two David at a depth of 11 from a set of 215,613,477 samples.
feet, resulting in a more accurate reconstruction. The reconstruction was computed in 1.9 hours and 5.2GB
of RAM, generating a 16,328,329 triangle model. Trying
to compute an equivalent reconstruction with methods such
5.3. Performance and Scalability
as the FFT approach would require constructing two voxel
Table 1 summarizes the temporal and spatial efficiency of grids at a resolution of 20483 and would require in excess of
our algorithm on the “dragon” model, and indicates that the 100GB of memory.

c The Eurographics Association 2006.
 Kazhdan et al. / Poisson Surface Reconstruction

Figure 8: Several images of the reconstruction of the head of Michelangelo’s David, obtained running our algorithm with a maximum tree
depth of 11. The ability to reconstruct the head at such a high resolution allows us to make out the fine features in the model such as the inset
iris, the drill marks in the hair, the chip on the eyelid, and the creases around the nose and mouth.

Method Time Peak Memory # of Tris. • Incorporate line-of-sight information from the scanning
Power Crust 380 2653 554,332 process into the solution process.
Robust Cocone 892 544 272,662 • Extend the system to allow out-of-core processing for
FastRBF 4919 796 1,798,154 huge datasets.
MPU 28 260 925,240
Hoppe et al 1992 70 330 950,562 Acknowledgements
VRIP 86 186 1,038,055
FFT 125 1684 910,320 The authors would like to express their thanks to the Stan-
Poisson 263 310 911,390 ford 3D Scanning Repository for their generosity in dis-
tributing their 3D models. The authors would also like to
Table 2: The running time (in seconds), the peak memory usage express particular gratitude to Szymon Rusinkiewicz and
(in megabytes), and the number of triangles in the reconstructed
Benedict Brown for sharing valuable experiences and ideas,
surface of the Stanford Bunny generated by the different methods.
and for providing non-rigid body aligned David data.
6. Conclusion
We have shown that surface reconstruction can be expressed References
as a Poisson problem, which seeks the indicator function that
[ABCO∗ 01] A LEXA M., B EHR J., C OHEN -O R D., F LEISHMAN
best agrees with a set of noisy, non-uniform observations,
S., L EVIN D., S ILVA C.: Point set surfaces. In Proc. of the
and we have demonstrated that this approach can robustly Conference on Visualization ’01 (2001), 21–28.
recover fine detail from noisy real-world scans.
[ABK98] A MENTA N., B ERN M., K AMVYSSELIS M.: A new
There are several avenues for future work: Voronoi-based surface reconstruction algorithm. Computer
• Extend the approach to exploit sample confidence values. Graphics (SIGGRAPH ’98) (1998), 415–21.

c The Eurographics Association 2006.
 Kazhdan et al. / Poisson Surface Reconstruction

[ACK01] A MENTA N., C HOI S., KOLLURI R.: The power crust, [NRDR05] N EHAB D., RUSINKIEWICZ S., DAVIS J., R A -
unions of balls, and the medial axis transform. Computational MAMOORTHI R.: Efficiently combining positions and normals
Geometry: Theory and Applications 19 (2001), 127–153. for precise 3D geometry. TOG (SIGGRAPH ’05) 24 (2005).
[BBX95] BAJAJ C., B ERNARDINI F., X U G.: Automatic recon- [OBA∗ 03] O HTAKE Y., B ELYAEV A., A LEXA M., T URK G.,
struction of surfaces and scalar fields from 3d scans. In SIG- S EIDEL H.: Multi-level partition of unity implicits. TOG (2003),
GRAPH (1995), 109–18. 463–470.
[BFGS03] B OLZ J., FARMER I., G RINSPUN E., S CHRÖDER P.: [Par62] PARZEN E.: On estimation of a probability density func-
Sparse matrix solvers on the GPU: Conjugate gradients and tion and mode. Ann. Math Stat. 33 (1962), 1065–1076.
multigrid. TOG 22 (2003), 917–924.
[PGB03] P ÉREZ P., G ANGNET M., B LAKE A.: Poisson image
[BMR∗ 99] B ERNARDINI F., M ITTLEMAN J., RUSHMEIER H., editing. TOG (SIGGRAPH ’03) 22 (2003), 313–318.
S ILVA C., TAUBIN G.: The ball-pivoting algorithm for surface
[SFYC96] S HEKHAR R., FAYYAD E., YAGEL R., C ORNHILL J.:
reconstruction. IEEE TVCG 5 (1999), 349–359.
Octree-based decimation of marching cubes surfaces. In IEEE
[Boi84] B OISSONNAT J.: Geometric structures for three dimen- Visualization (1996), 335–342.
sional shape representation. TOG (1984), 266–286.
[SOS04] S HEN C., O’B RIEN J., S HEWCHUK J.: Interpolating
[CBC∗ 01] C ARR J., B EATSON R., C HERRIE H., MITCHEL T., and approximating implicit surfaces from polygon soup. TOG
F RIGHT W., M C C ALLUM B., E VANS T.: Reconstruction and (SIGGRAPH ’04) 23 (2004), 896–904.
representation of 3D objects with radial basis functions. SIG-
[TO02] T URK G., O’B RIEN J.: Modelling with implicit surfaces
GRAPH (2001), 67–76.
that interpolate. In TOG (2002), 855–873.
[CL96] C URLESS B., L EVOY M.: A volumetric method for build-
ing complex models from range images. Computer Graphics [WG92] W ILHELMS J., G ELDER A. V.: Octrees for faster iso-
(SIGGRAPH ’96) (1996), 303–312. surface generation. TOG 11 (1992), 201–227.

[DG04] D EY T., G OSWAMI S.: Provable surface reconstruction [WKE99] W ESTERMANN R., KOBBELT L., E RTL T.: Real-time
from noisy samples. In Proc. of the Ann. Symp. Comp. Geom. exploration of regular volume data by adaptive reconstruction of
(2004), 428–438. iso-surfaces. The Visual Computer 15 (1999), 100–111.

[DMGL02] DAVIS J., M ARSCHNER S., G ARR M., L EVOY M.: [YZX∗ 04] Y U Y., Z HOU K., X U D., S HI X., BAO H., G UO B.,
Filling holes in complex surfaces using volumetric diffusion. In S HUM H.: Mesh editing with Poisson-based gradient field ma-
Int. Symp. 3DPVT (2002), 428–438. nipulation. TOG (SIGGRAPH ’04) 23 (2004), 641–648.

[EM94] E DELSBRUNNER H., M ÜCKE E.: Three-dimensional al-

pha shapes. TOG (1994), 43–72. Appendix A:
[FLW02] FATTAL R., L ISCHINKSI D., W ERMAN M.: Gradient The solution to surface reconstruction described in this paper
domain high dynamic range compression. In SIGGRAPH (2002),
approaches the problem in a manner similar to the solution
249–256.
of [Kaz05] in that the reconstructed surface is obtained by
[GKS02] G RINSPUN E., K RYSL P., S CHRÖDER P.: Charms: first computing the indicator function and then extracting the
a simple framework for adaptive simulation. In SIGGRAPH
appropriate isosurface.
(2002), 281–290.
[GWL∗ 03] G OODNIGHT N., W OOLLEY C., L EWIN G., L UE - While the two methods seem to approach the problem
BKE D., H UMPHREYS G.: A multigrid solver for boundary value of computing the indicator function in different manners
problems using programmable graphics hardware. In Graphics ( [Kaz05] uses Stokes’ Theorem to define the Fourier co-
Hardware (2003), 102–111. efficients of the indicator function while we use the Poisson
[HDD∗ 92] H OPPE H., D E ROSE T., D UCHAMP T., M C D ONALD equation), the two methods are in fact equivalent.
J., S TUETZLE W.: Surface reconstruction from unorganized
To show this, we use the fact that the Poisson equation
points. Computer Graphics 26 (1992), 71–78.
∆u = f where f is periodic can be solved using the Fourier
[Kaz05] K AZHDAN M.: Reconstruction of solid models from ori- transform. The Fourier series expansion is −|ζ |2 û(ζ ) =
ented point sets. SGP (2005), 73–82. fˆ(ζ ), or equivalently û(ζ ) = |−1 fˆ(ζ ).
ζ |2
[KSO04] KOLLURI R., S HEWCHUK J., O’B RIEN J.: Spectral
surface reconstruction from noisy point clouds. In SGP (2004), Thus, our Poisson equation ∆χ = ∇ ·
V can be solved us-
11–21.
ing χ̂ = −12 ∇
|ζ |
·
V . With the well known identity fˆ = −iζ fˆ
[LC87] L ORENSEN W., C LINE H.: Marching cubes: A high res-
olution 3d surface reconstruction algorithm. SIGGRAPH (1987), and its generalization ∇·
V = −iζ ·
Vˆ , we get χ̂ = i
|ζ |2
ζ ·
Vˆ ,
163–169. which is identical to [Kaz05].
[LGF04] L OSASSO F., G IBOU F., F EDKIW R.: Simulating water
and smoke with an octree data structure. TOG (SIGGRAPH ’04)
23 (2004), 457–462.
[Mur91] M URAKI S.: Volumetric shape description of range data
using “blobby model”. Computer Graphics 25 (1991), 227–235.

c The Eurographics Association 2006.
 PolyFit: Polygonal Surface Reconstruction from Point Clouds

Liangliang Nan Peter Wonka

Visual Computing Center, KAUST Visual Computing Center, KAUST
liangliang.nan@gmail.com pwonka@gmail.com

Abstract lightweight polygonal surface models of the objects.

We consider a method to be effective to this problem if it
We propose a novel framework for reconstructing meets the following requirements. First, the reconstructed
lightweight polygonal surfaces from point clouds1 . Unlike model should be an oriented 2-manifold and watertight,
traditional methods that focus on either extracting good ensuring physically valid models. This is essential for many
geometric primitives or obtaining proper arrangements of applications, e.g., simulation and fabrication. Second, it
primitives, the emphasis of this work lies in intersecting should be able to recover sharp features of the objects.
the primitives (planes only) and seeking for an appropriate Third, the method should not closely follow surface details
combination of them to obtain a manifold polygonal due to imperfections (i.e., noise, outliers, and missing data).
surface model without boundary. Instead, a lightweight representation is more preferred and
We show that reconstruction from point clouds can be beneficial for many applications (e.g., Augmented Real-
cast as a binary labeling problem. Our method is based on ity/Virtual Reality) that may involve many objects and
a hypothesizing and selection strategy. We first generate a large scenes. Besides, it is also helpful to provide a way
reasonably large set of face candidates by intersecting the to control the detail levels of the reconstructed models
extracted planar primitives. Then an optimal subset of the within the reconstruction pipeline. We proposed a novel
candidate faces is selected through optimization. Our opti- reconstruction framework that all the above requirements
mization is based on a binary linear programming formula- are satisfied in a single optimization process.
tion under hard constraints that enforce the final polygonal Instead of fitting dense smooth polygonal surfaces [8, 13,
surface model to be manifold and watertight. Experiments 22, 11], we seek for a more compact representation (i.e.,
on point clouds from various sources demonstrate that our simple polygonal surfaces) for piecewise planar objects.
method can generate lightweight polygonal surface models Our method is based on a hypothesizing and selection
of arbitrary piecewise planar objects. Besides, our method strategy. Specifically, we chose an optimal set of faces from
is capable of recovering sharp features and is robust to a large number of candidate faces to assemble a compact
noise, outliers, and missing data. polygonal surface model. The optimal set of faces is select-
ed through optimization that encourages good point support
and compactness of the final model. The manifold and
1. Introduction watertight requirements are enforced as hard constraints.
Reconstructing 3D models from sampled points has been Figure 1 shows an example of our reconstruction. Our work
a major problem in both computer vision and computer makes the following contributions:
graphics. Although it has been extensively researched in
the past few decades [8, 13, 22, 11, 3, 17, 27, 25, 12, 2, 15, • we cast polygonal surface reconstruction from point
20, 26, 16], obtaining faithful reconstructions of real-world clouds as a binary labeling problem based on a hypoth-
objects from unavoidable noisy and possibly incomplete esizing and selection strategy.
point clouds remains an open problem.
In this work, we focus on reconstructing piecewise • a surface reconstruction framework that is dedicated to
planar objects (i.e., man-made objects such as buildings). reconstructing piecewise planar objects.
Our inputs are point clouds sampled from real-world objects
that could be captured by various means (e.g., drones, hand- • a binary linear programming formulation for face se-
held scanners, and depth cameras). Our goal is to achieve lection that guarantees lightweight, manifold, and wa-
1 The source code of PolyFit is available under https://github. tertight reconstructions and meanwhile recovers sharp
com/LiangliangNan/PolyFit. features of the objects.

1
Figure 1. Pipeline. (a) Input point cloud. (b) Planar segments. (c) Supporting planes of the initial planar segments. (d) Supporting planes
of the refined planar segments. (e) Candidate faces. (f) Reconstructed model. All planar segments and faces are randomly colored.

2. Related Work from a sparse point cloud. Lin et al. [18] reconstruct
coarse building models by decomposing and fitting a set
In literature, a large body of research on polygonal of piecewise building blocks to the point clouds. By
surface reconstruction from point clouds aim at obtaining structuring and resampling planar primitives, Lafarge and
dense polygonal surfaces [8, 13, 22, 11]. In the last decades, Alliez [12] consolidate the point clouds and then obtain
extracting geometric primitives and identifying their com- surface models by solving a graph-cut problem. Using
bination to infer higher-level structures have become the a min-cut formulation, Verdie et. al. [26] reconstruct
most popular technique for reconstructing piecewise planar watertight buildings from an initial arrangement of planar
objects. In this section, we mainly review approaches that segments. These approaches specialize in reconstructing
focus on primitive extraction, primitive regularization, and urban buildings. Thus, they may not be suitable to handle
primitive- and hypothesis-based reconstruction. general piecewise planar objects. Based on space parti-
Primitive extraction. Researches falling in this cate- tioning and volumetric presentations, [3] and [2] obtain
gory aim at extracting high-quality instances of basic geo- piecewise planar reconstruction by determining if cells are
metric primitives (e.g., plane, cylinder) from point clouds occupied or not. These two techniques require visibility
corrupted by noise and outliers. The common practice information from the acquisition process as input while our
for this particular task is the Random Sample Consensus approach does not require this information.
(RANSAC) algorithm [6] and its variants [28, 24, 14]. Hypothesis-based reconstruction. By making stronger
These methods are reliable when the input data is contam- assumptions on the objects, researchers further regularize
inated by noise and outliers. So in our work, we use an the reconstruction problem and fit compound shapes (i.e.,
efficient implementation of the RANSAC algorithm [24] to a combination of multiple basic primitives) to the point
extract initial planar primitives. clouds [9, 16]. In Li et.al. [16], a set of axis-aligned boxes is
Primitive regularization. By exploiting the prior assembled to approximate the geometry of a building. Their
knowledge about the structure of an object, researchers strategy is generating a non-uniform grid and then selecting
further regularize the extracted primitives. Li et.al. [17] a subset of its cells that have good data support and are
discover global mutual relations between basic primitives smooth. In our work, we generalize this idea to reconstruct
and use such information as constraints to refine the general piecewise planar objects, and our reconstruction is
initial primitives base on local fitting and constrained based on optimization under hard constraints that guarantee
optimization. Monszpart et. al. [19] further exploit this idea manifold and watertight polygonal surface models.
to extract regular arrangements of planes. These methods
focus on inferring and regularizing potential relationship 3. Overview
between structures. However, obtaining manifold and
watertight surface models from the regularized primitives Our method takes as input a point cloud (possibly noisy,
remains unsolved. In our work, we provide a promising incomplete, and with outliers) of a piecewise planar object
framework to this challenging problem based on a and outputs a lightweight polygonal surface model. Our
hypothesizing and selection strategy. method consists of two main parts (see Figure 1):
Primitive-based reconstruction. In contrast with ap- Candidate face generation. We first extract a set of
proaches that focus on obtaining dense polygonal sur- planar segments from the point cloud using RANSAC [24].
faces, exploiting high-level primitives for man-made ob- Considering the detected planar segments may contain
jects became popular in the last years [25, 12, 15, 20, undesired elements due to noise, outliers, and missing data,
26, 16]. Arikan et al. [1] proposed an optimization-based we refine these planar segments by iteratively merging
interactive tool that can reconstruct architectural models plane pairs and fitting new planes. We then clip these
Figure 2. Two planes intersect with each other resulting in four parts (a). (b) to (g) show all possible combinations to ensure a 2-manifold
surface. The edge in (b) and (c) connects two co-planar parts, while in (d) to (g) it introduces sharp edges in the final model.

planes by an enlarged bounding box of the point cloud of the supporting planes for each pair of planar segments.
and compute pairwise intersections of the clipped planes to Then, starting from the pair (si , sj ) with the smallest angle,
hypothesize of the object’s faces. we test if the following two conditions are met. First, the
Face selection. We choose an optimal subset of the angle between the two planes is lower than a threshold,
candidate faces to assemble a manifold and watertight i.e., angle(si , sj ) < θt . Second, more than a specified
polygonal surface model. To do so, we formulate the face number (denoted as Nt ) of points lie on the supporting
selection as a binary linear programming problem. Our planes of both segments. If both conditions are satisfied,
objective function combines three terms that favor data we merge the two planar segments and fit a new supporting
fitting, point coverage, and model complexity, respectively. plane using PCA [10]. We iterate this process until no more
We also formulate hard constraints that ensure the final segment pair can be merged. In our experiments, we choose
model is manifold and watertight. θt = 10◦ and Nt = min(|si |, |sj |)/5, where |si | denotes
the number of supporting points of si . Figure 1 (c) and
4. Candidate Face Generation (d) show the extracted planes before and after refinement
respectively. It should be noted that an alternative approach,
The input to this step is a point cloud P of a piecewise
e.g., [17], can also be used to extract planar segments.
planar object, and the goal is to generate a reasonably large
Pairwise intersecting. To hypothesize the object’s
number of candidate faces.
faces, we crop the supporting planes of all planar segments
Plane extraction. We use the RANSAC-based primitive
by the bounding box of the point cloud. Then, the candidate
detection method proposed by Schnabel et al. [24] to detect
faces (i.e., a superset of the faces of the object) can be
a set of initial planar segments S = {si } from the point
obtained by intersecting the cropped planes. For simplicity,
cloud P . Here si is a set of points whose distances are
we compute pairwise intersections of the cropped planes
smaller than a threshold ε to a plane and each points can be
(see Figure 1 (e)). It should be noted that pairwise inter-
assigned to no more than one plane. This plane is denoted
sections introduce redundant candidate faces. Since most
as the supporting plane of si .
of the redundant faces do not represent actual structures of
Plane refinement. Due to the presence of noise and
an object, they are supported by no or very few (due to noise
outliers (especially for point clouds computed from Multi-
and outliers) point samples. The subsequent optimization-
view Stereo), RANSAC may produce a few undesired
based face selection is designed to favor choosing the most
planar segments. We observed that the undesired planar
confident faces while satisfying certain constraints.
segments usually have arbitrary orientations and poor point
The pairwise intersections maintain incidence informa-
support. Although the later optimization-based face se-
tion of the faces and edges. Each edge of a candidate face
lection procedure is designed to be capable of handling
is either connecting four neighboring candidate faces or
such outliers, there may still cause problems. First, these
representing a boundary. For example in Figure 2 (a), edge e
arbitrarily oriented planes may result in long but very
connects four faces while others are boundaries. We rely on
thin faces and small bumps in the final model. In some
such incidence information to formulate our manifold and
extreme cases, it may also introduce degenerate faces due
watertight constraints for face selection (see Section 5.2).
to the limit of floating point precision. This degeneracy
usually makes the manifold and watertight hard constraints
5. Face Selection
(see Section 5.2) impossible to be satisfied. Second, the
undesired candidate faces results in a larger optimization Given N candidate faces F = {fi |1 ≤ i ≤ N }
problem that is more expensive to be solved. generated in the previous step, we select a subset of these
To tackle this issue, we iteratively refine the initial planar candidate faces that can best describe the geometry of the
segments using an improved plane refinement algorithm object and ensure that the chosen faces form a manifold and
proposed in [15]. Specifically, we first compute the angle watertight polygonal surface. This is achieved through op-
timization. We define multiple energy terms that constitute
our objective function.
5.1. Energy terms
Let variable xi encode if a candidate face fi is chosen
(i.e., xi = 1) or not chosen (i.e., xi = 0), our objective
function consists of three energy terms: data-fitting, model
complexity, and point coverage.
Data-fitting. This term is intended to evaluate the fitting
quality of the faces to the point cloud while accounting for
an appropriate notion of confidence [21]. It is defined to Figure 3. Two examples of point coverage. The α-shape meshes
are in yellow and the candidate faces are in purple. The value
measure a confidence-weighted percentage of points that do
below each figure indicates the coverage ratio of each face.
not contribute to the final reconstruction
N
1 X Model complexity. Given incomplete point clouds
Ef = 1 − xi · support(fi ), (1)
|P | i=1 due to occlusions, the data-fitting term defined in Equa-
tion 1 tends to stubbornly comply with the incomplete
where |P | is the total number of points in P . support(fi ) data, resulting in gaps in the final model. Besides, noise
accounts for a notion of confidence and is defined at each and outliers also tend to introduce gaps and protrusions
point considering its local neighborhood in the reconstructed models. To avoid these defects, we
introduce a model complexity term to encourage simple
X dist(p, f )
support(f ) = (1 − ) · conf (p), structures (i.e., large planar regions). Considering gaps and
ε protrusions result in additional sharp edges, we define the
p,f |dist(p,f )<ε
(2) model complexity term as the ratio of sharp edges in the
where dist(p, f ) measures the Euclidean distance between model
a point p and a face f . We consider only points that are |E|
1 X
ε-close to f , i.e., points p satisfying dist(p, f ) < ε. The Em = corner(ei ), (4)
|E| i=1
confidence term conf (p) measures the local quality of the
point cloud at a point p. It is computed by examining the where |E| denotes the total number of pairwise intersections
local covariance matrices defined at p, as in [23]. In this in the candidate face set. corner(ei ) is an indicator function
work, we compute for p the covariance matrices of its local whose value is determined by the configuration of the two
neighbors at three static scales (i.e., different neighborhood selected faces connected by an edge ei . The intersecting
sizes). Then, conf (p) is defined as edges can be either planar or sharp. For example, the
3
intersecting edges e in Figure 2 (b) and (c) are planar edges
1X 3λ1i λ2 yielding larger planar polygons, but edges in (d) to (g) are
conf (p) = (1 − 1 2 3 ) · i3 , (3)
3 i=1 λi + λi + λi λi sharp edges that will introduce sharp features (if they are
selected in the final model). So corner(ei ) will have a value
where λ1i ≤ λ2i ≤ λ3i are the three eigenvalues of the of 1 if the faces associated with ei introduce a sharp edge
covariance matrix at scale i. conf (p) encode two geometric in the final model. Otherwise, corner(ei ) has a zero value
properties of the point samples near point p. The first meaning the two faces are coplanar.
property 1 − 3λ1 /(λ1 + λ2 + λ3 ) evaluates the quality of Point coverage. To handle missing data caused by oc-
fitting a local tangent plane at p, whose value of 0 indicates clusions, the unsupported regions (i.e., regions not covered
the worst point distribution and value of 1 indicates a by points) of the final model should be as small as possible.
perfectly fitted plane. The second property λ2 /λ3 evaluates To measure the coverage of a face fi , we first project all
the local sampling uniformity. Each eigenvalue ratio takes ε-close points onto fi . We then extract a mesh Miα by
on a value in the range [0, 1] with boundary values 0 and 1 constructing a 2D α-shape [5] from the projected points
corresponding to a perfect line distribution and a uniform (see Figure 3). We use only the points whose projections
disc distribution correspondingly. are inside fi for α-shape construction. We call Miα an α-
Intuitively, the data-fitting term favors choosing faces shape mesh. Intuitively, the α-shape mesh of a set of points
√
that are close to the input points and are supported by ensures that any three points with a circumradius rc ≤ α
densely sampled uniform regions. This term has a value are spanned by a triangle face. Given an appropriate value
in the range [0, 1] where a value of 1 indicates a noise-free of rc , the area of the α-shape mesh surface can provide us a
and outlier-free input data. reliable measure of the coverage of a candidate face by the
input points. Thus, our point coverage energy is defined as
the ratio of the uncovered regions in the model
N
1 X
Ec = xi · (area(fi ) − area(Miα )), (5)
area(M ) i=1

where area(M ), area(fi ), and area(Miα ) denote the

surface areas of the final model, a candidate face fi ,
and the α-shape mesh Miα of fi , respectively. In our
implementation, we chose the radius rc to be 5·density(P ),
where density(P ) denotes the average spacing of all the
points to their k nearest neighbors, kP
being set to 6.
N Figure 5. Reconstruction of a building with completely missing
Using the exact model area, i.e., i=1 xi · area(fi ), as
planes (top row). By adding an extra plane, model with different
the denominator ensures that the value of the point coverage structures can be obtained (bottom row).
term is within the range [0, 1]. However, this results in a
non-linear objective function that is difficult to optimize.
Considering the actual surface area of the final model is our method on various real-world objects of different com-
comparable to the area of the point cloud’s bounding box, plexity. Experiments demonstrated the advantages of our
we replace the exact model area with the area of the point hypothesizing and selection strategy.
cloud’s bounding box, i.e., area(M ) ≈ area(bbox(P )). Reconstruction results. Our method can reconstruct
piecewise planar objects such as buildings and other man-
5.2. Optimization
made objects. Figures 1 and 4 (a) - (d) show the reconstruc-
With the above energy terms, the optimal set of faces can tion of five buildings of different styles. The input of these
be obtained by minimizing a weighted sum of these terms buildings are point clouds computed from images using
under certain hard constraints enforcing the final model to Multi-View Stereo techniques. Although the point clouds
be manifold without boundary. are quite noisy and contain significant amounts of outliers
Remember that the candidate faces are obtained by the and missing data, our method faithfully reconstructed these
pairwise intersection of planes. Thus an edge is connected buildings. In Figures 4 (e) and (f), two indoor scenes
to either one (for boundary faces) or four faces (for inner consisting of multiple rooms are reconstructed. These two
faces). See Figure 2 (a) for an example. It is quite obvious scenes were captured by Google Tango tablets. Due to the
that the necessary and sufficient condition for manifold cluttered nature of indoor scenes, the datasets contains large
and watertight polygonal surfaces is that each edge of the holes. Our hypothesizing and selection strategy fills in the
model connects only two adjacent faces. Thus, the final missing regions and reconstructed the permanent structures
formulation for face selection can be written as (i.e., walls and roofs) of these scenes.
min λf · Ef + λm · Em + λc · Ec Besides the buildings, we also tested our method on
X other man-made objects. In Figure 4 (g), our method took
 X
 xj = 2 or 0, 1 ≤ i ≤ |E| (6) the laser scan of a small packing foam box as input and
s.t. j∈N (ei ) faithfully reconstructed all its structures with sharp features.
xi ∈ {0, 1}, 1≤i≤N We consider this as the first advantage of our approach. In

P (h) and (i), two sofas of different styles are reconstructed.
where j∈N (ei ) xj counts the number of faces connected
These pieces of furniture were scanned by PrimeSense
by an edge ei . This value is enforced to be either 0 or 2,
RGB-D cameras in the form of dense triangular meshes [4].
meaning none or two of the faces are selected (see Figure 2
Our method took the vertices of the meshes as input and
(b) - (g) for possible selections). These hard constraints
produced lightweight polygonal surface models.
guarantee that the final model is manifold and closed.
In rare cases where few planes are completely missing,
The problem defined in Equation 6 is a binary linear
our method may still generate plausible reconstructions. In
program. We solve it using the Gurobi solver [7]. After
Figure 5 (top), a building with its back roof completely
optimization, the union of the selected faces (i.e., faces
missing is reconstructed. It is interesting to see that by
having a corresponding variable xi = 1) comprises a
adding few extra planes, models with distinct structures can
polygonal surface model approximating the object.
also be obtained (bottom row). Thanks to the manifold
and watertight constraints, our method guarantees the re-
6. Results and Discussion
constructed models to be manifold without boundary. We
We implemented our method using C++. The α-shapes consider this as the second advantage of our approach.
are constructed using the CGAL library [5]. We tested In Table 1, we give statistics of our quantitative results.
Figure 4. Reconstruction of a set of piecewise planar objects from various data sources. The input for (a) - (d) are computed from images
using Multi-View Stereo; (e) and (f) are acquired by Google Tango tablets; (g) is captured by a laser scanner [17]; (h) and (i) are acquired
by PrimeSense Carmine RGB-D cameras [4]. From left to right: input point clouds, extracted planar segments, candidate faces (randomly
colored), reconstructed models, and models overlaid with the original point clouds (rendered with a smaller point size).
 Table 1. Statistics on the examples presented in Figure 4.
Model index in Figure 4 (a) (b) (c) (d) (e) (f) (g) (h) (i)
# points 129 K 101 K 73 K 577 K 186 K 176 K 382 K 756 K 911 K
# planar segments 36 21 22 22 28 33 52 17 19
# candidate faces 6239 472 946 959 1778 2770 7540 580 523
# faces of the final model 38 18 25 29 26 35 52 16 14
Primitive extraction (sec) 0.21 0.14 0.09 1.17 0.45 0.13 0.93 1.04 1.54
Candidate generation (sec) 0.05 0.01 0.01 0.02 0.01 0.02 0.06 0.02 0.03
Face selection (optimization) (sec) 2.73 1.46 1.17 8.26 2.90 3.05 7.39 13.26 16.01

Figure 6. Result without (a) and with (b) plane refinement step.
In (a), the top comprises ten faces originating from three different Figure 7. A building (a) reconstructed without (b) and with (c) the
planes. The red arrows indicate long but very thin polygonal faces. complexity term. The red arrows indicate bumps and gaps.

As can be seen from the number of faces in the final models, tion method lies in the optimization-based face selection
our method can produce lightweight 3D models. We that is designed to favor different aspects necessary for
consider this as the third main advantage of our approach. high-quality reconstructions. Both data fitting and coverage
Timings. Table 1 shows the running times of each are essential requirements. We observed that having only
step for the examples presented in Figure 4. We can see these two terms works perfectly for reasonably complete
that the face selection step is slower compared to primitive datasets, but it is still not sufficient for point clouds with
extraction and candidate face generation. It should be noted significant imperfections (i.e., noise, outliers, and missing
that the construction of the α-shape meshes dominated this data). Figure 7 (b) shows such an example. We can see that
step. Down-sampling of input point clouds may speed up the reconstructed model demonstrates some desired bumps
this process. and gaps. In contrast, with our model complexity term, the
Effect of plane refinement. We tested our algorithm reconstructed model is cleaner and compact (c). To further
on a few data sets by omitting the refinement step. One evaluate the behavior of the model complexity term, we
such example is shown in Figure 6. We can see that even ran our face selection on a dataset by gradually increasing
without the plane refinement step, our optimization still the weight (see Figure 9). Not surprisingly, increasing the
recovers the main structure of this building. However, influence of the model complexity term resulted in less
we observe that the two models have some differences in detailed 3D models. Thus, the model complexity term also
the details. First, the top-most face of the reconstructed provides control over the model details.
polyhedron in (a) is composed of ten candidate faces lying Parameters. The parameters for most examples are as
on three different planes, while in (b) it comprises fewer follows: λf = 0.46, λc = 0.27, and λm = 0.27 (Note
faces originating from the same plane. Second, some that weights in a wide range can produce the same results).
undesired bumpy structures can be avoided with the plane Slightly different weights (λf = 0.3, λc = 0.4, and λm =
refinement step. This is because some initial faces are 0.3) are used for the sofa example in Figure 4 (i), where the
quite close to each other, making the energy terms less background (ground plane) has a much higher density than
distinguishable. Another benefit of the plane refinement the object (sofa), thus the smaller data fitting weight.
step is the gain in efficiency. With plane refinement, Robustness and accuracy. We evaluated the robustness
the total number of candidate faces is reduced from 1185 of our approach by reconstructing from synthetic data
to 348. Accordingly, the run-time of the optimization with an increasing amount of noise (see Figure 10). In
decreased from 1.25 seconds to 0.31 seconds, i.e., more (c), the standard deviation of the noise distribution is
than four times faster. 0.6 meters high, our method still obtained topologically
Effect of the energy terms. The core of our reconstruc- accurate reconstruction. However, a much higher level
Figure 8. Comparison with four state-of-the-art methods on a building dataset. (a) Input point cloud. (b) Model reconstructed by the
Poisson surface reconstruction algorithm [11]. (c) The result of the 2.5D Dual Contouring approach [27]. (d) The result of [15]. (e) The
result of [16]. (f) Our result. The number under each sub-figure indicates the total number of faces in the corresponding model.

Limitations. Our hypothesizing and selection based

reconstruction strategy is intended for reconstructing simple
polygonal surfaces. It may encounter computation bot-
tlenecks for large complex objects (e.g., urban scenes of
many buildings). Running on such objects results in a huge
Figure 9. Reconstruction of the building shown in Figure 5 by number of candidate faces and the computation may not be
gradually increasing the influence of the model complexity term. affordable. In our experiments, we did not encounter such
The values are the weights used in the corresponding optimization. situations because we only tested our method on problems
of manageable scales.

7. Conclusions and Future Work

We introduced a novel framework that casts polygonal
surface reconstruction from point clouds as a binary label-
ing problem. Our framework is based on a hypothesizing
and selection strategy, and we proposed a novel optimiza-
tion formulation to obtain lightweight, watertight polygonal
surface models. We demonstrated the effectiveness of our
method on datasets captured by a range of devices. Since
our approach seeks to find an optimal combination of the in-
tersected planes under manifold and watertight constraints,
it is guaranteed to produce lightweight, manifold models
without boundary.
Figure 10. Reconstruction from synthetic data with increasing Future direction. Our current implementation exploit-
noise. Top: input. Middle: reconstruction. Bottom: reconstruction s planar primitives and it is suitable for reconstructing
error. σ indicates the standard deviation of the Gaussian noise. piecewise planar objects. However, the hypothesizing and
selection strategy is general and has the potential to handle
various types of basic primitives. In future work, we
of noise may further pollute the structures of the object,
would like to extend our method to incorporate more types
resulting in large errors in the final model.
of geometric primitives, such as cylinders and spheres.
Comparison. Figure 8 shows the comparison of our
To handle data with completely missing planes, we plan
approach with four other methods on a building dataset.
to infer the missing planes by exploiting structural priors
We can see that both Poisson [11] and 2.5D Dual Con-
(e.g., symmetry, parallelism, and orthogonality) to obtain
touring [27] generate dense surfaces with large numbers of
complete reconstructions.
bumps. The method in [15] uses only the roof information
and also produce a 2.5D reconstruction. By making the Acknowledgements
Manhattan-World assumption, the result of [16] is more
regularized. However, these two methods both produce We would like to thank Florent Lafarge, Michael Wim-
undesired structures. In contrast, our approach generates mer, and Pablo Speciale for providing us the data used in
the most compact and clean model. Besides, we can also Figures 4 (d), (a) , and (e)-(f), respectively. We also thank
observe that our result is less regular than that of [16]. Tina Smith for recording the voice-over for the video. This
However, this can be improved by a post-processing step research was supported by the KAUST Office of Sponsored
using the coarse model optimization technique proposed Research (award No. OCRF-2014-CGR3-62140401) and
in [20]. the Visual Computing Center (VCC) at KAUST.
References [19] A. Monszpart, N. Mellado, G. J. Brostow, and N. J.
Mitra. Rapter: Rebuilding man-made scenes with regular
[1] M. Arikan, M. Schwärzler, S. Flöry, M. Wimmer, and arrangements of planes. SIGGRAPH 2015, 34(4):103:1–
S. Maierhofer. O-snap: Optimization-based snapping for 103:12.
modeling architecture. ACM Transactions on Graphics,
[20] L. Nan, C. Jiang, B. Ghanem, and P. Wonka. Template
32:6:1–6:15, 2013.
assembly for detailed urban reconstruction. In Computer
[2] A. Boulch, M. de La Gorce, and R. Marlet. Piecewise-planar
Graphics Forum, volume 34, pages 217–228, 2015.
3d reconstruction with edge and corner regularization. In
[21] L. Nan, A. Sharf, H. Zhang, D. Cohen-Or, and B. Chen.
Computer Graphics Forum, volume 33, pages 55–64, 2014.
Smartboxes for interactive urban reconstruction. In
[3] A.-L. Chauve, P. Labatut, and J.-P. Pons. Robust piecewise-
SIGGRAPH 2010, volume 29, page 93.
planar 3d reconstruction and completion from large-scale
[22] Y. Ohtake, A. Belyaev, M. Alexa, G. Turk, and H.-P. Seidel.
unstructured point data. In CVPR 2010, pages 1261–1268.
Multi-level partition of unity implicits. In SIGGRAPH 2003,
[4] S. Choi, Q.-Y. Zhou, S. Miller, and V. Koltun. A large dataset
pages 463–470.
of object scans. arXiv preprint:1602.02481, 2016.
[23] M. Pauly, N. J. Mitra, J. Giesen, M. H. Gross, and L. J.
[5] T. K. F. Da. 2D alpha shapes. In CGAL User and Reference
Guibas. Example-based 3d scan completion. In Symposium
Manual. CGAL Editorial Board, 4.9 edition, 2016.
on Geometry Processing 2005, pages 23–32.
[6] M. A. Fischler and R. C. Bolles. Readings in computer
vision: Issues, problems, principles, and paradigms. chapter [24] R. Schnabel, R. Wahl, and R. Klein. Efficient ransac for
Random Sample Consensus: A Paradigm for Model point-cloud shape detection. In Computer graphics forum,
Fitting with Applications to Image Analysis and Automated volume 26, pages 214–226, 2007.
Cartography, pages 726–740. 1987. [25] C. A. Vanegas, D. G. Aliaga, and B. Benes. Automatic
[7] Gurobi. Gurobi optimization. http://www.gurobi. extraction of manhattan-world building masses from 3d
com/. laser range scans. IEEE transactions on visualization and
[8] H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, and computer graphics, 18(10):1627–1637, 2012.
W. Stuetzle. Surface reconstruction from unorganized points. [26] Y. Verdie, F. Lafarge, and P. Alliez. Lod generation for urban
In SIGGRAPH 1992, pages 71–78. scenes. ACM Transactions on Graphics, 34(3), 2015.
[9] H. Jiang and J. Xiao. A linear approach to matching cuboids [27] Q.-Y. Zhou and U. Neumann. 2.5d dual contouring: A robust
in rgbd images. In CVPR 2013. approach to creating building models from aerial lidar point
[10] I. Jolliffe. Principal component analysis. Encyclopedia of clouds. In ECCV 2010, pages 115–128.
Statistics in Behavioral Science. [28] M. Zuliani, C. S. Kenney, and B. Manjunath. The
[11] M. Kazhdan, M. Bolitho, and H. Hoppe. Poisson surface multiransac algorithm and its application to detect planar
reconstruction. Symposium on Geometry Processing 2006, homographies. In ICIP 2005, volume 3, pages III–153.
pages 61–70.
[12] F. Lafarge and P. Alliez. Surface reconstruction through
point set structuring. In Computer Graphics Forum,
volume 32, pages 225–234, 2013.
[13] M. Levoy, K. Pulli, B. Curless, S. Rusinkiewicz, D. Koller,
L. Pereira, M. Ginzton, S. Anderson, J. Davis, J. Ginsberg,
J. Shade, and D. Fulk. The digital michelangelo project: 3d
scanning of large statues. In SIGGRAPH 2000, pages 131–
144.
[14] B. Li, R. Schnabel, J. Shiyao, and R. Klein. Variational
surface approximation and model selection. Computer
Graphics Forum, 28(7), 2009.
[15] M. Li, L. Nan, N. Smith, and P. Wonka. Reconstructing
building mass models from uav images. Computers &
Graphics, 54(C):84–93, 2016.
[16] M. Li, P. Wonka, and L. Nan. Manhattan-world urban
reconstruction from point clouds. In ECCV 2016, pages 54–
69.
[17] Y. Li, X. Wu, Y. Chrysathou, A. Sharf, D. Cohen-Or,
and N. J. Mitra. Globfit: Consistently fitting primitives
by discovering global relations. In SIGGRAPH 2011,
volume 30, page 52.
[18] H. Lin, J. Gao, Y. Zhou, G. Lu, M. Ye, C. Zhang, L. Liu,
and R. Yang. Semantic decomposition and reconstruction
of residential scenes from lidar data. SIGGRAPH 2013,
32(4):66.