To average together quaternions in a well-defined way, the eigendecomposition method of Markley et al. may be used, from Averaging Quaternions, Journal of Guidance, Control, and Dynamics, 30(4):1193-1196, June 2007, Eqs. (12) and (13).

However, if a set of all quaternions are close to each other (meaning that they represent very similar rotations), then element-wise averaging of the quaternions followed by normalization may produce a sufficiently "central" quaternion. (Elementwise averaging is much faster than the eigendecomposition, which is important for some applications.)

However, the quaternions $\bf{q}$ and $\bf-{q}$ represent the same rotation (sometimes called the "double cover issue" of quaternions), so element-wise averaging cannot be applied without first somehow making sure that any quaternions that are to be averaged lie within the same "half" of the rotation group SO(3).

There are several possible methods for "standardizing" each quaternion in a set of quaternions so that the double-cover issue is not a problem, and I wrote about these in this answer, but I am not sure which of these methods is correct (or optimal, and under what assumptions). Some possible methods for standardizing all quaternions ${\bf q}_i \in Q$ (while ensuring that each quaternion still represents the same rotation) include the following:

- If the $w$ component is negative, negate the quaternion (i.e. replace ${\bf q}_i$ with $-{\bf q}_i$), so that the $w$ component is positive for all quaternions in the set $Q$.
- Take the dot product of ${\bf q}_1$ with all subsequent quaternions ${\bf q}_i$, for $2 \le i \le N$, and negate any of the subsequent quaternions whose dot product with ${\bf q}_i$ is negative.
- For each quaternion, measure the angle of rotation about the rotation axis of the quaternion, and normalize it so it always rotates the "short way around", i.e. such that $-\pi \le \theta \le \pi$. If it rotates the "long way around", i.e. $\theta \lt -\pi$ or $\theta \gt \pi$, then negate the quaternion.

These sometimes produce the same result, but they all produce different results in some cases (i.e. they all can negate different quaternions in a set of quaternions) -- therefore they are not equivalent.

**What is the best way to deal with quaternions in a standardized way in order to overcome the double cover issue in situations like this?**

Note that it is not just element-wise averaging of quaternions that can cause the double cover issue to affect the results. Another example is the swing-twist decomposition: in a naive implementation, the recovered rotation component around a given axis can represent either a rotation "the short way around" or a rotation "the long way around", which can lead to some unexpected or unstable results if you care only about the rotation about the axis, not the full quaternion.