arXiv:1711.10824v1 [cs.CV] 29 Nov 2017

Compression for Smooth Shape Analysis
V. Estellers
TUM
virginia.estellers@tum.de

F.R. Schmidt
TUM, BCAI
frank.schmidt12@de.bosch.com

Abstract

eigenfunctions of the Laplace-Beltrami operator necessary to estimate geodesics [7], compute shape descriptors [1, 28, 18, 19] or shape matches [25] with tractable
problem sizes. This creates a paradox, as shapes are first
discretized with fine meshes or large pointclouds to estimate high-dimensional shape operators that are mostly
discarded to analyze the shape. This paradox appears because we only need large meshes or dense pointclouds to
discretize accurate differential operators over the shape,
not for shape analysis.

Most 3D shape analysis methods use triangular meshes to
discretize both the shape and functions on it as piecewise
linear functions. With this representation, shape analysis requires fine meshes to represent smooth shapes and
geometric operators like normals, curvatures, or LaplaceBeltrami eigenfunctions at large computational and memory costs.
We avoid this bottleneck with a compression technique
that represents a smooth shape as subdivision surfaces
and exploits the subdivision scheme to parametrize smooth
functions on that shape with a few control parameters. This
compression does not affect the accuracy of the LaplaceBeltrami operator and its eigenfunctions and allow us to
compute shape descriptors and shape matchings at an accuracy comparable to triangular meshes but a fraction of
the computational cost.
Our framework can also compress surfaces represented
by point clouds to do shape analysis of 3D scanning data.

1

D. Cremers
TUM
cremers@tum.de

We propose a paradigm shift that avoids this paradox by
representing a shape with a subdivision surface [24] that
parametrizes the surface with a small set of smooth base
functions. Differentiability is then intrinsic to the shape
representation and does not require large bases to compute
principal curvatures, wave-kernel signatures, or LaplaceBeltrami eigenfunctions. Subdivision surfaces are a generalization of splines to surfaces of arbitrary topology that
provide set of base functions to compactly parametrize the
surface as well as smooth functions over it. As a result,
subdivision surfaces provide a compressed shape representation well suited to shape analysis.
Our first contribution (Section 3) is a method to estimate
a subdivision surface from the triangular meshes common in shape analysis[2]. This kind of manifold meshes,
however, are far from the meshes or pointclouds obtained
from range measurements or computer-vision models. As
these low-level representations proliferate, they create a
challenge for shape analysis to adapt to the characteristics of raw data (irregular surface sampling, noise, and
outliers) that cause havoc in the estimation of differential
operators. Up until now, raw measurements have been
largely ignored in shape analysis by estimating a mesh
representation from a pointcloud and performing analysis
on this representation. This avoids dealing with noise and
outliers but ignores the direct measurements and makes

Introduction

Many shape analysis tasks describe shapes as smooth manifolds and analyze them with respect to their geometry in
terms of normals, curvatures, and geodesics that require
accurate estimates of fist- and second-order derivatives
over the shape to compute tangent spaces, the Riemannian
metrics, or Laplace-Beltrami operators [28, 18, 34, 19].
When the surface is represented as a mesh or point cloud,
these differential operators can only be computed approximately at large memory and computational costs.
Once these differential operators have been used to compute normals or metrics, most of the high resolution information is discarded by operating only with the leading
1

input pointcloud
501007 points

triangular mesh from [15]
86246 vertices

30 geodesic level lines
cos(30πg)

5 geodesic 5 level lines
cos(5πg)

input pointcloud
501007 points

edge-collapsed mesh [15]
1500 vertices

30 geodesic level lines
cos(30πg)

5 geodesic 5 level lines
cos(5πg)

input pointcloud
501007 points

subdivision surface
1500 control vertices

30 geodesic level lines
cos(30πg)

5 geodesic 5 level lines
cos(5πg)

Figure 1: Computation of approximate geodesic g with different surface representations by the Heat method [7] .
Column 1: input Kinect pointcloud. Column 2: representation of the surface with a high-resolution mesh obtained
with Poisson reconstruction [15] (row 1), a low-resolution triangular mesh obtained by quadratic edge collapse [9] of
the high-resolution mesh (row 2), and with our subdivision surface (row 3). Columns 3 and 4: level lines of the geodesic
visualized as cos($g). Our subdivision surface is comparable to the surface obtained by Poisson reconstruction and
faithfully represent geodesics at low and high resolution with 2% of vertices. Compressing the Poisson mesh by edge
collapse looses all the small scale details of the surface and its geodesic, as highlighted by high-frequency level lines
of the geodesic.
shape analysis dependent on the model used to estimate
the surface.1.For this reasons, our second contribution is
a robust method that fits a subdivision surface to a noisy
pointcloud and estimates the geometric operators required
to do shape analysis on it. Our model is robust to both
noise and outliers in the pointcloud and to erroneous correspondences between the input points and the subdivision

surface. This is a critical aspect of fitting a parametric
surface to raw data without developing the non-manifold
artifacts that prevent shape analysis.

