3Dの計算にまつわるヨタ話し的なあれこれ

記述の正確さは保証しません。
もし誤りがありましたらご指摘いただけると幸いです。

Index


ゼロ近傍の闘い・または「ゼロの単位」

ゼロはどんな単位であっても「ゼロ」です。
このため、変数がゼロかどうかを判定する時にその単位を考えることを怠ってはいないでしょうか?
整数値であればゼロの単位が何か気にしなくとも全く問題ないのですが、浮動小数点値では問題になることがあります。 どういうことかといいますと、浮動小数点値がゼロかどうかを判定するのは

    f == 0.0

ではなくて、計算精度を考えて

    -0.00000001 < f && f < 0.00000001  ※1

のように、ゼロの近傍範囲でもって判定することになります。

そうです、ゼロに単位はなくとも、この判定用の「ゼロ近傍範囲」は単位に依存するのです。

図形に対する計算を行っていると この近傍範囲での浮動小数点値の判定を大量にやることになりますが、 「何の単位でもってゼロなのか」が揃っていないと、ある場合ではゼロと判定するのにある場合ではゼロと判定しない、 という不整合を産み、計算を間違えてしまう素になる恐れがあります。

単純な例で言うと、ミリメートルとキロメートルを混在して扱っている時に、ある時はミリメートル単位で※1の判定をし、 またある時はキロメートル単位に換算してから※1の判定をして、などということをやると、 ゼロ近傍範囲の大きさの意味が6ケタ分も変わってしまい、同じ長さなのにある時はゼロと判定するのにある時はゼロと判定しないといった 判定の不一致が生じて計算の不整合を産む要因になります。

これは非常に極端な例ですが、ちょっとした図形の計算でもこれと似たようなことが起り得ます。
以下、具体的な例で検証してみましょう。

ゼロ近傍の闘い1: 2つの2次元ベクトルが平行であるかどうかをチェックする

2次元のベクトルV1 (x1,y1) と、ベクトルV2 (x2,y2) が平行であるかどうかを調べる必要が出たとします。
いろいろな方法が考えられます。
  1. ベクトルの傾き y1/x1 と y2/x2 を計算し、同じなら平行
  2. ベクトルの各成分毎に一方からもう一方を割った値、x1/x2、y1/y2を求め、その値が同じなら平行
  3. 2つのベクトルを単位ベクトル化し、その各座標成分が全て同じ値であるか、絶対値が同じで全て符合が逆なら平行
  4. 片方のベクトルの単位化したベクトルを求め、もう片方のベクトルの単位化した直行ベクトルを求め、その内積を計算して0だったら平行
  5. 片方のベクトルの単位化した直行ベクトルを求め、もう片方のベクトルとの内積を計算して0だったら平行
  6. 片方のベクトルの直行ベクトル(長さは元と同じ)を求め、もう片方のベクトルとの内積を計算して0だったら平行 ( x1 * y2 - y1 * x2 == 0.0 )
  7. 2つのベクトルを単位ベクトル化し、その内積を計算して 1 か -1 だったら平行
たぶんまだまだあるでしょう。
4、5、6、は、直行する2つのベクトルの内積がゼロであることを利用してる点では同じですが、 両方を単位化するか、片方だけ単位化するか、あるいは両方とも単位化しないか、が異なります。

これらは全て数学的には正しい判定です。しかしながら、コンピュータでの計算上では非常に好ましくないものが含まれています。 あなたならどの方法を選ぶでしょうか?


上記7つの判定方法をひとつずつ見てみましょう。
以下、-0.0000001 < x && x < 0.0000001 のゼロ判定における '0.0000001' を「ゼロ近傍値」と呼ぶことにします。
1. ベクトルの傾き y1/x1 と y2/x2 を計算し、同じなら平行
中学校からお馴染みの「傾き」です。
これはこのままでは大きな問題を抱えています。まず、xが0だった場合は別途考えなければならないこと。 そしてより大きな問題は、傾きの大きさによってゼロ近傍の範囲が変わってしまうことです。

