Eurographics Symposium on Point-Based Graphics (2004)

M. Alexa, S. Rusinkiewicz, (Editors)

Meshing Point Clouds Using Spherical Parameterization

M. Zwicker † and C. Gotsman2‡

Massachusetts Institute of Technology

2 Harvard University

Abstract
We present a simple method for meshing a 3D point cloud to a manifold genus-0 mesh. Our approach is based on
recent methods for spherical embedding of planar graphs, where we use instead a k-nearest neighborhood graph
of the point cloud. Our approach proceeds in two steps: We first embed the neighborhood graph on a sphere using
an iterative procedure, minimizing the tangential Laplacian. Then we triangulate the embedded points and apply
the resulting mesh connectivity to the input points. Besides meshing, spherical embedding of point clouds may also
be used for other applications such as texture mapping or morphing.
Categories and Subject Descriptors (according to ACM CCS): I.3.5 [Computer Graphics]: Computational Geometry
and Object Modeling

1. Introduction and Reimers [FR01] for meshes with disk topologies. The
advantage of this method is its simplicity and robustness.
Contemporary 3D scanners produce sets of 3D data points,
It guarantees that the resulting mesh will always be closed
sampled from the surface of a 3D object. These points
manifold genus-0.
are frequently unorganized, and to use them in 3D appli-
cations requires computing a polygon (usually triangular)
mesh which best approximates the sampled surface. This 2. Previous Work
means associating a connectivity structure with the point set.
Floater and Reimers [FR01] proposed to mesh a 3D point
Many algorithms exist for meshing 3D point clouds set sampled from a manifold surface with disk topology and
(e.g. [ACK01, FR01, HR02]). The algorithms differ in tech- single boundary loop as shown in Figure 1:
nique, assumptions made on the input, complexity, robust-
ness and reproducability. By "reproducability", we mean that Function M = Mesh(point set V )
the algorithm, when fed as input the vertex set of a given 3D 1. Construct a graph G = hV, Ei by connecting each
mesh, will usually not reproduce the original connectivity vertex to its k nearest neighbors (k is a user param-
and not even the topology of the original 3D mesh. In fact, eter) in Euclidean space.
only a few of the algorithms, most notably the Crust family 2. Determine which sequence of vertices will form the
of algorithms [ACK01], provide conditions on the sampling boundary B of M.
density which guarantee that the output will be close to the 3. Embed G in the plane such that B forms a convex
original. shape and V − B are positioned at convex combi-
nations of their neighbors. Call the planar point set
In this paper we make use of recent results on spherical V 0.
embedding to mesh a 3D point cloud which has been sam- 4. Form G0 = hV 0 , E 0 i by triangulating V 0 (e.g. by De-
pled from a closed manifold genus-0 surface, i.e., a topolog- launay triangulation).
ical sphere. This generalizes a method proposed by Floater 5. Use E 0 to form M = hV, E 0 i. The triangles are the
faces of the mesh.
Figure 1: The Floater-Reimers algorithm to mesh a disk-like
† matthias@[Link]
point set.
‡ gotsman@[Link]

c The Eurographics Association 2004.

 M. Zwicker & C. Gotsman / Meshing Point Clouds Using Spherical Parameterization

