最准确的线交点坐标计算与浮点数？答案

【问题标题】：Most accurate line intersection ordinate computation with floats?最准确的线交点坐标计算与浮点数？
【发布时间】：2009-05-22 07:58:22
【问题描述】：

我正在计算给定横坐标 x 的线上的点的纵坐标 y。该线由其两个端点坐标 (x0,y0)(x1,y1) 定义。端点坐标是浮点数，必须以浮点精度进行计算才能在 GPU 中使用。

数学以及因此幼稚的实现都是微不足道的。

令 t = (x - x0)/(x1 - x0)，然后 y = (1 - t) * y0 + t * y1 = y0 + t * (y1 - y0)。

问题在于 x1 - x0 很小的时候。结果将引入取消错误。当与 x - x0 之一结合时，在除法中，我预计 t 会出现重大错误。

问题是是否存在另一种更准确地确定 y 的方法？

即我应该先计算 (x - x0)*(y1 - y0)，然后除以 (x1 - x0) 吗？

y1 - y0 的差异总是很大。

【问题讨论】：

你能使用三角函数比如 sin 或 cos 吗？ AFAIK，所有较新的 GPU 都将这些作为一条指令，所以它应该很快。
如果 x 在 x0 和 x1 之间，则没有大错误。在计算 t 时，您处理相同数量级的值：x-x0、x1-x0。

标签： c++ c optimization math numerical

【解决方案1】：

在很大程度上，您的根本问题是根本性的。当 (x1-x0) 较小时，表示 x1 和 x0 的尾数只有几位不同。并且通过扩展，在 x0 和 x1 之间只有有限数量的浮点数。例如。如果仅尾数的低4位不同，则它们之间最多有14个值。

在您的最佳算法中，t 术语代表这些低位。继续举例，如果 x0 和 x1 相差 4 位，那么 t 也只能取 16 个值。这些可能值的计算是相当稳健的。无论您是计算 3E0/14E0 还是 3E-12/14E-12，结果都将接近 3/14 的数学值。

您的公式具有 y0

（我假设您对浮点表示足够了解，因此“(x1-x0) 很小”实际上意味着“相对于 x1 和 x0 本身的值而言很小”。1E-1 的差异是x0=1E3 时小，x0=1E-6 时大）

【讨论】：

我完全同意你的分析。 x0 和 x1 将仅相差几个不太重要的位，并且 x 将是这个小浮点可表示范围内的一个值。我们可以从中得出什么结论？没有办法比使用简单的数学实现获得更好的结果精度吗？这可以解释观察到的结果。
是的，这就是我的结论。如果您对双精度数进行中间计算，您为 y 找到的值不会显着提高 - 当您将结果四舍五入为浮点数时，您会达到相同的限制。

【解决方案2】：

你可以看看 Qt 的“QLine”（如果我没记错的话）来源；他们已经实现了取自“Graphics Gems”一书中的交集确定算法（参考必须在代码 cmets 中，这本书几年前在 EDonkey 上），这反过来又对适用性有一些保证以给定的位宽执行计算时给定的屏幕分辨率（如果我没记错的话，他们使用定点算术）。

【讨论】：

【解决方案3】：

如果您有可能这样做，您可以在计算中引入两种情况，具体取决于 abs(x1-x0)

编辑。另一种可能性是使用二分搜索的变体逐位获得结果。这会比较慢，但在极端情况下可能会改善结果。

// Input is X
xmin = min(x0,x1);
xmax = max(x0,x1);
ymin = min(y0,y1);
ymax = max(y0,y1);
for (int i=0;i<20;i++) // get 20 bits in result
{
  xmid = (xmin+xmax)*0.5;
  ymid = (ymin+ymax)*0.5;
  if ( x < xmid ) { xmax = xmid; ymax = ymid; } // first half
  else { xmin = xmid; ymin = ymid; } // second half
}
// Output is some value in [ymin,ymax]
Y = ymin;

【讨论】：

很遗憾，这是不可能的。我必须为给定的 x 值找到 y。
如果直线平行于 y 轴，这甚至是不可能的。看我的帖子。
这是 Cohen-Sutherland 裁剪算法中描述的方法。它很聪明，但正如您所说，效率不高。我想把它转换成整数计算，这样我们就知道我们有多少准确的位了。

【解决方案4】：

我已经实现了一个基准程序来比较不同表达式的效果。

我使用双精度计算 y，然后使用具有不同表达式的单精度计算 y。

这里是测试的表达式：

inline double getYDbl( double x, double x0, double y0, double x1, double y1 )
{
    double const t = (x - x0)/(x1 - x0);
    return y0 + t*(y1 - y0);
} 

inline float getYFlt1( float x, float x0, float y0, float x1, float y1 )
{
    double const t = (x - x0)/(x1 - x0);
    return y0 + t*(y1 - y0);
} 

inline float getYFlt2( float x, float x0, float y0, float x1, float y1 )
{
    double const t = (x - x0)*(y1 - y0);
    return y0 + t/(x1 - x0);
} 

inline float getYFlt3( float x, float x0, float y0, float x1, float y1 )
{
    double const t = (y1 - y0)/(x1 - x0);
    return y0 + t*(x - x0);
} 

