【问题标题】：Is there an iterative way to calculate radii along a scanline?是否有一种迭代方法来计算沿扫描线的半径？
【发布时间】：2009-09-15 21:40:44
【问题描述】：

我正在处理一系列 Y 值相同但 X 值不同的点。我通过将 X 加一来遍历这些点。例如，我可能有 Y = 50，X 是从 -30 到 30 的整数。我的算法的一部分涉及找到从每个点到原点的距离，然后进行进一步处理。

分析后，我发现距离计算中的 sqrt 调用占用了我大量的时间。有没有迭代的方法来计算距离？

换句话说：

我想有效地计算：r[n] = sqrt(x[n]*x[n] + y*y))。我可以保存上一次迭代的信息。每次迭代都会通过增加 x 来改变，所以 x[n] = x[n-1] + 1。我不能使用 sqrt 或 trig 函数，因为它们太慢了，除了在每个扫描线的开头。

我可以使用近似值，只要它们足够好（误差小于 0.l%）并且引入的错误是平滑的（我无法将其放入预先计算的近似值表中）。

附加信息： x 和 y 总是介于 -150 和 150 之间的整数

我明天将尝试几个想法，并根据最快的答案标记最佳答案。

结果

我做了一些计时

距离公式：16 ms/迭代
Pete 的交叉解决方案：8 毫秒/迭代
wrang-wrang 预计算解：8ms / 迭代

我希望考试能在两者之间做出决定，因为我喜欢这两个答案。我会选择 Pete's，因为它使用的内存更少。

【问题讨论】：

我的第一个问题是：你真的需要计算 r，而不是 r^2，因为 (x+1)^2 = x^2 +2x + 1，这对更新。
是的，我确实需要计算 r，而不是 r^2。为什么人们一直问我这个？
...因为这就是让它变得困难的原因:-)您实际上是使用所有这些 r 值，还是只使用其中一些？如果我站在你的立场上，那我就会集中精力。
... 8 毫秒，60 点 = 133us / 点。这对我来说听起来很重要。
每次迭代都会处理数百条这样的扫描线。

标签： optimization mathematical-optimization

【解决方案1】：

只是为了感受一下，对于您的范围 y = 50，x = 0 给出 r = 50 和 y = 50，x = +/- 30 给出 r ~= 58.3。您需要 +/- 0.1% 或 +/- 0.05 绝对值的近似值。这比大多数库 sqrt 的准确度要低很多。

两种近似方法 - 您可以根据前一个值的插值计算 r，或者使用合适序列的一些项。

从前一个 r 插值

r = ( x² + y² ) ^1/2

博士/dx = 1/2 。 2倍。 ( x² + y² ) ^-1/2 = x/r

    double r = 50;
    
    for ( int x = 0; x <= 30; ++x ) {
        
        double r_true = Math.sqrt ( 50*50 + x*x );
        
        System.out.printf ( "x: %d r_true: %f r_approx: %f error: %f%%\n", x, r, r_true, 100 * Math.abs ( r_true - r ) / r );
        
        r = r + ( x + 0.5 ) / r; 
    }

给予：

x: 0 r_true: 50.000000 r_approx: 50.000000 error: 0.000000%
x: 1 r_true: 50.010000 r_approx: 50.009999 error: 0.000002%
....
x: 29 r_true: 57.825065 r_approx: 57.801384 error: 0.040953%
x: 30 r_true: 58.335225 r_approx: 58.309519 error: 0.044065%

这似乎符合0.1%误差的要求，所以我没有费心编写下一个，因为它需要更多的计算步骤。

截断序列

x 接近于零的 sqrt ( 1 + x ) 的泰勒级数是

sqrt ( 1 + x ) = 1 + 1/2 x - 1/8 x² ... + ( - 1 / 2 )ⁿ⁺¹ x ⁿ

使用 r = y sqrt ( 1 + (x/y)² ) 那么您正在寻找一个术语 t = ( - 1 / 2 )ⁿ⁺¹ 0.36ⁿ 幅度小于 0.001，log (0.002) > n log (0.18) 或 n > 3.6，因此对 x^4 取项应该没问题。

【讨论】：

我认为他的示例数据集是作为整个范围的一小部分给出的（否则 60 sqrts 应该不是问题）。如果他的 X 比你给出的 30 高得多，错误很快就会变得不可接受。也许稍作修改的版本是每 30 或 40 X 步重新计算 r_true 以控制错误。
哦，将 dr/dx 表示为函数 f(x,r_n-1) 是一个聪明的主意。这绝对应该足够好，至少可以将 sqrt() 的数量减少到 24 步中的 1 步。
我有大量（50,000/秒）相当短（100-200）的扫描线，因此在每条扫描线的开头重新启动近似值会很好。我正在考虑做积极的一面和消极的一面，都从 0 开始，这样误差是对称的。
我发现“从以前的 r 插值”解决方案存在一个小问题。当 Y 为 0 时，扫描线最终将到达 r 为 0 的点，当我尝试为 x = 1 计算 r 时，我将得到除以 0。我希望我已经完成了正面/负面的事情我之前想过。这将消除这个问题。
为 r