角度と傾きの関係は、θと tanθの関係になります。
で、タンジェントの値でもって同じ範囲で判定しても、角度で見た場合(左図で縦軸に対応する横軸の範囲を見た場合) 傾きが大きくなればなるほど角度ベースで見た場合の近傍範囲がどんどん小さくなってしまいます。これはうまくありません。
もしどうしても傾きベースでの判定をやりたいのであれば、ベクトルのx成分の絶対値とy成分の絶対値を比較し、 xの方が大きい場合は y/x を、yの方が大きい場合は x/y を考える、ということをやると 0除算問題も含めてかなり改善するかもしれません。

しかしながら、傾きは「長さ」を「長さ」で割るので「長さの比」、つまり無単位です。 無単位での判定は、ゼロ近傍値も単位としての意味を持たなくなります。
傾きでの判定はオススメしません。
傾きからは卒業しましょう。
2. ベクトルの各成分毎に一方からもう一方を割った値、x1/x2、y1/y2を求め、その値が同じなら平行
2つのベクトルが平行であれば各成分は同じ比率になる・・・ ので数学的には正しいのですが。
これもよろしくありません。

ベクトルV1がベクトルV2に比べてとてつもなく大きい場合を考えてみますと、成分比 x1/x2 は非常に大きい値になります。 ここでゼロ判定を考えてみましょう。x1/x2がとてつもなく大きな値の場合ゼロ近傍値は相対的にとてつもなく小さくなり、 ゼロかどうかの判定は非常にシビアになります。
逆に、ベクトルV1がベクトルV2に対してとてつもなく小さかったら成分比x1/x2は非常に小さな値になるのでゼロ近傍値は相対的に大きなものになり、 ゼロかどうかの判定はとてもゆるくなります。
これでは、ゼロかどうか(同じ値かどうか)の判定がものすごく曖昧なものになってしまいます。

さらにまた、長さを長さで割る「長さの比」、無単位に対してゼロ判定するので "ゼロ近傍値" は単位としての意味を持たなくなります。

「比での判定」はオススメしません。
ゼロと判定するかしないかが問題となるシビアな状況ではいろいろと厄介な問題を起こす元凶になります。
3. 2つのベクトルを単位ベクトル化し、その各座標成分が全て同じ値であるか、全て符合が逆なら平行
ベクトルで考えた場合の一番オーソドックスな解答がこれなんじゃないかと思います。

この方法では長さ1のベクトルに対して比較を行うため、ゼロ近傍の値の意味がベクトルの傾きや長さの比やその他諸々の要因によってほとんど変化しません。
で、その意味は何でしょう?
座標成分を調べてるから「距離」のように思うかもしれませんが、違います。 2つのベクトルは単位化されて長さは共に1。 2つのベクトルの長さが1で、座標成分がゼロ近傍値だけ離れた微妙な位置関係の場合、座標の差はベクトルの方向と直行する方に位置する2点の差になります。 これは、半径1の円の円周上の長さ、つまりラジアン角を考えるのに等しくなります。
実際には、円周方向の長さではなくx軸での範囲とy軸での範囲を考えてるので完全に角度ベースということではありませんが、意味的には「おおよそラジアン」と考えていいと思います。
x軸での近傍範囲とy軸での近傍範囲の作る矩形範囲で同一点判定を行うので、軸に平行な場合と対角線の45度の場合とで√2だけ判定の幅が変わりますが、 ケタが変わるとかのオーダーではないので許容範囲でしょう。

4. 片方のベクトルの単位化したベクトルを求め、もう片方のベクトルの単位化した直行ベクトルを求め、その内積を計算して0だったら平行
言葉で説明するとわかりずらいので図をみてください ^^;..

