2. Dual Quaternion Algebra

While quaternion algebra is very well suited for describing the orientation of a rigid body, it does not contain any information on the translational part of the motion. Translation and orientation are decoupled. Sometimes, however, one would like to have a unique entity that describes both. This is for example the case if one wants to nest Transform nodes within a script rather than only orientations. For this purpose, the VRML 97 specification contain the VrmlMatrix object which allows the designer to access the browser's mathematical description of the transformation within a transform node.

In this chapter we want to discuss a different, more sophisticated, description of rigid body transformations called dual quaternions. Dual quaternions are known since the early days of quaternion algebra in the 18th century. However, their use in Computer Animation is rather restricted. One reason for this is that the theory behind this concept is considered to be quite involved. Nevertheless, it will turn out that dual quaternion algebra does not only provide a useful way to address to many standard tasks in VRML, but also allows a relatively simple approach to the much more complex problem of Inverse Kinematics covered in the next chapter.

2.1. Example project: 3d slider puzzle

The project we are going to use for this chapter is the development of a 3d slider puzzle, also known as Sam Loyd grid puzzle, as shown in the right window. By clicking on one of the pieces next to the hole, the piece will change its position, thereby either sliding or rotating 90 degrees depending on wether the piece is on the same face as the hole or not.

Let us again first look at some of the VRML code. The puzzle consists of a total of 23 pieces, which are only distinguished by its texture and initial position. 

PROTO Element [ exposedField SFVec3f translation 0 0 0
                exposedField SFRotation rotation 0 0 1 0
                exposedField MFString textureUrl "dummyTexture.gif"
                exposedField SFBool enabled TRUE
                exposedField SFBool timeEnabled TRUE
                eventOut SFFloat fraction_changed
                eventOut SFBool isActive ] {
  Transform {
    translation IS translation
    rotation IS rotation
    children [
      DEF TS TouchSensor {
        enabled IS enabled
      }
      Shape {
        appearance Appearance {
          material Material {
            diffuseColor 0.8 0.4 0
          }
          texture ImageTexture {
            url "wood.jpg"
          }
        }
        geometry IndexedFaceSet {
          convex TRUE
          coord Coordinate {
            point [
              0.2 2 2, 0.2 2 -2, 0.2 -2 -2, 0.2 -2 2,
              -0.2 2 2, -0.2 2 -2, -0.2 -2 -2, -0.2 -2 2
            ]
          }
          coordIndex [
            4 5 6 7 -1
            0 1 5 4 -1
            1 2 6 5 -1
            2 3 7 6 -1
            3 0 4 7 -1
          ]
        }
      }
      Shape {
        appearance Appearance {
          material Material {
            diffuseColor 1 1 1
            emissiveColor 0.1 0.1 0.1
          }
          texture ImageTexture {
            url IS textureUrl
          }
        }
        geometry IndexedFaceSet {
          convex TRUE
          coord Coordinate {
            point [
              0.2 2 2, 0.2 2 -2, 0.2 -2 -2, 0.2 -2 2
            ]
          }
          coordIndex [
            3 2 1 0 -1
          ]
        }
      }
    ]
  }
  DEF TIMER TimeSensor {
    cycleInterval 2
    fraction_changed IS fraction_changed
    isActive IS isActive
    enabled IS timeEnabled
  }
  ROUTE TS.touchTime TO TIMER.startTime
}
The design is similar to the subcubes of the Rubic's cube. In this example, every gamepiece includes its own timesensor. In addition to the translation and rotation fields the definition contains a textureURL fields which holds the URL to the picture shown on the piece, an eventOut fraction_changed which will be used in order to move the piece and three additional fields for handling the sensor status. Initially only rotation, translation and textureURL have to be set. 
DEF E04 Element {
  translation 2.1 4.5 2.1
  rotation 0 0 1 1.5708
  textureUrl "plate04.jpg"
}
The overall situation is basically the same we encountered in the Rubic's Cube example. There is, however, one major difference. In this case the transformations we have to nest are not only orientation changes. They cause changes in both, the translation field as well as the rotation field of the individual elements. Therefore, we have to keep track of both.

