Home | History | Annotate | Download | only in transforms
      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 }
   1559