-
Notifications
You must be signed in to change notification settings - Fork 3
Expand file tree
/
Copy pathlinear-transformations-numpy.html
More file actions
313 lines (290 loc) · 18.7 KB
/
linear-transformations-numpy.html
File metadata and controls
313 lines (290 loc) · 18.7 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
<!DOCTYPE html>
<html lang="en" itemscope itemtype="http://schema.org/Article">
<head>
<title>Linear transformations in Numpy</title>
<meta charset="utf-8">
<meta property="og:title" content="Linear transformations in Numpy">
<meta property="og:site_name" content="Modesto Mas | Blog">
<meta property="og:image" content="https://mmas.github.io/images/profile.jpg">
<meta property="og:image:width" content="200">
<meta property="og:image:height" content="200">
<meta property="og:url" content="https://mmas.github.io/linear-transformations-numpy">
<meta property="og:locale" content="en_GB">
<meta name="twitter:image" content="https://mmas.github.io/images/profile.jpg">
<meta name="twitter:url" content="https://mmas.github.io/linear-transformations-numpy">
<meta name="twitter:card" content="summary">
<meta name="twitter:domain" content="mmas.github.io">
<meta name="twitter:title" content="Linear transformations in Numpy">
<meta name="description" content="A linear transformation of the plane \(\mathbb R^2\) is a geometric transformation of the form \[ f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}...">
<meta name="twitter:description" content="A linear transformation of the plane \(\mathbb R^2\) is a geometric transformation of the form \[ f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}...">
<meta property="og:description" content="A linear transformation of the plane \(\mathbb R^2\) is a geometric transformation of the form \[ f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}...">
<meta name="keywords" content="geometric-transformations,geometry,matplotlib,numpy,python">
<meta property="og:type" content="blog">
<meta name="viewport" content="width=device-width, initial-scale=1">
<meta property="og:type" content="article">
<meta property="article:author" content="https://github.com/mmas">
<meta property="article:section" content="geometric-transformations">
<meta property="article:tag" content="geometric-transformations,geometry,matplotlib,numpy,python">
<meta property="article:published_time" content="2016-06-11">
<meta property="article:modified_time" content="2016-06-11">
<link rel="stylesheet" type="text/css" href="/css/main.css">
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
CommonHTML: {
scale: 93,
showMathMenu: false
},
tex2jax: {
"inlineMath": [["$","$"], ["\\(","\\)"]]
}
});
</script>
<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-MML-AM_CHTML"></script>
</head>
<body class="entry-detail">
<header>
<div>
<img src="https://mmas.github.io/images/profile.jpg">
<a class="brand" href="/">Modesto Mas</a>
<span>Data/Python/DevOps Engineer</span>
<nav>
<ul>
<li><a href="/tags">Tags</a></li>
<li><a href="https://github.com/mmas/mmas.github.io/issues" target="_blank">Issues</a></li>
</ul>
</nav>
</div>
</header>
<section id="content" role="main">
<article>
<header>
<h1><a href="/linear-transformations-numpy">Linear transformations in Numpy</a></h1>
<time datetime="2016-06-11">Jun 11, 2016</time>
<a class="tag" href="/tags?tag=geometric-transformations">geometric-transformations</a>
<a class="tag" href="/tags?tag=geometry">geometry</a>
<a class="tag" href="/tags?tag=matplotlib">matplotlib</a>
<a class="tag" href="/tags?tag=numpy">numpy</a>
<a class="tag" href="/tags?tag=python">python</a>
</header>
<aside id="article-nav"></aside>
<section class="body">
<p>A linear transformation of the plane <span class="math">\(\mathbb R^2\)</span> is a geometric transformation of the form</p>
<p><span class="math">\[
f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}a&b\\c&d\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},
\]</span></p>
<p>where <span class="math">\(a\)</span>, <span class="math">\(b\)</span>, <span class="math">\(c\)</span> and <span class="math">\(d\)</span> are real constants.</p>
<p>Linear transformations leave the origin fixed and preserve parallelism. Scaling, shearing, rotation and reflexion of a plane are examples of linear transformations.</p>
<p>Applying a geometric transformation to a given matrix in Numpy requires applying the inverse of the transformation to the coordinates of the matrix, create a new matrix of indices from the coordinates and map the matrix to the new indices. Since this can be tricky, let's start with a simple example involving a matrix that represents the indices itself.</p>
<h2>A simple example. Indices transformation</h2>
<p>Start with a simple 3x4 matrix:</p>
<pre class="sourceCode python"><code class="sourceCode python">>>> <span class="ch">import</span> numpy <span class="ch">as</span> np
>>> M, N = <span class="dv">3</span>, <span class="dv">4</span>
>>> matrix = np.arange(M*N).reshape((M,N))
>>> matrix
array([[ <span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">2</span>, <span class="dv">3</span>],
[ <span class="dv">4</span>, <span class="dv">5</span>, <span class="dv">6</span>, <span class="dv">7</span>],
[ <span class="dv">8</span>, <span class="dv">9</span>, <span class="dv">10</span>, <span class="dv">11</span>]])</code></pre>
<p>Then, we need to obtain the indices pairs of the matrix in a matrix form. The new indices of the matrix will result from the product of the inverse of the transformation matrix and this matrix, therefore the indices pairs in this case need to be a 2x12 matrix as</p>
<p><span class="math">\[
\mathbf p =
\begin{pmatrix}
0 & 0 & 0 & 1 & 1 & 1 & 2 & 2 & 2 & 3 & 3 & 3 \\
0 & 1 & 2 & 0 & 1 & 2 & 0 & 1 & 2 & 0 & 1 & 2
\end{pmatrix},
\]</span></p>
<p>then:</p>
<pre class="sourceCode python"><code class="sourceCode python">>>> points = np.mgrid[<span class="dv">0</span>:N, <span class="dv">0</span>:M].reshape((<span class="dv">2</span>, M*N))
>>> points
array([[<span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">1</span>, <span class="dv">1</span>, <span class="dv">2</span>, <span class="dv">2</span>, <span class="dv">2</span>, <span class="dv">3</span>, <span class="dv">3</span>, <span class="dv">3</span>],
[<span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">2</span>, <span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">2</span>, <span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">2</span>, <span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">2</span>]])</code></pre>
<p>An alternative way to get the indices pairs, longer but maybe more obvious:</p>
<pre class="sourceCode python"><code class="sourceCode python">>>> x, y = np.mgrid[<span class="dv">0</span>:N, <span class="dv">0</span>:M]
>>> points = np.vstack([x.ravel(), y.ravel()])</code></pre>
<p>Now, apply the transformation to the indices pairs. The new indices pairs need to be integers to map the given matrix to the indices. In this example the result will be just casted to integers for simplicity, then it will be easy to spot the result of the transformation since it involve points halfway between two integers. Other methods require interpolations of the given matrix from the indices. As an example, the transformation matrix will be</p>
<p><span class="math">\[
\mathbf A = \begin{pmatrix}2&0\\0&1\end{pmatrix},
\]</span></p>
<p>which corresponds with scaling the plane in the <span class="math">\(x\)</span>-axis by a factor of <span class="math">\(2\)</span>. Hence, the new indices pairs are</p>
<p><span class="math">\[
\mathbf p' = \lfloor \mathbf A \mathbf p \rfloor =
\begin{pmatrix}
0&0&0&0&0&0&1&1&1&1&1&1\\
0&1&2&0&1&2&0&1&2&0&1&2
\end{pmatrix},
\]</span></p>
<p>then:</p>
<pre class="sourceCode python"><code class="sourceCode python">>>> a = np.array([[<span class="dv">2</span>, <span class="dv">0</span>],
[<span class="dv">0</span>, <span class="dv">1</span>]])
>>> new_points = np.linalg.inv(a).dot(points).astype(<span class="dt">int</span>)
>>> new_points
array([[<span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">1</span>, <span class="dv">1</span>, <span class="dv">1</span>, <span class="dv">1</span>, <span class="dv">1</span>],
[<span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">2</span>, <span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">2</span>, <span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">2</span>, <span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">2</span>]])</code></pre>
<p>To finish this example, convert the indices pairs to a matrix of indices which in this example corresponds with the resulting matrix.</p>
<p>The <span class="math">\(x\)</span> components is</p>
<p><span class="math">\[
\mathbf x =
\begin{pmatrix}
0 & 0 & 1 & 1 \\
0 & 0 & 1 & 1 \\
0 & 0 & 1 & 1
\end{pmatrix},
\]</span></p>
<p>the <span class="math">\(y\)</span> component is</p>
<p><span class="math">\[
\mathbf y =
\begin{pmatrix}
0 & 0 & 0 & 0 \\
1 & 1 & 1 & 1 \\
2 & 2 & 2 & 2
\end{pmatrix},
\]</span></p>
<p>and the resulting matrix from the <span class="math">\(x\)</span> and <span class="math">\(y\)</span> components is</p>
<p><span class="math">\[
\mathbf x + N \mathbf y =
\begin{pmatrix}
0 & 0 & 1 & 1 \\
4 & 4 & 5 & 5 \\
8 & 8 & 9 & 9
\end{pmatrix}.
\]</span></p>
<p>Here it is easy to spot the scaling of the values along the <span class="math">\(x\)</span>-axis in the matrix. These final steps correspond to:</p>
<pre class="sourceCode python"><code class="sourceCode python">>>> x, y = new_points.reshape((<span class="dv">2</span>, M, N), order=<span class="st">'F'</span>)
>>> x + N*y
array([[<span class="dv">0</span>, <span class="dv">0</span>, <span class="dv">1</span>, <span class="dv">1</span>],
[<span class="dv">4</span>, <span class="dv">4</span>, <span class="dv">5</span>, <span class="dv">5</span>],
[<span class="dv">8</span>, <span class="dv">8</span>, <span class="dv">9</span>, <span class="dv">9</span>]])</code></pre>
<h2>A more visual example. Matrix transformation</h2>
<p>In the following example we will use a bigger matrix, represented as an image for visual support. Once we calculate the new indices matrix we will map the original matrix to the new indices, wrapping the out-of-bounds indices to obtain a continuous plane using <a target="_blank" href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.take.html">numpy.take</a> with <code>mode='wrap'</code>.</p>
<pre class="sourceCode python"><code class="sourceCode python"><span class="ch">import</span> matplotlib <span class="ch">as</span> mpl
<span class="ch">import</span> matplotlib.pyplot <span class="ch">as</span> plt</code></pre>
<p>Some visual settings:</p>
<pre class="sourceCode python"><code class="sourceCode python">mpl.rcParams.update({<span class="st">'image.cmap'</span>: <span class="st">'Accent'</span>,
<span class="st">'image.interpolation'</span>: <span class="st">'none'</span>,
<span class="st">'xtick.major.width'</span>: <span class="dv">0</span>,
<span class="st">'xtick.labelsize'</span>: <span class="dv">0</span>,
<span class="st">'ytick.major.width'</span>: <span class="dv">0</span>,
<span class="st">'ytick.labelsize'</span>: <span class="dv">0</span>,
<span class="st">'axes.linewidth'</span>: <span class="dv">0</span>})</code></pre>
<p>Create a 200x200 matrix for this example:</p>
<pre class="sourceCode python"><code class="sourceCode python">aux = np.ones((<span class="dv">100</span>, <span class="dv">100</span>), dtype=<span class="dt">int</span>)
src = np.vstack([np.c_[aux, <span class="dv">2</span>*aux], np.c_[<span class="dv">3</span>*aux, <span class="dv">4</span>*aux]])
plt.imshow(src)
plt.show()</code></pre>
<div class="figure">
<img style="max-width:600px" alt="Linear transformation example in Numpy" src="https://mmas.github.io/images/figure_1_linear_transformation_numpy.png" />
</div>
<p>The linear transformation function, which includes the operations of the previous examples but rounding the new indices pairs and mapping the source matrix to the new indices might be written as follows:</p>
<pre class="sourceCode python"><code class="sourceCode python"><span class="kw">def</span> linear_transformation(src, a):
M, N = src.shape
points = np.mgrid[<span class="dv">0</span>:N, <span class="dv">0</span>:M].reshape((<span class="dv">2</span>, M*N))
new_points = np.linalg.inv(a).dot(points).<span class="dt">round</span>().astype(<span class="dt">int</span>)
x, y = new_points.reshape((<span class="dv">2</span>, M, N), order=<span class="st">'F'</span>)
indices = x + N*y
<span class="kw">return</span> np.take(src, indices, mode=<span class="st">'wrap'</span>)</code></pre>
<p>Let's see some linear transformations we can do.</p>
<p>Scaling the plane in the <span class="math">\(x\)</span>-axis by a factor of 1.5:</p>
<p><span class="math">\[
f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}1.5 & 0\\0 & 0\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},
\]</span></p>
<p>in Numpy:</p>
<pre class="sourceCode python"><code class="sourceCode python">a = np.array([[<span class="fl">1.5</span>, <span class="dv">0</span>],
[<span class="dv">0</span>, <span class="dv">1</span>]])
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()</code></pre>
<div class="figure">
<img style="max-width:600px" alt="Linear transformation example in Numpy" src="https://mmas.github.io/images/figure_2_linear_transformation_numpy.png" />
</div>
<p>Dilating the plane by a factor of 1.8:</p>
<p><span class="math">\[
f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}1.8 & 0\\0 & 1.8\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},
\]</span></p>
<p>in Numpy:</p>
<pre class="sourceCode python"><code class="sourceCode python">a = <span class="fl">1.8</span>*np.eye(<span class="dv">2</span>)
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()</code></pre>
<div class="figure">
<img style="max-width:600px" alt="Linear transformation example in Numpy" src="https://mmas.github.io/images/figure_3_linear_transformation_numpy.png" />
</div>
<p>Dilating the plane by a factor of 0.5:</p>
<p><span class="math">\[
f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}0.5 & 0\\0 & 0.5\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},
\]</span></p>
<p>in Numpy:</p>
<pre class="sourceCode python"><code class="sourceCode python">a = .<span class="dv">5</span>*np.eye(<span class="dv">2</span>)
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()</code></pre>
<div class="figure">
<img style="max-width:600px" alt="Linear transformation example in Numpy" src="https://mmas.github.io/images/figure_4_linear_transformation_numpy.png" />
</div>
<p>Scaling the plane in the <span class="math">\(y\)</span>-axis by a factor of 0.5:</p>
<p><span class="math">\[
f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}1 & 0\\0 & 0.5\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},
\]</span></p>
<p>in Numpy:</p>
<pre class="sourceCode python"><code class="sourceCode python">a = np.array([[<span class="dv">1</span>, <span class="dv">0</span>],
[<span class="dv">0</span>, .<span class="dv">5</span>]])
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()</code></pre>
<div class="figure">
<img style="max-width:600px" alt="Linear transformation example in Numpy" src="https://mmas.github.io/images/figure_5_linear_transformation_numpy.png" />
</div>
<p>Shearing about the <span class="math">\(y\)</span>-axis with a vertical displacement of <span class="math">\(+x/2\)</span>:</p>
<p><span class="math">\[
f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}1 & 0\\{1 \over 2} & 0\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},
\]</span></p>
<p>in Numpy:</p>
<pre class="sourceCode python"><code class="sourceCode python">a = np.array([[<span class="dv">1</span>, <span class="dv">0</span>],
[.<span class="dv">5</span>, <span class="dv">1</span>]])
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()</code></pre>
<div class="figure">
<img style="max-width:600px" alt="Linear transformation example in Numpy" src="https://mmas.github.io/images/figure_6_linear_transformation_numpy.png" />
</div>
<p>Rotation through <span class="math">\(45^{\circ}\)</span> about the origin:</p>
<p><span class="math">\[
f \begin{pmatrix}x\\y\end{pmatrix} =
\begin{pmatrix}
\cos{\pi \over 4} & -\sin{\pi \over 4} \\
\sin{\pi \over 4} & \cos{\pi \over 4}
\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},
\]</span></p>
<p>in Numpy:</p>
<pre class="sourceCode python"><code class="sourceCode python">alpha = np.pi/<span class="dv">4</span>
a = np.array([[np.cos(alpha), -np.sin(alpha)],
[np.sin(alpha), np.cos(alpha)]])
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()</code></pre>
<div class="figure">
<img style="max-width:600px" alt="Linear transformation example in Numpy" src="https://mmas.github.io/images/figure_7_linear_transformation_numpy.png" />
</div>
<p>Reflection in a line with inclination of <span class="math">\(45^{\circ}\)</span> through the origin:</p>
<p><span class="math">\[
f \begin{pmatrix}x\\y\end{pmatrix} =
\begin{pmatrix}
\cos{\pi \over 2} & \sin{\pi \over 2} \\
\sin{\pi \over 2} & -\cos{\pi \over 2}
\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},
\]</span></p>
<p>in Numpy:</p>
<pre class="sourceCode python"><code class="sourceCode python">alpha = np.pi/<span class="dv">4</span>
a = np.array([[np.cos(<span class="dv">2</span>*alpha), np.sin(<span class="dv">2</span>*alpha)],
[np.sin(<span class="dv">2</span>*alpha), -np.cos(<span class="dv">2</span>*alpha)]])
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()</code></pre>
<div class="figure">
<img style="max-width:600px" alt="Linear transformation example in Numpy" src="https://mmas.github.io/images/figure_8_linear_transformation_numpy.png" />
</div>
</section>
</article>
</section>
<footer></footer>
<script type="text/javascript" src="/js/article.js"></script>
</body>
</html>