EUROGRAPHICS 2020 / U. Assarsson and D.

Panozzo Volume 39 (2020), Number 2

(Guest Editors)

Spectral Mesh Simplification

Thibault Lescoat1 Hsueh-Ti Derek Liu2 Jean-Marc Thiery1 Alec Jacobson2 Tamy Boubekeur4,1 Maks Ovsjanikov3

1 LTCI, Télécom Paris, Institut Polytechnique de Paris, France 2 University of Toronto, Canada 3 École Polytechnique, France 4 Adobe

Ground truth Uniform Garland & Heckbert 1997 Ours

|| . ||L = 31.38 × 103 || . ||L = 39.35 × 103 || . ||L = 3.50 × 103

|| . ||D = 15.60 × 100 || . ||D = 15.13 × 100 || . ||D = 3.76 × 100

Figure 1: We propose to simplify a mesh using edge collapses while aiming to preserve the input eigenvectors and eigenvalues as much
as possible. While different strategies exist to reduce a mesh (here, from 25,727 vertices to 771 vertices, or 3% of its initial size), such as
enforcing uniform edge lengths or using the Quadric Error Metric [GH97], they do not focus on keeping the spectral properties of the mesh.
Reducing a mesh can be spectrally described using functional maps [OBCS∗ 12], shown here with the output meshes, and which should ideally
be diagonal. We also evaluate functional maps using two norms, the Laplacian commutativity k · kL and the orthogonality k · kD .
Abstract
The spectrum of the Laplace-Beltrami operator is instrumental for a number of geometric modeling applications, from processing
to analysis. Recently, multiple methods were developed to retrieve an approximation of a shape that preserves its eigenvectors
as much as possible, but these techniques output a subset of input points with no connectivity, which limits their potential
applications. Furthermore, the obtained Laplacian results from an optimization procedure, implying its storage alongside the
selected points. Focusing on keeping a mesh instead of an operator would allow to retrieve the latter using the standard cotangent
formulation, enabling easier processing afterwards. Instead, we propose to simplify the input mesh using a spectrum-preserving
mesh decimation scheme, so that the Laplacian computed on the simplified mesh is spectrally close to the one of the input mesh.
We illustrate the benefit of our approach for quickly approximating spectral distances and functional maps on low resolution
proxies of potentially high resolution input meshes.

1. Introduction The lack of a mesh limits the use of coarsening in many downstream
geometry processing tasks.
Triangle meshes remain a predominant representation of 3D sur-
We present the first mesh simplification method intentionally
faces. When the complexity of a given mesh exceeds computational
designed to preserve spectral properties. We propose adapting the
resources, we rely on mesh simplification methods to remove ver-
standard greedy edge-collapse mesh-simplification algorithm with a
tices, edges, and faces. In rendering, efficient simplification methods
novel cost function that measures spectral preservation of a given op-
can dramatically reduce the complexity of a mesh without affecting
erator (e.g., the cotangent Laplacian). Unlike algebraic methods that
its appearance. It is tempting to repurpose appearance-preserving
directly output a reduced operator (i.e., matrix), our method outputs
simplification methods for other geometry processing tasks.
a manifold triangle mesh with 3D vertex positions. Reconstructing
the operator on the output mesh will preserve both the eigenvalues
Unfortunately, appearance-based methods do not preserve the and eigenvectors of the operator on the input mesh.
spectral properties of the important differential operators upon which
much of modern geometry processing is built (see Figure 1). As a Confirmed by a series of experiments, our method preserves spec-
result, solutions computed on such a coarse mesh can be incorrect tral properties nearly as well as purely algebraic methods, while
or misleading. Alternatively, previous coarsening methods that do still outputting an embedded mesh like standard simplification algo-
preserve spectral properties work purely algebraically on the oper- rithms. We demonstrate our approach’s effectiveness for geodesic
ator matrices and do not produce a geometric mesh (see Figure 2). distance approximation and functional maps correspondence.

c 2020 The Author(s)

Computer Graphics Forum c 2020 The Eurographics Association and John
Wiley & Sons Ltd. Published by John Wiley & Sons Ltd.
 T. Lescoat, H.-T. Liu, J.-M. Thiery, A. Jacobson, T. Boubekeur, M. Ovsjanikov / Spectral Mesh Simplification

Ground truth Nasikun et al. 2018 Liu et al. 2019 Garland & Heckbert 1997 Ours

| |: 56,112

|| . ||L = 1.66 × 103 || . ||L = 8.02 × 103 || . ||L = 61.19 × 103 || . ||L = 1.71 × 103
|| . ||D = 2.62 × 100 || . ||D = 1.86 × 100 || . ||D = 42.25 × 100 || . ||D = 1.14 × 100

Figure 2: Reduction from 56,112 to 1,122 vertices (2% of input size). Nasikun et al. [NBH18] approximate the original Laplacian on a subset
of vertices obtained via Poisson-disk sampling. Liu et al. [LJO19] optimize the sampling and the operator, that they define on samples’ 3-rings.
Instead, by outputting a mesh, further processing can use a standard cotangent weighting scheme without knowledge of the reduction step.

2. Background An alternative approach to preserving properties of the opera-

