计算两个向量之间顺时针角度的直接方法

115

我想找出两个向量(二维或三维)之间的顺时针角度。

使用点积的经典方法可以给出内角(0-180度),但我需要使用一些if语句来确定结果是我需要的角度还是其补角。

有没有直接计算顺时针角度的方法?


9
为什么不使用 std::atan2() - user529758
5
在三维向量中,如何定义“顺时针角度”? - Martin R
1
@Felics请看我的回答。在三维空间中,顺时针方向没有明确定义。这是一个平面术语。 - kassak
5
@Felics说:“在二维中,“顺时针”有明确的定义,但在三维中没有。检查叉积的z坐标(如Nickolay O.的答案所述)意味着在三维中:‘对于从x/y平面上方向下观察的观察者来说是顺时针方向’。” - Martin R
2
@Felics 另外,我应该指出,由于毛球定理http://en.wikipedia.org/wiki/Hairy_ball_theorem,您无法连续定义3D顺时针角度。您始终会有一对向量,其中一个的epsilon运动将导致时钟方向的即时切换,从而导致角度符号的变化。 - kassak
显示剩余4条评论
10个回答

277

2D情况

就像点积与角度的余弦成正比一样,行列式与角度的正弦成正比。因此,你可以通过以下方式计算角度:

dot = x1*x2 + y1*y2      # Dot product between [x1, y1] and [x2, y2]
det = x1*y2 - y1*x2      # Determinant
angle = atan2(det, dot)  # atan2(y, x) or atan2(sin, cos)

这个角度的方向与坐标系的方向相匹配。在一个左手坐标系中,即常见于计算机图形的x向右,y向下,这意味着顺时针角度将得到正号。如果坐标系的方向是数学上的,y向上,那么得到的将是逆时针角度,这是数学中的约定。改变输入的顺序将改变符号,所以如果对符号不满意,只需交换输入即可。

三维情况

在三维空间中,任意放置的两个向量定义了它们自己的旋转轴,该轴垂直于两者。该旋转轴没有固定的方向,这意味着你无法唯一确定旋转角度的方向。一种常见的约定是让角度始终为正,并将轴定向为适合正角度的方式。在这种情况下,归一化向量的点积足以计算角度。
dot = x1*x2 + y1*y2 + z1*z2    # Between [x1, y1, z1] and [x2, y2, z2]
lenSq1 = x1*x1 + y1*y1 + z1*z1
lenSq2 = x2*x2 + y2*y2 + z2*z2
angle = acos(dot/sqrt(lenSq1 * lenSq2))

注意,一些评论和备选答案建议不要使用acos,原因是在测量的角度很小时会出现数值问题。
嵌入在三维空间中的平面
一个特殊情况是当你的向量不是随意放置的,而是位于一个已知法向量n所确定的平面上。那么旋转轴也将在n的方向上,并且n的方向将确定该轴的方向。在这种情况下,你可以将2D计算中的n包含到determinant中,使其大小为3×3。
dot = x1*x2 + y1*y2 + z1*z2
det = x1*y2*zn + x2*yn*z1 + xn*y1*z2 - z1*y2*xn - z2*yn*x1 - zn*y1*x2
angle = atan2(det, dot)

这个工作的一个条件是法向量n的长度必须为单位长度。如果不是,你需要对其进行归一化。
作为三重积
这个行列式也可以表示为三重积,正如@Excrubulent在建议编辑中指出的那样。
det = n · (v1 × v2)

这可能在某些API中更容易实现,并对这里发生的事情提供了不同的视角:叉积与角度的正弦成正比,且垂直于平面,因此是n的倍数。点积基本上测量了该向量的长度,但附加了正确的符号。

范围0-360°

