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