「直行するベクトルの内積を計算するとゼロになる」ということを利用します。
片方の直行ベクトルを求め、もう片方のベクトルとの内積を計算する。その値がゼロであれば、2つのベクトルは平行です。

2つのベクトルが直行している場合、それがどんな長さのベクトルでも内積を計算するとゼロになります。 あれ、じゃぁなんでわざわざ単位化するの? 意味あるの?  ・・と思われるかもしれませんが、ちゃんと意味があります。
内積は (V1の長さ)×(V2の長さ)× cosα ですから、2つが共に単位ベクトルで長さ1の場合、内積の値は単に cosα になります。 片方のベクトルの直行ベクトルに対して内積を計算しているので、元の2つのベクトルのなす角で考えると得られる値は sinθ (α+θ=90°)になります。

つまりこの計算は、“2つのベクトルのなす角θに対してsinθがゼロかどうか” の判定になります。 そして、θが限りなくゼロに近い場合 sinθ≒θが成り立ちますから、結局これは“2つのベクトルのなす角θがゼロかどうか”の判定です。 ゼロかどうかを判定するために -0.000001 < θ < 0.000001 とやった場合、この挟み込むゼロ近傍値の単位は「ラジアン」になります。

3の方法も「おおよそラジアン角ベース」でしたが、きっちりと角度ベースでの平行判定を行いたい場合はこっちの方がより正確です。
5. 片方のベクトルの単位化した直行ベクトルを求め、もう片方のベクトルとの内積を計算して0だったら平行
4.の場合と同じく、「直行するベクトルの内積を計算するとゼロになる」ということを利用します。 ただし、単位化するのは片方だけで、もう片方はそのままの長さで計算します。 片方だけ・・・ 単位化? ・・・ 意味あるの?・・・ と思われるかもしれませんが、これもちゃんと意味があります。図を見てみましょう。

内積を計算すると、(V1の長さ)×(V2の長さ)× cosθ はV1だけ単位化してるので (V2の長さ)×cosθ になります。
図をみるといっぱつですが、これはつまり、 「2つのベクトルの起点を揃えて、片方(単位化しなかった方)のベクトルからもう片方のベクトルへ降ろした垂線の長さ」です。
この長さがゼロであれば、平行、というワケです。
ここでゼロ判定 -0.000001 < l < 0.000001 をやると、この挟み込むゼロ近傍の単位は「距離」になります。

「直行ベクトルの内積を計算するとゼロ」という同じ手段を用いているのに、両方単位化してから計算するか、 片方だけ単位化してから計算するか、によって、"ゼロかどうか" を判定する時の単位が全く変わってしまいます。

さて、この4番と5番の2つの方法、どちらがより良い方法でしょうか?
それは「平行かどうか」の判定に何を求めているかによって変わります。
4番の角度ベースの判定では比較するベクトルがどんな長さであろうとも、ゼロ近傍の一定の角度内であれば平行と判定します。 これが望ましい場合もあるでしょう。また、同じ微小角度差であっても、ベクトルがどんどん長くなると隙間がどんどん開いてきます。 角度ではなく、この「開いた隙間」が問題となる場合は5の距離ベースでの判定の方がいいことになります。
6. 片方のベクトルの直行ベクトル(長さは元と同じ)を求め、もう片方のベクトルとの内積を計算して0だったら平行
さて、「直行するベクトルの内積を計算するとゼロになる」ということを利用した計算の3つ目。今度は、どちらも単位化しません。

この計算方法は非常に魅力的です。というのは、他の方法に比べてダントツに計算量が少ないのです。 単位化の計算は平方根の逆数を計算しなければなりませんし、1や2の比を使った計算でも割り算を含みます。 それに比べてこの6番目の方法は乗算と減算しか使いません。ベクトル V1 ( x1, y1 )、ベクトル V2 ( x2, y2 ) に対して

       x1 * y2 - y1 * x2

