Skip to content

Commit

Permalink
TransformationMatrix::Recompose() and Decompose() incorrectly transpo…
Browse files Browse the repository at this point in the history 
…se rotation.

https://bugs.webkit.org/show_bug.cgi?id=220856
<rdar://problem/73747851>

Reviewed by Dean Jackson.

The computation of quaternions during matrix decomposition appeared
to be based off of the matrix transpose. Indexing in the quaternion
calculation is consistent with row-major ordering; however, the "rows"
variable is actually column major order (goes back to the original
implementation of unmatrix in Graphics Gems II).  The rotation matrix
constructed from the quaternions also appeared to be assuming row-major
ordering. These discrepancies were mostly corrected for by flipping the
direction of the quaternion when extracting from the matrix.

This is a (partial) cherry-pick from Blink:
  - https://chromium-review.googlesource.com/c/chromium/src/+/2489727

* Source/WebCore/platform/graphics/transforms/RotateTransformOperation.cpp:
(WebCore::RotateTransformOperation::blend):
* Source/WebCore/platform/graphics/transforms/TransformationMatrix.cpp:
(WebCore::decompose4):
(WebCore::slerp):
(WebCore::accumulateQuaternion):
(WebCore::TransformationMatrix::recompose4):

Canonical link: https://commits.webkit.org/267424@main
  • Loading branch information
@mattwoodrow
mattwoodrow committed Aug 29, 2023
1 parent 616b80d commit a959acf
Show file tree
Hide file tree
Showing 2 changed files with 124 additions and 117 deletions.
6 changes: 3 additions & 3 deletions Source/WebCore/platform/graphics/transforms/RotateTransformOperation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,9 +114,9 @@ Ref<TransformOperation> RotateTransformOperation::blend(const TransformOperation
}

// Convert that to Axis/Angle form
double x = -decomp.quaternionX;
double y = -decomp.quaternionY;
double z = -decomp.quaternionZ;
double x = decomp.quaternionX;
double y = decomp.quaternionY;
double z = decomp.quaternionZ;
double length = std::hypot(x, y, z);
double angle = 0;