tor built from the input mesh is to directly simplify it as a matrix.
Our method builds on the long history of research in mesh simplifica- Recently, Liu et al. [LJO19] present a state-of-the-art method for
tion and recent developments in spectral coarsening for differential spectral preservation during algebraic coarsening of common op-
operators. We focus the attention of this related work section on erators used in geometry processing. We refer the reader to this
methods directly related in methodology or intention. recent work for a comprehensive review of previous algebraic and
numerical coarsening methods. Notably, Nasikun et al. [NBH18]
Classic methods for mesh simplification are based on preserving
aim at efficiently approximating eigenpairs. Both methods select
the rendered appearance of the geometric surface [SZL∗ 92, PH97,
a subset of points from the input; along with other recent devel-
GH97, HDD∗ 93, CSAD04]. These methods were extended to ac-
opments (e.g., [CBO∗ 19, OS19]), they share a common limitation
count for other signals stored on the mesh beyond geometry includ-
compared to our method: they do not produce a mesh.
ing texture coordinates and colors [CMO97,GH98,Hop99,LFJG17].
Similar to many of these methods, we adapt the per-edge cost We build upon the functional maps [OBCS∗ 12] machinery used
function of the basic greedy edge-collapse approach introduced by Liu et al. to evaluate how well a coarsening preserves spectral
by Garland & Heckbert [GH97]. Instead of optimizing perceptual properties. Importantly, unlike their computationally expensive alge-
metrics, we optimize a spectral metric. Spectral preservation re- braic optimization, our efficient edge-collapse algorithm maintains
lies on maintaining intrinsic properties of the surface (the metric). a manifold triangle mesh.
Previous methods have focused on maintaining extrinsic proper-
ties, such as keeping the coarse mesh within a small envelope
around the input [CVM∗ 96, ZG02] or strictly containing the in- 3. Method
put [SGG∗ 00, SVJ15]. In a rare previous example of spectral mesh Our spectrum-preserving simplification method is made of two
simplification, Li et al. [LFZ15] append modal displacement vec- main building blocks: a simplification algorithm based on edge-
tors for sound simulation as extra dimensions during greedy edge- collapses and a simplification metric which drives the algorithm.
collapse. However, this method preserves only the specific modes The metric associates a cost to any given edge-collapse, ordering
chosen and would not scale beyond a small number of frequen- them dynamically during the simplification.
cies. Our efficient method preserves a large span of low frequency
eigenmodes and corresponding values.
3.1. Input / output
Mesh simplification is closely related to graph reduction. In
this more general and less constrained context, recent works Our method takes as input a manifold triangle mesh M = (V, F ),
have investigated spectral preservation with theoretical guaran- which can optionally contain boundaries, and produces a simplified
tees [KS16, Lou19, LV18]. Liu et al. [LJO19] recently demonstrated mesh M f = (V,
e Fe), with a spectrum as close as possible to M when
superiority over [KS16] when applied to mesh Laplacians from evaluating it using the standard Laplacian operator. Optionally, a
geometry processing. Similarly, coreset selection algorithms aim coarse-to-fine restriction matrix can be produced and used when
to preserve statistics or properties of a larger point set or distribu- computing e.g., functional maps. We also take a unique parameter
tion [HCB16, CS18]. in the form of the number of eigenvectors to preserve.

c 2020 The Author(s)

Computer Graphics Forum c 2020 The Eurographics Association and John Wiley & Sons Ltd.
 T. Lescoat, H.-T. Liu, J.-M. Thiery, A. Jacobson, T. Boubekeur, M. Ovsjanikov / Spectral Mesh Simplification

input, | |: 20,212

|| . ||L = 2.57 × 103

|| . ||D = 2.58 × 100
ours, | |: 808

|| . ||L = 59.12 × 103 Garland & Heckbert φ5 φ6 ... φ35 φ36

|| . ||D = 71.03 × 100 1997, | |: 808

Figure 3: While the very first eigenvectors from both our method (middle row) and QSlim [GH97] (bottom row) are very close to those of the
input (top row), following eigenvectors are more sensitive; badly preserved eigenvectors exhibit patterns dissimilar to their input counterpart.

3.2. Simplification algorithm 3.3. Metric

