forked from alecjacobson/geometry-processing-deformation
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathREADME.html
807 lines (666 loc) · 40.8 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
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8"/>
<link type="text/css" rel="stylesheet" href="shared/css/style.css"/>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "all"} } });
</script>
<script type="text/javascript" src="shared/js/MathJax/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
</head>
<body>
<div style="display:none">
<span class="math">\(\newcommand{\A}{\mat{A}}\)</span>
<span class="math">\(\newcommand{\B}{\mat{B}}\)</span>
<span class="math">\(\newcommand{\C}{\mat{C}}\)</span>
<span class="math">\(\newcommand{\D}{\mat{D}}\)</span>
<span class="math">\(\newcommand{\E}{\mat{E}}\)</span>
<span class="math">\(\newcommand{\F}{\mat{F}}\)</span>
<span class="math">\(\newcommand{\G}{\mat{G}}\)</span>
<span class="math">\(\newcommand{\I}{\mat{I}}\)</span>
<span class="math">\(\newcommand{\K}{\mat{K}}\)</span>
<span class="math">\(\newcommand{\L}{\mat{L}}\)</span>
<span class="math">\(\newcommand{\M}{\mat{M}}\)</span>
<span class="math">\(\newcommand{\N}{\mat{N}}\)</span>
<span class="math">\(\newcommand{\One}{\mathbf{1}}\)</span>
<span class="math">\(\newcommand{\P}{\mat{P}}\)</span>
<span class="math">\(\newcommand{\Q}{\mat{Q}}\)</span>
<span class="math">\(\newcommand{\Rot}{\mat{R}}\)</span>
<span class="math">\(\newcommand{\R}{\mathbb{R}}\)</span>
<span class="math">\(\newcommand{\S}{\mathcal{S}}\)</span>
<span class="math">\(\newcommand{\T}{\mat{T}}\)</span>
<span class="math">\(\newcommand{\U}{\mat{U}}\)</span>
<span class="math">\(\newcommand{\V}{\mat{V}}\)</span>
<span class="math">\(\newcommand{\W}{\mat{W}}\)</span>
<span class="math">\(\newcommand{\X}{\mat{X}}\)</span>
<span class="math">\(\newcommand{\Y}{\mat{Y}}\)</span>
<span class="math">\(\newcommand{\argmax}{\mathop{\text{argmax}}}\)</span>
<span class="math">\(\newcommand{\argmin}{\mathop{\text{argmin}}}\)</span>
<span class="math">\(\newcommand{\c}{\vec{c}}\)</span>
<span class="math">\(\newcommand{\d}{\vec{d}}\)</span>
<span class="math">\(\newcommand{\e}{\vec{e}}\)</span>
<span class="math">\(\newcommand{\f}{\vec{f}}\)</span>
<span class="math">\(\newcommand{\g}{\vec{g}}\)</span>
<span class="math">\(\newcommand{\mat}[1]{\mathbf{#1}}\)</span>
<span class="math">\(\newcommand{\min}{\mathop{\text{min}}}\)</span>
<span class="math">\(\newcommand{\n}{\vec{n}}\)</span>
<span class="math">\(\newcommand{\p}{\vec{p}}\)</span>
<span class="math">\(\newcommand{\transpose}{{\mathsf T}}\)</span>
<span class="math">\(\newcommand{\tr}[1]{\mathop{\text{tr}}{\left(#1\right)}}\)</span>
<span class="math">\(\newcommand{\t}{\vec{t}}\)</span>
<span class="math">\(\newcommand{\u}{\vec{u}}\)</span>
<span class="math">\(\newcommand{\vec}[1]{\mathbf{#1}}\)</span>
<span class="math">\(\newcommand{\x}{\vec{x}}\)</span>
<span class="math">\(\newcommand{\y}{\vec{y}}\)</span>
<span class="math">\(\newcommand{\z}{\vec{z}}\)</span>
<span class="math">\(\renewcommand{\v}{\vec{v}}\)</span>
</div>
<h1 id="geometryprocessing–deformation">Geometry Processing – Deformation</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-deformation.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:</p>
<pre><code>./deformation [path to mesh.obj]
</code></pre>
<p>When the mesh is blue, the system is in “place handles” mode. Click on the mesh
to select vertex locations for control point handles. After pressing space to switch to
deformation mode, drag the handles. Pressing <code>m</code> will switch between different
deformation methods.</p>
<figure>
<img src="images/knight-deformation.gif" alt="" />
</figure>
<h2 id="background">Background</h2>
<p>In this assignment we explore smooth deformation of an existing shape. <a href="https://en.wikipedia.org/wiki/Deformation_(mechanics)">Shape
deformation</a> has many
applications in geometry process; we will focus on interactive <strong><em>handle-based
deformation</em></strong>. In this setup, the user repositions a sparse set of points and the
goal is to propagate the transformation at these “handles” to the rest of the
shape. To be interactive we should aim for computing the new deformation of the
shape at around 30 <a href="https://en.wikipedia.org/wiki/Frame_rate">frames per
second</a>.</p>
<p>Shape deformation is the transformation from its <em>rest shape</em> to a
new/current/deformed shape. If the position of a point on some 3D rest shape
given by <span class="math">\(\hat{\x} ∈ \R³\)</span> then we will write that the unknown position
on the deformed shape is given by <span class="math">\(\x ∈ \R³\)</span>. We can write this point’s
displacement vector as <span class="math">\(\x - \hat{\x} =: \d ∈ \R³\)</span>.</p>
<p>The propagation of the handles’ deformation can be thought of in two
complementary ways:</p>
<ol>
<li>as a <a href="https://en.wikipedia.org/wiki/Multivariate_interpolation#Irregular_grid_.28scattered_data.29">scattered data
interpolation</a>
problem, where handles provide sparse samples of an unknown <a href="https://en.wikipedia.org/wiki/Displacement_field_(mechanics)">displacement
field</a>, or</li>
<li>as a shape optimization problem, where we try to define a <em>new</em> shape that
retains the details of the old shape but fulfills the handle constraints.</li>
</ol>
<p>In the following discussion we will take advantage of the ability to switch
between thinking of the unknowns as the positions of the deformed shape (<span class="math">\(\x\)</span>)
and displacements (<span class="math">\(\d\)</span>). These views are equivalent, but often one or the
other provides a better intuitive understanding.</p>
<h4 id="continuity">Continuity</h4>
<p>We will limit ourselves to <em>continuous</em> deformations of shapes. That is, the
shape will not tear, crack or change its topological features. </p>
<p>If we represent our shape discretely as a <a href="https://en.wikipedia.org/wiki/Triangle_mesh">triangle
mesh</a> (e.g., with <em>rest</em> vertices
in <span class="math">\(\hat{\V}
∈ \R^{n × 3}\)</span> and faces in <span class="math">\(F ∈ \{1,…,n\}^{m × 3}\)</span>, then we can <em>trivially</em>
ensure a continuous deformation by determining new vertex positions <span class="math">\(\V\)</span>. The
topology (connectivity) of the mesh (<span class="math">\(F\)</span>) will not change.</p>
<h4 id="genericdistortionminimization">Generic Distortion Minimization</h4>
<p>A rest surface <span class="math">\(\hat{S}\)</span>
<a href="https://en.wikipedia.org/wiki/Immersion_(mathematics)">immersed</a> in <span class="math">\(\R³\)</span> can
be described as a mapping <span class="math">\(\hat{\x}\)</span> from <em>some</em> 2D parametric domain <span class="math">\(Ω\)</span>. For any
parameters <span class="math">\(u\)</span> and <span class="math">\(v\)</span>, <span class="math">\(\hat{\x}\)</span> describes the 3D position:</p>
<p><span class="math">\[
\hat{\x}(u,v) ∈ \R³.
\]</span></p>
<p>Similarly the deformed surface can be represented as a position function <span class="math">\(\x:
Ω → \R³\)</span>. The displacement vector field is thus a function taking the
difference: <span class="math">\(\d(u,v) = \x(u,v) - \hat{\x}(u,v)\)</span>.</p>
<p>For the handle-based deformation problem we would like to find a new surface
(defined by <span class="math">\(\x\)</span>) that:</p>
<ol>
<li>adds as little <em>distortion</em> as possible to the shape, and</li>
<li>satisfies the users constraints at selected handle positions.</li>
</ol>
<p>We can cast this as an energy optimization problem. Suppose we have
energy <a href="https://en.wikipedia.org/wiki/Functional_(mathematics)">functional</a>
<span class="math">\(E(\x)\)</span> that
measures the amount of distortion between the new shape (<span class="math">\(\x\)</span>) and the rest
shape <span class="math">\(\hat{x}\)</span>, then we could optimize for the best possible shape <span class="math">\(\x\)</span> by
minimizing <span class="math">\(E\)</span>:</p>
<p><span class="math">\[
\min_\x E(\x).
\]</span></p>
<p>To ensure that the user’s <span class="math">\(k\)</span> handle points are interpolated we add the
constraints:</p>
<p><span class="math">\[
\text{ subject to } \x(u_i, v_i) = \g_i \ ∀ i = \{1, … , k\},
\]</span>
where <span class="math">\(\g_i\)</span> is the position of the $i$th control point handle.</p>
<p>While the constraints are straightforward, we have many choices for how to
formulate the energy function <span class="math">\(E\)</span>. A natural choice is to measure distortion in
an egalitarian way by integrating a <em>local</em> measure of distortion at all points
on the surface:</p>
<p><span class="math">\[
\min_\x ∫_Ω ‖ e(\x) ‖² \ dA \quad \text{ subject to } \x(u_i, v_i) = \g_i \ ∀ i = \{1, … , k\},
\]</span></p>
<p>where <span class="math">\(e\)</span> is a vector- or scalar- valued function measuring local (unsquared)
distortion. We will now consider different choices for <span class="math">\(e\)</span>.</p>
<h3 id="linearmethods">Linear Methods</h3>
<p>If we assume that the deformation between the rest shape given by <span class="math">\(\hat{\x}\)</span>
and the new shape given by <span class="math">\(\x\)</span> is <em>small</em> then we can measure the distortion
of the deformation in terms of the smoothness of the displacement field. This
simplest methods will integrate the magnitude of derivatives of the
displacement field (<span class="math">\(\d\)</span>): if the displacement field has large variations or
sudden changes then it is inducing a lot of distortion.</p>
<h4 id="gradient-basedenergy">Gradient-based energy</h4>
<p>Let us first consider minimizing the integral of squared variation of the
displacement field:</p>
<p><span class="math">\[
\min_\d ∫_Ω ‖ ∇\d ‖_F^2 \ dA \quad \text{ subject to } \d_i =
\g_i-\hat{\x}_i \ ∀ i = \{1, … , k\},
\]</span>
where
<span class="math">\(∇\d = \left(\begin{array}{ccc}
\frac{∂d^x}{∂u} & \frac{∂d^y}{∂u} & \frac{∂d^z}{∂u} \\
\frac{∂d^x}{∂v} & \frac{∂d^y}{∂v} & \frac{∂d^z}{∂v}
\end{array}\right)\)</span>
the <a href="https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant">Jacobian</a>
matrix of the displacement field <span class="math">\(\d\)</span> </p>
<blockquote>
<h5 id="deformationgradient">Deformation Gradient</h5>
<p>If <span class="math">\(\I ∈ \R^{3 × 3}\)</span> is the identity matrix, then the quantity <span class="math">\(\F := \I +
∇\d\)</span> is referred to as the <a href="https://en.wikipedia.org/wiki/Finite_strain_theory#Deformation_gradient_tensor">deformation
gradient</a>
in the <a href="https://en.wikipedia.org/wiki/Continuum_mechanics">mechanics</a>
community.</p>
</blockquote>
<p>This is simply the familiar <a href="https://en.wikipedia.org/wiki/Dirichlet's_energy">Dirichlet
energy</a> applied to each
coordinate function of the displacement field <em>independently</em>.</p>
<p>We can discretize this over our triangle mesh surface the same way we have in
smoothing and parameterization assignments:</p>
<p><span class="math">\[
\min_{\D} \tr{\D^\transpose \L \D} \quad \text{subject to } \D_\text{handles} =
\g_\text{handles} - \hat{\V}_\text{handles},
\]</span>
where the rows of <span class="math">\(\g_\text{handles} ∈ \R^{k × 3}\)</span> contains the new positions
of the <span class="math">\(k\)</span> control point handles.</p>
<p>While easy to implement, this method suffers from a couple immediate problems:</p>
<ol>
<li>it is not smooth at constraints, and</li>
<li>then <em>influence</em> of handles dies off too quickly.</li>
</ol>
<figure>
<img src="images/bump-k-harmonic.jpg" alt="Two regions of selected points (purple) are constrained. The gradient-based
energy produces a continuous displacement but is not smooth at the inner
constrained points (sharp crease). The Laplacian-based energy produces a
displacement with continuous positions and
derivatives." />
<figcaption>Two regions of selected points (purple) are constrained. The gradient-based
energy produces a <em>continuous</em> displacement but is not <em>smooth</em> at the inner
constrained points (sharp crease). The Laplacian-based energy produces a
displacement with continuous positions and
derivatives.</figcaption>
</figure>
<p>By minimizing the Dirichlet energy, each coordinate of “displacement field” is
a <a href="https://en.wikipedia.org/wiki/Harmonic_function">harmonic function</a>.
Intuitively (however abstractly) we can think of each function as <em>diffusing</em>
the user’s constraints as if they were
<a href="https://en.wikipedia.org/wiki/Heat_equation">heat</a> values. As such, each
function diffuses quickly to an average, slowly varying value over most of the
domain. As a displacement field a constant value for each coordinate function
would mean a translation: most of the shape is simply translated.</p>
<figure>
<img src="images/bunny-harmonic.gif" alt="" />
</figure>
<p>The gradient operator (<span class="math">\(∇\)</span>) is a <a href="https://en.wikipedia.org/wiki/Linear_map"><em>linear</em>
operator</a>. We can alternatively view
our minimization above in terms of the unknown positions <span class="math">\(\x\)</span>:</p>
<p><span class="math">\[
\min_\d ∫_Ω ‖ ∇\d ‖_F^2 \ dA ⇒
\min_\x ∫_Ω ‖ ∇(\x - \hat{\x}) ‖_F^2 \ dA ⇒
\min_\x ∫_Ω ‖ \underbrace{∇\x}_\text{after} -
\underbrace{∇\hat{\x}}_\text{before} ‖_F^2 \ dA.
\]</span>
If we think of the gradient of the position function <span class="math">\(∇\x\)</span> (with respect to the
underlying parameterization <span class="math">\(u,v\)</span>) as a local geometric <a href="https://en.wikipedia.org/wiki/Feature_(computer_vision)">feature
descriptor</a> then this
energy can be re-understood as measuring the difference in this feature before
and after the deformation. This is very sensible as we are trying to measure
distortion. We would expect that a low-distortion deformation would correspond
with a small change to local features.</p>
<p>Unfortunately, the gradient of the position function <span class="math">\(\x\)</span> is a <em>poor</em>,
first-order local feature.</p>
<h4 id="laplacian-basedenergy">Laplacian-based energy</h4>
<p>If we model distortion as the change in a local feature descriptor, then a
natural local and <em>relative</em> descriptor would be one that compared the position
of some point on the shape to the average of its local neighborhood. We have
studied an operator that computes this in the smoothing assignment. The
<a href="https://en.wikipedia.org/wiki/Laplace_operator">Laplace(-Beltrami) operator</a>
can be derived as taking exactly the difference of a functions value at a point
and the average (i.e., <a href="https://en.wikipedia.org/wiki/Centroid">centroid</a>) of
an infinitesimal region around that point:</p>
<p><span class="math">\[
∆ f(\x) = \lim_{|B(\x)| → 0} \frac{1}{|B(\x))|} ∫_{B(\x)} f(\z) \;d\z - f(\x)
\]</span></p>
<p>(see, e.g., “Differential coordinates for local mesh morphing and deformation”
[Alexa et al. 2003] and expanded upon in “Laplacian Surface Editing”
[Sorkine et al. 2004]).</p>
<p>When applied to the embedding function <span class="math">\(\x\)</span> the Laplace operator computes the
difference in <em>position</em> between a point and its local neighborhood. This
<em>vector</em> points in the
<a href="https://en.wikipedia.org/wiki/Normal_(geometry)">normal</a> direction and its
magnitude corresponds inversely with <a href="https://en.wikipedia.org/wiki/Mean_curvature">how flat the surface is
locally</a>.</p>
<p>Let’s replace the gradient feature above with this <em>second-order</em> feature
descriptor and massage our optimization problem back in terms of
displacements:</p>
<p><span class="math">\[
\min_\x ∫_Ω ‖ \underbrace{∆\x}_\text{after} -
\underbrace{∆\hat{\x}}_\text{before} ‖^2 \ dA ⇒
\min_\x ∫_Ω ‖ ∆(\x - \hat{\x}) ‖^2 \ dA ⇒
\min_\d ∫_Ω ‖ ∆\d ‖^2 \ dA.
\]</span></p>
<p>Just as we can show that harmonic functions (<span class="math">\(∆\d = 0\)</span>) minimize the Dirichlet
energy, we can use <a href="https://en.wikipedia.org/wiki/Calculus_of_variations">calculus of
variations</a> apply
<a href="https://en.wikipedia.org/wiki/Green's_identities">Green’s identity</a> <em>twice</em> to
show that minimizers of the squared-Laplacian energy are _bi-_harmonic
functions (<span class="math">\(∆∆\d = 0\)</span> or <span class="math">\(∆²\d = 0\)</span>). Obviously all harmonic functions are also
biharmonic functions. This implies that the space of biharmonic functions is
strictly <em>larger</em>. In particular, this will allow use to ensure continuity of
first derivatives across constrained values at handles.</p>
<figure>
<img src="images/bunny-biharmonic.gif" alt="" />
</figure>
<p>The fact that each coordinate of the displacement field <span class="math">\(\d\)</span> will be a
bi-harmonic function is not so elucidating, but by minimizing the
squared-Laplacian energy integrated over the domain, we can say that it is
as-<em>harmonic</em>-as-possible. Harmonic functions include constant functions (i.e.,
translations) but also any displacement that <em>bends</em> in one direction by an
equal and opposite amount as it bends in the other direction: <span class="math">\(∆\d = 0 → ∂\d/∂u
+ ∂\d/∂v = 0 → ∂\d/∂u = -∂\d/∂v\)</span>.</p>
<p>To discretize this energy we can make use of our discrete Laplacian operator
<span class="math">\(\L\)</span>. This matrix computes the locally <em>integrated</em> Laplacian of a given
function specified by per-vertex values <span class="math">\(\f\)</span>. Now we would like to integrate
the square of the <em>point-wise</em> Laplacian. We can approximate the point-wise
Laplacian by the local <a href="https://en.wikipedia.org/wiki/Mean_of_a_function">integral
average</a> of the Laplacian
<span class="math">\(\M^{-1}\L \f\)</span>. Integrating this over the mesh we have our complete approximate
of the energy:</p>
<p><span class="math">\[
\begin{align}
∫_Ω ‖∆\d‖² \;dA &≈ \tr{ \D^\transpose \L^\transpose \M^{-\transpose} \M \M^{-1} \L
\D}\\
&= \tr{ \D^\transpose \underbrace{\L^\transpose \M^{-1} \L}_{\Q} \D },
\end{align}
\]</span>
where <span class="math">\(\M ∈ \R^{n × n}\)</span> is the mass-matrix for the given mesh and <span class="math">\(\Q ∈
\R^{n×n}\)</span> can be thought of as the bi-Laplacian matrix.</p>
<blockquote>
<h4 id="k-harmonic"><span class="math">\(k\)</span>-harmonic</h4>
<p>The logical continuation of harmonic and biharmonic deformation is to consider
<em>triharmonic</em> (<span class="math">\(∆³\d = 0\)</span>) and <em>tetraharmonic</em> (<span class="math">\(∆⁴\d = 0\)</span>) and so on. It’s
straightforward to implement these, though there are diminishing returns and
increasing costs.</p>
</blockquote>
<h5 id="precomputation">Precomputation</h5>
<p>With out loss of generality, assume that the rows of the unknown displacements
<span class="math">\(\D\)</span> have been sorted so that displacements corresponding to handle vertices
are in the bottom part:</p>
<p><span class="math">\[
\D = \left(\begin{array}{c}
\D_\text{u} \\
\D_\text{h}
\end{array}
\right)
\]</span></p>
<p>Since the displacements at handles are <em>known</em> before the optimization, we can
separate the knowns and unknowns in the energy:</p>
<p><span class="math">\[
\min_{\D_\text{u}}
\tr{(\D_\text{u}^\transpose \ \D_\text{h}^\transpose)
\left(\begin{array}{cc}
\Q_\text{u,u} & \Q_\text{u,h} \\
\Q_\text{h,u} & \Q_\text{h,h}
\end{array}\right)
\left(\begin{array}{c}
\D_\text{u} \\
\D_\text{h}
\end{array}
\right)} \\
\min_{\D_\text{u}}
\tr{\D_\text{u}^\transpose \Q_\text{u,u} \D_\text{u} +
2 \D_\text{u}^\transpose \Q_\text{u,h} \D_\text{h} +
\underbrace{\D_\text{h}^\transpose \Q_\text{h,h}
\D_\text{h}}_\text{constant}} \\
\min_{\D_\text{u}}
\tr{
\D_\text{u}^\transpose \Q_\text{u,u} \D_\text{u} +
2 \D_\text{u}^\transpose \Q_\text{u,h} \D_\text{h}}
\]</span>
where <span class="math">\(\Q_\text{u,u} ∈ \R^{(n-k) × (n-k)}\)</span> is the quadratic coefficients matrix
corresponding to the unknown displacements.</p>
<p>This quadratic optimization problem may solved by setting all partial
derivatives with respect to degrees of freedom in <span class="math">\(\D_\text{u}\)</span> to zero:</p>
<p><span class="math">\[
2 \Q_\text{u,u} \D_\text{u} + 2 \Q_\text{u,h} \D_\text{h}
= 0 → \D_\text{u} = \Q_\text{u,u}^{-1} \Q_\text{u,h} \D_\text{h}
\]</span></p>
<p>If we don’t change <em>which</em> vertices are handles, but only change the positions
of the selected handles, then only <span class="math">\(\D_\text{h}\)</span> changes above. In particular,
the matrix <span class="math">\(\Q_\text{u,u}\)</span> is unchanged. Therefore, we can
<a href="https://en.wikipedia.org/wiki/Cholesky_decomposition">prefactorize</a> it so that
<em>applying its inverse</em> is fast (<code>igl::min_quad_with_fixed</code> does this for you).</p>
<blockquote>
<p>Actually the entire term <span class="math">\(\Q_\text{u,u}^{-1} \Q_\text{u,h} =: \W ∈ \R^{(n-k)
× k}\)</span> does not change. The columns of <span class="math">\(\W\)</span> reveal how unknown displacements
respond to each handle point’s displacement (see also “An intuitive framework
for real-time freeform modeling” [Botsch & Kobbelt 2004]). This provides a
gateway to the relationship with linear blend skinning and automatic
(biharmonic) weighting functions (see “Bounded Biharmonic Weights for
Real-Time Deformation” [Jacobson et al. 2011]).</p>
</blockquote>
<h5 id="troubleinparadise">Trouble in paradise</h5>
<p>Biharmonic displacements work well for small deformations and deformations that
do not imply a large rotation. However, the Laplacian of the position function
<span class="math">\(∆\x\)</span> as a feature descriptor is <em>not</em> rotation invariant. This problem is true
for all linear differential features including the gradient of the embedding
function <span class="math">\(∇\x\)</span> considered above (see “On linear variational surface deformation
methods” [Botsch & Sorkine 2008]).</p>
<figure>
<img src="images/knight-biharmonic-large-rotation.gif" alt="" />
</figure>
<p>This means that if the user transforms all of the handle locations by a rigid
transformation <span class="math">\(\T\)</span> these energies will <em>not</em> measure zero for a displacement
equivalent to applying the rigid transformation <span class="math">\(\T\)</span> to the entire shape. We
would like <em>global</em> rotation invariance, but we would also like this property
to apply <em>locally</em> so that parts of the shape can easily rotate.</p>
<h3 id="as-rigid-as-possible">As-rigid-as-possible</h3>
<p>In the scenario where each handle
<span class="math">\(i\)</span> are perfectly transformed by a <a href="https://en.wikipedia.org/wiki/Rigid_transformation">rigid
transformation</a> <span class="math">\(\x_i =
\Rot \hat{\x}_i + \t\)</span>, where <span class="math">\(\Rot ∈ SO(3) ⊂ \R^{3×3}\)</span> is a rotation matrix and
<span class="math">\(\t∈\R^3\)</span> is a translation vector. If an
<a href="https://en.wikipedia.org/wiki/Oracle">oracle</a> could only tell us this
particular rigid transformation then we could repair
the gradient-based energy above by pre-rotating the rest shape by this
transformation:</p>
<p><span class="math">\[
\begin{align}
∫_Ω ‖ ∇ \x - ∇(\Rot \hat{\x} + \t) ‖² \;dA
&= ∫_Ω ‖ ∇ \x - ∇(\Rot \hat{\x}) - ∇\t‖² \;dA \\
&= ∫_Ω ‖ ∇ \x - \Rot ∇\hat{\x} ‖² \;dA,
\end{align}
\]</span></p>
<p>where the translation vector <span class="math">\(\t\)</span> falls out because a translation has constant
gradient.</p>
<p>We do not know the rotation <span class="math">\(\Rot\)</span> ahead of time, but we could be as generous
as possible and use the “best” rotation <span class="math">\(\Rot ← \argmin_\Rot ∫_Ω ‖ ∇\x - \Rot
∇\hat{\x} ‖² \;dA\)</span>:</p>
<p><span class="math">\[
∫_Ω \left\|∇\x - \left( \argmin_\Rot ∫_Ω ‖ ∇\x - \Rot ∇\hat{\x} ‖² \;dA \right)∇\hat{\x}\right\|² \;dA.
\]</span></p>
<p>If we treat <span class="math">\(\Rot\)</span> as a degree of freedom along with the unknown positions
<span class="math">\(\x\)</span>, we can unify this into an optimization over <span class="math">\(\x\)</span> and <span class="math">\(\Rot\)</span>:</p>
<p><span class="math">\[
\min_{\x,\Rot∈SO(3)} ∫_Ω \left\|∇\x - \Rot ∇\hat{\x}\right\|² \;dA.
\]</span></p>
<p>Optimizing this energy will ensure <em>global</em> rotation invariance. To ensure
<em>local</em> rotation invariance, we can replace <span class="math">\(\Rot ∈ SO(3)\)</span> with a spatially
varying <em>function</em> <span class="math">\(\Rot : Ω → SO(3)\)</span> that outputs a “best” rotation for any
point on the shape (see “A simple geometric model for elastic deformations”
[Chao et al. 2010]). In this way, the optimal rotation will be locally rigid
everywhere, or <em>as-rigid-as-possible</em> (ARAP).</p>
<figure>
<img src="images/knight-arap-large-rotation.gif" alt="" />
</figure>
<blockquote>
<p>For embedded solid shapes, we can take the rest shape given by <span class="math">\(\hat{\x}\)</span> as
the parameterization <span class="math">\(Ω\)</span>, so that <span class="math">\(∇\hat{\x} = \I\)</span>. This allows us to rewrite
the as-rigid-as-possible energy as the square of the difference between the
<a href="#deformationgradient">deformation gradient</a> and the closest rotation:</p>
<p><span class="math">\[
∫_Ω \left\|∇\x - \Rot ∇\hat{\x}\right\|² \;dA \\
∫_Ω \left\|(∇\x + \I - \I) - \Rot \I \right\|² \;dA \\
∫_Ω \left\|(\I + ∇\x - ∇\hat{x}) - \Rot \right\|² \;dA \\
∫_Ω \left\|(\I + ∇\d) - \Rot \right\|² \;dA \\
∫_Ω \left\|\F - \Rot \right\|² \;dA \\
\]</span></p>
<p>This form provides a bridge between the as-rigid-as-possible energy common in
geometry processing to <em>corotated linear elasticity</em> used in
mechanics/physically-based simulation (made explicit in “A simple geometric
model for elastic deformations” [Chao et al. 2010]). See Section 3.4 of “FEM
Simulation of 3D Deformable Solids” [Sifakis 2012] for a graphics-mechanics
perspective.</p>
</blockquote>
<h4 id="discreteas-rigid-as-possibleenergy">Discrete as-rigid-as-possible energy</h4>
<p>For a triangle mesh with displacing vertices, the gradient of the embedding
function is constant inside each triangle. In this way we can write the raw
gradient energy above as a double sum over all half-edges <span class="math">\(ij\)</span> of all faces <span class="math">\(f\)</span>
in the mesh:</p>
<p><span class="math">\[
½ ∫_Ω ‖ ∇ \x - ∇\hat{\x}‖² \;dA = ½ ∑\limits_{f ∈ F} ∑\limits_{ ij ∈ f} c_{ij} ‖
(\v_i-\v_j) - (\hat{\v}_i-\hat{\v}_j)‖²,
\]</span>
where <span class="math">\(c_{ij}\)</span> is cotangent of the angle opposite half-edge <span class="math">\(ij\)</span>.</p>
<p>To inject localized best fit rotations, we will assign an unknown rotation
matrix <span class="math">\(\Rot_k\)</span> to each vertex <span class="math">\(k\)</span> of the mesh and accounts for a third of the
energy integrated over incident triangles:
<span class="math">\[
½ ∫_Ω ‖ ∇ \x - \Rot ∇\hat{\x}‖² \;dA =
⅙ ∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)}
c_{ij} ‖ (\v_i-\v_j) - \Rot_k (\hat{\v}_i-\hat{\v}_j)‖²,
\]</span>
where <span class="math">\(F(k)\)</span> is the set of all faces incident on the <span class="math">\(k\)</span>-th vertex.</p>
<blockquote>
<h5 id="wheredorotationslive">Where do rotations <em>live</em>?</h5>
<p>We have assigned a rotation for each vertex: there are <span class="math">\(n\)</span> rotation matrices
as auxiliary degrees of freedom. This is in contrast to assigning
rotations per face (as in “A Local/Global Approach to Mesh Parameterization”
[Liu et al. 2008]). Per-face–or more generally per-element–rotations work
well for <a href="https://en.wikipedia.org/wiki/Codimension">codimension</a> zero
objects (triangle meshes in <span class="math">\(\R²\)</span> or tetrahedral meshes in <span class="math">\(\R³\)</span>). But for
triangle mesh surfaces in <span class="math">\(\R³\)</span> (i.e., codimension 1), per-face rotations
would lead to “crumpling”) because <em>bending</em> along edges would not be
measured.</p>
</blockquote>
<h4 id="optimization">Optimization</h4>
<p>The simplest method for optimizing the ARAP energy is by alternating between</p>
<ol>
<li>finding the optimal rotations <span class="math">\(\Rot_k\)</span> assuming the vertex positions <span class="math">\(\V\)</span> are
fixed, and</li>
<li>finding the optimal vertex positions <span class="math">\(\V\)</span> assuming all rotations <span class="math">\(\Rot_k\)</span> are
fixed.</li>
</ol>
<p>Each rotation <span class="math">\(\Rot_k\)</span> only affects the local energy and doesn’t interact with
the <em>other</em> rotations. So each can be optimized <em>locally</em>. In contrast, the
mesh vertex positions <span class="math">\(\V\)</span> depend on each other requiring a <em>global</em> solve. In
the geometry processing literature, this is known as a local-global
optimization (see “As-rigid-as-possible surface modeling” [Sorkine & Alexa
2007]). It is also known as “alternating block coordinate descent” because we
have separated the variables into disjoint sets (<span class="math">\(\V,\Rot_1,…,\Rot_n\)</span>) and taking
the optimal descent direction for each independently.</p>
<p>Observing the discrete energy above we can see that the energy is quadratic in
<span class="math">\(\V\)</span> and quadratic in each <span class="math">\(\Rot_k\)</span>. Let’s start by separating the terms that are
quadratic and linear in <span class="math">\(\V\)</span>:</p>
<p><span class="math">\[
⅙ \underbrace{∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)} c_{ij} (\v_i-\v_j)^\transpose(\v_i-\v_j)}_\text{quadratic}
+
⅙ \underbrace{∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)} c_{ij} (\v_i-\v_j)^\transpose \Rot_k (\hat{\v}_i-\hat{\v}_j)}_\text{linear}
\]</span></p>
<p>if we stack the rotation matrices <span class="math">\(\Rot_k\)</span> into large matrix <span class="math">\(\Rot ∈ \R^{3n ×
3}\)</span> then we can write this energy in matrix form as:</p>
<p><span class="math">\[
\tr{ \V^\transpose \L \V } + \tr{ \V^\transpose \K \Rot },
\]</span>
where <span class="math">\(\L ∈ \R^{n × n}\)</span> is the familiar cotangent discrete Laplacian matrix and
<span class="math">\(\K ∈ \R^{n × 3n}\)</span> sparse matrix containing cotangents multiplied against
differences across edges in the rest mesh (e.g., <span class="math">\(\hat{\v}_i - \hat{\v}_j\)</span>).</p>
<blockquote>
<h6 id="imsoconfused.whatsinthekmatrix">I’m so confused. What’s in the <span class="math">\(\K\)</span> matrix?</h6>
<p>Let’s take it slow. The <span class="math">\(\K\)</span> matrix is represents the
<a href="https://en.wikipedia.org/wiki/Bilinear_form">bilinear form</a> that combines unknown
vertex positions and unknown rotations. We have identified above that we can
write this in summation or matrix form:
<span class="math">\[
⅙ ∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)} c_{ij} (\v_i-\v_j)^\transpose \Rot_k (\hat{\v}_i-\hat{\v}_j)
= \tr{ \V^\transpose \K \Rot },
\]</span>
but how did we get here?</p>
<p>Let’s start with the summation form. The <em>constants</em> of this formula are the
<span class="math">\(c_{ij}\)</span> terms and the <span class="math">\((\hat{\v}_i-\hat{\v}_j)\)</span> terms. Since these always
appear together, let us merge them into weighted edge difference vectors
<span class="math">\(c_{ij} (\hat{\v}_i-\hat{\v}_j) =: \hat{\e}_{ij} ∈ \R³\)</span>:</p>
<p><span class="math">\[
⅙ ∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)} \underbrace{(\v_i-\v_j)^\transpose
\Rot_k \hat{\e}_{ij}}_{∈\R},
\]</span>
the inner term in the summation is an <a href="https://en.wikipedia.org/wiki/Inner_product_space">inner
product</a>; that is, a
<a href="https://en.wikipedia.org/wiki/Scalar_(mathematics)">scalar</a>. Let’s expose this
by expanding the matrix-vector products of the inner-product:
<span class="math">\[
⅙ ∑\limits_{k=1}^n ∑\limits_{ ij ∈ F(k)} ∑\limits_{α=1}^3 ∑\limits_{β=1}^3
(v_i^α - v_j^α)R_k^{αβ}\hat{e}_{ij}^β.
\]</span></p>
<p>If our mesh is stored as a vertex list and face list, it’s not easy/efficient
to loop over per-vertex rotations (outer sum) and then over all half-edges of
incident faces (second sum). Instead, let’s rearrange these sums to loop over
all faces first, then the half-edges of that face, and then over all
per-vertex rotations that involve this half-edge:
<span class="math">\[
⅙ ∑\limits_{f=1}^m ∑\limits_{ ij ∈ E(f)} \quad ∑\limits_{k | ij ∈F(k)} \quad ∑\limits_{α=1}^3 ∑\limits_{β=1}^3
(v_i^α - v_j^α)R_k^{αβ}\hat{e}_{ij}^β,
\]</span>
where the third sum is over all rotations <span class="math">\(k\)</span> such that the half-edge <span class="math">\(ij\)</span>
belongs to the half-edges of the faces incident on the <span class="math">\(k\)</span>-th vertex: <span class="math">\(k | ij
∈F(k)\)</span>. Well, this means <span class="math">\(k\)</span> can either be <span class="math">\(i\)</span> or <span class="math">\(j\)</span> or the third vertex of
the <span class="math">\(f\)</span>-th face.</p>
<p>Now let’s turn our attention back to the
<a href="https://en.wiktionary.org/wiki/summand">summand</a>. The terms indexed by <span class="math">\(α\)</span>
never <em>mix</em>. That is, we never add/multiply <span class="math">\(v_i^α\)</span>, <span class="math">\(v_j^γ\)</span>, and <span class="math">\(R_k^{δβ}\)</span>
unless <span class="math">\(α=γ=δ\)</span>. This implies that we can write this summation in matrix form
as:
<span class="math">\[
\V₁^\transpose \K₁ \Rot₁ +
\V₂^\transpose \K₂ \Rot₂ +
\V₃^\transpose \K₃ \Rot₃,
\]</span>
where <span class="math">\(\V_α ∈ \R^n\)</span> is <span class="math">\(α\)</span>-th column of <span class="math">\(\V\)</span>, <span class="math">\(\Rot_α ∈ \R^{3n}\)</span> is the <span class="math">\(α\)</span>-th column of
<span class="math">\(\Rot\)</span> and <span class="math">\(\K₁,\K₂,\K₃ ∈ \R^{n × 3n}\)</span> are sparse matrices.</p>
<p><em>Further</em>, the constant term <span class="math">\(\hat{e}^β_{ij}\)</span> in the summation acts the same on
<span class="math">\((v_i^α - v_j^α)R_k^{αβ}\)</span> for any <span class="math">\(α\)</span> value. This implies that <span class="math">\(\K₁=\K₂=\K₃\)</span>,
so we can reduce the matrix form to:
<span class="math">\[
\tr{\V \K \Rot}.
\]</span></p>
<p>Finally, we can answer what is in each entry <span class="math">\(K_{vw}\)</span>, where <span class="math">\(v\)</span> chooses the
row and <span class="math">\(w=3k+β\)</span> chooses the column of <span class="math">\(\K\)</span>. Treating our nested summation as
nested for-loops, we can increment entries in <span class="math">\(\K\)</span>
<span class="math">\[
\begin{align}
K_{i\ 3k+β} & \mathrel{+}= \hat{e}_{ij}^β, \\
K_{j\ 3k+β} & \mathrel{+}= -\hat{e}_{ij}^β, \\
\end{align}
\]</span>
for each half-edge encountered. </p>
</blockquote>
<h5 id="localstep">Local step</h5>
<p>Minimizing this energy with respect <span class="math">\(\Rot\)</span> corresponds to minimizing:</p>
<p><span class="math">\[
\tr{ \underbrace{\V^\transpose \K}_{\C^\transpose} \Rot },
\]</span>
where <span class="math">\(\C ∈ \R^{3n × 3}\)</span> stacks weighted covariance matrices <span class="math">\(\C_k ∈ \R^{3 ×
3}\)</span> for each region <em>covered</em> by the corresponding rotation <span class="math">\(\Rot_k\)</span>. We have
seen this problem before in the registration assignment. For each <span class="math">\(\C_k\)</span>,
<span class="math">\(\Rot_k\)</span> will be the closest rotation matrix solved via <a href="https://en.wikipedia.org/wiki/Singular_value_decomposition">singular value
decomposition</a>.</p>
<h5 id="globalstep">Global step</h5>
<p>Minimizing the energy above with respect to <span class="math">\(\V\)</span> corresponds to solving a
Dirichlet energy-like minimization problem:</p>
<p><span class="math">\[
\min_\V \tr{ \V^\transpose \L \V } + \tr{ \V^\transpose \B }
\]</span></p>
<p>where <span class="math">\(\K \Rot =: \B ∈ \R^{n × 3}\)</span> is a matrix of rotated vertex gradients.
Adding the handle constraints to the corresponding rows of <span class="math">\(\V\)</span> this is easily
minimized by setting all partial derivatives with respect to the unknowns in
<span class="math">\(\V\)</span> equal to zero (as in the linear methods above).</p>
<h5 id="implementation">Implementation</h5>
<p>In order to facilitate interactive deformation we would like our local and
global iterations to be computed as quickly as possible. Since the quadratic
form in the global step is the <em>same</em> regardless of the current rotations or
current handle positions, we can
<a href="https://en.wikipedia.org/wiki/Cholesky_decomposition">prefactorize</a> it (again,
as above). The matrix <span class="math">\(\K\)</span> also does not depend on the rotations, current
positions or handle positions, so we can pre-build this matrix. This way
computing <span class="math">\(\C\)</span> and <span class="math">\(\B\)</span> is a simple matrix-matrix multiplication.</p>
<blockquote>
<p><strong>Note:</strong> When constructing <span class="math">\(\K\)</span> it’s easiest to iterate over <em>all</em>
half-edges in the mesh (by iterating over all faces and then each of the
three edges). Each half-edge <span class="math">\(ij\)</span> <em>contributes</em> terms tying <span class="math">\(\v_i,\v_j\)</span> to
<em>each</em> of the (three) rotations <span class="math">\(\Rot_k\)</span> that apply against their difference
(see “Fast Automatic Skinning Transformations” [Jacobson et al. 2012]).</p>
</blockquote>
<h5 id="stilltroubleinparadise">Still trouble in paradise</h5>
<figure>
<img src="images/knight-arap-high-vs-low-resolution.gif" alt="" />
</figure>
<p>The as-rigid-as-possible deformation method for surfaces described above has a
number of remaining problems:</p>
<ol>
<li>like its gradient-based energy ancestor the deformation is not smooth at
constraints;</li>
<li>the energy punishes <em>bending</em> of the surface–which is good–but does so in
a way that diminishes as the mesh becomes higher and higher resolution, in
otherwords, the discrete energy is mesh-resolution dependent; and</li>
<li>the energy is biased by the original combinatorics of the mesh (even in
flat regions).</li>
</ol>
<h2 id="tasks">Tasks</h2>
<h3 id="blacklist">Blacklist</h3>
<ul>
<li><code>igl::arap</code></li>
<li><code>igl::arap_linear_block</code></li>
<li><code>igl::covariance_scatter_matrix</code></li>
<li><code>igl::harmonic</code></li>
</ul>
<h3 id="whitelist">Whitelist</h3>
<ul>
<li><code>igl::polar_svd3x3</code> (or your previous assignment’s <code>closest_rotation</code>)</li>
<li><code>igl::min_quad_with_fixed</code></li>
<li><code>igl::cotmatrix_entries</code></li>
<li><code>igl::cotmatrix</code> (or your previous implementation)</li>
<li><code>igl::massmatrix</code> (or your previous implementation)</li>
</ul>
<h3 id="srcbiharmonic_precompute.cpp"><code>src/biharmonic_precompute.cpp</code></h3>
<p>Precompute data needed to efficiently solve for a biharmonic deformation given
a mesh with vertices <code>V</code> and faces <code>F</code> and a list of selected vertices as
indices <code>b</code> into <code>V</code>. The output should be a prefacorized system using the
<code>data</code> struct employed by <code>igl::min_quad_with_fixed</code>.</p>
<h3 id="srcbiharmonic_solve.cpp"><code>src/biharmonic_solve.cpp</code></h3>
<p>Given precomputation <code>data</code> and a list of handle <em>displacements</em> determine
<em>displacements</em> for all vertices in the mesh.</p>
<h3 id="srcarap_precompute.cpp"><code>src/arap_precompute.cpp</code></h3>
<p>Precompute data needed to efficiently conduct local-global iterations for an
arap deformation. This includes the <code>data</code> struct employed by
<code>igl::min_quad_with_fixed</code> to solve the global step" and constructing the
bilinear form <code>K</code> that mixes rotation degrees of freedom with unknown
positions for preparing the covariance matrices of the local step and the
linear term of the global step.</p>
<h3 id="srcarap_single_iteration.cpp"><code>src/arap_single_iteration.cpp</code></h3>
<p>Given precomputed data (<code>data</code> and <code>K</code>), handle <em>positions</em> <code>bc</code> and current
positions of all vertices <code>U</code>, conduct a <em>single</em> iteration of the local-global
solver for minimizing the as-rigid-as-possible energy. Output the <em>positions</em>
of all vertices of the mesh (by overwriting <code>U</code>).</p>
</body>
</html>