平面幾何 [4/4]

前回に引き続き平面幾何の第4回です。今回扱っている話題は

直線と点の距離
線分と点の距離
線分の交差判定
線分の交点計算
直線の交点計算

です。

直線と点の距離

直線と点の距離は，点から「直線のどこか」への最短距離ですから，点から直線へ垂線を下ろせばその垂線の長さが求めたい「直線と点の距離」になります。それは，下図で示されているように，2点a,bを通る直線と点cとの距離は，ベクトル x=b-a と y=c-a を考えると距離は $|d| = |y| |\sin \theta|$ で求められます（sinはマイナスの場合もあるので絶対値を取っておきます）。

sinが入っているので外積の定義
$x \times y = |x||y| \sin \theta$
と見比べてみると，共通部分が多くあります。両辺を |x| で割ると，
$\frac{x \times y}{|x|} = \frac{|x||y| \sin \theta}{|x|} = |y| \sin \theta = d$
となるので，直線と点との距離は外積を使って $\frac{|x \times y|}{|x|}$ で求まります。コーディングは以下の通りです*1：

// 点a,bを通る直線と点cとの距離
double distance_l_p(P a, P b, P c) {
  return fabs(cross(b-a, c-a)) / abs(b-a);
}

線分と点の距離

今度は線分と点の距離を考えてみましょう。距離としてどのような値が欲しいのか，というのは問題依存なのですが，ここでは一般的な距離の定義に従って，点から「線分のどこか」への最短距離としてみます。そうすると，線分abに垂直な直線で点aを通る直線と点bを通る直線に囲まれた領域（下図の左の赤色領域に相当）にある点であれば，点から直線abへの垂線が最短距離になります。また，点cがこの赤色の領域の外であれば，点cがa側なら最短距離は2点a,cの距離と等しく，点cがb側なら最短距離は2点b,cの距離と等しくなることが分かります。

点cがaより内か外かを調べるには，前回の「点が線分上にあるかないか」の判定で用いた内積を利用した方法が使えます。すなわち，ベクトルx=b-aとy=c-aのcosが負であればaの外側ということが分かるので，その場合は点cと点aとの距離を最短距離とします（上図の右）。b側についても同様です。

コーディングは以下の通りです*2：

// 点a,bを端点とする線分と点cとの距離
double distance_ls_p(P a, P b, P c) {
  if ( dot(b-a, c-a) < EPS ) return abs(c-a);
  if ( dot(a-b, c-b) < EPS ) return abs(c-b);
  return fabs(cross(b-a, c-a)) / abs(b-a);
}

線分の交差判定

次に2線分が交差しているかどうかを考えます。

線分 a=a2-a1 と線分 b=b2-b1 が交差している場合，線分bの端点b1とb2は必ず線分aの両側にあります。また，線分aの端点a1とa2は必ず線分bの両側にあります（上図の左）。逆に交差していない場合，線分aか線分bのどちらかから見て，両点ともに片側にある場合が必ずあるといえます。上図の中央の例の場合，線分aから見ればb1,b2は線分aの両側にありますが，線分bから見ればa1,a2は両方とも線分bの右側にあります。これは交差していれば当然のことといえます。

ある点が線分（ベクトル）の左側にあるか右側にあるかというのは，外積を使えば判定できます（下図）。

sinは0度〜180度で正の値，180度〜360度で負の値を取るので，点がベクトルの左側なら外積は正，右側なら外積は負になります。従って，もし端点が両側にあるなら正と負，片側にあるなら正と正もしくは負と負になるので，2つの外積の値を掛けて正になれば非交差ということが分かります。また，線分aと線分bの両方について，相手線分の両端点との外積の値の積が負になれば交差ということが分かります。

以上より，2線分の交差判定のプログラムは次のようになります（外積の積が0であれば線上に点があることになりますが，これは交差しているとみなしています）：

// a1,a2を端点とする線分とb1,b2を端点とする線分の交差判定
int is_intersected_ls(P a1, P a2, P b1, P b2) {
  return ( cross(a2-a1, b1-a1) * cross(a2-a1, b2-a1) < EPS ) &&
         ( cross(b2-b1, a1-b1) * cross(b2-b1, a2-b1) < EPS );
}

線分の交点計算

次に，交差しているならばその交点を求めてみましょう。

