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