This procedure works reasonably well for surfaces with operator. This somewhat slow procedure usually converges
a boundary, the results depending on the specific recipe to a local minimum of the spring energy, which Gotsman et
chosen for the convex combination weights. However, al. speculate form a four to six dimensional subspace.
even a sophisticated choice of the weights (e.g. shape-
preserving [Flo97], harmonic [PP93], mean-value [Flo03])
3. Meshing using Spherical Embedding
will result in a significant metric distortion of the 2D em-
bedding relative to the original 3D mesh geometry if the We propose to use spherical embedding to mesh a point
mesh contains significantly curved regions. The results are cloud which is known to have spherical topology. This in-
also affected by the specific choice of the vertices forming volves eliminating step 2, replacing step 3 of the Floater-
the mesh boundary, which is usually heuristic, and the con- Reimers Algorithm (Figure 1) with a spherical embedding
vex shape which the boundary is mapped to in the embed- routine, and step 4 by a spherical triangulation routine. The
ding procedure. Obviously, if the 3D object has the topology result will, by definition, be a closed triangular genus-0 man-
of a sphere, this method is not the most natural choice. In a ifold mesh. The modified algorithm for meshing spherical
follow-up paper to [FR01] , Hormann and Reimers [HR02] point clouds is shown in Figure 2.
show how to mesh a spherical point cloud by segmenting it
to a number of disk-like subsets. This results in a somewhat Function M = Mesh(point set V )
complicated algorithm. 1. Construct a graph G = hV, Ei by connecting each
vertex to its k nearest neighbors (k is a user param-
Recently, researchers realized that for many mesh pro- eter) in Euclidean space.
cessing operations it is more natural to parameterize a closed 2. Embed G on the sphere. Call the spherical point set
manifold genus-0 mesh to a sphere, rather than cutting it in V 0.
various ways to reduce it to the case of a disk. The parame- 3. Form G0 = hV 0 , E 0 i by triangulating V 0 (e.g. by
terization operation is actually an embedding. However, em- spherical Delaunay triangulation).
bedding a graph on the sphere is much more difficult than 4. Use E 0 to form M = hV, E 0 i.
embedding in the plane, especially since the latter may be
done using essentially linear methods, and spherical em- Figure 2: Meshing a spherical point set.
bedding seems to be essentially non-linear. Parameterizing
a given closed manifold genus-0 3D mesh involves embed- The key step is the second one: spherical embedding. Al-
ding the 3D vertices on the sphere, such that the spherical though the graph G is not planar, it may be embedded on the
polygons induced by the mesh connectivity do not overlap. sphere using the same technique as proposed by [GGS03]
In the meshing application, where there is no given connec- and [GY02].
tivity to respect, the only requirement is that the distribution
of the points over the sphere is "similar" to their distribu- We now briefly describe our algorithm for spherical em-
tion in space. Essentially we would like the metric distortion bedding, but refer the reader to [GGS03] and [GY02] for
to be minimal, meaning that short-range distances should be more details and theoretical insights. Let us denote the po-
preserved as much as possible. Hence spherical embedding sitions of the input points in the set V by pi and the points
may be viewed as a "graph-drawing" operation in our con- embedded on the sphere by ui . Gu and Yauand [GY02] pro-
text, where the edge lengths are to be preserved as much as posed to embed a graph on a curved surface (in particular,
possible. the unit sphere) using the Laplace-Beltrami operator, which
is basically the tangential component of the Laplace oper-
Methods to embed a planar graph on the sphere were ator. Discrete approximations of the Laplace operator at a
proposed in [Ale00, DG97, GY02, KVLS, ST98]. These re- point ui have the form
quire a triangular graph as input. A simpler scheme was
Li = ∑ wi j ui − u j ,


proposed by Gu and Yau [GY02], and by Gotsman et (1)
j∈Ni
al. [GGS03]. This method is a generalization of the convex
combination method of Floater [Flo97] for planar graph em- where Ni is the index set of the neighbors of point i. In
bedding, but is non-linear. Essentially it means that instead our case, these neighbors are given by the k-nearest neigh-
of each vertex being located at some convex combination of bor graph of the input points pi . Typically, the weights wi j
its neighbor’s locations in the plane, the vector difference are chosen to be strictly positive. Popular options are unit
between this convex combination and the vertex location on weights, i.e., wi j = 1 for all j ∈ Ni , or inverse edge length
the sphere has only a radial component. This method has weights, i.e., wi j = 1/kpi − p j k for all j ∈ Ni . Other choices
the physical interpretation of a spring system (a zero-length are discussed in [Flo97, Flo03, FR01].
spring corresponding to each edge). The vertices are relaxed
On the unit sphere, the Laplace-Beltrami operator is sim-
to their minimal energy state, subject to the constraint that
ply approximated by the tangential component of the dis-
they are all located on the sphere. Computing the embedding
crete Laplacian:
involves solving a system of bilinear equations, which Gu
p
and Yau suggest to do iteratively using the Laplace-Beltrami Li = Li − (Li · Ni ) Ni , (2)

c The Eurographics Association 2004.

 M. Zwicker & C. Gotsman / Meshing Point Clouds Using Spherical Parameterization

where Ni = ui /kui k is the unit normal of point ui . Simi- in the spherical embedding algorithm (Figure 3). Inspired
lar as for planar embeddings [Flo97, Flo03, FR01], Gu and by [YBS04], we compute adaptive weights as
Yau [GY02] and Gotsman et al. [GGS03] showed that em-
wi j = kui − u j k/kpi − p j k for all j ∈ N − i. (4)
beddings on the sphere can be obtained by solving
p As shown by Yoshizawa et al. this minimizes the stretch or
Li = 0 for all i. (3)
distortion of the embedding.
This is a non-linear system of equations that we solve using To triangulate the points embedded on the sphere, we used
a simple iterative procedure shown in Figure 3, similar as the convex hull routine qconvex implemented in the qhull
proposed by [Ale00, GY02, KVLS]. package [qhu]. The resulting spherical triangulation is a De-
launay triangulation with relatively short edge lengths.
Function U = SphericalEmbedding(point set V )
1. Initialize ui = pi for all i
2. Translate the points ui such that ∑ ui = 0 4. Experimental Results
3. Project the ui onto the unit sphere We have run our meshing algorithm on a variety of 3D
4. For all i point datasets, three of which are shown in Figure 4: Tweety
p
- ui = ui + λLi (19,818 points), Max (52,809 points), and Igea (134,345
- ui = ui /kui k points).

