2

我们使用以下代码来计算两条线之间的锐角。 

def AcuteAngle2(line1,line2):
   ''':: line(x1,y1,x2,y2)'''
   u = (line1[2]-line1[0], line1[3]-line1[1])
   v = (line2[2]-line2[0], line2[3]-line2[1])
   return arccos(abs(dot(u,v)/(norm(u)*norm(v))))

它按预期工作。例如：

>>> AcuteAngle2([0,0,1,0],[0,0,0,1])
1.5707963267948966         #in rad = 90 degree

但是我们最近发现它在某些特殊情况下会失败！

>>> AcuteAngle2([0,0,1,0],[0,0,1,0])
0.0

这是正确的，但是：

>>> AcuteAngle2([0,0,1,1],[0,0,1,1])
2.1073424255447017e-08                #failed!

这是不正确的！它应该是 0.0。
有什么想法和解决方案吗？

更新 1：
在某些情况下， 使用Decimal答案中建议的包可能会有所帮助。然而，我们的问题仍未解决，因为 (1) 有很多代码需要大量时间来调整每个部分以供使用Decimal。此外，(2) 性能显着下降。此外，它需要 (3) 在处理numpy数组时进行大量更改。因此，它对我们的案例没有用。我们正在考虑某种装饰器等，而不改变事物并保持numpy性能。顺便说一句，有些人可能会建议使用多精度包，例如 gmpy 等，请注意，它们需要在代码中进行大量调整，这对我们的案例没有帮助。

4

3 回答 3

6

如果您关心准确性，则使用arccos锐角是一个坏主意。问题是对于接近 0 的角度的微小变化，该角度的余弦几乎不会改变。对于arccos情况相反 - 对于余弦角的非常非常小的变化变化更多。

在 2D 和 3D 中，更好的方法是使用atan2(crossproduct.length,scalarproduct)

在 2D 中，这变为atan2( dx1*dy2-dx2*dy1 , dx1*dy1+dx2*dy2 ). 请注意，您不需要对向量进行归一化，因此有两个改进：

  • 没有误差放大arccos
  • 没有额外的平方根误差
于 2013-06-08T04:44:28.353 回答
4

一种选择是使用小数模块来提高计算的精度：

from decimal import Decimal, getcontext

def AcuteAngle2(line1,line2):
   ''':: line(x1,y1,x2,y2)'''
   u = (Decimal(line1[2]-line1[0]), Decimal(line1[3]-line1[1]))
   v = (Decimal(line2[2]-line2[0]), Decimal(line2[3]-line2[1]))
   return arccos(float(abs(dot(u,v)/(norm(u)*norm(v)))))

看起来默认精度为 28 位，您将在这里得到预期的答案：

>>> getcontext().prec
28
>>> AcuteAngle2([0,0,1,1],[0,0,1,1])
0.0
于 2013-06-06T16:03:31.470 回答
2

在您的示例中，计算时引入了舍入误差dot(u,v)/(norm(u)*norm(v))。对于您的测试值，计算有效2/(sqrt(2)*sqrt(2))。sqrt(2) 的计算值被四舍五入为一个略大于无限精度值的值。

>>> math.sqrt(2)
1.4142135623730951
>>> math.sqrt(2)*math.sqrt(2)
2.0000000000000004
>>> 2/(math.sqrt(2)*math.sqrt(2))
0.9999999999999998
>>> math.acos(2/(math.sqrt(2)*math.sqrt(2)))
2.1073424255447017e-08

decimal@FJ的模块解决方案计算2/(sqrt(2)*sqrt(2))精度更高。当此值转换为浮点数（通过 arccos）时，它会四舍五入为 1.0。

>>> import decimal
>>> decimal.getcontext().sqrt(2)
Decimal('1.414213562373095048801688724')
>>> decimal.getcontext().sqrt(2)**2
Decimal('1.999999999999999999999999999')
>>> 2/decimal.getcontext().sqrt(2)**2
Decimal('1.000000000000000000000000001')
>>> float(2/decimal.getcontext().sqrt(2)**2)
1.0

2/(sqrt(2)*sqrt(2))使用和不同的精度计算decimal突出了另一个问题。

>>> for i in range(10,30):
...   decimal.getcontext().prec=i
...   print i,2/decimal.getcontext().sqrt(2)**2
... 
10 1.000000001
11 0.99999999995
12 1.00000000001
13 1
14 1
15 0.999999999999995
16 1
17 1.0000000000000001
18 1
19 0.9999999999999999995
20 1
21 1
22 0.9999999999999999999995
23 1
24 1
25 0.9999999999999999999999995
26 1.0000000000000000000000001
27 1.00000000000000000000000001
28 1.000000000000000000000000001
29 1

结果可能正好是 1、小于 1 或大于 1。如果您可以arccos不先四舍五入到浮点数，这可能会造成混淆。对于大于 1 的值，arccos未定义，因此结果为 a nan。如果在中间值超过 1 时看到这种类型的舍入错误会破坏纬度/经度计算。仅仅增加所有计算的精度，比如从 float64 到 float128，不会解决问题。它可能会将问题转移到一组不同的值，但仍会出现舍入错误。

更新 1

还有其他几个选项。您可以将公式重写为：

def AcuteAngle3(line1,line2):
   ''':: line(x1,y1,x2,y2)'''
   u = (line1[2]-line1[0], line1[3]-line1[1])
   v = (line2[2]-line2[0], line2[3]-line2[1])
   return arccos(sqrt(abs(dot(u,v)**2/(dot(u,u)*dot(v,v)))))

AcuteAngle3避免了你原来的问题，但它可能会dot(u,u)*dot(v,v)被舍入到一个比实际值略小的值，你可以尝试取大于 1 的值的 arccos。（但只对整个表达式使用 ROUND_UP 或 ROUND_DOWN不起作用；我在decimal示例中尝试使用不同的舍入模式，但仍然存在一些舍入“错误”。）

以下函数检查这些异常事件：

def AcuteAngle4(line1,line2):
   ''':: line(x1,y1,x2,y2)'''
   u = (line1[2]-line1[0], line1[3]-line1[1])
   v = (line2[2]-line2[0], line2[3]-line2[1])
   temp = sqrt(abs(dot(u,v)**2/(dot(u,u)*dot(v,v))))
   if temp > 1:
      return 0.0
   else:
      return arccos(temp)

计算中间值到更高的精度，然后向下舍入，或在计算表达式的每个分量时选择性地舍入或远离 0，是其他可能性。

于 2013-06-07T07:32:45.323 回答