| Flutter Impeller |

matrix.cc
Go to the documentation of this file.
1 // Copyright 2013 The Flutter Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style license that can be
3 // found in the LICENSE file.
4
5 #include "impeller/geometry/matrix.h"
6
7 #include <climits>
8 #include <sstream>
9
10 namespace impeller {
11
12 Matrix::Matrix(const MatrixDecomposition& d) : Matrix() {
13 /*
14 * Apply perspective.
15 */
16 for (int i = 0; i < 4; i++) {
17 e[i][3] = d.perspective.e[i];
18 }
19
20 /*
21 * Apply translation.
22 */
23 for (int i = 0; i < 3; i++) {
24 for (int j = 0; j < 3; j++) {
25 e[3][i] += d.translation.e[j] * e[j][i];
26 }
27 }
28
29 /*
30 * Apply rotation.
31 */
32
33 Matrix rotation;
34
35 const auto x = -d.rotation.x;
36 const auto y = -d.rotation.y;
37 const auto z = -d.rotation.z;
38 const auto w = d.rotation.w;
39
40 /*
41 * Construct a composite rotation matrix from the quaternion values.
42 */
43
44 rotation.e[0][0] = 1.0 - 2.0 * (y * y + z * z);
45 rotation.e[0][1] = 2.0 * (x * y - z * w);
46 rotation.e[0][2] = 2.0 * (x * z + y * w);
47 rotation.e[1][0] = 2.0 * (x * y + z * w);
48 rotation.e[1][1] = 1.0 - 2.0 * (x * x + z * z);
49 rotation.e[1][2] = 2.0 * (y * z - x * w);
50 rotation.e[2][0] = 2.0 * (x * z - y * w);
51 rotation.e[2][1] = 2.0 * (y * z + x * w);
52 rotation.e[2][2] = 1.0 - 2.0 * (x * x + y * y);
53
54 *this = *this * rotation;
55
56 /*
57 * Apply shear.
58 */
59 Matrix shear;
60
61 if (d.shear.e[2] != 0) {
62 shear.e[2][1] = d.shear.e[2];
63 *this = *this * shear;
64 }
65
66 if (d.shear.e[1] != 0) {
67 shear.e[2][1] = 0.0;
68 shear.e[2][0] = d.shear.e[1];
69 *this = *this * shear;
70 }
71
72 if (d.shear.e[0] != 0) {
73 shear.e[2][0] = 0.0;
74 shear.e[1][0] = d.shear.e[0];
75 *this = *this * shear;
76 }
77
78 /*
79 * Apply scale.
80 */
81 for (int i = 0; i < 3; i++) {
82 for (int j = 0; j < 3; j++) {
83 e[i][j] *= d.scale.e[i];
84 }
85 }
86 }
87
88 Matrix Matrix::operator+(const Matrix& o) const {
89 return Matrix(
90 m[0] + o.m[0], m[1] + o.m[1], m[2] + o.m[2], m[3] + o.m[3], //
91 m[4] + o.m[4], m[5] + o.m[5], m[6] + o.m[6], m[7] + o.m[7], //
92 m[8] + o.m[8], m[9] + o.m[9], m[10] + o.m[10], m[11] + o.m[11], //
93 m[12] + o.m[12], m[13] + o.m[13], m[14] + o.m[14], m[15] + o.m[15] //
94 );
95 }
96
97 Matrix Matrix::Invert() const {
98 Matrix tmp{
99 m[5] * m[10] * m[15] - m[5] * m[11] * m[14] - m[9] * m[6] * m[15] +
100 m[9] * m[7] * m[14] + m[13] * m[6] * m[11] - m[13] * m[7] * m[10],
101
102 -m[1] * m[10] * m[15] + m[1] * m[11] * m[14] + m[9] * m[2] * m[15] -
103 m[9] * m[3] * m[14] - m[13] * m[2] * m[11] + m[13] * m[3] * m[10],
104
105 m[1] * m[6] * m[15] - m[1] * m[7] * m[14] - m[5] * m[2] * m[15] +
106 m[5] * m[3] * m[14] + m[13] * m[2] * m[7] - m[13] * m[3] * m[6],
107
108 -m[1] * m[6] * m[11] + m[1] * m[7] * m[10] + m[5] * m[2] * m[11] -
109 m[5] * m[3] * m[10] - m[9] * m[2] * m[7] + m[9] * m[3] * m[6],
110
111 -m[4] * m[10] * m[15] + m[4] * m[11] * m[14] + m[8] * m[6] * m[15] -
112 m[8] * m[7] * m[14] - m[12] * m[6] * m[11] + m[12] * m[7] * m[10],
113
114 m[0] * m[10] * m[15] - m[0] * m[11] * m[14] - m[8] * m[2] * m[15] +
115 m[8] * m[3] * m[14] + m[12] * m[2] * m[11] - m[12] * m[3] * m[10],
116
117 -m[0] * m[6] * m[15] + m[0] * m[7] * m[14] + m[4] * m[2] * m[15] -
118 m[4] * m[3] * m[14] - m[12] * m[2] * m[7] + m[12] * m[3] * m[6],
119
120 m[0] * m[6] * m[11] - m[0] * m[7] * m[10] - m[4] * m[2] * m[11] +
121 m[4] * m[3] * m[10] + m[8] * m[2] * m[7] - m[8] * m[3] * m[6],
122
123 m[4] * m[9] * m[15] - m[4] * m[11] * m[13] - m[8] * m[5] * m[15] +
124 m[8] * m[7] * m[13] + m[12] * m[5] * m[11] - m[12] * m[7] * m[9],
125
126 -m[0] * m[9] * m[15] + m[0] * m[11] * m[13] + m[8] * m[1] * m[15] -
127 m[8] * m[3] * m[13] - m[12] * m[1] * m[11] + m[12] * m[3] * m[9],
128
129 m[0] * m[5] * m[15] - m[0] * m[7] * m[13] - m[4] * m[1] * m[15] +
130 m[4] * m[3] * m[13] + m[12] * m[1] * m[7] - m[12] * m[3] * m[5],
131
132 -m[0] * m[5] * m[11] + m[0] * m[7] * m[9] + m[4] * m[1] * m[11] -
133 m[4] * m[3] * m[9] - m[8] * m[1] * m[7] + m[8] * m[3] * m[5],
134
135 -m[4] * m[9] * m[14] + m[4] * m[10] * m[13] + m[8] * m[5] * m[14] -
136 m[8] * m[6] * m[13] - m[12] * m[5] * m[10] + m[12] * m[6] * m[9],
137
138 m[0] * m[9] * m[14] - m[0] * m[10] * m[13] - m[8] * m[1] * m[14] +
139 m[8] * m[2] * m[13] + m[12] * m[1] * m[10] - m[12] * m[2] * m[9],
140
141 -m[0] * m[5] * m[14] + m[0] * m[6] * m[13] + m[4] * m[1] * m[14] -
142 m[4] * m[2] * m[13] - m[12] * m[1] * m[6] + m[12] * m[2] * m[5],
143
144 m[0] * m[5] * m[10] - m[0] * m[6] * m[9] - m[4] * m[1] * m[10] +
145 m[4] * m[2] * m[9] + m[8] * m[1] * m[6] - m[8] * m[2] * m[5]};
146
147 Scalar det =
148 m[0] * tmp.m[0] + m[1] * tmp.m[4] + m[2] * tmp.m[8] + m[3] * tmp.m[12];
149
150 if (det == 0) {
151 return {};
152 }
153
154 det = 1.0 / det;
155
156 return {tmp.m[0] * det, tmp.m[1] * det, tmp.m[2] * det, tmp.m[3] * det,
157 tmp.m[4] * det, tmp.m[5] * det, tmp.m[6] * det, tmp.m[7] * det,
158 tmp.m[8] * det, tmp.m[9] * det, tmp.m[10] * det, tmp.m[11] * det,
159 tmp.m[12] * det, tmp.m[13] * det, tmp.m[14] * det, tmp.m[15] * det};
160 }
161
162 Scalar Matrix::GetDeterminant() const {
163 auto a00 = e[0][0];
164 auto a01 = e[0][1];
165 auto a02 = e[0][2];
166 auto a03 = e[0][3];
167 auto a10 = e[1][0];
168 auto a11 = e[1][1];
169 auto a12 = e[1][2];
170 auto a13 = e[1][3];
171 auto a20 = e[2][0];
172 auto a21 = e[2][1];
173 auto a22 = e[2][2];
174 auto a23 = e[2][3];
175 auto a30 = e[3][0];
176 auto a31 = e[3][1];
177 auto a32 = e[3][2];
178 auto a33 = e[3][3];
179
180 auto b00 = a00 * a11 - a01 * a10;
181 auto b01 = a00 * a12 - a02 * a10;
182 auto b02 = a00 * a13 - a03 * a10;
183 auto b03 = a01 * a12 - a02 * a11;
184 auto b04 = a01 * a13 - a03 * a11;
185 auto b05 = a02 * a13 - a03 * a12;
186 auto b06 = a20 * a31 - a21 * a30;
187 auto b07 = a20 * a32 - a22 * a30;
188 auto b08 = a20 * a33 - a23 * a30;
189 auto b09 = a21 * a32 - a22 * a31;
190 auto b10 = a21 * a33 - a23 * a31;
191 auto b11 = a22 * a33 - a23 * a32;
192
193 return b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
194 }
195
196 Scalar Matrix::GetMaxBasisLength() const {
197 Scalar max = 0;
198 for (int i = 0; i < 3; i++) {
199 max = std::max(max,
200 e[i][0] * e[i][0] + e[i][1] * e[i][1] + e[i][2] * e[i][2]);
201 }
202 return std::sqrt(max);
203 }
204
205 /*
206 * Adapted for Impeller from Graphics Gems:
207 * http://www.realtimerendering.com/resources/GraphicsGems/gemsii/unmatrix.c
208 */
209 std::optional<MatrixDecomposition> Matrix::Decompose() const {
210 /*
211 * Normalize the matrix.
212 */
213 Matrix self = *this;
214
215 if (self.e[3][3] == 0) {
216 return std::nullopt;
217 }
218
219 for (int i = 0; i < 4; i++) {
220 for (int j = 0; j < 4; j++) {
221 self.e[i][j] /= self.e[3][3];
222 }
223 }
224
225 /*
226 * `perspectiveMatrix` is used to solve for perspective, but it also provides
227 * an easy way to test for singularity of the upper 3x3 component.
228 */
229 Matrix perpectiveMatrix = self;
230 for (int i = 0; i < 3; i++) {
231 perpectiveMatrix.e[i][3] = 0;
232 }
233
234 perpectiveMatrix.e[3][3] = 1;
235
236 if (perpectiveMatrix.GetDeterminant() == 0.0) {
237 return std::nullopt;
238 }
239
240 MatrixDecomposition result;
241
242 /*
243 * ==========================================================================
244 * First, isolate perspective.
245 * ==========================================================================
246 */
247 if (self.e[0][3] != 0.0 || self.e[1][3] != 0.0 || self.e[2][3] != 0.0) {
248 /*
249 * prhs is the right hand side of the equation.
250 */
251 const Vector4 rightHandSide(self.e[0][3], //
252 self.e[1][3], //
253 self.e[2][3], //
254 self.e[3][3]);
255
256 /*
257 * Solve the equation by inverting `perspectiveMatrix` and multiplying
258 * prhs by the inverse.
259 */
260
261 result.perspective = perpectiveMatrix.Invert().Transpose() * rightHandSide;
262
263 /*
264 * Clear the perspective partition.
265 */
266 self.e[0][3] = self.e[1][3] = self.e[2][3] = 0;
267 self.e[3][3] = 1;
268 }
269
270 /*
271 * ==========================================================================
272 * Next, the translation.
273 * ==========================================================================
274 */
275 result.translation = {self.e[3][0], self.e[3][1], self.e[3][2]};
276 self.e[3][0] = self.e[3][1] = self.e[3][2] = 0.0;
277
278 /*
279 * ==========================================================================
280 * Next, the scale and shear.
281 * ==========================================================================
282 */
283 Vector3 row[3];
284 for (int i = 0; i < 3; i++) {
285 row[i].x = self.e[i][0];
286 row[i].y = self.e[i][1];
287 row[i].z = self.e[i][2];
288 }
289
290 /*
291 * Compute X scale factor and normalize first row.
292 */
293 result.scale.x = row[0].Length();
294 row[0] = row[0].Normalize();
295
296 /*
297 * Compute XY shear factor and make 2nd row orthogonal to 1st.
298 */
299 result.shear.xy = row[0].Dot(row[1]);
300 row[1] = Vector3::Combine(row[1], 1.0, row[0], -result.shear.xy);
301
302 /*
303 * Compute Y scale and normalize 2nd row.
304 */
305 result.scale.y = row[1].Length();
306 row[1] = row[1].Normalize();
307 result.shear.xy /= result.scale.y;
308
309 /*
310 * Compute XZ and YZ shears, orthogonalize 3rd row.
311 */
312 result.shear.xz = row[0].Dot(row[2]);
313 row[2] = Vector3::Combine(row[2], 1.0, row[0], -result.shear.xz);
314 result.shear.yz = row[1].Dot(row[2]);
315 row[2] = Vector3::Combine(row[2], 1.0, row[1], -result.shear.yz);
316
317 /*
318 * Next, get Z scale and normalize 3rd row.
319 */
320 result.scale.z = row[2].Length();
321 row[2] = row[2].Normalize();
322
323 result.shear.xz /= result.scale.z;
324 result.shear.yz /= result.scale.z;
325
326 /*
327 * At this point, the matrix (in rows[]) is orthonormal.
328 * Check for a coordinate system flip. If the determinant
329 * is -1, then negate the matrix and the scaling factors.
330 */
331 if (row[0].Dot(row[1].Cross(row[2])) < 0) {
332 result.scale.x *= -1;
333 result.scale.y *= -1;
334 result.scale.z *= -1;
335
336 for (int i = 0; i < 3; i++) {
337 row[i].x *= -1;
338 row[i].y *= -1;
339 row[i].z *= -1;
340 }
341 }
342
343 /*
344 * ==========================================================================
345 * Finally, get the rotations out.
346 * ==========================================================================
347 */
348 result.rotation.x =
349 0.5 * sqrt(fmax(1.0 + row[0].x - row[1].y - row[2].z, 0.0));
350 result.rotation.y =
351 0.5 * sqrt(fmax(1.0 - row[0].x + row[1].y - row[2].z, 0.0));
352 result.rotation.z =
353 0.5 * sqrt(fmax(1.0 - row[0].x - row[1].y + row[2].z, 0.0));
354 result.rotation.w =
355 0.5 * sqrt(fmax(1.0 + row[0].x + row[1].y + row[2].z, 0.0));
356
357 if (row[2].y > row[1].z) {
358 result.rotation.x = -result.rotation.x;
359 }
360 if (row[0].z > row[2].x) {
361 result.rotation.y = -result.rotation.y;
362 }
363 if (row[1].x > row[0].y) {
364 result.rotation.z = -result.rotation.z;
365 }
366
367 return result;
368 }
369
370 uint64_t MatrixDecomposition::GetComponentsMask() const {
371 uint64_t mask = 0;
372
373 Quaternion noRotation(0.0, 0.0, 0.0, 1.0);
374 if (rotation != noRotation) {
375 mask = mask | static_cast<uint64_t>(Component::kRotation);
376 }
377
378 Vector4 defaultPerspective(0.0, 0.0, 0.0, 1.0);
379 if (perspective != defaultPerspective) {
380 mask = mask | static_cast<uint64_t>(Component::kPerspective);
381 }
382
383 Shear noShear(0.0, 0.0, 0.0);
384 if (shear != noShear) {
385 mask = mask | static_cast<uint64_t>(Component::kShear);
386 }
387
388 Vector3 defaultScale(1.0, 1.0, 1.0);
389 if (scale != defaultScale) {
390 mask = mask | static_cast<uint64_t>(Component::kScale);
391 }
392
393 Vector3 defaultTranslation(0.0, 0.0, 0.0);
394 if (translation != defaultTranslation) {
395 mask = mask | static_cast<uint64_t>(Component::kTranslation);
396 }
397
398 return mask;
399 }
400
401 } // namespace impeller
impeller::Matrix::operator+
Matrix operator+(const Vector3 &t) const
Definition: matrix.h:394
impeller::Matrix::m
Scalar m[16]
Definition: matrix.h:39
impeller::Vector3::Dot
constexpr Scalar Dot(const Vector3 &other) const
Definition: vector.h:54
impeller::Matrix::Decompose
std::optional< MatrixDecomposition > Decompose() const
Definition: matrix.cc:209
impeller::Matrix::Matrix
constexpr Matrix()
Definition: matrix.h:47
impeller::Quaternion::z
Scalar z
Definition: quaternion.h:19
impeller::Scalar
float Scalar
Definition: scalar.h:18
impeller::Vector3::Combine
static constexpr Vector3 Combine(const Vector3 &a, Scalar aScale, const Vector3 &b, Scalar bScale)
Definition: vector.h:192
impeller::Quaternion::w
Scalar w
Definition: quaternion.h:20
impeller::MatrixDecomposition::Component::kPerspective
@ kPerspective
impeller::Vector4
Definition: vector.h:232
impeller::Matrix::GetDeterminant
Scalar GetDeterminant() const
Definition: matrix.cc:162
impeller::MatrixDecomposition::shear
Shear shear
Definition: matrix_decomposition.h:18
impeller::Vector3::x
Scalar x
Definition: vector.h:23
impeller::Quaternion::x
Scalar x
Definition: quaternion.h:17
impeller::Matrix::e
Scalar e[4][4]
Definition: matrix.h:40
impeller::MatrixDecomposition
Definition: matrix_decomposition.h:15
impeller::MatrixDecomposition::Component::kScale
@ kScale
matrix.h
impeller::Quaternion
Definition: quaternion.h:14
impeller::MatrixDecomposition::Component::kTranslation
@ kTranslation
impeller::MatrixDecomposition::perspective
Vector4 perspective
Definition: matrix_decomposition.h:19
impeller::Matrix::Transpose
constexpr Matrix Transpose() const
Definition: matrix.h:278
impeller::Matrix::GetMaxBasisLength
Scalar GetMaxBasisLength() const
Definition: matrix.cc:196
impeller::Vector3::z
Scalar z
Definition: vector.h:25
impeller::MatrixDecomposition::translation
Vector3 translation
Definition: matrix_decomposition.h:16
impeller::Vector3::Length
constexpr Scalar Length() const
Definition: vector.h:47
impeller::Vector3::e
Scalar e[3]
Definition: vector.h:27
impeller::Vector3::y
Scalar y
Definition: vector.h:24
impeller::Shear::yz
double yz
Definition: shear.h:17
impeller::MatrixDecomposition::Component::kShear
@ kShear
impeller::MatrixDecomposition::rotation
Quaternion rotation
Definition: matrix_decomposition.h:20
impeller::Matrix::Invert
Matrix Invert() const
Definition: matrix.cc:97
impeller::Vector3::Normalize
constexpr Vector3 Normalize() const
Definition: vector.h:49
impeller::MatrixDecomposition::scale
Vector3 scale
Definition: matrix_decomposition.h:17
impeller::MatrixDecomposition::Component::kRotation
@ kRotation
impeller::Shear::e
double e[3]
Definition: shear.h:19
impeller::Vector4::e
Scalar e[4]
Definition: vector.h:240
impeller::MatrixDecomposition::GetComponentsMask
uint64_t GetComponentsMask() const
Definition: matrix.cc:370
impeller::Shear::xz
double xz
Definition: shear.h:16
impeller::Quaternion::y
Scalar y
Definition: quaternion.h:18
impeller::Shear::xy
double xy
Definition: shear.h:15
impeller::Shear
Definition: shear.h:12
impeller
Definition: aiks_blur_unittests.cc:20
impeller::Matrix
A 4x4 matrix using column-major storage.
Definition: matrix.h:37
impeller::Vector3
Definition: vector.h:20