在不实际计算该角度由角度矢量进行排序最快的方式(Fastest way to sort vector

2019-09-03 02:13发布

站内文章 / 后端开发

77 0

女 | 书童

私信

许多算法（例如格雷厄姆扫描）规定点或向量可以通过角度进行排序（或许作为从其他点看到的，即，使用差矢量）。这个顺序是固有的循环，并在这个循环被打破计算线性值往往不管那么多了。但真正的角度值确实不太重要，只要循环顺序保持不变。所以，做一个atan2呼吁每一个点可能是一种浪费。什么更快的方法有计算的值是严格单调的角度，方式atan2是什么？这些功能显然已经被一些所谓的“pseudoangle”。

Answer 1:

我开始玩这个，实现该规范是一种不完整的。 atan2具有不连续性，因为作为dx和dy是多种多样的，有一个点atan2将-pi和+ PI之间跳转。下图显示了通过@MvG建议的两个公式，而事实上它们都具有在相对于不同的地方的不连续性atan2 。（NB：我添加3到第一公式和4替代，使得线不图表上的重叠）。如果我添加atan2到图形那么这将是直线Y = X。所以，在我看来，可能有不同的答案，不同的地方人愿意把间断。如果一个人真的想复制atan2 ，答案（在这一流派）会

# Input:  dx, dy: coordinates of a (difference) vector.
# Output: a number from the range [-2 .. 2] which is monotonic
#         in the angle this vector makes against the x axis.
#         and with the same discontinuity as atan2
def pseudoangle(dx, dy):
    p = dx/(abs(dx)+abs(dy)) # -1 .. 1 increasing with x
    if dy < 0: return p - 1  # -2 .. 0 increasing with x
    else:      return 1 - p  #  0 .. 2 decreasing with x

这意味着，如果你正在使用的语言有一个标志功能，你可以避开返回符号（DY）（1-P），其将0-5的答案在中断返回之间的作用分支-2和+2。而同样的伎俩将与@ MVG原来的工作方法，人们可以返回符号（DX）（P-1）。

更新在下方留言，@MvG表明一行C实现这一点，即

pseudoangle = copysign(1. - dx/(fabs(dx)+fabs(dy)),dy)

@MvG说，它工作得很好，它看起来对我好:-)。

Answer 2:

我知道，一个可能的这种功能，我将在这里描述。

# Input:  dx, dy: coordinates of a (difference) vector.
# Output: a number from the range [-1 .. 3] (or [0 .. 4] with the comment enabled)
#         which is monotonic in the angle this vector makes against the x axis.
def pseudoangle(dx, dy):
    ax = abs(dx)
    ay = abs(dy)
    p = dy/(ax+ay)
    if dx < 0: p = 2 - p
    # elif dy < 0: p = 4 + p
    return p

那么，为什么这项工作？有一点要注意的是，缩放所有输入的长度不会影响输出中。因此，向量的长度(dx, dy)是无关紧要的，只有其方向的事项。集中在第一象限，我们可以暂时假设dx == 1 。然后dy/(1+dy)单调从零增长dy == 0到一个无限dy （即，对于dx == 0 ）。现在，其他象限必须和处理。如果dy是负的，那么初始p 。因此，对于正dx我们已经有一个范围-1 <= p <= 1在角度单调的。对于dx < 0我们改变符号并添加两个。给出的范围1 <= p <= 3为dx < 0以及一系列的-1 <= p <= 3整体上。如果负数是由于某种原因不希望的，所述elif注释行可以被包括在内，这将在第四象限从移位-1…0至3…4 。

我不知道，如果上述功能有一个既定的名字，谁可能第一次出版它。我前一段得到它，并从一个项目复制到下一个。然而，我发现发生在网络上的这一点，所以我会考虑这个剪断公众足够的重复使用。

有一种方法，以获得在范围[0 ... 4]（真正角度[0 ...2π]）而不引入另外的情况下区分：

# Input:  dx, dy: coordinates of a (difference) vector.
# Output: a number from the range [0 .. 4] which is monotonic
#         in the angle this vector makes against the x axis.
def pseudoangle(dx, dy):
    p = dx/(abs(dx)+abs(dy)) # -1 .. 1 increasing with x
    if dy < 0: return 3 + p  #  2 .. 4 increasing with x
    else:      return 1 - p  #  0 .. 2 decreasing with x

Answer 3:

我有点像三角，所以我知道的角度映射到一定值，我们通常的最佳途径是切线。当然，如果我们希望以不具有比较麻烦的有限数量{符号（X），Y / X}，它变得更混乱一点。