大多数 atan2 实现会返回一个范围为 [-π, π] 的弧度值,对应于 [-180°, 180°] 的角度值。如果你需要正角度的范围是 [0, 2π] 或 [0°, 360°],你可以将任何负值结果加上 2π。或者你可以避免使用条件判断,无条件地使用 atan2(-det, -dot) + π。如果你在一个罕见的情况下需要相反的修正,即 atan2 返回非负的 [0, 2π],而你需要带符号的角度范围为 [-π, π],则可以使用 atan2(-det, -dot) - π。这个技巧实际上不仅适用于这个问题,而且可以应用于大多数使用 atan2 的情况。记得检查你的 atan2 函数是以度还是弧度为单位,并根据需要进行转换。

4
给你点个赞 - 我懒得去弄清楚其他答案是否正确,你的回答最清晰易读,因此是帮助了我的答案。 - Excrubulent
2
对于2D,我得到了(0,180)和(-180,0)。可以检查结果是否为负,并添加360以获得一个漂亮的顺时针角度(例如,如果是-180,则添加360的结果为180,对于-90添加360的结果为270等)。不知道是我的计算还是Qt框架中qAtan2(y,x)的实现,但如果有人和我遇到同样的问题,这可能会有所帮助。 - rbaleksandar
17
atan2通常的范围是[-180°,180°]。为了获得[0°,360°]的值而不需要进行情况判断,可以将 atan2(y,x) 替换为 atan2(-y,-x) + 180° - MvG
3
不要对点积取反余弦!那在数学上是正确的,但在实际中极不准确。你可以用另一个 atan2(det,dot) 来替换你的三维方法;在这种情况下,det 将是叉积的长度。 - Don Hatch
2
@N4ppeL 如果想了解acos点积的不良行为,请参考以下链接(这个问题在两个不同的网站上都有提问,有不同的答案和参考资料):https://math.stackexchange.com/questions/1143354/numerically-stable-method-for-angle-between-3d-vectors/1782769 https://scicomp.stackexchange.com/questions/27689/numerically-stable-way-of-computing-angles-between-vectors - Don Hatch
显示剩余14条评论

6

这个答案与MvG's相同,但是解释方式不同(这是我努力理解为什么MvG的解决方案有效的结果)。

逆时针角度theta从给定法线n||n|| = 1)的视点到xy之间,由以下公式给出:

atan2( dot(n, cross(x,y)), dot(x,y) )

(1) = atan2( ||x|| ||y|| sin(theta),  ||x|| ||y|| cos(theta) )

(2) = atan2( sin(theta), cos(theta) )

(3) = x轴与向量(cos(theta), sin(theta))之间的逆时针角度

(4) = theta

其中||x||表示x的大小。

步骤(1)通过注意到以下内容得出:

cross(x,y) = ||x|| ||y|| sin(theta) n,
因此
dot(n, cross(x,y)) = dot(n, ||x|| ||y|| sin(theta) n) = ||x|| ||y|| sin(theta) dot(n, n)
如果 ||n|| = 1,那么等于 ||x|| ||y|| sin(theta)
步骤(2)是根据 atan2 的定义得出的,注意到 atan2(cy, cx) = atan2(y,x),其中 c 是一个标量。步骤(3)是根据 atan2 的定义得出的。步骤(4)是从 cos 和 sin 的几何定义得出的。

5

为了计算角度,您只需要在二维情况下调用atan2(v1.s_cross(v2), v1.dot(v2))。其中s_cross是叉积的标量类比(平行四边形的有向面积)。

在二维情况下,这将是楔形积。

在三维情况下,您需要定义顺时针旋转,因为从平面的一侧顺时针方向是一个方向,而从平面的另一侧则是另一个方向=)

这是逆时针角度,而顺时针角度则相反。


