19

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)..

Hannesh
  • 735
  • Google gives http://en.wikipedia.org/wiki/Generalized_quaternion_interpolation Someone who knows more may be able to expand. – Peter Taylor Sep 01 '11 at 09:50
  • @PeterTaylor that Wikipedia page no longer exists, but probably the one you want is https://en.wikipedia.org/wiki/Slerp . See also my new answer below for more info. – Luke Hutchison Aug 25 '20 at 11:07
  • For dealing with the double cover issue specifically (i.e. ${\bf q}$ and $-{\bf q}$ representing the same rotation), see this question and answer: https://math.stackexchange.com/questions/3888504 . For how to apply this to quaternion averaging, see my answer to your question: https://math.stackexchange.com/a/3435296/365886 – Luke Hutchison Nov 08 '20 at 04:45
  • related: https://stackoverflow.com/questions/12374087/average-of-multiple-quaternions – Felipe G. Nievinski Apr 04 '22 at 07:40

6 Answers6

12

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 &gt; 0 &amp;&amp; quats[i].dot(quats[0]) &lt; 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();
}
Luke Hutchison
  • 391
  • 2
  • 12
4

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.

adrelino
  • 41
3

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.

Tim Seguine
  • 1,280
  • The solution by averaging canonical forms looks like jonathon's solution to his question : http://stackoverflow.com/a/12439567/1200764 – WillC Apr 07 '17 at 08:22
1

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

Community
  • 1
WillC
  • 111
1

I now provide the standard and weighted averaging as well as $L_p$ median of quaternions at: http://tbirdal.blogspot.com/

Tolga Birdal
  • 1,082
  • 9
  • 18
  • Wondering what are you thinking of this publication: "Means and Averaging in the Group of Rotations", Maher Moakher. Looks like Euclidean mean of the rotation matrix can also be an easy solution: https://github.com/lagadic/visp/blob/5cd5a377c510274d75e0b07a880dcad063316e17/modules/core/src/math/transformation/vpHomogeneousMatrix.cpp#L1050-L1083 – Catree Jul 01 '21 at 07:43
  • 1
    I think the difference is in the convexity region. To the best of my knowledge, chordal metrics are not convex beyond a ball of radius $\pi/2$ and strictly convex only in $\pi/4$. For a lot of applications this might suffice as well. But note that, computing proper mens is also not that hard or expensive. – Tolga Birdal Mar 04 '23 at 18:16
0

The Karcher/Fréchet mean algorithm can be used, which for quaternions boils down to:

# initial guess
curr = Quat.identity()

while True:
    # linearize samples at current guess and average
    mu = sum((curr.conjugate() * q).log() for q in qs) / len(qs)


# convergence check
if mu.norm() &lt; eps:
    return curr

# improve guess
curr = curr * Quat.exp(mu)

where log/exp are the quaternion logarithm/exponential.

If you can convert between rotation vectors (where the direction is the rotation axis and the magnitude is the angle) and quaternions you can use that as the exponential, and the logarithm is the reverse operation (both modulo a factor 2 that cancels out).

The above weighting scheme is uniform but other convex combinations can be used instead.

max
  • 331