用于测量向量之间角度的廉价算法

Question

找到两个向量之间的角度并不难 使用余弦法则。但是，由于我正在为资源非常有限的平台进行编程，因此我希望避免诸如 sqrt 和 arccos 之类的计算。即使是简单的划分也应该尽可能地受到限制。

幸运的是，我不需要角度本身，而只需要一些与所述角度成比例的值。

因此，我正在寻找一些计算成本低的算法来计算与两个向量之间的角度相关的量。到目前为止，我还没有找到符合要求的东西，我自己也无法想出一些东西。

Answer 1

如果您不需要实际的欧几里德角度，而是可以用作角度比较基础的东西，那么更改为出租车几何形状可能是一个选择，因为您可以放弃三角学，并且它的速度很慢，同时保持准确性（或至少准确性确实有很小的损失，见下文）。

在主要的现代浏览器引擎中，加速系数在 1.44 - 15.2 之间，并且精度几乎与 atan2 中的相同。计算菱形角度平均比 atan2 快 5.01 倍，并且在 Firefox 18 中使用内联代码，加速达到系数 15.2。速度比较：http://jsperf.com/diamond-angle-vs-atan2/2< /a>.

代码非常简单：

function DiamondAngle(y, x)
{
    if (y >= 0)
        return (x >= 0 ? y/(x+y) : 1-x/(-x+y)); 
    else
        return (x < 0 ? 2-y/(-x-y) : 3+x/(x-y)); 
}

上面的代码给出了 0 到 4 之间的角度，而 atan2 给出了 -PI 和 PI 之间的角度，如下表所示：

请注意，菱形角度始终为正且在 0-4 范围内，而 atan2 也给出负弧度。所以钻石角度更加标准化。另一个注意事项是atan2给出了更精确的结果，因为范围长度是2*pi（即6.283185307179586），而在菱形角度中它是4。在实践中这不是很重要，例如。 rad 2.3000000000000001 和 2.3000000000000002 都采用菱形角度 1.4718731421442295，但如果我们通过删除一个零来降低精度，则 rad 2.300000000000001 和 2.300000000000002 会给出两个不同的菱形角度。菱形角度的这种“精度损失”非常小，只有当距离很大时才会产生重大影响。您可以在 http://jsbin.com/bewodonase/1/edit?output（旧版本：http://jsbin.com/idoyon/1）：

上面的代码足以进行快速角度比较，但在许多情况下需要将菱形角度转换为弧度，反之亦然。如果你例如。有一些弧度角度的公差，然后你有一个 100,000 次循环，将此公差与其他角度进行比较，使用 atan2 进行比较是不明智的。相反，在循环之前，您将弧度公差更改为出租车（菱形角度）公差，并使用菱形公差进行循环内比较，这样您就不必在代码的速度关键部分使用慢速三角函数（= in循环）。

进行此转换的代码如下：

function RadiansToDiamondAngle(rad)
{
  var P = {"x": Math.cos(rad), "y": Math.sin(rad) };
  return DiamondAngle(P.y, P.x);
}  

正如您所注意到的，有 cos 和 sin。如您所知，它们很慢，但您不必在循环内进行转换，但在循环之前进行转换，并且加速是巨大的。

如果由于某种原因，您必须将菱形角度转换为弧度，例如。循环并进行角度比较后返回例如。比较的最小角度或其他弧度，代码如下：

function DiamondAngleToRadians(dia)
{
  var P = DiamondAngleToPoint(dia);
  return Math.atan2(P.y,P.x);
}

function DiamondAngleToPoint(dia)
{
  return {"x": (dia < 2 ? 1-dia : dia-3), 
  "y": (dia < 3 ? ((dia > 1) ? 2-dia : dia) : dia-4)};
}

这里使用 atan2，它很慢，但想法是在任何循环之外使用它。您不能通过简单地乘以某个因子来将菱形角转换为弧度，而是在出租车几何形状中找到一个点，该点与正 X 轴之间的菱形角就是所讨论的菱形角，然后使用 atan2 将此点转换为弧度。

这对于快速角度比较来说应该足够了。

当然，还有其他 atan2 加速技术（例如 CORDIC 和查找表），但据我所知，它们都失去了准确性，并且仍然可能比 atan2 慢。

背景：我已经测试了几种技术：点积、内积、余弦定律、单位圆、查找表等，但在速度和准确性都很重要的情况下，没有什么是足够的。最后我在 http:// www.freesteel.co.uk/wpblog/2009/06/encoding-2d-angles-without-trigonometry/ 具有所需的功能和原理。

我首先假设出租车距离也可以用于准确和快速的距离比较，因为欧几里得中的距离越大，出租车中的距离也越大。我意识到与欧几里德距离相反，起点和终点之间的角度会影响出租车距离。只有垂直和水平向量的长度可以在欧几里德和出租车之间轻松快速地转换，但在所有其他情况下，您都必须考虑角度，然后过程太慢（？）。

因此，作为结论，我认为在速度关键的应用程序中，角度和/或距离的多次比较的循环或递归，在出租车空间中比较角度和在欧几里得（平方，不使用 sqrt）空间中的距离比较速度更快。

Answer 2

您是否尝试过 CORDIC 算法？它是一个仅用加/减/位移+表格来解决极坐标 ↔ 矩形问题的通用框架，本质上是通过 tan^-1 (2^-n) 形式的角度进行旋转。您可以通过更改迭代次数来权衡准确性和执行时间。

