src/third_party/WebKit/Source/platform/transforms/TransformationMatrix.cpp

   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
  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
 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
 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
 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
 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
 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
 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
 262 static double v3Length(Vector3 a)
 263 {
 264     return sqrt((a[0] * a[0]) + (a[1] * a[1]) + (a[2] * a[2]));
 265 }
 266
 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
 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)
 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 */
 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
 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
 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
 523 TransformationMatrix::TransformationMatrix(const AffineTransform& t)
 524 {
 525     setMatrix(t.a(), t.b(), t.c(), t.d(), t.e(), t.f());
 526 }
 527
 528 TransformationMatrix& TransformationMatrix::scale(double s)
 529 {
 530     return scaleNonUniform(s, s);
 531 }
 532
 533 TransformationMatrix& TransformationMatrix::rotateFromVector(double x, double y)
 534 {
 535     return rotate(rad2deg(atan2(y, x)));
 536 }
 537
 538 TransformationMatrix& TransformationMatrix::flipX()
 539 {
 540     return scaleNonUniform(-1.0, 1.0);
 541 }
 542
 543 TransformationMatrix& TransformationMatrix::flipY()
 544 {
 545     return scaleNonUniform(1.0, -1.0);
 546 }
 547
 548 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
 597 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
 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
 628 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
 650 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
 672 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
 680 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
 690 IntRect TransformationMatrix::mapRect(const IntRect &rect) const
 691 {
 692     return enclosingIntRect(mapRect(FloatRect(rect)));
 693 }
 694
 695 LayoutRect TransformationMatrix::mapRect(const LayoutRect& r) const
 696 {
 697     return enclosingLayoutRect(mapRect(FloatRect(r)));
 698 }
 699
 700 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
 720 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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
 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.
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
1271 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
1282 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
1295 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
1308 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
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
1347 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
1353 static inline void blendFloat(double& from, double to, double progress)
1354 {
1355     if (from != to)
1356         from = from + (to - from) * progress;
1357 }
1358
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
1394 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
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
1464 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
1480 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
1488 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
1508 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
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 }