-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathblog_decomposition_polar.html
115 lines (83 loc) · 6.34 KB
/
blog_decomposition_polar.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
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="utf-8">
<title>Gabor Makes Games</title>
<meta name="author" content="Gabor Szauer">
<link rel="stylesheet" type="text/css" href="css/shared.css"><link rel="stylesheet" type="text/css" href="css/navigation.css"><link rel="stylesheet" type="text/css" href="css/font-raleway.css"><link rel="stylesheet" type="text/css" href="css/font-oxygen.css"><link rel="stylesheet" type="text/css" href="css/font-worksans.css"><link rel="stylesheet" type="text/css" href="css/codepretty/skins/desert.css"><script type="text/javascript" src="js/codepretty/prettify.js"></script><link rel="stylesheet" type="text/css" href="css/katex.min.css"><script type="text/javascript" src="js/katex.min.js"></script><script type="text/javascript" src="js/katex-autorender.min.js"></script><script type="text/javascript" src="js/navigation.js"></script><!-- Global site tag (gtag.js) - Google Analytics --><script async src="https://www.googletagmanager.com/gtag/js?id=UA-96941899-3"></script><script>window.dataLayer = window.dataLayer || [];function gtag(){dataLayer.push(arguments);}gtag('js', new Date());gtag('config', 'UA-96941899-3');</script> <link rel="stylesheet" type="text/css" href="css/blog.css">
<script type="text/javascript" src="js/blog_decomposition_polar.js"></script>
</head>
<body onload="MainNavOnLoad();renderMathInElement(document.body);PR.prettyPrint();Init();"> <div class="nav"> <ul class="menu"> <li class="logo"><a href="https://gabormakesgames.com">Gabor Makes Games</a></li> <li class="item"><a id="main-nav-active" href="blog.html">Blog</a></li> <li class="item"><a href="books.html">Books</a></li> <li class="item"><a href="https://github.com/gszauer/">Github</a></li> <li class="item"><a href="https://twitter.com/gszauer">@gszauer</a></li> <li class="toggle"><a href="#">Open Menu</a></li> </ul></div>
<div id="blog">
<div id="sidebar"><a class="sidebar-item sidebar-first sidebar-gap" href="blog.html">Back to blog</a></li><a class="sidebar-item sidebar-first" href="blog_decomposition_intro.html">Matrix Decomposition</a><a class="sidebar-item" href="blog_decomposition_code.html">Shoemake's Code</a><a class="sidebar-item" href="blog_decomposition_math.html">The Math</a><a class="sidebar-item" href="blog_decomposition_translation.html">Factoring Translation</a><a class="sidebar-item sidebar-active " href="blog_decomposition_polar.html">Polar Decomposition</a><a class="sidebar-item" href="blog_decomposition_qr.html">QR Decomposition</a><a class="sidebar-item" href="blog_decomposition_spect.html">Spectral Decomposition</a><a class="sidebar-item" href="blog_decomposition_adjust.html">Spectral Adjustment</a><a class="sidebar-item" href="blog_decomposition_affine.html">Affine Decomposition</a><a class="sidebar-item" href="blog_decomposition_summary.html">Summary</a></div>
<div id="content">
<!--?php OpenCodeTag(); ?-->
<!--?php CloseCodeTag(); ?-->
<!--img src="images/" alt=" " /-->
<!--p class="callout"></p-->
<h1>Polar Decomposition</h1>
<p>To decompose \(X\) so that \(X = QS\) we need set \(Q\) equal to \(X\) and then average \(Q\) with its inverse transpose repeatedly. Doing so over and over will eventually converge on a polar rotation, so long as the determinant of \(Q\) was positive. To find \(S\) you just need to multiply \(X\) by the inverse of \(Q\).</p>
<p>To actually implement this, set Q equal to X</p>
<p>$$ Q = X $$</p>
<p>Then, repeatedly average \(Q\) with its inverse transpose</p>
<p>$$ Q_{i + 1} = \frac{1}{2}(Q_{i} + (Q_{i}^{-1})^{T}) $$</p>
<p>Keep doing this until \(Q_{i + 1} - Q_{i} \approx 0\). Alternatley you can do a fixed number of iterations, 20 steps should be more than enough. After you know what \(Q\) is, finding \(S\) is trivial, just multiply the input matrix \(X\) by the inverse of \(Q\):</p>
<p>$$ S = Q^{-1} X $$</p>
<p>Below is an interactive example. The Red, Green and Blue lines represent the right, up and forward basis vectors of a transform matrix. This transform matrix contains both rotation and scale, causing the object to look skewed. The value of the slider in the upper left corner represents the current iteration out of a maximum of 9. Noice how past iteration 5 you can not visually tell that the matrix is being adjusted.</p>
<canvas id="polar_2d_canvas" width="800" height="80">
Canvas support required
</canvas>
<canvas id="polar_3d_canvas" width="800" height="400">
WebGL support required
</canvas>
<p>Understanding how the iterative algorithm works, we can actually implement polar decompositon in code.</p>
<pre class="prettyprint linenums">struct PolarDecompResult {
Matrix Q; // Q is a rotation, or the negative of a rotation
Matrix S; // Scale and skew matrix
}
bool PolarDecompositionEarlyOut(Matrix Q, Matrix Qprev) {
Matrix res = Q - Qprev;
for (var i = 0; i < 16; ++i) {
if (!EpsilonCompare(res[i], 0)) {
return false;
}
}
return true;
}
PolarDecompResult PolarDecomposition(Matrix X) {
PolarDecompResult result;
Matrix Q = X;
for (int i = 0; i < 20; ++i) {
Matrix Qprev = Q;
Matrix Qit = Inverse(Transpose(Q));
Q = MulMatFloat(AddMatMat(Q, Qit), 0.5)
if (PolarDecompositionEarlyOut(Q, Qprev)) {
break;
}
}
result.Q = Q;
result.S = Inverse(Q) * X;
return result;
}</pre>
<p>Remember, \(Q\) can be further factored into a flip and a rotation matrix. If the determinant of \(Q\) is negative, it contains a flip.</p>
<pre class="prettyprint linenums">struct FactorRotationResult {
Matrix F; // Positive or negative identity (flip if negative)
Matrix R; // Rotation matrix
}
FactorRotationResult FactorRotation(Matrix Q) {
FactorRotationResult result;
float f = 1.0;
if (Det(Q) < 0) {
f = -1;
Q = Q * -1; // Component wise multiplication
}
result.F[0] = f;
result.F[5] = f;
result.F[10] = f;
result.R = Q;
return result;
}</pre>
</div>
</div>
</body>
</html>