Averaging quaternions

Given multiple quaternions representing orientations, and I want to average them. Each one has a different weight, and they all sum up to one.

How can I get the average of them? Simple multiplication by weights and addition won't work, since it doesn't take into account that (qw, qx, qy, qz) = (-qw, -qx, -qy, -qz)..

I assume you are thinking of unit quaternions and you are using them to represent rotations? If that is the case then here is a paper on the related subject of means and averages in the rotation group. It might not be a very easy read though if you don't understand the notation.

Barring that, Here's what I might try: Pick a canonical form for your quaternions. Then convert each to the canonical form and finally perform your weighted linear combination.

There is a technical report from 2001 which states that the mean is actually quite a good approximation, provided that the quaternions lie close together. (for the case of -q=q, you could just flip the ones that point in the other direction by pre multiplying them by -1, so that all of the quaternions involved life in the same half sphere.

An even better approach is sketched in this paper from 2007, which involves using an SVD.

Averaging together pairs of quaternions

If you're trying to average just two quaternions, you can use quaternion interpolation (slerp) with an interpolation factor of 0.5 to find the "rotational midpoint" between the two quaternions.

(Note that Eq. 19 of the Markley et al. paper cited below also gives a direct / closed form solution for the average of two quaternions. This may or may not give the same result as slerp with an interpolation factor of 0.5, I haven't looked at it closely.)

Averaging together more than two quaternions

There are at least two methods for averaging together more than two quaternions.

1. Direct averaging (fast/approximate)

If quaternions represent similar rotations, and the quaternions are normalized, and a correction has been applied for the "double-cover problem", then the quaternions can be directly averaged and then the result normalized again, treating them as 4-dimensional vectors, to produce a quaternion representing a roughly-average rotation.

The double-cover problem is where $q$ and $-q \,$ represent the same rotation, as the original poster pointed out. This problem can be addressed by negating some subset of the quaternions to ensure that the dot product of any pair of quaternions is positive (or, w.l.o.g., negative).

Note that the more dissimilar the quaternions involved, the more inaccurate this method will be, since this is only an approximation to the eigendecomposition method -- so you only want to average together quaternions that are already fairly similar in terms of the rotations they represent.

I measured this method at 2.56x faster than the eigendecomposition method.

(Syntax used is for the JOML Java library)

public static Quaterniond averageApprox(Quaterniond[] quats, double[] weights) {
    if (weights != null && quats.length != weights.length) {
        throw new IllegalArgumentException("Args are of different length");
    }
    Quaterniond qAvg = new Quaterniond(0.0, 0.0, 0.0, 0.0);
    for (int i = 0; i < quats.length; i++) {
        Quaterniond q = quats[i];
        double weight = weights == null ? 1.0 : weights[i];
        
        // Correct for double cover, by ensuring that dot product
        // of quats[i] and quats[0] is positive
        if (i > 0 && quats[i].dot(quats[0]) < 0.0) {
            weight = -weight;
        }
        
        qAvg.add(weight * q.x, weight * q.y, weight * q.z, weight * q.w);
    }
    return qAvg.normalize();
}

2. Eigendecomposition of sum of outer products (maximum likelihood method)

This method is based on Markley et al., Averaging Quaternions, Journal of Guidance, Control, and Dynamics, 30(4):1193-1196, June 2007, Eqs. (12) and (13). The authors proved that the solution is unique, and that it corresponds to a maximum likelihood criterion. (tbirdal's implementation of this same algorithm, plus some more powerful variants, were linked from a couple of other answers.)

(Syntax used is for the Java JOML library and the Apache Commons Math Linear Algebra library)

public static Quaterniond average(Quaterniond[] quats, double[] weights) {
    if (weights != null && quats.length != weights.length) {
        throw new IllegalArgumentException("Args are of different length");
    }
    RealMatrix accum = MatrixUtils.createRealMatrix(4, 4);
    for (int i = 0; i < quats.length; i++) {
        Quaterniond q = quats[i];
        double weight = weights == null ? 1.0 : weights[i];
        RealVector qVec = MatrixUtils.createRealVector(
                new double[] { q.x, q.y, q.z, q.w });
        RealMatrix qOuterProdWeighted = qVec.outerProduct(qVec)
                .scalarMultiply(weight);
        accum = accum.add(qOuterProdWeighted);
    }
    EigenDecomposition eigDecomp = new EigenDecomposition(accum);
    double[] ev0 = eigDecomp.getEigenvector(0).toArray();
    return new Quaterniond(ev0[0], ev0[1], ev0[2], ev0[3]).normalize();
}

This answer by Gouda gives the working of the Matlab implementation by Tolga Birdal, which is well-commented.