logo标签列表

多个四元数的“平均值”是什么？

mathquaternions
70

我正在尝试在我的OpenGL程序中将骨架动画的矩阵转换为四元数，但是我遇到了一个问题：

给定一些单位四元数，我需要获得一个四元数，用它来变换一个向量，使得变换后的向量是每个四元数分别变换后向量的平均值。(对于矩阵，我只需要将矩阵相加并除以矩阵数量即可)

5
这不是一个可解决的问题。四元数无法编码三维空间中任意线性变换，它们只能编码正交变换。但是多个正交变换的“平均值”（即：对于每个向量，取其所有变换并求平均）通常不是正交的，因此不能通过四元数获得。可能存在一些更微妙的“平均”概念可以保持正交性，但你得不到你所想要的平均值。 - darij grinberg
1
所有关于四元数的可能算术运算都在这个页面中（http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/index.htm）。如果你在那里找不到，那么...祝你好运。 - Jav_Rock
你有没有考虑过查看开源骨骼动画代码，而不是重复造轮子？例如：http://home.gna.org/cal3d/ - JWWalker
相关：https://math.stackexchange.com/questions/61146/averaging-quaternions - Felipe G. Nievinski
16个回答
78
与计算机图形学界的普遍观点相反，解决这个问题有一个简单、稳健、准确的算法，它来自于航空航天工业。它的运行时间是线性的，与被平均的四元数数量成正比，再加上一个（相当大的）常数因子。
设 Q = [a1q1, a2q2, ..., anqn]，其中 ai 是第 i 个四元数的权重，qi 是要求平均值的第 i 个四元数，表示为列向量。因此，Q 是一个 4×N 的矩阵。
对应于 QQ^T 的最大特征值的归一化特征向量是加权平均值。由于 QQ^T 是自共轭的且至少半正定，因此可以使用快速而稳健的方法解决该特征值问题。计算矩阵乘积是唯一随着被平均元素数量增长而增加的步骤。
请参阅《2007 年导航、控制和动力学杂志》中的技术说明，其中总结了此方法和其他方法。在现代，我引用的方法在实现可靠性和稳健性方面取得了良好的折衷，而且早在 1978 年就已经出版在教科书上了！
15
@GabeHalsmer：这是*不正确的。如果四元数足够相似，则Unity代码是适当的（它仅是平均加负四元数检查）。所引用的论文提供了一种完全不同的解决方案，适用于所有四元数，但计算成本显着更高。 - Martin
1
权重 a_1，a_2，..，a_n 总和是否为 1.0？ - Maksym Ganenko
1
@Maksym Ganenko 我测试了一下，似乎是必要的。 - ntv1000
2
这个答案有错误吗？我会在你的公式中加入$a_i^{1/2}$，以便$QQ^\intercal$可以恢复$a_i$。我这么说是因为我在你的链接http://www.acsu.buffalo.edu/~johnc/ave_quat07.pdf的12号方程中没有看到权重的平方。 - Alexander Soare
1
在Python中，这里是：https://scikit-surgerycore.readthedocs.io/en/stable/_modules/sksurgerycore/algorithms/averagequaternions.html#average_quaternions - ttst
显示剩余3条评论
18

很不幸，这并不是非常简单的事情，但是它是可能实现的。这里有一篇白皮书解释了背后的数学原理：http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20070017872_2007014421.pdf

请查看Unity3D Wiki页面（下面的代码示例来自同一篇文章）：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;
    }
}

还有这篇文章：http://forum.unity3d.com/threads/86898-Average-quaternions

顺便提一下，这是一个 Matlab 实现 https://au.mathworks.com/matlabcentral/fileexchange/40098-tolgabirdal-averaging_quaternions - R.Falque
请问您能否解释一下“AreQuaternionsClose”的逻辑？ - GalDude33
@GalDude33 两个向量的点积也等于它们长度的余弦值乘以它们之间夹角的余弦值（请参阅维基百科点积文章的“几何定义”部分）。长度始终为正，因此点积的符号等于“向量”之间夹角的余弦值的符号，这本身是一个有些复杂的值，但它将告诉您这两个四元数是否... - Jason C
如果您将它们视为4维点而不是四元数，这些点会“指向”相同的方向（如果您考虑每个分量的符号是否被倒转）。这只是对问题“每个分量的符号是否反转”的近似回答。如果有意义的话。这很简单，但我在评论中解释起来有些困难，呵呵。 - Jason C
如果你把四元数想象成在4D空间中的普通4D向量，它们的点积将是负数，如果它们指向相反的方向（如果它们之间的角度大于90度）。是的，这是一个更简单的解释。或许吧。点积中隐藏了很多信息。 - Jason C
13

这里是我用于方位估计的四元数平均值的 MATLAB 函数实现。转换成其他语言很容易，但这种特定方法（Markley 2007）需要计算特征向量和特征值。有许多库（包括 Eigen C++）可以为您完成此操作。

您可以阅读文件的描述/标题以查看原始论文中的数学内容。

从http://www.mathworks.com/matlabcentral/fileexchange/40098-tolgabirdal-averaging-quaternions获得的 Matlab 文件:

