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