2

次のコードを使用して、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) パフォーマンスが大幅に低下します。numpyさらに、(3)配列を扱う際に大規模な変更が必要です。したがって、私たちの場合には役に立ちません。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 )ます。ベクトルを正規化する必要がないことに注意してください。そのため、2 つの改善点があります。

  • によるエラー増幅なしarccos
  • 平方根による追加誤差なし
于 2013-06-08T04:44:28.353 に答える
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

@FJdecimalによるモジュール ソリューションは2/(sqrt(2)*sqrt(2))、より高い精度で計算されます。この値が (arccos によって) float に変換されると、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は定義されていないため、結果は になり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 より大きい値のアークコスを取得しようとする可能性があります (ただし、式全体に 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 に答える