1/*
2 * Copyright (C) 2005, 2006, 2013 Apple 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 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 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 "FloatQuad.h"
32#include "FloatRect.h"
33#include "IntRect.h"
34#include "LayoutRect.h"
35#include <cmath>
36#include <wtf/Assertions.h>
37#include <wtf/MathExtras.h>
38#include <wtf/Optional.h>
39#include <wtf/text/TextStream.h>
40
41#if CPU(X86_64)
42#include <emmintrin.h>
43#endif
44
45namespace WebCore {
46
47//
48// Adapted from Matrix Inversion by Richard Carling, Graphics Gems <http://tog.acm.org/GraphicsGems/index.html>.
49
50// EULA: The Graphics Gems code is copyright-protected. In other words, you cannot claim the text of the code
51// as your own and resell it. Using the code is permitted in any program, product, or library, non-commercial
52// or commercial. Giving credit is not required, though is a nice gesture. The code comes as-is, and if there
53// are any flaws or problems with any Gems code, nobody involved with Gems - authors, editors, publishers, or
54// webmasters - are to be held responsible. Basically, don't be a jerk, and remember that anything free comes
55// with no guarantee.
56
57// A clarification about the storage of matrix elements
58//
59// This class uses a 2 dimensional array internally to store the elements of the matrix. The first index into
60// the array refers to the column that the element lies in; the second index refers to the row.
61//
62// In other words, this is the layout of the matrix:
63//
64// | m_matrix[0][0] m_matrix[1][0] m_matrix[2][0] m_matrix[3][0] |
65// | m_matrix[0][1] m_matrix[1][1] m_matrix[2][1] m_matrix[3][1] |
66// | m_matrix[0][2] m_matrix[1][2] m_matrix[2][2] m_matrix[3][2] |
67// | m_matrix[0][3] m_matrix[1][3] m_matrix[2][3] m_matrix[3][3] |
68
69typedef double Vector4[4];
70typedef double Vector3[3];
71
72const double SMALL_NUMBER = 1.e-8;
73
74const TransformationMatrix TransformationMatrix::identity { };
75
76// inverse(original_matrix, inverse_matrix)
77//
78// calculate the inverse of a 4x4 matrix
79//
80// -1
81// A = ___1__ adjoint A
82// det A
83
84// double = determinant2x2(double a, double b, double c, double d)
85//
86// calculate the determinant of a 2x2 matrix.
87
88static double determinant2x2(double a, double b, double c, double d)
89{
90 return a * d - b * c;
91}
92
93// double = determinant3x3(a1, a2, a3, b1, b2, b3, c1, c2, c3)
94//
95// Calculate the determinant of a 3x3 matrix
96// in the form
97//
98// | a1, b1, c1 |
99// | a2, b2, c2 |
100// | a3, b3, c3 |
101
102static double determinant3x3(double a1, double a2, double a3, double b1, double b2, double b3, double c1, double c2, double c3)
103{
104 return a1 * determinant2x2(b2, b3, c2, c3)
105 - b1 * determinant2x2(a2, a3, c2, c3)
106 + c1 * determinant2x2(a2, a3, b2, b3);
107}
108
109// double = determinant4x4(matrix)
110//
111// calculate the determinant of a 4x4 matrix.
112
113static double determinant4x4(const TransformationMatrix::Matrix4& m)
114{
115 // Assign to individual variable names to aid selecting
116 // correct elements
117
118 double a1 = m[0][0];
119 double b1 = m[0][1];
120 double c1 = m[0][2];
121 double d1 = m[0][3];
122
123 double a2 = m[1][0];
124 double b2 = m[1][1];
125 double c2 = m[1][2];
126 double d2 = m[1][3];
127
128 double a3 = m[2][0];
129 double b3 = m[2][1];
130 double c3 = m[2][2];
131 double d3 = m[2][3];
132
133 double a4 = m[3][0];
134 double b4 = m[3][1];
135 double c4 = m[3][2];
136 double d4 = m[3][3];
137
138 return a1 * determinant3x3(b2, b3, b4, c2, c3, c4, d2, d3, d4)
139 - b1 * determinant3x3(a2, a3, a4, c2, c3, c4, d2, d3, d4)
140 + c1 * determinant3x3(a2, a3, a4, b2, b3, b4, d2, d3, d4)
141 - d1 * determinant3x3(a2, a3, a4, b2, b3, b4, c2, c3, c4);
142}
143
144// adjoint( original_matrix, inverse_matrix )
145//
146// calculate the adjoint of a 4x4 matrix
147//
148// Let a denote the minor determinant of matrix A obtained by
149// ij
150//
151// deleting the ith row and jth column from A.
152//
153// i+j
154// Let b = (-1) a
155// ij ji
156//
157// The matrix B = (b ) is the adjoint of A
158// ij
159
160static void adjoint(const TransformationMatrix::Matrix4& matrix, TransformationMatrix::Matrix4& result)
161{
162 // Assign to individual variable names to aid
163 // selecting correct values
164 double a1 = matrix[0][0];
165 double b1 = matrix[0][1];
166 double c1 = matrix[0][2];
167 double d1 = matrix[0][3];
168
169 double a2 = matrix[1][0];
170 double b2 = matrix[1][1];
171 double c2 = matrix[1][2];
172 double d2 = matrix[1][3];
173
174 double a3 = matrix[2][0];
175 double b3 = matrix[2][1];
176 double c3 = matrix[2][2];
177 double d3 = matrix[2][3];
178
179 double a4 = matrix[3][0];
180 double b4 = matrix[3][1];
181 double c4 = matrix[3][2];
182 double d4 = matrix[3][3];
183
184 // Row column labeling reversed since we transpose rows & columns
185 result[0][0] = determinant3x3(b2, b3, b4, c2, c3, c4, d2, d3, d4);
186 result[1][0] = - determinant3x3(a2, a3, a4, c2, c3, c4, d2, d3, d4);
187 result[2][0] = determinant3x3(a2, a3, a4, b2, b3, b4, d2, d3, d4);
188 result[3][0] = - determinant3x3(a2, a3, a4, b2, b3, b4, c2, c3, c4);
189
190 result[0][1] = - determinant3x3(b1, b3, b4, c1, c3, c4, d1, d3, d4);
191 result[1][1] = determinant3x3(a1, a3, a4, c1, c3, c4, d1, d3, d4);
192 result[2][1] = - determinant3x3(a1, a3, a4, b1, b3, b4, d1, d3, d4);
193 result[3][1] = determinant3x3(a1, a3, a4, b1, b3, b4, c1, c3, c4);
194
195 result[0][2] = determinant3x3(b1, b2, b4, c1, c2, c4, d1, d2, d4);
196 result[1][2] = - determinant3x3(a1, a2, a4, c1, c2, c4, d1, d2, d4);
197 result[2][2] = determinant3x3(a1, a2, a4, b1, b2, b4, d1, d2, d4);
198 result[3][2] = - determinant3x3(a1, a2, a4, b1, b2, b4, c1, c2, c4);
199
200 result[0][3] = - determinant3x3(b1, b2, b3, c1, c2, c3, d1, d2, d3);
201 result[1][3] = determinant3x3(a1, a2, a3, c1, c2, c3, d1, d2, d3);
202 result[2][3] = - determinant3x3(a1, a2, a3, b1, b2, b3, d1, d2, d3);
203 result[3][3] = determinant3x3(a1, a2, a3, b1, b2, b3, c1, c2, c3);
204}
205
206// Returns false if the matrix is not invertible
207static bool inverse(const TransformationMatrix::Matrix4& matrix, TransformationMatrix::Matrix4& result)
208{
209 // Calculate the adjoint matrix
210 adjoint(matrix, result);
211
212 // Calculate the 4x4 determinant
213 // If the determinant is zero,
214 // then the inverse matrix is not unique.
215 double det = determinant4x4(matrix);
216
217 if (fabs(det) < SMALL_NUMBER)
218 return false;
219
220 // Scale the adjoint matrix to get the inverse
221
222 for (int i = 0; i < 4; i++)
223 for (int j = 0; j < 4; j++)
224 result[i][j] = result[i][j] / det;
225
226 return true;
227}
228
229// End of code adapted from Matrix Inversion by Richard Carling
230
231// Perform a decomposition on the passed matrix, return false if unsuccessful
232// From Graphics Gems: unmatrix.c
233
234// Transpose rotation portion of matrix a, return b
235static void transposeMatrix4(const TransformationMatrix::Matrix4& a, TransformationMatrix::Matrix4& b)
236{
237 for (int i = 0; i < 4; i++)
238 for (int j = 0; j < 4; j++)
239 b[i][j] = a[j][i];
240}
241
242// Multiply a homogeneous point by a matrix and return the transformed point
243static void v4MulPointByMatrix(const Vector4 p, const TransformationMatrix::Matrix4& m, Vector4 result)
244{
245 result[0] = (p[0] * m[0][0]) + (p[1] * m[1][0]) +
246 (p[2] * m[2][0]) + (p[3] * m[3][0]);
247 result[1] = (p[0] * m[0][1]) + (p[1] * m[1][1]) +
248 (p[2] * m[2][1]) + (p[3] * m[3][1]);
249 result[2] = (p[0] * m[0][2]) + (p[1] * m[1][2]) +
250 (p[2] * m[2][2]) + (p[3] * m[3][2]);
251 result[3] = (p[0] * m[0][3]) + (p[1] * m[1][3]) +
252 (p[2] * m[2][3]) + (p[3] * m[3][3]);
253}
254
255static double v3Length(Vector3 a)
256{
257 return sqrt((a[0] * a[0]) + (a[1] * a[1]) + (a[2] * a[2]));
258}
259
260static void v3Scale(Vector3 v, double desiredLength)
261{
262 double len = v3Length(v);
263 if (len != 0) {
264 double l = desiredLength / len;
265 v[0] *= l;
266 v[1] *= l;
267 v[2] *= l;
268 }
269}
270
271static double v3Dot(const Vector3 a, const Vector3 b)
272{
273 return (a[0] * b[0]) + (a[1] * b[1]) + (a[2] * b[2]);
274}
275
276// Make a linear combination of two vectors and return the result.
277// result = (a * ascl) + (b * bscl)
278static void v3Combine(const Vector3 a, const Vector3 b, Vector3 result, double ascl, double bscl)
279{
280 result[0] = (ascl * a[0]) + (bscl * b[0]);
281 result[1] = (ascl * a[1]) + (bscl * b[1]);
282 result[2] = (ascl * a[2]) + (bscl * b[2]);
283}
284
285// Return the cross product result = a cross b */
286static void v3Cross(const Vector3 a, const Vector3 b, Vector3 result)
287{
288 result[0] = (a[1] * b[2]) - (a[2] * b[1]);
289 result[1] = (a[2] * b[0]) - (a[0] * b[2]);
290 result[2] = (a[0] * b[1]) - (a[1] * b[0]);
291}
292
293static bool decompose2(const TransformationMatrix::Matrix4& matrix, TransformationMatrix::Decomposed2Type& result)
294{
295 double row0x = matrix[0][0];
296 double row0y = matrix[0][1];
297 double row1x = matrix[1][0];
298 double row1y = matrix[1][1];
299 result.translateX = matrix[3][0];
300 result.translateY = matrix[3][1];
301
302 // Compute scaling factors.
303 result.scaleX = sqrt(row0x * row0x + row0y * row0y);
304 result.scaleY = sqrt(row1x * row1x + row1y * row1y);
305
306 // If determinant is negative, one axis was flipped.
307 double determinant = row0x * row1y - row0y * row1x;
308 if (determinant < 0) {
309 // Flip axis with minimum unit vector dot product.
310 if (row0x < row1y)
311 result.scaleX = -result.scaleX;
312 else
313 result.scaleY = -result.scaleY;
314 }
315
316 // Renormalize matrix to remove scale.
317 if (result.scaleX) {
318 row0x *= 1 / result.scaleX;
319 row0y *= 1 / result.scaleX;
320 }
321 if (result.scaleY) {
322 row1x *= 1 / result.scaleY;
323 row1y *= 1 / result.scaleY;
324 }
325
326 // Compute rotation and renormalize matrix.
327 result.angle = atan2(row0y, row0x);
328
329 if (result.angle) {
330 // Rotate(-angle) = [cos(angle), sin(angle), -sin(angle), cos(angle)]
331 // = [row0x, -row0y, row0y, row0x]
332 // Thanks to the normalization above.
333 double sn = -row0y;
334 double cs = row0x;
335 double m11 = row0x, m12 = row0y;
336 double m21 = row1x, m22 = row1y;
337
338 row0x = cs * m11 + sn * m21;
339 row0y = cs * m12 + sn * m22;
340 row1x = -sn * m11 + cs * m21;
341 row1y = -sn * m12 + cs * m22;
342 }
343
344 result.m11 = row0x;
345 result.m12 = row0y;
346 result.m21 = row1x;
347 result.m22 = row1y;
348
349 // Convert into degrees because our rotation functions expect it.
350 result.angle = rad2deg(result.angle);
351
352 return true;
353}
354
355static bool decompose4(const TransformationMatrix::Matrix4& mat, TransformationMatrix::Decomposed4Type& result)
356{
357 TransformationMatrix::Matrix4 localMatrix;
358 memcpy(localMatrix, mat, sizeof(TransformationMatrix::Matrix4));
359
360 // Normalize the matrix.
361 if (localMatrix[3][3] == 0)
362 return false;
363
364 int i, j;
365 for (i = 0; i < 4; i++)
366 for (j = 0; j < 4; j++)
367 localMatrix[i][j] /= localMatrix[3][3];
368
369 // perspectiveMatrix is used to solve for perspective, but it also provides
370 // an easy way to test for singularity of the upper 3x3 component.
371 TransformationMatrix::Matrix4 perspectiveMatrix;
372 memcpy(perspectiveMatrix, localMatrix, sizeof(TransformationMatrix::Matrix4));
373 for (i = 0; i < 3; i++)
374 perspectiveMatrix[i][3] = 0;
375 perspectiveMatrix[3][3] = 1;
376
377 if (determinant4x4(perspectiveMatrix) == 0)
378 return false;
379
380 // First, isolate perspective. This is the messiest.
381 if (localMatrix[0][3] != 0 || localMatrix[1][3] != 0 || localMatrix[2][3] != 0) {
382 // rightHandSide is the right hand side of the equation.
383 Vector4 rightHandSide;
384 rightHandSide[0] = localMatrix[0][3];
385 rightHandSide[1] = localMatrix[1][3];
386 rightHandSide[2] = localMatrix[2][3];
387 rightHandSide[3] = localMatrix[3][3];
388
389 // Solve the equation by inverting perspectiveMatrix and multiplying
390 // rightHandSide by the inverse. (This is the easiest way, not
391 // necessarily the best.)
392 TransformationMatrix::Matrix4 inversePerspectiveMatrix, transposedInversePerspectiveMatrix;
393 inverse(perspectiveMatrix, inversePerspectiveMatrix);
394 transposeMatrix4(inversePerspectiveMatrix, transposedInversePerspectiveMatrix);
395
396 Vector4 perspectivePoint;
397 v4MulPointByMatrix(rightHandSide, transposedInversePerspectiveMatrix, perspectivePoint);
398
399 result.perspectiveX = perspectivePoint[0];
400 result.perspectiveY = perspectivePoint[1];
401 result.perspectiveZ = perspectivePoint[2];
402 result.perspectiveW = perspectivePoint[3];
403
404 // Clear the perspective partition
405 localMatrix[0][3] = localMatrix[1][3] = localMatrix[2][3] = 0;
406 localMatrix[3][3] = 1;
407 } else {
408 // No perspective.
409 result.perspectiveX = result.perspectiveY = result.perspectiveZ = 0;
410 result.perspectiveW = 1;
411 }
412
413 // Next take care of translation (easy).
414 result.translateX = localMatrix[3][0];
415 localMatrix[3][0] = 0;
416 result.translateY = localMatrix[3][1];
417 localMatrix[3][1] = 0;
418 result.translateZ = localMatrix[3][2];
419 localMatrix[3][2] = 0;
420
421 // Vector4 type and functions need to be added to the common set.
422 Vector3 row[3], pdum3;
423
424 // Now get scale and shear.
425 for (i = 0; i < 3; i++) {
426 row[i][0] = localMatrix[i][0];
427 row[i][1] = localMatrix[i][1];
428 row[i][2] = localMatrix[i][2];
429 }
430
431 // Compute X scale factor and normalize first row.
432 result.scaleX = v3Length(row[0]);
433 v3Scale(row[0], 1.0);
434
435 // Compute XY shear factor and make 2nd row orthogonal to 1st.
436 result.skewXY = v3Dot(row[0], row[1]);
437 v3Combine(row[1], row[0], row[1], 1.0, -result.skewXY);
438
439 // Now, compute Y scale and normalize 2nd row.
440 result.scaleY = v3Length(row[1]);
441 v3Scale(row[1], 1.0);
442 result.skewXY /= result.scaleY;
443
444 // Compute XZ and YZ shears, orthogonalize 3rd row.
445 result.skewXZ = v3Dot(row[0], row[2]);
446 v3Combine(row[2], row[0], row[2], 1.0, -result.skewXZ);
447 result.skewYZ = v3Dot(row[1], row[2]);
448 v3Combine(row[2], row[1], row[2], 1.0, -result.skewYZ);
449
450 // Next, get Z scale and normalize 3rd row.
451 result.scaleZ = v3Length(row[2]);
452 v3Scale(row[2], 1.0);
453 result.skewXZ /= result.scaleZ;
454 result.skewYZ /= result.scaleZ;
455
456 // At this point, the matrix (in rows[]) is orthonormal.
457 // Check for a coordinate system flip. If the determinant
458 // is -1, then negate the matrix and the scaling factors.
459 v3Cross(row[1], row[2], pdum3);
460 if (v3Dot(row[0], pdum3) < 0) {
461
462 result.scaleX *= -1;
463 result.scaleY *= -1;
464 result.scaleZ *= -1;
465
466 for (i = 0; i < 3; i++) {
467 row[i][0] *= -1;
468 row[i][1] *= -1;
469 row[i][2] *= -1;
470 }
471 }
472
473 // Now, get the rotations out, as described in the gem.
474
475 // FIXME - Add the ability to return either quaternions (which are
476 // easier to recompose with) or Euler angles (rx, ry, rz), which
477 // are easier for authors to deal with. The latter will only be useful
478 // when we fix https://bugs.webkit.org/show_bug.cgi?id=23799, so I
479 // will leave the Euler angle code here for now.
480
481 // ret.rotateY = asin(-row[0][2]);
482 // if (cos(ret.rotateY) != 0) {
483 // ret.rotateX = atan2(row[1][2], row[2][2]);
484 // ret.rotateZ = atan2(row[0][1], row[0][0]);
485 // } else {
486 // ret.rotateX = atan2(-row[2][0], row[1][1]);
487 // ret.rotateZ = 0;
488 // }
489
490 double s, t, x, y, z, w;
491
492 t = row[0][0] + row[1][1] + row[2][2] + 1.0;
493
494 if (t > 1e-4) {
495 s = 0.5 / sqrt(t);
496 w = 0.25 / s;
497 x = (row[2][1] - row[1][2]) * s;
498 y = (row[0][2] - row[2][0]) * s;
499 z = (row[1][0] - row[0][1]) * s;
500 } else if (row[0][0] > row[1][1] && row[0][0] > row[2][2]) {
501 s = sqrt(1.0 + row[0][0] - row[1][1] - row[2][2]) * 2.0; // S = 4 * qx.
502 x = 0.25 * s;
503 y = (row[0][1] + row[1][0]) / s;
504 z = (row[0][2] + row[2][0]) / s;
505 w = (row[2][1] - row[1][2]) / s;
506 } else if (row[1][1] > row[2][2]) {
507 s = sqrt(1.0 + row[1][1] - row[0][0] - row[2][2]) * 2.0; // S = 4 * qy.
508 x = (row[0][1] + row[1][0]) / s;
509 y = 0.25 * s;
510 z = (row[1][2] + row[2][1]) / s;
511 w = (row[0][2] - row[2][0]) / s;
512 } else {
513 s = sqrt(1.0 + row[2][2] - row[0][0] - row[1][1]) * 2.0; // S = 4 * qz.
514 x = (row[0][2] + row[2][0]) / s;
515 y = (row[1][2] + row[2][1]) / s;
516 z = 0.25 * s;
517 w = (row[1][0] - row[0][1]) / s;
518 }
519
520 result.quaternionX = x;
521 result.quaternionY = y;
522 result.quaternionZ = z;
523 result.quaternionW = w;
524
525 return true;
526}
527
528// Perform a spherical linear interpolation between the two
529// passed quaternions with 0 <= t <= 1.
530static void slerp(double qa[4], const double qb[4], double t)
531{
532 double ax, ay, az, aw;
533 double bx, by, bz, bw;
534 double cx, cy, cz, cw;
535 double angle;
536 double th, invth, scale, invscale;
537
538 ax = qa[0]; ay = qa[1]; az = qa[2]; aw = qa[3];
539 bx = qb[0]; by = qb[1]; bz = qb[2]; bw = qb[3];
540
541 angle = ax * bx + ay * by + az * bz + aw * bw;
542
543 if (angle < 0.0) {
544 ax = -ax; ay = -ay;
545 az = -az; aw = -aw;
546 angle = -angle;
547 }
548
549 if (angle + 1.0 > .05) {
550 if (1.0 - angle >= .05) {
551 th = acos (angle);
552 invth = 1.0 / sin (th);
553 scale = sin (th * (1.0 - t)) * invth;
554 invscale = sin (th * t) * invth;
555 } else {
556 scale = 1.0 - t;
557 invscale = t;
558 }
559 } else {
560 bx = -ay;
561 by = ax;
562 bz = -aw;
563 bw = az;
564 scale = sin(piDouble * (.5 - t));
565 invscale = sin (piDouble * t);
566 }
567
568 cx = ax * scale + bx * invscale;
569 cy = ay * scale + by * invscale;
570 cz = az * scale + bz * invscale;
571 cw = aw * scale + bw * invscale;
572
573 qa[0] = cx; qa[1] = cy; qa[2] = cz; qa[3] = cw;
574}
575
576// End of Supporting Math Functions
577
578TransformationMatrix::TransformationMatrix(const AffineTransform& t)
579{
580 setMatrix(t.a(), t.b(), t.c(), t.d(), t.e(), t.f());
581}
582
583TransformationMatrix& TransformationMatrix::scale(double s)
584{
585 return scaleNonUniform(s, s);
586}
587
588TransformationMatrix& TransformationMatrix::rotateFromVector(double x, double y)
589{
590 return rotate(rad2deg(atan2(y, x)));
591}
592
593TransformationMatrix& TransformationMatrix::flipX()
594{
595 return scaleNonUniform(-1.0, 1.0);
596}
597
598TransformationMatrix& TransformationMatrix::flipY()
599{
600 return scaleNonUniform(1.0, -1.0);
601}
602
603FloatPoint TransformationMatrix::projectPoint(const FloatPoint& p, bool* clamped) const
604{
605 // This is basically raytracing. We have a point in the destination
606 // plane with z=0, and we cast a ray parallel to the z-axis from that
607 // point to find the z-position at which it intersects the z=0 plane
608 // with the transform applied. Once we have that point we apply the
609 // inverse transform to find the corresponding point in the source
610 // space.
611 //
612 // Given a plane with normal Pn, and a ray starting at point R0 and
613 // with direction defined by the vector Rd, we can find the
614 // intersection point as a distance d from R0 in units of Rd by:
615 //
616 // d = -dot (Pn', R0) / dot (Pn', Rd)
617 if (clamped)
618 *clamped = false;
619
620 if (m33() == 0) {
621 // In this case, the projection plane is parallel to the ray we are trying to
622 // trace, and there is no well-defined value for the projection.
623 return FloatPoint();
624 }
625
626 double x = p.x();
627 double y = p.y();
628 double z = -(m13() * x + m23() * y + m43()) / m33();
629
630 // FIXME: use multVecMatrix()
631 double outX = x * m11() + y * m21() + z * m31() + m41();
632 double outY = x * m12() + y * m22() + z * m32() + m42();
633
634 double w = x * m14() + y * m24() + z * m34() + m44();
635 if (w <= 0) {
636 // Using int max causes overflow when other code uses the projected point. To
637 // represent infinity yet reduce the risk of overflow, we use a large but
638 // not-too-large number here when clamping.
639 const int largeNumber = 100000000 / kFixedPointDenominator;
640 outX = copysign(largeNumber, outX);
641 outY = copysign(largeNumber, outY);
642 if (clamped)
643 *clamped = true;
644 } else if (w != 1) {
645 outX /= w;
646 outY /= w;
647 }
648
649 return FloatPoint(static_cast<float>(outX), static_cast<float>(outY));
650}
651
652FloatQuad TransformationMatrix::projectQuad(const FloatQuad& q, bool* clamped) const
653{
654 FloatQuad projectedQuad;
655
656 bool clamped1 = false;
657 bool clamped2 = false;
658 bool clamped3 = false;
659 bool clamped4 = false;
660
661 projectedQuad.setP1(projectPoint(q.p1(), &clamped1));
662 projectedQuad.setP2(projectPoint(q.p2(), &clamped2));
663 projectedQuad.setP3(projectPoint(q.p3(), &clamped3));
664 projectedQuad.setP4(projectPoint(q.p4(), &clamped4));
665
666 if (clamped)
667 *clamped = clamped1 || clamped2 || clamped3 || clamped4;
668
669 // If all points on the quad had w < 0, then the entire quad would not be visible to the projected surface.
670 bool everythingWasClipped = clamped1 && clamped2 && clamped3 && clamped4;
671 if (everythingWasClipped)
672 return FloatQuad();
673
674 return projectedQuad;
675}
676
677static float clampEdgeValue(float f)
678{
679 ASSERT(!std::isnan(f));
680 return std::min<float>(std::max<float>(f, -LayoutUnit::max() / 2), LayoutUnit::max() / 2);
681}
682
683LayoutRect TransformationMatrix::clampedBoundsOfProjectedQuad(const FloatQuad& q) const
684{
685 FloatRect mappedQuadBounds = projectQuad(q).boundingBox();
686
687 float left = clampEdgeValue(floorf(mappedQuadBounds.x()));
688 float top = clampEdgeValue(floorf(mappedQuadBounds.y()));
689
690 float right;
691 if (std::isinf(mappedQuadBounds.x()) && std::isinf(mappedQuadBounds.width()))
692 right = LayoutUnit::max() / 2;
693 else
694 right = clampEdgeValue(ceilf(mappedQuadBounds.maxX()));
695
696 float bottom;
697 if (std::isinf(mappedQuadBounds.y()) && std::isinf(mappedQuadBounds.height()))
698 bottom = LayoutUnit::max() / 2;
699 else
700 bottom = clampEdgeValue(ceilf(mappedQuadBounds.maxY()));
701
702 return LayoutRect(LayoutUnit::clamp(left), LayoutUnit::clamp(top), LayoutUnit::clamp(right - left), LayoutUnit::clamp(bottom - top));
703}
704
705void TransformationMatrix::map4ComponentPoint(double& x, double& y, double& z, double& w) const
706{
707 if (isIdentityOrTranslation()) {
708 x += m_matrix[3][0];
709 y += m_matrix[3][1];
710 z += m_matrix[3][2];
711 return;
712 }
713
714 Vector4 input = { x, y, z, w };
715 Vector4 result;
716 v4MulPointByMatrix(input, m_matrix, result);
717
718 x = result[0];
719 y = result[1];
720 z = result[2];
721 w = result[3];
722}
723
724FloatPoint TransformationMatrix::mapPoint(const FloatPoint& p) const
725{
726 if (isIdentityOrTranslation())
727 return FloatPoint(p.x() + static_cast<float>(m_matrix[3][0]), p.y() + static_cast<float>(m_matrix[3][1]));
728
729 return internalMapPoint(p);
730}
731
732FloatPoint3D TransformationMatrix::mapPoint(const FloatPoint3D& p) const
733{
734 if (isIdentityOrTranslation())
735 return FloatPoint3D(p.x() + static_cast<float>(m_matrix[3][0]),
736 p.y() + static_cast<float>(m_matrix[3][1]),
737 p.z() + static_cast<float>(m_matrix[3][2]));
738
739 return internalMapPoint(p);
740}
741
742IntRect TransformationMatrix::mapRect(const IntRect &rect) const
743{
744 return enclosingIntRect(mapRect(FloatRect(rect)));
745}
746
747LayoutRect TransformationMatrix::mapRect(const LayoutRect& r) const
748{
749 return enclosingLayoutRect(mapRect(FloatRect(r)));
750}
751
752FloatRect TransformationMatrix::mapRect(const FloatRect& r) const
753{
754 if (isIdentityOrTranslation()) {
755 FloatRect mappedRect(r);
756 mappedRect.move(static_cast<float>(m_matrix[3][0]), static_cast<float>(m_matrix[3][1]));
757 return mappedRect;
758 }
759
760 FloatQuad result;
761
762 float maxX = r.maxX();
763 float maxY = r.maxY();
764 result.setP1(internalMapPoint(FloatPoint(r.x(), r.y())));
765 result.setP2(internalMapPoint(FloatPoint(maxX, r.y())));
766 result.setP3(internalMapPoint(FloatPoint(maxX, maxY)));
767 result.setP4(internalMapPoint(FloatPoint(r.x(), maxY)));
768
769 return result.boundingBox();
770}
771
772FloatQuad TransformationMatrix::mapQuad(const FloatQuad& q) const
773{
774 if (isIdentityOrTranslation()) {
775 FloatQuad mappedQuad(q);
776 mappedQuad.move(static_cast<float>(m_matrix[3][0]), static_cast<float>(m_matrix[3][1]));
777 return mappedQuad;
778 }
779
780 FloatQuad result;
781 result.setP1(internalMapPoint(q.p1()));
782 result.setP2(internalMapPoint(q.p2()));
783 result.setP3(internalMapPoint(q.p3()));
784 result.setP4(internalMapPoint(q.p4()));
785 return result;
786}
787
788TransformationMatrix& TransformationMatrix::scaleNonUniform(double sx, double sy)
789{
790 m_matrix[0][0] *= sx;
791 m_matrix[0][1] *= sx;
792 m_matrix[0][2] *= sx;
793 m_matrix[0][3] *= sx;
794
795 m_matrix[1][0] *= sy;
796 m_matrix[1][1] *= sy;
797 m_matrix[1][2] *= sy;
798 m_matrix[1][3] *= sy;
799 return *this;
800}
801
802TransformationMatrix& TransformationMatrix::scale3d(double sx, double sy, double sz)
803{
804 scaleNonUniform(sx, sy);
805
806 m_matrix[2][0] *= sz;
807 m_matrix[2][1] *= sz;
808 m_matrix[2][2] *= sz;
809 m_matrix[2][3] *= sz;
810 return *this;
811}
812
813TransformationMatrix& TransformationMatrix::rotate3d(double x, double y, double z, double angle)
814{
815 // Normalize the axis of rotation
816 double length = sqrt(x * x + y * y + z * z);
817 if (length == 0) {
818 // A direction vector that cannot be normalized, such as [0, 0, 0], will cause the rotation to not be applied.
819 return *this;
820 } else if (length != 1) {
821 x /= length;
822 y /= length;
823 z /= length;
824 }
825
826 // Angles are in degrees. Switch to radians.
827 angle = deg2rad(angle);
828
829 double sinTheta = sin(angle);
830 double cosTheta = cos(angle);
831
832 TransformationMatrix mat;
833
834 // Optimize cases where the axis is along a major axis
835 if (x == 1.0 && y == 0.0 && z == 0.0) {
836 mat.m_matrix[0][0] = 1.0;
837 mat.m_matrix[0][1] = 0.0;
838 mat.m_matrix[0][2] = 0.0;
839 mat.m_matrix[1][0] = 0.0;
840 mat.m_matrix[1][1] = cosTheta;
841 mat.m_matrix[1][2] = sinTheta;
842 mat.m_matrix[2][0] = 0.0;
843 mat.m_matrix[2][1] = -sinTheta;
844 mat.m_matrix[2][2] = cosTheta;
845 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
846 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
847 mat.m_matrix[3][3] = 1.0;
848 } else if (x == 0.0 && y == 1.0 && z == 0.0) {
849 mat.m_matrix[0][0] = cosTheta;
850 mat.m_matrix[0][1] = 0.0;
851 mat.m_matrix[0][2] = -sinTheta;
852 mat.m_matrix[1][0] = 0.0;
853 mat.m_matrix[1][1] = 1.0;
854 mat.m_matrix[1][2] = 0.0;
855 mat.m_matrix[2][0] = sinTheta;
856 mat.m_matrix[2][1] = 0.0;
857 mat.m_matrix[2][2] = cosTheta;
858 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
859 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
860 mat.m_matrix[3][3] = 1.0;
861 } else if (x == 0.0 && y == 0.0 && z == 1.0) {
862 mat.m_matrix[0][0] = cosTheta;
863 mat.m_matrix[0][1] = sinTheta;
864 mat.m_matrix[0][2] = 0.0;
865 mat.m_matrix[1][0] = -sinTheta;
866 mat.m_matrix[1][1] = cosTheta;
867 mat.m_matrix[1][2] = 0.0;
868 mat.m_matrix[2][0] = 0.0;
869 mat.m_matrix[2][1] = 0.0;
870 mat.m_matrix[2][2] = 1.0;
871 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
872 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
873 mat.m_matrix[3][3] = 1.0;
874 } else {
875 // This case is the rotation about an arbitrary unit vector.
876 //
877 // Formula is adapted from Wikipedia article on Rotation matrix,
878 // http://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
879 //
880 // An alternate resource with the same matrix: http://www.fastgraph.com/makegames/3drotation/
881 //
882 double oneMinusCosTheta = 1 - cosTheta;
883 mat.m_matrix[0][0] = cosTheta + x * x * oneMinusCosTheta;
884 mat.m_matrix[0][1] = y * x * oneMinusCosTheta + z * sinTheta;
885 mat.m_matrix[0][2] = z * x * oneMinusCosTheta - y * sinTheta;
886 mat.m_matrix[1][0] = x * y * oneMinusCosTheta - z * sinTheta;
887 mat.m_matrix[1][1] = cosTheta + y * y * oneMinusCosTheta;
888 mat.m_matrix[1][2] = z * y * oneMinusCosTheta + x * sinTheta;
889 mat.m_matrix[2][0] = x * z * oneMinusCosTheta + y * sinTheta;
890 mat.m_matrix[2][1] = y * z * oneMinusCosTheta - x * sinTheta;
891 mat.m_matrix[2][2] = cosTheta + z * z * oneMinusCosTheta;
892 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
893 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
894 mat.m_matrix[3][3] = 1.0;
895 }
896 multiply(mat);
897 return *this;
898}
899
900TransformationMatrix& TransformationMatrix::rotate3d(double rx, double ry, double rz)
901{
902 // Angles are in degrees. Switch to radians.
903 rx = deg2rad(rx);
904 ry = deg2rad(ry);
905 rz = deg2rad(rz);
906
907 TransformationMatrix mat;
908
909 double sinTheta = sin(rz);
910 double cosTheta = cos(rz);
911
912 mat.m_matrix[0][0] = cosTheta;
913 mat.m_matrix[0][1] = sinTheta;
914 mat.m_matrix[0][2] = 0.0;
915 mat.m_matrix[1][0] = -sinTheta;
916 mat.m_matrix[1][1] = cosTheta;
917 mat.m_matrix[1][2] = 0.0;
918 mat.m_matrix[2][0] = 0.0;
919 mat.m_matrix[2][1] = 0.0;
920 mat.m_matrix[2][2] = 1.0;
921 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
922 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
923 mat.m_matrix[3][3] = 1.0;
924
925 TransformationMatrix rmat(mat);
926
927 sinTheta = sin(ry);
928 cosTheta = cos(ry);
929
930 mat.m_matrix[0][0] = cosTheta;
931 mat.m_matrix[0][1] = 0.0;
932 mat.m_matrix[0][2] = -sinTheta;
933 mat.m_matrix[1][0] = 0.0;
934 mat.m_matrix[1][1] = 1.0;
935 mat.m_matrix[1][2] = 0.0;
936 mat.m_matrix[2][0] = sinTheta;
937 mat.m_matrix[2][1] = 0.0;
938 mat.m_matrix[2][2] = cosTheta;
939 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
940 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
941 mat.m_matrix[3][3] = 1.0;
942
943 rmat.multiply(mat);
944
945 sinTheta = sin(rx);
946 cosTheta = cos(rx);
947
948 mat.m_matrix[0][0] = 1.0;
949 mat.m_matrix[0][1] = 0.0;
950 mat.m_matrix[0][2] = 0.0;
951 mat.m_matrix[1][0] = 0.0;
952 mat.m_matrix[1][1] = cosTheta;
953 mat.m_matrix[1][2] = sinTheta;
954 mat.m_matrix[2][0] = 0.0;
955 mat.m_matrix[2][1] = -sinTheta;
956 mat.m_matrix[2][2] = cosTheta;
957 mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0;
958 mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0;
959 mat.m_matrix[3][3] = 1.0;
960
961 rmat.multiply(mat);
962
963 multiply(rmat);
964 return *this;
965}
966
967TransformationMatrix& TransformationMatrix::translate(double tx, double ty)
968{
969 m_matrix[3][0] += tx * m_matrix[0][0] + ty * m_matrix[1][0];
970 m_matrix[3][1] += tx * m_matrix[0][1] + ty * m_matrix[1][1];
971 m_matrix[3][2] += tx * m_matrix[0][2] + ty * m_matrix[1][2];
972 m_matrix[3][3] += tx * m_matrix[0][3] + ty * m_matrix[1][3];
973 return *this;
974}
975
976TransformationMatrix& TransformationMatrix::translate3d(double tx, double ty, double tz)
977{
978 m_matrix[3][0] += tx * m_matrix[0][0] + ty * m_matrix[1][0] + tz * m_matrix[2][0];
979 m_matrix[3][1] += tx * m_matrix[0][1] + ty * m_matrix[1][1] + tz * m_matrix[2][1];
980 m_matrix[3][2] += tx * m_matrix[0][2] + ty * m_matrix[1][2] + tz * m_matrix[2][2];
981 m_matrix[3][3] += tx * m_matrix[0][3] + ty * m_matrix[1][3] + tz * m_matrix[2][3];
982 return *this;
983}
984
985TransformationMatrix& TransformationMatrix::translateRight(double tx, double ty)
986{
987 if (tx != 0) {
988 m_matrix[0][0] += m_matrix[0][3] * tx;
989 m_matrix[1][0] += m_matrix[1][3] * tx;
990 m_matrix[2][0] += m_matrix[2][3] * tx;
991 m_matrix[3][0] += m_matrix[3][3] * tx;
992 }
993
994 if (ty != 0) {
995 m_matrix[0][1] += m_matrix[0][3] * ty;
996 m_matrix[1][1] += m_matrix[1][3] * ty;
997 m_matrix[2][1] += m_matrix[2][3] * ty;
998 m_matrix[3][1] += m_matrix[3][3] * ty;
999 }
1000
1001 return *this;
1002}
1003
1004TransformationMatrix& TransformationMatrix::translateRight3d(double tx, double ty, double tz)
1005{
1006 translateRight(tx, ty);
1007 if (tz != 0) {
1008 m_matrix[0][2] += m_matrix[0][3] * tz;
1009 m_matrix[1][2] += m_matrix[1][3] * tz;
1010 m_matrix[2][2] += m_matrix[2][3] * tz;
1011 m_matrix[3][2] += m_matrix[3][3] * tz;
1012 }
1013
1014 return *this;
1015}
1016
1017TransformationMatrix& TransformationMatrix::skew(double sx, double sy)
1018{
1019 // angles are in degrees. Switch to radians
1020 sx = deg2rad(sx);
1021 sy = deg2rad(sy);
1022
1023 TransformationMatrix mat;
1024 mat.m_matrix[0][1] = tan(sy); // note that the y shear goes in the first row
1025 mat.m_matrix[1][0] = tan(sx); // and the x shear in the second row
1026
1027 multiply(mat);
1028 return *this;
1029}
1030
1031TransformationMatrix& TransformationMatrix::applyPerspective(double p)
1032{
1033 TransformationMatrix mat;
1034 if (p != 0)
1035 mat.m_matrix[2][3] = -1/p;
1036
1037 multiply(mat);
1038 return *this;
1039}
1040
1041TransformationMatrix TransformationMatrix::rectToRect(const FloatRect& from, const FloatRect& to)
1042{
1043 ASSERT(!from.isEmpty());
1044 return TransformationMatrix(to.width() / from.width(),
1045 0, 0,
1046 to.height() / from.height(),
1047 to.x() - from.x(),
1048 to.y() - from.y());
1049}
1050
1051// this = mat * this.
1052TransformationMatrix& TransformationMatrix::multiply(const TransformationMatrix& mat)
1053{
1054#if CPU(ARM64) && defined(_LP64)
1055 double* leftMatrix = &(m_matrix[0][0]);
1056 const double* rightMatrix = &(mat.m_matrix[0][0]);
1057 asm volatile (
1058 // First, load the leftMatrix completely in memory. The leftMatrix is in v16-v23.
1059 "mov x4, %[leftMatrix]\n\t"
1060 "ld1 {v16.2d, v17.2d, v18.2d, v19.2d}, [%[leftMatrix]], #64\n\t"
1061 "ld1 {v20.2d, v21.2d, v22.2d, v23.2d}, [%[leftMatrix]]\n\t"
1062
1063 // First row.
1064 "ld4r {v24.2d, v25.2d, v26.2d, v27.2d}, [%[rightMatrix]], #32\n\t"
1065 "fmul v28.2d, v24.2d, v16.2d\n\t"
1066 "fmul v29.2d, v24.2d, v17.2d\n\t"
1067 "fmla v28.2d, v25.2d, v18.2d\n\t"
1068 "fmla v29.2d, v25.2d, v19.2d\n\t"
1069 "fmla v28.2d, v26.2d, v20.2d\n\t"
1070 "fmla v29.2d, v26.2d, v21.2d\n\t"
1071 "fmla v28.2d, v27.2d, v22.2d\n\t"
1072 "fmla v29.2d, v27.2d, v23.2d\n\t"
1073
1074 "ld4r {v0.2d, v1.2d, v2.2d, v3.2d}, [%[rightMatrix]], #32\n\t"
1075 "st1 {v28.2d, v29.2d}, [x4], #32\n\t"
1076
1077 // Second row.
1078 "fmul v30.2d, v0.2d, v16.2d\n\t"
1079 "fmul v31.2d, v0.2d, v17.2d\n\t"
1080 "fmla v30.2d, v1.2d, v18.2d\n\t"
1081 "fmla v31.2d, v1.2d, v19.2d\n\t"
1082 "fmla v30.2d, v2.2d, v20.2d\n\t"
1083 "fmla v31.2d, v2.2d, v21.2d\n\t"
1084 "fmla v30.2d, v3.2d, v22.2d\n\t"
1085 "fmla v31.2d, v3.2d, v23.2d\n\t"
1086
1087 "ld4r {v24.2d, v25.2d, v26.2d, v27.2d}, [%[rightMatrix]], #32\n\t"
1088 "st1 {v30.2d, v31.2d}, [x4], #32\n\t"
1089
1090 // Third row.
1091 "fmul v28.2d, v24.2d, v16.2d\n\t"
1092 "fmul v29.2d, v24.2d, v17.2d\n\t"
1093 "fmla v28.2d, v25.2d, v18.2d\n\t"
1094 "fmla v29.2d, v25.2d, v19.2d\n\t"
1095 "fmla v28.2d, v26.2d, v20.2d\n\t"
1096 "fmla v29.2d, v26.2d, v21.2d\n\t"
1097 "fmla v28.2d, v27.2d, v22.2d\n\t"
1098 "fmla v29.2d, v27.2d, v23.2d\n\t"
1099
1100 "ld4r {v0.2d, v1.2d, v2.2d, v3.2d}, [%[rightMatrix]], #32\n\t"
1101 "st1 {v28.2d, v29.2d}, [x4], #32\n\t"
1102
1103 // Fourth row.
1104 "fmul v30.2d, v0.2d, v16.2d\n\t"
1105 "fmul v31.2d, v0.2d, v17.2d\n\t"
1106 "fmla v30.2d, v1.2d, v18.2d\n\t"
1107 "fmla v31.2d, v1.2d, v19.2d\n\t"
1108 "fmla v30.2d, v2.2d, v20.2d\n\t"
1109 "fmla v31.2d, v2.2d, v21.2d\n\t"
1110 "fmla v30.2d, v3.2d, v22.2d\n\t"
1111 "fmla v31.2d, v3.2d, v23.2d\n\t"
1112
1113 "st1 {v30.2d, v31.2d}, [x4]\n\t"
1114
1115 : [leftMatrix]"+r"(leftMatrix), [rightMatrix]"+r"(rightMatrix)
1116 :
1117 : "memory", "x4", "v0", "v1", "v2", "v3", "v16", "v17", "v18", "v19", "v20", "v21", "v22", "v23", "v24", "v25", "v26", "v27", "v28", "v29", "v30", "v31");
1118#elif CPU(APPLE_ARMV7S)
1119 double* leftMatrix = &(m_matrix[0][0]);
1120 const double* rightMatrix = &(mat.m_matrix[0][0]);
1121 asm volatile (// First row of leftMatrix.
1122 "mov r3, %[leftMatrix]\n\t"
1123 "vld1.64 { d16-d19 }, [%[leftMatrix], :128]!\n\t"
1124 "vld1.64 { d0-d3}, [%[rightMatrix], :128]!\n\t"
1125 "vmul.f64 d4, d0, d16\n\t"
1126 "vld1.64 { d20-d23 }, [%[leftMatrix], :128]!\n\t"
1127 "vmla.f64 d4, d1, d20\n\t"
1128 "vld1.64 { d24-d27 }, [%[leftMatrix], :128]!\n\t"
1129 "vmla.f64 d4, d2, d24\n\t"
1130 "vld1.64 { d28-d31 }, [%[leftMatrix], :128]!\n\t"
1131 "vmla.f64 d4, d3, d28\n\t"
1132
1133 "vmul.f64 d5, d0, d17\n\t"
1134 "vmla.f64 d5, d1, d21\n\t"
1135 "vmla.f64 d5, d2, d25\n\t"
1136 "vmla.f64 d5, d3, d29\n\t"
1137
1138 "vmul.f64 d6, d0, d18\n\t"
1139 "vmla.f64 d6, d1, d22\n\t"
1140 "vmla.f64 d6, d2, d26\n\t"
1141 "vmla.f64 d6, d3, d30\n\t"
1142
1143 "vmul.f64 d7, d0, d19\n\t"
1144 "vmla.f64 d7, d1, d23\n\t"
1145 "vmla.f64 d7, d2, d27\n\t"
1146 "vmla.f64 d7, d3, d31\n\t"
1147 "vld1.64 { d0-d3}, [%[rightMatrix], :128]!\n\t"
1148 "vst1.64 { d4-d7 }, [r3, :128]!\n\t"
1149
1150 // Second row of leftMatrix.
1151 "vmul.f64 d4, d0, d16\n\t"
1152 "vmla.f64 d4, d1, d20\n\t"
1153 "vmla.f64 d4, d2, d24\n\t"
1154 "vmla.f64 d4, d3, d28\n\t"
1155
1156 "vmul.f64 d5, d0, d17\n\t"
1157 "vmla.f64 d5, d1, d21\n\t"
1158 "vmla.f64 d5, d2, d25\n\t"
1159 "vmla.f64 d5, d3, d29\n\t"
1160
1161 "vmul.f64 d6, d0, d18\n\t"
1162 "vmla.f64 d6, d1, d22\n\t"
1163 "vmla.f64 d6, d2, d26\n\t"
1164 "vmla.f64 d6, d3, d30\n\t"
1165
1166 "vmul.f64 d7, d0, d19\n\t"
1167 "vmla.f64 d7, d1, d23\n\t"
1168 "vmla.f64 d7, d2, d27\n\t"
1169 "vmla.f64 d7, d3, d31\n\t"
1170 "vld1.64 { d0-d3}, [%[rightMatrix], :128]!\n\t"
1171 "vst1.64 { d4-d7 }, [r3, :128]!\n\t"
1172
1173 // Third row of leftMatrix.
1174 "vmul.f64 d4, d0, d16\n\t"
1175 "vmla.f64 d4, d1, d20\n\t"
1176 "vmla.f64 d4, d2, d24\n\t"
1177 "vmla.f64 d4, d3, d28\n\t"
1178
1179 "vmul.f64 d5, d0, d17\n\t"
1180 "vmla.f64 d5, d1, d21\n\t"
1181 "vmla.f64 d5, d2, d25\n\t"
1182 "vmla.f64 d5, d3, d29\n\t"
1183
1184 "vmul.f64 d6, d0, d18\n\t"
1185 "vmla.f64 d6, d1, d22\n\t"
1186 "vmla.f64 d6, d2, d26\n\t"
1187 "vmla.f64 d6, d3, d30\n\t"
1188
1189 "vmul.f64 d7, d0, d19\n\t"
1190 "vmla.f64 d7, d1, d23\n\t"
1191 "vmla.f64 d7, d2, d27\n\t"
1192 "vmla.f64 d7, d3, d31\n\t"
1193 "vld1.64 { d0-d3}, [%[rightMatrix], :128]\n\t"
1194 "vst1.64 { d4-d7 }, [r3, :128]!\n\t"
1195
1196 // Fourth and last row of leftMatrix.
1197 "vmul.f64 d4, d0, d16\n\t"
1198 "vmla.f64 d4, d1, d20\n\t"
1199 "vmla.f64 d4, d2, d24\n\t"
1200 "vmla.f64 d4, d3, d28\n\t"
1201
1202 "vmul.f64 d5, d0, d17\n\t"
1203 "vmla.f64 d5, d1, d21\n\t"
1204 "vmla.f64 d5, d2, d25\n\t"
1205 "vmla.f64 d5, d3, d29\n\t"
1206
1207 "vmul.f64 d6, d0, d18\n\t"
1208 "vmla.f64 d6, d1, d22\n\t"
1209 "vmla.f64 d6, d2, d26\n\t"
1210 "vmla.f64 d6, d3, d30\n\t"
1211
1212 "vmul.f64 d7, d0, d19\n\t"
1213 "vmla.f64 d7, d1, d23\n\t"
1214 "vmla.f64 d7, d2, d27\n\t"
1215 "vmla.f64 d7, d3, d31\n\t"
1216 "vst1.64 { d4-d7 }, [r3, :128]\n\t"
1217 : [leftMatrix]"+r"(leftMatrix), [rightMatrix]"+r"(rightMatrix)
1218 :
1219 : "memory", "r3", "d0", "d1", "d2", "d3", "d4", "d5", "d6", "d7", "d16", "d17", "d18", "d19", "d20", "d21", "d22", "d23", "d24", "d25", "d26", "d27", "d28", "d29", "d30", "d31");
1220#elif CPU(ARM_VFP) && PLATFORM(IOS_FAMILY)
1221
1222#define MATRIX_MULTIPLY_ONE_LINE \
1223 "vldmia.64 %[rightMatrix]!, { d0-d3}\n\t" \
1224 "vmul.f64 d4, d0, d16\n\t" \
1225 "vmla.f64 d4, d1, d20\n\t" \
1226 "vmla.f64 d4, d2, d24\n\t" \
1227 "vmla.f64 d4, d3, d28\n\t" \
1228 \
1229 "vmul.f64 d5, d0, d17\n\t" \
1230 "vmla.f64 d5, d1, d21\n\t" \
1231 "vmla.f64 d5, d2, d25\n\t" \
1232 "vmla.f64 d5, d3, d29\n\t" \
1233 \
1234 "vmul.f64 d6, d0, d18\n\t" \
1235 "vmla.f64 d6, d1, d22\n\t" \
1236 "vmla.f64 d6, d2, d26\n\t" \
1237 "vmla.f64 d6, d3, d30\n\t" \
1238 \
1239 "vmul.f64 d7, d0, d19\n\t" \
1240 "vmla.f64 d7, d1, d23\n\t" \
1241 "vmla.f64 d7, d2, d27\n\t" \
1242 "vmla.f64 d7, d3, d31\n\t" \
1243 "vstmia.64 %[leftMatrix]!, { d4-d7 }\n\t"
1244
1245 double* leftMatrix = &(m_matrix[0][0]);
1246 const double* rightMatrix = &(mat.m_matrix[0][0]);
1247 // We load the full m_matrix at once in d16-d31.
1248 asm volatile("vldmia.64 %[leftMatrix], { d16-d31 }\n\t"
1249 :
1250 : [leftMatrix]"r"(leftMatrix)
1251 : "d16", "d17", "d18", "d19", "d20", "d21", "d22", "d23", "d24", "d25", "d26", "d27", "d28", "d29", "d30", "d31");
1252 for (unsigned i = 0; i < 4; ++i) {
1253 asm volatile(MATRIX_MULTIPLY_ONE_LINE
1254 : [leftMatrix]"+r"(leftMatrix), [rightMatrix]"+r"(rightMatrix)
1255 :
1256 : "memory", "d0", "d1", "d2", "d3", "d4", "d5", "d6", "d7");
1257 }
1258#undef MATRIX_MULTIPLY_ONE_LINE
1259
1260#elif defined(TRANSFORMATION_MATRIX_USE_X86_64_SSE2)
1261 // x86_64 has 16 XMM registers which is enough to do the multiplication fully in registers.
1262 __m128d matrixBlockA = _mm_load_pd(&(m_matrix[0][0]));
1263 __m128d matrixBlockC = _mm_load_pd(&(m_matrix[1][0]));
1264 __m128d matrixBlockE = _mm_load_pd(&(m_matrix[2][0]));
1265 __m128d matrixBlockG = _mm_load_pd(&(m_matrix[3][0]));
1266
1267 // First row.
1268 __m128d otherMatrixFirstParam = _mm_set1_pd(mat.m_matrix[0][0]);
1269 __m128d otherMatrixSecondParam = _mm_set1_pd(mat.m_matrix[0][1]);
1270 __m128d otherMatrixThirdParam = _mm_set1_pd(mat.m_matrix[0][2]);
1271 __m128d otherMatrixFourthParam = _mm_set1_pd(mat.m_matrix[0][3]);
1272
1273 // output00 and output01.
1274 __m128d accumulator = _mm_mul_pd(matrixBlockA, otherMatrixFirstParam);
1275 __m128d temp1 = _mm_mul_pd(matrixBlockC, otherMatrixSecondParam);
1276 __m128d temp2 = _mm_mul_pd(matrixBlockE, otherMatrixThirdParam);
1277 __m128d temp3 = _mm_mul_pd(matrixBlockG, otherMatrixFourthParam);
1278
1279 __m128d matrixBlockB = _mm_load_pd(&(m_matrix[0][2]));
1280 __m128d matrixBlockD = _mm_load_pd(&(m_matrix[1][2]));
1281 __m128d matrixBlockF = _mm_load_pd(&(m_matrix[2][2]));
1282 __m128d matrixBlockH = _mm_load_pd(&(m_matrix[3][2]));
1283
1284 accumulator = _mm_add_pd(accumulator, temp1);
1285 accumulator = _mm_add_pd(accumulator, temp2);
1286 accumulator = _mm_add_pd(accumulator, temp3);
1287 _mm_store_pd(&m_matrix[0][0], accumulator);
1288
1289 // output02 and output03.
1290 accumulator = _mm_mul_pd(matrixBlockB, otherMatrixFirstParam);
1291 temp1 = _mm_mul_pd(matrixBlockD, otherMatrixSecondParam);
1292 temp2 = _mm_mul_pd(matrixBlockF, otherMatrixThirdParam);
1293 temp3 = _mm_mul_pd(matrixBlockH, otherMatrixFourthParam);
1294
1295 accumulator = _mm_add_pd(accumulator, temp1);
1296 accumulator = _mm_add_pd(accumulator, temp2);
1297 accumulator = _mm_add_pd(accumulator, temp3);
1298 _mm_store_pd(&m_matrix[0][2], accumulator);
1299
1300 // Second row.
1301 otherMatrixFirstParam = _mm_set1_pd(mat.m_matrix[1][0]);
1302 otherMatrixSecondParam = _mm_set1_pd(mat.m_matrix[1][1]);
1303 otherMatrixThirdParam = _mm_set1_pd(mat.m_matrix[1][2]);
1304 otherMatrixFourthParam = _mm_set1_pd(mat.m_matrix[1][3]);
1305
1306 // output10 and output11.
1307 accumulator = _mm_mul_pd(matrixBlockA, otherMatrixFirstParam);
1308 temp1 = _mm_mul_pd(matrixBlockC, otherMatrixSecondParam);
1309 temp2 = _mm_mul_pd(matrixBlockE, otherMatrixThirdParam);
1310 temp3 = _mm_mul_pd(matrixBlockG, otherMatrixFourthParam);
1311
1312 accumulator = _mm_add_pd(accumulator, temp1);
1313 accumulator = _mm_add_pd(accumulator, temp2);
1314 accumulator = _mm_add_pd(accumulator, temp3);
1315 _mm_store_pd(&m_matrix[1][0], accumulator);
1316
1317 // output12 and output13.
1318 accumulator = _mm_mul_pd(matrixBlockB, otherMatrixFirstParam);
1319 temp1 = _mm_mul_pd(matrixBlockD, otherMatrixSecondParam);
1320 temp2 = _mm_mul_pd(matrixBlockF, otherMatrixThirdParam);
1321 temp3 = _mm_mul_pd(matrixBlockH, otherMatrixFourthParam);
1322
1323 accumulator = _mm_add_pd(accumulator, temp1);
1324 accumulator = _mm_add_pd(accumulator, temp2);
1325 accumulator = _mm_add_pd(accumulator, temp3);
1326 _mm_store_pd(&m_matrix[1][2], accumulator);
1327
1328 // Third row.
1329 otherMatrixFirstParam = _mm_set1_pd(mat.m_matrix[2][0]);
1330 otherMatrixSecondParam = _mm_set1_pd(mat.m_matrix[2][1]);
1331 otherMatrixThirdParam = _mm_set1_pd(mat.m_matrix[2][2]);
1332 otherMatrixFourthParam = _mm_set1_pd(mat.m_matrix[2][3]);
1333
1334 // output20 and output21.
1335 accumulator = _mm_mul_pd(matrixBlockA, otherMatrixFirstParam);
1336 temp1 = _mm_mul_pd(matrixBlockC, otherMatrixSecondParam);
1337 temp2 = _mm_mul_pd(matrixBlockE, otherMatrixThirdParam);
1338 temp3 = _mm_mul_pd(matrixBlockG, otherMatrixFourthParam);
1339
1340 accumulator = _mm_add_pd(accumulator, temp1);
1341 accumulator = _mm_add_pd(accumulator, temp2);
1342 accumulator = _mm_add_pd(accumulator, temp3);
1343 _mm_store_pd(&m_matrix[2][0], accumulator);
1344
1345 // output22 and output23.
1346 accumulator = _mm_mul_pd(matrixBlockB, otherMatrixFirstParam);
1347 temp1 = _mm_mul_pd(matrixBlockD, otherMatrixSecondParam);
1348 temp2 = _mm_mul_pd(matrixBlockF, otherMatrixThirdParam);
1349 temp3 = _mm_mul_pd(matrixBlockH, otherMatrixFourthParam);
1350
1351 accumulator = _mm_add_pd(accumulator, temp1);
1352 accumulator = _mm_add_pd(accumulator, temp2);
1353 accumulator = _mm_add_pd(accumulator, temp3);
1354 _mm_store_pd(&m_matrix[2][2], accumulator);
1355
1356 // Fourth row.
1357 otherMatrixFirstParam = _mm_set1_pd(mat.m_matrix[3][0]);
1358 otherMatrixSecondParam = _mm_set1_pd(mat.m_matrix[3][1]);
1359 otherMatrixThirdParam = _mm_set1_pd(mat.m_matrix[3][2]);
1360 otherMatrixFourthParam = _mm_set1_pd(mat.m_matrix[3][3]);
1361
1362 // output30 and output31.
1363 accumulator = _mm_mul_pd(matrixBlockA, otherMatrixFirstParam);
1364 temp1 = _mm_mul_pd(matrixBlockC, otherMatrixSecondParam);
1365 temp2 = _mm_mul_pd(matrixBlockE, otherMatrixThirdParam);
1366 temp3 = _mm_mul_pd(matrixBlockG, otherMatrixFourthParam);
1367
1368 accumulator = _mm_add_pd(accumulator, temp1);
1369 accumulator = _mm_add_pd(accumulator, temp2);
1370 accumulator = _mm_add_pd(accumulator, temp3);
1371 _mm_store_pd(&m_matrix[3][0], accumulator);
1372
1373 // output32 and output33.
1374 accumulator = _mm_mul_pd(matrixBlockB, otherMatrixFirstParam);
1375 temp1 = _mm_mul_pd(matrixBlockD, otherMatrixSecondParam);
1376 temp2 = _mm_mul_pd(matrixBlockF, otherMatrixThirdParam);
1377 temp3 = _mm_mul_pd(matrixBlockH, otherMatrixFourthParam);
1378
1379 accumulator = _mm_add_pd(accumulator, temp1);
1380 accumulator = _mm_add_pd(accumulator, temp2);
1381 accumulator = _mm_add_pd(accumulator, temp3);
1382 _mm_store_pd(&m_matrix[3][2], accumulator);
1383#else
1384 Matrix4 tmp;
1385
1386 tmp[0][0] = (mat.m_matrix[0][0] * m_matrix[0][0] + mat.m_matrix[0][1] * m_matrix[1][0]
1387 + mat.m_matrix[0][2] * m_matrix[2][0] + mat.m_matrix[0][3] * m_matrix[3][0]);
1388 tmp[0][1] = (mat.m_matrix[0][0] * m_matrix[0][1] + mat.m_matrix[0][1] * m_matrix[1][1]
1389 + mat.m_matrix[0][2] * m_matrix[2][1] + mat.m_matrix[0][3] * m_matrix[3][1]);
1390 tmp[0][2] = (mat.m_matrix[0][0] * m_matrix[0][2] + mat.m_matrix[0][1] * m_matrix[1][2]
1391 + mat.m_matrix[0][2] * m_matrix[2][2] + mat.m_matrix[0][3] * m_matrix[3][2]);
1392 tmp[0][3] = (mat.m_matrix[0][0] * m_matrix[0][3] + mat.m_matrix[0][1] * m_matrix[1][3]
1393 + mat.m_matrix[0][2] * m_matrix[2][3] + mat.m_matrix[0][3] * m_matrix[3][3]);
1394
1395 tmp[1][0] = (mat.m_matrix[1][0] * m_matrix[0][0] + mat.m_matrix[1][1] * m_matrix[1][0]
1396 + mat.m_matrix[1][2] * m_matrix[2][0] + mat.m_matrix[1][3] * m_matrix[3][0]);
1397 tmp[1][1] = (mat.m_matrix[1][0] * m_matrix[0][1] + mat.m_matrix[1][1] * m_matrix[1][1]
1398 + mat.m_matrix[1][2] * m_matrix[2][1] + mat.m_matrix[1][3] * m_matrix[3][1]);
1399 tmp[1][2] = (mat.m_matrix[1][0] * m_matrix[0][2] + mat.m_matrix[1][1] * m_matrix[1][2]
1400 + mat.m_matrix[1][2] * m_matrix[2][2] + mat.m_matrix[1][3] * m_matrix[3][2]);
1401 tmp[1][3] = (mat.m_matrix[1][0] * m_matrix[0][3] + mat.m_matrix[1][1] * m_matrix[1][3]
1402 + mat.m_matrix[1][2] * m_matrix[2][3] + mat.m_matrix[1][3] * m_matrix[3][3]);
1403
1404 tmp[2][0] = (mat.m_matrix[2][0] * m_matrix[0][0] + mat.m_matrix[2][1] * m_matrix[1][0]
1405 + mat.m_matrix[2][2] * m_matrix[2][0] + mat.m_matrix[2][3] * m_matrix[3][0]);
1406 tmp[2][1] = (mat.m_matrix[2][0] * m_matrix[0][1] + mat.m_matrix[2][1] * m_matrix[1][1]
1407 + mat.m_matrix[2][2] * m_matrix[2][1] + mat.m_matrix[2][3] * m_matrix[3][1]);
1408 tmp[2][2] = (mat.m_matrix[2][0] * m_matrix[0][2] + mat.m_matrix[2][1] * m_matrix[1][2]
1409 + mat.m_matrix[2][2] * m_matrix[2][2] + mat.m_matrix[2][3] * m_matrix[3][2]);
1410 tmp[2][3] = (mat.m_matrix[2][0] * m_matrix[0][3] + mat.m_matrix[2][1] * m_matrix[1][3]
1411 + mat.m_matrix[2][2] * m_matrix[2][3] + mat.m_matrix[2][3] * m_matrix[3][3]);
1412
1413 tmp[3][0] = (mat.m_matrix[3][0] * m_matrix[0][0] + mat.m_matrix[3][1] * m_matrix[1][0]
1414 + mat.m_matrix[3][2] * m_matrix[2][0] + mat.m_matrix[3][3] * m_matrix[3][0]);
1415 tmp[3][1] = (mat.m_matrix[3][0] * m_matrix[0][1] + mat.m_matrix[3][1] * m_matrix[1][1]
1416 + mat.m_matrix[3][2] * m_matrix[2][1] + mat.m_matrix[3][3] * m_matrix[3][1]);
1417 tmp[3][2] = (mat.m_matrix[3][0] * m_matrix[0][2] + mat.m_matrix[3][1] * m_matrix[1][2]
1418 + mat.m_matrix[3][2] * m_matrix[2][2] + mat.m_matrix[3][3] * m_matrix[3][2]);
1419 tmp[3][3] = (mat.m_matrix[3][0] * m_matrix[0][3] + mat.m_matrix[3][1] * m_matrix[1][3]
1420 + mat.m_matrix[3][2] * m_matrix[2][3] + mat.m_matrix[3][3] * m_matrix[3][3]);
1421
1422 memcpy(&m_matrix[0][0], &tmp[0][0], sizeof(Matrix4));
1423#endif
1424 return *this;
1425}
1426
1427void TransformationMatrix::multVecMatrix(double x, double y, double& resultX, double& resultY) const
1428{
1429 resultX = m_matrix[3][0] + x * m_matrix[0][0] + y * m_matrix[1][0];
1430 resultY = m_matrix[3][1] + x * m_matrix[0][1] + y * m_matrix[1][1];
1431 double w = m_matrix[3][3] + x * m_matrix[0][3] + y * m_matrix[1][3];
1432 if (w != 1 && w != 0) {
1433 resultX /= w;
1434 resultY /= w;
1435 }
1436}
1437
1438void TransformationMatrix::multVecMatrix(double x, double y, double z, double& resultX, double& resultY, double& resultZ) const
1439{
1440 resultX = m_matrix[3][0] + x * m_matrix[0][0] + y * m_matrix[1][0] + z * m_matrix[2][0];
1441 resultY = m_matrix[3][1] + x * m_matrix[0][1] + y * m_matrix[1][1] + z * m_matrix[2][1];
1442 resultZ = m_matrix[3][2] + x * m_matrix[0][2] + y * m_matrix[1][2] + z * m_matrix[2][2];
1443 double w = m_matrix[3][3] + x * m_matrix[0][3] + y * m_matrix[1][3] + z * m_matrix[2][3];
1444 if (w != 1 && w != 0) {
1445 resultX /= w;
1446 resultY /= w;
1447 resultZ /= w;
1448 }
1449}
1450
1451bool TransformationMatrix::isInvertible() const
1452{
1453 if (isIdentityOrTranslation())
1454 return true;
1455
1456 double det = WebCore::determinant4x4(m_matrix);
1457
1458 if (fabs(det) < SMALL_NUMBER)
1459 return false;
1460
1461 return true;
1462}
1463
1464Optional<TransformationMatrix> TransformationMatrix::inverse() const
1465{
1466 if (isIdentityOrTranslation()) {
1467 // identity matrix
1468 if (m_matrix[3][0] == 0 && m_matrix[3][1] == 0 && m_matrix[3][2] == 0)
1469 return TransformationMatrix();
1470
1471 // translation
1472 return TransformationMatrix(1, 0, 0, 0,
1473 0, 1, 0, 0,
1474 0, 0, 1, 0,
1475 -m_matrix[3][0], -m_matrix[3][1], -m_matrix[3][2], 1);
1476 }
1477
1478 TransformationMatrix invMat;
1479 // FIXME: Use LU decomposition to apply the inverse instead of calculating the inverse explicitly.
1480 // Calculating the inverse of a 4x4 matrix using cofactors is numerically unstable and unnecessary to apply the inverse transformation to a point.
1481 if (!WebCore::inverse(m_matrix, invMat.m_matrix))
1482 return WTF::nullopt;
1483
1484 return invMat;
1485}
1486
1487void TransformationMatrix::makeAffine()
1488{
1489 m_matrix[0][2] = 0;
1490 m_matrix[0][3] = 0;
1491
1492 m_matrix[1][2] = 0;
1493 m_matrix[1][3] = 0;
1494
1495 m_matrix[2][0] = 0;
1496 m_matrix[2][1] = 0;
1497 m_matrix[2][2] = 1;
1498 m_matrix[2][3] = 0;
1499
1500 m_matrix[3][2] = 0;
1501 m_matrix[3][3] = 1;
1502}
1503
1504AffineTransform TransformationMatrix::toAffineTransform() const
1505{
1506 return AffineTransform(m_matrix[0][0], m_matrix[0][1], m_matrix[1][0],
1507 m_matrix[1][1], m_matrix[3][0], m_matrix[3][1]);
1508}
1509
1510static inline void blendFloat(double& from, double to, double progress)
1511{
1512 if (from != to)
1513 from = from + (to - from) * progress;
1514}
1515
1516void TransformationMatrix::blend2(const TransformationMatrix& from, double progress)
1517{
1518 Decomposed2Type fromDecomp;
1519 Decomposed2Type toDecomp;
1520 if (!from.decompose2(fromDecomp) || !decompose2(toDecomp)) {
1521 if (progress < 0.5)
1522 *this = from;
1523 return;
1524 }
1525
1526 // If x-axis of one is flipped, and y-axis of the other, convert to an unflipped rotation.
1527 if ((fromDecomp.scaleX < 0 && toDecomp.scaleY < 0) || (fromDecomp.scaleY < 0 && toDecomp.scaleX < 0)) {
1528 fromDecomp.scaleX = -fromDecomp.scaleX;
1529 fromDecomp.scaleY = -fromDecomp.scaleY;
1530 fromDecomp.angle += fromDecomp.angle < 0 ? 180 : -180;
1531 }
1532
1533 // Don't rotate the long way around.
1534 if (!fromDecomp.angle)
1535 fromDecomp.angle = 360;
1536 if (!toDecomp.angle)
1537 toDecomp.angle = 360;
1538
1539 if (fabs(fromDecomp.angle - toDecomp.angle) > 180) {
1540 if (fromDecomp.angle > toDecomp.angle)
1541 fromDecomp.angle -= 360;
1542 else
1543 toDecomp.angle -= 360;
1544 }
1545
1546 blendFloat(fromDecomp.m11, toDecomp.m11, progress);
1547 blendFloat(fromDecomp.m12, toDecomp.m12, progress);
1548 blendFloat(fromDecomp.m21, toDecomp.m21, progress);
1549 blendFloat(fromDecomp.m22, toDecomp.m22, progress);
1550 blendFloat(fromDecomp.translateX, toDecomp.translateX, progress);
1551 blendFloat(fromDecomp.translateY, toDecomp.translateY, progress);
1552 blendFloat(fromDecomp.scaleX, toDecomp.scaleX, progress);
1553 blendFloat(fromDecomp.scaleY, toDecomp.scaleY, progress);
1554 blendFloat(fromDecomp.angle, toDecomp.angle, progress);
1555
1556 recompose2(fromDecomp);
1557}
1558
1559void TransformationMatrix::blend4(const TransformationMatrix& from, double progress)
1560{
1561 Decomposed4Type fromDecomp;
1562 Decomposed4Type toDecomp;
1563 if (!from.decompose4(fromDecomp) || !decompose4(toDecomp)) {
1564 if (progress < 0.5)
1565 *this = from;
1566 return;
1567 }
1568
1569 blendFloat(fromDecomp.scaleX, toDecomp.scaleX, progress);
1570 blendFloat(fromDecomp.scaleY, toDecomp.scaleY, progress);
1571 blendFloat(fromDecomp.scaleZ, toDecomp.scaleZ, progress);
1572 blendFloat(fromDecomp.skewXY, toDecomp.skewXY, progress);
1573 blendFloat(fromDecomp.skewXZ, toDecomp.skewXZ, progress);
1574 blendFloat(fromDecomp.skewYZ, toDecomp.skewYZ, progress);
1575 blendFloat(fromDecomp.translateX, toDecomp.translateX, progress);
1576 blendFloat(fromDecomp.translateY, toDecomp.translateY, progress);
1577 blendFloat(fromDecomp.translateZ, toDecomp.translateZ, progress);
1578 blendFloat(fromDecomp.perspectiveX, toDecomp.perspectiveX, progress);
1579 blendFloat(fromDecomp.perspectiveY, toDecomp.perspectiveY, progress);
1580 blendFloat(fromDecomp.perspectiveZ, toDecomp.perspectiveZ, progress);
1581 blendFloat(fromDecomp.perspectiveW, toDecomp.perspectiveW, progress);
1582
1583 slerp(&fromDecomp.quaternionX, &toDecomp.quaternionX, progress);
1584
1585 recompose4(fromDecomp);
1586}
1587
1588void TransformationMatrix::blend(const TransformationMatrix& from, double progress)
1589{
1590 if (from.isIdentity() && isIdentity())
1591 return;
1592
1593 if (from.isAffine() && isAffine())
1594 blend2(from, progress);
1595 else
1596 blend4(from, progress);
1597}
1598
1599bool TransformationMatrix::decompose2(Decomposed2Type& decomp) const
1600{
1601 if (isIdentity()) {
1602 memset(&decomp, 0, sizeof(decomp));
1603 decomp.scaleX = 1;
1604 decomp.scaleY = 1;
1605 decomp.m11 = 1;
1606 decomp.m22 = 1;
1607 return true;
1608 }
1609
1610 return WebCore::decompose2(m_matrix, decomp);
1611}
1612
1613bool TransformationMatrix::decompose4(Decomposed4Type& decomp) const
1614{
1615 if (isIdentity()) {
1616 memset(&decomp, 0, sizeof(decomp));
1617 decomp.perspectiveW = 1;
1618 decomp.scaleX = 1;
1619 decomp.scaleY = 1;
1620 decomp.scaleZ = 1;
1621 return true;
1622 }
1623
1624 return WebCore::decompose4(m_matrix, decomp);
1625}
1626
1627void TransformationMatrix::recompose2(const Decomposed2Type& decomp)
1628{
1629 makeIdentity();
1630
1631 m_matrix[0][0] = decomp.m11;
1632 m_matrix[0][1] = decomp.m12;
1633 m_matrix[1][0] = decomp.m21;
1634 m_matrix[1][1] = decomp.m22;
1635
1636 translate3d(decomp.translateX, decomp.translateY, 0);
1637 rotate(decomp.angle);
1638 scale3d(decomp.scaleX, decomp.scaleY, 1);
1639}
1640
1641void TransformationMatrix::recompose4(const Decomposed4Type& decomp)
1642{
1643 makeIdentity();
1644
1645 // First apply perspective.
1646 m_matrix[0][3] = decomp.perspectiveX;
1647 m_matrix[1][3] = decomp.perspectiveY;
1648 m_matrix[2][3] = decomp.perspectiveZ;
1649 m_matrix[3][3] = decomp.perspectiveW;
1650
1651 // Next, translate.
1652 translate3d(decomp.translateX, decomp.translateY, decomp.translateZ);
1653
1654 // Apply rotation.
1655 double xx = decomp.quaternionX * decomp.quaternionX;
1656 double xy = decomp.quaternionX * decomp.quaternionY;
1657 double xz = decomp.quaternionX * decomp.quaternionZ;
1658 double xw = decomp.quaternionX * decomp.quaternionW;
1659 double yy = decomp.quaternionY * decomp.quaternionY;
1660 double yz = decomp.quaternionY * decomp.quaternionZ;
1661 double yw = decomp.quaternionY * decomp.quaternionW;
1662 double zz = decomp.quaternionZ * decomp.quaternionZ;
1663 double zw = decomp.quaternionZ * decomp.quaternionW;
1664
1665 // Construct a composite rotation matrix from the quaternion values.
1666 TransformationMatrix rotationMatrix(1 - 2 * (yy + zz), 2 * (xy - zw), 2 * (xz + yw), 0,
1667 2 * (xy + zw), 1 - 2 * (xx + zz), 2 * (yz - xw), 0,
1668 2 * (xz - yw), 2 * (yz + xw), 1 - 2 * (xx + yy), 0,
1669 0, 0, 0, 1);
1670
1671 multiply(rotationMatrix);
1672
1673 // Apply skew.
1674 if (decomp.skewYZ) {
1675 TransformationMatrix tmp;
1676 tmp.setM32(decomp.skewYZ);
1677 multiply(tmp);
1678 }
1679
1680 if (decomp.skewXZ) {
1681 TransformationMatrix tmp;
1682 tmp.setM31(decomp.skewXZ);
1683 multiply(tmp);
1684 }
1685
1686 if (decomp.skewXY) {
1687 TransformationMatrix tmp;
1688 tmp.setM21(decomp.skewXY);
1689 multiply(tmp);
1690 }
1691
1692 // Finally, apply scale.
1693 scale3d(decomp.scaleX, decomp.scaleY, decomp.scaleZ);
1694}
1695
1696bool TransformationMatrix::isIntegerTranslation() const
1697{
1698 if (!isIdentityOrTranslation())
1699 return false;
1700
1701 // Check for translate Z.
1702 if (m_matrix[3][2])
1703 return false;
1704
1705 // Check for non-integer translate X/Y.
1706 if (static_cast<int>(m_matrix[3][0]) != m_matrix[3][0] || static_cast<int>(m_matrix[3][1]) != m_matrix[3][1])
1707 return false;
1708
1709 return true;
1710}
1711
1712bool TransformationMatrix::containsOnlyFiniteValues() const
1713{
1714 return std::isfinite(m_matrix[0][0]) && std::isfinite(m_matrix[0][1]) && std::isfinite(m_matrix[0][2]) && std::isfinite(m_matrix[0][3])
1715 && std::isfinite(m_matrix[1][0]) && std::isfinite(m_matrix[1][1]) && std::isfinite(m_matrix[1][2]) && std::isfinite(m_matrix[1][3])
1716 && std::isfinite(m_matrix[2][0]) && std::isfinite(m_matrix[2][1]) && std::isfinite(m_matrix[2][2]) && std::isfinite(m_matrix[2][3])
1717 && std::isfinite(m_matrix[3][0]) && std::isfinite(m_matrix[3][1]) && std::isfinite(m_matrix[3][2]) && std::isfinite(m_matrix[3][3]);
1718}
1719
1720TransformationMatrix TransformationMatrix::to2dTransform() const
1721{
1722 return TransformationMatrix(m_matrix[0][0], m_matrix[0][1], 0, m_matrix[0][3],
1723 m_matrix[1][0], m_matrix[1][1], 0, m_matrix[1][3],
1724 0, 0, 1, 0,
1725 m_matrix[3][0], m_matrix[3][1], 0, m_matrix[3][3]);
1726}
1727
1728auto TransformationMatrix::toColumnMajorFloatArray() const -> FloatMatrix4
1729{
1730 return { {
1731 float(m11()), float(m12()), float(m13()), float(m14()),
1732 float(m21()), float(m22()), float(m23()), float(m24()),
1733 float(m31()), float(m32()), float(m33()), float(m34()),
1734 float(m41()), float(m42()), float(m43()), float(m44()) } };
1735}
1736
1737bool TransformationMatrix::isBackFaceVisible() const
1738{
1739 // Back-face visibility is determined by transforming the normal vector (0, 0, 1) and
1740 // checking the sign of the resulting z component. However, normals cannot be
1741 // transformed by the original matrix, they require being transformed by the
1742 // inverse-transpose.
1743 //
1744 // Since we know we will be using (0, 0, 1), and we only care about the z-component of
1745 // the transformed normal, then we only need the m33() element of the
1746 // inverse-transpose. Therefore we do not need the transpose.
1747 //
1748 // Additionally, if we only need the m33() element, we do not need to compute a full
1749 // inverse. Instead, knowing the inverse of a matrix is adjoint(matrix) / determinant,
1750 // we can simply compute the m33() of the adjoint (adjugate) matrix, without computing
1751 // the full adjoint.
1752
1753 double determinant = WebCore::determinant4x4(m_matrix);
1754
1755 // If the matrix is not invertible, then we assume its backface is not visible.
1756 if (fabs(determinant) < SMALL_NUMBER)
1757 return false;
1758
1759 double cofactor33 = determinant3x3(m11(), m12(), m14(), m21(), m22(), m24(), m41(), m42(), m44());
1760 double zComponentOfTransformedNormal = cofactor33 / determinant;
1761
1762 return zComponentOfTransformedNormal < 0;
1763}
1764
1765TextStream& operator<<(TextStream& ts, const TransformationMatrix& transform)
1766{
1767 TextStream::IndentScope indentScope(ts);
1768 ts << "\n";
1769 ts << indent << "[" << transform.m11() << " " << transform.m12() << " " << transform.m13() << " " << transform.m14() << "]\n";
1770 ts << indent << "[" << transform.m21() << " " << transform.m22() << " " << transform.m23() << " " << transform.m24() << "]\n";
1771 ts << indent << "[" << transform.m31() << " " << transform.m32() << " " << transform.m33() << " " << transform.m34() << "]\n";
1772 ts << indent << "[" << transform.m41() << " " << transform.m42() << " " << transform.m43() << " " << transform.m44() << "]";
1773 return ts;
1774}
1775
1776}
1777