1 /*
2 * Copyright (C) 2005, 2006 Apple Computer, Inc. All rights reserved.
3 * Copyright (C) 2009 Torch Mobile, Inc.
4 * Copyright (C) 2013 Google Inc. All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 * 2. Redistributions in binary form must reproduce the above copyright
12 * notice, this list of conditions and the following disclaimer in the
13 * documentation and/or other materials provided with the distribution.
14 *
15 * THIS SOFTWARE IS PROVIDED BY APPLE COMPUTER, INC. ``AS IS'' AND ANY
16 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
18 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL APPLE COMPUTER, INC. OR
19 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
20 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
21 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
22 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
23 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 */
27
28 #include "config.h"
29 #include "platform/transforms/TransformationMatrix.h"
30
31 #include "platform/geometry/FloatBox.h"
32 #include "platform/geometry/FloatQuad.h"
33 #include "platform/geometry/FloatRect.h"
34 #include "platform/geometry/IntRect.h"
35 #include "platform/geometry/LayoutRect.h"
36 #include "platform/transforms/AffineTransform.h"
37
38 #include "wtf/Assertions.h"
39 #include "wtf/MathExtras.h"
40
41 #if CPU(X86_64)
42 #include <emmintrin.h>
43 #endif
44
45 using namespace std;
46
47 namespace WebCore {
48
49 //
50 // Supporting Math Functions
51 //
52 // This is a set of function from various places (attributed inline) to do things like
53 // inversion and decomposition of a 4x4 matrix. They are used throughout the code
54 //
55
56 //
57 // Adapted from Matrix Inversion by Richard Carling, Graphics Gems <http://tog.acm.org/GraphicsGems/index.html>.
58
59 // EULA: The Graphics Gems code is copyright-protected. In other words, you cannot claim the text of the code
60 // as your own and resell it. Using the code is permitted in any program, product, or library, non-commercial
61 // or commercial. Giving credit is not required, though is a nice gesture. The code comes as-is, and if there
62 // are any flaws or problems with any Gems code, nobody involved with Gems - authors, editors, publishers, or
63 // webmasters - are to be held responsible. Basically, don't be a jerk, and remember that anything free comes
64 // with no guarantee.
65
66 // A clarification about the storage of matrix elements
67 //
68 // This class uses a 2 dimensional array internally to store the elements of the matrix. The first index into
69 // the array refers to the column that the element lies in; the second index refers to the row.
70 //
71 // In other words, this is the layout of the matrix:
72 //
73 // | m_matrix[0][0] m_matrix[1][0] m_matrix[2][0] m_matrix[3][0] |
74 // | m_matrix[0][1] m_matrix[1][1] m_matrix[2][1] m_matrix[3][1] |
75 // | m_matrix[0][2] m_matrix[1][2] m_matrix[2][2] m_matrix[3][2] |
76 // | m_matrix[0][3] m_matrix[1][3] m_matrix[2][3] m_matrix[3][3] |
77
78 typedef double Vector4[4];
79 typedef double Vector3[3];
80
81 const double SMALL_NUMBER = 1.e-8;
82
83 // inverse(original_matrix, inverse_matrix)
84 //
85 // calculate the inverse of a 4x4 matrix
86 //
87 // -1
88 // A = ___1__ adjoint A
89 // det A
90
91 // double = determinant2x2(double a, double b, double c, double d)
92 //
93 // calculate the determinant of a 2x2 matrix.
94
determinant2x2(double a,double b,double c,double d)95 static double determinant2x2(double a, double b, double c, double d)
96 {
97 return a * d - b * c;
98 }
99
100 // double = determinant3x3(a1, a2, a3, b1, b2, b3, c1, c2, c3)
101 //
102 // Calculate the determinant of a 3x3 matrix
103 // in the form
104 //
105 // | a1, b1, c1 |
106 // | a2, b2, c2 |
107 // | a3, b3, c3 |
108
determinant3x3(double a1,double a2,double a3,double b1,double b2,double b3,double c1,double c2,double c3)109 static double determinant3x3(double a1, double a2, double a3, double b1, double b2, double b3, double c1, double c2, double c3)
110 {
111 return a1 * determinant2x2(b2, b3, c2, c3)
112 - b1 * determinant2x2(a2, a3, c2, c3)
113 + c1 * determinant2x2(a2, a3, b2, b3);
114 }
115
116 // double = determinant4x4(matrix)
117 //
118 // calculate the determinant of a 4x4 matrix.
119
determinant4x4(const TransformationMatrix::Matrix4 & m)120 static double determinant4x4(const TransformationMatrix::Matrix4& m)
121 {
122 // Assign to individual variable names to aid selecting
123 // correct elements
124
125 double a1 = m[0][0];
126 double b1 = m[0][1];
127 double c1 = m[0][2];
128 double d1 = m[0][3];
129
130 double a2 = m[1][0];
131 double b2 = m[1][1];
132 double c2 = m[1][2];
133 double d2 = m[1][3];
134
135 double a3 = m[2][0];
136 double b3 = m[2][1];
137 double c3 = m[2][2];
138 double d3 = m[2][3];
139
140 double a4 = m[3][0];
141 double b4 = m[3][1];
142 double c4 = m[3][2];
143 double d4 = m[3][3];
144
145 return a1 * determinant3x3(b2, b3, b4, c2, c3, c4, d2, d3, d4)
146 - b1 * determinant3x3(a2, a3, a4, c2, c3, c4, d2, d3, d4)
147 + c1 * determinant3x3(a2, a3, a4, b2, b3, b4, d2, d3, d4)
148 - d1 * determinant3x3(a2, a3, a4, b2, b3, b4, c2, c3, c4);
149 }
150
151 // adjoint( original_matrix, inverse_matrix )
152 //
153 // calculate the adjoint of a 4x4 matrix
154 //
155 // Let a denote the minor determinant of matrix A obtained by
156 // ij
157 //
158 // deleting the ith row and jth column from A.
159 //
160 // i+j
161 // Let b = (-1) a
162 // ij ji
163 //
164 // The matrix B = (b ) is the adjoint of A
165 // ij
166
adjoint(const TransformationMatrix::Matrix4 & matrix,TransformationMatrix::Matrix4 & result)167 static void adjoint(const TransformationMatrix::Matrix4& matrix, TransformationMatrix::Matrix4& result)
168 {
169 // Assign to individual variable names to aid
170 // selecting correct values
171 double a1 = matrix[0][0];
172 double b1 = matrix[0][1];
173 double c1 = matrix[0][2];
174 double d1 = matrix[0][3];
175
176 double a2 = matrix[1][0];
177 double b2 = matrix[1][1];
178 double c2 = matrix[1][2];
179 double d2 = matrix[1][3];
180
181 double a3 = matrix[2][0];
182 double b3 = matrix[2][1];
183 double c3 = matrix[2][2];
184 double d3 = matrix[2][3];
185
186 double a4 = matrix[3][0];
187 double b4 = matrix[3][1];
188 double c4 = matrix[3][2];
189 double d4 = matrix[3][3];
190
191 // Row column labeling reversed since we transpose rows & columns
192 result[0][0] = determinant3x3(b2, b3, b4, c2, c3, c4, d2, d3, d4);
193 result[1][0] = - determinant3x3(a2, a3, a4, c2, c3, c4, d2, d3, d4);
194 result[2][0] = determinant3x3(a2, a3, a4, b2, b3, b4, d2, d3, d4);
195 result[3][0] = - determinant3x3(a2, a3, a4, b2, b3, b4, c2, c3, c4);
196
197 result[0][1] = - determinant3x3(b1, b3, b4, c1, c3, c4, d1, d3, d4);
198 result[1][1] = determinant3x3(a1, a3, a4, c1, c3, c4, d1, d3, d4);
199 result[2][1] = - determinant3x3(a1, a3, a4, b1, b3, b4, d1, d3, d4);
200 result[3][1] = determinant3x3(a1, a3, a4, b1, b3, b4, c1, c3, c4);
201
202 result[0][2] = determinant3x3(b1, b2, b4, c1, c2, c4, d1, d2, d4);
203 result[1][2] = - determinant3x3(a1, a2, a4, c1, c2, c4, d1, d2, d4);
204 result[2][2] = determinant3x3(a1, a2, a4, b1, b2, b4, d1, d2, d4);
205 result[3][2] = - determinant3x3(a1, a2, a4, b1, b2, b4, c1, c2, c4);
206
207 result[0][3] = - determinant3x3(b1, b2, b3, c1, c2, c3, d1, d2, d3);
208 result[1][3] = determinant3x3(a1, a2, a3, c1, c2, c3, d1, d2, d3);
209 result[2][3] = - determinant3x3(a1, a2, a3, b1, b2, b3, d1, d2, d3);
210 result[3][3] = determinant3x3(a1, a2, a3, b1, b2, b3, c1, c2, c3);
211 }
212
213 // Returns false if the matrix is not invertible
inverse(const TransformationMatrix::Matrix4 & matrix,TransformationMatrix::Matrix4 & result)214 static bool inverse(const TransformationMatrix::Matrix4& matrix, TransformationMatrix::Matrix4& result)
215 {
216 // Calculate the adjoint matrix
217 adjoint(matrix, result);
218
219 // Calculate the 4x4 determinant
220 // If the determinant is zero,
221 // then the inverse matrix is not unique.
222 double det = determinant4x4(matrix);
223
224 if (fabs(det) < SMALL_NUMBER)
225 return false;
226
227 // Scale the adjoint matrix to get the inverse
228
229 for (int i = 0; i < 4; i++)
230 for (int j = 0; j < 4; j++)
231 result[i][j] = result[i][j] / det;
232
233 return true;
234 }
235
236 // End of code adapted from Matrix Inversion by Richard Carling
237
238 // Perform a decomposition on the passed matrix, return false if unsuccessful
239 // From Graphics Gems: unmatrix.c
240
241 // Transpose rotation portion of matrix a, return b
transposeMatrix4(const TransformationMatrix::Matrix4 & a,TransformationMatrix::Matrix4 & b)242 static void transposeMatrix4(const TransformationMatrix::Matrix4& a, TransformationMatrix::Matrix4& b)
243 {
244 for (int i = 0; i < 4; i++)
245 for (int j = 0; j < 4; j++)
246 b[i][j] = a[j][i];
247 }
248
249 // Multiply a homogeneous point by a matrix and return the transformed point
v4MulPointByMatrix(const Vector4 p,const TransformationMatrix::Matrix4 & m,Vector4 result)250 static void v4MulPointByMatrix(const Vector4 p, const TransformationMatrix::Matrix4& m, Vector4 result)
251 {
252 result[0] = (p[0] * m[0][0]) + (p[1] * m[1][0]) +
253 (p[2] * m[2][0]) + (p[3] * m[3][0]);
254 result[1] = (p[0] * m[0][1]) + (p[1] * m[1][1]) +
255 (p[2] * m[2][1]) + (p[3] * m[3][1]);
256 result[2] = (p[0] * m[0][2]) + (p[1] * m[1][2]) +
257 (p[2] * m[2][2]) + (p[3] * m[3][2]);
258 result[3] = (p[0] * m[0][3]) + (p[1] * m[1][3]) +
259 (p[2] * m[2][3]) + (p[3] * m[3][3]);
260 }
261
v3Length(Vector3 a)262 static double v3Length(Vector3 a)
263 {
264 return sqrt((a[0] * a[0]) + (a[1] * a[1]) + (a[2] * a[2]));
265 }
266
v3Scale(Vector3 v,double desiredLength)267 static void v3Scale(Vector3 v, double desiredLength)
268 {
269 double len = v3Length(v);
270 if (len != 0) {
271 double l = desiredLength / len;
272 v[0] *= l;
273 v[1] *= l;
274 v[2] *= l;
275 }
276 }
277
v3Dot(const Vector3 a,const Vector3 b)278 static double v3Dot(const Vector3 a, const Vector3 b)
279 {
280 return (a[0] * b[0]) + (a[1] * b[1]) + (a[2] * b[2]);
281 }
282
283 // Make a linear combination of two vectors and return the result.
284 // result = (a * ascl) + (b * bscl)
v3Combine(const Vector3 a,const Vector3 b,Vector3 result,double ascl,double bscl)285 static void v3Combine(const Vector3 a, const Vector3 b, Vector3 result, double ascl, double bscl)
286 {
287 result[0] = (ascl * a[0]) + (bscl * b[0]);
288 result[1] = (ascl * a[1]) + (bscl * b[1]);
289 result[2] = (ascl * a[2]) + (bscl * b[2]);
290 }
291
292 // Return the cross product result = a cross b */
v3Cross(const Vector3 a,const Vector3 b,Vector3 result)293 static void v3Cross(const Vector3 a, const Vector3 b, Vector3 result)
294 {
295 result[0] = (a[1] * b[2]) - (a[2] * b[1]);
296 result[1] = (a[2] * b[0]) - (a[0] * b[2]);
297 result[2] = (a[0] * b[1]) - (a[1] * b[0]);
298 }
299
decompose(const TransformationMatrix::Matrix4 & mat,TransformationMatrix::DecomposedType & result)300 static bool decompose(const TransformationMatrix::Matrix4& mat, TransformationMatrix::DecomposedType& result)
301 {
302 TransformationMatrix::Matrix4 localMatrix;
303 memcpy(localMatrix, mat, sizeof(TransformationMatrix::Matrix4));
304
305 // Normalize the matrix.
306 if (localMatrix[3][3] == 0)
307 return false;
308
309 int i, j;
310 for (i = 0; i < 4; i++)
311 for (j = 0; j < 4; j++)
312 localMatrix[i][j] /= localMatrix[3][3];
313
314 // perspectiveMatrix is used to solve for perspective, but it also provides
315 // an easy way to test for singularity of the upper 3x3 component.
316 TransformationMatrix::Matrix4 perspectiveMatrix;
317 memcpy(perspectiveMatrix, localMatrix, sizeof(TransformationMatrix::Matrix4));
318 for (i = 0; i < 3; i++)
319 perspectiveMatrix[i][3] = 0;
320 perspectiveMatrix[3][3] = 1;
321
322 if (determinant4x4(perspectiveMatrix) == 0)
323 return false;
324
325 // First, isolate perspective. This is the messiest.
326 if (localMatrix[0][3] != 0 || localMatrix[1][3] != 0 || localMatrix[2][3] != 0) {
327 // rightHandSide is the right hand side of the equation.
328 Vector4 rightHandSide;
329 rightHandSide[0] = localMatrix[0][3];
330 rightHandSide[1] = localMatrix[1][3];
331 rightHandSide[2] = localMatrix[2][3];
332 rightHandSide[3] = localMatrix[3][3];
333
334 // Solve the equation by inverting perspectiveMatrix and multiplying
335 // rightHandSide by the inverse. (This is the easiest way, not
336 // necessarily the best.)
337 TransformationMatrix::Matrix4 inversePerspectiveMatrix, transposedInversePerspectiveMatrix;
338 inverse(perspectiveMatrix, inversePerspectiveMatrix);
339 transposeMatrix4(inversePerspectiveMatrix, transposedInversePerspectiveMatrix);
340
341 Vector4 perspectivePoint;
342 v4MulPointByMatrix(rightHandSide, transposedInversePerspectiveMatrix, perspectivePoint);
343
344 result.perspectiveX = perspectivePoint[0];
345 result.perspectiveY = perspectivePoint[1];
346 result.perspectiveZ = perspectivePoint[2];
347 result.perspectiveW = perspectivePoint[3];
348
349 // Clear the perspective partition
350 localMatrix[0][3] = localMatrix[1][3] = localMatrix[2][3] = 0;
351 localMatrix[3][3] = 1;
352 } else {
353 // No perspective.
354 result.perspectiveX = result.perspectiveY = result.perspectiveZ = 0;
355 result.perspectiveW = 1;
356 }
357
358 // Next take care of translation (easy).
359 result.translateX = localMatrix[3][0];
360 localMatrix[3][0] = 0;
361 result.translateY = localMatrix[3][1];
362 localMatrix[3][1] = 0;
363 result.translateZ = localMatrix[3][2];
364 localMatrix[3][2] = 0;
365
366 // Vector4 type and functions need to be added to the common set.
367 Vector3 row[3], pdum3;
368
369 // Now get scale and shear.
370 for (i = 0; i < 3; i++) {
371 row[i][0] = localMatrix[i][0];
372 row[i][1] = localMatrix[i][1];
373 row[i][2] = localMatrix[i][2];
374 }
375
376 // Compute X scale factor and normalize first row.
377 result.scaleX = v3Length(row[0]);
378 v3Scale(row[0], 1.0);
379
380 // Compute XY shear factor and make 2nd row orthogonal to 1st.
381 result.skewXY = v3Dot(row[0], row[1]);
382 v3Combine(row[1], row[0], row[1], 1.0, -result.skewXY);
383
384 // Now, compute Y scale and normalize 2nd row.
385 result.scaleY = v3Length(row[1]);
386 v3Scale(row[1], 1.0);
387 result.skewXY /= result.scaleY;
388
389 // Compute XZ and YZ shears, orthogonalize 3rd row.
390 result.skewXZ = v3Dot(row[0], row[2]);
391 v3Combine(row[2], row[0], row[2], 1.0, -result.skewXZ);
392 result.skewYZ = v3Dot(row[1], row[2]);
393 v3Combine(row[2], row[1], row[2], 1.0, -result.skewYZ);
394
395 // Next, get Z scale and normalize 3rd row.
396 result.scaleZ = v3Length(row[2]);
397 v3Scale(row[2], 1.0);
398 result.skewXZ /= result.scaleZ;
399 result.skewYZ /= result.scaleZ;
400
401 // At this point, the matrix (in rows[]) is orthonormal.
402 // Check for a coordinate system flip. If the determinant
403 // is -1, then negate the matrix and the scaling factors.
404 v3Cross(row[1], row[2], pdum3);
405 if (v3Dot(row[0], pdum3) < 0) {
406
407 result.scaleX *= -1;
408 result.scaleY *= -1;
409 result.scaleZ *= -1;
410
411 for (i = 0; i < 3; i++) {
412 row[i][0] *= -1;
413 row[i][1] *= -1;
414 row[i][2] *= -1;
415 }
416 }
417
418 // Now, get the rotations out, as described in the gem.
419
420 // FIXME - Add the ability to return either quaternions (which are
421 // easier to recompose with) or Euler angles (rx, ry, rz), which
422 // are easier for authors to deal with. The latter will only be useful
423 // when we fix https://bugs.webkit.org/show_bug.cgi?id=23799, so I
424 // will leave the Euler angle code here for now.
425
426 // ret.rotateY = asin(-row[0][2]);
427 // if (cos(ret.rotateY) != 0) {
428 // ret.rotateX = atan2(row[1][2], row[2][2]);
429 // ret.rotateZ = atan2(row[0][1], row[0][0]);
430 // } else {
431 // ret.rotateX = atan2(-row[2][0], row[1][1]);
432 // ret.rotateZ = 0;
433 // }
434
435 double s, t, x, y, z, w;
436
437 t = row[0][0] + row[1][1] + row[2][2] + 1.0;
438
439 if (t > 1e-4) {
440 s = 0.5 / sqrt(t);
441 w = 0.25 / s;
442 x = (row[2][1] - row[1][2]) * s;
443 y = (row[0][2] - row[2][0]) * s;
444 z = (row[1][0] - row[0][1]) * s;
445 } else if (row[0][0] > row[1][1] && row[0][0] > row[2][2]) {
446 s = sqrt (1.0 + row[0][0] - row[1][1] - row[2][2]) * 2.0; // S=4*qx
447 x = 0.25 * s;
448 y = (row[0][1] + row[1][0]) / s;
449 z = (row[0][2] + row[2][0]) / s;
450 w = (row[2][1] - row[1][2]) / s;
451 } else if (row[1][1] > row[2][2]) {
452 s = sqrt (1.0 + row[1][1] - row[0][0] - row[2][2]) * 2.0; // S=4*qy
453 x = (row[0][1] + row[1][0]) / s;
454 y = 0.25 * s;
455 z = (row[1][2] + row[2][1]) / s;
456 w = (row[0][2] - row[2][0]) / s;
457 } else {
458 s = sqrt(1.0 + row[2][2] - row[0][0] - row[1][1]) * 2.0; // S=4*qz
459 x = (row[0][2] + row[2][0]) / s;
460 y = (row[1][2] + row[2][1]) / s;
461 z = 0.25 * s;
462 w = (row[1][0] - row[0][1]) / s;
463 }
464
465 result.quaternionX = x;
466 result.quaternionY = y;
467 result.quaternionZ = z;
468 result.quaternionW = w;
469
470 return true;
471 }
472
473 // Perform a spherical linear interpolation between the two
474 // passed quaternions with 0 <= t <= 1
slerp(double qa[4],const double qb[4],double t)475 static void slerp(double qa[4], const double qb[4], double t)
476 {
477 double ax, ay, az, aw;
478 double bx, by, bz, bw;
479 double cx, cy, cz, cw;
480 double angle;
481 double th, invth, scale, invscale;
482
483 ax = qa[0]; ay = qa[1]; az = qa[2]; aw = qa[3];
484 bx = qb[0]; by = qb[1]; bz = qb[2]; bw = qb[3];
485
486 angle = ax * bx + ay * by + az * bz + aw * bw;
487
488 if (angle < 0.0) {
489 ax = -ax; ay = -ay;
490 az = -az; aw = -aw;
491 angle = -angle;
492 }
493
494 if (angle + 1.0 > .05) {
495 if (1.0 - angle >= .05) {
496 th = acos (angle);
497 invth = 1.0 / sin (th);
498 scale = sin (th * (1.0 - t)) * invth;
499 invscale = sin (th * t) * invth;
500 } else {
501 scale = 1.0 - t;
502 invscale = t;
503 }
504 } else {
505 bx = -ay;
506 by = ax;
507 bz = -aw;
508 bw = az;
509 scale = sin(piDouble * (.5 - t));
510 invscale = sin (piDouble * t);
511 }
512
513 cx = ax * scale + bx * invscale;
514 cy = ay * scale + by * invscale;
515 cz = az * scale + bz * invscale;
516 cw = aw * scale + bw * invscale;
517
518 qa[0] = cx; qa[1] = cy; qa[2] = cz; qa[3] = cw;
519 }
520
521 // End of Supporting Math Functions
522
TransformationMatrix(const AffineTransform & t)523 TransformationMatrix::TransformationMatrix(const AffineTransform& t)
524 {
525 setMatrix(t.a(), t.b(), t.c(), t.d(), t.e(), t.f());
526 }
527
scale(double s)528 TransformationMatrix& TransformationMatrix::scale(double s)
529 {
530 return scaleNonUniform(s, s);
531 }
532
rotateFromVector(double x,double y)533 TransformationMatrix& TransformationMatrix::rotateFromVector(double x, double y)
534 {
535 return rotate(rad2deg(atan2(y, x)));
536 }
537
flipX()538 TransformationMatrix& TransformationMatrix::flipX()
539 {
540 return scaleNonUniform(-1.0, 1.0);
541 }
542
flipY()543 TransformationMatrix& TransformationMatrix::flipY()
544 {
545 return scaleNonUniform(1.0, -1.0);
546 }
547
projectPoint(const FloatPoint & p,bool * clamped) const548 FloatPoint TransformationMatrix::projectPoint(const FloatPoint& p, bool* clamped) const
549 {
550 // This is basically raytracing. We have a point in the destination
551 // plane with z=0, and we cast a ray parallel to the z-axis from that
552 // point to find the z-position at which it intersects the z=0 plane
553 // with the transform applied. Once we have that point we apply the
554 // inverse transform to find the corresponding point in the source
555 // space.
556 //
557 // Given a plane with normal Pn, and a ray starting at point R0 and
558 // with direction defined by the vector Rd, we can find the
559 // intersection point as a distance d from R0 in units of Rd by:
560 //
561 // d = -dot (Pn', R0) / dot (Pn', Rd)
562 if (clamped)
563 *clamped = false;
564
565 if (m33() == 0) {
566 // In this case, the projection plane is parallel to the ray we are trying to
567 // trace, and there is no well-defined value for the projection.
568 return FloatPoint();
569 }
570
571 double x = p.x();
572 double y = p.y();
573 double z = -(m13() * x + m23() * y + m43()) / m33();
574
575 // FIXME: use multVecMatrix()
576 double outX = x * m11() + y * m21() + z * m31() + m41();
577 double outY = x * m12() + y * m22() + z * m32() + m42();
578
579 double w = x * m14() + y * m24() + z * m34() + m44();
580 if (w <= 0) {
581 // Using int max causes overflow when other code uses the projected point. To
582 // represent infinity yet reduce the risk of overflow, we use a large but
583 // not-too-large number here when clamping.
584 const int largeNumber = 100000000 / kFixedPointDenominator;
585 outX = copysign(largeNumber, outX);
586 outY = copysign(largeNumber, outY);
587 if (clamped)
588 *clamped = true;
589 } else if (w != 1) {
590 outX /= w;
591 outY /= w;
592 }
593
594 return FloatPoint(static_cast<float>(outX), static_cast<float>(outY));
595 }
596
projectQuad(const FloatQuad & q,bool * clamped) const597 FloatQuad TransformationMatrix::projectQuad(const FloatQuad& q, bool* clamped) const
598 {
599 FloatQuad projectedQuad;
600
601 bool clamped1 = false;
602 bool clamped2 = false;
603 bool clamped3 = false;
604 bool clamped4 = false;
605
606 projectedQuad.setP1(projectPoint(q.p1(), &clamped1));
607 projectedQuad.setP2(projectPoint(q.p2(), &clamped2));
608 projectedQuad.setP3(projectPoint(q.p3(), &clamped3));
609 projectedQuad.setP4(projectPoint(q.p4(), &clamped4));
610
611 if (clamped)
612 *clamped = clamped1 || clamped2 || clamped3 || clamped4;
613
614 // If all points on the quad had w < 0, then the entire quad would not be visible to the projected surface.
615 bool everythingWasClipped = clamped1 && clamped2 && clamped3 && clamped4;
616 if (everythingWasClipped)
617 return FloatQuad();
618
619 return projectedQuad;
620 }
621
clampEdgeValue(float f)622 static float clampEdgeValue(float f)
623 {
624 ASSERT(!std::isnan(f));
625 return min<float>(max<float>(f, (-LayoutUnit::max() / 2).toFloat()), (LayoutUnit::max() / 2).toFloat());
626 }
627
clampedBoundsOfProjectedQuad(const FloatQuad & q) const628 LayoutRect TransformationMatrix::clampedBoundsOfProjectedQuad(const FloatQuad& q) const
629 {
630 FloatRect mappedQuadBounds = projectQuad(q).boundingBox();
631
632 float left = clampEdgeValue(floorf(mappedQuadBounds.x()));
633 float top = clampEdgeValue(floorf(mappedQuadBounds.y()));
634
635 float right;
636 if (std::isinf(mappedQuadBounds.x()) && std::isinf(mappedQuadBounds.width()))
637 right = (LayoutUnit::max() / 2).toFloat();
638 else
639 right = clampEdgeValue(ceilf(mappedQuadBounds.maxX()));
640
641 float bottom;
642 if (std::isinf(mappedQuadBounds.y()) && std::isinf(mappedQuadBounds.height()))
643 bottom = (LayoutUnit::max() / 2).toFloat();
644 else
645 bottom = clampEdgeValue(ceilf(mappedQuadBounds.maxY()));
646
647 return LayoutRect(LayoutUnit::clamp(left), LayoutUnit::clamp(top), LayoutUnit::clamp(right - left), LayoutUnit::clamp(bottom - top));
648 }
649
transformBox(FloatBox & box) const650 void TransformationMatrix::transformBox(FloatBox& box) const
651 {
652 FloatBox bounds;
653 bool firstPoint = true;
654 for (size_t i = 0; i < 2; ++i) {
655 for (size_t j = 0; j < 2; ++j) {
656 for (size_t k = 0; k < 2; ++k) {
657 FloatPoint3D point(box.x(), box.y(), box.z());
658 point += FloatPoint3D(i * box.width(), j * box.height(), k * box.depth());
659 point = mapPoint(point);
660 if (firstPoint) {
661 bounds.setOrigin(point);
662 firstPoint = false;
663 } else {
664 bounds.expandTo(point);
665 }
666 }
667 }
668 }
669 box = bounds;
670 }
671
mapPoint(const FloatPoint & p) const672 FloatPoint TransformationMatrix::mapPoint(const FloatPoint& p) const
673 {
674 if (isIdentityOrTranslation())
675 return FloatPoint(p.x() + static_cast<float>(m_matrix[3][0]), p.y() + static_cast<float>(m_matrix[3][1]));
676
677 return internalMapPoint(p);
678 }
679
mapPoint(const FloatPoint3D & p) const680 FloatPoint3D TransformationMatrix::mapPoint(const FloatPoint3D& p) const
681 {
682 if (isIdentityOrTranslation())
683 return FloatPoint3D(p.x() + static_cast<float>(m_matrix[3][0]),
684 p.y() + static_cast<float>(m_matrix[3][1]),
685 p.z() + static_cast<float>(m_matrix[3][2]));
686
687 return internalMapPoint(p);
688 }
689
mapRect(const IntRect & rect) const690 IntRect TransformationMatrix::mapRect(const IntRect &rect) const
691 {
692 return enclosingIntRect(mapRect(FloatRect(rect)));
693 }
694
mapRect(const LayoutRect & r) const695 LayoutRect TransformationMatrix::mapRect(const LayoutRect& r) const
696 {
697 return enclosingLayoutRect(mapRect(FloatRect(r)));
698 }
699
mapRect(const FloatRect & r) const700 FloatRect TransformationMatrix::mapRect(const FloatRect& r) const
701 {
702 if (isIdentityOrTranslation()) {
703 FloatRect mappedRect(r);
704 mappedRect.move(static_cast<float>(m_matrix[3][0]), static_cast<float>(m_matrix[3][1]));
705 return mappedRect;
706 }
707
708 FloatQuad result;
709
710 float maxX = r.maxX();
711 float maxY = r.maxY();
712 result.setP1(internalMapPoint(FloatPoint(r.x(), r.y())));
713 result.setP2(internalMapPoint(FloatPoint(maxX, r.y())));
714 result.setP3(internalMapPoint(FloatPoint(maxX, maxY)));
715 result.setP4(internalMapPoint(FloatPoint(r.x(), maxY)));
716
717 return result.boundingBox();
718 }
719
mapQuad(const FloatQuad & q) const720 FloatQuad TransformationMatrix::mapQuad(const FloatQuad& q) const
721 {
722 if (isIdentityOrTranslation()) {
723 FloatQuad mappedQuad(q);
724 mappedQuad.move(static_cast<float>(m_matrix[3][0]), static_cast<float>(m_matrix[3][1]));
725 return mappedQuad;
726 }
727
728 FloatQuad result;
729 result.setP1(internalMapPoint(q.p1()));
730 result.setP2(internalMapPoint(q.p2()));
731 result.setP3(internalMapPoint(q.p3()));
732 result.setP4(internalMapPoint(q.p4()));
733 return result;
734 }
735
scaleNonUniform(double sx,double sy)736 TransformationMatrix& TransformationMatrix::scaleNonUniform(double sx, double sy)
737 {
738 m_matrix[0][0] *= sx;
739 m_matrix[0][1] *= sx;
740 m_matrix[0][2] *= sx;
741 m_matrix[0][3] *= sx;
742
743 m_matrix[1][0] *= sy;
744 m_matrix[1][1] *= sy;
745 m_matrix[1][2] *= sy;
746 m_matrix[1][3] *= sy;
747 return *this;
748 }
749
scale3d(double sx,double sy,double sz)750 TransformationMatrix& TransformationMatrix::scale3d(double sx, double sy, double sz)
751 {
752 scaleNonUniform(sx, sy);
753
754 m_matrix[2][0] *= sz;
755 m_matrix[2][1] *= sz;
756 m_matrix[2][2] *= sz;
757 m_matrix[2][3] *= sz;
758 return *this;
759 }
760
rotate3d(double x,double y,double z,double angle)761 TransformationMatrix& TransformationMatrix::rotate3d(double x, double y, double z, double angle)
762 {
763 // Normalize the axis of rotation
764 double length = sqrt(x * x + y * y + z * z);
765 if (length == 0) {
766 // A direction vector that cannot be normalized, such as [0, 0, 0], will cause the rotation to not be applied.
767 return *this;
768 } else if (length != 1) {
769 x /= length;
770 y /= length;
771 z /= length;
772 }
773
774 // Angles are in degrees. Switch to radians.
775 angle = deg2rad(angle);
776
777 double sinTheta = sin(angle);
778 double cosTheta = cos(angle);
779
780 TransformationMatrix mat;
781
782 // Optimize cases where the axis is along a major axis
783 if (x == 1.0 && y == 0.0 && z == 0.0) {
784 mat.m_matrix[0][0] = 1.0;
785 mat.m_matrix[0][1] = 0.0;
786 mat.m_matrix[0][2] = 0.0;
787 mat.m_matrix[1][0] = 0.0;
788 mat.m_matrix[1][1] = cosTheta;
789 mat.m_matrix[1][2] = sinTheta;
790 mat.m_matrix[2][0] = 0.0;
791 mat.m_matrix[2][1] = -sinTheta;
792 mat.m_matrix[2][2] = cosTheta;
793 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
794 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
795 mat.m_matrix[3][3] = 1.0;
796 } else if (x == 0.0 && y == 1.0 && z == 0.0) {
797 mat.m_matrix[0][0] = cosTheta;
798 mat.m_matrix[0][1] = 0.0;
799 mat.m_matrix[0][2] = -sinTheta;
800 mat.m_matrix[1][0] = 0.0;
801 mat.m_matrix[1][1] = 1.0;
802 mat.m_matrix[1][2] = 0.0;
803 mat.m_matrix[2][0] = sinTheta;
804 mat.m_matrix[2][1] = 0.0;
805 mat.m_matrix[2][2] = cosTheta;
806 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
807 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
808 mat.m_matrix[3][3] = 1.0;
809 } else if (x == 0.0 && y == 0.0 && z == 1.0) {
810 mat.m_matrix[0][0] = cosTheta;
811 mat.m_matrix[0][1] = sinTheta;
812 mat.m_matrix[0][2] = 0.0;
813 mat.m_matrix[1][0] = -sinTheta;
814 mat.m_matrix[1][1] = cosTheta;
815 mat.m_matrix[1][2] = 0.0;
816 mat.m_matrix[2][0] = 0.0;
817 mat.m_matrix[2][1] = 0.0;
818 mat.m_matrix[2][2] = 1.0;
819 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
820 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
821 mat.m_matrix[3][3] = 1.0;
822 } else {
823 // This case is the rotation about an arbitrary unit vector.
824 //
825 // Formula is adapted from Wikipedia article on Rotation matrix,
826 // http://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
827 //
828 // An alternate resource with the same matrix: http://www.fastgraph.com/makegames/3drotation/
829 //
830 double oneMinusCosTheta = 1 - cosTheta;
831 mat.m_matrix[0][0] = cosTheta + x * x * oneMinusCosTheta;
832 mat.m_matrix[0][1] = y * x * oneMinusCosTheta + z * sinTheta;
833 mat.m_matrix[0][2] = z * x * oneMinusCosTheta - y * sinTheta;
834 mat.m_matrix[1][0] = x * y * oneMinusCosTheta - z * sinTheta;
835 mat.m_matrix[1][1] = cosTheta + y * y * oneMinusCosTheta;
836 mat.m_matrix[1][2] = z * y * oneMinusCosTheta + x * sinTheta;
837 mat.m_matrix[2][0] = x * z * oneMinusCosTheta + y * sinTheta;
838 mat.m_matrix[2][1] = y * z * oneMinusCosTheta - x * sinTheta;
839 mat.m_matrix[2][2] = cosTheta + z * z * oneMinusCosTheta;
840 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
841 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
842 mat.m_matrix[3][3] = 1.0;
843 }
844 multiply(mat);
845 return *this;
846 }
847
rotate3d(double rx,double ry,double rz)848 TransformationMatrix& TransformationMatrix::rotate3d(double rx, double ry, double rz)
849 {
850 // Angles are in degrees. Switch to radians.
851 rx = deg2rad(rx);
852 ry = deg2rad(ry);
853 rz = deg2rad(rz);
854
855 TransformationMatrix mat;
856
857 double sinTheta = sin(rz);
858 double cosTheta = cos(rz);
859
860 mat.m_matrix[0][0] = cosTheta;
861 mat.m_matrix[0][1] = sinTheta;
862 mat.m_matrix[0][2] = 0.0;
863 mat.m_matrix[1][0] = -sinTheta;
864 mat.m_matrix[1][1] = cosTheta;
865 mat.m_matrix[1][2] = 0.0;
866 mat.m_matrix[2][0] = 0.0;
867 mat.m_matrix[2][1] = 0.0;
868 mat.m_matrix[2][2] = 1.0;
869 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
870 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
871 mat.m_matrix[3][3] = 1.0;
872
873 TransformationMatrix rmat(mat);
874
875 sinTheta = sin(ry);
876 cosTheta = cos(ry);
877
878 mat.m_matrix[0][0] = cosTheta;
879 mat.m_matrix[0][1] = 0.0;
880 mat.m_matrix[0][2] = -sinTheta;
881 mat.m_matrix[1][0] = 0.0;
882 mat.m_matrix[1][1] = 1.0;
883 mat.m_matrix[1][2] = 0.0;
884 mat.m_matrix[2][0] = sinTheta;
885 mat.m_matrix[2][1] = 0.0;
886 mat.m_matrix[2][2] = cosTheta;
887 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
888 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
889 mat.m_matrix[3][3] = 1.0;
890
891 rmat.multiply(mat);
892
893 sinTheta = sin(rx);
894 cosTheta = cos(rx);
895
896 mat.m_matrix[0][0] = 1.0;
897 mat.m_matrix[0][1] = 0.0;
898 mat.m_matrix[0][2] = 0.0;
899 mat.m_matrix[1][0] = 0.0;
900 mat.m_matrix[1][1] = cosTheta;
901 mat.m_matrix[1][2] = sinTheta;
902 mat.m_matrix[2][0] = 0.0;
903 mat.m_matrix[2][1] = -sinTheta;
904 mat.m_matrix[2][2] = cosTheta;
905 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
906 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
907 mat.m_matrix[3][3] = 1.0;
908
909 rmat.multiply(mat);
910
911 multiply(rmat);
912 return *this;
913 }
914
translate(double tx,double ty)915 TransformationMatrix& TransformationMatrix::translate(double tx, double ty)
916 {
917 m_matrix[3][0] += tx * m_matrix[0][0] + ty * m_matrix[1][0];
918 m_matrix[3][1] += tx * m_matrix[0][1] + ty * m_matrix[1][1];
919 m_matrix[3][2] += tx * m_matrix[0][2] + ty * m_matrix[1][2];
920 m_matrix[3][3] += tx * m_matrix[0][3] + ty * m_matrix[1][3];
921 return *this;
922 }
923
translate3d(double tx,double ty,double tz)924 TransformationMatrix& TransformationMatrix::translate3d(double tx, double ty, double tz)
925 {
926 m_matrix[3][0] += tx * m_matrix[0][0] + ty * m_matrix[1][0] + tz * m_matrix[2][0];
927 m_matrix[3][1] += tx * m_matrix[0][1] + ty * m_matrix[1][1] + tz * m_matrix[2][1];
928 m_matrix[3][2] += tx * m_matrix[0][2] + ty * m_matrix[1][2] + tz * m_matrix[2][2];
929 m_matrix[3][3] += tx * m_matrix[0][3] + ty * m_matrix[1][3] + tz * m_matrix[2][3];
930 return *this;
931 }
932
translateRight(double tx,double ty)933 TransformationMatrix& TransformationMatrix::translateRight(double tx, double ty)
934 {
935 if (tx != 0) {
936 m_matrix[0][0] += m_matrix[0][3] * tx;
937 m_matrix[1][0] += m_matrix[1][3] * tx;
938 m_matrix[2][0] += m_matrix[2][3] * tx;
939 m_matrix[3][0] += m_matrix[3][3] * tx;
940 }
941
942 if (ty != 0) {
943 m_matrix[0][1] += m_matrix[0][3] * ty;
944 m_matrix[1][1] += m_matrix[1][3] * ty;
945 m_matrix[2][1] += m_matrix[2][3] * ty;
946 m_matrix[3][1] += m_matrix[3][3] * ty;
947 }
948
949 return *this;
950 }
951
translateRight3d(double tx,double ty,double tz)952 TransformationMatrix& TransformationMatrix::translateRight3d(double tx, double ty, double tz)
953 {
954 translateRight(tx, ty);
955 if (tz != 0) {
956 m_matrix[0][2] += m_matrix[0][3] * tz;
957 m_matrix[1][2] += m_matrix[1][3] * tz;
958 m_matrix[2][2] += m_matrix[2][3] * tz;
959 m_matrix[3][2] += m_matrix[3][3] * tz;
960 }
961
962 return *this;
963 }
964
skew(double sx,double sy)965 TransformationMatrix& TransformationMatrix::skew(double sx, double sy)
966 {
967 // angles are in degrees. Switch to radians
968 sx = deg2rad(sx);
969 sy = deg2rad(sy);
970
971 TransformationMatrix mat;
972 mat.m_matrix[0][1] = tan(sy); // note that the y shear goes in the first row
973 mat.m_matrix[1][0] = tan(sx); // and the x shear in the second row
974
975 multiply(mat);
976 return *this;
977 }
978
applyPerspective(double p)979 TransformationMatrix& TransformationMatrix::applyPerspective(double p)
980 {
981 TransformationMatrix mat;
982 if (p != 0)
983 mat.m_matrix[2][3] = -1/p;
984
985 multiply(mat);
986 return *this;
987 }
988
rectToRect(const FloatRect & from,const FloatRect & to)989 TransformationMatrix TransformationMatrix::rectToRect(const FloatRect& from, const FloatRect& to)
990 {
991 ASSERT(!from.isEmpty());
992 return TransformationMatrix(to.width() / from.width(),
993 0, 0,
994 to.height() / from.height(),
995 to.x() - from.x(),
996 to.y() - from.y());
997 }
998
999 // this = mat * this.
multiply(const TransformationMatrix & mat)1000 TransformationMatrix& TransformationMatrix::multiply(const TransformationMatrix& mat)
1001 {
1002 #if CPU(APPLE_ARMV7S)
1003 double* leftMatrix = &(m_matrix[0][0]);
1004 const double* rightMatrix = &(mat.m_matrix[0][0]);
1005 asm volatile (// First row of leftMatrix.
1006 "mov r3, %[leftMatrix]\n\t"
1007 "vld1.64 { d16-d19 }, [%[leftMatrix], :128]!\n\t"
1008 "vld1.64 { d0-d3}, [%[rightMatrix], :128]!\n\t"
1009 "vmul.f64 d4, d0, d16\n\t"
1010 "vld1.64 { d20-d23 }, [%[leftMatrix], :128]!\n\t"
1011 "vmla.f64 d4, d1, d20\n\t"
1012 "vld1.64 { d24-d27 }, [%[leftMatrix], :128]!\n\t"
1013 "vmla.f64 d4, d2, d24\n\t"
1014 "vld1.64 { d28-d31 }, [%[leftMatrix], :128]!\n\t"
1015 "vmla.f64 d4, d3, d28\n\t"
1016
1017 "vmul.f64 d5, d0, d17\n\t"
1018 "vmla.f64 d5, d1, d21\n\t"
1019 "vmla.f64 d5, d2, d25\n\t"
1020 "vmla.f64 d5, d3, d29\n\t"
1021
1022 "vmul.f64 d6, d0, d18\n\t"
1023 "vmla.f64 d6, d1, d22\n\t"
1024 "vmla.f64 d6, d2, d26\n\t"
1025 "vmla.f64 d6, d3, d30\n\t"
1026
1027 "vmul.f64 d7, d0, d19\n\t"
1028 "vmla.f64 d7, d1, d23\n\t"
1029 "vmla.f64 d7, d2, d27\n\t"
1030 "vmla.f64 d7, d3, d31\n\t"
1031 "vld1.64 { d0-d3}, [%[rightMatrix], :128]!\n\t"
1032 "vst1.64 { d4-d7 }, [r3, :128]!\n\t"
1033
1034 // Second row of leftMatrix.
1035 "vmul.f64 d4, d0, d16\n\t"
1036 "vmla.f64 d4, d1, d20\n\t"
1037 "vmla.f64 d4, d2, d24\n\t"
1038 "vmla.f64 d4, d3, d28\n\t"
1039
1040 "vmul.f64 d5, d0, d17\n\t"
1041 "vmla.f64 d5, d1, d21\n\t"
1042 "vmla.f64 d5, d2, d25\n\t"
1043 "vmla.f64 d5, d3, d29\n\t"
1044
1045 "vmul.f64 d6, d0, d18\n\t"
1046 "vmla.f64 d6, d1, d22\n\t"
1047 "vmla.f64 d6, d2, d26\n\t"
1048 "vmla.f64 d6, d3, d30\n\t"
1049
1050 "vmul.f64 d7, d0, d19\n\t"
1051 "vmla.f64 d7, d1, d23\n\t"
1052 "vmla.f64 d7, d2, d27\n\t"
1053 "vmla.f64 d7, d3, d31\n\t"
1054 "vld1.64 { d0-d3}, [%[rightMatrix], :128]!\n\t"
1055 "vst1.64 { d4-d7 }, [r3, :128]!\n\t"
1056
1057 // Third row of leftMatrix.
1058 "vmul.f64 d4, d0, d16\n\t"
1059 "vmla.f64 d4, d1, d20\n\t"
1060 "vmla.f64 d4, d2, d24\n\t"
1061 "vmla.f64 d4, d3, d28\n\t"
1062
1063 "vmul.f64 d5, d0, d17\n\t"
1064 "vmla.f64 d5, d1, d21\n\t"
1065 "vmla.f64 d5, d2, d25\n\t"
1066 "vmla.f64 d5, d3, d29\n\t"
1067
1068 "vmul.f64 d6, d0, d18\n\t"
1069 "vmla.f64 d6, d1, d22\n\t"
1070 "vmla.f64 d6, d2, d26\n\t"
1071 "vmla.f64 d6, d3, d30\n\t"
1072
1073 "vmul.f64 d7, d0, d19\n\t"
1074 "vmla.f64 d7, d1, d23\n\t"
1075 "vmla.f64 d7, d2, d27\n\t"
1076 "vmla.f64 d7, d3, d31\n\t"
1077 "vld1.64 { d0-d3}, [%[rightMatrix], :128]\n\t"
1078 "vst1.64 { d4-d7 }, [r3, :128]!\n\t"
1079
1080 // Fourth and last row of leftMatrix.
1081 "vmul.f64 d4, d0, d16\n\t"
1082 "vmla.f64 d4, d1, d20\n\t"
1083 "vmla.f64 d4, d2, d24\n\t"
1084 "vmla.f64 d4, d3, d28\n\t"
1085
1086 "vmul.f64 d5, d0, d17\n\t"
1087 "vmla.f64 d5, d1, d21\n\t"
1088 "vmla.f64 d5, d2, d25\n\t"
1089 "vmla.f64 d5, d3, d29\n\t"
1090
1091 "vmul.f64 d6, d0, d18\n\t"
1092 "vmla.f64 d6, d1, d22\n\t"
1093 "vmla.f64 d6, d2, d26\n\t"
1094 "vmla.f64 d6, d3, d30\n\t"
1095
1096 "vmul.f64 d7, d0, d19\n\t"
1097 "vmla.f64 d7, d1, d23\n\t"
1098 "vmla.f64 d7, d2, d27\n\t"
1099 "vmla.f64 d7, d3, d31\n\t"
1100 "vst1.64 { d4-d7 }, [r3, :128]\n\t"
1101 : [leftMatrix]"+r"(leftMatrix), [rightMatrix]"+r"(rightMatrix)
1102 :
1103 : "memory", "r3", "d0", "d1", "d2", "d3", "d4", "d5", "d6", "d7", "d16", "d17", "d18", "d19", "d20", "d21", "d22", "d23", "d24", "d25", "d26", "d27", "d28", "d29", "d30", "d31");
1104 #elif defined(TRANSFORMATION_MATRIX_USE_X86_64_SSE2)
1105 // x86_64 has 16 XMM registers which is enough to do the multiplication fully in registers.
1106 __m128d matrixBlockA = _mm_load_pd(&(m_matrix[0][0]));
1107 __m128d matrixBlockC = _mm_load_pd(&(m_matrix[1][0]));
1108 __m128d matrixBlockE = _mm_load_pd(&(m_matrix[2][0]));
1109 __m128d matrixBlockG = _mm_load_pd(&(m_matrix[3][0]));
1110
1111 // First row.
1112 __m128d otherMatrixFirstParam = _mm_set1_pd(mat.m_matrix[0][0]);
1113 __m128d otherMatrixSecondParam = _mm_set1_pd(mat.m_matrix[0][1]);
1114 __m128d otherMatrixThirdParam = _mm_set1_pd(mat.m_matrix[0][2]);
1115 __m128d otherMatrixFourthParam = _mm_set1_pd(mat.m_matrix[0][3]);
1116
1117 // output00 and output01.
1118 __m128d accumulator = _mm_mul_pd(matrixBlockA, otherMatrixFirstParam);
1119 __m128d temp1 = _mm_mul_pd(matrixBlockC, otherMatrixSecondParam);
1120 __m128d temp2 = _mm_mul_pd(matrixBlockE, otherMatrixThirdParam);
1121 __m128d temp3 = _mm_mul_pd(matrixBlockG, otherMatrixFourthParam);
1122
1123 __m128d matrixBlockB = _mm_load_pd(&(m_matrix[0][2]));
1124 __m128d matrixBlockD = _mm_load_pd(&(m_matrix[1][2]));
1125 __m128d matrixBlockF = _mm_load_pd(&(m_matrix[2][2]));
1126 __m128d matrixBlockH = _mm_load_pd(&(m_matrix[3][2]));
1127
1128 accumulator = _mm_add_pd(accumulator, temp1);
1129 accumulator = _mm_add_pd(accumulator, temp2);
1130 accumulator = _mm_add_pd(accumulator, temp3);
1131 _mm_store_pd(&m_matrix[0][0], accumulator);
1132
1133 // output02 and output03.
1134 accumulator = _mm_mul_pd(matrixBlockB, otherMatrixFirstParam);
1135 temp1 = _mm_mul_pd(matrixBlockD, otherMatrixSecondParam);
1136 temp2 = _mm_mul_pd(matrixBlockF, otherMatrixThirdParam);
1137 temp3 = _mm_mul_pd(matrixBlockH, otherMatrixFourthParam);
1138
1139 accumulator = _mm_add_pd(accumulator, temp1);
1140 accumulator = _mm_add_pd(accumulator, temp2);
1141 accumulator = _mm_add_pd(accumulator, temp3);
1142 _mm_store_pd(&m_matrix[0][2], accumulator);
1143
1144 // Second row.
1145 otherMatrixFirstParam = _mm_set1_pd(mat.m_matrix[1][0]);
1146 otherMatrixSecondParam = _mm_set1_pd(mat.m_matrix[1][1]);
1147 otherMatrixThirdParam = _mm_set1_pd(mat.m_matrix[1][2]);
1148 otherMatrixFourthParam = _mm_set1_pd(mat.m_matrix[1][3]);
1149
1150 // output10 and output11.
1151 accumulator = _mm_mul_pd(matrixBlockA, otherMatrixFirstParam);
1152 temp1 = _mm_mul_pd(matrixBlockC, otherMatrixSecondParam);
1153 temp2 = _mm_mul_pd(matrixBlockE, otherMatrixThirdParam);
1154 temp3 = _mm_mul_pd(matrixBlockG, otherMatrixFourthParam);
1155
1156 accumulator = _mm_add_pd(accumulator, temp1);
1157 accumulator = _mm_add_pd(accumulator, temp2);
1158 accumulator = _mm_add_pd(accumulator, temp3);
1159 _mm_store_pd(&m_matrix[1][0], accumulator);
1160
1161 // output12 and output13.
1162 accumulator = _mm_mul_pd(matrixBlockB, otherMatrixFirstParam);
1163 temp1 = _mm_mul_pd(matrixBlockD, otherMatrixSecondParam);
1164 temp2 = _mm_mul_pd(matrixBlockF, otherMatrixThirdParam);
1165 temp3 = _mm_mul_pd(matrixBlockH, otherMatrixFourthParam);
1166
1167 accumulator = _mm_add_pd(accumulator, temp1);
1168 accumulator = _mm_add_pd(accumulator, temp2);
1169 accumulator = _mm_add_pd(accumulator, temp3);
1170 _mm_store_pd(&m_matrix[1][2], accumulator);
1171
1172 // Third row.
1173 otherMatrixFirstParam = _mm_set1_pd(mat.m_matrix[2][0]);
1174 otherMatrixSecondParam = _mm_set1_pd(mat.m_matrix[2][1]);
1175 otherMatrixThirdParam = _mm_set1_pd(mat.m_matrix[2][2]);
1176 otherMatrixFourthParam = _mm_set1_pd(mat.m_matrix[2][3]);
1177
1178 // output20 and output21.
1179 accumulator = _mm_mul_pd(matrixBlockA, otherMatrixFirstParam);
1180 temp1 = _mm_mul_pd(matrixBlockC, otherMatrixSecondParam);
1181 temp2 = _mm_mul_pd(matrixBlockE, otherMatrixThirdParam);
1182 temp3 = _mm_mul_pd(matrixBlockG, otherMatrixFourthParam);
1183
1184 accumulator = _mm_add_pd(accumulator, temp1);
1185 accumulator = _mm_add_pd(accumulator, temp2);
1186 accumulator = _mm_add_pd(accumulator, temp3);
1187 _mm_store_pd(&m_matrix[2][0], accumulator);
1188
1189 // output22 and output23.
1190 accumulator = _mm_mul_pd(matrixBlockB, otherMatrixFirstParam);
1191 temp1 = _mm_mul_pd(matrixBlockD, otherMatrixSecondParam);
1192 temp2 = _mm_mul_pd(matrixBlockF, otherMatrixThirdParam);
1193 temp3 = _mm_mul_pd(matrixBlockH, otherMatrixFourthParam);
1194
1195 accumulator = _mm_add_pd(accumulator, temp1);
1196 accumulator = _mm_add_pd(accumulator, temp2);
1197 accumulator = _mm_add_pd(accumulator, temp3);
1198 _mm_store_pd(&m_matrix[2][2], accumulator);
1199
1200 // Fourth row.
1201 otherMatrixFirstParam = _mm_set1_pd(mat.m_matrix[3][0]);
1202 otherMatrixSecondParam = _mm_set1_pd(mat.m_matrix[3][1]);
1203 otherMatrixThirdParam = _mm_set1_pd(mat.m_matrix[3][2]);
1204 otherMatrixFourthParam = _mm_set1_pd(mat.m_matrix[3][3]);
1205
1206 // output30 and output31.
1207 accumulator = _mm_mul_pd(matrixBlockA, otherMatrixFirstParam);
1208 temp1 = _mm_mul_pd(matrixBlockC, otherMatrixSecondParam);
1209 temp2 = _mm_mul_pd(matrixBlockE, otherMatrixThirdParam);
1210 temp3 = _mm_mul_pd(matrixBlockG, otherMatrixFourthParam);
1211
1212 accumulator = _mm_add_pd(accumulator, temp1);
1213 accumulator = _mm_add_pd(accumulator, temp2);
1214 accumulator = _mm_add_pd(accumulator, temp3);
1215 _mm_store_pd(&m_matrix[3][0], accumulator);
1216
1217 // output32 and output33.
1218 accumulator = _mm_mul_pd(matrixBlockB, otherMatrixFirstParam);
1219 temp1 = _mm_mul_pd(matrixBlockD, otherMatrixSecondParam);
1220 temp2 = _mm_mul_pd(matrixBlockF, otherMatrixThirdParam);
1221 temp3 = _mm_mul_pd(matrixBlockH, otherMatrixFourthParam);
1222
1223 accumulator = _mm_add_pd(accumulator, temp1);
1224 accumulator = _mm_add_pd(accumulator, temp2);
1225 accumulator = _mm_add_pd(accumulator, temp3);
1226 _mm_store_pd(&m_matrix[3][2], accumulator);
1227 #else
1228 Matrix4 tmp;
1229
1230 tmp[0][0] = (mat.m_matrix[0][0] * m_matrix[0][0] + mat.m_matrix[0][1] * m_matrix[1][0]
1231 + mat.m_matrix[0][2] * m_matrix[2][0] + mat.m_matrix[0][3] * m_matrix[3][0]);
1232 tmp[0][1] = (mat.m_matrix[0][0] * m_matrix[0][1] + mat.m_matrix[0][1] * m_matrix[1][1]
1233 + mat.m_matrix[0][2] * m_matrix[2][1] + mat.m_matrix[0][3] * m_matrix[3][1]);
1234 tmp[0][2] = (mat.m_matrix[0][0] * m_matrix[0][2] + mat.m_matrix[0][1] * m_matrix[1][2]
1235 + mat.m_matrix[0][2] * m_matrix[2][2] + mat.m_matrix[0][3] * m_matrix[3][2]);
1236 tmp[0][3] = (mat.m_matrix[0][0] * m_matrix[0][3] + mat.m_matrix[0][1] * m_matrix[1][3]
1237 + mat.m_matrix[0][2] * m_matrix[2][3] + mat.m_matrix[0][3] * m_matrix[3][3]);
1238
1239 tmp[1][0] = (mat.m_matrix[1][0] * m_matrix[0][0] + mat.m_matrix[1][1] * m_matrix[1][0]
1240 + mat.m_matrix[1][2] * m_matrix[2][0] + mat.m_matrix[1][3] * m_matrix[3][0]);
1241 tmp[1][1] = (mat.m_matrix[1][0] * m_matrix[0][1] + mat.m_matrix[1][1] * m_matrix[1][1]
1242 + mat.m_matrix[1][2] * m_matrix[2][1] + mat.m_matrix[1][3] * m_matrix[3][1]);
1243 tmp[1][2] = (mat.m_matrix[1][0] * m_matrix[0][2] + mat.m_matrix[1][1] * m_matrix[1][2]
1244 + mat.m_matrix[1][2] * m_matrix[2][2] + mat.m_matrix[1][3] * m_matrix[3][2]);
1245 tmp[1][3] = (mat.m_matrix[1][0] * m_matrix[0][3] + mat.m_matrix[1][1] * m_matrix[1][3]
1246 + mat.m_matrix[1][2] * m_matrix[2][3] + mat.m_matrix[1][3] * m_matrix[3][3]);
1247
1248 tmp[2][0] = (mat.m_matrix[2][0] * m_matrix[0][0] + mat.m_matrix[2][1] * m_matrix[1][0]
1249 + mat.m_matrix[2][2] * m_matrix[2][0] + mat.m_matrix[2][3] * m_matrix[3][0]);
1250 tmp[2][1] = (mat.m_matrix[2][0] * m_matrix[0][1] + mat.m_matrix[2][1] * m_matrix[1][1]
1251 + mat.m_matrix[2][2] * m_matrix[2][1] + mat.m_matrix[2][3] * m_matrix[3][1]);
1252 tmp[2][2] = (mat.m_matrix[2][0] * m_matrix[0][2] + mat.m_matrix[2][1] * m_matrix[1][2]
1253 + mat.m_matrix[2][2] * m_matrix[2][2] + mat.m_matrix[2][3] * m_matrix[3][2]);
1254 tmp[2][3] = (mat.m_matrix[2][0] * m_matrix[0][3] + mat.m_matrix[2][1] * m_matrix[1][3]
1255 + mat.m_matrix[2][2] * m_matrix[2][3] + mat.m_matrix[2][3] * m_matrix[3][3]);
1256
1257 tmp[3][0] = (mat.m_matrix[3][0] * m_matrix[0][0] + mat.m_matrix[3][1] * m_matrix[1][0]
1258 + mat.m_matrix[3][2] * m_matrix[2][0] + mat.m_matrix[3][3] * m_matrix[3][0]);
1259 tmp[3][1] = (mat.m_matrix[3][0] * m_matrix[0][1] + mat.m_matrix[3][1] * m_matrix[1][1]
1260 + mat.m_matrix[3][2] * m_matrix[2][1] + mat.m_matrix[3][3] * m_matrix[3][1]);
1261 tmp[3][2] = (mat.m_matrix[3][0] * m_matrix[0][2] + mat.m_matrix[3][1] * m_matrix[1][2]
1262 + mat.m_matrix[3][2] * m_matrix[2][2] + mat.m_matrix[3][3] * m_matrix[3][2]);
1263 tmp[3][3] = (mat.m_matrix[3][0] * m_matrix[0][3] + mat.m_matrix[3][1] * m_matrix[1][3]
1264 + mat.m_matrix[3][2] * m_matrix[2][3] + mat.m_matrix[3][3] * m_matrix[3][3]);
1265
1266 setMatrix(tmp);
1267 #endif
1268 return *this;
1269 }
1270
multVecMatrix(double x,double y,double & resultX,double & resultY) const1271 void TransformationMatrix::multVecMatrix(double x, double y, double& resultX, double& resultY) const
1272 {
1273 resultX = m_matrix[3][0] + x * m_matrix[0][0] + y * m_matrix[1][0];
1274 resultY = m_matrix[3][1] + x * m_matrix[0][1] + y * m_matrix[1][1];
1275 double w = m_matrix[3][3] + x * m_matrix[0][3] + y * m_matrix[1][3];
1276 if (w != 1 && w != 0) {
1277 resultX /= w;
1278 resultY /= w;
1279 }
1280 }
1281
multVecMatrix(double x,double y,double z,double & resultX,double & resultY,double & resultZ) const1282 void TransformationMatrix::multVecMatrix(double x, double y, double z, double& resultX, double& resultY, double& resultZ) const
1283 {
1284 resultX = m_matrix[3][0] + x * m_matrix[0][0] + y * m_matrix[1][0] + z * m_matrix[2][0];
1285 resultY = m_matrix[3][1] + x * m_matrix[0][1] + y * m_matrix[1][1] + z * m_matrix[2][1];
1286 resultZ = m_matrix[3][2] + x * m_matrix[0][2] + y * m_matrix[1][2] + z * m_matrix[2][2];
1287 double w = m_matrix[3][3] + x * m_matrix[0][3] + y * m_matrix[1][3] + z * m_matrix[2][3];
1288 if (w != 1 && w != 0) {
1289 resultX /= w;
1290 resultY /= w;
1291 resultZ /= w;
1292 }
1293 }
1294
isInvertible() const1295 bool TransformationMatrix::isInvertible() const
1296 {
1297 if (isIdentityOrTranslation())
1298 return true;
1299
1300 double det = WebCore::determinant4x4(m_matrix);
1301
1302 if (fabs(det) < SMALL_NUMBER)
1303 return false;
1304
1305 return true;
1306 }
1307
inverse() const1308 TransformationMatrix TransformationMatrix::inverse() const
1309 {
1310 if (isIdentityOrTranslation()) {
1311 // identity matrix
1312 if (m_matrix[3][0] == 0 && m_matrix[3][1] == 0 && m_matrix[3][2] == 0)
1313 return TransformationMatrix();
1314
1315 // translation
1316 return TransformationMatrix(1, 0, 0, 0,
1317 0, 1, 0, 0,
1318 0, 0, 1, 0,
1319 -m_matrix[3][0], -m_matrix[3][1], -m_matrix[3][2], 1);
1320 }
1321
1322 TransformationMatrix invMat;
1323 bool inverted = WebCore::inverse(m_matrix, invMat.m_matrix);
1324 if (!inverted)
1325 return TransformationMatrix();
1326
1327 return invMat;
1328 }
1329
makeAffine()1330 void TransformationMatrix::makeAffine()
1331 {
1332 m_matrix[0][2] = 0;
1333 m_matrix[0][3] = 0;
1334
1335 m_matrix[1][2] = 0;
1336 m_matrix[1][3] = 0;
1337
1338 m_matrix[2][0] = 0;
1339 m_matrix[2][1] = 0;
1340 m_matrix[2][2] = 1;
1341 m_matrix[2][3] = 0;
1342
1343 m_matrix[3][2] = 0;
1344 m_matrix[3][3] = 1;
1345 }
1346
toAffineTransform() const1347 AffineTransform TransformationMatrix::toAffineTransform() const
1348 {
1349 return AffineTransform(m_matrix[0][0], m_matrix[0][1], m_matrix[1][0],
1350 m_matrix[1][1], m_matrix[3][0], m_matrix[3][1]);
1351 }
1352
blendFloat(double & from,double to,double progress)1353 static inline void blendFloat(double& from, double to, double progress)
1354 {
1355 if (from != to)
1356 from = from + (to - from) * progress;
1357 }
1358
blend(const TransformationMatrix & from,double progress)1359 void TransformationMatrix::blend(const TransformationMatrix& from, double progress)
1360 {
1361 if (from.isIdentity() && isIdentity())
1362 return;
1363
1364 // decompose
1365 DecomposedType fromDecomp;
1366 DecomposedType toDecomp;
1367 if (!from.decompose(fromDecomp) || !decompose(toDecomp)) {
1368 if (progress < 0.5)
1369 *this = from;
1370 return;
1371 }
1372
1373 // interpolate
1374 blendFloat(fromDecomp.scaleX, toDecomp.scaleX, progress);
1375 blendFloat(fromDecomp.scaleY, toDecomp.scaleY, progress);
1376 blendFloat(fromDecomp.scaleZ, toDecomp.scaleZ, progress);
1377 blendFloat(fromDecomp.skewXY, toDecomp.skewXY, progress);
1378 blendFloat(fromDecomp.skewXZ, toDecomp.skewXZ, progress);
1379 blendFloat(fromDecomp.skewYZ, toDecomp.skewYZ, progress);
1380 blendFloat(fromDecomp.translateX, toDecomp.translateX, progress);
1381 blendFloat(fromDecomp.translateY, toDecomp.translateY, progress);
1382 blendFloat(fromDecomp.translateZ, toDecomp.translateZ, progress);
1383 blendFloat(fromDecomp.perspectiveX, toDecomp.perspectiveX, progress);
1384 blendFloat(fromDecomp.perspectiveY, toDecomp.perspectiveY, progress);
1385 blendFloat(fromDecomp.perspectiveZ, toDecomp.perspectiveZ, progress);
1386 blendFloat(fromDecomp.perspectiveW, toDecomp.perspectiveW, progress);
1387
1388 slerp(&fromDecomp.quaternionX, &toDecomp.quaternionX, progress);
1389
1390 // recompose
1391 recompose(fromDecomp);
1392 }
1393
decompose(DecomposedType & decomp) const1394 bool TransformationMatrix::decompose(DecomposedType& decomp) const
1395 {
1396 if (isIdentity()) {
1397 memset(&decomp, 0, sizeof(decomp));
1398 decomp.perspectiveW = 1;
1399 decomp.scaleX = 1;
1400 decomp.scaleY = 1;
1401 decomp.scaleZ = 1;
1402 }
1403
1404 if (!WebCore::decompose(m_matrix, decomp))
1405 return false;
1406 return true;
1407 }
1408
recompose(const DecomposedType & decomp)1409 void TransformationMatrix::recompose(const DecomposedType& decomp)
1410 {
1411 makeIdentity();
1412
1413 // first apply perspective
1414 m_matrix[0][3] = decomp.perspectiveX;
1415 m_matrix[1][3] = decomp.perspectiveY;
1416 m_matrix[2][3] = decomp.perspectiveZ;
1417 m_matrix[3][3] = decomp.perspectiveW;
1418
1419 // now translate
1420 translate3d(decomp.translateX, decomp.translateY, decomp.translateZ);
1421
1422 // apply rotation
1423 double xx = decomp.quaternionX * decomp.quaternionX;
1424 double xy = decomp.quaternionX * decomp.quaternionY;
1425 double xz = decomp.quaternionX * decomp.quaternionZ;
1426 double xw = decomp.quaternionX * decomp.quaternionW;
1427 double yy = decomp.quaternionY * decomp.quaternionY;
1428 double yz = decomp.quaternionY * decomp.quaternionZ;
1429 double yw = decomp.quaternionY * decomp.quaternionW;
1430 double zz = decomp.quaternionZ * decomp.quaternionZ;
1431 double zw = decomp.quaternionZ * decomp.quaternionW;
1432
1433 // Construct a composite rotation matrix from the quaternion values
1434 TransformationMatrix rotationMatrix(1 - 2 * (yy + zz), 2 * (xy - zw), 2 * (xz + yw), 0,
1435 2 * (xy + zw), 1 - 2 * (xx + zz), 2 * (yz - xw), 0,
1436 2 * (xz - yw), 2 * (yz + xw), 1 - 2 * (xx + yy), 0,
1437 0, 0, 0, 1);
1438
1439 multiply(rotationMatrix);
1440
1441 // now apply skew
1442 if (decomp.skewYZ) {
1443 TransformationMatrix tmp;
1444 tmp.setM32(decomp.skewYZ);
1445 multiply(tmp);
1446 }
1447
1448 if (decomp.skewXZ) {
1449 TransformationMatrix tmp;
1450 tmp.setM31(decomp.skewXZ);
1451 multiply(tmp);
1452 }
1453
1454 if (decomp.skewXY) {
1455 TransformationMatrix tmp;
1456 tmp.setM21(decomp.skewXY);
1457 multiply(tmp);
1458 }
1459
1460 // finally, apply scale
1461 scale3d(decomp.scaleX, decomp.scaleY, decomp.scaleZ);
1462 }
1463
isIntegerTranslation() const1464 bool TransformationMatrix::isIntegerTranslation() const
1465 {
1466 if (!isIdentityOrTranslation())
1467 return false;
1468
1469 // Check for translate Z.
1470 if (m_matrix[3][2])
1471 return false;
1472
1473 // Check for non-integer translate X/Y.
1474 if (static_cast<int>(m_matrix[3][0]) != m_matrix[3][0] || static_cast<int>(m_matrix[3][1]) != m_matrix[3][1])
1475 return false;
1476
1477 return true;
1478 }
1479
to2dTransform() const1480 TransformationMatrix TransformationMatrix::to2dTransform() const
1481 {
1482 return TransformationMatrix(m_matrix[0][0], m_matrix[0][1], 0, m_matrix[0][3],
1483 m_matrix[1][0], m_matrix[1][1], 0, m_matrix[1][3],
1484 0, 0, 1, 0,
1485 m_matrix[3][0], m_matrix[3][1], 0, m_matrix[3][3]);
1486 }
1487
toColumnMajorFloatArray(FloatMatrix4 & result) const1488 void TransformationMatrix::toColumnMajorFloatArray(FloatMatrix4& result) const
1489 {
1490 result[0] = m11();
1491 result[1] = m12();
1492 result[2] = m13();
1493 result[3] = m14();
1494 result[4] = m21();
1495 result[5] = m22();
1496 result[6] = m23();
1497 result[7] = m24();
1498 result[8] = m31();
1499 result[9] = m32();
1500 result[10] = m33();
1501 result[11] = m34();
1502 result[12] = m41();
1503 result[13] = m42();
1504 result[14] = m43();
1505 result[15] = m44();
1506 }
1507
isBackFaceVisible() const1508 bool TransformationMatrix::isBackFaceVisible() const
1509 {
1510 // Back-face visibility is determined by transforming the normal vector (0, 0, 1) and
1511 // checking the sign of the resulting z component. However, normals cannot be
1512 // transformed by the original matrix, they require being transformed by the
1513 // inverse-transpose.
1514 //
1515 // Since we know we will be using (0, 0, 1), and we only care about the z-component of
1516 // the transformed normal, then we only need the m33() element of the
1517 // inverse-transpose. Therefore we do not need the transpose.
1518 //
1519 // Additionally, if we only need the m33() element, we do not need to compute a full
1520 // inverse. Instead, knowing the inverse of a matrix is adjoint(matrix) / determinant,
1521 // we can simply compute the m33() of the adjoint (adjugate) matrix, without computing
1522 // the full adjoint.
1523
1524 double determinant = WebCore::determinant4x4(m_matrix);
1525
1526 // If the matrix is not invertible, then we assume its backface is not visible.
1527 if (fabs(determinant) < SMALL_NUMBER)
1528 return false;
1529
1530 double cofactor33 = determinant3x3(m11(), m12(), m14(), m21(), m22(), m24(), m41(), m42(), m44());
1531 double zComponentOfTransformedNormal = cofactor33 / determinant;
1532
1533 return zComponentOfTransformedNormal < 0;
1534 }
1535
toSkMatrix44(const TransformationMatrix & matrix)1536 SkMatrix44 TransformationMatrix::toSkMatrix44(const TransformationMatrix& matrix)
1537 {
1538 SkMatrix44 ret(SkMatrix44::kUninitialized_Constructor);
1539 ret.setDouble(0, 0, matrix.m11());
1540 ret.setDouble(0, 1, matrix.m21());
1541 ret.setDouble(0, 2, matrix.m31());
1542 ret.setDouble(0, 3, matrix.m41());
1543 ret.setDouble(1, 0, matrix.m12());
1544 ret.setDouble(1, 1, matrix.m22());
1545 ret.setDouble(1, 2, matrix.m32());
1546 ret.setDouble(1, 3, matrix.m42());
1547 ret.setDouble(2, 0, matrix.m13());
1548 ret.setDouble(2, 1, matrix.m23());
1549 ret.setDouble(2, 2, matrix.m33());
1550 ret.setDouble(2, 3, matrix.m43());
1551 ret.setDouble(3, 0, matrix.m14());
1552 ret.setDouble(3, 1, matrix.m24());
1553 ret.setDouble(3, 2, matrix.m34());
1554 ret.setDouble(3, 3, matrix.m44());
1555 return ret;
1556 }
1557
1558 }
1559