在GLSL中实现强壮的atan(y,x)函数以将XY坐标转换为角度

Question

在GLSL中（具体来说是我使用的3.00版本），有两个版本的atan()函数：atan(y_over_x)只能返回从-PI/2到PI/2之间的角度，而atan(y/x)可以考虑所有四个象限，因此角度范围涵盖了从-PI到PI的所有内容，就像C++中的atan2()一样。

我想使用第二个atan函数将XY坐标转换为角度。然而，在GLSL中，atan()除了无法处理x = 0的情况外，还不太稳定。特别是当x接近零时，除法可能会溢出，导致得到一个相反的角度（你得到的角度接近于-PI/2，而应该是大约PI/2）。

有什么好的、简单的实现方法可以在GLSL的atan(y,x)之上构建，使其更加健壮？

- HuaTham

---

在GLSL中，atan()函数除了无法处理x = 0的情况之外... -- 这是不正确的：atan(y,x)被设计用来处理x或y（但不能同时为零）的情况。当y为正时，atan(y,0.)返回PI/2，当y为负时，返回-PI/2。 - Jonathan Lidbeck

根据规范 https://www.opengl.org/sdk/docs/man/html/atan.xhtml，在 x = 0 时结果未定义。 - HuaTham

HuaTham：你说得对。哎呀。看起来至少它是常见的实现方式——我在运行测试模式的任何设备上都没有看到atan(y,0)的情况失败，但显然这不是可以依赖的行为。感谢您提供规范链接。 - Jonathan Lidbeck

@Jonathan，我也感到困惑，但事实证明该网页不准确。 事实上，完整规范（§8.1，p143）说明了您预期的情况 - x或y可以为0，但不能同时为0。 不确定该网页是如何出现这种情况的。 - AnOccasionalCashew

从您的完整规范中，atan(y,x)返回“其正切为y / x的角度”。按原样使用，可以看出其规范不能保证在x = 0时返回可用结果，因为y/0未定义。这仍然与Web版本一致。 - HuaTham

Answer 1

我将回答自己的问题，分享我的知识。我们首先注意到不稳定性出现在x接近于零的情况下。然而，我们也可以将其转换为abs(x) << abs(y)。因此，我们首先将平面（假设我们在单位圆上）分成两个区域：一个是|x| <= |y|，另一个是|x| > |y|，如下所示：

我们知道，在绿色区域内，atan(x,y)更加稳定——当x接近于零时，我们得到了与atan(0.0)非常稳定的结果，而通常的atan(y,x)在橙色区域更加稳定。您还可以亲自验证以下这个关系：

atan(x,y) = PI/2 - atan(y,x)

对于所有非原点(x,y)，该函数都成立，即使在其未定义的情况下，我们也谈论能够返回从-PI到PI整个范围内的角度值的atan(y,x)，而不是仅返回在-PI/2到PI/2之间的角度值的atan(y_over_x)。因此，我们健壮的GLSL的atan2()例程非常简单：

float atan2(in float y, in float x)
{
    bool s = (abs(x) > abs(y));
    return mix(PI/2.0 - atan(x,y), atan(y,x), s);
}

顺便提一下，数学函数 atan(x) 的恒等式实际上是：

atan(x) + atan(1/x) = sgn(x) * PI/2

这是正确的，因为它的范围是(-π/2, π/2)。

- HuaTham

---

被接受的答案实际上会加剧我的三星Galaxy S4上的问题，因为它的atan(Y,X)实现在Y的小值上失败，而不是X。 - Jonathan Lidbeck

Answer 2

根据您的目标平台，这可能是一个已解决的问题。OpenGL规范中的atan(y,x)指定它应该在所有象限中工作，只有当x和y都为0时行为未定义。

因此，人们会期望任何体面的实现在所有轴附近都是稳定的，因为这是二参数atan（或atan2）的整个目的。

提问者/回答者正确地指出了一些实现会采取捷径。然而，接受的解决方案假设一个不好的实现总是在x接近零时不稳定：在某些硬件上（例如我的Galaxy S4），当x接近零时该值是稳定的，但当y接近零时不稳定。

要测试您的GLSL渲染器对atan(y,x)的实现，这是一个WebGL测试图案。点击下面的链接，只要您的OpenGL实现体面，您应该会看到类似于这样的东西：

使用本地atan(y,x)的测试模式： http://glslsandbox.com/e#26563.2

如果一切正常，您应该看到8种不同的颜色（忽略中心）。

链接的演示样本对几个x和y的值进行了atan(y,x)采样，包括0、非常大和非常小的值。中央框是atan(0.,0.) - 在数学上未定义，并且实现因硬件而异。在我测试的硬件上，我看到了0（红色）、PI/2（绿色）和NaN（黑色）。

