logo标签列表

寻找表示从一个向量旋转到另一个向量的四元数

mathvectorquaternions
132

我有两个向量u和v。是否有一种方法可以找到表示从u旋转到v的四元数？

请查看https://laustep.github.io/stlahblog/posts/ReorientTransformation2.html。 - Stéphane Laurent
11个回答
160
Quaternion q;
vector a = crossproduct(v1, v2);
q.xyz = a;
q.w = sqrt((v1.Length ^ 2) * (v2.Length ^ 2)) + dotproduct(v1, v2);

不要忘记对q进行归一化。

Richard关于不存在唯一旋转的说法是正确的，但以上方法应该会给出“最短弧”，这可能是你所需要的。

41
请注意，这种方法无法处理平行向量的情况（既指向同一方向，又指向相反方向的情况）。在这些情况下，crossproduct函数将不会返回有效的结果，因此您需要首先检查分别是否满足dot(v1, v2) > 0.999999和 dot(v1, v2) < -0.999999的条件，如果是，则应分别返回标识四元数或绕任意轴的180度旋转对于平行向量和相反向量。 - sinisterchipmunk
6
实际上，如果 v1 = v2，则叉积将为 (0,0,0)，并且 w 将为正数，这将归一化为单位。根据 http://www.gamedev.net/topic/429507-finding-the-quaternion-betwee-two-vectors/page__p__3856228#entry3856228 ，它也应该适用于 v1 = -v2 及其附近。 - jpa
4
有人如何使这项技术工作？首先，sqrt((v1.Length ^ 2) * (v2.Length ^ 2))可以简化为v1.Length * v2.Length。我无法通过任何这种变化来产生有意义的结果。 - Joseph Thomson
5
是的，这个方法可行。请查看源代码。第61行处理向量面向相反方向的情况（返回PI），否则按@jpa的评论返回identity。第67行处理平行向量：在数学上是不必要的，但更快。第72行是Polaris878的答案，假设两个向量都是单位长度（避免了sqrt运算）。还可以查看单元测试。 - sinisterchipmunk
4
通常获取v.Length^2比获取v.Length本身要快，这就是我认为它被表述成这样的原因。 - lvella
显示剩余8条评论
87

半程向量解法

我想出了一个解决方案，我相信这正是Imbrondir试图呈现的（尽管有一个小错误，这可能是为什么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，因此我们必须沿着与 u（或v）正交的任意向量旋转 180 度作为我们的特殊情况解。这可以通过使用 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);
}

四元数的中点解法

这实际上是被接受回答中提出的解决方案，与半程向量的解法相比似乎稍微快一些（按我的测量结果约快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)));
}
14
+1：太好了！这个方法很有效。应该被选为最佳答案。 - Rekin
1
四元数语法在某些示例中被切换（Quaternion（xyz，w）和Quaternion（w，xyz））。此外，在最后一个代码块中似乎混合了弧度和度来表示角度（180 vs. k_cos_theta + k）。 - Guillermo Blasco
1
Quaternion(float, Vector3) 是通过标量-向量构造的，而 Quaternion(Vector3, float) 则是通过轴角构造的。这可能会让人感到混淆，但我认为它是正确的。如果您仍然认为它不对，请纠正我！ - Joseph Thomson
2
成功了！谢谢！不过，我找到了另一个类似的，并且解释得很好的链接（http://lolengine.net/blog/2014/02/24/quaternion-from-two-vectors-final）来执行上述操作。我想分享一下记录 ;) - sinner
1
@JosephThomson 半路四元数解决方案似乎来自于这里：http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/transforms/halfAngle.htm。 - legends2k
显示剩余6条评论
10

问题陈述不够明确：对于给定的向量对，不存在唯一的旋转。例如，考虑当u = <1, 0, 0>和v = <0, 1, 0>时。从u到v的一个旋转是绕z轴旋转π / 2。另一个从u到v的旋转是绕向量<1, 1, 0>旋转π。

3
实际上，可能有无穷多个可能的答案，因为在你将“from”向量与“to”向量对齐之后，仍然可以自由地绕其轴旋转结果。您知道通常可以用什么额外信息来约束此选择并使问题定义良好吗？ - Doug McClean
6