Our third contribution is a fast optimization algorithm
(Section A) that combines a sequential quadratic approximation with a good initialization strategy. This is important because fitting a subdivision surface is not a convex problem and vanilla initializations [13] oversmooth
1Least-squares models [15, 4], for instance, oversmooth surface edges
and require less eigenfunctions to match than the `1 models [31, 8] that the surface to avoid small-scale local minima and prevent
emphasize corners.
self-intersections.
2

Our fourth contribution (Section 2) describes how to
accurately compute many shape operators with a subdivision surface. We focus on Laplace-Beltrami operator,
wave kernel signatures [1], approximate geodesics [7] and
prove its efficacy in the full-pipeline task of shape matching [25]. Our experiments show state-of-the-art performance at a fraction of the memory requirements of shape
representations with triangular meshes.

2

The reference quad  = [0, 1] × [0, 1] parametrizes each
quadrilateral patch while the equivalence relation ∼ connects neighboring faces along their common edge.
A parametrization of the surface is then obtained by
defining a piecewise bilinear function
i
BM
: TM → R
i
for each vertex in xi ∈ TM satisfying BM
(xj ) = δij . The
resulting parametrization

Smooth Shapes

ΦM : TM → SM ⊂ R3

In this paper, a 3D shape is a smooth2, compact oriented
surface S ⊂ R3 without boundaries.
For computational purposes, shapes are usually discretized by a mesh M = (V, F ) in terms of two matrices:
the i-th row of V ∈ Rn×3 stores the location of the i-th
vertex of the mesh, while each row of F ∈ Nm×k stores
the indices of the vertices of a face of M. We focus on
quad meshes with k = 4.
In shape analysis, quad meshes usually discretize a
smooth surfaces as a piecewise C 1 surface that is not
smooth. Subdivision surfaces, on the other hand, only
use meshes to discretize the domain of a smooth surface
parametrization. To formally introduce this parametrization, we first revisit the topological space associated to a
quad mesh and show how it can be used to define a smooth
subdivision parametrization.

2.1

defines the piecewise C 1 surface SM := ΦM (TM ) by
bilinear interpolation
Φ(u) =

n
X

i
BM
(u) · vi .

i=1

With such a piecewise C 1 paremetrization, Riemannian
quantities like the first or second fundamental form can
only be approximated and need large meshes with small
quadrilaterals for accurate approximations. This limits
their use in real-time applications and calls for other shape
representations.

2.2

Smooth Subdivision Surfaces

A subdivision surface represents a smooth surface S by
a coarse control mesh M0 = (V 0 , F 0 ) and a subdivision scheme. The subdivision scheme transforms a mesh
Mk = (V k , F k ) into a finer mesh Mk+1 such that if we
start with M0 := (V, F ), we obtain the smooth surface as
the limit of iterating subdivision ad infinitum:

Meshes – Piecewise Bilinear Surfaces

The regular quad mesh parametrization of a surface SM
consists of quadrilaterals “glued” together along their
edges. Each quadrilateral in SM is modeled by a bilinear
patch that fits the fours vertices of a facet in M and is glued
to its neighboring facets. This topological gluing extends
the independent parametrization of each quadrilateral to
the whole mesh and is formalized by the topological space,
TM , that defines the domain of the parametrization [14].
Formally, given a quad mesh M = (V, F ), the topological space TM is defined by a set of quads that is indexed
by the mesh facets IF = {1, . . . , m} and glued together
by the equivalence relation ∼, that is TM :=  × IF / ∼.

S = lim SMk .
k→∞

In this paper we follow the Catmull-Clark subdivision
scheme Mk 7→ Mk+1 consisting of two steps: First, the
topology of the new mesh Mk+1 is defined by subdividing
k
each quad in TM
into its four quadrants that define four
k+1
adjacent facets in TM
. Second, the location of each
vertex of the new mesh Mk+1 is computed as a linear
combination of the positions of the nearby vertices. Figure
2 illustrates this subdivision process.
1
2
2 We call a C object smooth if it is C except for a finite subset.
This produces a compressed surface representation beThis differs from the mathematical notion, which identifies smoothness
with C ∞ .
cause the topology of Mk+1 is completely determined by
3

the topology of Mk . Similarly, the space of bilinear functions over TMk+1 created by subdivision is determined by
the space of bilinear functions over TMk . This means
that we can parametrize the i-th function created by subi
k
division BM
k+1 over the domain TM , and by induction
k+1
over TM , instead than TM . We denote these functions
k,i
BM
: TM → R and use them to parametrize the surface
SMk = ΦkM (TM ) as a function of the vertices of Mk .
In the second step, the position of the vertices of Mk+1
are computed as linear combination of the positions of the
vertices of Mk . The weights of this linear combination
are organized into a subdivision matrix Âk ∈ Rnk × nk−1
that computes the position of the vertices of Mk as

This scheme is designed for control grids with quadrilateral connectivity, but our compression generalizes to any
subdivision schemes that parametrizes C 2 surfaces.

3

Robust Fitting Model

We formulate the problem of fitting a subdivision surface
S to a set of surface samples P = {p1 , . . . , pN } and a
set of normals T = {t1 , . . . , tN } as the minimization
problem:
min E(S) := dist(S, P) + α · T (S, T) + β · R(S) ,
S
{z
} | {z }
| {z } |

Vk = Âk . . . Â1 V = Ak V.
| {z }

Point Fit

Ak ∈Rnk ×n

Tangent Fit

Regularization

(1)
Using this linear dependency, we can parametrize a surface
approximation SMk resulting from k subdivisions directly where α, β ≥ 0 are model parameters. The different terms
over TM as follows:
take different information into account: dist(S, P) fits the
n
n
n
k
k
surface S to the points P, T (S, T) matches the tangent
X k,i
X X k,i
ΦkM (u) =
BM (u) (Ak V )i =
BM (u)Ak,ij vj plane of S to the normal information T, and R(S) regui=1
j=1 i=1
larizes the discrete representation of S to avoid degenerate
By iterating this process, the subdivision surface of the faces in M. If no normal information is provided, α = 0.
limit

S = lim SMk

k→∞
Point Fit
can be formulated as a limit of functions defined over the
topological space TM . Evaluating the limit surface at the Let DS (p) = mins∈S ks − pk denote the Euclidean dislocation of the vertex thus amounts to computing the base tance between a point p ∈ R3 and the compact surface S,
functions
we define the point fit energy as
nk
X
k,i
Φj := lim
BM
· Ak,ij .
N
N
X
X
k→∞
q
q
i=1
dist(S, P) :=
DS (pj ) =
min kΦ(uj ) − pj k .
uj ∈TM
j=1
j=1
The existence of this limit and the smoothness of the
limit functions is guaranteed by the subdivision scheme.
All the limit base functions Φj are compactly supported Using the q-th power of the distance with q ∈ (1, 2), inbecause the subdivision rules act locally and converge to stead of the common squared distance, makes our model
spline basis functions. This equivalence provides analytic robust to outliers and scanning artifacts in the input samderivatives and fast evaluation techniques [29] for smooth ples. The minimization in u1 , . . . , uN results from the
functions or tangent vectors defined over the surface and parametric representation of the surface S = Φ(TM )
makes subdivision surfaces well-suited for shape analysis. and introduces a large number of additional variables
We focuses on Catmull-Clark subdivision [6] for sim- over which the objective function is not convex. We
plicity. This scheme ensures the existence and C 2 smooth- avoid the optimization of the correspondence parameters
ness of the basis functions everywhere except at extraor- u1 , . . . , uN by using a second-order approximation of the
dinary vertices (vertices with valence different than four). squared distance function to a surface in Section A.

4

M0
1701 Control Points

M1
6798 Vertices

M2
27186 Vertices

M3
108738 Vertices

Figure 2: Different subdivision levels Mk of Catmull-Clark subdivision scheme applied to a control mesh and a
smooth function defined over the mesh (its 18-th ∆ eigenfunction). M3 is already a good approximation of the smooth
surface S = M∞ .
Tangent Fit

3.1

To obtain a surface that also matches the normals, we
include the tangent energy

A major challenge of fitting any explicit parametrization
to a set of points is the lack of convex formulations for
the optimization problem. Non-convexity introduces the
difficulty of finding a suitable initialization that leads to
a good local minimum. It is an important step of the
algorithm that conditions the topology of the surface.
By parameterizing S with a subdivision surface, our
optimization problem is formulated in terms of the control
mesh M = (V, F ). To avoid solving a combinatorial
problem, we fix the topology of the control mesh, and
thus the topology of the surface, and only optimize for the
location of the vertices V . We design the control mesh
to satisfy two properties: it determines a surface with the
same topology as the input data, and it has a small number
of vertices to represent the surface compactly. Depending
on the input data, we apply different methods to obtain a
good initialization.
If the point samples P are the vertices of a mesh M,
the topology of the surface is already encoded by the
mesh. Without loss of generality, we assume that M
is a triangular mesh without boundaries and use quadratic
edge-collapse [9] to reduce the number of vertices in the
mesh and preserve the surface topology. We then transform this triangular mesh into a quad mesh by solving a
perfect matching problem that pairs triangles of the collapsed mesh to create quads [27]. For each edge e, let
α1 (e), . . . , α4 (e) ∈ S1 be the angles of the quad that results from removing e and η(e) ∈ S1 the angle between
the normals of the triangles incident to e. We set the cost

T (S, T) =

N
X

q

q

|htj , ∂1 Φ(uj )i| + |htj , ∂2 Φ(uj )i| ,

j=1

(2)
where tj is the input surface normal at point pj and
∂i Φ(uj ) the i-th base vector of the tangent space TΦ(uj ) S.
Each term in T aligns the tangent space of the subdivision
surface in the direction orthogonal to the sample normal.
This penalty is designed to be independent of the orientation of the normals to account for noisy point clouds
where estimating consistent normal orientations is prone
to fail. The `q -norm makes this term robust to noise and
outliers in the normals and the correspondence parameters
u1 , . . . , uN .
Regularization
The regularizer R(S) penalizes the squared distance between the vertices of the mesh M incident to the same
quad. This keeps the size and shape of the quads regular and avoids skewed elements that cause instabilities
in finite-element computations. The regularizer is thus a
simple quadratic penalty that can be described by a sparse
matrix R that has a row for each edge in the mesh with
entries ±1 at the columns of the vertices incident to the
edge. We choose this regularizer for its simplicity.
5

Initialization and Surface Topology

of removing an edge e to
 1 P4
c(e) =

4

i=1

d αi (e),


π 2
2

+ tan(η(e))2

∞

if |η(e)| <
if |η(e)| ≥

π
2
π
.
2

q

|d| ≤

to assign infinity costs to edges that connect strongly bent
triangles (η(e) > π2 ) or edges that would create quads
with straight or reflex angles. The cost favors rectangular
quads to obtain a subdivision surface where finite-element
computations are accurate and numerically stable. We
solve the perfect matching and find the quad mesh with
minimum cost with the Blossom V method [16].
When the input is a point cloud, we use an implicit representation to estimate the topology of the surface from
the point samples and only then extract a control mesh.
Estimating first the topology of the surface with an implicit representation avoids explicitly handling topology
changes during the compression process. We use the popular Poisson reconstruction [15] for simplicity but other
techniques like voxel hashing work well. From this implicit representation, we extract a closed triangular mesh
with marching cubes [21] and create a compact quad mesh
with the collapse and blossom techniques described for
meshes.

4

E(S|S0 ) ≥ E(S) that coincides with E at S0 . We derive an upper envelop of the data and tangent fits from the
inequality
q
q
q−2
q
|d0 | d2 + (1 − )|d0 |
2
2

∀d0 6= 0, q ∈ [1, 2].

The second step of the MM algorithm minimizes the
upper envelop the upper envelop and drives the value of
the original function downwards by solving
S m+1 ← min E(S|S m ).
S

In our case, the envelop
E(S|S m ) =

n
X

wjm DS (pj )2 +

j=1

2
X

2
m
αm
ij htj , ∂i Φ (uj )i + βR(S)

i=1

defines a weighted least-squares problem where the
m
weights wjm and αij
are determined by the MM principle
wjm =

q
DS m (pj )q−2
2

m
αij
=

αq
q−2
|htj , ∂i Φm (uj )i| .
2

To simplify the minimization step of MM, we approximate
DS m (·)2 with a quadratic function that parametrizes the
squared distance to a surface in terms of its geometry as
follows:


>
d · τ2 τ2>
2
> d · τ1 τ1
>
DS (x) ≈ (x − p)
+
+ νν (x − p),
d + ρ1
d + ρ2