5. Iterate from step 2. until convergence In Figure 5, we show triangulations of the Max and Igea
models, and the corresponding embeddings on the sphere.
Figure 3: Embedding a point set on the unit sphere.
We used adaptive weights and k = 25 neighbors.
Figure 8 shows results for different choices of the weights.
Note that λ is a damping coefficient that we set to λ = 0.5 Figure 8(a) was generated using uniform weights, (b) using
in all our experiments. inverse edge length weights, and (c) using adaptive weights.
The close-ups in Figure 8(d,e,f) show the differences in the
resulting meshes. The large distortion produced by uniform
3.1. Algorithm Parameters and inverse edge length weights leads to very thin triangles,
while adaptive weights produce a more natural triangula-
Our meshing algorithm has a number of parameters: the
tion. Figure 8(g,h,i) show the corresponding spherical em-
number of neighbors k used to form the embedded graph
beddings. Clearly, uniform weights and inverse edge length
connectivity, the convex combination weight recipe used to
weights lead to much more nonuniform distributions (i.e.,
weight these edges, and the spherical triangulation method.
high distortion) of the vertices on the sphere than adaptive
The average number of neighbors (valence) in a spherical weights. As proposed by Hormann and Reimers [HR02],
triangle mesh is six, so the value of k should be at least this, standard mesh optimization tools may be applied as a post-
else we risk the resulting graph not being well-connected. process to improve the quality of the final meshes.
This can lead to degenerate solutions of the spherical embed-
ding iteration, i.e., all points collapse to the same location.
4.1. Complexity
For uniformly sampled surfaces, values bigger than nine usu-
ally work fine. To compute all the results shown below, we Constructing the k-nearest neighbors graph on n points may
used a more conservative value of k = 25. be done in O(nlogn) [Cla83]. Computing the spherical em-
bedding iteratively seems to run in O(n2 ). Computing the
We used weights which are either uniform (i.e. inde-
3D convex hull requires O(nlogn) [Ski97], so the entire al-
pendent of the input geometry), proportional to the inverse
gorithm complexity is no more than O(n2 ).
edge lengths, or adaptive weights similar to those proposed
by Yoshizawa et al. [YBS04]. Surprisingly, there did not Our experimental results suggest that the spherical em-
seem to be much difference in results when uniform and bedding procedure does not converge to a vanishing value of
geometry-dependent weights were used. This is probably be- the Laplace-Beltrami operator, probably because the embed-
cause the k nearest neighbors are more or less at the same ded graph is not planar. In practice the iterative procedure is
distance from the vertex, hence distance-dependent weights stopped when it reaches a local minimum. The convergence
will be close to uniform. of the iteration with the Tweety model for different values of
k is illustrated in Figure 6. Here we plot the average length
We were able to achieve significantly different results
of the tangential Laplacian over the number of iterations. We
with error adaptive weights [YBS04]. In this procedure, the
suspect that the residual error is larger for larger values of k
weights do not only depend on the geometry (points pi ), but
because the corresponding graph is “less planar”.
also on the current embedding (points ui ). Hence, we need
to update the weights after each iteration over all the points The residual error is also different for the various choices

c The Eurographics Association 2004.

 M. Zwicker & C. Gotsman / Meshing Point Clouds Using Spherical Parameterization

(a) (b) (c)

Figure 4: Sample input point clouds: (a) Tweety (19,818 points), (b) Max (52,809 points), and (c) Igea (134,345 points).

(a) (b) (c) (d)

Figure 5: (a) Triangular mesh of the Max model, (b) Embedding of the Max model on the sphere, (c) Triangular mesh of the
Igea model, (d) Embedding of the Igea model on the sphere.