inline float getYFlt4( float x, float x0, float y0, float x1, float y1 )
{
    double const t = (x1 - x0)/(y1 - y0);
    return y0 + (x - x0)/t;
}

我计算了双精度结果和单精度结果之间差异的平均值和标准差。

结果是平均没有超过 1000 个和 10K 个随机值集。我使用了带优化和不带优化的 icc 编译器以及 g++。

请注意，我必须使用 isnan() 函数来过滤掉虚假值。我怀疑这些是由差异或除法中的下溢引起的。

我不知道编译器是否重新排列了表达式。

无论如何，这个测试的结论是，上述重新排列的表达式对计算精度没有影响。误差保持不变（平均而言）。

【讨论】：

为了更好的比较，将所有输入和返回值作为浮点数；向上/向下适当加倍。您的问题是有限的精度是否是由于使用浮点数进行的计算。我和其他人认为，您看到的是输入本身的准确性有限，而不是计算的准确性。
另外，您使用了哪些优化设置？优化器可能会重新排序表达式评估顺序。
抱歉，无法跟踪所有 cmets。糟糕，查看我在此处复制的代码，我发现我将所有 t 声明为双精度值，这是一个错误。它应该全部为浮点数或全部为双精度。
将所有 double const t 更改为 float const t 后，计算超过 10K 值的平均值时没有区别。对于 icc，优化设置为 icc -xS -O3。对于 g++，我没有尝试过。没有区别。

【解决方案5】：

检查 x0 和 x1 之间的距离是否很小，即 fabs(x1 - x0) 取决于 x。您有无数个 y 值，因此您必须区别对待这种情况。

【讨论】：

你可以假设条件 x0
我想你不明白我的意思。如果线平行于 y 轴，“x0
看问题：dy总是很大，所以线永远不会平行于x轴。事实上，dx
我从来没有说过它会是 :) 我说平行于 y 轴。

【解决方案6】：

如果您的源数据已经是浮点数，那么您已经存在根本的不准确性。

为了进一步解释，想象一下如果您以图形方式执行此操作。您有一张 2D 方格纸，并标记了 2 个点。

案例 1：这些点非常准确，并且已经用非常锋利的铅笔做了标记。很容易画出连接它们的线，然后很容易在给定 x 的情况下得到 y（反之亦然）。

案例 2：这些点已经用一个大的粗毡尖笔做了标记，就像宾果记号笔一样。显然，您绘制的线将不太准确。你穿过斑点的中心吗？上边缘？底边？一个的顶部，另一个的底部？显然有许多不同的选择。如果两个点彼此靠近，那么变化会更大。

由于它们表示数字的方式，浮点数固有一定程度的不准确性，因此它们更多地对应于案例 2 而不是案例 1（有人可能认为这相当于使用任意精度库）。世界上没有任何算法可以弥补这一点。输入不精确的数据，输出不精确的数据

【讨论】：

您对浮点精度的看法是正确的，但最后一句话是错误的。见docs.sun.com/source/806-3568/ncg_goldberg.html。
糟糕。再次阅读后，您解决了表示精度的问题，而我解决了计算中的精度问题。在计算均值或方差时，即使在数学上是正确的，Knuth 算法也能提供比简单实现更好的准确度。
有时这种基本的不准确性会产生更大的影响。例如，当 x 为 4.002 与 4.001 与 4.0 以及当 x 为 5.002 与 5.001 与 5.0 时，考虑 5/(5-x) 的结果差异。越接近 5，误差影响越大。这里，“问题是当 x1 - x0 很小时。”。

【解决方案7】：

如何计算类似的东西：

t = sign * power2 ( sqrt (abs(x - x0))/ sqrt (abs(x1 - x0)))

这个想法是使用一个数学等效公式，其中低 (x1-x0) 的影响较小。（不确定我写的是否符合这个条件）

【讨论】：

我将此表达式添加到测试代码中，它在平均误差和 stdDev 中给出了与其他方法相同的结果。所以没有好处。除此之外，它似乎效率不高。除非编译器足够聪明，可以删除 sqrt。
谢谢，关于平均误差，您是否检查过误差除以实际值？ 5.1 5.0和0.0 0.1是有区别的，后者与数值成比例误差较大。
这如何改善问题？
我很想尝试一下。你能告诉我我应该计算什么吗？我有区别，我应该用什么来划分它？浮点值，双精度值，有区别吗？

【解决方案8】：

正如 MSalters 所说，问题已经在原始数据中。

插值/外插需要斜率，它在给定条件下的精度已经很低（对于远离原点的非常短的线段最差）。

算法的选择无法恢复这种准确性损失。我的直觉是不同的评估顺序不会改变事情，因为错误是由减法而不是除法引入的。

想法：
如果在生成线条时有更准确的数据，可以将表示形式从 ((x0, y0), (x1, y1)) 更改为 (x0,y0, angle, length)。您可以存储角度或坡度，坡度有一个极点，但角度需要三角函数......丑陋。

当然，如果您经常需要端点，并且您的行数太多以至于无法存储其他数据，那当然行不通，我不知道。但也许还有另一种可以很好地满足您的需求的表示形式。

双打在大多数情况下都有足够的分辨率，但这也会使工作集加倍。

【讨论】：

我同意您关于输入数据精度的观点。但问题是解决计算过程。所描述的计算方法会破坏我们降低的精度吗？