Я пытаюсь сделать переход от матриц к кватернионам для скелетной анимации в моей программе OpenGL, но я столкнулся с проблемой:
Учитывая ряд единичных кватернионов, мне нужно получить кватернион, который при использовании для преобразования вектора даст вектор, являющийся средним вектором, преобразованным каждым кватернионом индивидуально. (с матрицами я просто добавил бы матрицы вместе и разделил бы на число матриц)
Я попробовал Slerping кватернионов, как было предложено здесь, но это не помогло, что я пытаюсь сделать (модель была искажена), поэтому я просто закончил преобразование векторов каждым кватернионом, а затем выполнил средний (пока не смог найти лучшее решение).
Вопреки распространенному мнению в индустрии компьютерной графики, существует простой алгоритм для решения этой проблемы, которая является надежной, точной и простой, что происходит от аэрокосмической промышленности. Он работает во времени линейным по числу усредненных кватернионов плюс (довольно большой) постоянный коэффициент.
Пусть Q = [a_1 * q_1 a_2 * q_2... a_n * q_n]
Где a_i - вес i-го кватерниона, а q_i - i-й кватернион, усредненный, как вектор-столбец. Q поэтому является матрицей 4xN.
Нормализованный собственный вектор, соответствующий наибольшему собственному значению Q * Q ^ T, является средневзвешенным. Так как Q * Q ^ T является самосопряженным и по крайней мере положительным полуопределенным, то быстрые и надежные методы решения этой собственной задачи доступны. Вычисление матрично-матричного произведения является единственным шагом, который растет с усреднением числа элементов.
См. техническую заметку в "Журнале руководства, контроля и динамики" от 2007 г., которая является краткой статьей этого и других методов. В современную эпоху метод, приведенный выше, делает хороший компромисс для надежности и надежности реализации и уже опубликован в учебниках в 1978 году!
Вот реализация функции 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
К сожалению, это не очень просто сделать, но это возможно. Вот документ, объясняющий математику за ней: http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20070017872_2007014421.pdf
Просмотрите страницу Wiki Unity3D (с кодом): http://wiki.unity3d.com/index.php/Averaging_Quaternions_and_Vectors
Также этот пост: http://forum.unity3d.com/threads/86898-Average-quaternions
Это моя реализация в 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]
Вы не можете добавить кватернионы. Вы можете найти кватернион, который вращается непрерывно между двумя углами, включая полпути. Интерполяция Quaternion известна как "slerp" и имеет страницу wikipedia. Это очень полезный трюк для анимации. В некоторых отношениях slerp является основной причиной использования кватернионов в компьютерной графике.
Существует технический отчет от 2001 года, в котором говорится, что среднее на самом деле довольно хорошее приближение при условии, что кватернионы лежат близко друг к другу. (для случая -q = q вы можете просто перевернуть те, которые указывают в другом направлении, предварительно умножив их на -1, так что все кватернионы включали жизнь в одну и ту же половину сферы.
Еще лучший подход представлен в этой статье с 2007 года, в которой используется SVD. Это тот же документ, на который ссылался Натан. Я хотел бы добавить, что существует не только С++, но и реализация Matlab. От выполнения теста script, который поставляется с кодом 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. Отклонить кватернион до усреднения, если его точечный продукт с предыдущей суммой отрицателен. 2. Нормализовать средний кватернион, конец усреднения, если ваша библиотека работает с единичными кватернионами.
Средний кватернион будет представлять приблизительно среднее вращение (максимальная погрешность около 5 градусов).
ПРЕДУПРЕЖДЕНИЕ: Средняя матрица разной ориентации может быть нарушена, если вращение слишком отличается.
Кватернионы не являются идеальным набором степеней свободы для вращения при вычислении неограниченного среднего.
Вот что я использую большую часть времени (
[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;
}'
Поскольку здесь есть разные подходы, я написал скрипт Matlab для их сравнения. Эти результаты, по-видимому, позволяют предположить, что простого усреднения и нормализации кватернионов (подход из единственной вики, называемый здесь 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://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']);