| 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 Scalar Matrix::GetMaxBasisLengthXY() const {
206 Scalar max = 0;
207 for (int i = 0; i < 3; i++) {
208 max = std::max(max, e[i][0] * e[i][0] + e[i][1] * e[i][1]);
209 }
210 return std::sqrt(max);
211 }
212
213 /*
214 * Adapted for Impeller from Graphics Gems:
215 * http://www.realtimerendering.com/resources/GraphicsGems/gemsii/unmatrix.c
216 */
217 std::optional<MatrixDecomposition> Matrix::Decompose() const {
218 /*
219 * Normalize the matrix.
220 */
221 Matrix self = *this;
222
223 if (self.e[3][3] == 0) {
224 return std::nullopt;
225 }
226
227 for (int i = 0; i < 4; i++) {
228 for (int j = 0; j < 4; j++) {
229 self.e[i][j] /= self.e[3][3];
230 }
231 }
232
233 /*
234 * `perspectiveMatrix` is used to solve for perspective, but it also provides
235 * an easy way to test for singularity of the upper 3x3 component.
236 */
237 Matrix perpectiveMatrix = self;
238 for (int i = 0; i < 3; i++) {
239 perpectiveMatrix.e[i][3] = 0;
240 }
241
242 perpectiveMatrix.e[3][3] = 1;
243
244 if (perpectiveMatrix.GetDeterminant() == 0.0) {
245 return std::nullopt;
246 }
247
248 MatrixDecomposition result;
249
250 /*
251 * ==========================================================================
252 * First, isolate perspective.
253 * ==========================================================================
254 */
255 if (self.e[0][3] != 0.0 || self.e[1][3] != 0.0 || self.e[2][3] != 0.0) {
256 /*
257 * prhs is the right hand side of the equation.
258 */
259 const Vector4 rightHandSide(self.e[0][3], //
260 self.e[1][3], //
261 self.e[2][3], //
262 self.e[3][3]);
263
264 /*
265 * Solve the equation by inverting `perspectiveMatrix` and multiplying
266 * prhs by the inverse.
267 */
268
269 result.perspective = perpectiveMatrix.Invert().Transpose() * rightHandSide;
270
271 /*
272 * Clear the perspective partition.
273 */
274 self.e[0][3] = self.e[1][3] = self.e[2][3] = 0;
275 self.e[3][3] = 1;
276 }
277
278 /*
279 * ==========================================================================
280 * Next, the translation.
281 * ==========================================================================
282 */
283 result.translation = {self.e[3][0], self.e[3][1], self.e[3][2]};
284 self.e[3][0] = self.e[3][1] = self.e[3][2] = 0.0;
285
286 /*
287 * ==========================================================================
288 * Next, the scale and shear.
289 * ==========================================================================
290 */
291 Vector3 row[3];
292 for (int i = 0; i < 3; i++) {
293 row[i].x = self.e[i][0];
294 row[i].y = self.e[i][1];
295 row[i].z = self.e[i][2];
296 }
297
298 /*
299 * Compute X scale factor and normalize first row.
300 */
301 result.scale.x = row[0].Length();
302 row[0] = row[0].Normalize();
303
304 /*
305 * Compute XY shear factor and make 2nd row orthogonal to 1st.
306 */
307 result.shear.xy = row[0].Dot(row[1]);
308 row[1] = Vector3::Combine(row[1], 1.0, row[0], -result.shear.xy);
309
310 /*
311 * Compute Y scale and normalize 2nd row.
312 */
313 result.scale.y = row[1].Length();
314 row[1] = row[1].Normalize();
315 result.shear.xy /= result.scale.y;
316
317 /*
318 * Compute XZ and YZ shears, orthogonalize 3rd row.
319 */
320 result.shear.xz = row[0].Dot(row[2]);
321 row[2] = Vector3::Combine(row[2], 1.0, row[0], -result.shear.xz);
322 result.shear.yz = row[1].Dot(row[2]);
323 row[2] = Vector3::Combine(row[2], 1.0, row[1], -result.shear.yz);
324
325 /*
326 * Next, get Z scale and normalize 3rd row.
327 */
328 result.scale.z = row[2].Length();
329 row[2] = row[2].Normalize();
330
331 result.shear.xz /= result.scale.z;
332 result.shear.yz /= result.scale.z;
333
334 /*
335 * At this point, the matrix (in rows[]) is orthonormal.
336 * Check for a coordinate system flip. If the determinant
337 * is -1, then negate the matrix and the scaling factors.
338 */
339 if (row[0].Dot(row[1].Cross(row[2])) < 0) {
340 result.scale.x *= -1;
341 result.scale.y *= -1;
342 result.scale.z *= -1;
343
344 for (int i = 0; i < 3; i++) {
345 row[i].x *= -1;
346 row[i].y *= -1;
347 row[i].z *= -1;
348 }
349 }
350
351 /*
352 * ==========================================================================
353 * Finally, get the rotations out.
354 * ==========================================================================
355 */
356 result.rotation.x =
357 0.5 * sqrt(fmax(1.0 + row[0].x - row[1].y - row[2].z, 0.0));
358 result.rotation.y =
359 0.5 * sqrt(fmax(1.0 - row[0].x + row[1].y - row[2].z, 0.0));
360 result.rotation.z =
361 0.5 * sqrt(fmax(1.0 - row[0].x - row[1].y + row[2].z, 0.0));
362 result.rotation.w =
363 0.5 * sqrt(fmax(1.0 + row[0].x + row[1].y + row[2].z, 0.0));
364
365 if (row[2].y > row[1].z) {
366 result.rotation.x = -result.rotation.x;
367 }
368 if (row[0].z > row[2].x) {
369 result.rotation.y = -result.rotation.y;
370 }
371 if (row[1].x > row[0].y) {
372 result.rotation.z = -result.rotation.z;
373 }
374
375 return result;
376 }
377
378 uint64_t MatrixDecomposition::GetComponentsMask() const {
379 uint64_t mask = 0;
380
381 Quaternion noRotation(0.0, 0.0, 0.0, 1.0);
382 if (rotation != noRotation) {
383 mask = mask | static_cast<uint64_t>(Component::kRotation);
384 }
385
386 Vector4 defaultPerspective(0.0, 0.0, 0.0, 1.0);
387 if (perspective != defaultPerspective) {
388 mask = mask | static_cast<uint64_t>(Component::kPerspective);
389 }
390
391 Shear noShear(0.0, 0.0, 0.0);
392 if (shear != noShear) {
393 mask = mask | static_cast<uint64_t>(Component::kShear);
394 }
395
396 Vector3 defaultScale(1.0, 1.0, 1.0);
397 if (scale != defaultScale) {
398 mask = mask | static_cast<uint64_t>(Component::kScale);
399 }
400
401 Vector3 defaultTranslation(0.0, 0.0, 0.0);
402 if (translation != defaultTranslation) {
403 mask = mask | static_cast<uint64_t>(Component::kTranslation);
404 }
405
406 return mask;
407 }
408
409 } // namespace impeller
impeller::Matrix::operator+
Matrix operator+(const Vector3 &t) const
Definition: matrix.h:391
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:217
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:378
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_context.cc:10
impeller::Matrix
A 4x4 matrix using column-major storage.
Definition: matrix.h:37
impeller::Vector3
Definition: vector.h:20
impeller::Matrix::GetMaxBasisLengthXY
Scalar GetMaxBasisLengthXY() const
Definition: matrix.cc:205