Efficient Optimization

Our optimization algorithm exploits the properties of subdivision surfaces, namely, the compact support of the
basis functions and the ability to analytically evaluate
its geometry, with a sequential quadratic program that
is more efficient than the gradient-based algorithms proposed by [5, 12, 13]. We solve (9) as a sequence of convex
problems

where p is a point on the surface S close to x (but not
necessarily the projection of x onto S), d := kx − pk
is the signed distance from x to this point, and ρi , τi , ν
are the principal curvature radii, the principal curvature
directions and the normal to the surface at p. This approximation is due to [26] and coincides with the second order
Taylor approximation within the radius of curvature of the
v m+1 ← min v > (Qm v − b)
(3) surface at p, where it can be relaxed into a positive definite
v
quadratic form by taking the absolute values of d, ρ1 , ρ2 .
that approximate the original energy around the current
Incorporating this approximation into the upper envelop
surface estimate by a least-squares problem. The approx- corresponds to using a quasi-Newton algorithm for optiimation is derived from the surface geometry, while the mizing jointly the control vertices of the subdivision sursparsity of Qm results from compactness of the subdivi- face and the parametric coordinates of the surface point
sion basis.
closest to x [32, 20]. The approximation let us thus work
We follow a Majorize-Minimize (MM) principle to min- with approximate correspondences between surface samimize E(S) by iterating two steps until convergence. The ples and input points, instead of explicitly finding the corfirst step finds an upper envelop of the objective function respondences by solving a larger optimization problem.
6

It also makes our model robust to erroneous correspondences by taking into account their approximate nature in
the objective function. The approximation, however, can
affect the majorizing property and puts our algorithm into
the category of sequential quadratic programs instead of
MM algorithms.
The combination of the MM upper envelope with the
quadratic approximation of the squared distance to S reduces to weighted least-squares problems in the vertices
V (c.f. Appendix) that we solve efficiently with conjugate gradient solver warm-started with the solution of the
previous MM minimization step.
After optimizing V at iteration m, we update the parameter values U where we approximate the squared distance
to the surface by sampling the surface S m = Φm (TM ) uniformly and creating a kd-tree from these samples to find
the parameter of the sample closest to each input point.This
results in an efficient algorithm that avoids the slow convergence rates of coordinate descent algorithms and the
large optimization problems of the Levenberg-Marquard
solvers advocated by [5, 12, 13].

5

defines a weak formulation of ∆ for Sobolev functions
g ∈ H 1 (S). The Galerkin method uses this identity to
define a discretization of ∆ in the ansatz space Hn1 (S)
spanned by Sobolev functions Φ1 , . . . , Φn ∈ H 1 (S).
Using Hn1 (S) as test space, g, h and ∆g are represented
by vectors α, β, γ ∈ Rn as follows
g=

n
X

αi Φi

i=1

h=

n
X

∆g =

i=1

n
X

γ i Φi

i=1

and the weak formulation of ∆ in Hn1 (S) reads
hγ, βiD0 = hα, βiD1 ∀β ∈ Rn .

(4)

The Laplace-Beltrami operator in Hn1 (S) thus depends on
the mass and stiffness matrices D0 , D1 ∈ Rn×n whose
entries
Z
Z
(D0 )ij =
Φi · Φj dx (D1 )ij =
h∇Φi , ∇Φj i dx.
S

S

compute the scalar product of functions and tangent vectors in the surface S.
From (4), the Laplace-Beltrami operator in Hn1 (S) can
be parametrized by the matrix D2 := D0−1 D1 . This matrix has a non-negative spectrum and its eigenvector discretizes the eigenfunctions of the operator ∆ in Hn1 (S).
n
Indeed, if v ∈ R
Pnis an eigenvector of D2 with eigenvalue
λ, then gv = i=1 vi Φi is an eigenfunction of ∆ with
eigenvalue λ in Hn1 (S).
To particularize this construction to functions defined
on a subdivision surface S, we consider the mapping

Smooth Shape Analysis

We have showed how to obtain smooth shapes from nonsmooth input data like point clouds or meshes. We now
show how to do smooth shape analysis with a subdivision
surface S. To this purpose, we focus on the computation
of the Laplace-Beltrami operator and its eigenfunctions
and use them to match shapes, compute shape descriptors,
and approximate geodesics.

X(u, f ) =

5.1

β i Φi

n
X

Φi (u, f )vi

i=1

Laplace-Beltrami-Operator

that parametrizes the surface as a linear combination of
Given a smooth shape S, the Laplace-Beltrami-Operator n limit basis functions Φ1 , . . . , Φn and n mesh vertices
∆ maps a twice continuously differentiable function v1 , . . . , vn . Similarly, we can parametrize Sobolev funcctions over the surface g ∈ H 1 (S) by n scalar values
g ∈ C 2 (S) to a continuous function
g1 , . . . , gn ∈ R with a mapping
∆(g) : S → R,
g : S →R

n
X

n
X

(5)

If g satisfies Neumann boundary conditions, then Stokes’
Φi (u, f ) · vi 7→
Φi (u, f ) · gi
(6)
theorem
i=1
i=1
Z
Z
that depends on each basis function Φi twice, once to
∆ g · h dx =
h∇ g, ∇ hi dx
describe
the function g and once to describe the function
S
S
7

Mesh 3,400 vertices

Mesh 27,984 vertices

Subdiv. Surf. 1700 vertices

Wave-kernel signatures

Figure 3: 12th ∆-eigenfunction and wave-kernel signature of the a surface point in red computed with fine (blue line)
and coarse triangular meshes (red line) and a subdivision surface (yellow). The subdivision surface obtains a wake
signature comparable to the fine mesh while the signature of the coarse mesh lacks detail.
domain S. As a result, the gradient of a surface function
is
∇g(X(u, f )) = JX(u, f )(G−1 (u, f )

n
X

5.2

Geodesic Distances

The classic approach to computing the geodesic distance
χ to a surface point S(u0 ) is to solve the eikonal equation

∇Φi (u, f ) · gi ),

i=1

k∇χk = 1

where JX is the Jacobian of the surface mapping X with
respect to u.
When we use this parametrization to compute the mass
and stiffness matrix, we need to evaluate the first fundamental form of S at a point x ∈ S explicitly. Let f ∈ IF
be an arbitrary face index, then the first fundamental form
at surface point X(u, f ) is
G(u, f ) = JX(u, f )> JX(u, f ).

subject to boundary conditions χ(u0 ) = 0. This results
in a nonlinear hyperbolic PDE that is computationally expensive to solve. To alleviate the computational burden,
[7] proposes a method to approximately solve (8) in terms
of the Heat equation. This method first computes a function f ∈ H 1 (S) whose gradient is parallel to ∇χ by
integrating the Heat flow f˙ = ∆f . It then approximates
∇f
∇χ by the vector field X = − k∇f
k by observing that
∇χ is parallel to ∇f and has unit norm as a result of
(8). RThe method then finds an approximate χ by minimizing S k∇χ − Xk2 through its Euler-Lagrange equation
∆χ = ∇ · X. As all the computations are formulated in
terms of the Heat kernel ∆, the method can be directly
use to compute approximate geodesics in our subdivision
surface framework. We refer to [7] for the details of the
method.

(7)

We compute the surface integrals in the mass and stiffness
matrices by pulling-back the integrand to the reference
quad as follows:
Z
X Z
p
g(s) ds =
g ◦ X(u, f ) det G(u, f ) du
S

f ∈IF

(8)



We can compute the integral over the quad with arbitrary
accuracy with high-order quadrature approximations. Experimentally, we obtain reliable mass and stiffness matrices with simple 3 × 3 Gaussian quadrature.
Note that when the subdivision surface S is derived
from a control mesh with n vertices, D0 , D1 , D2 ∈ Rn×n .
An eigenfunction gv of ∆ is then obtained by solving an
n × nP
generalized eigen problem and by parameterizing
n
gv = i=1 vi Φi . In our figures, we use the refined Mk
mesh and the refined eigenfunction gvk for visualization.
That is, gvk is not derived from the huge mass and stiffness
matrices of Mk but from the projection of gv ∈ C 1 (S)
onto the mesh Mk .

6

Related Work

Subdivision surfaces are a standard representation of
shapes in computer graphics and animation, but they have
largely been neglected in vision due to the challenges of
working with noisy measurements that favor convex formulations and implicit representations.
Computer graphics techniques fit a subdivision surface
by minimize the sum of squared distances between the
8