v1.cross(v2) 是一个向量,不是标量,不能像这样使用。Nickolay O. 在他的回答中描述了如何找到角度的“方向”。获取 2D 角度的一种方法是:angle = atan2f(v2.x, v2.y) - atan2f(v1.x, v1.y)。 - Mircea Ispas
2
@Felics 在二维叉乘中,通常意味着楔积。这是平行四边形的有符号面积。 对于二维情况,该公式绝对正确,因为它的点积等于 |v1||v2|*cos,而叉积等于 |v1||v2|sin。这就是为什么 atan2 可以在整个圆周范围内给出正确的角度。正如我之前所说,对于三维情况,你需要做一些假设来扩展顺时针方向的定义。 - kassak
1
@Felics:请注意 atan2f 的第一个参数是 y 坐标,因此应该为 angle = atan2f(v2.y, v2.x) - atan2f(v1.y, v1.x) - Martin R
1
@kassak:在二维情况下,您可以通过显式公式替换crossdot,这将消除有关cross返回3D向量的所有疑虑(但这只是一个建议,您可以忽略)。-否则,我喜欢这个解决方案,因为它只需要一个atan2f函数调用。 - Martin R
@Martin R,感谢您的好建议。我进行了一些更正以使公式的含义更清晰。 - kassak
你的公式不适用于向量(-1,0)和(0,1)的情况。 - user2083364

4

由于最简单和最优雅的解决方案之一隐藏在评论中,我认为将其作为单独的答案发布可能会有用。

acos 对于非常小的角度可能会导致不准确性,因此通常更喜欢使用 atan2。对于三维情况:

dot = x1*x2 + y1*y2 + z1*z2
cross_x = (y1*z2 – z1*y2)
cross_y = (z1*x2 – x1*z2) 
cross_z = (x1*y2 – y1*x2)
det = sqrt(cross_x*cross_x + cross_y*cross_y + cross_z*cross_z)
angle = atan2(det, dot)

2
链接到相应的答案或至少引用用户是一个不错的添加。 - N4ppeL
由于存在EN DASHes,这将导致编译错误。类似于"someFile.c:44:1: error: stray ‘\342’ in program. someFile.c:44:1: error: stray ‘\200’ in program. someFile.c:44:1: error: stray ‘\223’ in program"。在大多数文本编辑器的正则表达式模式中,可以通过\x{2013}进行搜索EN DASH(Unicode代码点U+2013)。 - Peter Mortensen
这三个数字是U+2013(十六进制0xE2 0x80 0x93)的UTF-8字节序列的八进制数。 - Peter Mortensen
这种错误的规范是 *Compilation error: stray ‘\302’ in program, etc*。 - Peter Mortensen
这是从哪里复制来的?某个网页上吗? - Peter Mortensen

2

对于二维方法,您可以使用余弦定理和“方向”方法。

计算线段P3:P1顺时针扫过线段P3:P2的角度。

    P1     P2
P3
    double d = direction(x3, y3, x2, y2, x1, y1);

    // c
    int d1d3 = distanceSqEucl(x1, y1, x3, y3);

    // b
    int d2d3 = distanceSqEucl(x2, y2, x3, y3);

    // a
    int d1d2 = distanceSqEucl(x1, y1, x2, y2);

    //cosine A = (b^2 + c^2 - a^2)/2bc
    double cosA = (d1d3 + d2d3 - d1d2)
        / (2 * Math.sqrt(d1d3 * d2d3));

    double angleA = Math.acos(cosA);

    if (d > 0) {
        angleA = 2.*Math.PI - angleA;
    }

这个和以上所提到的建议具有同样数量的超越操作,只多出一个或更少的浮点操作。

它使用的方法是:

 public int distanceSqEucl(int x1, int y1,
    int x2, int y2) {

    int diffX = x1 - x2;
    int diffY = y1 - y2;
    return (diffX * diffX + diffY * diffY);
}

public int direction(int x1, int y1, int x2, int y2,
    int x3, int y3) {

    int d = ((x2 - x1)*(y3 - y1)) - ((y2 - y1)*(x3 - x1));

    return d;
}

2

两个向量的数量积(点积)可以让你得到它们之间夹角的余弦

为了获得夹角的“方向”,你还应该计算叉积。它将让你检查(通过 z 坐标)夹角是否顺时针旋转,以及是否需要从 360 度中减去它。


