Нахождение кватерниона, представляющего вращение от одного вектора к другому
У меня есть два вектора u и v. Есть ли способ найти кватернион, представляющий поворот от u к v?
10 ответов
Quaternion q;
vector a = crossproduct(v1, v2);
q.xyz = a;
q.w = sqrt((v1.Length ^ 2) * (v2.Length ^ 2)) + dotproduct(v1, v2);
Не забудьте нормализовать q.
Ричард прав в том, что не существует уникального вращения, но вышеприведенное должно дать "самую короткую дугу", которая, вероятно, вам и нужна.
Половинное векторное решение
Я нашел решение, которое, как мне кажется, пытался представить Имброндир (хотя и с незначительной ошибкой, которая, вероятно, была причиной того, что у sinisterchipmunk возникли проблемы с его проверкой).
Учитывая, что мы можем построить кватернион, представляющий вращение вокруг оси, вот так:
q.w == cos(angle / 2)
q.x == sin(angle / 2) * axis.x
q.y == sin(angle / 2) * axis.y
q.z == sin(angle / 2) * axis.z
И что точка и перекрестное произведение двух нормализованных векторов:
dot == cos(theta)
cross.x == sin(theta) * perpendicular.x
cross.y == sin(theta) * perpendicular.y
cross.z == sin(theta) * perpendicular.z
Видя, что вращение от u к v может быть достигнуто вращением на тета (угол между векторами) вокруг перпендикулярного вектора, выглядит так, как будто мы можем напрямую построить кватернион, представляющий такое вращение, из результатов точечной и перекрестной произведений; однако в том виде, в котором он стоит, theta = angle / 2, что означает, что это приведет к удвоенному желаемому повороту.
Одно из решений состоит в том, чтобы вычислить вектор на полпути между u и v и использовать точку и скрещенное произведение u и вектор на полпути, чтобы построить кватернион, представляющий вращение, в два раза превышающее угол между u и вектором на полпути, который ведет нас к v!
Существует особый случай, когда u == -v и уникальный вектор на полпути становится невозможным для вычисления. Это ожидаемо, учитывая бесконечно много поворотов "кратчайшей дуги", которые могут привести нас от u к v, и мы должны просто повернуть на 180 градусов вокруг любого вектора, ортогонального к u (или v), в качестве нашего решения в особом случае. Это сделано, беря нормализованное перекрестное произведение u с любым другим вектором, не параллельным u.
Далее следует псевдокод (очевидно, что в действительности особый случай должен учитывать неточности с плавающей запятой - возможно, путем проверки точечных произведений по некоторому порогу, а не по абсолютному значению).
Также обратите внимание, что нет особого случая, когда u == v (кватернион идентичности произведен - проверьте и убедитесь сами).
// N.B. the arguments are _not_ axis and angle, but rather the
// raw scalar-vector components.
Quaternion(float w, Vector3 xyz);
Quaternion get_rotation_between(Vector3 u, Vector3 v)
{
// It is important that the inputs are of equal length when
// calculating the half-way vector.
u = normalized(u);
v = normalized(v);
// Unfortunately, we have to check for when u == -v, as u + v
// in this case will be (0, 0, 0), which cannot be normalized.
if (u == -v)
{
// 180 degree rotation around any orthogonal vector
return Quaternion(0, normalized(orthogonal(u)));
}
Vector3 half = normalized(u + v);
return Quaternion(dot(u, half), cross(u, half));
}
orthogonal Функция возвращает любой вектор, ортогональный данному вектору. Эта реализация использует перекрестное произведение с наиболее ортогональным базисным вектором.
Vector3 orthogonal(Vector3 v)
{
float x = abs(v.x);
float y = abs(v.y);
float z = abs(v.z);
Vector3 other = x < y ? (x < z ? X_AXIS : Z_AXIS) : (y < z ? Y_AXIS : Z_AXIS);
return cross(v, other);
}
Half-Way Quaternion Solution
Это на самом деле решение, представленное в принятом ответе, и оно кажется немного более быстрым, чем векторное решение на полпути (по моим измерениям на 20% быстрее, хотя не верьте мне на слово). Я добавляю это здесь в случае, если другие как я заинтересованы в объяснении.
По сути, вместо того, чтобы вычислять кватернион, используя вектор на полпути, вы можете вычислить кватернион, который приводит к удвоенному требуемому вращению (как подробно описано в другом решении), и найти кватернион на полпути между этим и нулевыми градусами.
Как я объяснял ранее, кватернион для двойного необходимого вращения:
q.w == dot(u, v)
q.xyz == cross(u, v)
И кватернион для нулевого вращения:
q.w == 1
q.xyz == (0, 0, 0)
Вычисление полупериодического кватерниона - это просто вопрос суммирования кватернионов и нормализации результата, как с векторами. Однако, как и в случае с векторами, кватернионы должны иметь одинаковую величину, в противном случае результат будет перекошен в сторону кватерниона с большей величиной.
Кватернион, построенный из точечного и перекрестного произведения двух векторов, будет иметь ту же величину, что и эти произведения: length(u) * length(v), Вместо того, чтобы делить все четыре компонента на этот фактор, мы можем вместо этого увеличить единичный кватернион. И если вам интересно, почему принятый ответ, по-видимому, усложняет ситуацию с помощью sqrt(length(u) ^ 2 * length(v) ^ 2)это потому, что квадрат вектора длины вычисляется быстрее, чем длина, поэтому мы можем сохранить один sqrt расчет. Результат:
q.w = dot(u, v) + sqrt(length_2(u) * length_2(v))
q.xyz = cross(u, v)
А затем нормализуйте результат. Псевдокод следует:
Quaternion get_rotation_between(Vector3 u, Vector3 v)
{
float k_cos_theta = dot(u, v);
float k = sqrt(length_2(u) * length_2(v));
if (k_cos_theta / k == -1)
{
// 180 degree rotation around any orthogonal vector
return Quaternion(0, normalized(orthogonal(u)));
}
return normalized(Quaternion(k_cos_theta + k, cross(u, v)));
}
Поставленная задача не является четко определенной: для данной пары векторов нет единственного поворота. Рассмотрим, например, случай, когда u = <1, 0, 0> и v = <0, 1, 0>. Одно вращение от u до v будет вращением pi / 2 вокруг оси z. Другое вращение от u к v будет вращением pi вокруг вектора <1, 1, 0>.
Почему бы не представить вектор с использованием чистых кватернионов? Лучше, если вы сначала нормализуете их.
q1 = (0 ux uy uz) '
q2 = (0 vx vy vz) '
q1 qrot = q2
Предварительно умножить с q1-1
qrot = q1-1 q2
где q1-1 = q1кон / qнорма
Это можно рассматривать как "левое разделение". Правильное деление, а это не то, что вы хотите:
qrot, right = q2-1 q1
Я не очень хорош в Quaternion. Однако я часами боролся с этим и не мог заставить работать решение Polaris878. Я пытался предварительно нормализовать v1 и v2. Нормализация q. Нормализующий q.xyz. И все же я не понимаю. Результат все еще не дал мне правильный результат.
В конце концов, хотя я нашел решение, которое сделал. Если это поможет кому-то еще, вот мой рабочий (python) код:
def diffVectors(v1, v2):
""" Get rotation Quaternion between 2 vectors """
v1.normalize(), v2.normalize()
v = v1+v2
v.normalize()
angle = v.dot(v2)
axis = v.cross(v2)
return Quaternion( angle, *axis )
Особый случай должен быть сделан, если v1 и v2 параллельны, как v1 == v2 или v1 == -v2 (с некоторым допуском), где я считаю, что решения должны быть Quaternion(1, 0,0,0) (без вращения) или кватернион (0, *v1) (поворот на 180 градусов)
Обобщенное решение
function align(Q, u, v)
U = quat(0, ux, uy, uz)
V = quat(0, vx, vy, vz)
return normalize(length(U*V)*Q - V*Q*U)
Чтобы найти кватернион с наименьшим вращением, которое вращается до , используйте
align(quat(1, 0, 0, 0), u, v)
Почему это обобщение?
это кватернион, ближайший к которому будет вращаться до . Что еще более важно ,Rэто кватернион, ближайший к которому локальное направление указывает в том же направлении, что и .
Это можно использовать, чтобы дать вам все возможные повороты, которые вращаются от до , в зависимости от выбора . Если вы хотите минимальное вращение отuкv, как дают другие решения, используйтеQ = quat(1, 0, 0, 0).
Чаще всего я нахожу, что реальная операция, которую вы хотите сделать, — это общее выравнивание одной оси по другой.
// If you find yourself often doing something like
quatFromTo(toWorldSpace(Q, localFrom), worldTo)*Q
// you should instead consider doing
align(Q, localFrom, worldTo)
Пример
Скажем, вам нужен кватернион, который представляет толькоQ's yaw, чистое вращение вокруг оси Y. Мы можем вычислитьYсо следующим.
Y = align(quat(Qw, Qx, Qy, Qz), vec(0, 1, 0), vec(0, 1, 0))
// simplifies to
Y = normalize(quat(Qw, 0, Qy, 0))
Выравнивание как проекционная матрица 4x4
Если вы хотите повторно выполнить одну и ту же операцию выравнивания, поскольку эта операция аналогична проекции кватерниона на двумерную плоскость, встроенную в четырехмерное пространство, мы можем представить эту операцию как умножение с матрицей проекции 4x4,A*Q.
I = mat4(
1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1)
A = I - leftQ(V)*rightQ(U)/length(U*V)
// which expands to
A = mat4(
1 + ux*vx + uy*vy + uz*vz, uy*vz - uz*vy, uz*vx - ux*vz, ux*vy - uy*vx,
uy*vz - uz*vy, 1 + ux*vx - uy*vy - uz*vz, uy*vx + ux*vy, uz*vx + ux*vz,
uz*vx - ux*vz, uy*vx + ux*vy, 1 - ux*vx + uy*vy - uz*vz, uz*vy + uy*vz,
ux*vy - uy*vx, uz*vx + ux*vz, uz*vy + uy*vz, 1 - ux*vx - uy*vy + uz*vz)
// A can be applied to Q with the usual matrix-vector multiplication
R = normalize(A*Q)
//LeftQ is a 4x4 matrix which represents the multiplication on the left
//RightQ is a 4x4 matrix which represents the multiplication on the Right
LeftQ(w, x, y, z) = mat4(
w, -x, -y, -z,
x, w, -z, y,
y, z, w, -x,
z, -y, x, w)
RightQ(w, x, y, z) = mat4(
w, -x, -y, -z,
x, w, z, -y,
y, -z, w, x,
z, y, -x, w)
С точки зрения алгоритма самое быстрое решение выглядит в псевдокоде
Quaternion shortest_arc(const vector3& v1, const vector3& v2 )
{
// input vectors NOT unit
Quaternion q( cross(v1, v2), dot(v1, v2) );
// reducing to half angle
q.w += q.magnitude(); // 4 multiplication instead of 6 and more numerical stable
// handling close to 180 degree case
//... code skipped
return q.normalized(); // normalize if you need UNIT quaternion
}
Убедитесь, что вам нужны единицы кватернионов (как правило, это требуется для интерполяции).
ПРИМЕЧАНИЕ: неединичные кватернионы могут использоваться с некоторыми операциями быстрее, чем единичные.
Некоторые из ответов, по-видимому, не учитывают вероятность того, что кросс-произведение может быть равно 0. Ниже приведен фрагмент представления по угловой оси:
//v1, v2 are assumed to be normalized
Vector3 axis = v1.cross(v2);
if (axis == Vector3::Zero())
axis = up();
else
axis = axis.normalized();
return toQuaternion(axis, ang);
toQuaternion может быть реализовано следующим образом:
static Quaternion toQuaternion(const Vector3& axis, float angle)
{
auto s = std::sin(angle / 2);
auto u = axis.normalized();
return Quaternion(std::cos(angle / 2), u.x() * s, u.y() * s, u.z() * s);
}
Если вы используете библиотеку Eigen, вы также можете просто сделать:
Quaternion::FromTwoVectors(from, to)
В соответствии с выводом вращения кватерниона между двумя углами можно повернуть вектор u на вектор v с
function fromVectors(u, v) {
d = dot(u, v)
w = cross(u, v)
return Quaternion(d + sqrt(d * d + dot(w, w)), w).normalize()
}
Если известно, что векторы от u до вектора v являются единичными векторами, функция сводится к
function fromUnitVectors(u, v) {
return Quaternion(1 + dot(u, v), cross(u, v)).normalize()
}
В зависимости от вашего варианта использования может потребоваться обработка случаев, когда скалярное произведение равно 1 (параллельные векторы) и -1 (векторы, указывающие в противоположных направлениях).
Работая только с нормализованными кватернионами, мы можем выразить ответ Джозефа Томпсона в следующих терминах.
Пусть q_v = (0, u_x, v_y, v_z) и q_w = (0, v_x, v_y, v_z) и рассмотрим
q = q_v * q_w = (-u точка v, uxv).
Таким образом, представляя q как q (q_0, q_1, q_2, q_3), мы имеем
q_r = (1 - q_0, q_1, q_2, q_3)