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