がゼロかどうか、を計算するだけです。ブラボー!

しかし、計算量を節約できたと喜んでいられるのはゼロ近傍問題が顕在化するまでの間だけです。
いや、顕在化したのならまだいい方で、多くの場合潜在的に問題を起こし、「数学的には正しい」アルゴリズムに潜むバグに頭を抱えることになります。

さて、この方法の何が問題なのでしょうか。

この方法の計算で得られる数値は2つのベクトルが作る平行四辺形の面積です。 2つのベクトルが平行であればどんなに長かろうと短かろうとこの面積はゼロになりますから、「ゼロだったら平行」の判定は数学的に何も間違っていません。

しかしながら、 -0.000001 < sq < 0.000001 の近傍値で挟み込む判定は面積の単位で行われることになります。

4の角度ベース、5の距離ベースでの判定が必要になることは多々あれど、この面積ベースでの判定が必要になるケースは滅多にありません。 そして「ゼロ近傍チェックの単位」を深く考えずにこの判定を使った場合、この「面積ベース」が悪さをします。面積ということはつまり、2乗のオーダーで値が変化するのです。

ほとんど平行だけど微小角だけ方向の異なる2つの1センチのベクトル、および、その同じ微小角だけ方向の異なる2つの100センチのベクトルを考えて、 それぞれをこの面積ベースでの平行判定を行った場合、後者の方の判定する面積は前者の1万倍になります。
これはつまりゼロ近傍値の意味が相対的に4ケタも変わってしまうということです。

大きなベクトルどうしの平行チェックではありえない勢いで判定が厳しくなり、小さなベクトルどうしの平行チェックではありえないくらいに判定がユルくなります。

図形の計算を行っていて全体としての整合性が重要な場合、このような判定基準の不一致や曖昧さが不整合を産む原因になります。

7. 2つのベクトルを単位ベクトル化し、その内積を計算して 1 か -1 だったら平行
これで最後です。
「平行な2つの単位ベクトルの内積は1(逆向きの場合 -1)である」というのを利用します。
これも数学的になら全く問題ありません。
なんとなく、4番の方法と似ています。

しかし、この方法は激しくおすすめしません。

この判定は、cosθが1かどうか、を判定することと同じです。
これがイケナイ。

サインカーブとコサインカーブの角度ゼロ付近を見てみると一発でわかります。

sinθがゼロかどうかの判定

  -0.000001 < sinθ < 0.000001

の場合、これに対応する角度θの範囲はおおよそ

  -0.000001 < θ < 0.000001

となります。角度ゼロ付近ではサインカーブの傾きは1で sinθ≒θが成り立ちます。
しかしながら、コサインカーブは角度ゼロ付近では傾きゼロです。
ゼロから角度θが増えていっても、cosθは暫くは1からなかなか減らないのです。
同じ精度で cosθが1かどうかの判定を行った場合、

  0.999999 < cosθ

に対応する角度θの範囲は大きく広がって

  -0.004472 < θ < 0.004472

となり、sinθ=0の判定と比べて3ケタもユルい角度で判定することになります。
cosθ=1の判定は一応は角度ベースでの判定なのですが、挟み込む1近傍範囲と、対応する角度の範囲に非常に大きな差があり、 全体的な精度を揃えるうえで面倒なことが起きる要因になり得ます。オススメしません。



さて、7通りの判定方法をみてきましたが、どうだったでしょうか。
数学的に考えるとどれでもOKなようにみえても、浮動小数点の同値判定問題が絡むおかげで平行かどうかの判定に大きな差が出てしまうことが わかっていただけたら幸いです。

ゼロ近傍の闘い2: 点が直線上に乗っているかどうかをチェックする

この問題は、上記の「闘い1:2つの2次元ベクトルが平行であるかどうかをチェックする」の応用です。
ですので、「距離でゼロ判定するか」「角度でゼロ判定するか」の2点に絞って簡単に説明します。