subdivision surface and the vertices of a clean and regular
mesh [11, 22]. Research has focused into efficient ways
of minimizing this non-convex energy by gradient-based
techniques [35], quasi-Newton methods [23, 32, 20], and
sequential convex programing [23, 17]. A very effective
approach iteratively approximates the squared distance to
the surface with a quadratic penalty [26] and minimizes
the resulting convex energy. Our method method adapts
this strategy to the robust energies necessary to fit a surface
to noisy data and outperform the squared distance penalty
of [23] in three points: 1) We use robust data terms to
account for noise and outliers in the input data. A robust
energy also reduces the effects of errors in the correspondences (between input samples and the surface point closest them) and the quadratic distance approximation. This
is improves the accuracy of the surface independently of
the noise present in the data, as our experiments show.
2) We use a second-order approximation to the squared
distance function, instead of the first order of [23], to take
into account curvature information and fit better the smallscale structures of the surface. 3) We penalize deviations
of the tangent space to improve the appearance of sharp
creases.
Computer Vision, on the other hand, has focused on
developing energy models to fit subdivision surfaces to
noisy data [30, 5, 12, 13], but has ignored the impact of
outliers and the optimization on the accuracy of the surface. [12] fits a sub vision surface to range data by solving
a weighted least-squares problem that assumes optimal
correspondences between input point and domain parameters. This assumption is removed in [30, 5, 13] by explicitly optimizing the correspondences with least-squares
models that preprocess the data to eliminate outliers and
can only handle gross outliers. These methods also assume the topology of the surface to be know and initialize
it with a manually-designed control mesh [30, 5] or a
sphere [13]. The sphere initialization is flexible but oversmooths the surface because to prevents self intersections
as surface evolves. In terms of optimization, all these
techniques use a Levenberg-Marquardt (LM) algorithm
that estimates jointly the vertices of the control mesh and
the correspondences between the input points and surface
samples. This increases the size of the problem and slows
the optimization because gradient-based techniques, like
LM, only update the surface correspondences to adjacent
faces of the control mesh. Our majorize-minimize (MM)

algorithm places no limits on these updates and can jump
to different valleys of the energy landscape.
Our application of subdivision surfaces to shape analysis builds on the work of [14, 3, 33, 10] that investigate
the accuracy of computing mass and stiffness matrices on
a subdivision surface, but do not consider the estimation
of wave-kernel signatures, geodesics, or shape matches in
the subdivision framework.

7

Experiments

Experiments on Surface Compression
Our experiments compare our subdivision representation
to triangular meshes with two different settings: reconstructing a surface from noisy Kinect pointclouds, and
approximating a fine artifact-free mesh with small-scale
structures. For all our experiments we set q = 1.3 and
manually choose α, β to fit each shape and model ( either
our model or the state-of-the-art methods we compare to ).
When the input is a mesh, we set these parameters to minimize the Haussdorff distance between the subdivsision
surface and the input mesh.
Figures 1, 4 show that our subdivision surface reconstructs a surface from a noisy pointcloud as accurately as
the mesh obtained by the state-of-the-art reconstruction
technique [15] but reduces the number of representation
variables to 2% of the mesh size. Our subdivision surfaces
preserve the smooth regions and sharp corners (torso and
hands of Figure 4), and texture details (human hair in Figure 1) of the fine mesh that are lost in a triangular mesh of
the same size obtained by quadratic edge collapse.
For all our experiments with meshes, we use the highresolution versions of the TOSCA dataset and fitted a subdivision surface with 1700 and 1000 vertices. Tables 1–2
report the execution time and the Hausdorff distance between the input mesh and our compressed subdivision representation for 72 high-resolution TOSCA meshes. Figures 5–7 show the subdivision surfaces fitted to the input
meshes with the different models for a qualitative comparison.
Comparison to State-of-the-Art Models When the input is a noisy pointcloud, our surfaces avoid the artifacts of
[23] because our energy model is robust to outliers (Figure
4). Robustness also guards us against erroneous corre9

Point cloud
139553 points

Poisson surface
20932 vertices

Point cloud
107534 points

collapsed Poisson
250 vertices

Poisson surface
12928 vertices

[23] subdiv. surf.
250 vertices

collapsed Poisson
250 vertices

our subdiv. surf.
250 vertices

[23] subdiv. surf.
250 vertices

recognizable
picture

our subdiv. surf.
250 vertices

Figure 4: Surfaces reconstructed from a Kinect point cloud (columns 1) with a state of the art Poisson surface
reconstruction method [15] (column 2) and with the proposed subdivision surface (columns 5). Our technique
reconstructs and represents the surface with an accuracy comparable to a high-resolution mesh while compressing the
representation to less than 2%. We reconstruct the nozzle of the bottle with higher accuracy thanks to our robust emery
model. The subdivision surface of [23] (columns 4) shrinks the statue’s V-shaped torso while it oversmooths its hands
and produces artifacts at the base of the bottle where samples are irregularly spaced and overlap.
spondences in thin structures of clean meshes. We also fit
better small curved structures because we use a secondorder distance approximation, instead of the first-order of
[23], that takes into account curvature information. This
is highlighted in the representation of mouths, eyes, and
ears of Figures 5–7.

Our fitting technique is moderately slower than the
quadratic model of [23] as a result of our robust energy, as
expected. This speed loss is compensated by our improved
surface quality, see Table 1.
Effects of the different terms in our energy model:
Our model always benefits from the robust term while the
tangent-fit is slightly detrimental in rare cases, see columns
Comparison to State-of-the-Art Optimization We
1–3 in Table 1.
compare the subdivision surfaces obtained by optimizing our energy model with the proposed algorithm and
by employing Levenberg-Marquardt to minimize the MM Experiments on Shape Analsysis
majorizer with respect to the mesh vertices and the correspondence parameters. Our optimization is an order of Figure 3 compares the wave kernel signature [1] at a
magnitude faster and less sensitive to local minima than surface point wx with three different discretizations: a
the LM approach proposed by [30, 5, 12, 13]. As a result, fine triangular mesh, a coarse triangular mesh, and our
our surfaces evolve further away from their initialization subdivision surface. This signature is a shape descriptor
to reproduce the small-scale details of the input meshes that assigns to each point x ∈ S a function wx that depends
on the value of the ∆-eigenfunctions and their eigenvalues
(ears and mouths of Figure 5).
10

11

TOSCA mesh, 52565 vertices

[23] subdiv. surface

LM solver

our subdiv. surface

zoom on TOSCA mesh

zoom on [23] surface

zoom on LM solver

zoom on our surface

TOSCA mesh, 15768 vertices

[23] subdiv. surface

LM solver

our subdiv. surface

zoom on TOSCA mesh

zoom on [23] surface

zoom on LM solver

zoom on our surface

Figure 5: Comparison to state-of-the-art: high-resolution input mesh (column 1) mesh and subdivision surfaces fitted
with [23] (column 2) and our proposed model optimized with Levenberg-Marquard (column 3) and our proposed
optimization algorithm (column 4). Comparing columns 2 and 4, shows the benefits of using a robust q-power distance
function and a second-order, instead of a first-order, approximation of the squared distance to the surface to fit thin
surface structures like the ears or lips. Comparing columns 3 and 4 show the benefits of quadratic sequential program
to avoid small-scale local minima close to the initialization that do not capture the shape details of the nose, lips, or
mesh ears. All the subdivision surfaces have 1700 control vertices and are visualized by refining M0 3 times.

12

TOSCA mesh, 19248 vertices

[23] subdiv. surface

LM solver

our subdiv. surface

zoom on TOSCA mesh

zoom on [23] surface

zoom on LM solver

zoom on our surface

TOSCA mesh, 45659 vertices

[23] subdiv. surface

LM solver

our subdiv. surface

zoom on TOSCA mesh

zoom on [23] surface

zoom on LM solver

zoom on our surface

Figure 6: Comparison to state-of-the-art: high-resolution input mesh (column 1) mesh and subdivision surfaces fitted
with [23] (column 2) and our proposed model optimized with Levenberg-Marquard (column 3) and our proposed
optimization algorithm (column 4). Comparing columns 2 and 4, shows the benefits of using a robust q-power distance
function to fit thin surface structures like the lips or mesh eyelids without artifacts. Comparing columns 3 and 4 show
the benefits of quadratic sequential program to avoid small-scale local minima close to the initialization that do not
capture the shape details of the nose, lips, or mesh ears. All the subdivision surfaces have 1700 control vertices and
are visualized by refining M0 3 times.

Table 1: Hausdorff distance (as a ‰of the bounding box diagonal) between the input mesh and a subdivision surface
with 1000 and 1700 control vertices. Columns 1–3 and 6–8 show the importance of each term in the energy model,
while columns 4–5 and 9–10 compare our algorithm to state-of-the-art model of [23] and to the Levenberg-Marquardt
(LM) optimization advocated by [30, 5, 13].
mesh
cat0
cat1
cat10
cat2
cat3
cat4
cat5
cat6
cat7
cat8
cat9
centaur0
centaur1
centaur2
centaur3
centaur4
centaur5
david0
david1
david10
david11
david12
david13
david6
dog0
dog1
dog10
dog2
dog3
dog5
dog6
dog7
dog8
gorilla1
gorilla14
gorilla5
gorilla8
horse0
horse10
horse14
horse15
horse17
horse5
horse6
horse7
michael0
michael1
michael10
michael11
michael12
michael13
michael15
michael16
michael18
michael19
michael2
michael3
michael4
michael5
michael6
michael7
victoria1
victoria10
victoria17
victoria2
victoria21
victoria24
victoria25
victoria4
wolf0
wolf1
wolf2