2.2. Dual numbers

It might be surprising, but in order to prepare dual quaternion algebra, we have to make a rather big stretch. The following is not exactly a shortcut, but it will help to develop the theory properly and will also improve the understanding of our approach to Inverse Kinematics.

Before we can do anything related to geometry we first have to introduce the concept of dual numbers. Basically, dual numbers are similar to complex numbers. They consist of two real numbers, the real part and the dual part which are connected by a new unit. In case of complex numbers this is the complex unit i with the property

i2 = -1,

in case of dual numbers it is the so-called dual unit e satisfying

e2 = 0.

In order to distinguish between dual and real numbers, we are going to use the following notation

a = a + e ae

where a is the real part and ae the dual part of the dual number a.

It is important to know, how dual numbers influence the way we are have to handle functions. In particular we will need to know how to calulate sin(a) and cos(a) for a dual number a. It is not difficult to prove, but goes beyond the scope of this tutorial, that any analytic function f(t) satisfies

f ( a ) = f ( a ) + e f' ( ae ),

where f' denotes the first derivative of f. We therefore obtain

sin ( a ) = sin ( a ) + e cos ( ae ) , and
cos ( a ) = cos ( a ) - e sin ( ae ).

Before we are able to define dual quaternions, we need one more thing.

2.3. Line coordinates

There are several ways to describe lines in 3d. One can for instance use the parametric form of a line, or one can describe it as the intersection of two planes, which in mathematical terms is the set of all points satifying two linear equations. In the history of mathematics, it has always been considered an advantage if one can describe a geometric entity with one single set of coodrinates. The question therefore is, wether it is possible to describe a straight line, which for the sake of simplicity may be considered infinitely long, can be described by one simple set of coordinate. The answer is yes and these coordinates are called line coordinates.

The definition is rather simple. Recall that we used homogeneous coordinates to describe points. Let us consider two points O=( o,ow ) and P=( p,pw ). The line coordinates of the line L connecting O and P are a set of six coordinates collected in a pair of vecors v and ve. They are defined by

L = ( v , ve ) = ( ow p - pw o , o × p ).

Notice that the line coordinates satisfy Plücker's identity

v · ve = 0,

which means that these vectors are perpendicular to each other. It can be shown that two vectors represent the coordinates of a line if and only if they satify this identity. Obviously, v is the direction vector of the line. The vector ve is the so-called moment vector.

The line coordinates of a line are not unique, since the homogeneous coordinates we started with were not unique either. Similar to the homogeneous coordinates of a point the line coordinates can be multiplied by any real number. Therefore they are often called homogeneous line coordinates. It is common to switch to normalized line coordinates by dividing both vectors by the length of the direction vector v. Normalized line coordinates are distinguished by the property that the direction vector has length 1. If v=(0,0,0) the line is at infinity.

The final trick is to collect the two vectors describing the line coordinates of a line in a single dual vector

v = v + e ve.

Using dual vectors in order to represent lines will have a significant effect on how we derive dual quaternion algebra.

2.4. Dual quaternion algebra

Recall that we used the rotation vector v and the rotation angle a to define the quaternion. In analogy to orientations where each orientation can be represented by a rotation, the displacement of a rigid body as implied by a Transform node can be represented by a screw motion. Let v be the normalized dual vector which describes the screw axis. Furthermore let p be the pitch of this screw motion and a the rotation angle as in the rotation case. Using dual numbers we first define the so-called dual angle

a = a + e a p.

Notice that ap is the amount of translation in direction of the screw axis. In case the rotation angle is zero, we have a translation in the direction of v. If d denotes the amount of the translation, the corresponding dual angle simply reads

a = e d.

We are finally ready to define the dual quaternion Q representing a displacement of a rigid body by

