"Average" of multiple quaternions?

I'm trying to make the switch from matrices to quaternions for skeletal animation in my OpenGL program, but I've encountered a problem:

Given a number of unit quaternions, I need to get a quaternion that when used to transform a vector will give a vector that is the average of the vector transformed by each quaternion individually. (with matrices I would simply add the matrices together and divide by the number of matrices)

Solution 1:

Contrary to popular belief in the computer graphics industry, there is a straightforward algorithm to solve this problem which is robust, accurate, and simple that comes from the aerospace industry. It runs in time linear in the number of quaternions being averaged plus a (largish) constant factor.

Let Q = [a₁q₁, a₂q₂, ..., a_nq_n],

where a_i is the weight of the ith quaternion, and q_i is the ith quaternion being averaged, as a column vector. Q is therefore a 4×N matrix.

The normalized eigenvector corresponding to the largest eigenvalue of QQ^T is the weighted average. Since QQ^T is self-adjoint and at least positive semi-definite, fast and robust methods of solving that eigenproblem are available. Computing the matrix-matrix product is the only step that grows with the number of elements being averaged.

See this technical note in the Journal of Guidance, Control, and Dynamics from 2007, which is a summary paper of this and other methods. In the modern era, the method I quoted above makes a good tradeoff for implementation reliability and robustness, and was already published in textbooks in 1978!

Solution 2:

Unfortunately it isn't terribly simple to do, but it is possible. Here's a whitepaper explaining the math behind it: http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20070017872_2007014421.pdf

Check out the Unity3D Wiki page (The below code sample is from the same article): http://wiki.unity3d.com/index.php/Averaging_Quaternions_and_Vectors

//Get an average (mean) from more then two quaternions (with two, slerp would be used).
//Note: this only works if all the quaternions are relatively close together.
//Usage: 
//-Cumulative is an external Vector4 which holds all the added x y z and w components.
//-newRotation is the next rotation to be added to the average pool
//-firstRotation is the first quaternion of the array to be averaged
//-addAmount holds the total amount of quaternions which are currently added
//This function returns the current average quaternion
public static Quaternion AverageQuaternion(ref Vector4 cumulative, Quaternion newRotation, Quaternion firstRotation, int addAmount){

    float w = 0.0f;
    float x = 0.0f;
    float y = 0.0f;
    float z = 0.0f;

    //Before we add the new rotation to the average (mean), we have to check whether the quaternion has to be inverted. Because
    //q and -q are the same rotation, but cannot be averaged, we have to make sure they are all the same.
    if(!Math3d.AreQuaternionsClose(newRotation, firstRotation)){

        newRotation = Math3d.InverseSignQuaternion(newRotation);    
    }

    //Average the values
    float addDet = 1f/(float)addAmount;
    cumulative.w += newRotation.w;
    w = cumulative.w * addDet;
    cumulative.x += newRotation.x;
    x = cumulative.x * addDet;
    cumulative.y += newRotation.y;
    y = cumulative.y * addDet;
    cumulative.z += newRotation.z;
    z = cumulative.z * addDet;      

    //note: if speed is an issue, you can skip the normalization step
    return NormalizeQuaternion(x, y, z, w);
}

public static Quaternion NormalizeQuaternion(float x, float y, float z, float w){

    float lengthD = 1.0f / (w*w + x*x + y*y + z*z);
    w *= lengthD;
    x *= lengthD;
    y *= lengthD;
    z *= lengthD;

    return new Quaternion(x, y, z, w);
}

//Changes the sign of the quaternion components. This is not the same as the inverse.
public static Quaternion InverseSignQuaternion(Quaternion q){

    return new Quaternion(-q.x, -q.y, -q.z, -q.w);
}

//Returns true if the two input quaternions are close to each other. This can
//be used to check whether or not one of two quaternions which are supposed to
//be very similar but has its component signs reversed (q has the same rotation as
//-q)
public static bool AreQuaternionsClose(Quaternion q1, Quaternion q2){

    float dot = Quaternion.Dot(q1, q2);

    if(dot < 0.0f){

        return false;                   
    }

    else{

        return true;
    }
}

Also this post: http://forum.unity3d.com/threads/86898-Average-quaternions