1
即使这是正确的,我想避免的是计算某个值并确定计算出的值是否表示我的角度或我的角度补角。 - Mircea Ispas
我想知道这是否可能 :) 如果有(也许!)更好的方法,为什么要使用一些低效的方法来做事情。如果没有更好的方法,我会使用“标准”方法,但是询问更好的方法总是好的! - Mircea Ispas
实际上,标准方法并不总是高效的。 - Nickolay Olshevsky
@NickolayOlshevsky,您所说的“通过z坐标检查”是什么意思?我该如何去做呢? - Ogen
你应该检查z坐标的符号,就我记得的而言。 - Nickolay Olshevsky

0
对于二维情况,atan2 可以轻松计算一个向量与 (1, 0) 向量(即 x 轴)之间的夹角。
公式如下:
Atan2(y, x)

因此,您可以轻松地计算相对于x轴的两个角度之间的差异:

angle = -(atan2(y2, x2) - atan2(y1, x1))

为什么它不被用作默认解决方案?atan2 不够高效。顶部答案的解决方案更好。在 C# 上的测试表明,这种方法的性能比 atan2 低19.6%(100,000,000 次迭代)。

虽然不是关键问题,但仍然令人不愉快。

因此,其他有用的信息包括:

内外角度之间的最小角度:

abs(angle * 180 / PI)

角度制的完整角:

angle = angle * 180 / PI
angle = angle > 0 ? angle : 360 - angle

或者

angle = angle * 180 / PI
if (angle < 0)
    angle = 360 - angle;

0
如果“直接方式”是指避免使用if语句,那么我认为没有真正通用的解决方案。
然而,如果你的问题允许在角度离散化上失去一些精度,并且你可以接受在类型转换上失去一些时间,那么你可以将phi角度的[-pi,pi]允许范围映射到某个有符号整数类型的允许范围。然后你就可以免费获得补充性。然而,我实际上并没有在实践中使用这个技巧。很可能,浮点数到整数和整数到浮点数的转换开销会抵消任何直接性的好处。当这个角度计算被频繁执行时,最好将你的优先事项放在编写自动可向量化或并行化代码上。
此外,如果你的问题细节是这样的,即对于角度方向存在明确更可能的结果,那么你可以使用编译器内置函数将这些信息提供给编译器,以便它更高效地优化分支。例如,在GCC的情况下,那就是__builtin_expect函数。当你将它包装成这样的likelyunlikely宏(就像在Linux内核中一样)时,它会更加方便。
#define likely(x)      __builtin_expect(!!(x), 1)
#define unlikely(x)    __builtin_expect(!!(x), 0)

-1

计算二维情况下,两个向量 (xa,ya)(xb,yb) 之间的顺时针夹角公式。

Angle(vec.a-vec,b) =
  pi()/2*((1 + sign(ya))*
  (1 - sign(xa^2)) - (1 + sign(yb))*
  (1 - sign(xb^2))) + pi()/4*
  ((2 + sign(ya))*sign(xa) - (2 + sign(yb))*
  sign(xb)) + sign(xa*ya)*
  atan((abs(ya) - abs(xa))/(abs(ya) + abs(xa))) - sign(xb*yb)*
  atan((abs(yb) - abs(xb))/(abs(yb) + abs(xb)))

1
这是什么编程语言? - gdbdable
为什么 "vec.a-vec,b" 看起来非常不对称?你是指 "vec.a-vec.b" 吗? - Peter Mortensen
好的,OP已经“离开了大厦”: “上次出现超过3年前”。 - Peter Mortensen

-4

只需复制并粘贴此内容:

angle = (acos((v1.x * v2.x + v1.y * v2.y)/((sqrt(v1.x*v1.x + v1.y*v1.y) * sqrt(v2.x*v2.x + v2.y*v2.y))))/pi*180);

7
尽管这段代码可能回答了问题,但是加入解释“为什么”和“如何”帮助解决问题可以提高您的回答质量和持久性,特别是对于像这样的老问题。请参阅《如何撰写好答案》(https://stackoverflow.com/help/how-to-answer)。 - slothiful

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