我正在尝试在我的 OpenGL 程序中将骨骼动画从矩阵切换到四元数，但我遇到了一个问题：
给定许多单位四元数，我需要得到一个四元数，当它用于转换向量时，将给出一个向量，该向量是每个四元数单独转换的向量的平均值。（对于矩阵，我只需将矩阵相加并除以矩阵的数量）
问问题
36445 次
14 回答
63
与计算机图形行业的普遍看法相反，有一种直接的算法可以解决这个问题，该算法来自航空航天工业，它强大、准确且简单。它在被平均的四元数加上一个（较大的）常数因子的时间线性运行。
令Q = [ a 1 q 1 , a 2 q 2 , ..., a n q n ],
其中a i是第i个四元数的权重，q i是被平均的第i个四元数，作为列向量。 Q因此是一个 4× N矩阵。
QQ T的最大特征值对应的归一化特征向量就是加权平均。由于QQ T是自伴随的并且至少是半正定的，因此可以使用快速且稳健的方法来解决该特征问题。计算矩阵-矩阵乘积是唯一随着被平均的元素数量而增长的步骤。
请参阅2007 年的《制导、控制和动力学杂志》中的此技术说明，这是此方法和其他方法的摘要论文。在现代，我上面引用的方法在实现可靠性和健壮性方面做出了很好的权衡，并且已经在 1978 年的教科书中发表！
于 2014-12-10T21:11:26.017 回答
16
不幸的是，这并不是非常简单，但它是可能的。这是一份解释其背后数学原理的白皮书：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
于 2013-11-27T17:57:16.740 回答
11
这是我用来平均四元数以进行方向估计的 MATLAB 函数的实现。将 MATLAB 转换为任何其他语言都很简单，除了这种特殊方法 (Markley 2007) 需要计算特征向量和特征值。有许多库（包括 Eigen C++）可以为您做到这一点。
您可以阅读文件的描述/标题以查看原始论文中的数学。
matlab 文件取自http://www.mathworks.com/matlabcentral/fileexchange/40098-tolgabirdal-averaging-quaternions：
% 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
于 2015-03-28T10:05:41.380 回答
4
我尝试按照这里的建议对四元数进行 Slerping，但这对我想要做的事情不起作用（模型失真），所以我只是最终通过每个四元数转换向量，然后做一个平均值（直到我能找到一个更好的解决方案）。
于 2012-09-15T17:27:56.083 回答
3
您不能添加四元数。你可以做的是找到一个在两个角度之间连续旋转的四元数，包括中途。四元数插值被称为“slerp”并且有一个维基百科页面。这是一个非常有用的动画技巧。在某些方面，slerp 是在计算机图形学中使用四元数的主要原因。
于 2012-09-14T19:34:47.163 回答
3
这是我在 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]
于 2018-04-06T10:39:06.290 回答
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;
于 2014-10-20T13:16:07.373 回答
1
使用四元数，您可以做同样的事情，但需要进行少量修正： 1. 如果四元数与先前总和的点积为负，则在平均之前取反四元数。2. 如果您的库使用单位四元数，则归一化平均四元数，即平均的结束。
平均四元数将代表近似平均旋转（最大误差约 5 度）。
警告：如果旋转太不同，则可能会破坏不同方向的平均矩阵。
于 2014-03-26T15:29:41.620 回答
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;
}`
于 2018-09-30T06:54:27.550 回答
1
由于这里有不同的方法，我编写了一个 Matlab 脚本来比较它们。这些结果似乎表明，对于simple_average四元数足够相似且可以接受小偏差的情况，简单地对四元数进行平均和归一化（来自统一 wiki 的方法，称为此处）可能就足够了。
这是输出：
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://stackoverflow.com/a/29315869/1221661
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']);
于 2019-02-05T16:47:55.967 回答
1
Markley 的解决方案最简单的实现（在 Python 中使用 Numpy>=3.6）是：
import numpy as np
def q_average(Q, W=None):
if W is not None:
Q *= W[:, None]
eigvals, eigvecs = np.linalg.eig(Q.T@Q)
return eigvecs[:, eigvals.argmax()]
其中Q大小为 N×4。得到的四元数已经标准化。
在这种情况下，权重默认等于 1。否则，您可以给出大小为 N 的权重列表（每个四元数一个。）
就是这样。
于 2020-04-03T14:10:40.663 回答
1
在此处查看我对加权平均以及四元数 Lp 中位数的解决方案。
于 2020-04-17T21:50:31.547 回答
0
这是一个 GitHub 存储库，其中实现了这个建议的算法 :) https://github.com/christophhagen/averaging-quaternions
感谢和感谢 christophhagen ofc ;)
于 2020-07-24T13:56:52.907 回答
0
还有另一种基于 slerp 的方法，这可能是您在平均四元数时真正想要的。
让我们首先将其与基于特征分析的平均进行比较：
考虑两个四元数 A 和 B 的平均值，对应于围绕 X 轴旋转 0° 和 90° 度，权重 w_A = 2 和 w_B = 1。预期的加权平均值应对应于围绕 X 轴旋转 30° . 基于 Slerp 的加权平均值将返回预期值。基于特征分析的加权平均将返回 26.56° 的旋转。
基于特征分析的方法将返回最小化相应旋转矩阵的 Frobenius 范数的四元数。基于 slerp 的方法将改为根据四元数计算 3D 旋转空间中的平均值。
import math
import numpy as np
import quaternion # (pip install numpy-quaternion)
d2r = math.pi/180
r2d = 180/math.pi
def recover_angle(mat):
return np.arctan2(mat[1,0], mat[0,0])
# ground truth
angles = np.array([0,90])
weights = np.array([2,1])
mean_angle = np.sum(angles*(weights/weights.sum()))
quaternion_mean = quaternion.from_euler_angles(mean_angle*d2r,0,0)
# eigen analysis
Q = quaternion.as_float_array(
[
(weight**0.5) * quaternion.from_euler_angles(0,0,angle*d2r)
for angle,weight
in zip(angles,weights)
]
).T
QQT = Q @ Q.T
eigval,eigvec = np.linalg.eig(QQT)
quaternion_mean_eig = quaternion.from_float_array( eigvec[:,np.argmax(eigval)] )
# slerp based
def slerp_avg(quaternions, weights):
# welford's mean in terms of linear mix operations:
toqt = quaternion.from_float_array
mix = lambda a,b,k: quaternion.slerp(toqt(a),toqt(b),0,1,k)
wmean, wsum, num = quaternions[0], weights[0], len(weights)
for i in range(1,num):
wsum += weights[i]
k = (weights[i]/wsum)
# wmean := wmean*(1-k) + quaternions[i]*k
wmean = mix(wmean,quaternions[i],k)
return wmean
quaternion_mean_slerp = slerp_avg(quaternion.as_float_array(
[quaternion.from_euler_angles(0,0,angle*d2r) for angle in angles]
), weights)
mat_mean = quaternion.as_rotation_matrix(quaternion_mean)
mat_mean_eig = quaternion.as_rotation_matrix(quaternion_mean_eig)
mat_mean_slerp = quaternion.as_rotation_matrix(quaternion_mean_slerp)
print("expected", recover_angle(mat_mean)*r2d)
print("eigen", recover_angle(mat_mean_eig)*r2d)
print("slerp", recover_angle(mat_mean_slerp)*r2d)
输出：
expected 29.999999999999996
eigen 26.565051177077994
slerp 30.00000000000001
您可以在 https://github.com/xaedes/average_affine_transform_mat/找到一个零依赖 C++ 库
于 2021-11-13T21:08:16.087 回答