我对四元数不是很擅长。然而，我为此苦苦挣扎了数小时，也无法使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)（没有旋转）或者Quaternion(0，*v1) （180度旋转）。
我已经有一个可行的实现，但是你的更漂亮，所以我真的希望它能够工作。不幸的是，它未通过我的所有测试用例。我的所有测试看起来都像这样：quat = diffVectors(v1, v2); assert quat * v1 == v2。 - sinisterchipmunk
由于“角度”取决于点积，因此这很不可能起作用。 - sam hocevar
1
Quaternion() 函数在哪里？ - June Wang
我没有尝试过这个，但是看起来你可能只需要删除 v.normalize()。因此答案的标量部分将是 v.dot(v2) = (v1+v2).dot(v2) = 1 + v1.dot(v2)，而向量部分将是 v.cross(v2) = (v1+v2).cross(v2) = v1.cross(v2)。 - Don Hatch
4

要找到将u旋转到v的最小旋转四元数，请使用以下方法：

align(quat(1, 0, 0, 0), u, v)

通用解决方案

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)

为什么要做这种泛化？

R = align(Q, u, v)

最重要的是，R 是最接近 Q 的四元数，其本地 u 方向在全局指向 v 方向。因此，R 是最接近 Q 的四元数，它将旋转 u 到 v。

这可以用来给出所有可能的旋转，这些旋转从 u 到 v，取决于选择的 Q。如果您想要从 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)

例子

假如您想要得到四元数 Y，该四元数仅代表了绕 y 轴的旋转（也就是 Q 的偏航角），我们可以使用以下方法计算出 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投影矩阵作为对齐

如果您想要重复执行相同的对齐操作，因为此操作与将四元数投影到嵌入4D空间中的2D平面相同，我们可以将此操作表示为与4x4投影矩阵A*Q相乘。

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)

I = mat4(
    1, 0, 0, 0,
    0, 1, 0, 0,
    0, 0, 1, 0,
    0, 0, 0, 1)
1
这真是太神奇了！只是有点可悲的是我在十年前没有学到这个，如果那时候学会了就可以帮助很多IMU传感器数据处理了。有趣的是Madgwick过滤器使用梯度下降优化步骤解决相同的问题，但现在我知道你可以在一次对准中解决这个问题。真的非常有帮助！ - xaedes
4
为什么不使用纯四元数来表示向量？或许在归一化之后会更好。
q1 = (0 ux uy uz)'
q2 = (0 vx vy vz)'
q1 qrot = q2
用q1-1作为前置因子：
qrot = q1-1 q2
其中，q1-1 = q1conj / qnorm
这可以被看作是“左除”。 而“右除”则不是你想要的：
qrot,right = q2-1 q1
4
我迷失了，从q1旋转到q2的计算不是为q_2 = q_rot q_1 q_rot^-1 吗？ - yota
你是对的。我已经尝试过了，它不起作用。 - Paweł Iwaneczko
3

从算法角度来看，最快的解决方案在伪代码中。

 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
 }

请确认您需要使用单位四元数（通常，这是插值所需的）。

注意： 某些操作可以比单位四元数更快地使用非单位四元数。

2

有些答案似乎没有考虑到叉积可能为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)
toQuaternion(axis, ang) -> you forgot to specify what is ang - Maksym Ganenko
第二个参数是“角度”，它是四元数的轴角表示法的一部分，以弧度为单位测量。 - Shital Shah
你被要求获取四元数以从一个向量旋转到另一个向量。你没有角度，你必须先计算它。你的答案应该包含角度的计算。干杯！ - Maksym Ganenko
这是 C++ 吗？u.x() 是什么？ - June Wang
是的，这是C++。如果您正在使用Eigen库，则u是向量类型。 - Shital Shah
2
根据两个角度之间四元数旋转的推导，可以使用以下方法将向量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（指向相反方向的向量）的情况。
老实说，检查这种情况的最好方法就是在归一化之前检查幅度。只要不为0，就应该没问题。计算幅度的方式需要平方值。因此，任何需要比机器精度更接近一半的东西都具有0的幅度。对于浮点数，10^-20 * 10^-20 -> 0。 - Trey Reynolds
0

仅使用归一化四元数，我们可以用以下术语表达Joseph Thompson的答案。

设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）。normalize（）

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