-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathindex.html
More file actions
383 lines (300 loc) · 20.3 KB
/
Copy pathindex.html
File metadata and controls
383 lines (300 loc) · 20.3 KB
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
<!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>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']], skipTags: ["script","noscript","style","textarea"]}});
</script>
</head>
<body>
<div style="display:none">
<span class="math">\[\newcommand{\mat}[1]{\mathbf{#1}}\]</span>
<span class="math">\[\newcommand{\vec}[1]{\mathbf{#1}}\]</span>
<span class="math">\[\newcommand{\x}{\vec{x}}\]</span>
<span class="math">\[\newcommand{\R}{\mathbb{R}}\]</span>
<span class="math">\[\renewcommand{\v}{\vec{v}}\]</span>
<span class="math">\[\newcommand{\C}{\mat{C}}\]</span>
</div>
<h1 id="geometry_processing_project">Conformal Quaternion Based Mesh (Re)Construction </h1>
<p><strong>To get started:</strong> Clone this repo to your local machine</p>
<h2 id="installation_layout_and_compilation">Installation, Layout, and Compilation</h2>
<ol>
<li>Create a folder called build</li>
<li>Run CMake at the same level as the build folder, with the build folder as the target. Be sure to use an x64 compiler</li>
</ol>
<h2 id="execution">Execution</h2>
<p>Once built, run deformation. It takes up a lot of RAM so be prepared for that. The default is deforming a sphere into a head. To change this you can change the paths in main.cpp</p>
<h2 id="background">Background</h2>
<p>
This project attempts to construct a mesh given a point cloud. To do so we begin with a seed mesh and deform it in such as way that it becomes the mesh implied by the point cloud.
</p>
<img src="./images/front_sphere_to_head_overnight.PNG" style="width:200px; padding-left: 130px; padding-right: 50px"/>
<img src="./images/side_sphere_to_head_overnight.PNG" style="width:200px; padding-right: 50px"/>
<img src="./images/back_sphere_to_head_overnight.PNG" style="width:200px" />
<p style="padding-left:250px"> A sphere deformed to match a point cloud sampling of a human head </p>
<h2 id="spin_transformations">Spin Transformations</h2>
<p>
To do so we must first go over some of the results from the paper <a href="https://www.cs.cmu.edu/~kmcrane/Projects/SpinTransformations/paper.pdf">Spin Transformations of Discrete Surfaces</a> by Keenan et al
</p>
<p>
Keenan et al build on previous results detailing so called "spin transformations" of surfaces. These are transformations of a surface $f$ by a <a href="https://en.wikipedia.org/wiki/Quaternion">quaternion</a> $\lambda$ defined by:
</p>
<p>
$d\tilde{f} = \bar{\lambda} df \lambda$
</p>
<p>
where $df$ and $d\tilde{f}$ are the differential <a href="https://en.wikipedia.org/wiki/One-form#Differential_of_a_function">one forms</a> associated with the orignal surface $f$ and transformed surface $\tilde{f}$.
</p>
<p>
This equation is saying that on a infinitely small scale the relationship between $f$ and $\tilde{f}$ is a <a href="http://run.usc.edu/cs520-s12/quaternions/quaternions-cs520.pdf">rotation and scale</a>.
</p>
<p>
In their paper Keenan et al develop a discrete version of this equation by assigning each vertex $v_{i}$ a quaternion $\lambda_{i}$. The relationship between the original vertices
and the transformed vertices is constructed as a <a href="https://en.wikipedia.org/wiki/Poisson%27s_equation">Poisson problem</a> over the edges of the meshes:
</p>
<p>
$\Delta \tilde{f} = \nabla \cdot \tilde{e}$
</p>
<p>
Discretizing this equation leads to:
</p>
<p>
$\mat{L}\vec{\tilde{f}} = \mat{D}\vec{\tilde{e}}$
</p>
<p>
Where:
<ul>
<li>$\vec{\tilde{e}}$ is a $\mathbb{H}^{|E| \times 1}$ vector of quaternions representing each edge in the transformed mesh</li>
<li>$\vec{\tilde{f}}$ is a $\mathbb{H}^{|V| \times 1}$ vector of quaternions representing each vertex in the transformed mesh</li>
<li>
$\mat{L}$ is the $\mathbb{R}^{|V| \times |V|}$ discrete <a href="http://www.ctralie.com/Teaching/LapMesh/">cotangent laplacian</a>
</li>
<li>
$\mat{D}$ is the $\mathbb{R}^{|V| \times |E|}$ discrete
<a href="https://books.google.ca/books?id=-0T1CAAAQBAJ&pg=PA122&lpg=PA122&dq=cotangent+divergence+operator&source=bl&ots=mI5b5C2Z8R&sig=gED6erLH1TRDn_6
q4dM9KZ-uCAw&hl=en&sa=X&ved=0ahUKEwjqrKnn2ZfTAhXL4IMKHbmaD-IQ6AEIIzAB#v=onepage&q=cotangent%20divergence%20operator&f=false">cotangent divergence</a>
</li>
</ul>
</p>
<p>
This system can be solved in a <a href="https://en.wikipedia.org/wiki/Least_squares">least squares</a> sense for $\tilde{f}$
</p>
<blockquote>
<h5 id="quaternion_vectors">Hold the phone. Quaternion vectors?</h5>
<p>
So you might be wondering how exactly $\tilde{e}$ and $\tilde{f}$ represent edges and vertices. After all each one is just a vector of quaternions. Keenan et al use the following embedding of $\mathbb{R}^{3}$ into $\mathbb{H}$.
Each vector $\vec{p}=(x,y,z)$ in $\mathbb{R}^{3}$ is mapped to a purely imaginary quaternion $(0,x,y,z)$ So the vector $\vec{\tilde{e}} \in \mathbb{H}^{|E| \times 1}$ is represented in $\mathbb{R}$
as an expanded vector $\vec{\tilde{e}}_{\mathbb{R}} \in \mathbb{R}^{4|E| \times 1}$. Similarly for $\vec{\tilde{f}}$.
</p>
<h5 id="quaternion_matrix">That's great, but what about $\mat{L}$ and $\mat{D}$?</h5>
<p>Ah, good question! One way to represent any quaterion $q = (a,b,c,d)$ is with a $4 \times 4$ matrix $Q$</p>
<p>
$\begin{bmatrix}
a & -b & -c & -d\\
b & a & -c & d \\
c & d & a & -b \\
d & -c & -b & a \\
\end{bmatrix}$
</p>
<p>
This representation has many useful properties such as $q_{1}q_{2} = q_{3} \iff Q_{1}Q_{2} = Q_{3}$ and $\bar{q} = Q^{T}$. For our purposes we consider the entries of $\mat{L}$ and $\mat{D}$ as purely real quaternions.
As such each entry gets blown up to a diagonal $4 \times 4$ block in the real versions of $\mat{L}$ and $\mat{D}$
</p>
</blockquote>
<p>
So if we know the transformed edges $\vec{\tilde{e}}$ we can solve for the transformed vertex positions $\vec{\tilde{f}}$. So what's the relationship between $\vec{e}$ and $\vec{\tilde{e}}$?
The details are in the paper so I'll just state the final result
</p>
<p>
$\tilde{e}_{ij} = \frac{1}{3}\bar{\lambda}_{i}e_{ij}\lambda_{i} + \frac{1}{6}\bar{\lambda}_{i}e_{ij}\lambda_{j} + \frac{1}{6}\bar{\lambda}_{j}e_{ij}\lambda_{i} + \frac{1}{3}\bar{\lambda}_{j}e_{ij}\lambda_{j}$
</p>
<p>
Recalling the quaternion multiplication is not commutative. Here $e_{ij}$ is the quaternion representation of the edge vector connecting vertex $i$ to vertex $j$
</p>
<h2 id="goal">Goal</h2>
<p>
Our goal is to find an appropriate vector of quaternions $\vec{\lambda} = (\lambda_{1},\lambda_{2},\ldots,\lambda_{n})$ that will deform the seed mesh into the mesh implied by the point cloud.
To do that we will need some way to measuring how "good" the resuling mesh is. Therefore we will cast this problem as an <a href="https://en.wikipedia.org/wiki/Mathematical_optimization#Optimization_problems">energy minimization problem</a>. The first step in this direction is defining our energy
</p>
<h2 id="energy">Energy</h2>
<p>
There are two easy to implement energies to consider for this kind of problem: point-to-point and point-to-plane. For the following discussion we will be putting aside
the fact that we are dealing with a point cloud sampling of $\tilde{f}$ and instead imagine it as a manifold
</p>
<h3 id="point_to_point">Point to point</h3>
<p>
Point to point energy measures the distance between points on the seed mesh and their closest point on the target mesh (or vice versea).
Denoting by $P_{\tilde{f}}(x)$ the closest point on the target mesh $\tilde{f}$ to a given point $x$ on the seed mesh $f$ we can define the energy as follows:
</p>
<p>
$E = \int\limits_{x \in f}||x - p(x)||^{2}df$
</p>
<p>
Where $p = P_{\tilde{f}}(x)$. Since we will be deforming the seed mesh according to the vector $\vec{\lambda}$ we seek to miniminze this energy with respect to $\lambda$.
In order to make the dependence on $\lambda$ more explicit lets rewrite the energy minimzation problem as
</p>
<p>
$\min\limits_{\vec{\lambda}} E(\vec{\lambda}) = \min\limits_{\vec{\lambda}} \int\limits_{x \in f}||x(\vec{\lambda}) - p(x(\vec{\lambda}))||^{2}df$
</p>
<p>Here $x(\vec{\lambda})$ represents the position of the point $x \in f$ after $f$ is transformed by $\vec{\lambda}$.</p>
<h3 id="point_to_plane">Point to plane</h3>
<p>
Point to plane energy requires some extra information. If in addition to the spacial position of each point $z \in \tilde{f}$ we also have a normal vector $\vec{n}_{z}$
at that point then we can use point to plane energy:
</p>
<p>
$E(\vec{\lambda}) = \int\limits_{x \in f}((x(\vec{\lambda}) - p(x(\vec{\lambda})))\cdot \vec{n}_{p})^{2}df$
</p>
<p>
This energy is measuring the distance between a point $x \in f$ and the plane defined by $p \in \tilde{f}$ and $n_{p}$. It is a looser restriction since the points $x$ and $p$
need only lie in the same plane for their contribution to the energy to be zero
</p>
<h2 id="minimizing_the_energy">Minimizing the energy</h2>
<h3 id="ICP"> Iterative Closest Point </h3>
<p>
We will be using an energy minimization technique called <a href="https://en.wikipedia.org/wiki/Iterative_closest_point">Iterative Closest Point</a> (ICP) whereby we use samples from the meshes involved to approximate $E(\vec{\lambda})$
These are two approches to doing this:
</p>
<ul>
<li>Sample the point cloud and for each point in the sample project it onto the seed mesh to find its closest point</li>
<li>Sample the seed mesh and for each point in the sample loop over the points in the point cloud to find its closest point</li>
</ul>
<p>
We will refer to this as "cloud to mesh" and "mesh to cloud" respectively. In practice it works best to alternate between both methods. In the code from the repo this switch happens if the energy seems to have plateaued between multiple ICP iterations
</p>
<p>Once we have our $k$ pairs $(\vec{x}_{i},\vec{p}_{i})_{i=1,...,k}$ our energies become:</p>
<p>
$E(\vec{\lambda}) \approx \frac{1}{k}\sum\limits_{i=1}^{k}||\vec{x}_{i}(\vec{\lambda}) - \vec{p}_{i}||^{2} \\
E(\vec{\lambda}) \approx \frac{1}{k}\sum\limits_{i=1}^{k}((\vec{x}_{i}(\vec{\lambda}) - \vec{p}_{i}) \cdot \vec{n}_{i})^{2}$
</p>
<p>
Geometrically the first energy is assuming that in the neighbourhood of each point $\vec{p}_{i}$ the target mesh $\tilde{f}$ is approximately a point. Imagine perhaps a large spike with the point at the tip.
This is obviously an unrealistic assumption but it serves as a good base line and works okay in practice.
</p>
<p>
The second energy can be seen as instead assuming that in the neighbourhood of each point $\vec{p}_{i}$ the target mesh $\tilde{f}$ is approximately a plane.
This is much more realistic assumption and tends to work better as long as your points don't "slide" too much as they are moved close to eachother. Meaning you do want $\vec{x}_{i}$ and $\vec{p}_{i}$ to end up
close to each other, not just on the same infinite plane.
</p>
<h3 id="gradient_decent"> Gradient decent </h3>
<p>
Even after making this approxiation we still are faced with a major issue. The relationship between $\vec{x}_{i}$ and $\vec{x}_{i}(\lambda)$ is a complicated one involving solving a least squares problem. As such finding a closed form minimizer
$\tilde{\vec{\lambda}}$ is difficult. Instead we will use <a href="https://en.wikipedia.org/wiki/Gradient_descent">gradient decent</a>. Gradient decent is an iterative algorithm that alternates between 2 steps:
</p>
<ul>
<li>Compute $\vec{\beta} = \nabla E(\vec{\lambda})|_{\lambda = \eta}$</li>
<li>Update $\vec{\eta}$ by $\vec{\eta} \rightarrow \vec{\eta} - \alpha||\vec{\beta}||$</li>
</ul>
<p>
Where $\alpha > 0$ is a parameter that we set to control how far we step in the gradient direction. This works because the gradient always points in the direction of steepest change.
The negative sign in the second step ensures that we will be stepping "downhill" towards a minimum instea of "uphill" towards a maximum.
All that's required to start the algorithm is an inital guess at the minimizer $\vec{\tilde{\lambda}} \approx \vec{\eta}_{0} \in \mathbb{H}^{|V| \times 1}$. For our purposes we use $\vec{\eta}_{0} = (\mathbb{I},\mathbb{I},...,\mathbb{I})$ where $\mathbb{I}$ is
the identity transformation quaternion $\mathbb{I}=(1,0,0,0)$; which represents no rotation or scale.
</p>
<h2 id="bringing_it_together">Bringing it all together</h2>
<p>
We are now in a place to dicuss our algorithm as a whole.
</p>
<code>
1. Define $\vec{\lambda} = (\mathbb{I},\mathbb{I},...,\mathbb{I})$ and $\alpha > 0$ </br>
2. LOOP FOREVER: </br>
3. Select sampling method based on current energy value $E(\vec{\lambda})$</br>
4. Sample seed mesh and point cloud to generate $k$ pairs $(\vec{x}_{i},\vec{p}_{i})$ </br>
5. REPEAT UNTIL CONVERGENCE OR UNTIL MAX ITERATION NUMBER IS REACHED </br>
6. Compute $\vec{x}_{i}(\vec{\lambda})$ for each $\vec{x}_{i}$ </br>
7. Compute $E(\vec{\lambda})$ and $\nabla E(\vec{\lambda})$ </br>
8. $\vec{\lambda} \rightarrow \vec{\lambda} - \alpha||\nabla E(\vec{\lambda})||$ </br>
</code>
<br />
<br />
<img src="./images/trimmed.gif" />
<p style="padding-left:150px">Sped up animation of a deforming a sphere into a human head. The pauses occur when the algorithm is resampling</p>
<h2 id="technical_details">Technical details</h2>
<h3 id="tech_details_sampling">Sampling</h3>
<h4> Mesh to cloud</h4>
<p>
The sampling proceedure for a triangle mesh is described in the readme for <a href="https://github.com/AdamSturge/geometry-processing-registration">registration</a> under the heading "Uniform random sampling of a triangle mesh".
It is repreated here for completeness.
</p>
<blockquote>
<h3 id="uniformrandomsamplingofasingletriangle">Uniform random sampling of a single triangle</h3>
<p>In order to pick a point uniformly randomly in a triangle with corners <span class="math">\(\v_1,
\v_2, \v_3 ∈ \R^3\)</span> we will <em>first</em> pick a point uniformly randomly in the
<a href="https://en.wikipedia.org/wiki/Parallelogram">parallelogram</a> formed by
reflecting <span class="math">\(\v_1\)</span> across the line <span class="math">\(\overline{\v_2\v_3}\)</span>:</p>
<p><span class="math">\[
\x = \v_1 + α (\v_2-\v_1) + β (\v_3 - \v_1)
\]</span></p>
<p>where <span class="math">\(α,β\)</span> are uniformly sampled from the unit interval <span class="math">\([0,1]\)</span>. If <span class="math">\(α+β > 1\)</span>
then the point <span class="math">\(\x\)</span> above will lie in the reflected triangle rather than the
original one. In this case, preprocess <span class="math">\(α\)</span> and <span class="math">\(β\)</span> by setting <span class="math">\(α←1-α\)</span> and
<span class="math">\(β←1-β\)</span> to reflect the point <span class="math">\(\x\)</span> back into the original triangle.</p>
<h3 id="area-weightedrandomsamplingoftriangles">Area-weighted random sampling of triangles</h3>
<p>Assuming we know how to draw a <em>continuous</em> uniform random variable <span class="math">\(γ\)</span> from
the unit interval <span class="math">\([0,1]\)</span>, we would now like to draw a <em>discrete</em> random
triangle index <span class="math">\(T\)</span> from the sequence <span class="math">\({1,‥,m}\)</span> with likelihood proportional to
the relative area of each triangle in the mesh.</p>
<p>We can achieve this by first computing the <a href="https://en.wikipedia.org/wiki/Running_total">cumulative
sum</a> <span class="math">\(\C ∈ \R^{m}\)</span> of the relative
areas:</p>
<p><span class="math">\[
C_i = ∑_{j=1}^m \frac{A_j}{A_X},
\]</span></p>
<p>Then our random index is found by identifying the first entry in <span class="math">\(\C\)</span> whose
value is greater than a uniform random variable <span class="math">\(γ\)</span>. Since <span class="math">\(\C\)</span> is sorted,
locating this entry can be done in <span class="math">\(O(\log m)\)</span>
<a href="https://en.wikipedia.org/wiki/Big_O_notation">time</a>.</p>
</blockquote>
<p>
Once the sampling of the seed mesh is complete it is straightforward to loop over the target point cloud and find the closest point for each sample.
</p>
<h4> Cloud to mesh</h4>
<p>
Given a list of $N$ points we want to sample $k$ points uniformly without replacement. This can be accomplished with the following algorithm:
</p>
<code>
1. FOR i = 0,...,k </br>
2. N = len(point_cloud)</br>
3. j = round(rand(0,1)*(N-1)) % round up if decimal porition is greater than 0.5</br>
4. samples[i] = point_cloud[j]</br>
5. remove point_cloud[j] from point_cloud</br>
</code>
<p>
Once the sampling of the target cloud is complete the closest points on the seed mesh can be found in the following manner: For each point loop over the triangle mesh faces to determine point to face distance.
Once the closest face is found project the point onto the plane defined by the triangle. Once on the plane find the closest point to the triangle by computing
the <a href="https://en.wikipedia.org/wiki/Barycentric_coordinate_system">Barycentric coordinates</a> of the projected point defined by the face vertices.
For more details see <a href="http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.104.4264&rep=rep1&type=pdf">3D Distance from a Point to a Triangle</a>
</p>
<h3 id="tech_details_solving_lease_sqaures">Solving the least squares problem</h3>
<p>
Each iteration of gradient decent requires solving the least squares problem $\mat{L}\vec{\tilde{f}} = \mat{D}\vec{\tilde{e}}$ for $\vec{\tilde{f}}$. As such it is imperative that this system be solved as fast as possible.
To that end a <a href="https://en.wikipedia.org/wiki/Cholesky_decomposition#LDL_decomposition">LDL</a> decomposition is used to solve this system. In practice I've found the solving the system
$\mat{L}^{T}\mat{L}\vec{\tilde{f}} = \mat{L}^{T}\mat{D}\vec{\tilde{e}}$ produces better results. Even with this fast solver each iteration of gradient decent can take between one and two seconds.
</p>
<p>
Once the new vertices $\vec{\tilde{f}}$ have been computed $\vec{x}_{i}(\vec{\lambda})$ can be found by a taking the barycentric coordinates for $\vec{x}_{i}$ and applying them to the appropriate transformed vertices.
</p>
<h3 id="tech_details_gradient_decent">Gradient decent</h3>
<p>
Selecting the ideal $\alpha$ value is a matter of trail and error. One common technique is to start out with a fairly high $\alpha$, say $\alpha=1$, and decrease it every so many iterations.
This way you make large jumps when you are most likely far away from your local minima and become more fine grained as time goes on.
</p>
<p>
In practice I found it better to run a relatiely small number of gradient decent iterations, around $10$, before resampling. This helped jiggle the decent out of any high energy local minimia it might get stuck in.
As a note on this point, since our optimization problem is not convex there is no gaurentee that gradient decent will converge to a global minimum. However a local minimum is better than a random guess.
</p>
<p>
Finding a closed form for the gradient is part of what makes this problem so difficult. However in recent years there has been an explosion of numerical automatic differentiation software packages. In my code I make use of
Stan (<a href="https://arxiv.org/abs/1509.07164">paper</a>, <a href="https://github.com/stan-dev/math">code</a>). Stan uses <a href="https://en.wikipedia.org/wiki/Automatic_differentiation#Reverse_accumulation">reverse mode automatic differentiation</a>
to enable efficent computation of both $E(\vec{\lambda})$ and $\nabla E(\vec{\lambda})$.
</p>
<h3 id="tech_details_ICP">Iterative Closest Point</h3>
<p>
In classic ICP you update your mesh after each iteration. I found it was generally better to always deform the smooth geometry of the initial seed mesh over the potentially jagged geometry of an intermediate mesh.
Therefore my code uses the deformed mesh for resampling and computing closests points, but attempts to find a $\vec{\lambda}$ that will deform the sphere to minimize the energy between those points.
Cast in this light the algorithm defined above is for finding an optimal pairing of points $(\vec{x}_{i},\vec{p}_{i})$ as well as an optimal $\vec{\lambda}$ to bring those points together.
All this means in practice is that you do not reassign your mesh after each step of ICP.
</p>
</body>