Split LayoutTypes.h to improve compile time
[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) 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 all points on the quad had w < 0, then the entire quad would not be visible to the projected surface.
607     bool everythingWasClipped = clamped1 && clamped2 && clamped3 && clamped4;
608     if (everythingWasClipped)
609         return FloatQuad();
610
611     return projectedQuad;
612 }
613
614 static float clampEdgeValue(float f)
615 {
616     ASSERT(!isnan(f));
617     return min<float>(max<float>(f, -FractionalLayoutUnit::max() / 2), FractionalLayoutUnit::max() / 2);
618 }
619
620 FractionalLayoutRect TransformationMatrix::clampedBoundsOfProjectedQuad(const FloatQuad& q) const
621 {
622     FloatRect mappedQuadBounds = projectQuad(q).boundingBox();
623
624     float left = clampEdgeValue(floorf(mappedQuadBounds.x()));
625     float top = clampEdgeValue(floorf(mappedQuadBounds.y()));
626
627     float right;
628     if (isinf(mappedQuadBounds.x()) && isinf(mappedQuadBounds.width()))
629         right = FractionalLayoutUnit::max() / 2;
630     else
631         right = clampEdgeValue(ceilf(mappedQuadBounds.maxX()));
632
633     float bottom;
634     if (isinf(mappedQuadBounds.y()) && isinf(mappedQuadBounds.height()))
635         bottom = FractionalLayoutUnit::max() / 2;
636     else
637         bottom = clampEdgeValue(ceilf(mappedQuadBounds.maxY()));
638
639     return FractionalLayoutRect(FractionalLayoutUnit::clamp(left), FractionalLayoutUnit::clamp(top), 
640                                 FractionalLayoutUnit::clamp(right - left), FractionalLayoutUnit::clamp(bottom - top));
641 }
642
643 FloatPoint TransformationMatrix::mapPoint(const FloatPoint& p) const
644 {
645     if (isIdentityOrTranslation())
646         return FloatPoint(p.x() + static_cast<float>(m_matrix[3][0]), p.y() + static_cast<float>(m_matrix[3][1]));
647
648     double x, y;
649     multVecMatrix(p.x(), p.y(), x, y);
650     return FloatPoint(static_cast<float>(x), static_cast<float>(y));
651 }
652
653 FloatPoint3D TransformationMatrix::mapPoint(const FloatPoint3D& p) const
654 {
655     if (isIdentityOrTranslation())
656         return FloatPoint3D(p.x() + static_cast<float>(m_matrix[3][0]),
657                             p.y() + static_cast<float>(m_matrix[3][1]),
658                             p.z() + static_cast<float>(m_matrix[3][2]));
659
660     double x, y, z;
661     multVecMatrix(p.x(), p.y(), p.z(), x, y, z);
662     return FloatPoint3D(static_cast<float>(x), static_cast<float>(y), static_cast<float>(z));
663 }
664
665 IntRect TransformationMatrix::mapRect(const IntRect &rect) const
666 {
667     return enclosingIntRect(mapRect(FloatRect(rect)));
668 }
669
670 FractionalLayoutRect TransformationMatrix::mapRect(const FractionalLayoutRect& r) const
671 {
672     return enclosingFractionalLayoutRect(mapRect(FloatRect(r)));
673 }
674
675 FloatRect TransformationMatrix::mapRect(const FloatRect& r) const
676 {
677     if (isIdentityOrTranslation()) {
678         FloatRect mappedRect(r);
679         mappedRect.move(static_cast<float>(m_matrix[3][0]), static_cast<float>(m_matrix[3][1]));
680         return mappedRect;
681     }
682
683     FloatQuad resultQuad = mapQuad(FloatQuad(r));
684     return resultQuad.boundingBox();
685 }
686
687 FloatQuad TransformationMatrix::mapQuad(const FloatQuad& q) const
688 {
689     if (isIdentityOrTranslation()) {
690         FloatQuad mappedQuad(q);
691         mappedQuad.move(static_cast<float>(m_matrix[3][0]), static_cast<float>(m_matrix[3][1]));
692         return mappedQuad;
693     }
694
695     FloatQuad result;
696     result.setP1(mapPoint(q.p1()));
697     result.setP2(mapPoint(q.p2()));
698     result.setP3(mapPoint(q.p3()));
699     result.setP4(mapPoint(q.p4()));
700     return result;
701 }
702
703 TransformationMatrix& TransformationMatrix::scaleNonUniform(double sx, double sy)
704 {
705     m_matrix[0][0] *= sx;
706     m_matrix[0][1] *= sx;
707     m_matrix[0][2] *= sx;
708     m_matrix[0][3] *= sx;
709     
710     m_matrix[1][0] *= sy;
711     m_matrix[1][1] *= sy;
712     m_matrix[1][2] *= sy;
713     m_matrix[1][3] *= sy;
714     return *this;
715 }
716
717 TransformationMatrix& TransformationMatrix::scale3d(double sx, double sy, double sz)
718 {
719     scaleNonUniform(sx, sy);
720     
721     m_matrix[2][0] *= sz;
722     m_matrix[2][1] *= sz;
723     m_matrix[2][2] *= sz;
724     m_matrix[2][3] *= sz;
725     return *this;
726 }
727
728 TransformationMatrix& TransformationMatrix::rotate3d(double x, double y, double z, double angle)
729 {
730     // Angles are in degrees. Switch to radians.
731     angle = deg2rad(angle);
732
733     double sinTheta = sin(angle);
734     double cosTheta = cos(angle);
735     
736     // Normalize the axis of rotation
737     double length = sqrt(x * x + y * y + z * z);
738     if (length == 0) {
739         // bad vector, just use something reasonable
740         x = 0;
741         y = 0;
742         z = 1;
743     } else if (length != 1) {
744         x /= length;
745         y /= length;
746         z /= length;
747     }
748     
749     TransformationMatrix mat;
750
751     // Optimize cases where the axis is along a major axis
752     if (x == 1.0f && y == 0.0f && z == 0.0f) {
753         mat.m_matrix[0][0] = 1.0f;
754         mat.m_matrix[0][1] = 0.0f;
755         mat.m_matrix[0][2] = 0.0f;
756         mat.m_matrix[1][0] = 0.0f;
757         mat.m_matrix[1][1] = cosTheta;
758         mat.m_matrix[1][2] = sinTheta;
759         mat.m_matrix[2][0] = 0.0f;
760         mat.m_matrix[2][1] = -sinTheta;
761         mat.m_matrix[2][2] = cosTheta;
762         mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
763         mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
764         mat.m_matrix[3][3] = 1.0f;
765     } else if (x == 0.0f && y == 1.0f && z == 0.0f) {
766         mat.m_matrix[0][0] = cosTheta;
767         mat.m_matrix[0][1] = 0.0f;
768         mat.m_matrix[0][2] = -sinTheta;
769         mat.m_matrix[1][0] = 0.0f;
770         mat.m_matrix[1][1] = 1.0f;
771         mat.m_matrix[1][2] = 0.0f;
772         mat.m_matrix[2][0] = sinTheta;
773         mat.m_matrix[2][1] = 0.0f;
774         mat.m_matrix[2][2] = cosTheta;
775         mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
776         mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
777         mat.m_matrix[3][3] = 1.0f;
778     } else if (x == 0.0f && y == 0.0f && z == 1.0f) {
779         mat.m_matrix[0][0] = cosTheta;
780         mat.m_matrix[0][1] = sinTheta;
781         mat.m_matrix[0][2] = 0.0f;
782         mat.m_matrix[1][0] = -sinTheta;
783         mat.m_matrix[1][1] = cosTheta;
784         mat.m_matrix[1][2] = 0.0f;
785         mat.m_matrix[2][0] = 0.0f;
786         mat.m_matrix[2][1] = 0.0f;
787         mat.m_matrix[2][2] = 1.0f;
788         mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
789         mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
790         mat.m_matrix[3][3] = 1.0f;
791     } else {
792         // This case is the rotation about an arbitrary unit vector.
793         //
794         // Formula is adapted from Wikipedia article on Rotation matrix,
795         // http://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
796         //
797         // An alternate resource with the same matrix: http://www.fastgraph.com/makegames/3drotation/
798         //
799         double oneMinusCosTheta = 1 - cosTheta;
800         mat.m_matrix[0][0] = cosTheta + x * x * oneMinusCosTheta;
801         mat.m_matrix[0][1] = y * x * oneMinusCosTheta + z * sinTheta;
802         mat.m_matrix[0][2] = z * x * oneMinusCosTheta - y * sinTheta;
803         mat.m_matrix[1][0] = x * y * oneMinusCosTheta - z * sinTheta;
804         mat.m_matrix[1][1] = cosTheta + y * y * oneMinusCosTheta;
805         mat.m_matrix[1][2] = z * y * oneMinusCosTheta + x * sinTheta;
806         mat.m_matrix[2][0] = x * z * oneMinusCosTheta + y * sinTheta;
807         mat.m_matrix[2][1] = y * z * oneMinusCosTheta - x * sinTheta;
808         mat.m_matrix[2][2] = cosTheta + z * z * oneMinusCosTheta;
809         mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
810         mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
811         mat.m_matrix[3][3] = 1.0f;
812     }
813     multiply(mat);
814     return *this;
815 }
816
817 TransformationMatrix& TransformationMatrix::rotate3d(double rx, double ry, double rz)
818 {
819     // Angles are in degrees. Switch to radians.
820     rx = deg2rad(rx);
821     ry = deg2rad(ry);
822     rz = deg2rad(rz);
823     
824     TransformationMatrix mat;
825     
826     double sinTheta = sin(rz);
827     double cosTheta = cos(rz);
828     
829     mat.m_matrix[0][0] = cosTheta;
830     mat.m_matrix[0][1] = sinTheta;
831     mat.m_matrix[0][2] = 0.0f;
832     mat.m_matrix[1][0] = -sinTheta;
833     mat.m_matrix[1][1] = cosTheta;
834     mat.m_matrix[1][2] = 0.0f;
835     mat.m_matrix[2][0] = 0.0f;
836     mat.m_matrix[2][1] = 0.0f;
837     mat.m_matrix[2][2] = 1.0f;
838     mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
839     mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
840     mat.m_matrix[3][3] = 1.0f;
841     
842     TransformationMatrix rmat(mat);
843     
844     sinTheta = sin(ry);
845     cosTheta = cos(ry);
846     
847     mat.m_matrix[0][0] = cosTheta;
848     mat.m_matrix[0][1] = 0.0f;
849     mat.m_matrix[0][2] = -sinTheta;
850     mat.m_matrix[1][0] = 0.0f;
851     mat.m_matrix[1][1] = 1.0f;
852     mat.m_matrix[1][2] = 0.0f;
853     mat.m_matrix[2][0] = sinTheta;
854     mat.m_matrix[2][1] = 0.0f;
855     mat.m_matrix[2][2] = cosTheta;
856     mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
857     mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
858     mat.m_matrix[3][3] = 1.0f;
859
860     rmat.multiply(mat);
861
862     sinTheta = sin(rx);
863     cosTheta = cos(rx);
864     
865     mat.m_matrix[0][0] = 1.0f;
866     mat.m_matrix[0][1] = 0.0f;
867     mat.m_matrix[0][2] = 0.0f;
868     mat.m_matrix[1][0] = 0.0f;
869     mat.m_matrix[1][1] = cosTheta;
870     mat.m_matrix[1][2] = sinTheta;
871     mat.m_matrix[2][0] = 0.0f;
872     mat.m_matrix[2][1] = -sinTheta;
873     mat.m_matrix[2][2] = cosTheta;
874     mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
875     mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
876     mat.m_matrix[3][3] = 1.0f;
877
878     rmat.multiply(mat);
879
880     multiply(rmat);
881     return *this;
882 }
883
884 TransformationMatrix& TransformationMatrix::translate(double tx, double ty)
885 {
886     m_matrix[3][0] += tx * m_matrix[0][0] + ty * m_matrix[1][0];
887     m_matrix[3][1] += tx * m_matrix[0][1] + ty * m_matrix[1][1];
888     m_matrix[3][2] += tx * m_matrix[0][2] + ty * m_matrix[1][2];
889     m_matrix[3][3] += tx * m_matrix[0][3] + ty * m_matrix[1][3];
890     return *this;
891 }
892
893 TransformationMatrix& TransformationMatrix::translate3d(double tx, double ty, double tz)
894 {
895     m_matrix[3][0] += tx * m_matrix[0][0] + ty * m_matrix[1][0] + tz * m_matrix[2][0];
896     m_matrix[3][1] += tx * m_matrix[0][1] + ty * m_matrix[1][1] + tz * m_matrix[2][1];
897     m_matrix[3][2] += tx * m_matrix[0][2] + ty * m_matrix[1][2] + tz * m_matrix[2][2];
898     m_matrix[3][3] += tx * m_matrix[0][3] + ty * m_matrix[1][3] + tz * m_matrix[2][3];
899     return *this;
900 }
901
902 TransformationMatrix& TransformationMatrix::translateRight(double tx, double ty)
903 {
904     if (tx != 0) {
905         m_matrix[0][0] +=  m_matrix[0][3] * tx;
906         m_matrix[1][0] +=  m_matrix[1][3] * tx;
907         m_matrix[2][0] +=  m_matrix[2][3] * tx;
908         m_matrix[3][0] +=  m_matrix[3][3] * tx;
909     }
910
911     if (ty != 0) {
912         m_matrix[0][1] +=  m_matrix[0][3] * ty;
913         m_matrix[1][1] +=  m_matrix[1][3] * ty;
914         m_matrix[2][1] +=  m_matrix[2][3] * ty;
915         m_matrix[3][1] +=  m_matrix[3][3] * ty;
916     }
917
918     return *this;
919 }
920
921 TransformationMatrix& TransformationMatrix::translateRight3d(double tx, double ty, double tz)
922 {
923     translateRight(tx, ty);
924     if (tz != 0) {
925         m_matrix[0][2] +=  m_matrix[0][3] * tz;
926         m_matrix[1][2] +=  m_matrix[1][3] * tz;
927         m_matrix[2][2] +=  m_matrix[2][3] * tz;
928         m_matrix[3][2] +=  m_matrix[3][3] * tz;
929     }
930
931     return *this;
932 }
933
934 TransformationMatrix& TransformationMatrix::skew(double sx, double sy)
935 {
936     // angles are in degrees. Switch to radians
937     sx = deg2rad(sx);
938     sy = deg2rad(sy);
939     
940     TransformationMatrix mat;
941     mat.m_matrix[0][1] = tan(sy); // note that the y shear goes in the first row
942     mat.m_matrix[1][0] = tan(sx); // and the x shear in the second row
943
944     multiply(mat);
945     return *this;
946 }
947
948 TransformationMatrix& TransformationMatrix::applyPerspective(double p)
949 {
950     TransformationMatrix mat;
951     if (p != 0)
952         mat.m_matrix[2][3] = -1/p;
953
954     multiply(mat);
955     return *this;
956 }
957
958 TransformationMatrix TransformationMatrix::rectToRect(const FloatRect& from, const FloatRect& to)
959 {
960     ASSERT(!from.isEmpty());
961     return TransformationMatrix(to.width() / from.width(),
962                                 0, 0,
963                                 to.height() / from.height(),
964                                 to.x() - from.x(),
965                                 to.y() - from.y());
966 }
967
968 //
969 // *this = mat * *this
970 //
971 TransformationMatrix& TransformationMatrix::multiply(const TransformationMatrix& mat)
972 {
973     Matrix4 tmp;
974     
975     tmp[0][0] = (mat.m_matrix[0][0] * m_matrix[0][0] + mat.m_matrix[0][1] * m_matrix[1][0]
976                + mat.m_matrix[0][2] * m_matrix[2][0] + mat.m_matrix[0][3] * m_matrix[3][0]);
977     tmp[0][1] = (mat.m_matrix[0][0] * m_matrix[0][1] + mat.m_matrix[0][1] * m_matrix[1][1]
978                + mat.m_matrix[0][2] * m_matrix[2][1] + mat.m_matrix[0][3] * m_matrix[3][1]);
979     tmp[0][2] = (mat.m_matrix[0][0] * m_matrix[0][2] + mat.m_matrix[0][1] * m_matrix[1][2]
980                + mat.m_matrix[0][2] * m_matrix[2][2] + mat.m_matrix[0][3] * m_matrix[3][2]);
981     tmp[0][3] = (mat.m_matrix[0][0] * m_matrix[0][3] + mat.m_matrix[0][1] * m_matrix[1][3]
982                + mat.m_matrix[0][2] * m_matrix[2][3] + mat.m_matrix[0][3] * m_matrix[3][3]);
983
984     tmp[1][0] = (mat.m_matrix[1][0] * m_matrix[0][0] + mat.m_matrix[1][1] * m_matrix[1][0]
985                + mat.m_matrix[1][2] * m_matrix[2][0] + mat.m_matrix[1][3] * m_matrix[3][0]);
986     tmp[1][1] = (mat.m_matrix[1][0] * m_matrix[0][1] + mat.m_matrix[1][1] * m_matrix[1][1]
987                + mat.m_matrix[1][2] * m_matrix[2][1] + mat.m_matrix[1][3] * m_matrix[3][1]);
988     tmp[1][2] = (mat.m_matrix[1][0] * m_matrix[0][2] + mat.m_matrix[1][1] * m_matrix[1][2]
989                + mat.m_matrix[1][2] * m_matrix[2][2] + mat.m_matrix[1][3] * m_matrix[3][2]);
990     tmp[1][3] = (mat.m_matrix[1][0] * m_matrix[0][3] + mat.m_matrix[1][1] * m_matrix[1][3]
991                + mat.m_matrix[1][2] * m_matrix[2][3] + mat.m_matrix[1][3] * m_matrix[3][3]);
992
993     tmp[2][0] = (mat.m_matrix[2][0] * m_matrix[0][0] + mat.m_matrix[2][1] * m_matrix[1][0]
994                + mat.m_matrix[2][2] * m_matrix[2][0] + mat.m_matrix[2][3] * m_matrix[3][0]);
995     tmp[2][1] = (mat.m_matrix[2][0] * m_matrix[0][1] + mat.m_matrix[2][1] * m_matrix[1][1]
996                + mat.m_matrix[2][2] * m_matrix[2][1] + mat.m_matrix[2][3] * m_matrix[3][1]);
997     tmp[2][2] = (mat.m_matrix[2][0] * m_matrix[0][2] + mat.m_matrix[2][1] * m_matrix[1][2]
998                + mat.m_matrix[2][2] * m_matrix[2][2] + mat.m_matrix[2][3] * m_matrix[3][2]);
999     tmp[2][3] = (mat.m_matrix[2][0] * m_matrix[0][3] + mat.m_matrix[2][1] * m_matrix[1][3]
1000                + mat.m_matrix[2][2] * m_matrix[2][3] + mat.m_matrix[2][3] * m_matrix[3][3]);
1001
1002     tmp[3][0] = (mat.m_matrix[3][0] * m_matrix[0][0] + mat.m_matrix[3][1] * m_matrix[1][0]
1003                + mat.m_matrix[3][2] * m_matrix[2][0] + mat.m_matrix[3][3] * m_matrix[3][0]);
1004     tmp[3][1] = (mat.m_matrix[3][0] * m_matrix[0][1] + mat.m_matrix[3][1] * m_matrix[1][1]
1005                + mat.m_matrix[3][2] * m_matrix[2][1] + mat.m_matrix[3][3] * m_matrix[3][1]);
1006     tmp[3][2] = (mat.m_matrix[3][0] * m_matrix[0][2] + mat.m_matrix[3][1] * m_matrix[1][2]
1007                + mat.m_matrix[3][2] * m_matrix[2][2] + mat.m_matrix[3][3] * m_matrix[3][2]);
1008     tmp[3][3] = (mat.m_matrix[3][0] * m_matrix[0][3] + mat.m_matrix[3][1] * m_matrix[1][3]
1009                + mat.m_matrix[3][2] * m_matrix[2][3] + mat.m_matrix[3][3] * m_matrix[3][3]);
1010     
1011     setMatrix(tmp);
1012     return *this;
1013 }
1014
1015 void TransformationMatrix::multVecMatrix(double x, double y, double& resultX, double& resultY) const
1016 {
1017     resultX = m_matrix[3][0] + x * m_matrix[0][0] + y * m_matrix[1][0];
1018     resultY = m_matrix[3][1] + x * m_matrix[0][1] + y * m_matrix[1][1];
1019     double w = m_matrix[3][3] + x * m_matrix[0][3] + y * m_matrix[1][3];
1020     if (w != 1 && w != 0) {
1021         resultX /= w;
1022         resultY /= w;
1023     }
1024 }
1025
1026 void TransformationMatrix::multVecMatrix(double x, double y, double z, double& resultX, double& resultY, double& resultZ) const
1027 {
1028     resultX = m_matrix[3][0] + x * m_matrix[0][0] + y * m_matrix[1][0] + z * m_matrix[2][0];
1029     resultY = m_matrix[3][1] + x * m_matrix[0][1] + y * m_matrix[1][1] + z * m_matrix[2][1];
1030     resultZ = m_matrix[3][2] + x * m_matrix[0][2] + y * m_matrix[1][2] + z * m_matrix[2][2];
1031     double w = m_matrix[3][3] + x * m_matrix[0][3] + y * m_matrix[1][3] + z * m_matrix[2][3];
1032     if (w != 1 && w != 0) {
1033         resultX /= w;
1034         resultY /= w;
1035         resultZ /= w;
1036     }
1037 }
1038
1039 bool TransformationMatrix::isInvertible() const
1040 {
1041     if (isIdentityOrTranslation())
1042         return true;
1043
1044     double det = WebCore::determinant4x4(m_matrix);
1045
1046     if (fabs(det) < SMALL_NUMBER)
1047         return false;
1048
1049     return true;
1050 }
1051
1052 TransformationMatrix TransformationMatrix::inverse() const 
1053 {
1054     if (isIdentityOrTranslation()) {
1055         // identity matrix
1056         if (m_matrix[3][0] == 0 && m_matrix[3][1] == 0 && m_matrix[3][2] == 0)
1057             return TransformationMatrix();
1058         
1059         // translation
1060         return TransformationMatrix(1, 0, 0, 0,
1061                                     0, 1, 0, 0,
1062                                     0, 0, 1, 0,
1063                                     -m_matrix[3][0], -m_matrix[3][1], -m_matrix[3][2], 1);
1064     }
1065     
1066     TransformationMatrix invMat;
1067     bool inverted = WebCore::inverse(m_matrix, invMat.m_matrix);
1068     if (!inverted)
1069         return TransformationMatrix();
1070     
1071     return invMat;
1072 }
1073
1074 void TransformationMatrix::makeAffine()
1075 {
1076     m_matrix[0][2] = 0;
1077     m_matrix[0][3] = 0;
1078     
1079     m_matrix[1][2] = 0;
1080     m_matrix[1][3] = 0;
1081     
1082     m_matrix[2][0] = 0;
1083     m_matrix[2][1] = 0;
1084     m_matrix[2][2] = 1;
1085     m_matrix[2][3] = 0;
1086     
1087     m_matrix[3][2] = 0;
1088     m_matrix[3][3] = 1;
1089 }
1090
1091 AffineTransform TransformationMatrix::toAffineTransform() const
1092 {
1093     return AffineTransform(m_matrix[0][0], m_matrix[0][1], m_matrix[1][0],
1094                            m_matrix[1][1], m_matrix[3][0], m_matrix[3][1]);
1095 }
1096
1097 static inline void blendFloat(double& from, double to, double progress)
1098 {
1099     if (from != to)
1100         from = from + (to - from) * progress;
1101 }
1102
1103 void TransformationMatrix::blend(const TransformationMatrix& from, double progress)
1104 {
1105     if (from.isIdentity() && isIdentity())
1106         return;
1107         
1108     // decompose
1109     DecomposedType fromDecomp;
1110     DecomposedType toDecomp;
1111     from.decompose(fromDecomp);
1112     decompose(toDecomp);
1113
1114     // interpolate
1115     blendFloat(fromDecomp.scaleX, toDecomp.scaleX, progress);
1116     blendFloat(fromDecomp.scaleY, toDecomp.scaleY, progress);
1117     blendFloat(fromDecomp.scaleZ, toDecomp.scaleZ, progress);
1118     blendFloat(fromDecomp.skewXY, toDecomp.skewXY, progress);
1119     blendFloat(fromDecomp.skewXZ, toDecomp.skewXZ, progress);
1120     blendFloat(fromDecomp.skewYZ, toDecomp.skewYZ, progress);
1121     blendFloat(fromDecomp.translateX, toDecomp.translateX, progress);
1122     blendFloat(fromDecomp.translateY, toDecomp.translateY, progress);
1123     blendFloat(fromDecomp.translateZ, toDecomp.translateZ, progress);
1124     blendFloat(fromDecomp.perspectiveX, toDecomp.perspectiveX, progress);
1125     blendFloat(fromDecomp.perspectiveY, toDecomp.perspectiveY, progress);
1126     blendFloat(fromDecomp.perspectiveZ, toDecomp.perspectiveZ, progress);
1127     blendFloat(fromDecomp.perspectiveW, toDecomp.perspectiveW, progress);
1128     
1129     slerp(&fromDecomp.quaternionX, &toDecomp.quaternionX, progress);
1130         
1131     // recompose
1132     recompose(fromDecomp);
1133 }
1134
1135 bool TransformationMatrix::decompose(DecomposedType& decomp) const
1136 {
1137     if (isIdentity()) {
1138         memset(&decomp, 0, sizeof(decomp));
1139         decomp.perspectiveW = 1;
1140         decomp.scaleX = 1;
1141         decomp.scaleY = 1;
1142         decomp.scaleZ = 1;
1143     }
1144     
1145     if (!WebCore::decompose(m_matrix, decomp))
1146         return false;
1147     return true;
1148 }
1149
1150 void TransformationMatrix::recompose(const DecomposedType& decomp)
1151 {
1152     makeIdentity();
1153     
1154     // first apply perspective
1155     m_matrix[0][3] = (float) decomp.perspectiveX;
1156     m_matrix[1][3] = (float) decomp.perspectiveY;
1157     m_matrix[2][3] = (float) decomp.perspectiveZ;
1158     m_matrix[3][3] = (float) decomp.perspectiveW;
1159     
1160     // now translate
1161     translate3d((float) decomp.translateX, (float) decomp.translateY, (float) decomp.translateZ);
1162     
1163     // apply rotation
1164     double xx = decomp.quaternionX * decomp.quaternionX;
1165     double xy = decomp.quaternionX * decomp.quaternionY;
1166     double xz = decomp.quaternionX * decomp.quaternionZ;
1167     double xw = decomp.quaternionX * decomp.quaternionW;
1168     double yy = decomp.quaternionY * decomp.quaternionY;
1169     double yz = decomp.quaternionY * decomp.quaternionZ;
1170     double yw = decomp.quaternionY * decomp.quaternionW;
1171     double zz = decomp.quaternionZ * decomp.quaternionZ;
1172     double zw = decomp.quaternionZ * decomp.quaternionW;
1173     
1174     // Construct a composite rotation matrix from the quaternion values
1175     TransformationMatrix rotationMatrix(1 - 2 * (yy + zz), 2 * (xy - zw), 2 * (xz + yw), 0, 
1176                            2 * (xy + zw), 1 - 2 * (xx + zz), 2 * (yz - xw), 0,
1177                            2 * (xz - yw), 2 * (yz + xw), 1 - 2 * (xx + yy), 0,
1178                            0, 0, 0, 1);
1179     
1180     multiply(rotationMatrix);
1181     
1182     // now apply skew
1183     if (decomp.skewYZ) {
1184         TransformationMatrix tmp;
1185         tmp.setM32((float) decomp.skewYZ);
1186         multiply(tmp);
1187     }
1188     
1189     if (decomp.skewXZ) {
1190         TransformationMatrix tmp;
1191         tmp.setM31((float) decomp.skewXZ);
1192         multiply(tmp);
1193     }
1194     
1195     if (decomp.skewXY) {
1196         TransformationMatrix tmp;
1197         tmp.setM21((float) decomp.skewXY);
1198         multiply(tmp);
1199     }
1200     
1201     // finally, apply scale
1202     scale3d((float) decomp.scaleX, (float) decomp.scaleY, (float) decomp.scaleZ);
1203 }
1204
1205 bool TransformationMatrix::isIntegerTranslation() const
1206 {
1207     if (!isIdentityOrTranslation())
1208         return false;
1209
1210     // Check for translate Z.
1211     if (m_matrix[3][2])
1212         return false;
1213
1214     // Check for non-integer translate X/Y.
1215     if (static_cast<int>(m_matrix[3][0]) != m_matrix[3][0] || static_cast<int>(m_matrix[3][1]) != m_matrix[3][1])
1216         return false;
1217
1218     return true;
1219 }
1220
1221 TransformationMatrix TransformationMatrix::to2dTransform() const
1222 {
1223     return TransformationMatrix(m_matrix[0][0], m_matrix[0][1], 0, m_matrix[0][3],
1224                                 m_matrix[1][0], m_matrix[1][1], 0, m_matrix[1][3],
1225                                 0, 0, 1, 0,
1226                                 m_matrix[3][0], m_matrix[3][1], 0, m_matrix[3][3]);
1227 }
1228
1229 void TransformationMatrix::toColumnMajorFloatArray(FloatMatrix4& result) const
1230 {
1231     result[0] = m11();
1232     result[1] = m12();
1233     result[2] = m13();
1234     result[3] = m14();
1235     result[4] = m21();
1236     result[5] = m22();
1237     result[6] = m23();
1238     result[7] = m24();
1239     result[8] = m31();
1240     result[9] = m32();
1241     result[10] = m33();
1242     result[11] = m34();
1243     result[12] = m41();
1244     result[13] = m42();
1245     result[14] = m43();
1246     result[15] = m44();
1247 }
1248
1249 bool TransformationMatrix::isBackFaceVisible() const
1250 {
1251     // Back-face visibility is determined by transforming the normal vector (0, 0, 1) and
1252     // checking the sign of the resulting z component. However, normals cannot be
1253     // transformed by the original matrix, they require being transformed by the
1254     // inverse-transpose.
1255     //
1256     // Since we know we will be using (0, 0, 1), and we only care about the z-component of
1257     // the transformed normal, then we only need the m33() element of the
1258     // inverse-transpose. Therefore we do not need the transpose.
1259     //
1260     // Additionally, if we only need the m33() element, we do not need to compute a full
1261     // inverse. Instead, knowing the inverse of a matrix is adjoint(matrix) / determinant,
1262     // we can simply compute the m33() of the adjoint (adjugate) matrix, without computing
1263     // the full adjoint.
1264
1265     double determinant = WebCore::determinant4x4(m_matrix);
1266
1267     // If the matrix is not invertible, then we assume its backface is not visible.
1268     if (fabs(determinant) < SMALL_NUMBER)
1269         return false;
1270
1271     double cofactor33 = determinant3x3(m11(), m12(), m14(), m21(), m22(), m24(), m41(), m42(), m44());
1272     double zComponentOfTransformedNormal = cofactor33 / determinant;
1273
1274     return zComponentOfTransformedNormal < 0;
1275 }
1276
1277 }