• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*
2  * Copyright (C) 2005, 2006 Apple Computer, Inc.  All rights reserved.
3  * Copyright (C) 2009 Torch Mobile, Inc.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions
7  * are met:
8  * 1. Redistributions of source code must retain the above copyright
9  *    notice, this list of conditions and the following disclaimer.
10  * 2. Redistributions in binary form must reproduce the above copyright
11  *    notice, this list of conditions and the following disclaimer in the
12  *    documentation and/or other materials provided with the distribution.
13  *
14  * THIS SOFTWARE IS PROVIDED BY APPLE COMPUTER, INC. ``AS IS'' AND ANY
15  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
17  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL APPLE COMPUTER, INC. OR
18  * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
19  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
20  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
21  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
22  * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
24  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25  */
26 
27 #include "config.h"
28 #include "TransformationMatrix.h"
29 
30 #include "AffineTransform.h"
31 #include "FloatPoint3D.h"
32 #include "FloatRect.h"
33 #include "FloatQuad.h"
34 #include "IntRect.h"
35 
36 #include <wtf/Assertions.h>
37 #include <wtf/MathExtras.h>
38 
39 namespace WebCore {
40 
41 //
42 // Supporting Math Functions
43 //
44 // This is a set of function from various places (attributed inline) to do things like
45 // inversion and decomposition of a 4x4 matrix. They are used throughout the code
46 //
47 
48 //
49 // Adapted from Matrix Inversion by Richard Carling, Graphics Gems <http://tog.acm.org/GraphicsGems/index.html>.
50 
51 // EULA: The Graphics Gems code is copyright-protected. In other words, you cannot claim the text of the code
52 // as your own and resell it. Using the code is permitted in any program, product, or library, non-commercial
53 // or commercial. Giving credit is not required, though is a nice gesture. The code comes as-is, and if there
54 // are any flaws or problems with any Gems code, nobody involved with Gems - authors, editors, publishers, or
55 // webmasters - are to be held responsible. Basically, don't be a jerk, and remember that anything free comes
56 // with no guarantee.
57 
58 // A clarification about the storage of matrix elements
59 //
60 // This class uses a 2 dimensional array internally to store the elements of the matrix.  The first index into
61 // the array refers to the column that the element lies in; the second index refers to the row.
62 //
63 // In other words, this is the layout of the matrix:
64 //
65 // | m_matrix[0][0] m_matrix[1][0] m_matrix[2][0] m_matrix[3][0] |
66 // | m_matrix[0][1] m_matrix[1][1] m_matrix[2][1] m_matrix[3][1] |
67 // | m_matrix[0][2] m_matrix[1][2] m_matrix[2][2] m_matrix[3][2] |
68 // | m_matrix[0][3] m_matrix[1][3] m_matrix[2][3] m_matrix[3][3] |
69 
70 typedef double Vector4[4];
71 typedef double Vector3[3];
72 
73 const double SMALL_NUMBER = 1.e-8;
74 
75 // inverse(original_matrix, inverse_matrix)
76 //
77 // calculate the inverse of a 4x4 matrix
78 //
79 // -1
80 // A  = ___1__ adjoint A
81 //       det A
82 
83 //  double = determinant2x2(double a, double b, double c, double d)
84 //
85 //  calculate the determinant of a 2x2 matrix.
86 
determinant2x2(double a,double b,double c,double d)87 static double determinant2x2(double a, double b, double c, double d)
88 {
89     return a * d - b * c;
90 }
91 
92 //  double = determinant3x3(a1, a2, a3, b1, b2, b3, c1, c2, c3)
93 //
94 //  Calculate the determinant of a 3x3 matrix
95 //  in the form
96 //
97 //      | a1,  b1,  c1 |
98 //      | a2,  b2,  c2 |
99 //      | a3,  b3,  c3 |
100 
determinant3x3(double a1,double a2,double a3,double b1,double b2,double b3,double c1,double c2,double c3)101 static double determinant3x3(double a1, double a2, double a3, double b1, double b2, double b3, double c1, double c2, double c3)
102 {
103     return a1 * determinant2x2(b2, b3, c2, c3)
104          - b1 * determinant2x2(a2, a3, c2, c3)
105          + c1 * determinant2x2(a2, a3, b2, b3);
106 }
107 
108 //  double = determinant4x4(matrix)
109 //
110 //  calculate the determinant of a 4x4 matrix.
111 
determinant4x4(const TransformationMatrix::Matrix4 & m)112 static double determinant4x4(const TransformationMatrix::Matrix4& m)
113 {
114     // Assign to individual variable names to aid selecting
115     // correct elements
116 
117     double a1 = m[0][0];
118     double b1 = m[0][1];
119     double c1 = m[0][2];
120     double d1 = m[0][3];
121 
122     double a2 = m[1][0];
123     double b2 = m[1][1];
124     double c2 = m[1][2];
125     double d2 = m[1][3];
126 
127     double a3 = m[2][0];
128     double b3 = m[2][1];
129     double c3 = m[2][2];
130     double d3 = m[2][3];
131 
132     double a4 = m[3][0];
133     double b4 = m[3][1];
134     double c4 = m[3][2];
135     double d4 = m[3][3];
136 
137     return a1 * determinant3x3(b2, b3, b4, c2, c3, c4, d2, d3, d4)
138          - b1 * determinant3x3(a2, a3, a4, c2, c3, c4, d2, d3, d4)
139          + c1 * determinant3x3(a2, a3, a4, b2, b3, b4, d2, d3, d4)
140          - d1 * determinant3x3(a2, a3, a4, b2, b3, b4, c2, c3, c4);
141 }
142 
143 // adjoint( original_matrix, inverse_matrix )
144 //
145 //   calculate the adjoint of a 4x4 matrix
146 //
147 //    Let  a   denote the minor determinant of matrix A obtained by
148 //         ij
149 //
150 //    deleting the ith row and jth column from A.
151 //
152 //                  i+j
153 //   Let  b   = (-1)    a
154 //        ij            ji
155 //
156 //  The matrix B = (b  ) is the adjoint of A
157 //                   ij
158 
adjoint(const TransformationMatrix::Matrix4 & matrix,TransformationMatrix::Matrix4 & result)159 static void adjoint(const TransformationMatrix::Matrix4& matrix, TransformationMatrix::Matrix4& result)
160 {
161     // Assign to individual variable names to aid
162     // selecting correct values
163     double a1 = matrix[0][0];
164     double b1 = matrix[0][1];
165     double c1 = matrix[0][2];
166     double d1 = matrix[0][3];
167 
168     double a2 = matrix[1][0];
169     double b2 = matrix[1][1];
170     double c2 = matrix[1][2];
171     double d2 = matrix[1][3];
172 
173     double a3 = matrix[2][0];
174     double b3 = matrix[2][1];
175     double c3 = matrix[2][2];
176     double d3 = matrix[2][3];
177 
178     double a4 = matrix[3][0];
179     double b4 = matrix[3][1];
180     double c4 = matrix[3][2];
181     double d4 = matrix[3][3];
182 
183     // Row column labeling reversed since we transpose rows & columns
184     result[0][0]  =   determinant3x3(b2, b3, b4, c2, c3, c4, d2, d3, d4);
185     result[1][0]  = - determinant3x3(a2, a3, a4, c2, c3, c4, d2, d3, d4);
186     result[2][0]  =   determinant3x3(a2, a3, a4, b2, b3, b4, d2, d3, d4);
187     result[3][0]  = - determinant3x3(a2, a3, a4, b2, b3, b4, c2, c3, c4);
188 
189     result[0][1]  = - determinant3x3(b1, b3, b4, c1, c3, c4, d1, d3, d4);
190     result[1][1]  =   determinant3x3(a1, a3, a4, c1, c3, c4, d1, d3, d4);
191     result[2][1]  = - determinant3x3(a1, a3, a4, b1, b3, b4, d1, d3, d4);
192     result[3][1]  =   determinant3x3(a1, a3, a4, b1, b3, b4, c1, c3, c4);
193 
194     result[0][2]  =   determinant3x3(b1, b2, b4, c1, c2, c4, d1, d2, d4);
195     result[1][2]  = - determinant3x3(a1, a2, a4, c1, c2, c4, d1, d2, d4);
196     result[2][2]  =   determinant3x3(a1, a2, a4, b1, b2, b4, d1, d2, d4);
197     result[3][2]  = - determinant3x3(a1, a2, a4, b1, b2, b4, c1, c2, c4);
198 
199     result[0][3]  = - determinant3x3(b1, b2, b3, c1, c2, c3, d1, d2, d3);
200     result[1][3]  =   determinant3x3(a1, a2, a3, c1, c2, c3, d1, d2, d3);
201     result[2][3]  = - determinant3x3(a1, a2, a3, b1, b2, b3, d1, d2, d3);
202     result[3][3]  =   determinant3x3(a1, a2, a3, b1, b2, b3, c1, c2, c3);
203 }
204 
205 // Returns false if the matrix is not invertible
inverse(const TransformationMatrix::Matrix4 & matrix,TransformationMatrix::Matrix4 & result)206 static bool inverse(const TransformationMatrix::Matrix4& matrix, TransformationMatrix::Matrix4& result)
207 {
208     // Calculate the adjoint matrix
209     adjoint(matrix, result);
210 
211     // Calculate the 4x4 determinant
212     // If the determinant is zero,
213     // then the inverse matrix is not unique.
214     double det = determinant4x4(matrix);
215 
216     if (fabs(det) < SMALL_NUMBER)
217         return false;
218 
219     // Scale the adjoint matrix to get the inverse
220 
221     for (int i = 0; i < 4; i++)
222         for (int j = 0; j < 4; j++)
223             result[i][j] = result[i][j] / det;
224 
225     return true;
226 }
227 
228 // End of code adapted from Matrix Inversion by Richard Carling
229 
230 // Perform a decomposition on the passed matrix, return false if unsuccessful
231 // From Graphics Gems: unmatrix.c
232 
233 // Transpose rotation portion of matrix a, return b
transposeMatrix4(const TransformationMatrix::Matrix4 & a,TransformationMatrix::Matrix4 & b)234 static void transposeMatrix4(const TransformationMatrix::Matrix4& a, TransformationMatrix::Matrix4& b)
235 {
236     for (int i = 0; i < 4; i++)
237         for (int j = 0; j < 4; j++)
238             b[i][j] = a[j][i];
239 }
240 
241 // Multiply a homogeneous point by a matrix and return the transformed point
v4MulPointByMatrix(const Vector4 p,const TransformationMatrix::Matrix4 & m,Vector4 result)242 static void v4MulPointByMatrix(const Vector4 p, const TransformationMatrix::Matrix4& m, Vector4 result)
243 {
244     result[0] = (p[0] * m[0][0]) + (p[1] * m[1][0]) +
245                 (p[2] * m[2][0]) + (p[3] * m[3][0]);
246     result[1] = (p[0] * m[0][1]) + (p[1] * m[1][1]) +
247                 (p[2] * m[2][1]) + (p[3] * m[3][1]);
248     result[2] = (p[0] * m[0][2]) + (p[1] * m[1][2]) +
249                 (p[2] * m[2][2]) + (p[3] * m[3][2]);
250     result[3] = (p[0] * m[0][3]) + (p[1] * m[1][3]) +
251                 (p[2] * m[2][3]) + (p[3] * m[3][3]);
252 }
253 
v3Length(Vector3 a)254 static double v3Length(Vector3 a)
255 {
256     return sqrt((a[0] * a[0]) + (a[1] * a[1]) + (a[2] * a[2]));
257 }
258 
v3Scale(Vector3 v,double desiredLength)259 static void v3Scale(Vector3 v, double desiredLength)
260 {
261     double len = v3Length(v);
262     if (len != 0) {
263         double l = desiredLength / len;
264         v[0] *= l;
265         v[1] *= l;
266         v[2] *= l;
267     }
268 }
269 
v3Dot(const Vector3 a,const Vector3 b)270 static double v3Dot(const Vector3 a, const Vector3 b)
271 {
272     return (a[0] * b[0]) + (a[1] * b[1]) + (a[2] * b[2]);
273 }
274 
275 // Make a linear combination of two vectors and return the result.
276 // result = (a * ascl) + (b * bscl)
v3Combine(const Vector3 a,const Vector3 b,Vector3 result,double ascl,double bscl)277 static void v3Combine(const Vector3 a, const Vector3 b, Vector3 result, double ascl, double bscl)
278 {
279     result[0] = (ascl * a[0]) + (bscl * b[0]);
280     result[1] = (ascl * a[1]) + (bscl * b[1]);
281     result[2] = (ascl * a[2]) + (bscl * b[2]);
282 }
283 
284 // Return the cross product result = a cross b */
v3Cross(const Vector3 a,const Vector3 b,Vector3 result)285 static void v3Cross(const Vector3 a, const Vector3 b, Vector3 result)
286 {
287     result[0] = (a[1] * b[2]) - (a[2] * b[1]);
288     result[1] = (a[2] * b[0]) - (a[0] * b[2]);
289     result[2] = (a[0] * b[1]) - (a[1] * b[0]);
290 }
291 
decompose(const TransformationMatrix::Matrix4 & mat,TransformationMatrix::DecomposedType & result)292 static bool decompose(const TransformationMatrix::Matrix4& mat, TransformationMatrix::DecomposedType& result)
293 {
294     TransformationMatrix::Matrix4 localMatrix;
295     memcpy(localMatrix, mat, sizeof(TransformationMatrix::Matrix4));
296 
297     // Normalize the matrix.
298     if (localMatrix[3][3] == 0)
299         return false;
300 
301     int i, j;
302     for (i = 0; i < 4; i++)
303         for (j = 0; j < 4; j++)
304             localMatrix[i][j] /= localMatrix[3][3];
305 
306     // perspectiveMatrix is used to solve for perspective, but it also provides
307     // an easy way to test for singularity of the upper 3x3 component.
308     TransformationMatrix::Matrix4 perspectiveMatrix;
309     memcpy(perspectiveMatrix, localMatrix, sizeof(TransformationMatrix::Matrix4));
310     for (i = 0; i < 3; i++)
311         perspectiveMatrix[i][3] = 0;
312     perspectiveMatrix[3][3] = 1;
313 
314     if (determinant4x4(perspectiveMatrix) == 0)
315         return false;
316 
317     // First, isolate perspective.  This is the messiest.
318     if (localMatrix[0][3] != 0 || localMatrix[1][3] != 0 || localMatrix[2][3] != 0) {
319         // rightHandSide is the right hand side of the equation.
320         Vector4 rightHandSide;
321         rightHandSide[0] = localMatrix[0][3];
322         rightHandSide[1] = localMatrix[1][3];
323         rightHandSide[2] = localMatrix[2][3];
324         rightHandSide[3] = localMatrix[3][3];
325 
326         // Solve the equation by inverting perspectiveMatrix and multiplying
327         // rightHandSide by the inverse.  (This is the easiest way, not
328         // necessarily the best.)
329         TransformationMatrix::Matrix4 inversePerspectiveMatrix, transposedInversePerspectiveMatrix;
330         inverse(perspectiveMatrix, inversePerspectiveMatrix);
331         transposeMatrix4(inversePerspectiveMatrix, transposedInversePerspectiveMatrix);
332 
333         Vector4 perspectivePoint;
334         v4MulPointByMatrix(rightHandSide, transposedInversePerspectiveMatrix, perspectivePoint);
335 
336         result.perspectiveX = perspectivePoint[0];
337         result.perspectiveY = perspectivePoint[1];
338         result.perspectiveZ = perspectivePoint[2];
339         result.perspectiveW = perspectivePoint[3];
340 
341         // Clear the perspective partition
342         localMatrix[0][3] = localMatrix[1][3] = localMatrix[2][3] = 0;
343         localMatrix[3][3] = 1;
344     } else {
345         // No perspective.
346         result.perspectiveX = result.perspectiveY = result.perspectiveZ = 0;
347         result.perspectiveW = 1;
348     }
349 
350     // Next take care of translation (easy).
351     result.translateX = localMatrix[3][0];
352     localMatrix[3][0] = 0;
353     result.translateY = localMatrix[3][1];
354     localMatrix[3][1] = 0;
355     result.translateZ = localMatrix[3][2];
356     localMatrix[3][2] = 0;
357 
358     // Vector4 type and functions need to be added to the common set.
359     Vector3 row[3], pdum3;
360 
361     // Now get scale and shear.
362     for (i = 0; i < 3; i++) {
363         row[i][0] = localMatrix[i][0];
364         row[i][1] = localMatrix[i][1];
365         row[i][2] = localMatrix[i][2];
366     }
367 
368     // Compute X scale factor and normalize first row.
369     result.scaleX = v3Length(row[0]);
370     v3Scale(row[0], 1.0);
371 
372     // Compute XY shear factor and make 2nd row orthogonal to 1st.
373     result.skewXY = v3Dot(row[0], row[1]);
374     v3Combine(row[1], row[0], row[1], 1.0, -result.skewXY);
375 
376     // Now, compute Y scale and normalize 2nd row.
377     result.scaleY = v3Length(row[1]);
378     v3Scale(row[1], 1.0);
379     result.skewXY /= result.scaleY;
380 
381     // Compute XZ and YZ shears, orthogonalize 3rd row.
382     result.skewXZ = v3Dot(row[0], row[2]);
383     v3Combine(row[2], row[0], row[2], 1.0, -result.skewXZ);
384     result.skewYZ = v3Dot(row[1], row[2]);
385     v3Combine(row[2], row[1], row[2], 1.0, -result.skewYZ);
386 
387     // Next, get Z scale and normalize 3rd row.
388     result.scaleZ = v3Length(row[2]);
389     v3Scale(row[2], 1.0);
390     result.skewXZ /= result.scaleZ;
391     result.skewYZ /= result.scaleZ;
392 
393     // At this point, the matrix (in rows[]) is orthonormal.
394     // Check for a coordinate system flip.  If the determinant
395     // is -1, then negate the matrix and the scaling factors.
396     v3Cross(row[1], row[2], pdum3);
397     if (v3Dot(row[0], pdum3) < 0) {
398         for (i = 0; i < 3; i++) {
399             result.scaleX *= -1;
400             row[i][0] *= -1;
401             row[i][1] *= -1;
402             row[i][2] *= -1;
403         }
404     }
405 
406     // Now, get the rotations out, as described in the gem.
407 
408     // FIXME - Add the ability to return either quaternions (which are
409     // easier to recompose with) or Euler angles (rx, ry, rz), which
410     // are easier for authors to deal with. The latter will only be useful
411     // when we fix https://bugs.webkit.org/show_bug.cgi?id=23799, so I
412     // will leave the Euler angle code here for now.
413 
414     // ret.rotateY = asin(-row[0][2]);
415     // if (cos(ret.rotateY) != 0) {
416     //     ret.rotateX = atan2(row[1][2], row[2][2]);
417     //     ret.rotateZ = atan2(row[0][1], row[0][0]);
418     // } else {
419     //     ret.rotateX = atan2(-row[2][0], row[1][1]);
420     //     ret.rotateZ = 0;
421     // }
422 
423     double s, t, x, y, z, w;
424 
425     t = row[0][0] + row[1][1] + row[2][2] + 1.0;
426 
427     if (t > 1e-4) {
428         s = 0.5 / sqrt(t);
429         w = 0.25 / s;
430         x = (row[2][1] - row[1][2]) * s;
431         y = (row[0][2] - row[2][0]) * s;
432         z = (row[1][0] - row[0][1]) * s;
433     } else if (row[0][0] > row[1][1] && row[0][0] > row[2][2]) {
434         s = sqrt (1.0 + row[0][0] - row[1][1] - row[2][2]) * 2.0; // S=4*qx
435         x = 0.25 * s;
436         y = (row[0][1] + row[1][0]) / s;
437         z = (row[0][2] + row[2][0]) / s;
438         w = (row[2][1] - row[1][2]) / s;
439     } else if (row[1][1] > row[2][2]) {
440         s = sqrt (1.0 + row[1][1] - row[0][0] - row[2][2]) * 2.0; // S=4*qy
441         x = (row[0][1] + row[1][0]) / s;
442         y = 0.25 * s;
443         z = (row[1][2] + row[2][1]) / s;
444         w = (row[0][2] - row[2][0]) / s;
445     } else {
446         s = sqrt(1.0 + row[2][2] - row[0][0] - row[1][1]) * 2.0; // S=4*qz
447         x = (row[0][2] + row[2][0]) / s;
448         y = (row[1][2] + row[2][1]) / s;
449         z = 0.25 * s;
450         w = (row[1][0] - row[0][1]) / s;
451     }
452 
453     result.quaternionX = x;
454     result.quaternionY = y;
455     result.quaternionZ = z;
456     result.quaternionW = w;
457 
458     return true;
459 }
460 
461 // Perform a spherical linear interpolation between the two
462 // passed quaternions with 0 <= t <= 1
slerp(double qa[4],const double qb[4],double t)463 static void slerp(double qa[4], const double qb[4], double t)
464 {
465     double ax, ay, az, aw;
466     double bx, by, bz, bw;
467     double cx, cy, cz, cw;
468     double angle;
469     double th, invth, scale, invscale;
470 
471     ax = qa[0]; ay = qa[1]; az = qa[2]; aw = qa[3];
472     bx = qb[0]; by = qb[1]; bz = qb[2]; bw = qb[3];
473 
474     angle = ax * bx + ay * by + az * bz + aw * bw;
475 
476     if (angle < 0.0) {
477         ax = -ax; ay = -ay;
478         az = -az; aw = -aw;
479         angle = -angle;
480     }
481 
482     if (angle + 1.0 > .05) {
483         if (1.0 - angle >= .05) {
484             th = acos (angle);
485             invth = 1.0 / sin (th);
486             scale = sin (th * (1.0 - t)) * invth;
487             invscale = sin (th * t) * invth;
488         } else {
489             scale = 1.0 - t;
490             invscale = t;
491         }
492     } else {
493         bx = -ay;
494         by = ax;
495         bz = -aw;
496         bw = az;
497         scale = sin(piDouble * (.5 - t));
498         invscale = sin (piDouble * t);
499     }
500 
501     cx = ax * scale + bx * invscale;
502     cy = ay * scale + by * invscale;
503     cz = az * scale + bz * invscale;
504     cw = aw * scale + bw * invscale;
505 
506     qa[0] = cx; qa[1] = cy; qa[2] = cz; qa[3] = cw;
507 }
508 
509 // End of Supporting Math Functions
510 
scale(double s)511 TransformationMatrix& TransformationMatrix::scale(double s)
512 {
513     return scaleNonUniform(s, s);
514 }
515 
rotateFromVector(double x,double y)516 TransformationMatrix& TransformationMatrix::rotateFromVector(double x, double y)
517 {
518     return rotate(rad2deg(atan2(y, x)));
519 }
520 
flipX()521 TransformationMatrix& TransformationMatrix::flipX()
522 {
523     return scaleNonUniform(-1.0f, 1.0f);
524 }
525 
flipY()526 TransformationMatrix& TransformationMatrix::flipY()
527 {
528     return scaleNonUniform(1.0f, -1.0f);
529 }
530 
projectPoint(const FloatPoint & p) const531 FloatPoint TransformationMatrix::projectPoint(const FloatPoint& p) const
532 {
533     // This is basically raytracing. We have a point in the destination
534     // plane with z=0, and we cast a ray parallel to the z-axis from that
535     // point to find the z-position at which it intersects the z=0 plane
536     // with the transform applied. Once we have that point we apply the
537     // inverse transform to find the corresponding point in the source
538     // space.
539     //
540     // Given a plane with normal Pn, and a ray starting at point R0 and
541     // with direction defined by the vector Rd, we can find the
542     // intersection point as a distance d from R0 in units of Rd by:
543     //
544     // d = -dot (Pn', R0) / dot (Pn', Rd)
545 
546     double x = p.x();
547     double y = p.y();
548     double z = -(m13() * x + m23() * y + m43()) / m33();
549 
550     double outX = x * m11() + y * m21() + z * m31() + m41();
551     double outY = x * m12() + y * m22() + z * m32() + m42();
552 
553     double w = x * m14() + y * m24() + z * m34() + m44();
554     if (w != 1 && w != 0) {
555         outX /= w;
556         outY /= w;
557     }
558 
559     return FloatPoint(static_cast<float>(outX), static_cast<float>(outY));
560 }
561 
projectQuad(const FloatQuad & q) const562 FloatQuad TransformationMatrix::projectQuad(const FloatQuad& q) const
563 {
564     FloatQuad projectedQuad;
565     projectedQuad.setP1(projectPoint(q.p1()));
566     projectedQuad.setP2(projectPoint(q.p2()));
567     projectedQuad.setP3(projectPoint(q.p3()));
568     projectedQuad.setP4(projectPoint(q.p4()));
569     return projectedQuad;
570 }
571 
mapPoint(const FloatPoint & p) const572 FloatPoint TransformationMatrix::mapPoint(const FloatPoint& p) const
573 {
574     if (isIdentityOrTranslation())
575         return FloatPoint(p.x() + static_cast<float>(m_matrix[3][0]), p.y() + static_cast<float>(m_matrix[3][1]));
576 
577     double x, y;
578     multVecMatrix(p.x(), p.y(), x, y);
579     return FloatPoint(static_cast<float>(x), static_cast<float>(y));
580 }
581 
mapPoint(const FloatPoint3D & p) const582 FloatPoint3D TransformationMatrix::mapPoint(const FloatPoint3D& p) const
583 {
584     if (isIdentityOrTranslation())
585         return FloatPoint3D(p.x() + static_cast<float>(m_matrix[3][0]),
586                             p.y() + static_cast<float>(m_matrix[3][1]),
587                             p.z() + static_cast<float>(m_matrix[3][2]));
588 
589     double x, y, z;
590     multVecMatrix(p.x(), p.y(), p.z(), x, y, z);
591     return FloatPoint3D(static_cast<float>(x), static_cast<float>(y), static_cast<float>(z));
592 }
593 
mapRect(const IntRect & rect) const594 IntRect TransformationMatrix::mapRect(const IntRect &rect) const
595 {
596     return enclosingIntRect(mapRect(FloatRect(rect)));
597 }
598 
mapRect(const FloatRect & r) const599 FloatRect TransformationMatrix::mapRect(const FloatRect& r) const
600 {
601     if (isIdentityOrTranslation()) {
602         FloatRect mappedRect(r);
603         mappedRect.move(static_cast<float>(m_matrix[3][0]), static_cast<float>(m_matrix[3][1]));
604         return mappedRect;
605     }
606 
607     FloatQuad resultQuad = mapQuad(FloatQuad(r));
608     return resultQuad.boundingBox();
609 }
610 
mapQuad(const FloatQuad & q) const611 FloatQuad TransformationMatrix::mapQuad(const FloatQuad& q) const
612 {
613     if (isIdentityOrTranslation()) {
614         FloatQuad mappedQuad(q);
615         mappedQuad.move(static_cast<float>(m_matrix[3][0]), static_cast<float>(m_matrix[3][1]));
616         return mappedQuad;
617     }
618 
619     FloatQuad result;
620     result.setP1(mapPoint(q.p1()));
621     result.setP2(mapPoint(q.p2()));
622     result.setP3(mapPoint(q.p3()));
623     result.setP4(mapPoint(q.p4()));
624     return result;
625 }
626 
scaleNonUniform(double sx,double sy)627 TransformationMatrix& TransformationMatrix::scaleNonUniform(double sx, double sy)
628 {
629     TransformationMatrix mat;
630     mat.m_matrix[0][0] = sx;
631     mat.m_matrix[1][1] = sy;
632 
633     multiply(mat);
634     return *this;
635 }
636 
scale3d(double sx,double sy,double sz)637 TransformationMatrix& TransformationMatrix::scale3d(double sx, double sy, double sz)
638 {
639     TransformationMatrix mat;
640     mat.m_matrix[0][0] = sx;
641     mat.m_matrix[1][1] = sy;
642     mat.m_matrix[2][2] = sz;
643 
644     multiply(mat);
645     return *this;
646 }
647 
rotate3d(double x,double y,double z,double angle)648 TransformationMatrix& TransformationMatrix::rotate3d(double x, double y, double z, double angle)
649 {
650     // angles are in degrees. Switch to radians
651     angle = deg2rad(angle);
652 
653     angle /= 2.0f;
654     double sinA = sin(angle);
655     double cosA = cos(angle);
656     double sinA2 = sinA * sinA;
657 
658     // normalize
659     double length = sqrt(x * x + y * y + z * z);
660     if (length == 0) {
661         // bad vector, just use something reasonable
662         x = 0;
663         y = 0;
664         z = 1;
665     } else if (length != 1) {
666         x /= length;
667         y /= length;
668         z /= length;
669     }
670 
671     TransformationMatrix mat;
672 
673     // optimize case where axis is along major axis
674     if (x == 1.0f && y == 0.0f && z == 0.0f) {
675         mat.m_matrix[0][0] = 1.0f;
676         mat.m_matrix[0][1] = 0.0f;
677         mat.m_matrix[0][2] = 0.0f;
678         mat.m_matrix[1][0] = 0.0f;
679         mat.m_matrix[1][1] = 1.0f - 2.0f * sinA2;
680         mat.m_matrix[1][2] = 2.0f * sinA * cosA;
681         mat.m_matrix[2][0] = 0.0f;
682         mat.m_matrix[2][1] = -2.0f * sinA * cosA;
683         mat.m_matrix[2][2] = 1.0f - 2.0f * sinA2;
684         mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
685         mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
686         mat.m_matrix[3][3] = 1.0f;
687     } else if (x == 0.0f && y == 1.0f && z == 0.0f) {
688         mat.m_matrix[0][0] = 1.0f - 2.0f * sinA2;
689         mat.m_matrix[0][1] = 0.0f;
690         mat.m_matrix[0][2] = -2.0f * sinA * cosA;
691         mat.m_matrix[1][0] = 0.0f;
692         mat.m_matrix[1][1] = 1.0f;
693         mat.m_matrix[1][2] = 0.0f;
694         mat.m_matrix[2][0] = 2.0f * sinA * cosA;
695         mat.m_matrix[2][1] = 0.0f;
696         mat.m_matrix[2][2] = 1.0f - 2.0f * sinA2;
697         mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
698         mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
699         mat.m_matrix[3][3] = 1.0f;
700     } else if (x == 0.0f && y == 0.0f && z == 1.0f) {
701         mat.m_matrix[0][0] = 1.0f - 2.0f * sinA2;
702         mat.m_matrix[0][1] = 2.0f * sinA * cosA;
703         mat.m_matrix[0][2] = 0.0f;
704         mat.m_matrix[1][0] = -2.0f * sinA * cosA;
705         mat.m_matrix[1][1] = 1.0f - 2.0f * sinA2;
706         mat.m_matrix[1][2] = 0.0f;
707         mat.m_matrix[2][0] = 0.0f;
708         mat.m_matrix[2][1] = 0.0f;
709         mat.m_matrix[2][2] = 1.0f;
710         mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
711         mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
712         mat.m_matrix[3][3] = 1.0f;
713     } else {
714         double x2 = x*x;
715         double y2 = y*y;
716         double z2 = z*z;
717 
718         mat.m_matrix[0][0] = 1.0f - 2.0f * (y2 + z2) * sinA2;
719         mat.m_matrix[0][1] = 2.0f * (x * y * sinA2 + z * sinA * cosA);
720         mat.m_matrix[0][2] = 2.0f * (x * z * sinA2 - y * sinA * cosA);
721         mat.m_matrix[1][0] = 2.0f * (y * x * sinA2 - z * sinA * cosA);
722         mat.m_matrix[1][1] = 1.0f - 2.0f * (z2 + x2) * sinA2;
723         mat.m_matrix[1][2] = 2.0f * (y * z * sinA2 + x * sinA * cosA);
724         mat.m_matrix[2][0] = 2.0f * (z * x * sinA2 + y * sinA * cosA);
725         mat.m_matrix[2][1] = 2.0f * (z * y * sinA2 - x * sinA * cosA);
726         mat.m_matrix[2][2] = 1.0f - 2.0f * (x2 + y2) * sinA2;
727         mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
728         mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
729         mat.m_matrix[3][3] = 1.0f;
730     }
731     multiply(mat);
732     return *this;
733 }
734 
rotate3d(double rx,double ry,double rz)735 TransformationMatrix& TransformationMatrix::rotate3d(double rx, double ry, double rz)
736 {
737     // angles are in degrees. Switch to radians
738     rx = deg2rad(rx);
739     ry = deg2rad(ry);
740     rz = deg2rad(rz);
741 
742     TransformationMatrix mat;
743 
744     rz /= 2.0f;
745     double sinA = sin(rz);
746     double cosA = cos(rz);
747     double sinA2 = sinA * sinA;
748 
749     mat.m_matrix[0][0] = 1.0f - 2.0f * sinA2;
750     mat.m_matrix[0][1] = 2.0f * sinA * cosA;
751     mat.m_matrix[0][2] = 0.0f;
752     mat.m_matrix[1][0] = -2.0f * sinA * cosA;
753     mat.m_matrix[1][1] = 1.0f - 2.0f * sinA2;
754     mat.m_matrix[1][2] = 0.0f;
755     mat.m_matrix[2][0] = 0.0f;
756     mat.m_matrix[2][1] = 0.0f;
757     mat.m_matrix[2][2] = 1.0f;
758     mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
759     mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
760     mat.m_matrix[3][3] = 1.0f;
761 
762     TransformationMatrix rmat(mat);
763 
764     ry /= 2.0f;
765     sinA = sin(ry);
766     cosA = cos(ry);
767     sinA2 = sinA * sinA;
768 
769     mat.m_matrix[0][0] = 1.0f - 2.0f * sinA2;
770     mat.m_matrix[0][1] = 0.0f;
771     mat.m_matrix[0][2] = -2.0f * sinA * cosA;
772     mat.m_matrix[1][0] = 0.0f;
773     mat.m_matrix[1][1] = 1.0f;
774     mat.m_matrix[1][2] = 0.0f;
775     mat.m_matrix[2][0] = 2.0f * sinA * cosA;
776     mat.m_matrix[2][1] = 0.0f;
777     mat.m_matrix[2][2] = 1.0f - 2.0f * sinA2;
778     mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
779     mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
780     mat.m_matrix[3][3] = 1.0f;
781 
782     rmat.multiply(mat);
783 
784     rx /= 2.0f;
785     sinA = sin(rx);
786     cosA = cos(rx);
787     sinA2 = sinA * sinA;
788 
789     mat.m_matrix[0][0] = 1.0f;
790     mat.m_matrix[0][1] = 0.0f;
791     mat.m_matrix[0][2] = 0.0f;
792     mat.m_matrix[1][0] = 0.0f;
793     mat.m_matrix[1][1] = 1.0f - 2.0f * sinA2;
794     mat.m_matrix[1][2] = 2.0f * sinA * cosA;
795     mat.m_matrix[2][0] = 0.0f;
796     mat.m_matrix[2][1] = -2.0f * sinA * cosA;
797     mat.m_matrix[2][2] = 1.0f - 2.0f * sinA2;
798     mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
799     mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
800     mat.m_matrix[3][3] = 1.0f;
801 
802     rmat.multiply(mat);
803 
804     multiply(rmat);
805     return *this;
806 }
807 
translate(double tx,double ty)808 TransformationMatrix& TransformationMatrix::translate(double tx, double ty)
809 {
810     m_matrix[3][0] += tx * m_matrix[0][0] + ty * m_matrix[1][0];
811     m_matrix[3][1] += tx * m_matrix[0][1] + ty * m_matrix[1][1];
812     m_matrix[3][2] += tx * m_matrix[0][2] + ty * m_matrix[1][2];
813     m_matrix[3][3] += tx * m_matrix[0][3] + ty * m_matrix[1][3];
814     return *this;
815 }
816 
translate3d(double tx,double ty,double tz)817 TransformationMatrix& TransformationMatrix::translate3d(double tx, double ty, double tz)
818 {
819     m_matrix[3][0] += tx * m_matrix[0][0] + ty * m_matrix[1][0] + tz * m_matrix[2][0];
820     m_matrix[3][1] += tx * m_matrix[0][1] + ty * m_matrix[1][1] + tz * m_matrix[2][1];
821     m_matrix[3][2] += tx * m_matrix[0][2] + ty * m_matrix[1][2] + tz * m_matrix[2][2];
822     m_matrix[3][3] += tx * m_matrix[0][3] + ty * m_matrix[1][3] + tz * m_matrix[2][3];
823     return *this;
824 }
825 
translateRight(double tx,double ty)826 TransformationMatrix& TransformationMatrix::translateRight(double tx, double ty)
827 {
828     if (tx != 0) {
829         m_matrix[0][0] +=  m_matrix[0][3] * tx;
830         m_matrix[1][0] +=  m_matrix[1][3] * tx;
831         m_matrix[2][0] +=  m_matrix[2][3] * tx;
832         m_matrix[3][0] +=  m_matrix[3][3] * tx;
833     }
834 
835     if (ty != 0) {
836         m_matrix[0][1] +=  m_matrix[0][3] * ty;
837         m_matrix[1][1] +=  m_matrix[1][3] * ty;
838         m_matrix[2][1] +=  m_matrix[2][3] * ty;
839         m_matrix[3][1] +=  m_matrix[3][3] * ty;
840     }
841 
842     return *this;
843 }
844 
translateRight3d(double tx,double ty,double tz)845 TransformationMatrix& TransformationMatrix::translateRight3d(double tx, double ty, double tz)
846 {
847     translateRight(tx, ty);
848     if (tz != 0) {
849         m_matrix[0][2] +=  m_matrix[0][3] * tz;
850         m_matrix[1][2] +=  m_matrix[1][3] * tz;
851         m_matrix[2][2] +=  m_matrix[2][3] * tz;
852         m_matrix[3][2] +=  m_matrix[3][3] * tz;
853     }
854 
855     return *this;
856 }
857 
skew(double sx,double sy)858 TransformationMatrix& TransformationMatrix::skew(double sx, double sy)
859 {
860     // angles are in degrees. Switch to radians
861     sx = deg2rad(sx);
862     sy = deg2rad(sy);
863 
864     TransformationMatrix mat;
865     mat.m_matrix[0][1] = tan(sy); // note that the y shear goes in the first row
866     mat.m_matrix[1][0] = tan(sx); // and the x shear in the second row
867 
868     multiply(mat);
869     return *this;
870 }
871 
applyPerspective(double p)872 TransformationMatrix& TransformationMatrix::applyPerspective(double p)
873 {
874     TransformationMatrix mat;
875     if (p != 0)
876         mat.m_matrix[2][3] = -1/p;
877 
878     multiply(mat);
879     return *this;
880 }
881 
rectToRect(const FloatRect & from,const FloatRect & to)882 TransformationMatrix TransformationMatrix::rectToRect(const FloatRect& from, const FloatRect& to)
883 {
884     ASSERT(!from.isEmpty());
885     return TransformationMatrix(to.width() / from.width(),
886                                 0, 0,
887                                 to.height() / from.height(),
888                                 to.x() - from.x(),
889                                 to.y() - from.y());
890 }
891 
892 //
893 // *this = mat * *this
894 //
multiply(const TransformationMatrix & mat)895 TransformationMatrix& TransformationMatrix::multiply(const TransformationMatrix& mat)
896 {
897     Matrix4 tmp;
898 
899     tmp[0][0] = (mat.m_matrix[0][0] * m_matrix[0][0] + mat.m_matrix[0][1] * m_matrix[1][0]
900                + mat.m_matrix[0][2] * m_matrix[2][0] + mat.m_matrix[0][3] * m_matrix[3][0]);
901     tmp[0][1] = (mat.m_matrix[0][0] * m_matrix[0][1] + mat.m_matrix[0][1] * m_matrix[1][1]
902                + mat.m_matrix[0][2] * m_matrix[2][1] + mat.m_matrix[0][3] * m_matrix[3][1]);
903     tmp[0][2] = (mat.m_matrix[0][0] * m_matrix[0][2] + mat.m_matrix[0][1] * m_matrix[1][2]
904                + mat.m_matrix[0][2] * m_matrix[2][2] + mat.m_matrix[0][3] * m_matrix[3][2]);
905     tmp[0][3] = (mat.m_matrix[0][0] * m_matrix[0][3] + mat.m_matrix[0][1] * m_matrix[1][3]
906                + mat.m_matrix[0][2] * m_matrix[2][3] + mat.m_matrix[0][3] * m_matrix[3][3]);
907 
908     tmp[1][0] = (mat.m_matrix[1][0] * m_matrix[0][0] + mat.m_matrix[1][1] * m_matrix[1][0]
909                + mat.m_matrix[1][2] * m_matrix[2][0] + mat.m_matrix[1][3] * m_matrix[3][0]);
910     tmp[1][1] = (mat.m_matrix[1][0] * m_matrix[0][1] + mat.m_matrix[1][1] * m_matrix[1][1]
911                + mat.m_matrix[1][2] * m_matrix[2][1] + mat.m_matrix[1][3] * m_matrix[3][1]);
912     tmp[1][2] = (mat.m_matrix[1][0] * m_matrix[0][2] + mat.m_matrix[1][1] * m_matrix[1][2]
913                + mat.m_matrix[1][2] * m_matrix[2][2] + mat.m_matrix[1][3] * m_matrix[3][2]);
914     tmp[1][3] = (mat.m_matrix[1][0] * m_matrix[0][3] + mat.m_matrix[1][1] * m_matrix[1][3]
915                + mat.m_matrix[1][2] * m_matrix[2][3] + mat.m_matrix[1][3] * m_matrix[3][3]);
916 
917     tmp[2][0] = (mat.m_matrix[2][0] * m_matrix[0][0] + mat.m_matrix[2][1] * m_matrix[1][0]
918                + mat.m_matrix[2][2] * m_matrix[2][0] + mat.m_matrix[2][3] * m_matrix[3][0]);
919     tmp[2][1] = (mat.m_matrix[2][0] * m_matrix[0][1] + mat.m_matrix[2][1] * m_matrix[1][1]
920                + mat.m_matrix[2][2] * m_matrix[2][1] + mat.m_matrix[2][3] * m_matrix[3][1]);
921     tmp[2][2] = (mat.m_matrix[2][0] * m_matrix[0][2] + mat.m_matrix[2][1] * m_matrix[1][2]
922                + mat.m_matrix[2][2] * m_matrix[2][2] + mat.m_matrix[2][3] * m_matrix[3][2]);
923     tmp[2][3] = (mat.m_matrix[2][0] * m_matrix[0][3] + mat.m_matrix[2][1] * m_matrix[1][3]
924                + mat.m_matrix[2][2] * m_matrix[2][3] + mat.m_matrix[2][3] * m_matrix[3][3]);
925 
926     tmp[3][0] = (mat.m_matrix[3][0] * m_matrix[0][0] + mat.m_matrix[3][1] * m_matrix[1][0]
927                + mat.m_matrix[3][2] * m_matrix[2][0] + mat.m_matrix[3][3] * m_matrix[3][0]);
928     tmp[3][1] = (mat.m_matrix[3][0] * m_matrix[0][1] + mat.m_matrix[3][1] * m_matrix[1][1]
929                + mat.m_matrix[3][2] * m_matrix[2][1] + mat.m_matrix[3][3] * m_matrix[3][1]);
930     tmp[3][2] = (mat.m_matrix[3][0] * m_matrix[0][2] + mat.m_matrix[3][1] * m_matrix[1][2]
931                + mat.m_matrix[3][2] * m_matrix[2][2] + mat.m_matrix[3][3] * m_matrix[3][2]);
932     tmp[3][3] = (mat.m_matrix[3][0] * m_matrix[0][3] + mat.m_matrix[3][1] * m_matrix[1][3]
933                + mat.m_matrix[3][2] * m_matrix[2][3] + mat.m_matrix[3][3] * m_matrix[3][3]);
934 
935     setMatrix(tmp);
936     return *this;
937 }
938 
multVecMatrix(double x,double y,double & resultX,double & resultY) const939 void TransformationMatrix::multVecMatrix(double x, double y, double& resultX, double& resultY) const
940 {
941     resultX = m_matrix[3][0] + x * m_matrix[0][0] + y * m_matrix[1][0];
942     resultY = m_matrix[3][1] + x * m_matrix[0][1] + y * m_matrix[1][1];
943     double w = m_matrix[3][3] + x * m_matrix[0][3] + y * m_matrix[1][3];
944     if (w != 1 && w != 0) {
945         resultX /= w;
946         resultY /= w;
947     }
948 }
949 
multVecMatrix(double x,double y,double z,double & resultX,double & resultY,double & resultZ) const950 void TransformationMatrix::multVecMatrix(double x, double y, double z, double& resultX, double& resultY, double& resultZ) const
951 {
952     resultX = m_matrix[3][0] + x * m_matrix[0][0] + y * m_matrix[1][0] + z * m_matrix[2][0];
953     resultY = m_matrix[3][1] + x * m_matrix[0][1] + y * m_matrix[1][1] + z * m_matrix[2][1];
954     resultZ = m_matrix[3][2] + x * m_matrix[0][2] + y * m_matrix[1][2] + z * m_matrix[2][2];
955     double w = m_matrix[3][3] + x * m_matrix[0][3] + y * m_matrix[1][3] + z * m_matrix[2][3];
956     if (w != 1 && w != 0) {
957         resultX /= w;
958         resultY /= w;
959         resultZ /= w;
960     }
961 }
962 
isInvertible() const963 bool TransformationMatrix::isInvertible() const
964 {
965     if (isIdentityOrTranslation())
966         return true;
967 
968     double det = WebCore::determinant4x4(m_matrix);
969 
970     if (fabs(det) < SMALL_NUMBER)
971         return false;
972 
973     return true;
974 }
975 
inverse() const976 TransformationMatrix TransformationMatrix::inverse() const
977 {
978     if (isIdentityOrTranslation()) {
979         // identity matrix
980         if (m_matrix[3][0] == 0 && m_matrix[3][1] == 0 && m_matrix[3][2] == 0)
981             return TransformationMatrix();
982 
983         // translation
984         return TransformationMatrix(1, 0, 0, 0,
985                                     0, 1, 0, 0,
986                                     0, 0, 1, 0,
987                                     -m_matrix[3][0], -m_matrix[3][1], -m_matrix[3][2], 1);
988     }
989 
990     TransformationMatrix invMat;
991     bool inverted = WebCore::inverse(m_matrix, invMat.m_matrix);
992     if (!inverted)
993         return TransformationMatrix();
994 
995     return invMat;
996 }
997 
makeAffine()998 void TransformationMatrix::makeAffine()
999 {
1000     m_matrix[0][2] = 0;
1001     m_matrix[0][3] = 0;
1002 
1003     m_matrix[1][2] = 0;
1004     m_matrix[1][3] = 0;
1005 
1006     m_matrix[2][0] = 0;
1007     m_matrix[2][1] = 0;
1008     m_matrix[2][2] = 1;
1009     m_matrix[2][3] = 0;
1010 
1011     m_matrix[3][2] = 0;
1012     m_matrix[3][3] = 1;
1013 }
1014 
toAffineTransform() const1015 AffineTransform TransformationMatrix::toAffineTransform() const
1016 {
1017     return AffineTransform(m_matrix[0][0], m_matrix[0][1], m_matrix[1][0],
1018                            m_matrix[1][1], m_matrix[3][0], m_matrix[3][1]);
1019 }
1020 
blendFloat(double & from,double to,double progress)1021 static inline void blendFloat(double& from, double to, double progress)
1022 {
1023     if (from != to)
1024         from = from + (to - from) * progress;
1025 }
1026 
blend(const TransformationMatrix & from,double progress)1027 void TransformationMatrix::blend(const TransformationMatrix& from, double progress)
1028 {
1029     if (from.isIdentity() && isIdentity())
1030         return;
1031 
1032     // decompose
1033     DecomposedType fromDecomp;
1034     DecomposedType toDecomp;
1035     from.decompose(fromDecomp);
1036     decompose(toDecomp);
1037 
1038     // interpolate
1039     blendFloat(fromDecomp.scaleX, toDecomp.scaleX, progress);
1040     blendFloat(fromDecomp.scaleY, toDecomp.scaleY, progress);
1041     blendFloat(fromDecomp.scaleZ, toDecomp.scaleZ, progress);
1042     blendFloat(fromDecomp.skewXY, toDecomp.skewXY, progress);
1043     blendFloat(fromDecomp.skewXZ, toDecomp.skewXZ, progress);
1044     blendFloat(fromDecomp.skewYZ, toDecomp.skewYZ, progress);
1045     blendFloat(fromDecomp.translateX, toDecomp.translateX, progress);
1046     blendFloat(fromDecomp.translateY, toDecomp.translateY, progress);
1047     blendFloat(fromDecomp.translateZ, toDecomp.translateZ, progress);
1048     blendFloat(fromDecomp.perspectiveX, toDecomp.perspectiveX, progress);
1049     blendFloat(fromDecomp.perspectiveY, toDecomp.perspectiveY, progress);
1050     blendFloat(fromDecomp.perspectiveZ, toDecomp.perspectiveZ, progress);
1051     blendFloat(fromDecomp.perspectiveW, toDecomp.perspectiveW, progress);
1052 
1053     slerp(&fromDecomp.quaternionX, &toDecomp.quaternionX, progress);
1054 
1055     // recompose
1056     recompose(fromDecomp);
1057 }
1058 
decompose(DecomposedType & decomp) const1059 bool TransformationMatrix::decompose(DecomposedType& decomp) const
1060 {
1061     if (isIdentity()) {
1062         memset(&decomp, 0, sizeof(decomp));
1063         decomp.perspectiveW = 1;
1064         decomp.scaleX = 1;
1065         decomp.scaleY = 1;
1066         decomp.scaleZ = 1;
1067     }
1068 
1069     if (!WebCore::decompose(m_matrix, decomp))
1070         return false;
1071     return true;
1072 }
1073 
recompose(const DecomposedType & decomp)1074 void TransformationMatrix::recompose(const DecomposedType& decomp)
1075 {
1076     makeIdentity();
1077 
1078     // first apply perspective
1079     m_matrix[0][3] = (float) decomp.perspectiveX;
1080     m_matrix[1][3] = (float) decomp.perspectiveY;
1081     m_matrix[2][3] = (float) decomp.perspectiveZ;
1082     m_matrix[3][3] = (float) decomp.perspectiveW;
1083 
1084     // now translate
1085     translate3d((float) decomp.translateX, (float) decomp.translateY, (float) decomp.translateZ);
1086 
1087     // apply rotation
1088     double xx = decomp.quaternionX * decomp.quaternionX;
1089     double xy = decomp.quaternionX * decomp.quaternionY;
1090     double xz = decomp.quaternionX * decomp.quaternionZ;
1091     double xw = decomp.quaternionX * decomp.quaternionW;
1092     double yy = decomp.quaternionY * decomp.quaternionY;
1093     double yz = decomp.quaternionY * decomp.quaternionZ;
1094     double yw = decomp.quaternionY * decomp.quaternionW;
1095     double zz = decomp.quaternionZ * decomp.quaternionZ;
1096     double zw = decomp.quaternionZ * decomp.quaternionW;
1097 
1098     // Construct a composite rotation matrix from the quaternion values
1099     TransformationMatrix rotationMatrix(1 - 2 * (yy + zz), 2 * (xy - zw), 2 * (xz + yw), 0,
1100                            2 * (xy + zw), 1 - 2 * (xx + zz), 2 * (yz - xw), 0,
1101                            2 * (xz - yw), 2 * (yz + xw), 1 - 2 * (xx + yy), 0,
1102                            0, 0, 0, 1);
1103 
1104     multiply(rotationMatrix);
1105 
1106     // now apply skew
1107     if (decomp.skewYZ) {
1108         TransformationMatrix tmp;
1109         tmp.setM32((float) decomp.skewYZ);
1110         multiply(tmp);
1111     }
1112 
1113     if (decomp.skewXZ) {
1114         TransformationMatrix tmp;
1115         tmp.setM31((float) decomp.skewXZ);
1116         multiply(tmp);
1117     }
1118 
1119     if (decomp.skewXY) {
1120         TransformationMatrix tmp;
1121         tmp.setM21((float) decomp.skewXY);
1122         multiply(tmp);
1123     }
1124 
1125     // finally, apply scale
1126     scale3d((float) decomp.scaleX, (float) decomp.scaleY, (float) decomp.scaleZ);
1127 }
1128 
1129 }
1130