LP0と長さ1の法線ベクトルLNで表される直線、および、点Pがあり、点Pが直線上にあるかどうかを判定します。

まず、距離で判定する方から見てみましょう。

直線上の点LP0から点Pへのベクトル (P-LP0) を求めます。これは単位化しません。
これと、単位化されてる直線の法線ベクトル LN との内積を求めますと、出てくる答えは
「点Pから直線に降ろした垂線の長さ h 」になります。
点が直線上にある場合この h はゼロです。
ここで h がゼロかどうかの  -0.000001 < h < 0.000001 を行うと、つまり、

  直線から両側に 0.000001 の幅の帯の内側に点がある場合、その点は直線上とする

という判定を行うことになります。点が図の黄色い帯の中に点が入っていたら「直線上」というカンジで、これは非常に良い判定です。

次に角度での判定をみてみましょう。

同じく ベクトル (P-LP0) を求めて、今度は単位化します。
これと、単位化されてる直線の法線ベクトル LN との内積を求めますと、図での sinθに相当する値になります。
この sinθ がゼロであれば直線上であるとみなして、-0.000001 < sinθ < 0.000001 の判定を行いますと、θが限りなくゼロに近い場合はsinθ≒θが成り立って、-0.000001 < θ < 0.000001、 つまり、LP0 を中心として、直線から 0.000001ラジアンの角度内にある場合「直線上」という判定をすることになります。

これは図をみてすぐ判りますが、非常にマズい判定です。
点LP0の近傍ではどんなに直線に近くても「直線上でない」と判定してしまいます。
逆に点LP0から限りなく離れたところではかなり判定が緩くなってしまいます。


この2つはどちらも「法線ベクトルとの内積を計算してゼロかどうか」という判定をしています。
単位化するかしないか、ただそれだけの違いです。
しかしこの違いによってゼロかどうかを判定する数値の単位が全く変わってしまい、後者の方は非常に問題のある判定になってしまいます。


バウンダリBOXに関するかなりどうでもいい考察

図形全体を囲うバウンダリボックスを用意しといて、図形との細かな計算に入る前にバウンダリボックスと粗い計算をして計算対象かどうかチェックする のはよくやると思います。ところでこのバウンダリBOXってどういうデータ構造で持ちますか?

たぶん、以下の例のように、min点 max点の2点で表すのが一般的なんじゃないかと思います。
  struct POINT2{	// 2次元の点を表す構造体
    double x;
    double y;
  };

  struct BOUNDARY_BOX{	// 2次元のバウンダリを表す構造体
    POINT2 min;
    POINT2 max;
  };

「各軸の最小値を集めた点」と「各軸の最大値を集めた点」の2点で表現します。
しかしこれは、xの最小値と最大値、yの最小値と最大値、の4つの値がが欲しいだけなんだけど、それを2つの点の中に押し込めた、という感が拭えないのは私だけでしょうか。 min点とmax点という2つの「点」は、バウンダリボックスを構成する4頂点のうち2点という、それ以上の意味もそれ以下の意味もありません。

一方、私のプログラムではバウンダリは以下のような構造のものを用いています。
  // 1次元の範囲を表す構造体 RANGE
  struct RANGE{
    double min;
    double max;
  };

  // 2次元の範囲を表す構造体 BOUNDARY_BOX
  struct BOUNDARY_BOX{
    RANGE xRange;	// X軸方向の範囲
    RANGE yRange;	// Y軸方向の範囲
  };
一つの軸毎に「最小値と最大値が作る範囲」を考えて、それを軸の数だけ用意します。

なんだ、bounds.min.x ってのが bounds.xRange.min って「軸」と「最大最少」の順序が変わるだけじゃん、って思うかもしれません。
このちょっとした差が大きな差を産みます。

この違いを図にすると下図のようになります。

