REGRESSION (tile cache layers): bits of tiled layers are missing with certain 3D...
[WebKit-https.git] / Source / WebCore / platform / graphics / transforms / TransformationMatrix.cpp
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 "FractionalLayoutRect.h"
35 #include "IntRect.h"
36
37 #include <wtf/Assertions.h>
38 #include <wtf/MathExtras.h>
39
40 using namespace std;
41
42 namespace WebCore {
43
44 //
45 // Supporting Math Functions
46 //
47 // This is a set of function from various places (attributed inline) to do things like
48 // inversion and decomposition of a 4x4 matrix. They are used throughout the code
49 //
50
51 //
52 // Adapted from Matrix Inversion by Richard Carling, Graphics Gems <http://tog.acm.org/GraphicsGems/index.html>.
53
54 // EULA: The Graphics Gems code is copyright-protected. In other words, you cannot claim the text of the code 
55 // as your own and resell it. Using the code is permitted in any program, product, or library, non-commercial 
56 // or commercial. Giving credit is not required, though is a nice gesture. The code comes as-is, and if there 
57 // are any flaws or problems with any Gems code, nobody involved with Gems - authors, editors, publishers, or 
58 // webmasters - are to be held responsible. Basically, don't be a jerk, and remember that anything free comes 
59 // with no guarantee.
60
61 // A clarification about the storage of matrix elements
62 //
63 // This class uses a 2 dimensional array internally to store the elements of the matrix.  The first index into 
64 // the array refers to the column that the element lies in; the second index refers to the row.
65 //
66 // In other words, this is the layout of the matrix:
67 //
68 // | m_matrix[0][0] m_matrix[1][0] m_matrix[2][0] m_matrix[3][0] |
69 // | m_matrix[0][1] m_matrix[1][1] m_matrix[2][1] m_matrix[3][1] |
70 // | m_matrix[0][2] m_matrix[1][2] m_matrix[2][2] m_matrix[3][2] |
71 // | m_matrix[0][3] m_matrix[1][3] m_matrix[2][3] m_matrix[3][3] |
72
73 typedef double Vector4[4];
74 typedef double Vector3[3];
75
76 const double SMALL_NUMBER = 1.e-8;
77
78 // inverse(original_matrix, inverse_matrix)
79 // 
80 // calculate the inverse of a 4x4 matrix
81 // 
82 // -1     
83 // A  = ___1__ adjoint A
84 //       det A
85
86 //  double = determinant2x2(double a, double b, double c, double d)
87 //  
88 //  calculate the determinant of a 2x2 matrix.
89
90 static double determinant2x2(double a, double b, double c, double d)
91 {
92     return a * d - b * c;
93 }
94
95 //  double = determinant3x3(a1, a2, a3, b1, b2, b3, c1, c2, c3)
96 //  
97 //  Calculate the determinant of a 3x3 matrix
98 //  in the form
99 // 
100 //      | a1,  b1,  c1 |
101 //      | a2,  b2,  c2 |
102 //      | a3,  b3,  c3 |
103
104 static double determinant3x3(double a1, double a2, double a3, double b1, double b2, double b3, double c1, double c2, double c3)
105 {
106     return a1 * determinant2x2(b2, b3, c2, c3)
107          - b1 * determinant2x2(a2, a3, c2, c3)
108          + c1 * determinant2x2(a2, a3, b2, b3);
109 }
110
111 //  double = determinant4x4(matrix)
112 //  
113 //  calculate the determinant of a 4x4 matrix.
114
115 static double determinant4x4(const TransformationMatrix::Matrix4& m)
116 {
117     // Assign to individual variable names to aid selecting
118     // correct elements
119
120     double a1 = m[0][0];
121     double b1 = m[0][1]; 
122     double c1 = m[0][2];
123     double d1 = m[0][3];
124
125     double a2 = m[1][0];
126     double b2 = m[1][1]; 
127     double c2 = m[1][2];
128     double d2 = m[1][3];
129
130     double a3 = m[2][0]; 
131     double b3 = m[2][1];
132     double c3 = m[2][2];
133     double d3 = m[2][3];
134
135     double a4 = m[3][0];
136     double b4 = m[3][1]; 
137     double c4 = m[3][2];
138     double d4 = m[3][3];
139
140     return a1 * determinant3x3(b2, b3, b4, c2, c3, c4, d2, d3, d4)
141          - b1 * determinant3x3(a2, a3, a4, c2, c3, c4, d2, d3, d4)
142          + c1 * determinant3x3(a2, a3, a4, b2, b3, b4, d2, d3, d4)
143          - d1 * determinant3x3(a2, a3, a4, b2, b3, b4, c2, c3, c4);
144 }
145
146 // adjoint( original_matrix, inverse_matrix )
147 //
148 //   calculate the adjoint of a 4x4 matrix
149 //
150 //    Let  a   denote the minor determinant of matrix A obtained by
151 //         ij
152 //
153 //    deleting the ith row and jth column from A.
154 // 
155 //                  i+j
156 //   Let  b   = (-1)    a
157 //        ij            ji
158 // 
159 //  The matrix B = (b  ) is the adjoint of A
160 //                   ij
161
162 static void adjoint(const TransformationMatrix::Matrix4& matrix, TransformationMatrix::Matrix4& result)
163 {
164     // Assign to individual variable names to aid
165     // selecting correct values
166     double a1 = matrix[0][0];
167     double b1 = matrix[0][1]; 
168     double c1 = matrix[0][2];
169     double d1 = matrix[0][3];
170
171     double a2 = matrix[1][0];
172     double b2 = matrix[1][1]; 
173     double c2 = matrix[1][2];
174     double d2 = matrix[1][3];
175
176     double a3 = matrix[2][0];
177     double b3 = matrix[2][1];
178     double c3 = matrix[2][2];
179     double d3 = matrix[2][3];
180
181     double a4 = matrix[3][0];
182     double b4 = matrix[3][1]; 
183     double c4 = matrix[3][2];
184     double d4 = matrix[3][3];
185
186     // Row column labeling reversed since we transpose rows & columns
187     result[0][0]  =   determinant3x3(b2, b3, b4, c2, c3, c4, d2, d3, d4);
188     result[1][0]  = - determinant3x3(a2, a3, a4, c2, c3, c4, d2, d3, d4);
189     result[2][0]  =   determinant3x3(a2, a3, a4, b2, b3, b4, d2, d3, d4);
190     result[3][0]  = - determinant3x3(a2, a3, a4, b2, b3, b4, c2, c3, c4);
191         
192     result[0][1]  = - determinant3x3(b1, b3, b4, c1, c3, c4, d1, d3, d4);
193     result[1][1]  =   determinant3x3(a1, a3, a4, c1, c3, c4, d1, d3, d4);
194     result[2][1]  = - determinant3x3(a1, a3, a4, b1, b3, b4, d1, d3, d4);
195     result[3][1]  =   determinant3x3(a1, a3, a4, b1, b3, b4, c1, c3, c4);
196         
197     result[0][2]  =   determinant3x3(b1, b2, b4, c1, c2, c4, d1, d2, d4);
198     result[1][2]  = - determinant3x3(a1, a2, a4, c1, c2, c4, d1, d2, d4);
199     result[2][2]  =   determinant3x3(a1, a2, a4, b1, b2, b4, d1, d2, d4);
200     result[3][2]  = - determinant3x3(a1, a2, a4, b1, b2, b4, c1, c2, c4);
201         
202     result[0][3]  = - determinant3x3(b1, b2, b3, c1, c2, c3, d1, d2, d3);
203     result[1][3]  =   determinant3x3(a1, a2, a3, c1, c2, c3, d1, d2, d3);
204     result[2][3]  = - determinant3x3(a1, a2, a3, b1, b2, b3, d1, d2, d3);
205     result[3][3]  =   determinant3x3(a1, a2, a3, b1, b2, b3, c1, c2, c3);
206 }
207
208 // Returns false if the matrix is not invertible
209 static bool inverse(const TransformationMatrix::Matrix4& matrix, TransformationMatrix::Matrix4& result)
210 {
211     // Calculate the adjoint matrix
212     adjoint(matrix, result);
213
214     // Calculate the 4x4 determinant
215     // If the determinant is zero, 
216     // then the inverse matrix is not unique.
217     double det = determinant4x4(matrix);
218
219     if (fabs(det) < SMALL_NUMBER)
220         return false;
221
222     // Scale the adjoint matrix to get the inverse
223
224     for (int i = 0; i < 4; i++)
225         for (int j = 0; j < 4; j++)
226             result[i][j] = result[i][j] / det;
227
228     return true;
229 }
230
231 // End of code adapted from Matrix Inversion by Richard Carling
232
233 // Perform a decomposition on the passed matrix, return false if unsuccessful
234 // From Graphics Gems: unmatrix.c
235
236 // Transpose rotation portion of matrix a, return b
237 static void transposeMatrix4(const TransformationMatrix::Matrix4& a, TransformationMatrix::Matrix4& b)
238 {
239     for (int i = 0; i < 4; i++)
240         for (int j = 0; j < 4; j++)
241             b[i][j] = a[j][i];
242 }
243
244 // Multiply a homogeneous point by a matrix and return the transformed point
245 static void v4MulPointByMatrix(const Vector4 p, const TransformationMatrix::Matrix4& m, Vector4 result)
246 {
247     result[0] = (p[0] * m[0][0]) + (p[1] * m[1][0]) +
248                 (p[2] * m[2][0]) + (p[3] * m[3][0]);
249     result[1] = (p[0] * m[0][1]) + (p[1] * m[1][1]) +
250                 (p[2] * m[2][1]) + (p[3] * m[3][1]);
251     result[2] = (p[0] * m[0][2]) + (p[1] * m[1][2]) +
252                 (p[2] * m[2][2]) + (p[3] * m[3][2]);
253     result[3] = (p[0] * m[0][3]) + (p[1] * m[1][3]) +
254                 (p[2] * m[2][3]) + (p[3] * m[3][3]);
255 }
256
257 static double v3Length(Vector3 a)
258 {
259     return sqrt((a[0] * a[0]) + (a[1] * a[1]) + (a[2] * a[2]));
260 }
261
262 static void v3Scale(Vector3 v, double desiredLength) 
263 {
264     double len = v3Length(v);
265     if (len != 0) {
266         double l = desiredLength / len;
267         v[0] *= l;
268         v[1] *= l;
269         v[2] *= l;
270     }
271 }
272
273 static double v3Dot(const Vector3 a, const Vector3 b) 
274 {
275     return (a[0] * b[0]) + (a[1] * b[1]) + (a[2] * b[2]);
276 }
277
278 // Make a linear combination of two vectors and return the result.
279 // result = (a * ascl) + (b * bscl)
280 static void v3Combine(const Vector3 a, const Vector3 b, Vector3 result, double ascl, double bscl)
281 {
282     result[0] = (ascl * a[0]) + (bscl * b[0]);
283     result[1] = (ascl * a[1]) + (bscl * b[1]);
284     result[2] = (ascl * a[2]) + (bscl * b[2]);
285 }
286
287 // Return the cross product result = a cross b */
288 static void v3Cross(const Vector3 a, const Vector3 b, Vector3 result)
289 {
290     result[0] = (a[1] * b[2]) - (a[2] * b[1]);
291     result[1] = (a[2] * b[0]) - (a[0] * b[2]);
292     result[2] = (a[0] * b[1]) - (a[1] * b[0]);
293 }
294
295 static bool decompose(const TransformationMatrix::Matrix4& mat, TransformationMatrix::DecomposedType& result)
296 {
297     TransformationMatrix::Matrix4 localMatrix;
298     memcpy(localMatrix, mat, sizeof(TransformationMatrix::Matrix4));
299
300     // Normalize the matrix.
301     if (localMatrix[3][3] == 0)
302         return false;
303
304     int i, j;
305     for (i = 0; i < 4; i++)
306         for (j = 0; j < 4; j++)
307             localMatrix[i][j] /= localMatrix[3][3];
308
309     // perspectiveMatrix is used to solve for perspective, but it also provides
310     // an easy way to test for singularity of the upper 3x3 component.
311     TransformationMatrix::Matrix4 perspectiveMatrix;
312     memcpy(perspectiveMatrix, localMatrix, sizeof(TransformationMatrix::Matrix4));
313     for (i = 0; i < 3; i++)
314         perspectiveMatrix[i][3] = 0;
315     perspectiveMatrix[3][3] = 1;
316
317     if (determinant4x4(perspectiveMatrix) == 0)
318         return false;
319
320     // First, isolate perspective.  This is the messiest.
321     if (localMatrix[0][3] != 0 || localMatrix[1][3] != 0 || localMatrix[2][3] != 0) {
322         // rightHandSide is the right hand side of the equation.
323         Vector4 rightHandSide;
324         rightHandSide[0] = localMatrix[0][3];
325         rightHandSide[1] = localMatrix[1][3];
326         rightHandSide[2] = localMatrix[2][3];
327         rightHandSide[3] = localMatrix[3][3];
328
329         // Solve the equation by inverting perspectiveMatrix and multiplying
330         // rightHandSide by the inverse.  (This is the easiest way, not
331         // necessarily the best.)
332         TransformationMatrix::Matrix4 inversePerspectiveMatrix, transposedInversePerspectiveMatrix;
333         inverse(perspectiveMatrix, inversePerspectiveMatrix);
334         transposeMatrix4(inversePerspectiveMatrix, transposedInversePerspectiveMatrix);
335
336         Vector4 perspectivePoint;
337         v4MulPointByMatrix(rightHandSide, transposedInversePerspectiveMatrix, perspectivePoint);
338  
339         result.perspectiveX = perspectivePoint[0];
340         result.perspectiveY = perspectivePoint[1];
341         result.perspectiveZ = perspectivePoint[2];
342         result.perspectiveW = perspectivePoint[3];
343         
344         // Clear the perspective partition
345         localMatrix[0][3] = localMatrix[1][3] = localMatrix[2][3] = 0;
346         localMatrix[3][3] = 1;
347     } else {
348         // No perspective.
349         result.perspectiveX = result.perspectiveY = result.perspectiveZ = 0;
350         result.perspectiveW = 1;
351     }
352     
353     // Next take care of translation (easy).
354     result.translateX = localMatrix[3][0];
355     localMatrix[3][0] = 0;
356     result.translateY = localMatrix[3][1];
357     localMatrix[3][1] = 0;
358     result.translateZ = localMatrix[3][2];
359     localMatrix[3][2] = 0;
360
361     // Vector4 type and functions need to be added to the common set.
362     Vector3 row[3], pdum3;
363
364     // Now get scale and shear.
365     for (i = 0; i < 3; i++) {
366         row[i][0] = localMatrix[i][0];
367         row[i][1] = localMatrix[i][1];
368         row[i][2] = localMatrix[i][2];
369     }
370
371     // Compute X scale factor and normalize first row.
372     result.scaleX = v3Length(row[0]);
373     v3Scale(row[0], 1.0);
374
375     // Compute XY shear factor and make 2nd row orthogonal to 1st.
376     result.skewXY = v3Dot(row[0], row[1]);
377     v3Combine(row[1], row[0], row[1], 1.0, -result.skewXY);
378
379     // Now, compute Y scale and normalize 2nd row.
380     result.scaleY = v3Length(row[1]);
381     v3Scale(row[1], 1.0);
382     result.skewXY /= result.scaleY;
383
384     // Compute XZ and YZ shears, orthogonalize 3rd row.
385     result.skewXZ = v3Dot(row[0], row[2]);
386     v3Combine(row[2], row[0], row[2], 1.0, -result.skewXZ);
387     result.skewYZ = v3Dot(row[1], row[2]);
388     v3Combine(row[2], row[1], row[2], 1.0, -result.skewYZ);
389
390     // Next, get Z scale and normalize 3rd row.
391     result.scaleZ = v3Length(row[2]);
392     v3Scale(row[2], 1.0);
393     result.skewXZ /= result.scaleZ;
394     result.skewYZ /= result.scaleZ;
395  
396     // At this point, the matrix (in rows[]) is orthonormal.
397     // Check for a coordinate system flip.  If the determinant
398     // is -1, then negate the matrix and the scaling factors.
399     v3Cross(row[1], row[2], pdum3);
400     if (v3Dot(row[0], pdum3) < 0) {
401
402         result.scaleX *= -1;
403         result.scaleY *= -1;
404         result.scaleZ *= -1;
405
406         for (i = 0; i < 3; i++) {
407             row[i][0] *= -1;
408             row[i][1] *= -1;
409             row[i][2] *= -1;
410         }
411     }
412  
413     // Now, get the rotations out, as described in the gem.
414     
415     // FIXME - Add the ability to return either quaternions (which are
416     // easier to recompose with) or Euler angles (rx, ry, rz), which
417     // are easier for authors to deal with. The latter will only be useful
418     // when we fix https://bugs.webkit.org/show_bug.cgi?id=23799, so I
419     // will leave the Euler angle code here for now.
420
421     // ret.rotateY = asin(-row[0][2]);
422     // if (cos(ret.rotateY) != 0) {
423     //     ret.rotateX = atan2(row[1][2], row[2][2]);
424     //     ret.rotateZ = atan2(row[0][1], row[0][0]);
425     // } else {
426     //     ret.rotateX = atan2(-row[2][0], row[1][1]);
427     //     ret.rotateZ = 0;
428     // }
429     
430     double s, t, x, y, z, w;
431
432     t = row[0][0] + row[1][1] + row[2][2] + 1.0;
433
434     if (t > 1e-4) {
435         s = 0.5 / sqrt(t);
436         w = 0.25 / s;
437         x = (row[2][1] - row[1][2]) * s;
438         y = (row[0][2] - row[2][0]) * s;
439         z = (row[1][0] - row[0][1]) * s;
440     } else if (row[0][0] > row[1][1] && row[0][0] > row[2][2]) { 
441         s = sqrt (1.0 + row[0][0] - row[1][1] - row[2][2]) * 2.0; // S=4*qx 
442         x = 0.25 * s;
443         y = (row[0][1] + row[1][0]) / s; 
444         z = (row[0][2] + row[2][0]) / s; 
445         w = (row[2][1] - row[1][2]) / s;
446     } else if (row[1][1] > row[2][2]) { 
447         s = sqrt (1.0 + row[1][1] - row[0][0] - row[2][2]) * 2.0; // S=4*qy
448         x = (row[0][1] + row[1][0]) / s; 
449         y = 0.25 * s;
450         z = (row[1][2] + row[2][1]) / s; 
451         w = (row[0][2] - row[2][0]) / s;
452     } else { 
453         s = sqrt(1.0 + row[2][2] - row[0][0] - row[1][1]) * 2.0; // S=4*qz
454         x = (row[0][2] + row[2][0]) / s;
455         y = (row[1][2] + row[2][1]) / s; 
456         z = 0.25 * s;
457         w = (row[1][0] - row[0][1]) / s;
458     }
459
460     result.quaternionX = x;
461     result.quaternionY = y;
462     result.quaternionZ = z;
463     result.quaternionW = w;
464     
465     return true;
466 }
467
468 // Perform a spherical linear interpolation between the two
469 // passed quaternions with 0 <= t <= 1
470 static void slerp(double qa[4], const double qb[4], double t)
471 {
472     double ax, ay, az, aw;
473     double bx, by, bz, bw;
474     double cx, cy, cz, cw;
475     double angle;
476     double th, invth, scale, invscale;
477
478     ax = qa[0]; ay = qa[1]; az = qa[2]; aw = qa[3];
479     bx = qb[0]; by = qb[1]; bz = qb[2]; bw = qb[3];
480
481     angle = ax * bx + ay * by + az * bz + aw * bw;
482
483     if (angle < 0.0) {
484         ax = -ax; ay = -ay;
485         az = -az; aw = -aw;
486         angle = -angle;
487     }
488
489     if (angle + 1.0 > .05) {
490         if (1.0 - angle >= .05) {
491             th = acos (angle);
492             invth = 1.0 / sin (th);
493             scale = sin (th * (1.0 - t)) * invth;
494             invscale = sin (th * t) * invth;
495         } else {
496             scale = 1.0 - t;
497             invscale = t;
498         }
499     } else {
500         bx = -ay;
501         by = ax;
502         bz = -aw;
503         bw = az;
504         scale = sin(piDouble * (.5 - t));
505         invscale = sin (piDouble * t);
506     }
507
508     cx = ax * scale + bx * invscale;
509     cy = ay * scale + by * invscale;
510     cz = az * scale + bz * invscale;
511     cw = aw * scale + bw * invscale;
512
513     qa[0] = cx; qa[1] = cy; qa[2] = cz; qa[3] = cw;
514 }
515
516 // End of Supporting Math Functions
517
518 TransformationMatrix::TransformationMatrix(const AffineTransform& t)
519 {
520     setMatrix(t.a(), t.b(), t.c(), t.d(), t.e(), t.f());
521 }
522
523 TransformationMatrix& TransformationMatrix::scale(double s)
524 {
525     return scaleNonUniform(s, s);
526 }
527
528 TransformationMatrix& TransformationMatrix::rotateFromVector(double x, double y)
529 {
530     return rotate(rad2deg(atan2(y, x)));
531 }
532
533 TransformationMatrix& TransformationMatrix::flipX()
534 {
535     return scaleNonUniform(-1.0f, 1.0f);
536 }
537
538 TransformationMatrix& TransformationMatrix::flipY()
539 {
540     return scaleNonUniform(1.0f, -1.0f);
541 }
542
543 FloatPoint TransformationMatrix::projectPoint(const FloatPoint& p, bool* clamped) const
544 {
545     // This is basically raytracing. We have a point in the destination
546     // plane with z=0, and we cast a ray parallel to the z-axis from that
547     // point to find the z-position at which it intersects the z=0 plane
548     // with the transform applied. Once we have that point we apply the
549     // inverse transform to find the corresponding point in the source
550     // space.
551     // 
552     // Given a plane with normal Pn, and a ray starting at point R0 and
553     // with direction defined by the vector Rd, we can find the
554     // intersection point as a distance d from R0 in units of Rd by:
555     // 
556     // d = -dot (Pn', R0) / dot (Pn', Rd)
557     if (clamped)
558         *clamped = false;
559
560     if (m33() == 0) {
561         // In this case, the projection plane is parallel to the ray we are trying to
562         // trace, and there is no well-defined value for the projection.
563         return FloatPoint();
564     }
565     
566     double x = p.x();
567     double y = p.y();
568     double z = -(m13() * x + m23() * y + m43()) / m33();
569
570     // FIXME: use multVecMatrix()
571     double outX = x * m11() + y * m21() + z * m31() + m41();
572     double outY = x * m12() + y * m22() + z * m32() + m42();
573
574     double w = x * m14() + y * m24() + z * m34() + m44();
575     if (w <= 0) {
576         // Using int max causes overflow when other code uses the projected point. To
577         // represent infinity yet reduce the risk of overflow, we use a large but
578         // not-too-large number here when clamping.
579         const int largeNumber = 100000000 / kFixedPointDenominator;
580         outX = copysign(largeNumber, outX);
581         outY = copysign(largeNumber, outY);
582         if (clamped)
583             *clamped = true;
584     } else if (w != 1) {
585         outX /= w;
586         outY /= w;
587     }
588
589     return FloatPoint(static_cast<float>(outX), static_cast<float>(outY));
590 }
591
592 FloatQuad TransformationMatrix::projectQuad(const FloatQuad& q, bool* clamped) const
593 {
594     FloatQuad projectedQuad;
595
596     bool clamped1 = false;
597     bool clamped2 = false;
598     bool clamped3 = false;
599     bool clamped4 = false;
600
601     projectedQuad.setP1(projectPoint(q.p1(), &clamped1));
602     projectedQuad.setP2(projectPoint(q.p2(), &clamped2));
603     projectedQuad.setP3(projectPoint(q.p3(), &clamped3));
604     projectedQuad.setP4(projectPoint(q.p4(), &clamped4));
605
606     if (clamped)
607         *clamped = clamped1 || clamped2 || clamped3 || clamped4;
608         
609     // If all points on the quad had w < 0, then the entire quad would not be visible to the projected surface.
610     bool everythingWasClipped = clamped1 && clamped2 && clamped3 && clamped4;
611     if (everythingWasClipped)
612         return FloatQuad();
613
614     return projectedQuad;
615 }
616
617 static float clampEdgeValue(float f)
618 {
619     ASSERT(!isnan(f));
620     return min<float>(max<float>(f, -FractionalLayoutUnit::max() / 2), FractionalLayoutUnit::max() / 2);
621 }
622
623 FractionalLayoutRect TransformationMatrix::clampedBoundsOfProjectedQuad(const FloatQuad& q) const
624 {
625     FloatRect mappedQuadBounds = projectQuad(q).boundingBox();
626
627     float left = clampEdgeValue(floorf(mappedQuadBounds.x()));
628     float top = clampEdgeValue(floorf(mappedQuadBounds.y()));
629
630     float right;
631     if (isinf(mappedQuadBounds.x()) && isinf(mappedQuadBounds.width()))
632         right = FractionalLayoutUnit::max() / 2;
633     else
634         right = clampEdgeValue(ceilf(mappedQuadBounds.maxX()));
635
636     float bottom;
637     if (isinf(mappedQuadBounds.y()) && isinf(mappedQuadBounds.height()))
638         bottom = FractionalLayoutUnit::max() / 2;
639     else
640         bottom = clampEdgeValue(ceilf(mappedQuadBounds.maxY()));
641
642     return FractionalLayoutRect(FractionalLayoutUnit::clamp(left), FractionalLayoutUnit::clamp(top), 
643                                 FractionalLayoutUnit::clamp(right - left), FractionalLayoutUnit::clamp(bottom - top));
644 }
645
646 FloatPoint TransformationMatrix::mapPoint(const FloatPoint& p) const
647 {
648     if (isIdentityOrTranslation())
649         return FloatPoint(p.x() + static_cast<float>(m_matrix[3][0]), p.y() + static_cast<float>(m_matrix[3][1]));
650
651     double x, y;
652     multVecMatrix(p.x(), p.y(), x, y);
653     return FloatPoint(static_cast<float>(x), static_cast<float>(y));
654 }
655
656 FloatPoint3D TransformationMatrix::mapPoint(const FloatPoint3D& p) const
657 {
658     if (isIdentityOrTranslation())
659         return FloatPoint3D(p.x() + static_cast<float>(m_matrix[3][0]),
660                             p.y() + static_cast<float>(m_matrix[3][1]),
661                             p.z() + static_cast<float>(m_matrix[3][2]));
662
663     double x, y, z;
664     multVecMatrix(p.x(), p.y(), p.z(), x, y, z);
665     return FloatPoint3D(static_cast<float>(x), static_cast<float>(y), static_cast<float>(z));
666 }
667
668 IntRect TransformationMatrix::mapRect(const IntRect &rect) const
669 {
670     return enclosingIntRect(mapRect(FloatRect(rect)));
671 }
672
673 FractionalLayoutRect TransformationMatrix::mapRect(const FractionalLayoutRect& r) const
674 {
675     return enclosingFractionalLayoutRect(mapRect(FloatRect(r)));
676 }
677
678 FloatRect TransformationMatrix::mapRect(const FloatRect& r) const
679 {
680     if (isIdentityOrTranslation()) {
681         FloatRect mappedRect(r);
682         mappedRect.move(static_cast<float>(m_matrix[3][0]), static_cast<float>(m_matrix[3][1]));
683         return mappedRect;
684     }
685
686     FloatQuad resultQuad = mapQuad(FloatQuad(r));
687     return resultQuad.boundingBox();
688 }
689
690 FloatQuad TransformationMatrix::mapQuad(const FloatQuad& q) const
691 {
692     if (isIdentityOrTranslation()) {
693         FloatQuad mappedQuad(q);
694         mappedQuad.move(static_cast<float>(m_matrix[3][0]), static_cast<float>(m_matrix[3][1]));
695         return mappedQuad;
696     }
697
698     FloatQuad result;
699     result.setP1(mapPoint(q.p1()));
700     result.setP2(mapPoint(q.p2()));
701     result.setP3(mapPoint(q.p3()));
702     result.setP4(mapPoint(q.p4()));
703     return result;
704 }
705
706 TransformationMatrix& TransformationMatrix::scaleNonUniform(double sx, double sy)
707 {
708     m_matrix[0][0] *= sx;
709     m_matrix[0][1] *= sx;
710     m_matrix[0][2] *= sx;
711     m_matrix[0][3] *= sx;
712     
713     m_matrix[1][0] *= sy;
714     m_matrix[1][1] *= sy;
715     m_matrix[1][2] *= sy;
716     m_matrix[1][3] *= sy;
717     return *this;
718 }
719
720 TransformationMatrix& TransformationMatrix::scale3d(double sx, double sy, double sz)
721 {
722     scaleNonUniform(sx, sy);
723     
724     m_matrix[2][0] *= sz;
725     m_matrix[2][1] *= sz;
726     m_matrix[2][2] *= sz;
727     m_matrix[2][3] *= sz;
728     return *this;
729 }
730
731 TransformationMatrix& TransformationMatrix::rotate3d(double x, double y, double z, double angle)
732 {
733     // Normalize the axis of rotation
734     double length = sqrt(x * x + y * y + z * z);
735     if (length == 0) {
736         // A direction vector that cannot be normalized, such as [0, 0, 0], will cause the rotation to not be applied. 
737         return *this;
738     } else if (length != 1) {
739         x /= length;
740         y /= length;
741         z /= length;
742     }
743
744     // Angles are in degrees. Switch to radians.
745     angle = deg2rad(angle);
746
747     double sinTheta = sin(angle);
748     double cosTheta = cos(angle);
749     
750     TransformationMatrix mat;
751
752     // Optimize cases where the axis is along a major axis
753     if (x == 1.0f && y == 0.0f && z == 0.0f) {
754         mat.m_matrix[0][0] = 1.0f;
755         mat.m_matrix[0][1] = 0.0f;
756         mat.m_matrix[0][2] = 0.0f;
757         mat.m_matrix[1][0] = 0.0f;
758         mat.m_matrix[1][1] = cosTheta;
759         mat.m_matrix[1][2] = sinTheta;
760         mat.m_matrix[2][0] = 0.0f;
761         mat.m_matrix[2][1] = -sinTheta;
762         mat.m_matrix[2][2] = cosTheta;
763         mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
764         mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
765         mat.m_matrix[3][3] = 1.0f;
766     } else if (x == 0.0f && y == 1.0f && z == 0.0f) {
767         mat.m_matrix[0][0] = cosTheta;
768         mat.m_matrix[0][1] = 0.0f;
769         mat.m_matrix[0][2] = -sinTheta;
770         mat.m_matrix[1][0] = 0.0f;
771         mat.m_matrix[1][1] = 1.0f;
772         mat.m_matrix[1][2] = 0.0f;
773         mat.m_matrix[2][0] = sinTheta;
774         mat.m_matrix[2][1] = 0.0f;
775         mat.m_matrix[2][2] = cosTheta;
776         mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
777         mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
778         mat.m_matrix[3][3] = 1.0f;
779     } else if (x == 0.0f && y == 0.0f && z == 1.0f) {
780         mat.m_matrix[0][0] = cosTheta;
781         mat.m_matrix[0][1] = sinTheta;
782         mat.m_matrix[0][2] = 0.0f;
783         mat.m_matrix[1][0] = -sinTheta;
784         mat.m_matrix[1][1] = cosTheta;
785         mat.m_matrix[1][2] = 0.0f;
786         mat.m_matrix[2][0] = 0.0f;
787         mat.m_matrix[2][1] = 0.0f;
788         mat.m_matrix[2][2] = 1.0f;
789         mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
790         mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
791         mat.m_matrix[3][3] = 1.0f;
792     } else {
793         // This case is the rotation about an arbitrary unit vector.
794         //
795         // Formula is adapted from Wikipedia article on Rotation matrix,
796         // http://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
797         //
798         // An alternate resource with the same matrix: http://www.fastgraph.com/makegames/3drotation/
799         //
800         double oneMinusCosTheta = 1 - cosTheta;
801         mat.m_matrix[0][0] = cosTheta + x * x * oneMinusCosTheta;
802         mat.m_matrix[0][1] = y * x * oneMinusCosTheta + z * sinTheta;
803         mat.m_matrix[0][2] = z * x * oneMinusCosTheta - y * sinTheta;
804         mat.m_matrix[1][0] = x * y * oneMinusCosTheta - z * sinTheta;
805         mat.m_matrix[1][1] = cosTheta + y * y * oneMinusCosTheta;
806         mat.m_matrix[1][2] = z * y * oneMinusCosTheta + x * sinTheta;
807         mat.m_matrix[2][0] = x * z * oneMinusCosTheta + y * sinTheta;
808         mat.m_matrix[2][1] = y * z * oneMinusCosTheta - x * sinTheta;
809         mat.m_matrix[2][2] = cosTheta + z * z * oneMinusCosTheta;
810         mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
811         mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
812         mat.m_matrix[3][3] = 1.0f;
813     }
814     multiply(mat);
815     return *this;
816 }
817
818 TransformationMatrix& TransformationMatrix::rotate3d(double rx, double ry, double rz)
819 {
820     // Angles are in degrees. Switch to radians.
821     rx = deg2rad(rx);
822     ry = deg2rad(ry);
823     rz = deg2rad(rz);
824     
825     TransformationMatrix mat;
826     
827     double sinTheta = sin(rz);
828     double cosTheta = cos(rz);
829     
830     mat.m_matrix[0][0] = cosTheta;
831     mat.m_matrix[0][1] = sinTheta;
832     mat.m_matrix[0][2] = 0.0f;
833     mat.m_matrix[1][0] = -sinTheta;
834     mat.m_matrix[1][1] = cosTheta;
835     mat.m_matrix[1][2] = 0.0f;
836     mat.m_matrix[2][0] = 0.0f;
837     mat.m_matrix[2][1] = 0.0f;
838     mat.m_matrix[2][2] = 1.0f;
839     mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
840     mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
841     mat.m_matrix[3][3] = 1.0f;
842     
843     TransformationMatrix rmat(mat);
844     
845     sinTheta = sin(ry);
846     cosTheta = cos(ry);
847     
848     mat.m_matrix[0][0] = cosTheta;
849     mat.m_matrix[0][1] = 0.0f;
850     mat.m_matrix[0][2] = -sinTheta;
851     mat.m_matrix[1][0] = 0.0f;
852     mat.m_matrix[1][1] = 1.0f;
853     mat.m_matrix[1][2] = 0.0f;
854     mat.m_matrix[2][0] = sinTheta;
855     mat.m_matrix[2][1] = 0.0f;
856     mat.m_matrix[2][2] = cosTheta;
857     mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
858     mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
859     mat.m_matrix[3][3] = 1.0f;
860
861     rmat.multiply(mat);
862
863     sinTheta = sin(rx);
864     cosTheta = cos(rx);
865     
866     mat.m_matrix[0][0] = 1.0f;
867     mat.m_matrix[0][1] = 0.0f;
868     mat.m_matrix[0][2] = 0.0f;
869     mat.m_matrix[1][0] = 0.0f;
870     mat.m_matrix[1][1] = cosTheta;
871     mat.m_matrix[1][2] = sinTheta;
872     mat.m_matrix[2][0] = 0.0f;
873     mat.m_matrix[2][1] = -sinTheta;
874     mat.m_matrix[2][2] = cosTheta;
875     mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
876     mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
877     mat.m_matrix[3][3] = 1.0f;
878
879     rmat.multiply(mat);
880
881     multiply(rmat);
882     return *this;
883 }
884
885 TransformationMatrix& TransformationMatrix::translate(double tx, double ty)
886 {
887     m_matrix[3][0] += tx * m_matrix[0][0] + ty * m_matrix[1][0];
888     m_matrix[3][1] += tx * m_matrix[0][1] + ty * m_matrix[1][1];
889     m_matrix[3][2] += tx * m_matrix[0][2] + ty * m_matrix[1][2];
890     m_matrix[3][3] += tx * m_matrix[0][3] + ty * m_matrix[1][3];
891     return *this;
892 }
893
894 TransformationMatrix& TransformationMatrix::translate3d(double tx, double ty, double tz)
895 {
896     m_matrix[3][0] += tx * m_matrix[0][0] + ty * m_matrix[1][0] + tz * m_matrix[2][0];
897     m_matrix[3][1] += tx * m_matrix[0][1] + ty * m_matrix[1][1] + tz * m_matrix[2][1];
898     m_matrix[3][2] += tx * m_matrix[0][2] + ty * m_matrix[1][2] + tz * m_matrix[2][2];
899     m_matrix[3][3] += tx * m_matrix[0][3] + ty * m_matrix[1][3] + tz * m_matrix[2][3];
900     return *this;
901 }
902
903 TransformationMatrix& TransformationMatrix::translateRight(double tx, double ty)
904 {
905     if (tx != 0) {
906         m_matrix[0][0] +=  m_matrix[0][3] * tx;
907         m_matrix[1][0] +=  m_matrix[1][3] * tx;
908         m_matrix[2][0] +=  m_matrix[2][3] * tx;
909         m_matrix[3][0] +=  m_matrix[3][3] * tx;
910     }
911
912     if (ty != 0) {
913         m_matrix[0][1] +=  m_matrix[0][3] * ty;
914         m_matrix[1][1] +=  m_matrix[1][3] * ty;
915         m_matrix[2][1] +=  m_matrix[2][3] * ty;
916         m_matrix[3][1] +=  m_matrix[3][3] * ty;
917     }
918
919     return *this;
920 }
921
922 TransformationMatrix& TransformationMatrix::translateRight3d(double tx, double ty, double tz)
923 {
924     translateRight(tx, ty);
925     if (tz != 0) {
926         m_matrix[0][2] +=  m_matrix[0][3] * tz;
927         m_matrix[1][2] +=  m_matrix[1][3] * tz;
928         m_matrix[2][2] +=  m_matrix[2][3] * tz;
929         m_matrix[3][2] +=  m_matrix[3][3] * tz;
930     }
931
932     return *this;
933 }
934
935 TransformationMatrix& TransformationMatrix::skew(double sx, double sy)
936 {
937     // angles are in degrees. Switch to radians
938     sx = deg2rad(sx);
939     sy = deg2rad(sy);
940     
941     TransformationMatrix mat;
942     mat.m_matrix[0][1] = tan(sy); // note that the y shear goes in the first row
943     mat.m_matrix[1][0] = tan(sx); // and the x shear in the second row
944
945     multiply(mat);
946     return *this;
947 }
948
949 TransformationMatrix& TransformationMatrix::applyPerspective(double p)
950 {
951     TransformationMatrix mat;
952     if (p != 0)
953         mat.m_matrix[2][3] = -1/p;
954
955     multiply(mat);
956     return *this;
957 }
958
959 TransformationMatrix TransformationMatrix::rectToRect(const FloatRect& from, const FloatRect& to)
960 {
961     ASSERT(!from.isEmpty());
962     return TransformationMatrix(to.width() / from.width(),
963                                 0, 0,
964                                 to.height() / from.height(),
965                                 to.x() - from.x(),
966                                 to.y() - from.y());
967 }
968
969 //
970 // *this = mat * *this
971 //
972 TransformationMatrix& TransformationMatrix::multiply(const TransformationMatrix& mat)
973 {
974     Matrix4 tmp;
975     
976     tmp[0][0] = (mat.m_matrix[0][0] * m_matrix[0][0] + mat.m_matrix[0][1] * m_matrix[1][0]
977                + mat.m_matrix[0][2] * m_matrix[2][0] + mat.m_matrix[0][3] * m_matrix[3][0]);
978     tmp[0][1] = (mat.m_matrix[0][0] * m_matrix[0][1] + mat.m_matrix[0][1] * m_matrix[1][1]
979                + mat.m_matrix[0][2] * m_matrix[2][1] + mat.m_matrix[0][3] * m_matrix[3][1]);
980     tmp[0][2] = (mat.m_matrix[0][0] * m_matrix[0][2] + mat.m_matrix[0][1] * m_matrix[1][2]
981                + mat.m_matrix[0][2] * m_matrix[2][2] + mat.m_matrix[0][3] * m_matrix[3][2]);
982     tmp[0][3] = (mat.m_matrix[0][0] * m_matrix[0][3] + mat.m_matrix[0][1] * m_matrix[1][3]
983                + mat.m_matrix[0][2] * m_matrix[2][3] + mat.m_matrix[0][3] * m_matrix[3][3]);
984
985     tmp[1][0] = (mat.m_matrix[1][0] * m_matrix[0][0] + mat.m_matrix[1][1] * m_matrix[1][0]
986                + mat.m_matrix[1][2] * m_matrix[2][0] + mat.m_matrix[1][3] * m_matrix[3][0]);
987     tmp[1][1] = (mat.m_matrix[1][0] * m_matrix[0][1] + mat.m_matrix[1][1] * m_matrix[1][1]
988                + mat.m_matrix[1][2] * m_matrix[2][1] + mat.m_matrix[1][3] * m_matrix[3][1]);
989     tmp[1][2] = (mat.m_matrix[1][0] * m_matrix[0][2] + mat.m_matrix[1][1] * m_matrix[1][2]
990                + mat.m_matrix[1][2] * m_matrix[2][2] + mat.m_matrix[1][3] * m_matrix[3][2]);
991     tmp[1][3] = (mat.m_matrix[1][0] * m_matrix[0][3] + mat.m_matrix[1][1] * m_matrix[1][3]
992                + mat.m_matrix[1][2] * m_matrix[2][3] + mat.m_matrix[1][3] * m_matrix[3][3]);
993
994     tmp[2][0] = (mat.m_matrix[2][0] * m_matrix[0][0] + mat.m_matrix[2][1] * m_matrix[1][0]
995                + mat.m_matrix[2][2] * m_matrix[2][0] + mat.m_matrix[2][3] * m_matrix[3][0]);
996     tmp[2][1] = (mat.m_matrix[2][0] * m_matrix[0][1] + mat.m_matrix[2][1] * m_matrix[1][1]
997                + mat.m_matrix[2][2] * m_matrix[2][1] + mat.m_matrix[2][3] * m_matrix[3][1]);
998     tmp[2][2] = (mat.m_matrix[2][0] * m_matrix[0][2] + mat.m_matrix[2][1] * m_matrix[1][2]
999                + mat.m_matrix[2][2] * m_matrix[2][2] + mat.m_matrix[2][3] * m_matrix[3][2]);
1000     tmp[2][3] = (mat.m_matrix[2][0] * m_matrix[0][3] + mat.m_matrix[2][1] * m_matrix[1][3]
1001                + mat.m_matrix[2][2] * m_matrix[2][3] + mat.m_matrix[2][3] * m_matrix[3][3]);
1002
1003     tmp[3][0] = (mat.m_matrix[3][0] * m_matrix[0][0] + mat.m_matrix[3][1] * m_matrix[1][0]
1004                + mat.m_matrix[3][2] * m_matrix[2][0] + mat.m_matrix[3][3] * m_matrix[3][0]);
1005     tmp[3][1] = (mat.m_matrix[3][0] * m_matrix[0][1] + mat.m_matrix[3][1] * m_matrix[1][1]
1006                + mat.m_matrix[3][2] * m_matrix[2][1] + mat.m_matrix[3][3] * m_matrix[3][1]);
1007     tmp[3][2] = (mat.m_matrix[3][0] * m_matrix[0][2] + mat.m_matrix[3][1] * m_matrix[1][2]
1008                + mat.m_matrix[3][2] * m_matrix[2][2] + mat.m_matrix[3][3] * m_matrix[3][2]);
1009     tmp[3][3] = (mat.m_matrix[3][0] * m_matrix[0][3] + mat.m_matrix[3][1] * m_matrix[1][3]
1010                + mat.m_matrix[3][2] * m_matrix[2][3] + mat.m_matrix[3][3] * m_matrix[3][3]);
1011     
1012     setMatrix(tmp);
1013     return *this;
1014 }
1015
1016 void TransformationMatrix::multVecMatrix(double x, double y, double& resultX, double& resultY) const
1017 {
1018     resultX = m_matrix[3][0] + x * m_matrix[0][0] + y * m_matrix[1][0];
1019     resultY = m_matrix[3][1] + x * m_matrix[0][1] + y * m_matrix[1][1];
1020     double w = m_matrix[3][3] + x * m_matrix[0][3] + y * m_matrix[1][3];
1021     if (w != 1 && w != 0) {
1022         resultX /= w;
1023         resultY /= w;
1024     }
1025 }
1026
1027 void TransformationMatrix::multVecMatrix(double x, double y, double z, double& resultX, double& resultY, double& resultZ) const
1028 {
1029     resultX = m_matrix[3][0] + x * m_matrix[0][0] + y * m_matrix[1][0] + z * m_matrix[2][0];
1030     resultY = m_matrix[3][1] + x * m_matrix[0][1] + y * m_matrix[1][1] + z * m_matrix[2][1];
1031     resultZ = m_matrix[3][2] + x * m_matrix[0][2] + y * m_matrix[1][2] + z * m_matrix[2][2];
1032     double w = m_matrix[3][3] + x * m_matrix[0][3] + y * m_matrix[1][3] + z * m_matrix[2][3];
1033     if (w != 1 && w != 0) {
1034         resultX /= w;
1035         resultY /= w;
1036         resultZ /= w;
1037     }
1038 }
1039
1040 bool TransformationMatrix::isInvertible() const
1041 {
1042     if (isIdentityOrTranslation())
1043         return true;
1044
1045     double det = WebCore::determinant4x4(m_matrix);
1046
1047     if (fabs(det) < SMALL_NUMBER)
1048         return false;
1049
1050     return true;
1051 }
1052
1053 TransformationMatrix TransformationMatrix::inverse() const 
1054 {
1055     if (isIdentityOrTranslation()) {
1056         // identity matrix
1057         if (m_matrix[3][0] == 0 && m_matrix[3][1] == 0 && m_matrix[3][2] == 0)
1058             return TransformationMatrix();
1059         
1060         // translation
1061         return TransformationMatrix(1, 0, 0, 0,
1062                                     0, 1, 0, 0,
1063                                     0, 0, 1, 0,
1064                                     -m_matrix[3][0], -m_matrix[3][1], -m_matrix[3][2], 1);
1065     }
1066     
1067     TransformationMatrix invMat;
1068     bool inverted = WebCore::inverse(m_matrix, invMat.m_matrix);
1069     if (!inverted)
1070         return TransformationMatrix();
1071     
1072     return invMat;
1073 }
1074
1075 void TransformationMatrix::makeAffine()
1076 {
1077     m_matrix[0][2] = 0;
1078     m_matrix[0][3] = 0;
1079     
1080     m_matrix[1][2] = 0;
1081     m_matrix[1][3] = 0;
1082     
1083     m_matrix[2][0] = 0;
1084     m_matrix[2][1] = 0;
1085     m_matrix[2][2] = 1;
1086     m_matrix[2][3] = 0;
1087     
1088     m_matrix[3][2] = 0;
1089     m_matrix[3][3] = 1;
1090 }
1091
1092 AffineTransform TransformationMatrix::toAffineTransform() const
1093 {
1094     return AffineTransform(m_matrix[0][0], m_matrix[0][1], m_matrix[1][0],
1095                            m_matrix[1][1], m_matrix[3][0], m_matrix[3][1]);
1096 }
1097
1098 static inline void blendFloat(double& from, double to, double progress)
1099 {
1100     if (from != to)
1101         from = from + (to - from) * progress;
1102 }
1103
1104 void TransformationMatrix::blend(const TransformationMatrix& from, double progress)
1105 {
1106     if (from.isIdentity() && isIdentity())
1107         return;
1108         
1109     // decompose
1110     DecomposedType fromDecomp;
1111     DecomposedType toDecomp;
1112     from.decompose(fromDecomp);
1113     decompose(toDecomp);
1114
1115     // interpolate
1116     blendFloat(fromDecomp.scaleX, toDecomp.scaleX, progress);
1117     blendFloat(fromDecomp.scaleY, toDecomp.scaleY, progress);
1118     blendFloat(fromDecomp.scaleZ, toDecomp.scaleZ, progress);
1119     blendFloat(fromDecomp.skewXY, toDecomp.skewXY, progress);
1120     blendFloat(fromDecomp.skewXZ, toDecomp.skewXZ, progress);
1121     blendFloat(fromDecomp.skewYZ, toDecomp.skewYZ, progress);
1122     blendFloat(fromDecomp.translateX, toDecomp.translateX, progress);
1123     blendFloat(fromDecomp.translateY, toDecomp.translateY, progress);
1124     blendFloat(fromDecomp.translateZ, toDecomp.translateZ, progress);
1125     blendFloat(fromDecomp.perspectiveX, toDecomp.perspectiveX, progress);
1126     blendFloat(fromDecomp.perspectiveY, toDecomp.perspectiveY, progress);
1127     blendFloat(fromDecomp.perspectiveZ, toDecomp.perspectiveZ, progress);
1128     blendFloat(fromDecomp.perspectiveW, toDecomp.perspectiveW, progress);
1129     
1130     slerp(&fromDecomp.quaternionX, &toDecomp.quaternionX, progress);
1131         
1132     // recompose
1133     recompose(fromDecomp);
1134 }
1135
1136 bool TransformationMatrix::decompose(DecomposedType& decomp) const
1137 {
1138     if (isIdentity()) {
1139         memset(&decomp, 0, sizeof(decomp));
1140         decomp.perspectiveW = 1;
1141         decomp.scaleX = 1;
1142         decomp.scaleY = 1;
1143         decomp.scaleZ = 1;
1144     }
1145     
1146     if (!WebCore::decompose(m_matrix, decomp))
1147         return false;
1148     return true;
1149 }
1150
1151 void TransformationMatrix::recompose(const DecomposedType& decomp)
1152 {
1153     makeIdentity();
1154     
1155     // first apply perspective
1156     m_matrix[0][3] = (float) decomp.perspectiveX;
1157     m_matrix[1][3] = (float) decomp.perspectiveY;
1158     m_matrix[2][3] = (float) decomp.perspectiveZ;
1159     m_matrix[3][3] = (float) decomp.perspectiveW;
1160     
1161     // now translate
1162     translate3d((float) decomp.translateX, (float) decomp.translateY, (float) decomp.translateZ);
1163     
1164     // apply rotation
1165     double xx = decomp.quaternionX * decomp.quaternionX;
1166     double xy = decomp.quaternionX * decomp.quaternionY;
1167     double xz = decomp.quaternionX * decomp.quaternionZ;
1168     double xw = decomp.quaternionX * decomp.quaternionW;
1169     double yy = decomp.quaternionY * decomp.quaternionY;
1170     double yz = decomp.quaternionY * decomp.quaternionZ;
1171     double yw = decomp.quaternionY * decomp.quaternionW;
1172     double zz = decomp.quaternionZ * decomp.quaternionZ;
1173     double zw = decomp.quaternionZ * decomp.quaternionW;
1174     
1175     // Construct a composite rotation matrix from the quaternion values
1176     TransformationMatrix rotationMatrix(1 - 2 * (yy + zz), 2 * (xy - zw), 2 * (xz + yw), 0, 
1177                            2 * (xy + zw), 1 - 2 * (xx + zz), 2 * (yz - xw), 0,
1178                            2 * (xz - yw), 2 * (yz + xw), 1 - 2 * (xx + yy), 0,
1179                            0, 0, 0, 1);
1180     
1181     multiply(rotationMatrix);
1182     
1183     // now apply skew
1184     if (decomp.skewYZ) {
1185         TransformationMatrix tmp;
1186         tmp.setM32((float) decomp.skewYZ);
1187         multiply(tmp);
1188     }
1189     
1190     if (decomp.skewXZ) {
1191         TransformationMatrix tmp;
1192         tmp.setM31((float) decomp.skewXZ);
1193         multiply(tmp);
1194     }
1195     
1196     if (decomp.skewXY) {
1197         TransformationMatrix tmp;
1198         tmp.setM21((float) decomp.skewXY);
1199         multiply(tmp);
1200     }
1201     
1202     // finally, apply scale
1203     scale3d((float) decomp.scaleX, (float) decomp.scaleY, (float) decomp.scaleZ);
1204 }
1205
1206 bool TransformationMatrix::isIntegerTranslation() const
1207 {
1208     if (!isIdentityOrTranslation())
1209         return false;
1210
1211     // Check for translate Z.
1212     if (m_matrix[3][2])
1213         return false;
1214
1215     // Check for non-integer translate X/Y.
1216     if (static_cast<int>(m_matrix[3][0]) != m_matrix[3][0] || static_cast<int>(m_matrix[3][1]) != m_matrix[3][1])
1217         return false;
1218
1219     return true;
1220 }
1221
1222 TransformationMatrix TransformationMatrix::to2dTransform() const
1223 {
1224     return TransformationMatrix(m_matrix[0][0], m_matrix[0][1], 0, m_matrix[0][3],
1225                                 m_matrix[1][0], m_matrix[1][1], 0, m_matrix[1][3],
1226                                 0, 0, 1, 0,
1227                                 m_matrix[3][0], m_matrix[3][1], 0, m_matrix[3][3]);
1228 }
1229
1230 void TransformationMatrix::toColumnMajorFloatArray(FloatMatrix4& result) const
1231 {
1232     result[0] = m11();
1233     result[1] = m12();
1234     result[2] = m13();
1235     result[3] = m14();
1236     result[4] = m21();
1237     result[5] = m22();
1238     result[6] = m23();
1239     result[7] = m24();
1240     result[8] = m31();
1241     result[9] = m32();
1242     result[10] = m33();
1243     result[11] = m34();
1244     result[12] = m41();
1245     result[13] = m42();
1246     result[14] = m43();
1247     result[15] = m44();
1248 }
1249
1250 bool TransformationMatrix::isBackFaceVisible() const
1251 {
1252     // Back-face visibility is determined by transforming the normal vector (0, 0, 1) and
1253     // checking the sign of the resulting z component. However, normals cannot be
1254     // transformed by the original matrix, they require being transformed by the
1255     // inverse-transpose.
1256     //
1257     // Since we know we will be using (0, 0, 1), and we only care about the z-component of
1258     // the transformed normal, then we only need the m33() element of the
1259     // inverse-transpose. Therefore we do not need the transpose.
1260     //
1261     // Additionally, if we only need the m33() element, we do not need to compute a full
1262     // inverse. Instead, knowing the inverse of a matrix is adjoint(matrix) / determinant,
1263     // we can simply compute the m33() of the adjoint (adjugate) matrix, without computing
1264     // the full adjoint.
1265
1266     double determinant = WebCore::determinant4x4(m_matrix);
1267
1268     // If the matrix is not invertible, then we assume its backface is not visible.
1269     if (fabs(determinant) < SMALL_NUMBER)
1270         return false;
1271
1272     double cofactor33 = determinant3x3(m11(), m12(), m14(), m21(), m22(), m24(), m41(), m42(), m44());
1273     double zComponentOfTransformedNormal = cofactor33 / determinant;
1274
1275     return zComponentOfTransformedNormal < 0;
1276 }
1277
1278 }