Q = Q + e Qe = ( q , qw ) = ( sin(a/2) v , cos(a/2) ).

It can be split into a real part Q and a dual part Qe. This definition includes pure rotations and pure translations. In the first case the dual part of the dual angle is zero, in case of a translation the real part is zero.

If we compare this definition to the definition of a quaternion as used in the previous section, we see that the only difference is that we exchanged a vector by a dual vector and a real number by a dual number. In fact, most formulas that can be derived for rotations can be translated into the spatial case simply by making this change from real variables to dual ones. In particlar the nesting of two dual quaternions reads

O = P * Q = ( pw q + qw p + p × q , pw qw - p · q ).

Here the dual quaternion O describes the transformation that results from nesting the transformations associated with P and Q. The quaternion multiplication is exactly the same we used in the previous chapter. Only this time all numbers are dual instead of real.

2.5. Implementing dual quaternion algebra in VRML

Since a dual quaternion consists of 8 parameters we are going to use an MFFloat fields with 8 elements for storage. The following function converts the data from a Transform node into a dual quaternion. The rotation field is assumed to be stored in the first four and the translation field in the following three elements of array. The formula that is used in the return statement can simply be derived by storing rotation and translation in different dual quaternions and applying the quaternion product. Note that the rotation is the first transformation, the translation is applied afterwards. 
function toDualQuaternion(array) {
  len=Math.sqrt(array[0]*array[0]+array[1]*array[1]+array[2]*array[2]);
  if(len<0.01) {
    ax=0.0; ay=0.0; az=0.0; s=0.0; c=1.0;
  } else {
    ax=array[0]/len; ay=array[1]/len; az=array[2]/len;
    s=Math.sin(array[3]/2); c=Math.cos(array[3]/2);
  }
  vx=array[4]/2; vy=array[5]/2; vz=array[6]/2;
  return new MFFloat(ax*s,ay*s,az*s,c,
                     vx*c+(vy*az-vz*ay)*s,vy*c+(vz*ax-vx*az)*s,
                     vz*c+(vx*ay-vy*ax)*s,-(ax*vx+ay*vy+az*vz)*s);
}
The second function transforms a dual quaternion into a MFFloat field where rotation and translation are stored in the above manner. In order to make our live easier we first determine the rotation and apply the inverse rotation in order to extract the translation field. 
function fromDualQuaternion(quat) {
  s=Math.sqrt(quat[0]*quat[0]+quat[1]*quat[1]+quat[2]*quat[2]+quat[3]*quat[3]);
  if(s<0.01) {
    return new MFFloat(0,0,1,0,0,0,0,0);
  } else {
    qu=new MFFloat(quat[0]/s,quat[1]/s,quat[2]/s,quat[3]/s,
                   quat[4]/s,quat[5]/s,quat[6]/s,quat[7]/s);
    rotPart=new MFFloat(-qu[0],-qu[1],-qu[2],qu[3],0,0,0,0);
    transPart=multiply(qu,rotPart);
    len=Math.sqrt(qu[0]*qu[0]+qu[1]*qu[1]+qu[2]*qu[2]);
    if(len<0.01) {
      return new MFFloat(0,0,1,0,2*transPart[4],2*transPart[5],2*transPart[6],0);
    } else {
      return new MFFloat(qu[0]/len,qu[1]/len,qu[2]/len,2*Math.acos(qu[3]),
                         2*transPart[4],2*transPart[5],2*transPart[6],0);
    }
  }
}
Finally, we need a function for multiplying dual quaternions. 
function multiply(a,b) {
  return new MFFloat(a[3]*b[0]+a[0]*b[3]+a[1]*b[2]-a[2]*b[1],
                     a[3]*b[1]+a[1]*b[3]+a[2]*b[0]-a[0]*b[2],
                     a[3]*b[2]+a[2]*b[3]+a[0]*b[1]-a[1]*b[0],
                     a[3]*b[3]-a[0]*b[0]-a[1]*b[1]-a[2]*b[2],
                     a[7]*b[0]+a[4]*b[3]+a[5]*b[2]-a[6]*b[1]+
                       a[3]*b[4]+a[0]*b[7]+a[1]*b[6]-a[2]*b[5],
                     a[7]*b[1]+a[5]*b[3]+a[6]*b[0]-a[4]*b[2]+
                       a[3]*b[5]+a[1]*b[7]+a[2]*b[4]-a[0]*b[6],
                     a[7]*b[2]+a[6]*b[3]+a[4]*b[1]-a[5]*b[0]+
                       a[3]*b[6]+a[2]*b[7]+a[0]*b[5]-a[1]*b[4],
                     a[7]*b[3]-a[4]*b[0]-a[5]*b[1]-a[6]*b[2]+
                       a[3]*b[7]-a[0]*b[4]-a[1]*b[5]-a[2]*b[6]);
}
Since the current specification of ECMAScript does not provide objects we have to code the multiplication formula explicitly.