of the weights. As shown in Figure 7, uniform weights lead Model Number of Points Runtime (sec)
to a smaller error than inverse edge length and adaptive Tweety 19,818 125
weigths. However, adaptive weights converge faster, which Max 52,809 424
is further illustrated in Figure 9. This figure shows how the Igea 134,345 843
tail of the Tweety, which generates a fold-over in the ini- Table 1: Runtime of the spherical embedding iteration (1600
tial projection, is unfolded during the iteration. Adaptive iteration steps) for different models.
weights lead to a faster unfolding because they do not only
depend on the geometry, but also on the distortion produced
by the embedding [YBS04].

The current runtimes for the spherical embedding on a Preliminary results suggest that the spherical embedding
state-of-the-art PC (2.8 GHz P4 processor with 1 GB RAM) procedure can be accelerated significantly using hierarchi-
are reported in Table 1. For all three models, we performed cal algebraic multigrid methods, similarly to Aksoylu et
1600 iterations for these measurements. al. [AKS04]. Triangulation of the convex hull of the embed-

c The Eurographics Association 2004.

 M. Zwicker & C. Gotsman / Meshing Point Clouds Using Spherical Parameterization
2
10 input. This could be done in a pre-process where points too
far from the average of their neighbors are removed.
Average length of tangential Laplacian

We tested our algorithm on a noisy point cloud that was

obtained by randomly perturbing the points of a smooth
3 model along their normal direction. We chose a range of
10

perturbation b relative to the average distance to the k = 25

(a)
nearest neighbors. Figure 10 shows a mesh of the Max data
set with b = 1. We also applied smoothing after meshing to
4
improve the quality of the result, as shown in Figure 10(c).
10
(b) Figure 10(d,e,f) shows close-ups of results for b = 0, b = 1,
(c)
and b = 3.

5. Conclusion and Discussion

5
10 0 1 2 3 4
10 10 10 10 10 We have described an algorithm for spherical meshing of a
Iterations
3D point cloud which is simple and efficient. As its main
Figure 6: Convergence of the embedding iteration with uni- advantage, it guarantees a closed manifold genus-0 result,
form weights for different parameters k: (a) k=100, (b) even for very noisy inputs. As with every other meshing al-
k=25, (c) k=9. gorithm, a number of independent pre-processing techniques
may be applied to the input point set and post-processing
techniques to the output mesh.
2
10
Embedding the point set on the sphere has an important
side-effect: the coordinates of the embedded points on the
Average length of tangential Laplacian

sphere may be used as spherical texture coordinates, or con-

verted to planar texture coordinates in a variety of standard
10
3 ways, to allow texture mapping onto the model. The embed-
ding, when viewed as a parameterization, can also be used to
establish a one-to-one correspondence between different ob-
(a) jects, which is useful for applications such as morphing. Fur-
ther, the spherical parameterization could be used to recon-
4 (b)
10
(c) struct a continuous surface similar as described by Zwicker
et al. [ZPKG] for points parameterized to a planar domain.
In our experiments we observed that we do not obtain
spherical embeddings with completely vanishing tangential
5
10 0
10 10
1
10
2
10
3
10
4 Laplacians. While it has been proven that it is possible to
Iterations embed planar 3-connected graphs on the sphere with vanish-
ing Laplacians [GGS03], we do not know of any result for
Figure 7: Convergence of the embedding iteration for k = 25 non-planar k-nearest neighbor graphs. An in-depth analysis
and different choices of the weights: (a) adaptive weights, (b) of this problem seems interesting. In the future, we would
inverse edge length weights, (c) uniform weights. also like to investigate hierarchical techniques such as alge-
braic multigrid or geometric clustering (as in [ZPKG] for
planar embedding) to accelerate the spherical embedding it-
eration. We also plan to experiment with a graph given by lo-
ded points using qconvex [qhu] took less than five seconds cal Delaunay neighborhoods as proposed by Floater [FR01]
for all models. instead of using k-nearest neighbors. This might improve the
convergence of the embedding procedure.
4.2. Sensitivity to Noise
Acknowledgements
The algorithm in its basic form will always produce a closed
manifold genus-0 mesh on the entire input data set, which is The work of Craig Gotsman was partially supported by Is-
a major advantage of this particular algorithm. In the event rael Ministry of Science grant 01-01-01509, German-Israel
of noisy input, however, this mesh may intersect itself. Note Fund (GIF) Grant I-627-45.6/1999, and European FP6 NoE
that this does not contradict the fact that the mesh is man- grant 506766 (AIM@SHAPE). Matthias Zwicker is partially
ifold. It can be rectified by a post-processing smoothing supported by a stipend of the Swiss National Science Foun-
stage. It would also be desirable to eliminate outliers in the dation.

