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 = [a1q1, a2q2, ..., anqn],
where ai is the weight of the ith quaternion, and qi 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 QQT is the weighted average. Since QQT 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!
与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…