| 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 #include "flutter/fml/logging.h"
11
12 namespace impeller {
13
14 Matrix::Matrix(const MatrixDecomposition& d) : Matrix() {
15 /*
16 * Apply perspective.
17 */
18 for (int i = 0; i < 4; i++) {
19 e[i][3] = d.perspective.e[i];
20 }
21
22 /*
23 * Apply translation.
24 */
25 for (int i = 0; i < 3; i++) {
26 for (int j = 0; j < 3; j++) {
27 e[3][i] += d.translation.e[j] * e[j][i];
28 }
29 }
30
31 /*
32 * Apply rotation.
33 */
34
35 Matrix rotation;
36
37 const auto x = -d.rotation.x;
38 const auto y = -d.rotation.y;
39 const auto z = -d.rotation.z;
40 const auto w = d.rotation.w;
41
42 /*
43 * Construct a composite rotation matrix from the quaternion values.
44 */
45
46 rotation.e[0][0] = 1.0 - 2.0 * (y * y + z * z);
47 rotation.e[0][1] = 2.0 * (x * y - z * w);
48 rotation.e[0][2] = 2.0 * (x * z + y * w);
49 rotation.e[1][0] = 2.0 * (x * y + z * w);
50 rotation.e[1][1] = 1.0 - 2.0 * (x * x + z * z);
51 rotation.e[1][2] = 2.0 * (y * z - x * w);
52 rotation.e[2][0] = 2.0 * (x * z - y * w);
53 rotation.e[2][1] = 2.0 * (y * z + x * w);
54 rotation.e[2][2] = 1.0 - 2.0 * (x * x + y * y);
55
56 *this = *this * rotation;
57
58 /*
59 * Apply shear.
60 */
61 Matrix shear;
62
63 if (d.shear.e[2] != 0) {
64 shear.e[2][1] = d.shear.e[2];
65 *this = *this * shear;
66 }
67
68 if (d.shear.e[1] != 0) {
69 shear.e[2][1] = 0.0;
70 shear.e[2][0] = d.shear.e[1];
71 *this = *this * shear;
72 }
73
74 if (d.shear.e[0] != 0) {
75 shear.e[2][0] = 0.0;
76 shear.e[1][0] = d.shear.e[0];
77 *this = *this * shear;
78 }
79
80 /*
81 * Apply scale.
82 */
83 for (int i = 0; i < 3; i++) {
84 for (int j = 0; j < 3; j++) {
85 e[i][j] *= d.scale.e[i];
86 }
87 }
88 }
89
90 Matrix Matrix::operator+(const Matrix& o) const {
91 return Matrix(
92 m[0] + o.m[0], m[1] + o.m[1], m[2] + o.m[2], m[3] + o.m[3], //
93 m[4] + o.m[4], m[5] + o.m[5], m[6] + o.m[6], m[7] + o.m[7], //
94 m[8] + o.m[8], m[9] + o.m[9], m[10] + o.m[10], m[11] + o.m[11], //
95 m[12] + o.m[12], m[13] + o.m[13], m[14] + o.m[14], m[15] + o.m[15] //
96 );
97 }
98
99 Matrix Matrix::Invert() const {
100 Matrix tmp{
101 m[5] * m[10] * m[15] - m[5] * m[11] * m[14] - m[9] * m[6] * m[15] +
102 m[9] * m[7] * m[14] + m[13] * m[6] * m[11] - m[13] * m[7] * m[10],
103
104 -m[1] * m[10] * m[15] + m[1] * m[11] * m[14] + m[9] * m[2] * m[15] -
105 m[9] * m[3] * m[14] - m[13] * m[2] * m[11] + m[13] * m[3] * m[10],
106
107 m[1] * m[6] * m[15] - m[1] * m[7] * m[14] - m[5] * m[2] * m[15] +
108 m[5] * m[3] * m[14] + m[13] * m[2] * m[7] - m[13] * m[3] * m[6],
109
110 -m[1] * m[6] * m[11] + m[1] * m[7] * m[10] + m[5] * m[2] * m[11] -
111 m[5] * m[3] * m[10] - m[9] * m[2] * m[7] + m[9] * m[3] * m[6],
112
113 -m[4] * m[10] * m[15] + m[4] * m[11] * m[14] + m[8] * m[6] * m[15] -
114 m[8] * m[7] * m[14] - m[12] * m[6] * m[11] + m[12] * m[7] * m[10],
115
116 m[0] * m[10] * m[15] - m[0] * m[11] * m[14] - m[8] * m[2] * m[15] +
117 m[8] * m[3] * m[14] + m[12] * m[2] * m[11] - m[12] * m[3] * m[10],
118
119 -m[0] * m[6] * m[15] + m[0] * m[7] * m[14] + m[4] * m[2] * m[15] -
120 m[4] * m[3] * m[14] - m[12] * m[2] * m[7] + m[12] * m[3] * m[6],
121
122 m[0] * m[6] * m[11] - m[0] * m[7] * m[10] - m[4] * m[2] * m[11] +
123 m[4] * m[3] * m[10] + m[8] * m[2] * m[7] - m[8] * m[3] * m[6],
124
125 m[4] * m[9] * m[15] - m[4] * m[11] * m[13] - m[8] * m[5] * m[15] +
126 m[8] * m[7] * m[13] + m[12] * m[5] * m[11] - m[12] * m[7] * m[9],
127
128 -m[0] * m[9] * m[15] + m[0] * m[11] * m[13] + m[8] * m[1] * m[15] -
129 m[8] * m[3] * m[13] - m[12] * m[1] * m[11] + m[12] * m[3] * m[9],
130
131 m[0] * m[5] * m[15] - m[0] * m[7] * m[13] - m[4] * m[1] * m[15] +
132 m[4] * m[3] * m[13] + m[12] * m[1] * m[7] - m[12] * m[3] * m[5],
133
134 -m[0] * m[5] * m[11] + m[0] * m[7] * m[9] + m[4] * m[1] * m[11] -
135 m[4] * m[3] * m[9] - m[8] * m[1] * m[7] + m[8] * m[3] * m[5],
136
137 -m[4] * m[9] * m[14] + m[4] * m[10] * m[13] + m[8] * m[5] * m[14] -
138 m[8] * m[6] * m[13] - m[12] * m[5] * m[10] + m[12] * m[6] * m[9],
139
140 m[0] * m[9] * m[14] - m[0] * m[10] * m[13] - m[8] * m[1] * m[14] +
141 m[8] * m[2] * m[13] + m[12] * m[1] * m[10] - m[12] * m[2] * m[9],
142
143 -m[0] * m[5] * m[14] + m[0] * m[6] * m[13] + m[4] * m[1] * m[14] -
144 m[4] * m[2] * m[13] - m[12] * m[1] * m[6] + m[12] * m[2] * m[5],
145
146 m[0] * m[5] * m[10] - m[0] * m[6] * m[9] - m[4] * m[1] * m[10] +
147 m[4] * m[2] * m[9] + m[8] * m[1] * m[6] - m[8] * m[2] * m[5]};
148
149 Scalar det =
150 m[0] * tmp.m[0] + m[1] * tmp.m[4] + m[2] * tmp.m[8] + m[3] * tmp.m[12];
151
152 if (det == 0) {
153 return {};
154 }
155
156 det = 1.0 / det;
157
158 return {tmp.m[0] * det, tmp.m[1] * det, tmp.m[2] * det, tmp.m[3] * det,
159 tmp.m[4] * det, tmp.m[5] * det, tmp.m[6] * det, tmp.m[7] * det,
160 tmp.m[8] * det, tmp.m[9] * det, tmp.m[10] * det, tmp.m[11] * det,
161 tmp.m[12] * det, tmp.m[13] * det, tmp.m[14] * det, tmp.m[15] * det};
162 }
163
164 Scalar Matrix::GetDeterminant() const {
165 auto a00 = e[0][0];
166 auto a01 = e[0][1];
167 auto a02 = e[0][2];
168 auto a03 = e[0][3];
169 auto a10 = e[1][0];
170 auto a11 = e[1][1];
171 auto a12 = e[1][2];
172 auto a13 = e[1][3];
173 auto a20 = e[2][0];
174 auto a21 = e[2][1];
175 auto a22 = e[2][2];
176 auto a23 = e[2][3];
177 auto a30 = e[3][0];
178 auto a31 = e[3][1];
179 auto a32 = e[3][2];
180 auto a33 = e[3][3];
181
182 auto b00 = a00 * a11 - a01 * a10;
183 auto b01 = a00 * a12 - a02 * a10;
184 auto b02 = a00 * a13 - a03 * a10;
185 auto b03 = a01 * a12 - a02 * a11;
186 auto b04 = a01 * a13 - a03 * a11;
187 auto b05 = a02 * a13 - a03 * a12;
188 auto b06 = a20 * a31 - a21 * a30;
189 auto b07 = a20 * a32 - a22 * a30;
190 auto b08 = a20 * a33 - a23 * a30;
191 auto b09 = a21 * a32 - a22 * a31;
192 auto b10 = a21 * a33 - a23 * a31;
193 auto b11 = a22 * a33 - a23 * a32;
194
195 return b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
196 }
197
198 /*
199 * Adapted for Impeller from Graphics Gems:
200 * http://www.realtimerendering.com/resources/GraphicsGems/gemsii/unmatrix.c
201 */
202 std::optional<MatrixDecomposition> Matrix::Decompose() const {
203 /*
204 * Normalize the matrix.
205 */
206 Matrix self = *this;
207
208 if (self.e[3][3] == 0) {
209 return std::nullopt;
210 }
211
212 for (int i = 0; i < 4; i++) {
213 for (int j = 0; j < 4; j++) {
214 self.e[i][j] /= self.e[3][3];
215 }
216 }
217
218 /*
219 * `perspectiveMatrix` is used to solve for perspective, but it also provides
220 * an easy way to test for singularity of the upper 3x3 component.
221 */
222 Matrix perpectiveMatrix = self;
223 for (int i = 0; i < 3; i++) {
224 perpectiveMatrix.e[i][3] = 0;
225 }
226
227 perpectiveMatrix.e[3][3] = 1;
228
229 if (!perpectiveMatrix.IsInvertible()) {
230 return std::nullopt;
231 }
232
233 MatrixDecomposition result;
234
235 /*
236 * ==========================================================================
237 * First, isolate perspective.
238 * ==========================================================================
239 */
240 if (self.e[0][3] != 0.0 || self.e[1][3] != 0.0 || self.e[2][3] != 0.0) {
241 /*
242 * prhs is the right hand side of the equation.
243 */
244 const Vector4 rightHandSide(self.e[0][3], //
245 self.e[1][3], //
246 self.e[2][3], //
247 self.e[3][3]);
248
249 /*
250 * Solve the equation by inverting `perspectiveMatrix` and multiplying
251 * prhs by the inverse.
252 */
253
254 result.perspective = perpectiveMatrix.Invert().Transpose() * rightHandSide;
255
256 /*
257 * Clear the perspective partition.
258 */
259 self.e[0][3] = self.e[1][3] = self.e[2][3] = 0;
260 self.e[3][3] = 1;
261 }
262
263 /*
264 * ==========================================================================
265 * Next, the translation.
266 * ==========================================================================
267 */
268 result.translation = {self.e[3][0], self.e[3][1], self.e[3][2]};
269 self.e[3][0] = self.e[3][1] = self.e[3][2] = 0.0;
270
271 /*
272 * ==========================================================================
273 * Next, the scale and shear.
274 * ==========================================================================
275 */
276 Vector3 row[3];
277 for (int i = 0; i < 3; i++) {
278 row[i].x = self.e[i][0];
279 row[i].y = self.e[i][1];
280 row[i].z = self.e[i][2];
281 }
282
283 /*
284 * Compute X scale factor and normalize first row.
285 */
286 result.scale.x = row[0].GetLength();
287 row[0] = row[0].Normalize();
288
289 /*
290 * Compute XY shear factor and make 2nd row orthogonal to 1st.
291 */
292 result.shear.xy = row[0].Dot(row[1]);
293 row[1] = Vector3::Combine(row[1], 1.0, row[0], -result.shear.xy);
294
295 /*
296 * Compute Y scale and normalize 2nd row.
297 */
298 result.scale.y = row[1].GetLength();
299 row[1] = row[1].Normalize();
300 result.shear.xy /= result.scale.y;
301
302 /*
303 * Compute XZ and YZ shears, orthogonalize 3rd row.
304 */
305 result.shear.xz = row[0].Dot(row[2]);
306 row[2] = Vector3::Combine(row[2], 1.0, row[0], -result.shear.xz);
307 result.shear.yz = row[1].Dot(row[2]);
308 row[2] = Vector3::Combine(row[2], 1.0, row[1], -result.shear.yz);
309
310 /*
311 * Next, get Z scale and normalize 3rd row.
312 */
313 result.scale.z = row[2].GetLength();
314 row[2] = row[2].Normalize();
315
316 result.shear.xz /= result.scale.z;
317 result.shear.yz /= result.scale.z;
318
319 /*
320 * At this point, the matrix (in rows[]) is orthonormal.
321 * Check for a coordinate system flip. If the determinant
322 * is -1, then negate the matrix and the scaling factors.
323 */
324 if (row[0].Dot(row[1].Cross(row[2])) < 0) {
325 result.scale.x *= -1;
326 result.scale.y *= -1;
327 result.scale.z *= -1;
328
329 for (int i = 0; i < 3; i++) {
330 row[i].x *= -1;
331 row[i].y *= -1;
332 row[i].z *= -1;
333 }
334 }
335
336 /*
337 * ==========================================================================
338 * Finally, get the rotations out.
339 * ==========================================================================
340 */
341 result.rotation.x =
342 0.5 * sqrt(fmax(1.0 + row[0].x - row[1].y - row[2].z, 0.0));
343 result.rotation.y =
344 0.5 * sqrt(fmax(1.0 - row[0].x + row[1].y - row[2].z, 0.0));
345 result.rotation.z =
346 0.5 * sqrt(fmax(1.0 - row[0].x - row[1].y + row[2].z, 0.0));
347 result.rotation.w =
348 0.5 * sqrt(fmax(1.0 + row[0].x + row[1].y + row[2].z, 0.0));
349
350 if (row[2].y > row[1].z) {
351 result.rotation.x = -result.rotation.x;
352 }
353 if (row[0].z > row[2].x) {
354 result.rotation.y = -result.rotation.y;
355 }
356 if (row[1].x > row[0].y) {
357 result.rotation.z = -result.rotation.z;
358 }
359
360 return result;
361 }
362
363 std::optional<std::pair<Scalar, Scalar>> Matrix::GetScales2D() const {
364 if (HasPerspective2D()) {
365 return std::nullopt;
366 }
367
368 // We only operate on the uppermost 2x2 matrix since those are the only
369 // values that can induce a scale on 2D coordinates.
370 // [ a b ]
371 // [ c d ]
372 double a = m[0];
373 double b = m[1];
374 double c = m[4];
375 double d = m[5];
376
377 if (b == 0.0f && c == 0.0f) {
378 return {{std::abs(a), std::abs(d)}};
379 }
380
381 if (a == 0.0f && d == 0.0f) {
382 return {{std::abs(b), std::abs(c)}};
383 }
384
385 // Compute eigenvalues for the matrix (transpose(A) * A):
386 // [ a2 b2 ] == [ a b ] [ a c ] == [ aa + bb ac + bd ]
387 // [ c2 d2 ] [ c d ] [ b d ] [ ac + bd cc + dd ]
388 // (note the reverse diagonal entries in the answer are identical)
389 double a2 = a * a + b * b;
390 double b2 = a * c + b * d;
391 double c2 = b2;
392 double d2 = c * c + d * d;
393
394 //
395 // If L is an eigenvalue, then
396 // det(this - L*Identity) == 0
397 // det([ a - L b ]
398 // [ c d - L ]) == 0
399 // (a - L) * (d - L) - bc == 0
400 // ad - aL - dL + L^2 - bc == 0
401 // L^2 + (-a + -d)L + ad - bc == 0
402 //
403 // Using quadratic equation for (Ax^2 + Bx + C):
404 // A == 1
405 // B == -(a2 + d2)
406 // C == a2d2 - b2c2
407 //
408 // (We use -B for calculations because the square is the same as B and we
409 // need -B for the final quadratic equation computations anyway.)
410 double minus_B = a2 + d2;
411 double C = a2 * d2 - b2 * c2;
412 double B_squared_minus_4AC = minus_B * minus_B - 4 * 1.0f * C;
413
414 double quadratic_sqrt;
415 if (B_squared_minus_4AC <= 0.0f) {
416 // This test should never fail, but we might be slightly negative
417 FML_DCHECK(B_squared_minus_4AC + kEhCloseEnough >= 0.0f);
418 // Uniform scales (possibly rotated) would tend to end up here
419 // in which case both eigenvalues are identical
420 quadratic_sqrt = 0.0f;
421 } else {
422 quadratic_sqrt = std::sqrt(B_squared_minus_4AC);
423 }
424
425 // Since this is returning the sqrt of the values, we can guarantee that
426 // the returned scales are non-negative.
427 FML_DCHECK(minus_B - quadratic_sqrt >= 0.0f);
428 FML_DCHECK(minus_B + quadratic_sqrt >= 0.0f);
429 return {{std::sqrt((minus_B - quadratic_sqrt) / 2.0f),
430 std::sqrt((minus_B + quadratic_sqrt) / 2.0f)}};
431 }
432
433 uint64_t MatrixDecomposition::GetComponentsMask() const {
434 uint64_t mask = 0;
435
436 Quaternion noRotation(0.0, 0.0, 0.0, 1.0);
437 if (rotation != noRotation) {
438 mask = mask | static_cast<uint64_t>(Component::kRotation);
439 }
440
441 Vector4 defaultPerspective(0.0, 0.0, 0.0, 1.0);
442 if (perspective != defaultPerspective) {
443 mask = mask | static_cast<uint64_t>(Component::kPerspective);
444 }
445
446 Shear noShear(0.0, 0.0, 0.0);
447 if (shear != noShear) {
448 mask = mask | static_cast<uint64_t>(Component::kShear);
449 }
450
451 Vector3 defaultScale(1.0, 1.0, 1.0);
452 if (scale != defaultScale) {
453 mask = mask | static_cast<uint64_t>(Component::kScale);
454 }
455
456 Vector3 defaultTranslation(0.0, 0.0, 0.0);
457 if (translation != defaultTranslation) {
458 mask = mask | static_cast<uint64_t>(Component::kTranslation);
459 }
460
461 return mask;
462 }
463
464 } // namespace impeller
x
int32_t x
Definition: dl_golden_unittests.cc:458
matrix.h
impeller::saturated::b
SI b
Definition: saturated_math.h:87
impeller
Definition: allocation.cc:12
impeller::Scalar
float Scalar
Definition: scalar.h:19
impeller::kEhCloseEnough
constexpr float kEhCloseEnough
Definition: constants.h:57
c
Point c
Definition: round_superellipse_param.cc:315
impeller::MatrixDecomposition
Definition: matrix_decomposition.h:15
impeller::MatrixDecomposition::shear
Shear shear
Definition: matrix_decomposition.h:18
impeller::MatrixDecomposition::scale
Vector3 scale
Definition: matrix_decomposition.h:17
impeller::MatrixDecomposition::Component::kTranslation
@ kTranslation
impeller::MatrixDecomposition::Component::kScale
@ kScale
impeller::MatrixDecomposition::Component::kRotation
@ kRotation
impeller::MatrixDecomposition::Component::kPerspective
@ kPerspective
impeller::MatrixDecomposition::Component::kShear
@ kShear
impeller::MatrixDecomposition::perspective
Vector4 perspective
Definition: matrix_decomposition.h:19
impeller::MatrixDecomposition::GetComponentsMask
uint64_t GetComponentsMask() const
Definition: matrix.cc:433
impeller::MatrixDecomposition::rotation
Quaternion rotation
Definition: matrix_decomposition.h:20
impeller::MatrixDecomposition::translation
Vector3 translation
Definition: matrix_decomposition.h:16
impeller::Matrix
A 4x4 matrix using column-major storage.
Definition: matrix.h:37
impeller::Matrix::Matrix
constexpr Matrix()
Definition: matrix.h:47
impeller::Matrix::m
Scalar m[16]
Definition: matrix.h:39
impeller::Matrix::IsInvertible
bool IsInvertible() const
Definition: matrix.h:321
impeller::Matrix::operator+
Matrix operator+(const Vector3 &t) const
Definition: matrix.h:559
impeller::Matrix::Invert
Matrix Invert() const
Definition: matrix.cc:99
impeller::Matrix::Decompose
std::optional< MatrixDecomposition > Decompose() const
Definition: matrix.cc:202
impeller::Matrix::HasPerspective2D
constexpr bool HasPerspective2D() const
Definition: matrix.h:414
impeller::Matrix::e
Scalar e[4][4]
Definition: matrix.h:40
impeller::Matrix::GetScales2D
std::optional< std::pair< Scalar, Scalar > > GetScales2D() const
Compute the two non-negative scales applied by this matrix to 2D coordinates and return them as an op...
Definition: matrix.cc:363
impeller::Matrix::GetDeterminant
Scalar GetDeterminant() const
Definition: matrix.cc:164
impeller::Matrix::Transpose
constexpr Matrix Transpose() const
Definition: matrix.h:306
impeller::Quaternion
Definition: quaternion.h:14
impeller::Quaternion::y
Scalar y
Definition: quaternion.h:18
impeller::Quaternion::z
Scalar z
Definition: quaternion.h:19
impeller::Quaternion::w
Scalar w
Definition: quaternion.h:20
impeller::Quaternion::x
Scalar x
Definition: quaternion.h:17
impeller::Shear
Definition: shear.h:12
impeller::Shear::yz
double yz
Definition: shear.h:17
impeller::Shear::xy
double xy
Definition: shear.h:15
impeller::Shear::xz
double xz
Definition: shear.h:16
impeller::Shear::e
double e[3]
Definition: shear.h:19
impeller::Vector3
Definition: vector.h:20
impeller::Vector3::Normalize
Vector3 Normalize() const
Definition: vector.h:49
impeller::Vector3::x
Scalar x
Definition: vector.h:23
impeller::Vector3::e
Scalar e[3]
Definition: vector.h:27
impeller::Vector3::y
Scalar y
Definition: vector.h:24
impeller::Vector3::Combine
static constexpr Vector3 Combine(const Vector3 &a, Scalar aScale, const Vector3 &b, Scalar bScale)
Definition: vector.h:192
impeller::Vector3::Dot
constexpr Scalar Dot(const Vector3 &other) const
Definition: vector.h:54
impeller::Vector3::GetLength
Scalar GetLength() const
Definition: vector.h:47
impeller::Vector3::z
Scalar z
Definition: vector.h:25
impeller::Vector4
Definition: vector.h:232
impeller::Vector4::e
Scalar e[4]
Definition: vector.h:240