2.6. Implementing the slider puzzle

The actual implementation of the slider puzzle is similar to the Rubic's cube. For each of the pieces we have to store their current translation and rotation fields in two arrays called transArray and rotArray. Furthermore we again use a lookup table to determine the position of the piece on the cube. The following function is reponsible for the motion of the piece with index pieceInd. 
function calcCurPos(value) {
  if(onSameFace==TRUE) {
    curPos[0]=transArray[pieceInd][0]+4.2*value*transVec[0];
    curPos[1]=transArray[pieceInd][1]+4.2*value*transVec[1];
    curPos[2]=transArray[pieceInd][2]+4.2*value*transVec[2];
    curRot=rotArray[pieceInd];
  } else {
    baseArray=new MFFloat(rotArray[pieceInd][0],rotArray[pieceInd][1],
                          rotArray[pieceInd][2],rotArray[pieceInd][3],
                          transArray[pieceInd][0],transArray[pieceInd][1],
                          transArray[pieceInd][2],0);
    relArrayA=new MFFloat(0,0,1,0,
                          -2.1*transVec[0],-2.1*transVec[1],-2.1*transVec[2],0);
    relArrayB=new MFFloat(rotVec[0],rotVec[1],rotVec[2],1.5708*value,
                          2.1*transVec[0],2.1*transVec[1],2.1*transVec[2],0);
    qu=toDualQuaternion(baseArray);
    baseQuat=new MFFloat(qu[0],qu[1],qu[2],qu[3],qu[4],qu[5],qu[6],qu[7]);
    qu=toDualQuaternion(relArrayA);
    relQuatA=new MFFloat(qu[0],qu[1],qu[2],qu[3],qu[4],qu[5],qu[6],qu[7]);
    qu=toDualQuaternion(relArrayB);
    relQuatB=new MFFloat(qu[0],qu[1],qu[2],qu[3],qu[4],qu[5],qu[6],qu[7]);
    qu=multiply(relQuatB,relQuatA);
    relQuat=new MFFloat(qu[0],qu[1],qu[2],qu[3],qu[4],qu[5],qu[6],qu[7]);
    // CP1.0.2/IRIX workaround!
    res=fromDualQuaternion(multiply(relQuat,baseQuat));
    curRot=new SFRotation(res[0],res[1],res[2],res[3]);
    curPos=new SFVec3f(res[4],res[5],res[6]);
  }
}
The field onSameFace is set true, if the chosen piece and the hole are on the same face. In this case the motion is a pure translation in direction of a vector stored in transVec. Otherwise, the motion is a rotation around an axis which does not pass through the center of the coordinate frame. We therefore first have to translate the piece back such that the center is on the rotation axis. This is done by dual Quaternion relQuatA. Then the piece is rotated and moved back by the dual quaternion relQuatB. Finally, the dual quaternions are multiplied with the quaternion describing the initial transformation and rotational and translational part are extracted.

The details can be found in the source code. To play with the puzzle, click here.