Given a cost for each edge of M, our simplification algorithm (see While the Quadric Error Metric, introduced in the original work of
Algorithm 1) follows the seminal idea of Garland and Heckbert Garland & Heckbert [GH97], maintain the visual appearance of the
[GH97]: each input edge is pushed to a priority-queue, where the original mesh as much as possible, we propose a new alternative
priority is dictated by our specific cost metric. At all time, the edge metric focused on spectral preservation.
located at the head of the queue is the next best edge to collapse
i.e., with minimum cost. Once popped from the queue, the edge Let L, M ∈ R|V|×|V| be the Laplacian and the diagonal mass ma-
e e
is collapsed, effectively removing one vertex from the mesh and eM
trix of M respectively. Similarly, L, e ∈ R|V|×| V|
denote the Lapla-
resulting in a merged vertex positioned to optimize the metric. Last, e
cian and the diagonal mass matrix of M. We note P ∈ R|V|×|V| the
f
the cost of the incident edges are updated, reordering the queue to fine-to-coarse restriction matrix and Ni (v) the i-ring of vertex v. We
respect the priority measure. This atomic mesh reduction step is also use the following weighted norm:
f removing
iterated until reaching the desired output resolution for M,
one vertex at a time. kXk2Me = tr(X ⊤ MX)
e (1)
Algorithm 1: Edge-collapse progressive simplification
Input: mesh M = (V, F ), target size N, metric We formulate the preservation of the eigenvectors of the Lapla-
m : V × V 7→ R cian by the commutativity of the Laplacian and the reduction. More
f = (V,
Output: simplified mesh M e Fe) generally, for a signal f ∈ R|V| , we aim for M
e −1 LP
e f = PM −1 L f .
Thus, given K signals to preserve (here, eigenvectors of the Lapla-
e ← V ; Fe ← F ; queue ← {} ;
V
cian), represented as a matrix F ∈ R|V|×K , the reduction metric is:
for edge e ∈ M do
add (e, m(e)) to queue ;
−1 e −1 e 2
e > N and queue not empty do
while |V| | {zLF} −M LPFkMe
E = kP M (2)
(e, c) ← pop edge e with lowest cost c from queue ; Z
e and F)
collapse e (this changes V e ; F and Z are computed only once, at the very beginning. This metric
for n ∈ e’s neighbors do will also preserve eigenvalues, as shown in Appendix C of [LJO19].
update n in queue ;
M −1 L
Optionally, we can augment this reduction algorithm to generate f −−−−−−−→ •
a restriction matrix P to be used for the computation of functional  
maps for instance. More precisely, when collapsing the edge (u, v), Py yP
we generate the restriction matrix Q such that V a f ter = QV be f ore . • −−−−−−−→ fe
Note that all its coefficients are positive, and its rows sum to 1. Me −1 L
e
Then, with Qi the restriction matrix of the i-th operation, the global
restriction matrix is formed as follow: P = Qn Qn−1 ...Q2 Q1 . e is a diagonal matrix, the weighted norm of Equation (1)
Since M

c 2020 The Author(s)

Computer Graphics Forum c 2020 The Eurographics Association and John Wiley & Sons Ltd.
 T. Lescoat, H.-T. Liu, J.-M. Thiery, A. Jacobson, T. Boubekeur, M. Ovsjanikov / Spectral Mesh Simplification

can be decomposed as follow: α=1

α=0 t
cos
kXk2Me = tr(X ⊤ MX) ev k rowv (X)k2
e ⊤ diag(XX ⊤ ) = ∑ M
e = diag(M) α*
v
u v
Thus we can define the metric E in Equation (2) for each vertex v: optimal position

e −1 LPF)k
ev k rowv (PZ − M e 2
E = ∑ Ev and Ev = M .
v

Figure 5: Our vertex optimization scheme finds α∗ ∈ [0, 1] which

We can observe that only the minimizes the cost to determine the merged position (green).
1-ring of v matters, as the cost
will not change for vertices fur- e
ther away from this region (Fig- u v 4. Evaluation
ure 4). More precisely, noting
H = {u, v} ∪ N1 (u, v) when col- We evaluated the performance of our simplification method using a
lapsing edge e = (u, v), the met- variety of criteria. We implemented our technique in C++ with the
ric only changes for vertices of Figure 4: When collapsing help of Spectra, and tested it on a workstation with an Intel Xeon
H. Therefore, we only need the 2- e (blue), only the 1-ring en- 3.0 GHz CPU, 32 GB of RAM. We used the same dataset as Liu et
ring of {u, v} to compute the cost tries of Ew (black) change. al. [LJO19]. In all figures, functional maps are 100 × 100, and we
of a given edge collapse: aim to preserve 100 eigenvectors (K = 100).

cost(e) = E a f ter − E be f ore

4.1. Functional maps
a f ter a f ter be f ore be f ore
= ∑ Ew + ∑ Ew − ∑ Ew − ∑ Ew First, several quantities used in this evaluation are from the func-
w∈H w∈H
/ w∈H w∈H
/
tional maps field; let us briefly introduce the concept of functional
a f ter be f ore
= ∑ Ew − ∑ Ew maps [OBCS∗ 12] here. These are maps between two shapes, but
w∈H w∈H instead of having a point-to-point map, we use a linear mapping
As a result, we only need to compute Ew for w ∈ H, allowing between function spaces. Given a base of functions on each shape
e j , we can map one base function (say, Φi ) on the base functions
Φi , Φ
to track the global cost via local updates. Note that the signals to
preserve, F, are very important when determining which frequencies e j ). This allows to map any function
of the other shape (here, ∑ j Ci j Φ
to keep during the decimation process. Consequently, we use the that we can decompose on these bases from one shape to the other.
K first eigenvectors of L, to focus only on the low-frequencies of The functional map can be represented by the matrix C = (Ci j ).
interest. Finally, once the edge is collapsed, with Q the restriction
Here we consider the functional map C ∈ RK×K between the
matrix associated to this operation, we update P, F, and Z as follow:
input and output shapes, so we don’t take high frequencies into
P ← QP, F ← QF, Z ← QZ.
account. Noting Φ the matrix whose columns are the K first eigen-
e this matrix can be defined as:
vectors of L (and idem for Φ),
e ⊤ MPΦ
C=Φ e
3.4. Merged vertex optimization
and given a function f on M ( f = Φx), its corresponding function
When collapsing an edge, we need to reposition the resulting merged f is (g = ΦCx).
e
e and L,
e in a non- on M To enable fair comparison with other methods,
vertex. This position impacts our metric via P, M e j k e = 1.
we normalize all eigenvectors: ∀i, kΦi kM = 1 and ∀ j, kΦ
linear fashion. Thus, we use an approximation of the metric to find M
e = 1.
We also scale the meshes to have a unit area: tr(M) = tr(M)
its minimizer. Several strategies are possible: (i) always merge at
the edge center, (ii) use a 1D quadratic approximation, restricted on
the edge, or (iii) use a 3D quadratic approximation, unrestricted. Ex- 4.2. Norms
perimentally, we found that (ii) provides the best trade-off between
Ideally, the functional map C between input and output should be as
accuracy and computational cost (see Appendix A for more details).
close as possible to the identity. In order not to rely only on visual
More precisely, let e = (u, v) be the edge to collapse, the resulting inspection, we use two norms on the functional maps to quantify the
position is denoted w(α ∈ [0, 1]) = (1 − α)u + αv, with cost(e, α) result of the simplification:
the collapse cost. We construct the quadratic polynomial p such that e 2
kCΛ − ΛCk
p(α) = cost(e, α) for α ∈ {0, 0.5, 1}, and optimize α∗ ∈ [0, 1] to Laplacian commutativity: kCk2L = (3)
kCk2
minimize p, yielding the optimal merged position (Figure 5).
Orthonormality: kCk2D = kC⊤C − Idk2 (4)
The restriction matrix associated with the collapse is Q ∈ Rn−1×n
with n the number of vertices prior to the collapse. Let ŵ be the where Λ and Λ e are the diagonal matrices of the eigenvalues of the
index of vertex w post collapse, and v the removed vertex: the only Laplacian operator on the input and output mesh, respectively. While
non-zeros are Qûu = 1 − α, Qûv = α, and Qŵw = 1. If needed by the the ideal functional map is often the identity, multiplicity of eigen-
application, we can force Q to be binary by rounding 1 − α and α. values may occur (e.g., in spheres or the bumpy cube in Figure 6).

c 2020 The Author(s)

Computer Graphics Forum c 2020 The Eurographics Association and John Wiley & Sons Ltd.
 T. Lescoat, H.-T. Liu, J.-M. Thiery, A. Jacobson, T. Boubekeur, M. Ovsjanikov / Spectral Mesh Simplification

Both norms are valid in these cases. While Equations (2) & (4) are Nasikun et al. Liu et al. Garland & ~
Ours
not related, we have that if kCk 6= 0 then E = 0 ⇐⇒ kCkL = 0 2018 2019 Heckbert 1997 | |: 622
(proof in Appendix B).
Our goal is to show that C is orthonormal and commutes with
the Laplacians in the reduced basis if and only if it preserves corre-
sponding eigenfunctions and eigenvalues exactly. We do not assume
any constraints on C (e.g., association with a pointwise map).

Theorem 1. For a square functional map, the following statements

are equivalent:

(1) C⊤C = Id and CΛ = ΛC e

f and solves
e is orthonormal on M
(2) the set of functions Y = ΦC
the eigenvalue problem LYe = MYe Λe and moreover Λ = Λ
e
|| . ||L = 4.59 × 103 || . ||L = 5.58 × 103 || . ||L = 92.10 × 103 || . ||L = 23.99 × 103
(3) the set of functions X = ΦC⊤ is orthonormal on M and satis- || . ||D = 7.90 × 100 || . ||D = 1.74 × 100 || . ||D = 18.36 × 100 || . ||D = 10.31 × 100
fies LX = MXΛ and moreover Λ = Λ e
eigenvalues
Intuitively condition (2) (respectively (3)) above implies that Nasikun et al. 2018
C (respectively C⊤ ) preserves the given set of eigenpairs of the Liu et al. 2019
Laplacian. Then, the theorem can also be stated simply as follows: Ours
C is both orthonormal and commutes with the diagonal matrices of Reference
eigenvalues, if and only if it preserves the eigenfunctions and their
corresponding eigenvalues of the Laplacians. We provide a proof of Garland & Heckbert 1997
this theorem as Appendix C. index

Figure 7: To spectrally preserve the Laplacian, both coarse eigen-

4.3. Analysis vectors and eigenvalues should be close to their fine equivalent.
Type of operator. The kind of Laplacian operator on the output
shape differs between the compared methods. The operator retrieved
by the method of Nasikun et al. [NBH18] uses geodesic distances
to determine the sparsity of the Laplacian. This could yield a denser
Norms. As shown before, the ideal case is when both the Laplacian
Laplacian compared to the cotangent Laplacian on a triangle mesh.
commutativity and the orthogonality are zero, as it means the spec-
The method of Liu et al. [LJO19] use the 3-ring of vertices for the
trum is exactly preserved. We observed that in general, the method
Laplacian. Both methods specifically optimize for the operator.
of Nasikun et al. [NBH18] and of Liu et al. [LJO19] outperform our
On the contrary, our method and the method of Garland & Heck- method for small ratios of output size on input size, and the method
bert [GH97] output a mesh from which we can compute the Lapla- of Nasikun et al. is still at least on par when increasing the ratio. The
cian operator using standard formulations on the 1-ring; we use operator coarsening of Liu et al., however, becomes less accurate
the usual cotangent Laplacian for the evaluations. However, this as the ratio increases, notably because more vertices means more
formulation is more constrained, and while usually the Laplacian coefficients in the operator to optimize, making the optimization
operator is easily derived from the input mesh, our problem is an in- more difficult. The method of Garland & Heckbert [GH97] is usu-
verse problem: finding the best reduced mesh that respects a specific ally worse than our method for both metrics, and often create slivers
Laplacian operator (given by its eigenpairs). in the output shape that will heavily impact the Laplacian operator.
We show typical examples of the preservation of eigenvectors in
Figure 3, of eigenvalues in Figure 7, and of the norms in Figure 8.
~
| |: 44,954 | |: 900
Storage. The method of Nasikun et al. [NBH18] aims at generating
coarse eigenpairs, which are dense matrices and thus are very costly
to store. Instead, as with the method of Liu et al. [LJO19], one could
store only the resulting sparse operator, along with the selected
vertices. This is still a lot larger than simply storing a coarse mesh.
We measured the storage size as a function of the target size, in
|| . ||L = 1.54 × 103 bytes/vertex (B/v), yielding a median of 228 B/v for the method
|| . ||D = 1.18 × 100 of Nasikun et al. [NBH18], 262 B/v for the method of Liu et al.
[LJO19], and 48 B/v for our method (more details in Appendix D).
Figure 6: The functional map from the reduction should be block- We require the same storage as the method of Garland & Heckbert
diagonal following the multiplicity of the eigenvalues. [GH97], but with a higher quality Laplacian.

c 2020 The Author(s)

Computer Graphics Forum c 2020 The Eurographics Association and John Wiley & Sons Ltd.
 T. Lescoat, H.-T. Liu, J.-M. Thiery, A. Jacobson, T. Boubekeur, M. Ovsjanikov / Spectral Mesh Simplification

|| . ||L seconds
Garland & Heckbert 1997 Liu et al. 2019
40000
40
Liu et al. 2019
20000

Ours
Nasikun et al. 2018 20
Ours
0
~
0% 20% 40% 60% | |/| | Nasikun et al. 2018
|| . ||D Garland & Heckbert 1997
20 0
~
Garland & Heckbert 1997 0% 20% 40% 60% | |/| |

Figure 10: Typical reduction time (seconds) in function of the ratio

10 output size / input size. Here, |V| = 20685, the method of Liu et
al. [LJO19] was run using MATLAB, and the method of Nasikun et
Nasikun et al. 2018 al. [NBH18] ran out of memory past 21% of input size.
Liu et al. 2019
Ours
0 ~
0% 20% 40% 60% | |/| |

Figure 8: One can see from these plots, k · kL and k · kD over the
output size relative to the input size, that most methods have the sim-
ilar behaviors, except for the method of Liu et al. [LJO19] for which
the optimization becomes harder with more output vertices. The
simplification of Garland & Heckbert [GH97] distort the spectrum
and thus exhibit higher norms.

Liu et al. 2019 Nasikun et al. 2018

mean || . ||L = 2.72 × 103 4.65 × 103
variance || . ||L = 2.99 × 105 4.76 × 105
mean || . ||D = 3.69 × 10-1 3.30 × 100
variance || . ||D = 1.32 × 10-2 1.72 × 10-2

Figure 9: Both [LJO19] and [NBH18] are not deterministic meth-

ods. Although the average of the commutativity and the orthogonal-
ity out of ten runs are small, they still have high variance.

Figure 11: By looking at the norms in function of number of eigen-

Determinism. Previous methods from Nasikun et al. [NBH18] and
vectors relative to the output size, we can see that the output size
Liu et al. [LJO19] have a vertex selection step, following some regu-
should be 3 times the number of eigenvectors for a correct spectral
larity metric. For both of these methods, this step is initialized with
preservation of the Laplacian.
a random selection, making these methods non-deterministic (Fig-
ure 9). This alters not only the final result, but also the time needed
to compute the output. Mirroring Garland & Heckbert [GH97], our
method is deterministic as it only depends on the input. can be quite fast for extreme reductions but timings quickly increase
following the number of vertices. Similarly, the method of Nasikun
Timings. Similarly to the method of Liu et al. [LJO19], the initial et al. [NBH18] can be really fast for small output sizes, notably due
eigenvectors computation depends heavily on the input size and can to GPU usage, and is slower when reducing less, while consuming
be accelerated via a faster eigen solver. Then the reduction time is much more memory. Although Figure 10 show the timings for the
linear in the number of removed vertices, as can be seen in Figure 10. deer head, these behaviors are typical across all of our test dataset,
This behavior is similar for the method of Garland & Heckbert as can be seen in Figure 20 in the appendix.
[GH97], although the latter is much faster as not only the metric is
defined per vertex instead of per edge, but also their setup step is less Number of eigenvectors. Increasing the number of eigenvectors
intensive. As the operator coarsening method [LJO19] optimizes in the metric of Equation (2) leads to a better preservation of high
for a Laplacian operator whose size depends on the output size, it frequencies, up to a certain point. We observed from our test dataset

c 2020 The Author(s)

Computer Graphics Forum c 2020 The Eurographics Association and John Wiley & Sons Ltd.
 T. Lescoat, H.-T. Liu, J.-M. Thiery, A. Jacobson, T. Boubekeur, M. Ovsjanikov / Spectral Mesh Simplification

Garland & eigenvalue, these distances are detailed in the following:

Jakob et. al. 2015 Ours
Heckbert 1997
ddiffusion (u, v,t) = ∑(φi (u) − φi (v))2 e−2λi t
i
dbiharmonic (u, v) = ∑(φi (u) − φi (v))2 /λ2i
i
Z tmax
WKS(u,t) − WKS(v,t)
dWKS (u, v) = dt
tmin WKS(u,t) + WKS(v,t)
dcommute (u, v) = ∑(φi (u) − φi (v))2 /λi
i
where WKS is the wave kernel signature [ASC11]; we also evaluated
the heat kernel signature (HKS) [SOG09]:
|| . ||L = 3.49 × 103 || . ||L = 3.20 × 103 || . ||L = 3.07 × 103 (t−log λi )2 (t−log λi )2
WKS(v,t) = ∑ φ2i (v)e
− −
|| . ||D = 2.01 × 100 || . ||D = 1.88 × 100 || . ||D = 1.23 × 100 2σ2 /∑e 2σ2

i i
Figure 12: Both the shape (by removing small details) and the heat_kernel(u, v,t) = ∑ φi (u)φi (v)e −λi t
discretization (by avoiding slivers) impact the preservation of the i
spectrum (here we remove 90% of the vertices).
HKS(v,t) = heat_kernel(v, v,t) = ∑ φ2i (v)e−λi t
i
While heavily reducing the input mesh will decrease the quality of
that the output size should be at least 3x the number of eigenvectors these distances and signatures, our method preserves them better
to preserve (Figure 11), for a correct spectral preservation of the than the method of Garland & Heckbert [GH97]. In particular, dis-
Laplacian. From Figure 11, we also observe that Laplacian commu- tances on meshes reduced using their method exhibit a lot more
tativity has less variability than orthogonality for high ratios. spurious local optimums than on meshes reduced with our method
(see Figure 13). The Heat Kernel Signature tends to be more resis-
Factors of spectral properties. Both the shape and its discretiza- tant at small values of t, but is less conserved at large values of t
tion impact the spectral properties, especially when the number of (see Figure 14).
vertices is limited. Ideally, faces should be regular, to get a high
quality Laplacian operator [WMKG07]. As the method of Garland &
Heckbert [GH97] focuses on the shape, we also evaluated a remesh- Garland & Heckbert
Reference Ours (-96%) 1997 (-96%)
ing method that focuses on the discretization: Instant Field-aligned
Meshes (IFM) [JTPSH15]. On average, for the same mesh and
diffusion

target size, functional maps from IFM have better orthogonality

while those from our method show better laplacian commutativ-
ity. However, IFM often produced non-manifold meshes, especially
in presence of thin features. As with the method of Garland & | |: 25,214
Heckbert [GH97], we have the guarantee to stay manifold. When
wave kernel biharnmonic

remeshing, and for a given vertex budget, it seems the spectrum

is better preserved when focusing on the discretization [JTPSH15]
than on the shape [GH97] (see Figure 12).
| |: 10,044

5. Applications
Our simplification method was designed to enable faster computa-
tionally expensive shape analysis tasks, by replacing dense input | |: 23,570
meshes with coarser substitutes, yet optimized to carry on as much
as possible the original spectral properties onto which these task
commute

build upon. We illustrate this behavior for two applications: spectral

distance computation and functional map generation.

| |: 25,727

5.1. Spectral distances Figure 13: Spectral distance comparison: the source point is de-
We now show the preservation of spectral distances between vertices, picted in blue and the iso-lines in black, with t set to 0.01 for the
and evaluated several different distances. Noting φi the i-th eigen- diffusion distance. With 25x less vertices in these reduced meshes,
vector of the Laplacian operator on the mesh and λi its associated computing spectral distances is on average 18x faster.

c 2020 The Author(s)

Computer Graphics Forum c 2020 The Eurographics Association and John Wiley & Sons Ltd.
 T. Lescoat, H.-T. Liu, J.-M. Thiery, A. Jacobson, T. Boubekeur, M. Ovsjanikov / Spectral Mesh Simplification

Garland & Heckbert correspondences (%) correspondences (%)

Reference | |: 38,485 Ours (-96%)
1997 (-96%) 100 100 Ours
Ours
HKS (t = 0.1)

50 50
Garland & Garland &
Heckbert 1997 Heckbert 1997

TOSCA isometric TOSCA nonisometric

(90 pairs) (90 pairs)
0 0
0 0.4 0.8 1.2 0 0.5 1.0 1.5
HKS (t = 10000.0)

geodesic error geodesic error

Figure 16: Reducing meshes from the TOSCA dataset from around
30K vertices to 600 vertices allowed to use BCICP [RPWO18] as it
would run out of memory on the fine meshes, and our method pro-
vides better correspondences than the simplification of Garland &
Heckbert [GH97]. On the reduced meshes, the BCICP computation
took on average 67s, from 16s to 142s (variance: 802).
Figure 14: Despite the high reduction, the heat kernel signature is
still similar between coarse and fine mesh when using our method.
Figure 15). This simplification is unaware of the spectrum to pre-
serve and can distort it, limiting the accuracy of the match. We show
5.2. Faster functional maps that we can use our method to enable robust matching, while still
As briefly shown before, functional maps are a powerful tool for find- being faster than without the reduction. This hierarchical scheme
ing correspondences between shapes, and are complemented with can be written in matrix form:
approaches to retrieve a point-to-point mapping. One can perform CX ,Ye = CY,Ye CX ,Y = CXe ,Ye CX ,Xe
shape matching using Product Manifold Filter (PMF) [VLB∗ 17] or
Bijective and Continuous ICP (BCICP) [RPWO18], with excellent where CX ,Y is the functional map from shape X to shape Y. We
results. Both methods are iterative: PMF solves a linear assignment retrieve CX ,Xe and CY,Ye from our mesh simplification method, and
problem at each iteration to determine a bijection between shapes, CXe ,Ye from either PMF [VLB∗ 17] or BCICP [RPWO18]. Then,
which is has a high algorithmic complexity and also needs the shapes CX ,Y can be computed by solving a least squares system. While
to have the same number of vertices. BCICP instead does not need the method of Liu et al. [LJO19] is also suitable for PMF but not
the shapes to have an equal number of vertices, and refines both for BCICP since the latter require a connectivity, we can use our
the functional map and the point-to-point maps at each iteration. method with both PMF and BCICP. We evaluated using BCICP
However these methods do not scale performance-wise when the in the hierarchical functional maps, and compare the results to a
number of vertices increases. reduction with the method of Garland & Heckbert [GH97]. As can
This performance problem can be circumvented by simplifying be observed in Figure 16, our method enable better matching, and is
the meshes prior to the matching, usually using the method of Gar- significantly faster than running BCICP on the fine shape. Indeed,
land & Heckbert [LRR∗ 17], yielding a hierarchical scheme (see while running BCICP on meshes with 600 vertices to 1 minute per
shape pair, BCICP on meshes with 1000 vertices took on average 2
minutes, and ran out of memory for mesh around 30K vertices. For
this application, we used binary restriction matrices.

6. Discussion
Edge flips. The cost of an edge col-
lapse and the preservation metric are
defined in such a way that we can
easily extend it to edge flips, keeping
only edge flips with a negative cost
in order to avoid cycles (Figure 17).
We tried this during the reduction pro- Figure 17: An edge flip
cess, but this always overfitted and will not change the Ew
Figure 15: One can accelerate the computation of functional maps generated poor results (Figure 18). corresponding to the
between detailed meshes by performing shape matching on simpli- Doing a post-process optimization us- 1-ring vertices (black)
fied shapes before upscaling the resulting functional maps. ing edge flips does indeed lead to a and beyond.

c 2020 The Author(s)

Computer Graphics Forum c 2020 The Eurographics Association and John Wiley & Sons Ltd.
 T. Lescoat, H.-T. Liu, J.-M. Thiery, A. Jacobson, T. Boubekeur, M. Ovsjanikov / Spectral Mesh Simplification

result of better quality, but this is quantifiably marginal and is often

not worth the time and complexity.

Figure 19: For large meshes, most of the time is spent in the eigen-
e = 10%|V|). After the setup step, the reduction
solver (left, with |V|
Figure 18: Allowing edge flips often results in whole parts missing. time (right) is linear in the number of collapsed edges.

Limitations. There are at least two main areas where our approach References
can further be developed, both related to time performance. First, [ASC11] AUBRY M., S CHLICKEWEI U., C REMERS D.: The wave kernel
the eigenvectors need to be computed at the beginning of the algo- signature: A quantum mechanical approach to shape analysis. In Proc.
rithm, while one may seek computing them at query time, when ICCV (2011). 7
needed. Second, our cost evaluation involves fairly large matrices, at [CBO∗ 19] C HEN J., B UDNINSKIY M., OWHADI H., BAO H., H UANG
each evaluation which could be optimized using an approximation J., D ESBRUN M.: Material-adapted refinable basis functions for elasticity
simulation. ACM Trans. on Graphics (2019). 2
scheme.
[CMO97] C OHEN J., M ANOCHA D., O LANO M.: Simplifying polygonal
models using successive mappings. In Proc. IEEE Vis (1997). 2
Scalability. By far, the most limiting factor when scaling to large
[CS18] C LAICI S., S OLOMON J.: Wasserstein coresets for Lipschitz costs.
meshes is the eigen solver, which can be seamlessly replaced by arXiv preprint arXiv:1805.07412 (2018). 2
a faster one (see Figure 19). Interestingly, one could use a fast
[CSAD04] C OHEN -S TEINER D., A LLIEZ P., D ESBRUN M.: Variational
approximation [NBH18] for this step. Lowering the number of shape approximation. ACM Trans. on Graphics (2004). 2
eigenvectors K would help for both the setup and the reduction, as
[CVM∗ 96] C OHEN J., VARSHNEY A., M ANOCHA D., T URK G., W E -
the matrices to manipulate would be smaller. It is also possible to BER H., AGARWAL P., B ROOKS F., W RIGHT W.: Simplification en-
partially reduce the input using the QEM [GH97] (e.g., down to velopes. In Proc. SIGGRAPH (1996). 2
100K or 200K vertices) before using our method: the impact of such [GH97] G ARLAND M., H ECKBERT P. S.: Surface simplification using
pre-processing on the spectrum would be limited. quadric error metrics. In Proc. SIGGRAPH (1997). 1, 2, 3, 5, 6, 7, 8, 9
[GH98] G ARLAND M., H ECKBERT P. S.: Simplifying surfaces with color
and texture using quadric error metrics. In Proc. IEEE Vis (1998). 2
7. Conclusion [HCB16] H UGGINS J., C AMPBELL T., B RODERICK T.: Coresets for
scalable Bayesian logistic regression. In Proc. NeurIPS (2016). 2
We have introduced a new mesh simplification algorithm which is
[HDD∗ 93] H OPPE H., D E ROSE T., D UCHAMP T., M C D ONALD J.,
designed to preserve, as much possible, the spectral properties of S TUETZLE W.: Mesh optimization. In Proc. SIGGRAPH (1993). 2
the input surface. Our method is built on a standard graph reduction
[Hop99] H OPPE H.: New quadric metric for simplifying meshes with
algorithm for which we introduced a custom metric driving a cost appearance attributes. In Proc. IEEE Vis (1999). 2
designed specifically to preserve the spectrum, together with a repo-
[JTPSH15] JAKOB W., TARINI M., PANOZZO D., S ORKINE -H ORNUNG
sitioning strategy for merged vertices. We illustrated the superior O.: Instant field-aligned meshes. ACM Trans. on Graphics (2015). 7
behavior of our decimation scheme compared to appearance preserv-
[KS16] K YNG R., S ACHDEVA S.: Approximate Gaussian elimination for
ing methods, for spectral distance computation and functional map Laplacians-fast, sparse, and simple. In Proc. FOCS (2016). 2
generation. Yet, we believe more spectrum-dependent applications
[LFJG17] L IU S., F ERGUSON Z., JACOBSON A., G INGOLD Y. I.: Seam-
may find immediate benefit from our approach. less: seam erasure and seam-aware decoupling of shape from mesh reso-
lution. ACM Trans. on Graphics (2017). 2
[LFZ15] L I D., F EI Y., Z HENG C.: Interactive acoustic transfer approxi-
Acknowledgments mation for modal sound. ACM Trans. on Graphics (2015). 2
[LJO19] L IU H.-T. D., JACOBSON A., OVSJANIKOV M.: Spectral coars-
Parts of this work were supported by the KAUST OSR Award ening of geometric operators. ACM Trans. on Graphics (2019). 2, 3, 4, 5,
No. CRG-2017-3426, the ERC Starting Grant No. 758800 6, 8, 10
(EXPROTEA), NSERC Discovery (RGPIN2017–05235, RG- [Lou19] L OUKAS A.: Graph reduction with spectral and cut guarantees.
PAS–2017–507938), the Ontario Early Research Award program, J. Machine Learning Research (2019). 2
the Canada Research Chairs Program, the Fields Centre for Quantita- [LRR∗ 17] L ITANY O., R EMEZ T., RODOLÀ E., B RONSTEIN A. M.,
tive Analysis and Modelling and gifts by Adobe Systems, Autodesk, B RONSTEIN M. M.: Deep functional maps: Structured prediction for
and MESH Inc. dense shape correspondence. In Proc. ICCV (2017). 8

c 2020 The Author(s)

Computer Graphics Forum c 2020 The Eurographics Association and John Wiley & Sons Ltd.
 T. Lescoat, H.-T. Liu, J.-M. Thiery, A. Jacobson, T. Boubekeur, M. Ovsjanikov / Spectral Mesh Simplification

[LV18] L OUKAS A., VANDERGHEYNST P.: Spectrally approximating

large graphs with smaller graphs. In Proc. ICML (2018). 2
[NBH18] NASIKUN A., B RANDT C., H ILDEBRANDT K.: Fast approxi-
mation of Laplace-Beltrami eigenproblems. Comp. Graph. Forum (2018).
2, 5, 6, 9
[OBCS∗ 12] OVSJANIKOV M., B EN -C HEN M., S OLOMON J.,
B UTSCHER A., G UIBAS L.: Functional maps: A flexible repre-
sentation of maps between shapes. ACM Trans. on Graphics (2012). 1, 2,
4
[OS19] OWHADI H., S COVEL C.: Operator-Adapted Wavelets, Fast
Solvers, and Numerical Homogenization: From a Game Theoretic Ap-
proach to Numerical Approximation and Algorithm Design. 2019. 2
[PH97] P OPOVI Ć J., H OPPE H.: Progressive simplicial complexes. In e with |V| fixed per curve.
Figure 20: Timings as a function of |V|
Proc. SIGGRAPH (1997). 2 The behavior shown in Figure 10 is also visible here.
[RPWO18] R EN J., P OULENARD A., W ONKA P., OVSJANIKOV M.:
Continuous and orientation-preserving correspondences via functional e ⊤ M,
side is equal to PΦΛ. We left-multiply both sides by Φ e yielding
maps. ACM Trans. on Graphics (2018). 8 CΛ = Φ ⊤
e LPΦ
e since C = Φ ⊤
e MPΦ.
e eΦ
Recalling that L e =MeΦ e the
e Λ,
[SGG∗ 00] S ANDER P. V., G U X., G ORTLER S. J., H OPPE H., S NYDER right-hand side is equal to (M e ⊤ PΦ, thus ΛC.
e Λ)
eΦ e This means that
J.: Silhouette clipping. In Proc. SIGGRAPH (2000). 2 e proving that E = 0 ⇐⇒ kCkL = 0.
CΛ = ΛC,
[SOG09] S UN J., OVSJANIKOV M., G UIBAS L.: A concise and provably
informative multi-scale signature based on heat diffusion. In Proc. SGP
(2009). 7
C. Proof of Theorem 1
[SVJ15] S ACHT L., VOUGA E., JACOBSON A.: Nested cages. ACM
Trans. on Graphics (2015). 2 Proof. We will prove the equivalence between (1) and (2). The
[SZL∗ 92] S CHROEDER W. J., Z ARGE J. A., L ORENSEN W. E., ET AL .: equivalence between (1) and (3) is proved identically. Suppose (2)
Decimation of triangle meshes. In Proc. SIGGRAPH (1992). 2 holds. To show that (1) must hold, first note that from orthonormality
of Y on Mf by definition we get Y ⊤ MY e = Id, i.e. C⊤ Φe ⊤M e = Id,
e ΦC
[VLB∗ 17] V ESTNER M., L ÄHNER Z., B OYARSKI A., L ITANY O.,
e f ⊤
S LOSSBERG R., R EMEZ T., RODOLA E., B RONSTEIN A., B RONSTEIN and since Φ is orthonormal on M this implies C C = Id. Moreover,
M., K IMMEL R., C REMERS D.: Efficient deformable shape correspon- by assumption LY e = MY e Λ,e and Λ = Λ.e Thus, L e =M
eΦC e
e ΦCΛ. Since
dence via kernel matching. 8 by definition LeΦe =M eΦ eΛ e we get M eΦ e =M
e ΛC e
e ΦCΛ which implies
[WMKG07] WARDETZKY M., M ATHUR S., K AELBERER F., G RINSPUN e = CΛ.
ΛC
E.: Discrete Laplace operators: No free lunch. In Proc. SGP (2007). 7
Conversely, suppose that (1) holds. If Y = ΦC, e then C⊤C =
[ZG02] Z ELINKA S., G ARLAND M.: Permission grids: Practical, error- ⊤ e ⊤ e⊤ e e e is orthonormal
bounded simplification. ACM Trans. on Graphics (2002). 2 Id implies: Y MY = C Φ M ΦC = Id, since Φ
with respect to M.e Now, LY e =L e =M
eΦC eΦ e Since ΛC
e ΛC. e = CΛ by
assumption, we get LY e =M e
e ΦCΛ e Λ. Therefore, Y solves the
= MY
Appendix eigenvalue problem of (L, e M)
e with the eigenvalues Λ. It remains
e
to prove that Λ = Λ. For this, note that ΛC e = CΛ implies Ci2j (eλi −
A. Merged position
2 th
λ j ) = 0 ∀i, j, where eλi is the i eigenvalue of L
e (idem for L). Using
Strategy k · kL k · kD Time this and C⊤C = Id implies that C must be block orthonormal with
(i) middle 1.0 1.0 1.0 blocks corresponding to the equal eigenvalues (to see this, note that
(ii) on edge 0.7 0.6 2.1 for each eλi there must be an equal λ j otherwise a row or column of
(iii) unrestricted 0.8 0.9 4.9 C would have to be zero). Moreover the blocks must be square by
e
orthonormality of C, so that Λ = Λ.
Table 1: Median norms and time relative to strategy (i), for the same
mesh and parameters.
D. Storage sizes
Let cost(e, x) be the cost of collapsing e with x the merged posi- Let n be the target size and c the number of coefficient per row for
tion. As α is needed for P (and thus impacts the cost), we define it e We only consider 64bit floats here. We can store the selected
L.
(x−u)·(v−u)
as kv−uk2 and clamp it in [0, 1]. For strategy (iii), we sample e which is diagonal,
subset of vertices via a list of n 32bit indices. M,
e is symmetric so we need to store only n(c + 1)/2
requires n floats. L
cost(e, x) for x in a sphere around e to determine the approximating
e never
coefficients, each one taking one float and two 16bit integers (L
quadric (in R4×4 ), which we then minimize to get x∗ (and α∗ ). We
refer to Table 1 for a comparison relative to strategy (i). exceeds 65,535 rows in our tests). Overall, this gives (18 + 6c)n
bytes. Considering the observed median size for the method of Liu
et al. [LJO19] (about 262 B/v), we find c ≈ 41, close to what the
B. Relation between Equations (2) & (3) authors report just before their Section 3.2.
Proof. Having E = 0 (from Equation (2)) is equivalent to kCkL = 0 For meshes, we need 3n floats for vertices. Assuming an average
e are strictly positive,
(from Equation (3)). Since all coefficients of M valence of 6, there are 2n triangles, needing 3 × 2n 32bit integers.
E = 0 ⇐⇒ PM −1 LΦ = M e −1 LPΦ.
e As LΦ = MΦΛ, the left-hand This yields a total of 48n bytes, which we experimentally observe.

c 2020 The Author(s)

Computer Graphics Forum c 2020 The Eurographics Association and John Wiley & Sons Ltd.