这是已接受解决方案的测试页面。注意：主机的WebGL版本缺少mix(float,float,bool)，因此我添加了一个与规范匹配的实现。

使用已接受答案的atan2(y,x)的测试模式： http://glslsandbox.com/e#26666.0

- Jonathan Lidbeck

---

根据GLSL规范，atan函数会找到y/x上的一个角度。很自然地，我们可以看出当x越接近于0时，y/x变得数值不稳定。因此我认为三星S4在y = 0附近发生不稳定意味着其实现与规范不兼容。例如，在第一象限中，它可能会在内部计算atan(x/y)后从pi/2中减去。我无法访问三星设备，但今天（2019年）在Mac / iPhone上两种实现的渲染结果看起来非常好。 - HuaTham

不完全正确。规范说明atan“返回一个角度，其正切为y/x”。这并不需要实际计算y/x。x/y同样有用，但是只有糟糕的实现才会依赖于其中之一（例如S4，它只使用x/y）。 - Jonathan Lidbeck

Answer 3

您提出的解决方案在 x=y=0 的情况下仍然无法解决。在这种情况下，atan() 函数都返回 NaN。

此外，我不会依赖于 mix 在两种情况之间切换。我不确定这是如何实现/编译的，但是 IEEE 浮点数规则对于 x*NaN 和 x+NaN 的结果再次为 NaN。因此，如果您的编译器确实使用了 mix/插值，则 x=0 或 y=0 的结果应该是 NaN。

以下是另一个修复方法，对我来说解决了这个问题：

float atan2(in float y, in float x)
{
    return x == 0.0 ? sign(y)*PI/2 : atan(y, x);
}

当 x=0 时，角度可以是 ±π/2。其中一个取决于仅有的 y。如果 y=0，则角度可以是任意值（向量长度为0）。在这种情况下，sign(y) 返回 0，这是可以接受的。

- Johannes Jendersie

---

数学上来说，在 (x, y) = (0, 0) 处是不确定的。你期望在原点处的角度是多少？它可以是任意值。另外，请仔细阅读我的帖子。我的 mix() 版本使用的是布尔插值器，而不是浮点数插值器。 - HuaTham

我并没有说角度在（x，y）=（0,0）处被定义。我只是将其设置为零作为有效选择。这里可以是任何角度，但NaN使我的程序膨胀了（NaN不是有效的角度，即使可以任意选择）我不知道有布尔混合操作（但你是对的）。 - Johannes Jendersie

Answer 4

有时候，提高代码性能的最佳方法是避免首次调用。例如，您可能想确定向量的角度是为了使用该角度构建旋转矩阵，其中组合了角度的正弦和余弦。然而，相对于原点的向量的正弦和余弦已经隐藏在向量本身中。您只需要通过将每个向量坐标除以向量的总长度来创建一个归一化版本的向量。以下是计算二维向量 [ x y ] 角度的正弦和余弦的示例：

double length = sqrt(x*x + y*y);
double cos = x / length;
double sin = y / length;

一旦您获取了正弦和余弦值，就可以直接使用这些值填充一个旋转矩阵，通过相同的角度顺时针或逆时针旋转任意向量，或者您可以连接第二个旋转矩阵以旋转到非零角度。在这种情况下，您可以将旋转矩阵视为将任意向量“标准化”为零角度。该方法也可拓展到三维（或N维）情形，不过例如在3D旋转中，需要计算三个角度和六组正弦/余弦（每平面一个角度）。

在能够使用该方法的情况下，您可以完全跳过atan的计算，从而大幅提高效率，因为您想要确定角度的唯一原因是为了计算正弦和余弦值。通过跳过角度空间的转换，您不仅避免了分母为零的问题，还改善了接近极点的角度的精度，否则这些角度会被大数乘除所影响。我曾经在一个GLSL程序中成功地使用这种方法将一个场景旋转到零度，以简化计算。

有时候，我们会陷入当前问题迷雾之中，而忘记为什么需要这些信息。虽然这种方法不适用于所有情况，但思考不同角度可能会帮助你解决问题。

- Steve Phillips

Answer 5

一个公式，可以为任何坐标值x和y给出四个象限中的角度。当x=y=0时，结果未定义。

f(x,y)=pi()-pi()/2*(1+sign(x))* (1-sign(y^2))-pi()/4*(2+sign(x))*sign(y)

   -sign(x*y)*atan((abs(x)-abs(y))/(abs(x)+abs(y))) 

- theodore panagos