Aの一般的な構造の場合、先ず最初にバウンダリBOXという矩形を考えて、その矩形を表すためのmin点とmax点を指定するカンジです。 一方Bの場合は、X軸上の範囲とY軸上の範囲を考えると、その両方を満たすエリアが "結果的に" バウンダリボックスなる、という構造です。
このBによるバウンダリの考え方は、Cの「点P(px,py) はx軸上の値px とy軸上の値py の交差する点」という点の考え方と似ていて 「軸の直行性」をうまく表現しています。

バウンダリBOXに対する計算の多くは、軸毎に独立しています。
簡単な例を2つ挙げましょう。

「点がバウンダリBOXの中にあるか」を判定する場合、
  • 点のx座標がバウンダリBOXのX軸方向の範囲内にあり、かつ、
  • 点のy座標がバウンダリBOXのY軸方向の範囲内にある
という、X軸方向の判定とY軸方向の判定という2つの軸毎に独立した計算の AND が求める答えになります。

また、「2つのバウンダリBOXが共通部分を持つか」を判定する場合も同じく、
  • 2つのバウンダリの X軸方向の範囲が共通部分を持ち、かつ、
  • 2つのバウンダリの Y軸方向の範囲が共通部分を持つ
と、やはりX軸方向の判定とY軸方向の判定という2つの軸毎に独立した計算の AND が求める答えになります。
X軸の判定を行ってるときはY座標は全く関係ありませんし、その逆も同じです。
バウンダリBOXがこのような特性を持っているので、バウンダリBOXの構造も min点 max点で表すよりも x軸の範囲と y軸の範囲で表した方がよっぽど便利なことが多いです。

具体的に、「2つのバウンダリBOXが共通部分を持つかどうか」を判定するコードを比較してみましょう。

まず一般的な「min点 max点」によるバウンダリBOXのコードはこんなカンジでしょう。
// 2つのバウンダリBOXが共通部分を持つ場合に true
bool BoundsSectBounds( const BOUNDARY_BOX &b1, const BOUNDARY_BOX &b2 ){
	return (b1.min.x < b2.max.x) && (b2.min.x < b1.max.x) &&	// X軸方向で共通区間をもち、かつ
		(b1.min.y < b2.max.y) && (b2.min.y < b1.max.y);		// Y軸方向で共通区間をもつ
}
軸毎に独立した計算をしているのに、「各軸の最少値を集めた点」と「各軸の最大値を集めた点」という 直行性を全く活かしてない構造で持っているためこれ以上コードをスッキリさせることができません。
min.x、min.y、max.x、max.y .... と、ごちゃごちゃして紛らわしいコードになりがちです。

一方、「X軸の範囲、Y軸の範囲」によるバウンダリBOXではこうなります。
// 2つの “一次元の範囲” が共通部分を持つ場合に true
bool RangeSectRange( const RANGE &r1, const RANGE &r2 ){
	return r1.min < r2.max && r2.min < r1.max;
}

// 2つのバウンダリBOXが共通部分を持つ場合に true
bool BoundsSectBounds( const BOUNDARY_BOX &b1, const BOUNDARY_BOX &b2 ){
	return RangeSectRange( b1.xRange, b2.xRange ) && // X軸方向の2つの範囲が共通区間をもち、かつ、
		RangeSectRange( b1.yRange, b2.yRange );	 // Y軸方向の2つの範囲が共通区間をもつ
}
各軸毎に同じ計算をして、さらに各軸毎に構造体がまとまっているので、各軸毎の計算を行う部分を下請け関数にまとめることができました。コードが読みやすくなったとおもいませんか?

さらにz軸を加えて3次元のバウンダリCUBEを考えたり、境界判定を厳しくして「境界を含む場合」と「境界を含めない」場合を区別したり・・・   などと、計算量が増えたり判定パターンが増えたりする度に、このちょっとした差がコード全体の可読性の差になって表れます。





copyright©2008 dendrocopos All rights reserved.