Re: [css3-transitions] Interpolating transforms: a proposal to avoid Euler angles in favor of using quaternions

On 21/10/2010, at 5:19 AM, Boris Zbarsky wrote:

> This is a message that was sent to me with the intent of it being forwarded on to this list.  It describes some problems with the existing setup for interpolating transforms and suggests alternatives.
> 
> -------- Original Message --------
> 
> http://dev.w3.org/csswg/css3-2d-transforms/#animation describes the
> current method.
> 
> As I see it, there are two main problems with the current approach: the
> interpolation of rotations and the handling of matrices with negative
> determinant.
> 
> Rotations are currently decomposed into a series of Euler angles, and
> the angle values interpolated directly. Such a situation is notoriously
> susceptible to Gimbal lock: http://en.wikipedia.org/wiki/Gimbal_lock
> This results in numeric instability near the poles: the value of one of
> the angles derived from a rotation matrix becomes almost arbitrary,
> depending on numerical rounding error, which can result in macroscopic
> deviations in the animation. A much better approach is to use
> quaternions: http://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation
> 
> Simple ways to interpolate rotations with quaternions have long been
> known: http://portal.acm.org/citation.cfm?doid=325334.325242
> They completely avoid Gimbal lock, and provide a number of other
> benefits, like "binvariance": given three rotations Q, R, and S, and a
> path between Q and R, denoted P(t), then the path between S.Q and S.R is
> given by S.P(t) (left invariance), and the path between R.S and Q.S is
> given by P(t).S (right invariance). Direct interpolation via Euler
> angles does not have this property.

It turns out this is what we implemented in WebKit, so we support making such a change. Here are the details of the WebKit implementation:

Run unmatrix as specified, but don't use the rotate[] values. Instead, calculate four quaternions using

    double quaternionX, quaternionY, quaternionZ, quaternionW;
    double s, t;

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

    if (t > 1e-4) {
        s = 0.5 / sqrt(t);
        quaternionX = (row[2][1] - row[1][2]) * s;
        quaternionY = (row[0][2] - row[2][0]) * s;
        quaternionZ = (row[1][0] - row[0][1]) * s;
        quaternionW = 0.25 / 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 
        quaternionX = 0.25 * s;
        quaternionY = (row[0][1] + row[1][0]) / s; 
        quaternionZ = (row[0][2] + row[2][0]) / s; 
        quaternionW = (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
        quaternionX = (row[0][1] + row[1][0]) / s; 
        quaternionY = 0.25 * s;
        quaternionZ = (row[1][2] + row[2][1]) / s; 
        quaternionW = (row[0][2] - row[2][0]) / s;
    } else { 
        s = sqrt(1.0 + row[2][2] - row[0][0] - row[1][1]) * 2.0; // S=4*qz
        quaternionX = (row[0][2] + row[2][0]) / s;
        quaternionY = (row[1][2] + row[2][1]) / s; 
        quaternionZ = 0.25 * s;
        quaternionW = (row[1][0] - row[0][1]) / s;
    }

Then, when we interpolate between matrices, do the scale, skew, translate and perspective parts as specified, but use slerp (http://en.wikipedia.org/wiki/Slerp) for the interpolation of the quaternions. Then recompose a resulting matrix using the interpolated values.

This still leaves the shorthand properties susceptible to Gimbal lock, but I think that's fine. The author can always use a matrix if that's an issue for them.

Now, I'm going to wimp out here and ask that Boris, the original author, or someone else check this to make sure it is acceptable. If so, I'll add it to the specification.

There is one more thing below, the suggestion to  matrix[2][2]=sign(det(matrix2d)) when converting from 2d to 3d. I propose we accept that as well.

Dean



> 
> The best approach I've found for converting a 3x3 rotation matrix into a  quaternion is:
> 
> #define MAX(a,b) ((a)>(b)?(a):(b))
> 
>  w=0.5*sqrt(MAX(1+row[0][0]+row[1][1]+row[2][2],0));
>  x=0.5*sqrt(MAX(1+row[0][0]-row[1][1]-row[2][2],0));
>  y=0.5*sqrt(MAX(1-row[0][0]+row[1][1]-row[2][2],0));
>  z=0.5*sqrt(MAX(1-row[0][0]-row[1][1]+row[2][2],0));
>  if(row[2][1]<row[1][2])x=-x;
>  if(row[0][2]<row[2][0])y=-y;
>  if(row[1][0]<row[0][1])z=-z;
> 
> The result, {w,x,y,z}, is a unit quaternion. You may need to flip the
> direction of the comparisons in the last three rows if you're applying
> your rotation matrices from the right instead of the left. This is the
> only approach I've seen that requires no divisions at all. It also
> requires no epsilons or other fudge factors that would complicate the
> definition of a standard. It is numerically stable for rotations near 0, 90, and 180 degrees around any axis, unlike many other approaches. The definition of the macro MAX here is very specific. The second argument must be the default when the two compare equal. Otherwise the arithmetic  in the first argument can produce -0 as a result, causing sqrt to return NaN even with the check in place.
> 
> This entire discussion only applies to 3D rotations: for 2D rotations
> quaternions are equivalent to Euler angles, because all rotations can be represented by a single angle describing the rotation around the Z axis.
> 
> However, it becomes relevant when coupled with the next issue.
> 
> David Baron outlined some of the problems with matrices with a negative
> determinant in
> http://lists.w3.org/Archives/Public/www-style/2010Jun/0602.html
> He comes up with two possibilities for decomposing them into a series of 2D transforms:
> 
> a) Invert scaleY and XYshear (i.e., multiply by
> {{1,0,0},{0,-1,0},{0,0,1}} on the right).
> 
> b) Invert scaleX, XYshear (i.e., multiply by {{-1,0},{0,1}} on the
> left), and A, B, C, and D (so that the rotation computed from A and B
> ends up in the opposite quadrant).
> 
> There is yet a third option that I will argue makes much more intuitive
> sense, which is this: when converting from a 2D matrix to a 3D matrix,
> set matrix[2][2]=sign(det(matrix2d)), i.e., force the determinant of the 3D matrix to be positive.
> 
> The reasoning behind this is as follows. The space of invertible
> matrices forms a mathematical structure known as the "General Linear
> Group", which is differentiable (a Lie group) and has smooth manifold
> structure. As a smooth manifold, it has two, disjoint pieces: those with positive determinant, and those with negative determinant. Thus there can be _no_ smooth path that moves from a matrix with a positive
> determinant to one with a negative determinant. Any path between those
> two matrices must pass through a singular matrix at some point.
> 
> There are many possible paths, but one that has some intuitive
> conceptual meaning is to view a 2D matrix as a projection of a 3D
> matrix. This is how the current interpolation scheme is defined. As
> stated above, any invertible 2D matrix can be extended to a 3D matrix
> with positive determinant. Then, there is a smooth path between such
> matrices, and its projection back down to 2D gives the "natural"
> appearance of a flat object being flipped over. Compare to the current
> behavior of http://dbaron.org/css/test/2010/transition-negative-determinant
> 
> When combined with the quaternion approach for interpolating rotations,
> the result is very intuitive: the first example in the "no skew,
> negative determinant line" simply flips the rectangle around the X axis
> (compared to the current behavior of some kind of combined rotation and
> flip). The behavior of the second example remains the same: a simple
> flip around the Y axis. The next two examples are simple flips across
> the two 45-degree diagonals, instead of more combined rotation-flips.
> 
> There has been some complaint that all of the 3D math is
> over-complicated for 2D transforms, and there is a proposal to simplify
> it: http://lists.w3.org/Archives/Public/www-style/2010Sep/0697.html
> I believe this approach of a guaranteeing a positive determinant of the
> 3D transformation and interpolating rotations with quaternions has a
> similar, simple 2D form. One finds that the projection of just the
> rotation component from 3D to 2D during the animation always has
> det(rot(t))=cos(t*pi) (or cos((1-t)*pi) if going from negative to
> positive determinant), and that the scaling is along a single axis. I
> haven't worked through the details of the algebra to derive an
> expression for the axis for the case of more general rotations.
> 
> Such a simplification for 2D could probably produce a generalization for 3D as well. This would avoid the need to go to 4D rotations in order to handle 3D matrices with a negative determinant, which is a much harder thing to argue for. Unlike 2D, there's no obvious intuitive advantage of the result that would justify the complexity, and a solution described analogous generality to the one presented here becomes much more complex for 4D than it is for 3D.
>

Received on Tuesday, 30 August 2011 01:19:35 UTC