[css-transforms] Quaternion calculation for 3D matrix decomposition #3710
Description
Decomposition of a 3D matrix:
https://drafts.csswg.org/css-transforms-2/#decomposing-a-3d-matrix
The pseudocode for quaternion calculations contains the following:
// Now, get the rotations out
quaternion[0] = 0.5 * sqrt(max(1 + row[0][0] - row[1][1] - row[2][2], 0))
quaternion[1] = 0.5 * sqrt(max(1 - row[0][0] + row[1][1] - row[2][2], 0))
quaternion[2] = 0.5 * sqrt(max(1 - row[0][0] - row[1][1] + row[2][2], 0))
quaternion[3] = 0.5 * sqrt(max(1 + row[0][0] + row[1][1] + row[2][2], 0))
if (row[2][1] > row[1][2])
quaternion[0] = -quaternion[0]
if (row[0][2] > row[2][0])
quaternion[1] = -quaternion[1]
if (row[1][0] > row[0][1])
quaternion[2] = -quaternion[2]
This implementation hits a degenerate case when off diagonal elements in the orthonormal matrix used for determining the sign are equal.
Consider the decomposition of matrix(1, 1, 1, 0, 0, 0)
Following the decomposition steps, the orthonormal matrix for computing the quaternion is:
[-0.7071 -0.7071 0]
[-0.7071 0.7071 0]
[ 0 0 -1]
The off-diagonal elements are equal and we cannot determine, which elements in the quaternion should be negative.
Based on the algorithm in the spec, (x, y, z, w) = (0.382683, 0.92388, 0, 0) when it should be (x, y, z, w) = (-0.382683, 0.92388, 0, 0)
A more robust algorithm is:
t = q_xx + q_yy + q_zz
if (t > 0)
r = sqrt(1.0 + t)
s = 0.5 / r
w = 0.5 * r
x = (q_zy - q_yz) * s
y = (q_xz - q_zx) * s
z = (q_yx - q_xy) * s
else if (q_xx > q_yy && q_xx > q_zz)
r = sqrt(1.0 + q_xx - q_yy - q_zz)
s = 0.5 / r
x = 0.5 * r
y = (q_xy + q_yx) * s
z = (q_xz + q_zx) * s
w = (q_zy - q_yz) * s
else if (q_yy > q_zz)
r = sqrt(1.0 - q_xx + q_yy - q_zz)
s = 0.5 / r
x = (q_xy + q_yx) * s
y = 0.5 * r
z = (q_yz + q_zy) * s
w = (q_xz - q_zx) * s
else
r = sqrt(1.0 - q_xx - q_yy + q_zz)
s = 0.5 / r
x = (q_xz + q_zx) * s
y = (q_yz + q_zy) * s
z = 0.5 * r
w = (q_yx - q_xy) * s