在您的示例中,计算时引入了舍入误差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,是其他可能性。