q=2
0.180
0.251
0.213
0.436
0.210
0.254
0.219
0.390
0.437
0.233
0.283
0.448
0.508
0.632
0.443
0.545
0.474
0.274
0.403
0.396
0.477
0.409
0.349
0.297
0.302
0.433
0.440
0.360
0.384
0.336
0.365
0.444
0.458
0.497
0.854
0.711
0.793
0.372
0.351
0.370
0.352
0.379
0.416
0.353
0.351
0.309
0.376
0.417
0.665
0.436
0.438
0.440
0.403
0.431
0.417
0.578
0.426
0.414
0.366
0.384
0.474
0.321
0.278
0.370
0.302
0.537
0.311
0.283
0.296
0.429
0.410
0.443

α=0
0.176
0.237
0.201
0.390
0.200
0.247
0.208
0.390
0.391
0.209
0.263
0.432
0.449
0.578
0.446
0.488
0.466
0.250
0.377
0.329
0.434
0.363
0.341
0.295
0.291
0.407
0.378
0.333
0.377
0.327
0.335
0.429
0.435
0.379
0.688
0.568
0.650
0.366
0.323
0.348
0.331
0.356
0.374
0.350
0.333
0.276
0.340
0.373
0.576
0.365
0.371
0.387
0.366
0.371
0.381
0.476
0.356
0.364
0.326
0.363
0.394
0.267
0.275
0.305
0.267
0.449
0.293
0.258
0.280
0.427
0.419
0.463

1000 control vertices
full model
[23]
0.158
0.209
0.218
0.295
0.019
0.246
0.352
0.475
0.183
0.259
0.232
0.303
0.201
0.254
0.359
0.454
0.373
0.468
0.197
0.251
0.248
0.336
0.400
0.490
0.424
0.507
0.553
0.653
0.417
0.489
0.478
0.550
0.441
0.521
0.232
0.273
0.370
0.445
0.324
0.378
0.421
0.479
0.360
0.391
0.321
0.381
0.272
0.333
0.272
0.349
0.398
0.470
0.369
0.469
0.311
0.425
0.362
0.453
0.309
0.391
0.328
0.408
0.413
0.496
0.408
0.523
0.398
0.457
0.719
0.885
0.546
0.689
0.621
0.775
0.355
0.431
0.315
0.375
0.343
0.422
0.336
0.399
0.352
0.414
0.364
0.430
0.340
0.411
0.327
0.398
0.266
0.320
0.333
0.396
0.346
0.418
0.547
0.681
0.344
0.431
0.364
0.422
0.369
0.449
0.349
0.416
0.346
0.422
0.367
0.433
0.470
0.537
0.346
0.413
0.350
0.415
0.313
0.384
0.349
0.411
0.365
0.450
0.281
0.330
0.252
0.314
0.309
0.357
0.243
0.322
0.457
0.532
0.277
0.354
0.233
0.312
0.267
0.317
0.414
0.477
0.389
0.455
0.427
0.501

LM solver
0.237
0.327
0.251
0.572
0.263
0.345
0.286
0.443
0.483
0.261
0.356
0.589
0.526
0.733
0.597
0.651
0.563
0.402
0.645
0.453
1.004
0.465
0.647
0.343
0.424
0.624
0.539
0.421
0.463
0.449
0.455
0.628
0.598
0.753
1.348
0.880
1.080
1.252
1.223
0.457
0.478
0.448
0.484
0.424
0.435
0.392
0.480
0.618
0.879
0.507
0.698
0.639
0.552
0.594
0.517
0.647
0.564
0.600
0.454
0.453
0.645
0.365
0.567
0.436
0.363
1.168
0.416
0.377
0.390
0.510
0.517
0.570

q=2
0.113
0.160
0.132
0.272
0.128
0.168
0.139
0.260
0.268
0.156
0.179
0.305
0.323
0.402
0.306
0.336
0.319
0.174
0.275
0.218
0.313
0.257
0.238
0.185
0.195
0.269
0.244
0.214
0.273
0.209
0.221
0.288
0.288
0.340
0.599
0.446
0.530
0.250
0.226
0.246
0.230
0.257
0.260
0.233
0.238
0.185
0.246
0.267
0.443
0.268
0.283
0.283
0.263
0.266
0.259
0.358
0.262
0.264
0.222
0.248
0.287
0.194
0.165
0.203
0.164
0.322
0.188
0.152
0.174
0.331
0.298
0.330

α=0
0.103
0.155
0.123
0.256
0.127
0.166
0.131
0.257
0.253
0.145
0.171
0.303
0.321
0.403
0.310
0.346
0.325
0.171
0.262
0.216
0.293
0.233
0.227
0.186
0.184
0.264
0.246
0.207
0.253
0.205
0.206
0.280
0.287
0.274
0.486
0.380
0.449
0.246
0.225
0.248
0.230
0.250
0.264
0.237
0.233
0.189
0.226
0.256
0.409
0.262
0.267
0.263
0.259
0.255
0.254
0.325
0.249
0.249
0.220
0.234
0.273
0.176
0.166
0.190
0.157
0.293
0.173
0.144
0.167
0.334
0.315
0.345

1700 control vertices
full model
[23]
0.103
0.123
0.148
0.192
0.121
0.155
0.240
0.300
0.126
0.155
0.154
0.195
0.129
0.166
0.240
0.313
0.246
0.304
0.139
0.168
0.166
0.204
0.285
0.347
0.298
0.361
0.381
0.459
0.289
0.352
0.320
0.393
0.304
0.371
0.163
0.187
0.255
0.296
0.199
0.251
0.284
0.351
0.219
0.272
0.216
0.253
0.177
0.211
0.180
0.225
0.250
0.323
0.232
0.289
0.202
0.249
0.242
0.297
0.193
0.244
0.206
0.247
0.268
0.331
0.279
0.340
0.274
0.312
0.489
0.572
0.383
0.444
0.451
0.538
0.233
0.290
0.216
0.274
0.240
0.298
0.212
0.273
0.239
0.297
0.250
0.314
0.222
0.283
0.224
0.284
0.176
0.210
0.214
0.256
0.240
0.290
0.412
0.475
0.255
0.297
0.255
0.296
0.258
0.294
0.249
0.296
0.245
0.285
0.243
0.294
0.318
0.365
0.233
0.279
0.239
0.282
0.206
0.251
0.227
0.265
0.258
0.311
0.166
0.205
0.156
0.202
0.174
0.222
0.150
0.189
0.281
0.337
0.169
0.205
0.137
0.174
0.160
0.197
0.320
0.355
0.293
0.335
0.327
0.363

LM solver
0.141
0.192
0.184
0.366
0.165
0.209
0.165
0.313
0.444
0.189
0.238
0.357
0.374
0.571
0.361
0.394
0.422
0.203
0.502
0.283
0.410
0.381
0.304
0.288
0.286
0.381
0.295
0.303
0.402
0.266
0.391
0.389
0.362
0.449
0.776
0.706
0.863
0.317
0.294
0.346
0.281
0.342
0.370
0.309
0.284
0.257
0.322
0.384
0.904
0.394
0.382
0.427
0.459
0.449
0.363
0.527
0.278
0.386
0.342
0.325
0.387
0.298
0.239
0.310
0.240
0.449
0.329
0.236
0.243
0.431
0.367
0.471

Table 2: Execution time (in seconds) to fit a subdivision surface with 1000 and 1700 control vertices. Columns 1–3
and 6–8 show how each additional term in the energy model slows the optimization, while columns 4–5 and 9–10
compare our technique to state-of-the-art model of [23] and to the Levenberg-Marquardt (LM) optimization algorithm
advocated by [30, 5, 13].
14
1000 control vertices
1700 control vertices
mesh
cat0
cat1
cat10
cat2
cat3
cat4
cat5
cat6
cat7
cat8
cat9
centaur0
centaur1
centaur2
centaur3
centaur4
centaur5
david0
david1
david10
david11
david12
david13
david6
dog0
dog1
dog10
dog2
dog3
dog5
dog6
dog7
dog8
gorilla1
gorilla14
gorilla5
gorilla8
horse0
horse10
horse14
horse15
horse17
horse5
horse6
horse7
michael0
michael1
michael10
michael11
michael12
michael13
michael15
michael16
michael18
michael19
michael2
michael3
michael4
michael5
michael6
michael7
victoria1
victoria10
victoria17
victoria2
victoria21
victoria24
victoria25
victoria4
wolf0
wolf1
wolf2