上図より，線分aの始点a1からのパラメータtを求めれば，交点 p は $p = a1 + at$ で求めることができます。

パラメータtは，点a1から直線bの距離d1と点a2から直線bの距離d2との比 (t:1-t = d1:d2) から，
$t = \frac{d1}{d1 + d2}$
で計算できます。点と直線の距離は前述の通りで，外積を用いて求めます。

以上より，2線分の交点を求めるプログラムは次のようになります*3（前提：2線分は交差している）：

// a1,a2を端点とする線分とb1,b2を端点とする線分の交点計算
P intersection_ls(P a1, P a2, P b1, P b2) {
  P b = b2-b1;
  double d1 = fabs(cross(b, a1-b1));
  double d2 = fabs(cross(b, a2-b1));
  double t = d1 / (d1 + d2);

  return a1 + (a2-a1) * t;
}

なお，このプログラムはa1,a2を端点とする線分とb1,b2を通る直線の交点を求めることもできます（先に交差判定が必要です）。

直線の交点計算

最後に直線の交点計算を考えます。まず，2直線は平行でなければ必ず交差しており，交点を持ちます。平行判定は外積の性質である
$a \times b = 0 \Leftrightarrow a // b$
を使用すればOKです。一応プログラムで書けば

// a1,a2を通る直線とb1,b2を通る直線の交差判定
int is_intersected_l(P a1, P a2, P b1, P b2) {
  return !EQ( cross(a1-a2, b1-b2), 0.0 );
}

となりますが，平行時の例外処理は問題依存です。以降は，平行でない（＝どこかで交わっている）という前提で話をします。

さて，点a1,a2を通る直線と点b1,b2を通る直線の交点を求めるのですが，ここで提示する基本的な戦略は

ベクトル b=b2-b1 は直線として扱う
ベクトル a=a2-a1 を伸び縮みさせて，交点まで到達するベクトルを作る（結果，そのベクトルは交点座標を持つ）

ことです。

ベクトルbは直線として扱うので1通り（向きはありますが）に対して，ベクトルaは下図の3通りが考えられます。

まず直線bに垂直な軸を考えて，その軸上で上図に示す d1, d2 の長さを求めます。この直線bとその法線が成す直交座標系における d1:d2 の比は，元の座標系における目的のベクトル:元々のベクトルa の比と同じになります（目的のベクトルはa1を始点，交点を終点とするベクトル）。従って，交点を指す位置ベクトルは $a (d1 / d2)$ で求めることができます。

言い換えると，ベクトルaの大きさを d2 で割って正規化し，それを長さ d1 倍するとも考えられます。正規化するための係数 d2 は，上図より $d2 = b \times a$ で求めることができます。また，長さを表す係数 d1 は $d1 = b \times (b1-a1)$ で求めることができます。

なお，係数 d1 および d2 の正負は検証しておく必要があるでしょう。d1に関しては，左・中央の図では正，右の図では負になり，正常に動作しそうです（外積の正負に関しては「線分の交差判定」の図を参照）。d2に関しては，上図では全て正となるのでOKです。図を上下反転させる，すなわちa=a2-a1ではなくa=a1-a2ととった場合や，図の左右反転，すなわちb=b1-b2とする場合でも，d1およびd2は期待通りです。

以上より，次のプログラムが導かれます（前提：2直線は平行でない）：

// a1,a2を通る直線とb1,b2を通る直線の交点計算
P intersection_l(P a1, P a2, P b1, P b2) {
  P a = a2 - a1; P b = b2 - b1;
  return a1 + a * cross(b, b1-a1) / cross(b, a);
}

線分の交点計算は上図の中央のケースですから，当然このアルゴリズムは線分の交差判定でも使用できます。

以上で，平面幾何の解説を終わります。最後に，次のエントリで今回記述したソースコードを載せておきます。

*1:C++ではfabsの替わりにabsでも問題ありません。Cではabsはintに対してのみ有効なので気をつけてください。

*2:ここでは比較に EPS を使わずに 0.0 を用いても構わない（計算結果の距離はせいぜい計算誤差の範囲内になる）のですが，このようにEPSを使って垂線上の点も含めた方がこの関数の計算全体がより早く終了します。

*3:距離を求めるときに |b| で割っていないのは，d1 / (d1 + d2) で結局打ち消されるからです。