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 
6 
7 #include <climits>
8 #include <sstream>
9 
10 namespace impeller {
11 
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 
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 
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.IsInvertible()) {
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 
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
int32_t x
float Scalar
Definition: scalar.h:19
uint64_t GetComponentsMask() const
Definition: matrix.cc:361
A 4x4 matrix using column-major storage.
Definition: matrix.h:37
constexpr Matrix()
Definition: matrix.h:47
Scalar m[16]
Definition: matrix.h:39
bool IsInvertible() const
Definition: matrix.h:321
Matrix operator+(const Vector3 &t) const
Definition: matrix.h:492
Matrix Invert() const
Definition: matrix.cc:97
std::optional< MatrixDecomposition > Decompose() const
Definition: matrix.cc:200
Scalar e[4][4]
Definition: matrix.h:40
Scalar GetDeterminant() const
Definition: matrix.cc:162
constexpr Matrix Transpose() const
Definition: matrix.h:306
double yz
Definition: shear.h:17
double xy
Definition: shear.h:15
double xz
Definition: shear.h:16
double e[3]
Definition: shear.h:19
Vector3 Normalize() const
Definition: vector.h:49
Scalar e[3]
Definition: vector.h:27
static constexpr Vector3 Combine(const Vector3 &a, Scalar aScale, const Vector3 &b, Scalar bScale)
Definition: vector.h:192
constexpr Scalar Dot(const Vector3 &other) const
Definition: vector.h:54
Scalar GetLength() const
Definition: vector.h:47
Scalar e[4]
Definition: vector.h:240