Expand Down
235 changes: 121 additions & 114 deletions Source/WebCore/platform/graphics/transforms/TransformationMatrix.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -419,103 +419,131 @@ static bool decompose4(const TransformationMatrix::Matrix4& mat, TransformationM
result.translateZ = localMatrix[3][2];
localMatrix[3][2] = 0;

// Vector4 type and functions need to be added to the common set.
Vector3 row[3], pdum3;
// Note: Deviating from the spec in terms of variable naming. The matrix is
// stored on column major order and not row major. Using the variable 'row'
// instead of 'column' in the spec pseudocode has been the source of
// confusion, specifically in sorting out rotations.
Vector3 column[3], pdum3;

// Now get scale and shear.
for (i = 0; i < 3; i++) {
row[i][0] = localMatrix[i][0];
row[i][1] = localMatrix[i][1];
row[i][2] = localMatrix[i][2];
column[i][0] = localMatrix[i][0];
column[i][1] = localMatrix[i][1];
column[i][2] = localMatrix[i][2];
}

// Compute X scale factor and normalize first row.
result.scaleX = v3Length(row[0]);
v3Scale(row[0], 1.0);
// Compute X scale factor and normalize the first column.
result.scaleX = v3Length(column[0]);
v3Scale(column[0], 1.0);

// Compute XY shear factor and make 2nd row orthogonal to 1st.
result.skewXY = v3Dot(row[0], row[1]);
v3Combine(row[1], row[0], row[1], 1.0, -result.skewXY);
// Compute XY shear factor and make 2nd column orthogonal to 1st.
result.skewXY = v3Dot(column[0], column[1]);
v3Combine(column[1], column[0], column[1], 1.0, -result.skewXY);

// Now, compute Y scale and normalize 2nd row.
result.scaleY = v3Length(row[1]);
v3Scale(row[1], 1.0);
// Now, compute Y scale and normalize 2nd column.
result.scaleY = v3Length(column[1]);
v3Scale(column[1], 1.0);
result.skewXY /= result.scaleY;

// Compute XZ and YZ shears, orthogonalize 3rd row.
result.skewXZ = v3Dot(row[0], row[2]);
v3Combine(row[2], row[0], row[2], 1.0, -result.skewXZ);
result.skewYZ = v3Dot(row[1], row[2]);
v3Combine(row[2], row[1], row[2], 1.0, -result.skewYZ);
// Compute XZ and YZ shears, orthogonalize 3rd column.
result.skewXZ = v3Dot(column[0], column[2]);
v3Combine(column[2], column[0], column[2], 1.0, -result.skewXZ);
result.skewYZ = v3Dot(column[1], column[2]);
v3Combine(column[2], column[1], column[2], 1.0, -result.skewYZ);

// Next, get Z scale and normalize 3rd row.
result.scaleZ = v3Length(row[2]);
v3Scale(row[2], 1.0);
// Next, get Z scale and normalize 3rd column.
result.scaleZ = v3Length(column[2]);
v3Scale(column[2], 1.0);
result.skewXZ /= result.scaleZ;
result.skewYZ /= result.scaleZ;

// At this point, the matrix (in rows[]) is orthonormal.
// At this point, the matrix (in column[]) is orthonormal.
// Check for a coordinate system flip. If the determinant
// is -1, then negate the matrix and the scaling factors.
v3Cross(row[1], row[2], pdum3);
if (v3Dot(row[0], pdum3) < 0) {
v3Cross(column[1], column[2], pdum3);
if (v3Dot(column[0], pdum3) < 0) {

result.scaleX *= -1;
result.scaleY *= -1;
result.scaleZ *= -1;

for (i = 0; i < 3; i++) {
row[i][0] *= -1;
row[i][1] *= -1;
row[i][2] *= -1;
column[i][0] *= -1;
column[i][1] *= -1;
column[i][2] *= -1;
}
}

// Now, get the rotations out, as described in the gem.

// FIXME - Add the ability to return either quaternions (which are
// easier to recompose with) or Euler angles (rx, ry, rz), which
// are easier for authors to deal with. The latter will only be useful
// when we fix https://bugs.webkit.org/show_bug.cgi?id=23799, so I
// will leave the Euler angle code here for now.

// ret.rotateY = asin(-row[0][2]);
// if (cos(ret.rotateY) != 0) {
// ret.rotateX = atan2(row[1][2], row[2][2]);
// ret.rotateZ = atan2(row[0][1], row[0][0]);
// } else {
// ret.rotateX = atan2(-row[2][0], row[1][1]);
// ret.rotateZ = 0;
// }

double s, t, x, y, z, w;

t = row[0][0] + row[1][1] + row[2][2] + 1.0;

if (t > 1e-4) {
s = 0.5 / sqrt(t);
w = 0.25 / s;
x = (row[2][1] - row[1][2]) * s;
y = (row[0][2] - row[2][0]) * s;
z = (row[1][0] - row[0][1]) * s;
} else if (row[0][0] > row[1][1] && row[0][0] > row[2][2]) {
s = sqrt(1.0 + row[0][0] - row[1][1] - row[2][2]) * 2.0; // S = 4 * qx.
x = 0.25 * s;
y = (row[0][1] + row[1][0]) / s;
z = (row[0][2] + row[2][0]) / s;
w = (row[2][1] - row[1][2]) / s;
} else if (row[1][1] > row[2][2]) {
s = sqrt(1.0 + row[1][1] - row[0][0] - row[2][2]) * 2.0; // S = 4 * qy.
x = (row[0][1] + row[1][0]) / s;
y = 0.25 * s;
z = (row[1][2] + row[2][1]) / s;
w = (row[0][2] - row[2][0]) / s;
// Lastly, compute the quaternions.
// See https://en.wikipedia.org/wiki/Rotation_matrix#Quaternion.
// Note: deviating from spec (http://www.w3.org/TR/css3-transforms/)
// which has a degenerate case when the trace (t) of the orthonormal matrix
// (Q) approaches -1. In the Wikipedia article, Q_ij is indexing on row then
// column. Thus, Q_ij = column[j][i].

// The following are equivalent represnetations of the rotation matrix:
//
// Axis-angle form:
//
// [ c+(1-c)x^2 (1-c)xy-sz (1-c)xz+sy ] c = cos theta
// R = [ (1-c)xy+sz c+(1-c)y^2 (1-c)yz-sx ] s = sin theta
// [ (1-c)xz-sy (1-c)yz+sx c+(1-c)z^2 ] [x,y,z] = axis or rotation
//
// The sum of the diagonal elements (trace) is a simple function of the cosine
// of the angle. The w component of the quaternion is cos(theta/2), and we
// make use of the double angle formula to directly compute w from the
// trace. Differences between pairs of skew symmetric elements in this matrix
// isolate the remaining components. Since w can be zero (also numerically
// unstable if near zero), we cannot rely solely on this approach to compute
// the quaternion components.
//
// Quaternion form:
//
// [ 1-2(y^2+z^2) 2(xy-zw) 2(xz+yw) ]
// r = [ 2(xy+zw) 1-2(x^2+z^2) 2(yz-xw) ] q = (x,y,y,w)
// [ 2(xz-yw) 2(yz+xw) 1-2(x^2+y^2) ]
//
// Different linear combinations of the diagonal elements isolates x, y or z.
// Sums or differences between skew symmetric elements isolate the remainder.

double r, s, t, x, y, z, w;

t = column[0][0] + column[1][1] + column[2][2];

// https://en.wikipedia.org/wiki/Rotation_matrix#Quaternion
if (1 + t > 0.001) {
// Numerically stable as long as 1+t is not close to zero. Otherwise use the
// diagonal element with the greatest value to compute the quaternions.
r = std::sqrt(1.0 + t);
s = 0.5 / r;
w = 0.5 * r;
x = (column[1][2] - column[2][1]) * s;
y = (column[2][0] - column[0][2]) * s;
z = (column[0][1] - column[1][0]) * s;
} else if (column[0][0] > column[1][1] && column[0][0] > column[2][2]) {
// Q_xx is largest.
r = std::sqrt(1.0 + column[0][0] - column[1][1] - column[2][2]);
s = 0.5 / r;
x = 0.5 * r;
y = (column[1][0] - column[0][1]) * s;
z = (column[2][0] + column[0][2]) * s;
w = (column[1][2] - column[2][1]) * s;
} else if (column[1][1] > column[2][2]) {
// Q_yy is largest.
r = std::sqrt(1.0 - column[0][0] + column[1][1] - column[2][2]);
s = 0.5 / r;
x = (column[1][0] + column[0][1]) * s;
y = 0.5 * r;
z = (column[2][1] + column[1][2]) * s;
w = (column[2][0] - column[0][2]) * s;
} else {
s = sqrt(1.0 + row[2][2] - row[0][0] - row[1][1]) * 2.0; // S = 4 * qz.
x = (row[0][2] + row[2][0]) / s;
y = (row[1][2] + row[2][1]) / s;
z = 0.25 * s;
w = (row[1][0] - row[0][1]) / s;
// Q_zz is largest.
r = std::sqrt(1.0 - column[0][0] - column[1][1] + column[2][2]);
s = 0.5 / r;
x = (column[2][0] + column[0][2]) * s;
y = (column[2][1] + column[1][2]) * s;
z = 0.5 * r;
w = (column[0][1] - column[1][0]) * s;
}

result.quaternionX = x;
Expand All @@ -528,44 +556,38 @@ static bool decompose4(const TransformationMatrix::Matrix4& mat, TransformationM

// Perform a spherical linear interpolation between the two
// passed quaternions with 0 <= t <= 1.

static void slerp(double qa[4], const double qb[4], double t)
{
const double kEpsilon = 1e-5;
double ax, ay, az, aw;
double bx, by, bz, bw;
double cx, cy, cz, cw;
double angle;
double th, invth, scale, invscale;

ax = qa[0]; ay = qa[1]; az = qa[2]; aw = qa[3];
bx = qb[0]; by = qb[1]; bz = qb[2]; bw = qb[3];

angle = ax * bx + ay * by + az * bz + aw * bw;
double cosHalfAngle = ax * bx + ay * by + az * bz + aw * bw;

if (angle < 0.0) {
if (cosHalfAngle < 0.0) {
ax = -ax; ay = -ay;
az = -az; aw = -aw;
angle = -angle;
cosHalfAngle = -cosHalfAngle;
}

if (angle + 1.0 > .05) {
if (1.0 - angle >= .05) {
th = acos (angle);
invth = 1.0 / sin (th);
scale = sin (th * (1.0 - t)) * invth;
invscale = sin (th * t) * invth;
} else {
scale = 1.0 - t;
invscale = t;
}
} else {
bx = -ay;
by = ax;
bz = -aw;
bw = az;
scale = sin(piDouble * (.5 - t));
invscale = sin (piDouble * t);
if (cosHalfAngle > 1)
cosHalfAngle = 1;

double sinHalfAngle = std::sqrt(1.0 - cosHalfAngle * cosHalfAngle);
if (sinHalfAngle < kEpsilon) {
// Quaternions share common axis and angle.
return;
}

double halfAngle = std::acos(cosHalfAngle);
double scale = std::sin((1 - t) * halfAngle) / sinHalfAngle;
double invscale = std::sin(t * halfAngle) / sinHalfAngle;

cx = ax * scale + bx * invscale;
cy = ay * scale + by * invscale;
cz = az * scale + bz * invscale;
Expand All @@ -581,7 +603,6 @@ TransformationMatrix::TransformationMatrix(const AffineTransform& t)
setMatrix(t.a(), t.b(), t.c(), t.d(), t.e(), t.f());
}

// FIXME: Once https://bugs.webkit.org/show_bug.cgi?id=220856 is addressed we can reuse this function in TransformationMatrix::recompose4().
TransformationMatrix TransformationMatrix::fromQuaternion(double qx, double qy, double qz, double qw)
{
double xx = qx * qx;
Expand Down Expand Up @@ -1730,10 +1751,10 @@ void TransformationMatrix::blend2(const TransformationMatrix& from, double progr
// Compute quaternion multiplication
static void accumulateQuaternion(double qa[4], const double qb[4])
{
auto qx = (qb[3] * qa[0]) + (qb[0] * qa[3]) + (qb[1] * qa[2]) - (qb[2] * qa[1]);
auto qy = (qb[3] * qa[1]) + (qb[1] * qa[3]) + (qb[2] * qa[0]) - (qb[0] * qa[2]);
auto qz = (qb[3] * qa[2]) + (qb[2] * qa[3]) + (qb[0] * qa[1]) - (qb[1] * qa[0]);
auto qw = (qb[3] * qa[3]) - (qb[0] * qa[0]) - (qb[1] * qa[1]) - (qb[2] * qa[2]);
auto qx = qa[3] * qb[0] + qa[0] * qb[3] + qa[1] * qb[2] - qa[2] * qb[1];
auto qy = qa[3] * qb[1] - qa[0] * qb[2] + qa[1] * qb[3] + qa[2] * qb[0];
auto qz = qa[3] * qb[2] + qa[0] * qb[1] - qa[1] * qb[0] + qa[2] * qb[3];
auto qw = qa[3] * qb[3] - qa[0] * qb[0] - qa[1] * qb[1] - qa[2] * qb[2];
qa[0] = qx; qa[1] = qy; qa[2] = qz; qa[3] = qw;
}

Expand Down Expand Up @@ -1859,21 +1880,7 @@ void TransformationMatrix::recompose4(const Decomposed4Type& decomp)
translate3d(decomp.translateX, decomp.translateY, decomp.translateZ);

// Apply rotation.
double xx = decomp.quaternionX * decomp.quaternionX;
double xy = decomp.quaternionX * decomp.quaternionY;
double xz = decomp.quaternionX * decomp.quaternionZ;
double xw = decomp.quaternionX * decomp.quaternionW;
double yy = decomp.quaternionY * decomp.quaternionY;
double yz = decomp.quaternionY * decomp.quaternionZ;
double yw = decomp.quaternionY * decomp.quaternionW;
double zz = decomp.quaternionZ * decomp.quaternionZ;
double zw = decomp.quaternionZ * decomp.quaternionW;

// Construct a composite rotation matrix from the quaternion values.
TransformationMatrix rotationMatrix(1 - 2 * (yy + zz), 2 * (xy - zw), 2 * (xz + yw), 0,
2 * (xy + zw), 1 - 2 * (xx + zz), 2 * (yz - xw), 0,
2 * (xz - yw), 2 * (yz + xw), 1 - 2 * (xx + yy), 0,
0, 0, 0, 1);
TransformationMatrix rotationMatrix = TransformationMatrix::fromQuaternion(decomp.quaternionX, decomp.quaternionY, decomp.quaternionZ, decomp.quaternionW);

multiply(rotationMatrix);

Expand Down

0 comments on commit a959acf

Please sign in to comment.