c The Eurographics Association 2004.

 M. Zwicker & C. Gotsman / Meshing Point Clouds Using Spherical Parameterization

(a) (b) (c)

(d) (e) (f)

(g) (h) (i)

Figure 8: Sample results using k = 25 neighbors in the graph: (a,d) use uniform weights, (b,e) inverse edge length weights, and
(c,f) were computed using the adaptive weights; (g,h,i) depict the corresponding embeddings on the sphere.

References P.: Multilevel solvers for unstructured surface

meshes. To appear in SIAM Journal of Scientific
[ACK01] A MENTA N., C HOI S., KOLLURI R.: The power Computing (2004).
crust. In Proceedings of 6th ACM Sympo-sium on
Solid Modeling (2001), pp. 249–260.
[Ale00] A LEXA M.: Merging polyhedral shapes with
[AKS04] A KSOYLU B., K HODAKOVSKY A., S CHRÖDER scattered features. The Visual Computer 16, 1

c The Eurographics Association 2004.

 M. Zwicker & C. Gotsman / Meshing Point Clouds Using Spherical Parameterization

(a) (b) (c) (d)

(e) (f) (g) (h)

Figure 9: Unfolding of the tail of the tweety model. (a-d) uniform weights, (e-h) adaptive weights, (a,e) initial embedding, (b,f)
after 100 iterations, (c,g) after 200 iterations, (d,h) after 800 iterations.

(2000), 26–37. ceedings of Curve and Surface Design (2002),

pp. 215–224.
[Cla83] C LARKSON K.: Fast algorithms for the all nearest
neighbors problem. In Proceedings of 24th IEEE [KVLS] KOBBELT L., VORSATZ J., L ABISK U., S EIDEL
FOCS (1983), pp. 226–232. H.-P.: A shrink-wrapping approach to remeshing
polygonal surfaces. In Proceedings of Eurograph-
[DG97] DAS G., G OODRICH M.: On the complexity of ics 1999, pp. 119–130.
optimization problems for 3-dimensional convex
polyhedra and decision trees. Computational Ge- [PP93] P INKALL U., P OLTHIER K.: Computing discrete
ometry 8 (1997), 123–137. minimal surfaces and their conjugates. Experi-
mental Mathematics 2 (1993), 15–36.
[Flo97] F LOATER M.: Parameterization and smooth ap-
[qhu] : [Link].
proximation of surface triangulations. Computer
Aided Geometric Design 14 (1997), 231–250. [Ski97] S KIENA S.: The algorithm design manual. Telos,
1997.
[Flo03] F LOATER M. S.: Mean-value coordinates. Com-
puter Aided Geometric Design 20 (2003), 19–27. [ST98] S HAPIRO A., TAL A.: Polygon realization for
shape transformation. The Visual Computer 14,
[FR01] F LOATER M. S., R EIMERS M.: Meshless param- 8-9 (1998), 429–444.
eterization and surface reconstruction. Computer
Aided Geometric Design 18 (2001), 77–92. [YBS04] YOSHIZAWA S., B ELYAEV A., S EIDEL H.-P.: A
fast and simple stretch-minimizing mesh param-
[GGS03] G OTSMAN C., G U X., S HEFFER A.: Fundamen- eterization. In Proceedings of SMI’04, to appear
tals of spherical parameterization for 3d meshes. (2004).
In ACM Transactions on Graphics (Proceedings
[ZPKG] Z WICKER M., PAULY M., K NOLL T., G ROSS
of SIGGRAPH 2003) (2003), pp. 358–363.
M.: Pointshop3d: An interactive system for point-
[GY02] G U X., YAU S.-T.: Computing conformal struc- based surface editing. In Proceedings of SIG-
tures of surfaces. Communications in Information GRAPH 2002, pp. 322–329.
and Systems 2, 2 (2002), 121–146.

[HR02] H ORMANN K., R EIMERS M.: Triangulating

point clouds with spherical topology. In Pro-

c The Eurographics Association 2004.

 M. Zwicker & C. Gotsman / Meshing Point Clouds Using Spherical Parameterization

(a) (b) (c)

(d) (e) (f)

Figure 10: Triangular mesh of a noisy point cloud with b = 1: (a) Mesh, (b) Flat shaded version of (a), and (c) flat shaded mesh
after smoothing. Close-ups of some results on inputs with different amounts of noise: (d) original with b = 0, (e) b = 1, and (f)
b = 3. We used k = 25 and adaptive weights. While the entire mesh is guaranteed to be manifold genus-0, self-intersections may
occur.

c The Eurographics Association 2004.