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