q=2
32.7
37.1
24.7
23.5
28.1
30.6
18.7
24.9
30.1
34.7
22.6
20.7
12.0
16.0
18.3
16.9
14.3
59.5
58.8
44.7
68.1
63.4
61.0
50.9
20.7
22.9
21.4
28.1
21.1
31.7
35.4
21.0
20.4
36.8
32.8
37.9
46.5
29.6
23.1
17.5
17.5
16.9
17.2
23.2
28.0
63.1
64.5
54.4
44.1
79.7
53.7
66.4
46.8
59.0
53.6
55.3
51.8
48.6
77.6
63.6
55.4
49.7
41.8
39.2
53.5
35.1
40.7
49.9
46.5
7.4
6.5
7.0

α=0
34.7
41.2
36.5
33.9
36.6
43.4
38.8
43.9
23.9
42.1
35.7
16.5
21.9
15.9
14.1
24.0
17.0
59.4
56.4
47.1
58.2
52.0
45.7
42.6
25.4
27.4
30.5
27.2
34.0
28.2
31.3
28.3
27.4
33.0
41.3
47.2
44.5
22.5
26.7
26.2
21.6
22.6
28.9
34.8
24.0
54.9
55.0
63.1
80.9
75.2
70.3
61.5
50.9
60.1
52.5
54.4
57.4
47.5
59.3
39.5
65.5
49.4
36.8
64.7
38.0
69.0
46.2
52.9
43.7
12.1
8.5
8.0

full model
56.0
44.0
34.9
43.6
47.4
46.7
41.6
40.7
40.9
35.9
56.6
21.5
16.5
26.0
25.9
26.9
28.1
86.9
70.1
72.1
128.5
61.7
86.8
71.3
48.8
30.1
39.6
38.9
32.3
44.8
53.0
38.6
41.8
61.2
52.9
55.1
59.6
46.4
47.0
31.8
34.1
21.4
30.5
25.4
38.6
74.1
86.9
108.5
115.8
100.4
103.5
121.4
65.4
98.9
58.6
55.5
95.9
76.3
106.4
68.5
96.4
52.5
40.3
55.6
45.8
80.2
84.5
57.7
60.4
11.5
9.1
9.7

[23]
20.8
12.9
18.2
21.3
16.4
16.5
17.4
18.0
14.3
17.4
19.1
8.2
8.9
9.9
8.5
8.8
8.8
39.3
34.9
37.6
37.3
38.8
37.4
30.7
14.0
15.7
12.8
14.3
16.1
20.4
18.8
15.2
14.8
23.6
21.0
38.3
19.9
12.9
12.6
13.4
12.0
13.9
13.2
11.2
11.1
37.4
47.2
38.4
29.3
37.2
39.0
34.5
47.8
32.9
34.6
35.4
36.2
33.7
27.4
26.5
35.0
29.3
24.4
27.7
27.6
27.0
22.7
23.6
27.5
6.8
8.2
5.7

LM solver
408.0
345.2
403.3
450.0
333.7
344.1
547.1
365.1
402.6
267.5
361.4
165.4
279.3
209.2
192.4
267.1
254.2
602.6
447.5
589.8
683.6
528.8
807.0
592.7
325.0
281.2
283.4
251.6
303.0
271.4
293.2
427.1
304.5
531.9
519.9
565.8
436.8
98.6
120.0
261.4
489.1
233.1
322.4
250.3
194.6
632.1
588.5
476.4
583.0
607.9
553.5
690.1
495.6
311.1
473.3
764.9
452.0
569.7
468.1
389.0
508.0
439.4
527.8
582.8
448.4
321.3
366.5
308.4
393.3
64.2
78.5
52.9

q=2
24.6
34.5
30.3
31.3
41.5
28.6
29.6
28.5
26.6
27.4
35.5
16.2
23.8
19.2
18.8
24.2
18.9
64.8
60.1
59.0
59.2
63.8
58.8
54.2
23.2
39.2
31.0
29.3
28.3
34.3
24.8
44.0
28.5
45.8
41.6
28.3
37.8
22.0
34.9
30.9
20.1
24.9
29.3
30.2
19.7
60.2
56.3
69.4
54.3
77.8
52.2
62.5
80.9
55.9
71.0
61.2
52.8
76.4
57.1
68.5
62.0
30.3
59.3
45.0
48.0
37.5
27.6
47.3
41.6
10.6
17.0
8.0

α=0
44.5
43.3
50.0
44.3
51.8
34.8
47.3
48.3
50.9
48.7
29.7
22.1
22.0
24.9
21.9
28.3
25.0
56.1
55.2
59.8
76.5
53.2
58.4
60.7
31.1
39.6
31.4
45.7
27.1
48.8
45.5
47.4
45.4
48.3
72.1
55.3
64.6
42.0
36.7
45.8
24.3
35.5
35.0
44.5
29.8
54.9
56.5
51.8
56.9
61.6
51.5
53.3
69.7
69.6
65.9
47.2
49.6
52.2
53.4
67.0
57.8
38.9
60.4
43.9
57.0
42.9
47.2
79.1
46.4
16.7
22.8
11.6

full model
36.0
52.6
38.7
42.3
32.4
54.7
33.4
41.1
51.4
44.1
47.2
24.0
33.7
24.1
26.4
27.9
25.9
63.4
63.0
79.9
100.4
74.9
108.8
63.2
39.6
41.8
32.8
48.0
34.1
34.8
49.4
45.0
37.6
64.4
65.6
63.8
69.5
34.0
32.1
45.9
40.4
36.8
37.5
38.8
28.0
76.2
63.7
76.5
84.3
118.6
71.6
68.3
72.0
78.6
60.5
68.9
60.7
70.8
57.7
65.8
62.4
52.5
57.3
72.0
67.9
68.9
76.1
65.5
58.7
9.1
15.8
16.7

[23]
18.6
15.1
37.4
21.0
18.6
20.6
17.1
16.1
20.3
20.8
17.0
15.8
11.7
10.4
10.2
11.3
9.8
40.2
36.0
40.2
33.1
45.6
43.0
33.7
14.6
16.8
15.0
18.8
15.4
17.5
18.5
22.2
15.7
34.1
25.3
24.3
26.9
13.0
11.8
13.2
16.5
13.1
16.5
13.1
12.7
40.9
44.5
29.3
31.3
40.5
39.4
38.2
56.1
34.9
38.0
35.7
38.3
36.3
38.9
37.6
40.5
27.5
30.3
30.6
27.6
24.5
26.8
28.5
23.7
5.7
16.6
8.1

LM solver
346.5
367.5
544.6
399.9
382.5
384.0
372.9
332.3
468.0
378.9
468.4
185.1
374.1
171.3
258.7
307.5
197.1
473.5
799.9
611.4
713.5
622.2
832.3
551.9
321.8
368.4
348.5
323.8
311.2
207.0
378.7
355.7
366.5
416.0
620.7
837.3
695.7
234.2
302.6
406.2
252.8
316.4
346.2
254.3
290.6
599.8
909.7
517.2
762.0
573.2
546.4
443.0
571.5
472.5
535.3
476.2
658.7
576.4
1117.8
439.6
489.5
381.9
405.8
534.7
584.8
444.0
549.3
412.3
478.4
114.0
137.9
106.7

TOSCA mesh, 52565 vertices

[23] subdiv. surface

LM solver

our subdiv. surface

zoom on TOSCA mesh

zoom on [23] surface

zoom on LM solver

zoom on our surface

Figure 7: Comparison to state-of-the-art: high-resolution input mesh (column 1) mesh and subdivision surfaces fitted
with [23] (column 2) and our proposed model optimized with Levenberg-Marquard (column 3) and our proposed
optimization algorithm (column 4). Comparing columns 2 and 4, shows the benefits of using the second-order, instead
of a first-order, approximation of the squared distance to the surface to capture the curvature details of small structures
like the mesh ears. Comparing columns 3 and 4 show the benefits of quadratic sequential program to avoid small-scale
local minima close to the initialization that do not capture the shape details of the nose, lips, or mesh eyes. All the
subdivision surfaces have 1700 control vertices and are visualized by refining M0 3 times.
at x. The point signature of coarse mesh is different from racy at the tails of some cats, but this is compensated by
the signature of fine mesh, which is well matched by the a surface representation and matching problem of 10%
subdivision surface.
of the original size. To visualize the results, we show a
template shape with a color map that interprets the three
The functional map framework [25] formulates the spatial coordinates of every point as RGB color informaproblem of matching points of two shapes as matching tion and transfer this color map to each target shape with
smooth functions over them. The problem is invariant the matching estimated with functional maps [25]. In this
to isometric transformations by representing smooth func- experiment, using a coarse triangular mesh to do the shape
tion with ∆-eigenfunctions and by finding a linear map be- matching simply fails.
tween these eigenfunctions. Figure 8 shows the computed
matching with respect to two different shape representaFigures 1–10 compare geodesics computed by the Heat
tions: a fine mesh with 27894 vertices and a subdivision method [7] with different surface representations fitted to
surface with 1700 control vertices. Our experiments show the Kinect pointclouds of Figures 1 and 4. As a proxy
comparable matching results with a slight loss of accu- for ground-truth, we use a high-resolution triangular mesh
15

