forked from alecjacobson/geometry-processing-smoothing
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathREADME.html
644 lines (519 loc) · 29.4 KB
/
README.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8"/>
<link type="text/css" rel="stylesheet" href="shared/css/style.css"/>
<script type="text/javascript" src="shared/js/MathJax/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
</head>
<body>
<div style="display:none">
$$\newcommand{\mat}[1]{\mathbf{#1}}$$
$$\newcommand{\vec}[1]{\mathbf{#1}}$$
$$\newcommand{\A}{\mat{A}}$$
$$\newcommand{\B}{\mat{B}}$$
$$\newcommand{\C}{\mat{C}}$$
$$\newcommand{\D}{\mat{D}}$$
$$\newcommand{\I}{\mat{I}}$$
$$\newcommand{\G}{\mat{G}}$$
$$\newcommand{\N}{\mat{N}}$$
$$\newcommand{\P}{\mat{P}}$$
$$\newcommand{\Rot}{\mat{R}}$$
$$\newcommand{\R}{\mathbb{R}}$$
$$\newcommand{\One}{\mathbf{1}}$$
$$\newcommand{\S}{\mathcal{S}}$$
$$\newcommand{\M}{\mat{M}}$$
$$\newcommand{\U}{\mat{U}}$$
$$\newcommand{\L}{\mat{L}}$$
$$\newcommand{\V}{\mat{V}}$$
$$\newcommand{\W}{\mat{W}}$$
$$\newcommand{\X}{\mat{X}}$$
$$\newcommand{\Y}{\mat{Y}}$$
$$\newcommand{\c}{\vec{c}}$$
$$\newcommand{\f}{\vec{f}}$$
$$\newcommand{\g}{\vec{g}}$$
$$\newcommand{\n}{\vec{n}}$$
$$\newcommand{\e}{\vec{e}}$$
$$\newcommand{\p}{\vec{p}}$$
$$\newcommand{\tr}[1]{\mathop{\text{tr}}{\left(#1\right)}}$$
$$\newcommand{\t}{\vec{t}}$$
$$\newcommand{\x}{\vec{x}}$$
$$\newcommand{\y}{\vec{y}}$$
$$\newcommand{\u}{\vec{u}}$$
$$\newcommand{\z}{\vec{z}}$$
$$\renewcommand{\v}{\vec{v}}$$
$$\newcommand{\transpose}{{\mathsf T}}$$
$$\newcommand{\argmin}{\mathop{\text{argmin}}}$$
$$\newcommand{\argmax}{\mathop{\text{argmax}}}$$
</div>
<h1 id="geometryprocessing–smoothing">Geometry Processing – Smoothing</h1>
<blockquote>
<p><strong>To get started:</strong> Fork this repository then issue</p>
<pre><code>git clone --recursive http://github.com/[username]/geometry-processing-smoothing.git
</code></pre>
</blockquote>
<h2 id="installationlayoutandcompilation">Installation, Layout, and Compilation</h2>
<p>See
<a href="http://github.com/alecjacobson/geometry-processing-introduction">introduction</a>.</p>
<h2 id="execution">Execution</h2>
<p>Once built, you can execute the assignment from inside the <code>build/</code> by running
on a given mesh with given scalar field (in
<a href="http://libigl.github.io/libigl/file-formats/dmat.html">dmat</a> format).</p>
<pre><code>./smoothing [path to mesh.obj] [path to data.dmat]
</code></pre>
<p>or to load a mesh with phony noisy data use:</p>
<pre><code>./smoothing [path to mesh.obj] n
</code></pre>
<p>or to load a mesh with smooth z-values as data (for mesh smoothing only):</p>
<pre><code>./smoothing [path to mesh.obj]
</code></pre>
<h2 id="background">Background</h2>
<p>In this assignment we will explore how to smooth a data <em>signal</em> defined over a
curved surface. The data <em>signal</em> may be a scalar field defined on a static
surface: for example, noisy temperatures on the surface of an airplane.
Smoothing in this context can be understood as <a href="https://en.wikipedia.org/wiki/Noise_reduction">data
denoising</a>:</p>
<figure>
<img src="images/plane-smooth-data.gif" alt="" />
</figure>
<p>The <em>signal</em> could also be the geometric positions of the surface itself. In
this context, smoothing acts also affects the underlying geometry of the
domain. We can understand this operation as <a href="https://en.wikipedia.org/wiki/Surface_fairng">surface
fairing</a>:</p>
<figure>
<img src="images/plane-smooth-geometry.gif" alt="" />
</figure>
<h3 id="flow-basedformulation">Flow-based formulation</h3>
<p>In both cases, the underlying mathematics of both operations will be very
similar. If we think of the signal as undergoing a <em>flow</em> toward a smooth
solution over some phony notion of “time”, then the governing partial
differential equation we will start with sets the change in signal value <span class="math">\(u\)</span>
over time <a href="https://en.wikipedia.org/wiki/Proportionality_(mathematics)">proportional
to</a> the
<a href="https://en.wikipedia.org/wiki/Laplace_operator">Laplacian</a> of the signal <span class="math">\(∆u\)</span>
(for now, roughly the second derivative of the signal as we move <em>on</em> the
surface):</p>
<p><span class="math">\[
\frac{∂ u}{∂ t} = λ ∆ u,
\]</span></p>
<p>where the scalar parameter <span class="math">\(λ\)</span> controls the rate of smoothing.</p>
<p>When the signal is the surface geometry, we call this a <a href="https://en.wikipedia.org/wiki/Geometric_flow">geometric
flow</a>.</p>
<p>There are various ways to motivate this choice of flow for
data-/geometry-smoothing. Let us consider one way that will introduce the
Laplacian as a form of local averaging.</p>
<p>Given a noisy signal <span class="math">\(f\)</span>, intuitively we can <em>smooth</em> <span class="math">\(f\)</span> by averaging every
value with its neighbors’ values. In continuous math, we might write that the
smoothed value <span class="math">\(u(\x)\)</span> at any point on our surface <span class="math">\(\x ∈ \S\)</span> should be equal to
the average value of some small
<a href="https://en.wikipedia.org/wiki/Ball_(mathematics)">ball</a> of nearby points:</p>
<p><span class="math">\[
u(\x) = \frac{1}{|B(\x))|} ∫_{B(\x)} f(\z) \;d\z,
\]</span></p>
<p>If the ball <span class="math">\(B(\x)\)</span> is small, then we will have to repeat this averaging many
times to see a global smoothing effect. Hence, we can write that the current
value <span class="math">\(u^t\)</span> <em>flows</em> toward smooth solution by small steps <span class="math">\(δt\)</span> in time:</p>
<p><span class="math">\[
u^{t+δt}(\x) = \frac{1}{|B(\x))|} ∫_{B(\x)} u^t(\z) \;d\z.
\]</span></p>
<p>Subtracting the current value <span class="math">\(u^t(\x)\)</span> from both sides and introducing a
flow-speed parameter <span class="math">\(λ\)</span> we have a flow equation
describing the change in value as an integral of relative values:</p>
<p><span class="math">\[
\frac{∂ u}{∂ t}
= λ \frac{1}{|B(\x))|} ∫_{B(\x)} (u(\z)-u(\x)) \;d\z.
\]</span></p>
<p>For harmonic functions, <span class="math">\(∆u =0\)</span>, this integral becomes zero in the limit as the
radius of the ball shrinks to zero via satisfaction of the
<a href="https://en.wikipedia.org/wiki/Harmonic_function#The_mean_value_property">mean value
theorem</a>.
It follows for a non-harmonic <span class="math">\(∆u ≠ 0\)</span> this integral is equal to the Laplacian
of the <span class="math">\(u\)</span>, so we have arrived at our flow equation:</p>
<p><span class="math">\[
\frac{∂ u}{∂ t}
= \lim_{|B(\x)| → 0} λ \frac{1}{|B(\x))|} ∫_{B(\x)} (u(\z)-u(\x)) \;d\z
= λ ∆ u.
\]</span></p>
<h3 id="energy-basedformulation">Energy-based formulation</h3>
<p>Alternatively, we can think of a single smoothing operation as the solution to
an energy minimization problem. If <span class="math">\(f\)</span> is our noisy signal over the surface,
then we want to find a signal <span class="math">\(u\)</span> such that it simultaneously minimizes its
difference with <span class="math">\(f\)</span> and minimizes its variation over the surface:</p>
<p><span class="math">\[
u^*
= \argmin_u E_(u)
=
\argmin_u ½∫_\S ( \underbrace{(f-u)²}_\text{data} +
\underbrace{λ‖∇u‖²}_\text{smoothness} )\;dA,
\]</span></p>
<p>where again the scalar parameter <span class="math">\(λ\)</span> controls the rate of smoothing. This
energy-based formulation is equivalent to the flow-based formulation.
Minimizing these energies is identical to stepping forward one temporal unit in
the flow.</p>
<h4 id="calculusofvariations">Calculus of variations</h4>
<p>In the smooth setting, our energy <span class="math">\(E\)</span> is a function that measures scalar value
of a given function <em>u</em>, making it a
<a href="https://en.wikipedia.org/wiki/Functional_(mathematics)">functional</a>. To
understand how to <em>minimize</em> a functional with respect to an unknown function,
we will need concepts from the <a href="https://en.wikipedia.org/wiki/Calculus_of_variations">calculus of
variations</a>.</p>
<p>We are used to working with minimizing quadratic <em>functions</em> with respect to a
discrete set of variables, where the minimum is obtained when the gradient of
the energy with respect to the variables is zero.</p>
<p>In our case, the functional <span class="math">\(E(u)\)</span> is quadratic in <span class="math">\(u\)</span> (recall that the
<a href="https://en.wikipedia.org/wiki/Gradient">gradient operator</a> <span class="math">\(∇\)</span> is a linear
operator). The function <span class="math">\(u\)</span> that minimizes <span class="math">\(E(u)\)</span> will be obtained when any
small change or <em>variation</em> in <span class="math">\(u\)</span> has no change on the energy values. To
create a small change in a function <span class="math">\(u\)</span> we will add another function <span class="math">\(v\)</span> times
a infinitesimal scalar <span class="math">\(ε\)</span>. If <span class="math">\(E(u)\)</span> is minimized for a function <span class="math">\(w\)</span> and we
are given another arbitrary function <span class="math">\(v\)</span>, then let us define a function new
function </p>
<p><span class="math">\[
Φ(ε) = E(w+εv) = ½∫_\S ((f-w+εv)² + λ ‖∇w + ε∇v‖²) \;dA,
\]</span>
where we observe that <span class="math">\(Φ\)</span> is quadratic in <span class="math">\(ε\)</span>.</p>
<p>Since <span class="math">\(E(w)\)</span> is minimal then <span class="math">\(Φ\)</span> is minimized when <span class="math">\(ε\)</span> is zero, and if <span class="math">\(Φ\)</span> is
minimal at <span class="math">\(ε=0\)</span>, then the derivative of <span class="math">\(Φ\)</span> with respect <span class="math">\(ε\)</span> must be zero:</p>
<p><span class="math">\[
\begin{align}
0 & = \left.\frac{∂Φ}{∂ε} \right|_{ε = 0},\\
& = \left.\frac{∂}{∂ε} ½∫_\S ( (f-w-εv)² + λ ‖∇w + ε∇v‖²)\;dA, \right|_{ε = 0} \\
& = \left.\frac{∂}{∂ε} ½∫_\S (
f^2 - 2wf - 2εfv +w²+2εvw +ε²v² + λ ‖∇w‖² + λ2ε∇v⋅∇w + λ ε²‖∇w‖²) \;dA \right|_{ε = 0}\\
& = \left.∫_\S (-fv + vw +2εvw + λ∇v⋅∇w + λ ε‖∇w‖²) \;dA \right|_{ε = 0}\\
& = ∫_\S (v(w-f) + λ∇v⋅∇w )\;dA.
\end{align}
\]</span></p>
<p>The choice of “test” function <span class="math">\(v\)</span> was arbitrary, so this must hold for any
(reasonable) choice of <span class="math">\(v\)</span>:</p>
<p><span class="math">\[
0 = ∫_\S (v(w-f) + λ∇v⋅∇w) \;dA \quad ∀ v.
\]</span></p>
<p>It is difficult to claim much about <span class="math">\(w\)</span> from this equation directly because
derivatives of <span class="math">\(v\)</span> are still involved. We can <em>move</em> a derivative from <span class="math">\(v\)</span> to a
<span class="math">\(w\)</span> by applying <a href="https://en.wikipedia.org/wiki/Green's_identities">Green’s first
identity</a>:</p>
<p><span class="math">\[
0 = ∫_\S (v(w-f) - λv∆w )\;dA \quad (+ \text{boundary term} )\quad ∀ v,
\]</span>
where we choose to <em>ignore</em> the boundary terms (for now) or equivalently we
agree to work on <em>closed</em> surfaces <span class="math">\(\S\)</span>.</p>
<p>Since this equality must hold of <em>any</em> <span class="math">\(v\)</span> let us consider functions that are
little <a href="https://en.wikipedia.org/wiki/Bump_function">“blips”</a> centered at any
arbitrary point <span class="math">\(\x ∈ \S\)</span>. A function <span class="math">\(v\)</span> that is one at <span class="math">\(\x\)</span> and quickly
decays to zero everywhere else. To satisfy the equation above at <span class="math">\(\x\)</span> with this
blip <span class="math">\(v\)</span> we must have that:</p>
<p><span class="math">\[
w(\x)-f(\x) = λ∆w(\x).
\]</span>
The choice of <span class="math">\(\x\)</span> was arbitrary so this must hold <em>everywhere</em>.</p>
<p>Because we invoke <em>variations</em> to arrive at this equation, we call the
<em>energy-based</em> formulation a <em>variational formulation</em>.</p>
<h3 id="implicitsmoothingiteration">Implicit smoothing iteration</h3>
<p>Now we understand that the flow-based formulation and the variational
formulation lead to the same system, let us concretely write out the implicit
smoothing step. </p>
<p>Letting <span class="math">\(u^0 = f\)</span> we compute a new smoothed function <span class="math">\(u^{t+1}\)</span> given the
current solution <span class="math">\(u^t\)</span> by solving the <em>linear</em> system of equations:</p>
<p><span class="math">\[
u^t(\x) = (\text{id}-λ∆)u^{t+1}(\x), \quad ∀ \x ∈ \S
\]</span>
where <span class="math">\(\text{id}\)</span> is the <a href="https://en.wikipedia.org/wiki/Identity_function">identity
operator</a>. In the discrete
case, we will need discrete approximations of the <span class="math">\(\text{id}\)</span> and <span class="math">\(∆\)</span>
operators.</p>
<h2 id="discretelaplacian">Discrete Laplacian</h2>
<p>There are many ways to derive a discrete
approximation of the Laplacian <span class="math">\(∆\)</span> operator on a triangle mesh using:</p>
<ul>
<li><a href="https://en.wikipedia.org/wiki/Finite_volume_method">finite volume method</a>,
<ul>
<li>“The solution of partial differential equations by means of electrical
networks” [MacNeal 1949, pp. 68],</li>
<li>“Discrete differential-geometry operators for triangulated
2-manifolds” [Meyer et al. 2002],</li>
<li><em>Polygon mesh processing</em> [Botsch et al. 2010],</li>
</ul></li>
<li><a href="https://en.wikipedia.org/wiki/Finite_element_method">finite element
method</a>,
<ul>
<li>“Variational methods for the solution of problems of equilibrium and
vibrations” [Courant 1943],</li>
<li><em>Algorithms and Interfaces for Real-Time Deformation of 2D and 3D Shapes</em>
[Jacobson 2013, pp. 9]</li>
</ul></li>
<li><a href="https://en.wikipedia.org/wiki/Discrete_exterior_calculus">discrete exterior
calculus</a>
<ul>
<li><em>Discrete Exterior Calculus</em> [Hirani 2003, pp. 69]</li>
<li><em>Discrete Differential Geometry: An Applied Introduction</em> [Crane 2013,
pp. 71]</li>
</ul></li>
<li><a href="https://en.wikipedia.org/wiki/Mean_curvature_flow">gradient of surface area → mean curvature
flow</a>
<ul>
<li>“Computing Discrete Minimal Surfaces and Their Conjugates” [Pinkall &
Polthier 1993]</li>
</ul></li>
</ul>
<p>All of these techniques will produce the <em>same</em> sparse <em>Laplacian matrix</em> <span class="math">\(\L ∈ \R^{n×n}\)</span>
for a mesh with <span class="math">\(n\)</span> vertices. </p>
<h3 id="finiteelementderivationofthediscretelaplacian">Finite element derivation of the discrete Laplacian</h3>
<p>We want to approximate the Laplacian of a function <span class="math">\(∆u\)</span>. Let us consider <span class="math">\(u\)</span> to
be <a href="https://en.wikipedia.org/wiki/Piecewise_linear_function">piecewise-linear</a>
represented by scalar values at each vertex, collected in <span class="math">\(\u ∈ \R^n\)</span>.</p>
<p>Any piecewise-linear function can be expressed as a sum of values at mesh
vertices times corresponding piecewise-linear basis functions (a.k.a hat
functions, <span class="math">\(φ_i\)</span>):</p>
<p><span class="math">\[
\begin{align}
u(\x) &= ∑_{i=1}^n u_i φ_i(\x), \\
φ(\x) &= \begin{cases}
1 & \text{if $\x = \v_i$}, \\
\frac{\text{Area($\x$,$\v_j$,$\v_k$)}}{\text{Area($\v_i$,$\v_j$,$\v_k$)}}
& \text{if $\x ∈ \text{triangle}(i,j,k)$}, \\
0 & \text{otherwise}.
\end{cases}
\end{align}
\]</span></p>
<figure>
<img src="images/hat-function.png" alt="" />
</figure>
<p>By plugging this definition into our smoothness energy above, we have discrete
energy that is quadratic in the values at each mesh vertex:</p>
<p><span class="math">\[
\begin{align}
∫_\S ‖∇u(\x)‖² \;dA
&= ∫_\S \left\|∇\left(∑_{i=1}^n u_i φ_i(\x)\right)\right\|^2 \;dA \\
&= ∫_\S \left(∑_{i=1}^n u_i ∇φ_i(\x)\right)⋅\left(∑_{i=1}^n u_i ∇φ_i(\x)\right) \;dA \\
&= ∫_\S ∑_{i=1}^n ∑_{j=1}^n ∇φ_i⋅∇φ_j u_i u_j \;dA \\
&= \u^\transpose \L \u, \quad \text{where } L_{ij} = ∫_\S ∇φ_i⋅∇φ_j \;dA.
\end{align}
\]</span></p>
<p>By defining <span class="math">\(φ_i\)</span> as piecewise-linear hat functions, the values in the system
matrix <span class="math">\(L_{ij}\)</span> are uniquely determined by the geometry of the underlying mesh.
These values are famously known as <em>cotangent weights</em>. “Cotangent”
because, as we will shortly see, of their trigonometric formulae and “weights”
because as a matrix <span class="math">\(\L\)</span> they define a weighted <a href="https://en.wikipedia.org/wiki/Laplacian_matrix">graph
Laplacian</a> for the given mesh.
Graph Laplacians are employed often in geometry processing, and often in
discrete contexts ostensibly disconnected from FEM. The choice or manipulation
of Laplacian weights and subsequent use as a discrete Laplace operator has been
a point of controversy in geometry processing research (see “Discrete laplace
operators: no free lunch” [Wardetzky et al. 2007]).</p>
<p>We first notice that <span class="math">\(∇φ_i\)</span> are constant on each triangle, and only nonzero on
triangles incident on node <span class="math">\(i\)</span>. For such a triangle, <span class="math">\(T_α\)</span>, this <span class="math">\(∇φ_i\)</span> points
perpendicularly from the opposite edge <span class="math">\(e_i\)</span> with inverse magnitude equal to
the height <span class="math">\(h\)</span> of the triangle treating that opposite edge as base:
\begin{equation} |∇φ_i| = \frac{1}{h} = \frac{|\e_i|}{2A}, \end{equation}
where <span class="math">\(\e_i\)</span> is the edge <span class="math">\(e_i\)</span> as a vector and <span class="math">\(A\)</span> is the area of the triangle.</p>
<figure>
<img src="images/hat-function-gradient.png" alt="Leftpoints perpendicular to opposite edges. Right" />
<figcaption> <meta name="Left" content="the gradient $∇ φ_i$ of a hat function $φ_i$ is piecewise-constant and"/>
<meta name="points perpendicular to opposite edges. Right" content="hat function gradients $∇ φ_i$
and $∇ φ_j$ of neighboring nodes meet at angle $θ = π -
α_{ij}$."/>
</figcaption>
</figure>
<p>Now, consider two neighboring nodes <span class="math">\(i\)</span> and <span class="math">\(j\)</span>, connected by some edge
<span class="math">\(\e_{ij}\)</span>. Then <span class="math">\(∇φ_i\)</span> points toward node <span class="math">\(i\)</span> perpendicular to <span class="math">\(\e_i\)</span> and
likewise <span class="math">\(∇ φ_j\)</span> points toward node <span class="math">\(j\)</span> perpendicular to <span class="math">\(\e_j\)</span>. Call the angle
formed between these two vectors <span class="math">\(θ\)</span>. So we may write:</p>
<p><span class="math">\[
∇ φ_i ⋅ ∇ φ_j = \|∇ φ_i\| \|∇ φ_j\| \cos θ =
\frac{\|\e_j\|}{2A}\frac{\|\e_i\|}{2A} \cos θ.
\]</span></p>
<p>Now notice that the angle between <span class="math">\(\e_i\)</span> and <span class="math">\(\e_j\)</span>, call it <span class="math">\(α_{ij}\)</span>,
is <span class="math">\(π - θ\)</span>, but more importantly that:
<span class="math">\[
\cos θ = - \cos \left(π - θ\right) = -\cos α_{ij}.
\]</span>
So, we can rewrite equation the cosine law equation above into:
<span class="math">\[
-\frac{\|\e_j\|}{2A}\frac{\|\e_i\|}{2A} \cos
α_{ij}.
\]</span>
Now, apply the definition of sine for right triangles:
<span class="math">\[
\sin α_{ij} = \frac{h_j}{\|\e_i\|} = \frac{h_i}{\|\e_j\|},
\]</span>
where <span class="math">\(h_i\)</span> is the height of the triangle treating <span class="math">\(\e_i\)</span> as base, and
likewise for <span class="math">\(h_j\)</span>. Rewriting the equation above, replacing one of the edge norms,
e.g.\ <span class="math">\(\|\e_i\|\)</span>:
<span class="math">\[
-\frac{\|\e_j\|}{2A} \frac{\frac{h_j}{\sinα_{ij}}}{2A} \cos α_{ij}.
\]</span></p>
<p>Combine the cosine and sine terms:
<span class="math">\[
-\frac{\|\e_j\|}{2A} \frac{h_j}{2A} \cot α_{ij}.
\]</span></p>
<p>Finally, since <span class="math">\(\|\e_j\|h_j=2A\)</span>, our constant dot product of these
gradients in our triangle is:
<span class="math">\[
∇ φ_i ⋅ ∇ φ_j = -\frac{\cot α_{ij}}{2A}.
\]</span></p>
<p>Similarly, inside the other triangle <span class="math">\(T_β\)</span> incident
on nodes <span class="math">\(i\)</span> and <span class="math">\(j\)</span> with angle <span class="math">\(β_{ij}\)</span> we have a constant dot
product:
<span class="math">\[
∇ φ_i ⋅ ∇ φ_j = -\frac{\cot β_{ij}}{2B},
\]</span>
where <span class="math">\(B\)</span> is the area <span class="math">\(T_β\)</span>.</p>
<p>Recall that <span class="math">\(φ_i\)</span> and <span class="math">\(φ_j\)</span> are only both nonzero inside these two
triangles, <span class="math">\(T_α\)</span> and <span class="math">\(T_β\)</span>. So, since these constants are inside an
integral over area we may write:
<span class="math">\[
\int\limits_\S ∇ φ_i ⋅ ∇ φ_j \;dA =
\left.A∇ φ_i ⋅ ∇ φ_j \right|_{T_α} + \left.B∇ φ_i ⋅ ∇ φ_j \right|_{T_β}
=
-\frac{1}{2} \left( \cot α_{ij} + \cot β_{ij} \right).
\]</span></p>
<h2 id="massmatrix">Mass matrix</h2>
<p>Treated as an <em>operator</em> (i.e., when used multiplied against a vector <span class="math">\(\L\u\)</span>),
the Laplacian matrix <span class="math">\(\L\)</span> computes the local integral of the Laplacian of a
function <span class="math">\(u\)</span>. In the energy-based formulation of the smoothing problem this is
not an issue. If we used a similar FEM derivation for the <em>data term</em> we would
get another sparse matrix <span class="math">\(\M ∈ \R^{n × n}\)</span>:</p>
<p><span class="math">\[
∫_\S (u-f)² \;dA
= ∫_\S ∑_{i=1}^n ∑_{j=1}^n φ_i⋅φ_j (u_i-f_i) (u_j-f_j) \;dA =
(\u-\f)^\transpose \M (\u-\f),
\]</span>
where <span class="math">\(\M\)</span> as an operator computes the local integral of a function’s value
(i.e., <span class="math">\(\M\u\)</span>).</p>
<p>This matrix <span class="math">\(\M\)</span> is often <em>diagonalized</em> or <em>lumped</em> into a diagonal matrix,
even in the context of FEM. So often we will simply set:</p>
<p><span class="math">\[
M_{ij} =
\begin{cases}
⅓ ∑_{t=1}^m \begin{cases}
\text{Area}(t) & \text{if triangle $t$ contains vertex $i$} \\
0 & \text{otherwise}
\end{cases}
& \text{if $i=j$}\\
0 & \text{otherwise},
\end{cases}
\]</span>
for a mesh with <span class="math">\(m\)</span> triangles.</p>
<p>If we start directly with the continuous smoothing iteration equation, then we
have a point-wise equality. To fit in our integrated Laplacian <span class="math">\(\L\)</span> we should
convert it to a point-wise quantity. From a units perspective, we need to
divide by the local area. This would result in a discrete smoothing iteration
equation:</p>
<p><span class="math">\[
\u^t = (\I - λ\M^{-1} \L)\u^{t+1},
\]</span></p>
<p>where <span class="math">\(\I ∈ \R^{n×n}\)</span> is the identity matrix. This equation is <em>correct</em> but
the resulting matrix <span class="math">\(\A := \I - λ\M^{-1} \L\)</span> is not symmetric and thus slower
to solve against.</p>
<p>Instead, we could take the healthier view of requiring our smoothing iteration
equation to hold in a locally integrated sense. In this case, we replace mass
matrices on either side:</p>
<p><span class="math">\[
\M \u^t = (\M - λ\L)\u^{t+1}.
\]</span></p>
<p>Now the system matrix <span class="math">\(\A := \M + λ\L\)</span> will be symmetric and we can use
<a href="https://en.wikipedia.org/wiki/Cholesky_decomposition">Cholesky factorization</a>
to solve with it.</p>
<h3 id="laplaceoperatorisintrinsic">Laplace Operator is Intrinsic</h3>
<p>The discrete Laplacian operator and its accompanying mass matrix are
<em>intrinsic</em> operators in the sense that they <em>only</em> depend on lengths. In
practical terms, this means we do not need to know <em>where</em> vertices are
actually positioned in space (i.e., <span class="math">\(\V\)</span>). Rather we only need to know the
relative distances between neighboring vertices (i.e., edge lengths). We do not
even need to know which dimension this mesh is <a href="https://en.wikipedia.org/wiki/Embedding">living
in</a>.</p>
<p>This also means that applying a transformation to a shape that does not change
any lengths on the surface (e.g., bending a sheet of paper) will have no affect
on the Laplacian.</p>
<h3 id="datadenoising">Data denoising</h3>
<p>For the data denoising application, our geometry of the domain is not changing
only the scalar function living upon it. We can build our discrete Laplacian
<span class="math">\(\L\)</span> and mass matrix <span class="math">\(\M\)</span> and apply the above formula with a chosen <span class="math">\(λ\)</span>
parameter.</p>
<h3 id="geometricsmoothing">Geometric smoothing</h3>
<p>For geometric smoothing, the Laplacian operator (both <span class="math">\(∆\)</span> in the continuous
setting and <span class="math">\(\L,\M\)</span> in the discrete setting) depend on the geometry of the
surface <span class="math">\(\S\)</span>. So if the signal <span class="math">\(u\)</span> is replaced with the positions of points on
the surface (say, <span class="math">\(\V ∈ \R^{n×3}\)</span> in the discrete case), then the smoothing
iteration update rule is a <em>non-linear</em> function if we write it as:</p>
<p><span class="math">\[
\M^{t+1} \V^t = (\M^{t+1} - λ\L^{t+1})\V^{t+1}.
\]</span></p>
<p>However, if we assume that small changes in <span class="math">\(\V\)</span> have a negligible effect on
<span class="math">\(\L\)</span> and <span class="math">\(\M\)</span> then we can discretize <em>explicitly</em> by computing <span class="math">\(\L\)</span> and <span class="math">\(\M\)</span>
<em>before</em> performing the update:</p>
<p><span class="math">\[
\M^{t} \V^t = (\M^{t} - λ\L^{t}) \V^{t+1}.
\]</span></p>
<h3 id="whydidmymeshdisappear">Why did my mesh disappear?</h3>
<p>Repeated application of geometric smoothing may cause the mesh to “disappear”.
Actually the updated vertex values are being set to
<a href="https://en.wikipedia.org/wiki/NaN">NaNs</a> due to degenerate numerics. We are
rebuilding the discrete Laplacian at every new iteration, regardless of the
“quality” of the mesh’s triangles. In particular, if a triangle tends to become
skinnier and skinnier during smoothing, what will happen to the cotangents of
its angles?</p>
<p>In “Can Mean-Curvature Flow Be Made Non-Singular?”, Kazhdan et al. derive a new
type of geometric flow that is stable (so long as the mesh at time <span class="math">\(t=0\)</span> is
reasonable). Their change is remarkably simple: do not update <span class="math">\(\L\)</span>, only update
<span class="math">\(\M\)</span>.</p>
<h2 id="tasks">Tasks</h2>
<h3 id="learnanalternativederivationofcotangentlaplacian">Learn an alternative derivation of cotangent Laplacian</h3>
<p>The “cotangent Laplacian” by far the most important tool in geometry
processing. It comes up everywhere. It is important to understand where it
comes from and be able to derive it (in one way or another).</p>
<p>The background section above contains a FEM derivation of the discrete
“cotangnet Laplacian”. For this (unmarked) task, read and understand one of the
<em>other</em> derivations listed above.</p>
<blockquote>
<p><strong>Hint:</strong> The finite-volume method used in [Botsch et al. 2010] is perhaps
the most accessible alternative.</p>
</blockquote>
<h3 id="whitelist">White list</h3>
<ul>
<li><code>igl::doublearea</code></li>
<li><code>igl::edge_lengths</code></li>
</ul>
<h3 id="blacklist">Black list</h3>
<ul>
<li><code>igl::cotmatrix_entries</code></li>
<li><code>igl::cotmatrix</code></li>
<li><code>igl::massmatrix</code></li>
<li>Trig functions <code>sin</code>, <code>cos</code>, <code>tan</code> etc. (e.g., from <code>#include <cmath></code>)
<em>See background notes about “intrinisic”-ness</em></li>
</ul>
<h3 id="srccotmatrix.cpp"><code>src/cotmatrix.cpp</code></h3>
<p>Construct the “cotangent Laplacian” for a mesh with edge lengths <code>l</code>. Each
entry in the output sparse, symmetric matrix <code>L</code> is given by:</p>
<p><span class="math">\[
L_{ij} = \begin{cases}
½ \cot{α_{ij}} + ½ \cot{β_{ij}} & \text{if edge $ij$ exists} \\
∑_{j≠i} L_{ij} & \text{if $i = j$} \\
0 & \text{otherwise}
\end{cases}
\]</span></p>
<blockquote>
<p>Hint: Review the <a href="https://en.wikipedia.org/wiki/Law_of_sines">law of sines</a>
and <a href="https://en.wikipedia.org/wiki/Law_of_cosines">law of cosines</a> and
<a href="https://en.wikipedia.org/wiki/Heron's_formula">Heron’s ancient formula</a> to
derive a formula for the cotangent of each triangle angle that <em>only</em> uses
edge lengths.</p>
</blockquote>
<h3 id="srcmassmatrix.cpp"><code>src/massmatrix.cpp</code></h3>
<p>Construct the diagonal(ized) mass matrix <code>M</code> for a mesh with given face indices
in <code>F</code> and edge lengths <code>l</code>.</p>
<h3 id="srcsmooth.cpp"><code>src/smooth.cpp</code></h3>
<p>Given a mesh (<code>V</code>,<code>F</code>) and data specified per-vertex (<code>G</code>), smooth this data
using a single implicit Laplacian smoothing step.</p>
<p>This data could be a scalar field on the surface and smoothing corresponds to
data denoising.</p>
<figure>
<img src="images/beetle-data-denoising.gif" alt="" />
</figure>
<p>Or the data could be the vector field of the surface’s own geometry. This
corresponds to geometric smoothing.</p>
<figure>
<img src="images/sphere-geometric-smoothing.gif" alt="" />
</figure>
</body>
</html>