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