在您的情况下，将一个向量作为固定参考，并将另一个向量复制到临时向量，使用弦线角度将其旋转到第一个向量（大致二等分），直到达到所需的角度精度。

（编辑：使用点积的符号来确定每一步是否向前或向后旋转。尽管如果乘法足够便宜以允许使用点积，那么就不必担心 CORDIC，也许可以使用角度 π/2ⁿ 的旋转矩阵的正弦/余弦对表，以解决二分问题。）

（编辑：我喜欢 Eric Bainville 在评论中的建议：旋转两个向量都趋于零并跟踪角度差。）

Answer 3

回到几 K RAM 和数学能力有限的机器的时代，我使用查找表和线性插值。基本思想很简单：创建一个具有所需分辨率的数组（更多元素可以减少插值产生的误差）。然后在查找值之间进行插值。

这是一个处理中的示例 （原始死链接）。

您也可以使用其他三角函数来执行此操作。在 6502 处理器上，这使得计算完整 3D 线框图形的速度提高了一个数量级。

Answer 4

在此，我仍然没有发表评论的特权（尽管我在 math.se 上有评论），所以这实际上是对 Timo 关于钻石角度的帖子的回复。

基于 L1 范数的菱形角的整个概念是最有趣的，如果它只是比较哪个向量相对于正 X 轴具有更大/更小，那就足够了。然而，OP确实提到了两个通用向量之间的角度，我猜想OP想要将它与一些公差进行比较，以找到平滑度/角状态或类似的东西，但不幸的是，似乎只有jsperf.com上提供的公式或者 freesteel.co.uk （上面的链接）似乎不可能使用菱形角来做到这一点。

观察我的公式 Asymptote 实现的以下输出：

Vectors : 50,20 and -40,40
Angle diff found by acos      : 113.199
Diff of angles found by atan2 : 113.199
Diamond minus diamond         : 1.21429
Convert that to degrees       : 105.255
Rotate same vectors by 30 deg.
Vectors : 33.3013,42.3205 and -54.641,14.641
Angle diff found by acos      : 113.199
Diff of angles found by atan2 : 113.199
Diamond minus diamond         : 1.22904
Convert that to degrees       : 106.546
Rotate same vectors by 30 deg.
Vectors : 7.67949,53.3013 and -54.641,-14.641
Angle diff found by acos      : 113.199
Diff of angles found by atan2 : 113.199
Diamond minus diamond         : 1.33726
Convert that to degrees       : 116.971

因此，重点是您不能执行 Diamond(alpha)-diamond(beta) 并将其与某个容差进行比较，这与您可以对 atan2 的输出进行的操作不同。如果您想做的只是钻石（alpha）>钻石（beta）那么我认为钻石就可以了。

Answer 5

叉积与两个向量之间的角度成正比，当向量归一化且角度较小时，由于角度近似值较小，叉积非常接近以弧度表示的实际角度。

具体来说：

I1Q2-I2Q1 与 I1Q1 和 I2Q2 之间的角度成正比。

Answer 6

如果使用极坐标而不是笛卡尔坐标（或者“以及”使用笛卡尔坐标）来定义/存储向量，那么解决方案将是微不足道的。

Answer 7

两个向量 (x1, y1) 和 (x2, y2) 的点积等于

x1 * x2 + y1 * y2 

且 等于两个向量的长度乘以它们之间角度的余弦的乘积。

因此，如果您首先对两个向量进行归一化（将坐标除以长度），

Where length of V1 L1 = sqrt(x1^2 + y1^2),  
  and length of V2 L2 = sqrt(x2^2 + y2^2),

那么归一化向量是

(x1/L1, y1/L1),  and (x2/L2, y2/L2),  

归一化向量的点积（与向量之间角度的余弦相同）当然

 (x1*x2 + y1*y2)
 -----------------
     (L1*L2)

这可能就像计算上一样与计算余弦一样困难

Answer 8

如果您需要计算平方根，请考虑使用 invsqrt hack。

acos((x1*x2 + y1*y2) * invsqrt((x1*x1+y1*y1)*(x2*x2+y2*y2)));

Answer 9

点积可能适合您的情况。它与角度不成正比，而是“相关”。

Answer 10

为了更便宜地计算向量差角的度量，我制作了“全弧度”方程的简化版本，它仅使用 1 次除法和多次乘法（可以轻松地放入 SIMD 中以并行处理多对向量）：

也仅一个条件（仅符号检查） - 也可以将其放入带有符号位提取的 SIMD 中并用作掩码（或与零比较）。

没有使用 sqrt() 和 arccos()。

它将向量差值返回为固定范围 0..2（浮点）内的唯一正度量。其中 0 是完全相干矢量（零角度差），2 是完全反向（最大可能角度差）。

当前的 C 函数是：

float fDiffAngleVect(int x1, int y1, int x2, int y2)
{
  float fResult = 0.0f;
  // check if any of 2 input vectors is zero vector - return 0
  if ((x1 == 0) && (y1 == 0) || (x2 == 0) && (y2 == 0))
    return 0.0f;

  int iUpper = x1 * x2 + y1 * y2;

  if (iUpper > 0)
  {
    fResult = 1.0f - ((float)(iUpper * iUpper) / (float((x1 * x1 + y1 * y1) * (x2 * x2 + y2 * y2))));
  }
  else
  {
    fResult = 1.0f + ((float)(iUpper * iUpper) / (float((x1 * x1 + y1 * y1) * (x2 * x2 + y2 * y2))));
  }

  return fResult;

}

我的用例是要比较的仅整数 x,y 向量，以便您可以将其重写为全浮点输入和输出。