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
45namespace 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
76typedef double Vector4[4];
77typedef double Vector3[3];
78
79const 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
93static 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
107static 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
118static 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
165static 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
212static 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
240static 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
248static 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
260static double v3Length(Vector3 a)
261{
262    return std::sqrt((a[0] * a[0]) + (a[1] * a[1]) + (a[2] * a[2]));
263}
264
265static 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
276static 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)
283static 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 */
291static 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
298static 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
473static 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
521TransformationMatrix::TransformationMatrix(const AffineTransform& t)
522{
523    setMatrix(t.a(), t.b(), t.c(), t.d(), t.e(), t.f());
524}
525
526TransformationMatrix& TransformationMatrix::scale(double s)
527{
528    return scaleNonUniform(s, s);
529}
530
531TransformationMatrix& TransformationMatrix::rotateFromVector(double x, double y)
532{
533    return rotate(rad2deg(atan2(y, x)));
534}
535
536TransformationMatrix& TransformationMatrix::flipX()
537{
538    return scaleNonUniform(-1.0, 1.0);
539}
540
541TransformationMatrix& TransformationMatrix::flipY()
542{
543    return scaleNonUniform(1.0, -1.0);
544}
545
546FloatPoint 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
595FloatQuad 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
620static 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
626LayoutRect 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
648void 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
670FloatPoint 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
678FloatPoint3D 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
688IntRect TransformationMatrix::mapRect(const IntRect &rect) const
689{
690    return enclosingIntRect(mapRect(FloatRect(rect)));
691}
692
693LayoutRect TransformationMatrix::mapRect(const LayoutRect& r) const
694{
695    return enclosingLayoutRect(mapRect(FloatRect(r)));
696}
697
698FloatRect 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
718FloatQuad 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
734TransformationMatrix& 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
748TransformationMatrix& 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
759TransformationMatrix& 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
846TransformationMatrix& 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
913TransformationMatrix& 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
922TransformationMatrix& 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
931TransformationMatrix& 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
950TransformationMatrix& 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
963TransformationMatrix& 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
977TransformationMatrix& 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
987TransformationMatrix 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.
998TransformationMatrix& 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
1330void 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
1341void 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
1354bool 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
1367TransformationMatrix 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
1389void 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
1406AffineTransform 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
1412static inline void blendFloat(double& from, double to, double progress)
1413{
1414    if (from != to)
1415        from = from + (to - from) * progress;
1416}
1417
1418void 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
1453bool 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
1468void 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
1523bool 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
1539TransformationMatrix 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
1547void 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
1567SkMatrix44 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}
1590