【解决方案2】：

Y=10000
Y2=Y*Y
for x=0..Y2 do
  D[x]=sqrt(Y2+x*x)

norm(x,y)=
  if (y==0) x
  else if (x>y) norm(y,x) 
  else {
     s=Y/y
     D[round(x*s)]/s
  }

如果你的坐标是平滑的，那么这个想法可以通过线性插值来扩展。要获得更高的精度，请增加 Y。

这个想法是 s*(x,y) 位于 y=Y 线上，您已经为其预先计算了距离。得到距离，然后除以 s。

我假设你真的确实需要距离而不是它的平方。

您也许还可以找到一个通用的 sqrt 实现，该实现会牺牲一些准确性来提高速度，但我很难想象它会击败 FPU 可以做什么。

通过线性插值，我的意思是将D[round(x)]改为：

f=floor(x)
a=x-f
D[f]*(1-a)+D[f+1]*a

【讨论】：

【解决方案3】：

这并不能真正回答您的问题，但可能会有所帮助...

我要问的第一个问题是：

“我需要 sqrt 吗？”。
“如果没有，如何减少sqrts的数量？”
然后是你的：“我可以用聪明的计算替换剩余的 sqrts 吗？”

所以我要开始：

您是否需要精确的半径，或者半径平方是否可以接受？ sqrt 有快速的近似值，但对于您的规范可能不够准确。
能否使用镜像象限或八分之一处理图像？通过批量处理相同半径值的所有像素，可以将计算次数减少 8 倍。
你能预先计算半径值吗？您只需要一个大小为您正在处理的图像大小的四分之一（或八分之一）的表格，并且该表格只需预先计算一次，然后在多次运行算法时重复使用。

如此聪明的数学可能不是最快的解决方案。

【讨论】：

【解决方案4】：

嗯，总是在尝试优化你的 sqrt，我见过的最快的是旧的 carmack quake 3 sqrt：

http://betterexplained.com/articles/understanding-quakes-fast-inverse-square-root/

也就是说，由于 sqrt 是非线性的，因此您将无法沿直线进行简单的线性插值以获得结果。最好的办法是使用表查找，因为它可以让您快速访问数据。而且，由于您似乎是按整数进行迭代，因此表查找应该非常准确。

【讨论】：

Carmack 的 fastsqrt 不再优于硬件，我从最近的时间中发现。在蹩脚的 x87 FPU 时代，它确实很有用，但 SSE 引入了一个具有 1/1 吞吐量 (RCPSS) 的标量 recip 平方根估计操作码。
编译器不能保证使用它。无论如何，值得一试，如果速度不快，可能值得去 SSE。再说一次，它可能不会。
Carmack 技术的另一个问题是现代管道非常深，浮点寄存器和 int 寄存器之间没有直接连线，因此从 int 转换为 float 会导致很大的 load-hit-store 停顿。您可以通过将 /arch:sse2 指定给 MSVC 或使用内部函数来让编译器使用 SSE 操作。
使用 FPU x87 可能会更快，但操作没有说明他正在研究的架构，我在 8 位 PIC 上使用了这种方法来获得更快的 sqrt 计算并使用了这篇论文：lomont.org/Math/Papers/2003/InvSqrt.pdf（也链接在上面的文章中）为我的应用程序找到最好的初始猜测
作为记录，我使用的是双核和四核英特尔 CPU。

【解决方案5】：

好吧，您可以从 x=0 开始镜像（您只需要计算 n>=0，并将这些结果复制到相应的 n

如果这还不够准确，我可以指出这对于 SIMD 来说是一项非常好的工作，它将为您提供 SSE 和 VMX（以及着色器模型 2）上的倒数平方根运算。

【讨论】：

【解决方案6】：

这有点像HAKMEM item：

第 149 项（明斯基）：循环算法这是一种优雅的绘制方式点绘图显示上的圆圈：
NEW X = OLD X - epsilon * OLD Y
NEW Y = OLD Y + epsilon * NEW(!) X
这是一个非常圆的椭圆以原点为中心，大小为由初始点决定。 epsilon 确定角度循环点的速度，和稍微影响偏心度。如果 epsilon 是 2 的幂，那么我们不甚至需要乘法，更不用说平方根、正弦和余弦！这 “圆”将完全稳定因为积分很快就变成了定期。

圆算法是由当我试图保存一个错误时在显示黑客中注册！本·格利有一个惊人的显示黑客只使用大约六七条指令，以及这是一个伟大的奇迹。但它是基本上是面向线的。发生了对我来说，这将是令人兴奋的有曲线，我试图得到一个最小的曲线显示技巧说明。

【讨论】：