% by Tolga Birdal
% Q is an Mx4 matrix of quaternions. weights is an Mx1 vector, a weight for
% each quaternion.
% Qavg is the weightedaverage quaternion
% This function is especially useful for example when clustering poses
% after a matching process. In such cases a form of weighting per rotation
% is available (e.g. number of votes), which can guide the trust towards a
% specific pose. weights might then be interpreted as the vector of votes 
% per pose.
% Markley, F. Landis, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman. 
% "Averaging quaternions." Journal of Guidance, Control, and Dynamics 30, 
% no. 4 (2007): 1193-1197.
function [Qavg]=quatWAvgMarkley(Q, weights)

% Form the symmetric accumulator matrix
A=zeros(4,4);
M=size(Q,1);
wSum = 0;

for i=1:M
    q = Q(i,:)';
    w_i = weights(i);
    A=w_i.*(q*q')+A; % rank 1 update
    wSum = wSum + w_i;
end

% scale
A=(1.0/wSum)*A;

% Get the eigenvector corresponding to largest eigen value
[Qavg, ~]=eigs(A,1);

end
2
为什么需要使用 for i = 1:M 循环？你不能只是加权 Q 的值然后计算 A=Q'*Q 吗？ - fdermishin
2
@fdermishin：是的，我测试过了，它产生的值相同，误差大约为10^-15。而且它可能会更快。这是代码：Q = (weights .* Q) ./ sum(weights); A = transpose(Q) * Q; - Fritz
为什么要应用缩放，它只影响特征值的大小，而不影响特征向量。 ％缩放 A =（1.0 / wSum）* A; - minorlogic
5
这是我用Python实现的Tolga Birdal算法：
import numpy as np

def quatWAvgMarkley(Q, weights):
    '''
    Averaging Quaternions.

    Arguments:
        Q(ndarray): an Mx4 ndarray of quaternions.
        weights(list): an M elements list, a weight for each quaternion.
    '''

    # Form the symmetric accumulator matrix
    A = np.zeros((4, 4))
    M = Q.shape[0]
    wSum = 0

    for i in range(M):
        q = Q[i, :]
        w_i = weights[i]
        A += w_i * (np.outer(q, q)) # rank 1 update
        wSum += w_i

    # scale
    A /= wSum

    # Get the eigenvector corresponding to largest eigen value
    return np.linalg.eigh(A)[1][:, -1]
7
补充一下，对于大量四元数数据，您可以使用np.linalg.eigh（np.einsum（'ij，ik，i->... jk'，q，q，w））[1] [:，-1]来实现此操作，这比原来的方法快几百倍。累加器的缩放全部在特征值中（不包括规范正交特征向量），因此可以省略。 - reve_etrange
1
@reve_etrange：那行代码是否意味着要替换整个函数，以便您的 q 是上面答案中的 Q（而w 是 weights）？ - ketza
1
好的，现在可以确认它确实替换了整个函数，并且产生了相同的结果，精确到第七位小数左右。 - ketza
@ketz没错！很高兴能帮到你，我认为这是正确的方法。我认为高精度上的差异只是不同的浮点误差，而不是算法上的差异。 - reve_etrange
5

我尝试按照这里的建议对四元数进行旋转，但对于我要做的事情（模型变形）没有效果，因此我最终决定通过每个四元数来变换向量，然后求平均值（直到我能找到更好的解决方案）。

你提供的链接中的slerping方法对我确实有效。然而四元数相对较近。 - Eran
3

您不能直接相加四元数。但是，您可以找到一个四元数，它可以在两个角度之间连续旋转，包括中间位置。四元数插值通常被称为“slerp”，并且有一个维基页面来介绍它。这是动画制作中非常有用的技巧。在某些方面，四元数插值（slerp）是使用计算机图形学中四元数的主要原因。

8
你可以对四元数进行加法运算，但是两个单位四元数的和不会得到一个单位四元数。 - JWWalker
2
但是你可以对它们求平均值，就像Gouda在上面展示的那样 - https://dev59.com/N2ct5IYBdhLWcg3wEJcm#29315869 - WillC
2

有一份来自2001年的技术报告指出，平均数实际上是一个相当好的近似值，只要四元数彼此接近。（对于-q=q的情况，您可以通过在它们之前乘以-1来翻转指向另一个方向的那些四元数，以便涉及的所有四元数都处于同一半球中。）

更好的方法在这篇2007年的论文中概述，其中涉及使用奇异值分解（SVD）。这正是Nathan提到的同一篇论文。我想补充说明的是，不仅有C++实现，还有Matlab实现。通过执行与Matlab代码一起提供的测试脚本，我可以说它对涉及的四元数的小扰动（0.004 *均匀噪声）给出了相当不错的结果：

qinit=rand(4,1);
Q=repmat(qinit,1,10);

% apply small perturbation to the quaternions 
perturb=0.004;
Q2=Q+rand(size(Q))*perturb;
1
这个关于四元数平均的问题在math.stackexchange.com上有相同的提问。 - adrelino
2

点击这里查看我对加权平均和四元数Lp中位数的解决方案。

2
由于这里有不同的方法，我编写了一个Matlab脚本来比较它们。这些结果似乎表明，对四元数进行简单平均和归一化（来自Unity Wiki的方法，在此称为“simple_average”）可能足够处理四元数相似且小偏差可接受的情况。
以下是输出结果：
everything okay, max angle offset == 9.5843
qinit to average: 0.47053 degrees
qinit to simple_average: 0.47059 degrees
average to simple_average: 0.00046228 degrees
loop implementation to matrix implementation: 3.4151e-06 degrees

以下是代码：
%% Generate random unity quaternion
rng(42); % set arbitrary seed for random number generator
M = 100;
qinit=rand(1,4) - 0.5;
qinit=qinit/norm(qinit);
Qinit=repmat(qinit,M,1);

%% apply small perturbation to the quaternions 
perturb=0.05; % 0.05 => +- 10 degrees of rotation (see angles_deg)
Q = Qinit + 2*(rand(size(Qinit)) - 0.5)*perturb;
Q = Q ./ vecnorm(Q, 2, 2); % Normalize perturbed quaternions
Q_inv = Q * diag([1 -1 -1 -1]); % calculated inverse perturbed rotations

%% Test if everything worked as expected: assert(Q2 * Q2_inv = unity)
unity = quatmultiply(Q, Q_inv);
Q_diffs = quatmultiply(Qinit, Q_inv);
angles = 2*acos(Q_diffs(:,1));
angles_deg = wrapTo180(rad2deg(angles));
if sum(sum(abs(unity - repmat([1 0 0 0], M, 1)))) > 0.0001
    disp('error, quaternion inversion failed for some reason');
else
    disp(['everything okay, max angle offset == ' num2str(max(angles_deg))])
end

%% Calculate average using matrix implementation of eigenvalues algorithm
[average,~] = eigs(transpose(Q) * Q, 1);
average = transpose(average);
diff = quatmultiply(qinit, average * diag([1 -1 -1 -1]));
diff_angle = 2*acos(diff(1));

%% Calculate average using algorithm from https://dev59.com/N2ct5IYBdhLWcg3wEJcm#29315869
average2 = quatWAvgMarkley(Q, ones(M,1));
diff2 = quatmultiply(average, average2 * diag([1 -1 -1 -1]));
diff2_angle = 2*acos(diff2(1));

%% Simply add coefficients and normalize the result
simple_average = sum(Q) / norm(sum(Q));
simple_diff = quatmultiply(qinit, simple_average * diag([1 -1 -1 -1]));
simple_diff_angle = 2*acos(simple_diff(1));
simple_to_complex = quatmultiply(simple_average, average * diag([1 -1 -1 -1]));
simple_to_complex_angle = 2*acos(simple_to_complex(1));

%% Compare results
disp(['qinit to average: ' num2str(wrapTo180(rad2deg(diff_angle))) ' degrees']);
disp(['qinit to simple_average: ' num2str(wrapTo180(rad2deg(simple_diff_angle))) ' degrees']);
disp(['average to simple_average: ' num2str(wrapTo180(rad2deg(simple_to_complex_angle))) ' degrees']);
disp(['loop implementation to matrix implementation: ' num2str(wrapTo180(rad2deg(diff2_angle))) ' degrees']);
1

当计算无约束平均值时，四元数不是用于旋转的理想自由度集。

以下是我大部分时间使用的内容（

[MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static Vector3 ToAngularVelocity( this Quaternion q )
    {
        if ( abs(q.w) > 1023.5f / 1024.0f)
            return new Vector3();
            var angle = acos( abs(q.w) );
            var gain = Sign(q.w)*2.0f * angle / Sin(angle);

        return new Vector3(q.x * gain, q.y * gain, q.z * gain);
    }


    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static Quaternion FromAngularVelocity( this Vector3 w )
    {
        var mag = w.magnitude;
        if (mag <= 0)
            return Quaternion.identity;
        var cs = cos(mag * 0.5f);
        var siGain = sin(mag * 0.5f) / mag;
        return new Quaternion(w.x * siGain, w.y * siGain, w.z * siGain, cs);

    }

    internal static Quaternion Average(this Quaternion refence, Quaternion[] source)
    {

        var refernceInverse = refence.Inverse();
        Assert.IsFalse(source.IsNullOrEmpty());
        Vector3 result = new Vector3();
        foreach (var q in source)
        {
            result += (refernceInverse*q).ToAngularVelocity();
        }

        return reference*((result / source.Length).FromAngularVelocity());
    }
     internal static Quaternion Average(Quaternion[] source)
    {
        Assert.IsFalse(source.IsNullOrEmpty());
        Vector3 result = new Vector3();
        foreach (var q in source)
        {
            result += q.ToAngularVelocity();
        }

        return (result / source.Length).FromAngularVelocity();
    }
     internal static Quaternion Average(Quaternion[] source, int iterations)
    {
        Assert.IsFalse(source.IsNullOrEmpty());
        var reference = Quaternion.identity;
        for(int i = 0;i < iterations;i++)
        {
            reference = Average(reference,source);

        }
        return reference;

    }`
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接