Fine mesh

Matched shape

Matched shape

Matched shape

Subdivision surface

Matched shape

Matched shape

Matched shape

Figure 8: Shape matching results of the functional map framework [25]. After the matching is determined, the color
information of the target shape is transferred to the matched shape with the mapping estimated by functional maps[25].
obtained by the Poisson reconstruction method [15] and
compare it to two compressed surface representations of
the same size: a low-resolution triangular mesh obtained
by quadratic edge collapse [9] of the high-resolution mesh
and our subdivision representation. Our subdivision surface is comparable to the high-resolution mesh and can
be used to compressed both the shape and its geodesics
functions. In contrast, the low-resolution triangular mesh
provides a discretization too coarse to represent the small
scale details of the surface and its geodesic. To highlight the loss of small-scale structure on the geodesic, we
represent the level lines of the geodesic at high and low
resolution and demonstrate how the geodesic estimated
from the coarse triangular mesh is unable to capture the
high-frequency information of the geodesic (as estimated
from the high-resolution mesh) while our surface representation captures both low and high frequency information
accurately.

representations and develop standard shape analysis operators with it.
Our surface fitting model is robust to noise and outliers and can handle both clean input meshes and noisy
point cloud obtained from real data. The resulting subdivision surfaces let us represent smooth shapes at high
accuracy with a fraction of the variables used in triangular meshes and use the same subdivision schemes to
parametrize functions over the surface and compute tangent spaces or curvatures analytically. Our experiments
show how the smoothness of the representation reduces
the size of the Laplace-Beltrami operator and its eigenfunctions without loss of accuracy and allows us to compute wave-kernel signatures, approximate geodesics, and
shape matches in a compact representation.

8

We optimize the objective energy

Conclusions

A

We have presented a method to compress surfaces with
a representation well suited for shape analysis. To this
purpose, we fit a subdivision surface to low-level surface

Appendix: Quadratic Solver

min E(S) = dist(S, P) + αT (S, T) + βR(S) , (9)
S
| {z } | {z }
| {z }
Point Fit

16

Tangent Fit

Regularization

triangular mesh from [15]
12928 vertices

approximate geodesic
g

30 geodesic level lines
cos(30πg)

5 geodesic 5 level lines
cos(5πg)

edge-collapsed mesh [15]
250 vertices

approximate geodesic
g

30 geodesic level lines
cos(30πg)

5 geodesic 5 level lines
cos(5πg)

subdivision surface
250 control vertices

approximate geodesic
g

30 geodesic level lines
cos(30πg)

5 geodesic 5 level lines
cos(5πg)

Figure 9: Computation of approximate geodesics by the Heat method [7] with different surface representations fitted
to a Kinect pointcloud. Column 1: representation of the surface with a high-resolution mesh obtained with Poisson
reconstruction [15] (row 1), a low-resolution triangular mesh obtained by edge collapse of the high-resolution mesh
(row 2), and with our subdivision surface (row 3). Column 2: approximate geodesic g estimated with the Heat method
[7] on each surface representation. Columns 3 and 4: level lines of the geodesics visualized as cos($g), with $ = 30π
in column 3 and $ = 5π in column 4. Our subdivision surface is comparable to the surface obtained with the Poisson
reconstruction method [15] and can faithfully represent geodesics at both low and high resolution even though it has
2% of its vertices. In contrast, compressing the state-of-the-art mesh with the quadratic edge collapse method of [9]
looses all the small scale details of the surface and its geodesics. This is highlighted by representing the level lines of
the geodesic on the mesh at high resolution (column 3) where we show how the discretization associated with a coarse
triangular mesh representation cannot estimate the high-frequency information of the geodesic accurately.
that fits a subdivision surface S to a set of point P and eration m we substitute our robust energy by the weighted
normal data T by solving a sequence of convex problems least-squares problem
V m+1 ← min V > (Qm V − bm )
V

(10)
min
U,V

N
X

2

wjm kpj − Φ(uj )k +

2
X

2

m
αij
htj , ∂i Φ(uj )i

j=1
i=1
over the vertices V = [v1 , . . . , vn ] of the control mesh M.
2
In this section, we show how to obtain the quadratic ob+ βkRV k
(11)
jective function V > (Qm V − bm ) that approximates the
original energy at iteration m.
over the mesh vertices V and the correspondence paramFollowing a majorize-minimize (MM) principle, at it- eters u1 , . . . , uN ∈ TM . The parametrization of the sub-

17

triangular mesh from [15]
20932 vertices

approximate geodesic
g

30 geodesic level lines
cos(30πg)

5 geodesic 5 level lines
cos(5πg)

edge-collapsed mesh [15]
250 vertices

approximate geodesic
g

30 geodesic level lines
cos(30πg)

5 geodesic 5 level lines
cos(5πg)

subdivision surface
250 control vertices

approximate geodesic
g

30 geodesic level lines
cos(30πg)

5 geodesic 5 level lines
cos(5πg)

Figure 10: Computation of approximate geodesics by the Heat method [7] with different surface representations fitted
to a Kinect pointcloud. Column 1: representation of the surface with a high-resolution mesh obtained with Poisson
reconstruction [15] (row 1), a low-resolution triangular mesh obtained by edge collapse of the high-resolution mesh
(row 2), and with our subdivision surface (row 3). Column 2: approximate geodesic g estimated with the Heat method
[7] on each surface representation. Columns 3 and 4: level lines of the geodesics visualized as cos($g), with $ = 30π
in column 3 and $ = 5π in column 4. Our subdivision surface is comparable to the surface obtained with the Poisson
reconstruction method [15] and can faithfully represent geodesics at both low and high resolution even though it has
2% of its vertices. In contrast, compressing the state-of-the-art mesh with the quadratic edge collapse method of [9]
looses all the small scale details of the surface and its geodesics. This is highlighted by representing the level lines of
the geodesic on the mesh at high resolution (column 3) where we show how the discretization associated with a coarse
triangular mesh representation cannot estimate the high-frequency information of the geodesic accurately.
division surface

Φi (uj ) that samples the surface as
Φ(u) =

n
X

vi Φi (u),

xj =

i=1

n
X

vi Φi (uj ) j = 1, . . . , N

X = L · V.

i=1

depends linearly on the mesh vertices and can be described
compactly with an N × n sparse matrix L with Lji = Analogously, we define L1 and L2 for the partial deriva18

tives of the basis functions

After optimizing V at iteration m, we update the parameter
values u1 , . . . , un where we approximate the squared
n
X
distance
to the surface with the matrices D1m , . . . , Dnm .
∂i Φ(uj ) =
Li,jk vk
(12)
In particular, we sample the surface S m = Φm (TM ) unik=1
formly and create a kd-tree from these samples to find
and introduce the auxiliary variables x0ij = ∂i Φ(uj ) that quickly the parameters u ∈ T of the surface sample
j
M
sample the tangent plane of the surface by
closest to each input point pj ∈ P. Although this update
strategy might seem reminiscent of coordinate descent,
Xi0 =Li · V ∈ RN ×3 .
our combination of coordinate updates with the quadratic
We remove the variables uj from the optimization by approximation corresponds to applying a quasi-Newton
approximating the squared distance with the quadratic ex- algorithm [32, 20] to the majorizer to avoid the slow convergence rates of coordinate descent. This removes the
pression in V
correspondence parameters from the optimization but can
n
n
X
X
2
2
affect the majorizing property of (11) and puts our algom
m
wj kpj − Φ(uj )k ≈
wj kpj − (L · V )j kDm ,
j
rithm into the category of quadratic sequential programs
j=1
j=1
instead of majorize-minimize algorithms. For this reason,
n
X
m
> m
wj (pj − xj ) Dj (pj − xj ) we stop the process as soon as the original energy E(·)
≈
does not decrease any more. Our experiments indicate
j=1
that this method provides better results than state-of-thewhere the symmetric positive definite matrix
art methods.
Djm =

d
d
τ1 τ > +
τ2 τ > + νν >
d + ρ1 1
d + ρ2 2

(13)