但是存在一个映射函数[1，+ INF [到[1,0 [称为逆，这将使我们能够具有有限的范围，这是我们将映射角度。正切的倒数是众所周知的余切，从而X / Y（是的，它的那样简单）。

一点说明，示出了切线和余切单位圆上的值：

您看到的值相同，当| X | = | Y |，你也看到，如果我们的颜色，输出之间[-1,1]上的两个圆，我们设法颜色的整圆的值的部分。有价值观的这种映射是连续的，单调的，我们可以做两本：

使用余切的相对具有相同的单调切线
添加2〜-cotan，以具有值一致，其中黄褐色= 1
添加4到圆的一半（比方说，下面在x = -y对角线）具有值装配在不连续的所述一个。

这给出了下面的分段函数，它是角度的连续和单调函数，只有一个不连续性（即最小）：

double pseudoangle(double dx, double dy) 
{
    // 1 for above, 0 for below the diagonal/anti-diagonal
    int diag = dx > dy; 
    int adiag = dx > -dy;

    double r = !adiag ? 4 : 0;

    if (dy == 0)
        return r;

    if (diag ^ adiag)
        r += 2 - dx / dy; 
    else
        r += dy / dx; 

    return r;
}

请注意，这是非常接近福勒的角度，具有相同的属性。形式上， pseudoangle(dx,dy) + 1 % 8 == Fowler(dx,dy)

谈性能，它比福勒的代码枝少得多（通常较不复杂的IMO）。编译-O3上的gcc 6.1.1，上述功能产生具有4个分支，其中两个都来自汇编代码dy == 0 （一个检查所述两个操作数都是“无序”，因此，如果DY是NaN ，并其他检查，如果他们是平等的）。

我认为这个版本比别人更准确，因为它只是使用尾数保留操作，直到结果向右移动的时间间隔。这应该是特别明显，当| X | << | Y | 或| Y | >> | X |，则操作| X | + | Y | 失去颇有些精度。

正如你在图上看到的角度pseudoangle关系也很好地接近线性的。

寻找其中的分支从何而来，我们可以做出以下回应：

我的代码不依赖于abs ，也没有copysign ，这使得它看起来更加自足。然而，随着符号位玩浮点值实际上是非常不重要的，因为它只是翻转一个单独的位（无分支！），所以这是比较吃亏的。
此外其他的解决方案建议在这里不检查是否abs(dx) + abs(dy) == 0将它划分之前，但是这个版本会失败，只要只有一个组件（DY）为0 -所以，在一个分支抛出（或在我的情况2）。

如果我们选择得到大致相同的结果（最多舍入误差），但没有分支机构，我们可能滥用copsign写：

double pseudoangle(double dx, double dy) 
{
    double s = dx + dy; 
    double d = dx - dy; 
    double r = 2 * (1.0 - copysign(1.0, s)); 
    double xor_sign = copysign(1.0, d) * copysign(1.0, s); 

    r += (1.0 - xor_sign);
    r += (s - xor_sign * d) / (d + xor_sign * s);

    return r;
}

如果dx和dy是绝对值接近更大的错误可能发生比以前的实现，由于消除在d或秒。有没有检查是否被零除要与提出的其他实现可比性，因为这既时dx和dy为0只发生。

Answer 4:

如果你能在排序时喂原有载体的角度，而不是进入一个比较功能，你可以把它一起工作：

只是一个单一的分支。
只有浮点比较和乘法。

避免加减使得数字更强大。双实际上可以随时精确地表示两个浮点数的产品，但不一定是他们的总和。这意味着单精度输入你能保证毫不费力地完美无瑕的结果。

这基本上是CIMBALI的解决方案重复两个向量，用树枝消除分歧乘了。它返回一个整数，符号匹配的比较结果（正，负或零）：

signed int compare(double x1, double y1, double x2, double y2) {
    unsigned int d1 = x1 > y1;
    unsigned int d2 = x2 > y2;
    unsigned int a1 = x1 > -y1;
    unsigned int a2 = x2 > -y2;

    // Quotients of both angles.
    unsigned int qa = d1 * 2 + a1;
    unsigned int qb = d2 * 2 + a2;

    if(qa != qb) return((0x6c >> qa * 2 & 6) - (0x6c >> qb * 2 & 6));

    d1 ^= a1;

    double p = x1 * y2;
    double q = x2 * y1;

    // Numerator of each remainder, multiplied by denominator of the other.
    double na = q * (1 - d1) - p * d1;
    double nb = p * (1 - d1) - q * d1;

    // Return signum(na - nb)
    return((na > nb) - (na < nb));
}

Answer 5:

不错..这里是返回-Pi，皮像许多arctan2功能的变体光盘。

编辑注：改变了我的pseudoscode适当的蟒蛇.. ARG变化，以与蟒蛇数学模块ATAN2兼容性（）。 EDIT2打扰更多代码来捕获的情况下DX = 0。

def pseudoangle( dy , dx ):
  """ returns approximation to math.atan2(dy,dx)*2/pi"""
  if dx == 0 :
      s = cmp(dy,0)
  else::
      s = cmp(dx*dy,0)  # cmp == "sign" in many other languages.
  if s == 0 : return 0 # doesnt hurt performance much.but can omit if 0,0 never happens
  p = dy/(dx+s*dy)
  if dx < 0: return p-2*s
  return  p

在这种形式下，最大误差仅为0.07〜为所有角度的弧度。（当然离开了π/ 2，如果你不关心的大小。）

现在的坏消息 - 使用python math.atan2我的系统上约25％的速度显然更换一个简单的解释代码犯规打编译intrisic。