is determined by the approximation to the squared distance
to the surface S m at pj . In particular, xj is a surface point
close to pj (but not necessarily its projection onto S m ), d
is the distance from pj to this point, and ρ1 , ρ2 , τ1 , τ2 , ν
are the absolute values of the principal curvature radii,
the principal curvature directions, and the normal to the
surface S m at xj . Considering Djm independent of xj at
iteration m, the data fit term defines a quadratic problem.
The tangent fit energy

References

[1] M. Aubry, U. Schlickewei, and D. Cremers. The wave
kernel signature: A quantum mechanical approach to shape
analysis. In IEEE ICCV, pages 1626–1633, 2011. 1, 3, 10
[2] A. Bronstein, M. Bronstein, and R. Kimmel. Numerical
Geometry of Non-Rigid Shapes. Springer, 2008. 1
[3] D. Burkhart, B. Hamann, and G. Umlauf. Iso-geometric
finite element analysis based on Catmull-Clark subdivision
solids. Eurographics Symp. Geom. Process., 29(5):1575–
1584, 2010. 9
N
2
2
N
D
E
XX
XX
2 [4] F. Calakli and G. Taubin. SSD: Smooth Signed Dis2
m
m
αij
tj , (Lm
αij
htj , ∂i Φm (uj )i =
.
tance Surface Reconstruction. Comput. Graph. Forum,
i · V )j
30(7):1993–2002, sep 2011. 2
j=1 i=1
j=1 i=1
[5]
T. J. Cashman and A. W. Fitzgibbon. What shape are
m
is also quadratic in V . Due to non-negative αij
this exdolphins? building 3D morphable models from 2D images.
pression is positive semi-definite. As the regularization
IEEE Trans. Pattern Anal. Mach. Intell., 35(1):232–244,
2
kRV k = V > R> RV is also quadratic and positive semi2013. 6, 7, 9, 10, 13, 14
definite, we can efficiently solve the linear system that [6] E. Catmull and J. Clark. Recursively generated b-spline
characterizes its unique minimum with the conjugate grasurfaces on arbitrary topological meshes. Computer-Aided
dient method. This method is well suited to the problem at
Design, 10(6):350 – 355, 1978. 4
hand, because the matrices that define our quadratic objec- [7] K. Crane, C. Weischedel, and M. Wardetzky. Geodesics in
tive are sparse and we can use the MM iterations to warm
heat: A new approach to computing distance based on heat
start the optimization at iteration m with the solution at
flow. ACM Trans. Graph., 32(3):10, 2013. 1, 2, 3, 8, 15,
iteration m − 1.
17, 18

19

[8] V. Estellers, M. Scott, and S. Soatto. Robust Surface Reconstruction. SIAM J. Imaging Sci., 9(4):2073–2098, jan
2016. 2
[9] M. Garland and P. S. Heckbert. Surface simplification using
quadric error metrics. SIGGRAPH, pages 209–216, 1997.
2, 5, 16, 17, 18
[10] F. Goes, M. Desbrun, M. Meyer, and T. Derose. Subdivision
Exterior Calculus for Geometry Processing. ACM Trans.
Graph. Artic., 35(10), 2016. 9
[11] H. Hoppe, T. DeRose, T. Duchamp, M. Halstead, H. Jin,
J. McDonald, J. Schweitzer, and W. Struetzle. Piecewise
smooth surface reconstruction. Comput. Graph. (ACM).,
pages 294–302, 1994. 9
[12] S. Ili. Using Subdivision Surfaces for 3 âĂŞ D Reconstruction from Noisy Data. Deform. Image Regist. Deform.
Environ., pages 1—-10, 2006. 6, 7, 9, 10
[13] M. Jaimez, T. J. Cashman, A. Fitzgibbon, J. Gonzalezjimenez, and D. Cremers. An Efficient Background Term
for 3D Reconstruction and Tracking with Smooth Surface
Models. IEEE Conf. Comput. Vis. Pattern Recognit., 2017.
2, 6, 7, 9, 10, 13, 14
[14] B. Jüttler, A. Mantzaflaris, R. Perl, and M. Rumpf. On
numerical integration in isogeometric subdivision methods
for {PDEs} on surfaces. Computer Methods in Applied
Mechanics and Engineering, 302:131 – 146, 2016. 3, 9
[15] M. Kazhdan and H. Hoppe. Screened poisson surface reconstruction. ACM Trans. Graph., 32(3):1–13, jun 2013.
2, 6, 9, 10, 15, 17, 18
[16] V. Kolmogorov. Blossom v: a new implementation of a
minimum cost perfect matching algorithm. Mathematical
Programming Computation, 1(1):43–67, July 2009. 6
[17] G. Lavoué, F. Dupont, and A. Baskurt. Subdivision surface
fitting for efficient compression and coding of 3D models.
Vis. Commun. Image Process., 5960(33):59603F–59603F–
12, 2005. 9
[18] B. Levy. Laplace-Beltrami Eigenfunctions: Towards an
Algorithm that Understand s Geometry. IEEE Int. Conf.
Shape Model. Appl. inv ited talk, 2006, 2006. 1
[19] R. Litman and A. M. Bronstein. Learning spectral descriptors for deformable shape correspondence. IEEE Trans.
Pattern Anal. Mach. Intell., 36(1):171–180, 2014. 1
[20] Y. Liu and W. Wang. A revisit to least squares orthogonal
distance fitting of parametric curves and surfaces. Lect.
Notes Comput. Sci., 4975 LNCS(1):384–397, 2008. 6, 9,
19
[21] W. E. Lorensen and H. E. Cline. Marching cubes: A
high resolution 3D surface construction algorithm. ACM
SIGGRAPH Comput. Graph., 21(4):163–169, 1987. 6

20

[22] W. Ma, X. Ma, S. K. Tso, and Z. Pan. Subdivision surface
fitting from a dense triangle mesh. Proc. - Geom. Model.
Process. Theory Appl. GMP 2002, 36:94–103, 2002. 9
[23] M. Marinov and L. Kobbelt. Optimization methods for scattered data approximation with subdivision surfaces. Graph.
Models, 67(5):452–473, 2005. 9, 10, 11, 12, 13, 14, 15
[24] C. Notes and D. Zorin. Subdivision for Modeling and
Animation. New York, 98:1 – 203, 2000. 1
[25] M. Ovsjanikov, M. Ben-Chen, J. Solomon, A. Butscher,
and L. Guibas. Functional maps. ACM Trans. Graph.,
31(4):1–11, jul 2012. 1, 3, 15, 16
[26] H. Pottmann and S. Leopoldseder. A concept for parametric
surface fitting which avoids the parametrization problem.
Comput. Aided Geom. Des., 20(6):343–362, 2003. 6, 9
[27] J. F. Remacle, J. Lambrechts, B. Seny, E. Marchandise,
A. Johnen, and C. Geuzainet. Blossom-Quad: A nonuniform quadrilateral mesh generator using a minimumcost perfect-matching algorithm. Int. J. Numer. Methods
Eng., 89(9):1102–1119, 2012. 5
[28] M. Reuter, F.-E. Wolter, and N. Peinecke. {LaplaceBeltrami} Spectra As ’Shape-{DNA}’ of Surfaces and
Solids. Comput. Aided Des., 38(4):342–366, apr 2006.
1
[29] J. Stam. Exact evaluation of Catmull-Clark subdivision
surfaces at arbitrary parameter values. Proc. 25th Annu.
Conf. Comput. Graph. Interact. Tech. - SIGGRAPH ’98,
pages 395–404, 1998. 4
[30] J. Taylor, J. Valentin, B. Luff, R. Banks, A. Fitzgibbon,
T. Sharp, E. Wood, and J. Shotton. Efficient and Precise Interactive Hand Tracking Through Joint , Continuous
Optimization of Pose and Correspondences. ACM Trans.
Graph. ACM SIGGRAPH, 35(4):143, 2016. 9, 10, 13, 14
[31] B. Ummenhofer and T. Brox. Global , Dense Multiscale
Reconstruction for a Billion Points. Int. Conf. Comput. Vis.
Pattern Recongition, 2015. 2
[32] W. Wang, H. Pottmann, and Y. Liu. Fitting B-spline curves
to point clouds by curvature-based squared distance minimization. ACM Trans. Graph., 25(2):214–238, 2006. 6, 9,
19
[33] A. Wawrzinek and K. Polthier. Computer-Aided Design
Integration of generalized B-spline functions on Catmull
âĂŞ Clark. Comput. Des., 78:60–70, 2016. 9
[34] T. Windheuser, M. Vestner, E. Rodola, R. Triebel, and
D. Cremers. Optimal Intrinsic Descriptors for Non-Rigid
Shape Analysis. Br. Mach. Vis. Conf., page 4, 2014. 1
[35] W. Zheng, P. Bo, Y. Liu, and W. Wang. Fast B-spline curve
fitting by L-BFGS. Comput. Aided Geom. Des., 29(7):448–
462, 2012. 9

