平面幾何 [3/4]

ACM/ICPC

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

  • 内積・外積
  • 2直線の直交判定・平行判定
  • 点が線上にあるかないかの判定

です。

内積・外積

ベクトルの内積 (inner product, dot product, scalar product) と外積 (outer product, cross product, vector product) という演算を用いると幾何の問題を解く考え方が簡単になります。幾何学における内積や外積はもともと3次元空間上で定義されるものなので，まずは3次元空間上で幾何学的な内積・外積を導入し，それらが線形代数的なベクトル演算と等価であることを利用し，内積・外積を2次元平面上に拡張（縮小？）します。

3次元空間上において，ベクトルの内積（ドット積）は a \cdot b で表され，以下の式で定義されます：
a \cdot b = |a||b| \cos \theta
内積はcosを使って定義されている点と，内積の結果は単一の値＝スカラーになる点に注意してください。

また，ベクトルの外積（クロス積）は a \times b で表され，その大きさ |a \times b| は
|a \times b| = |a||b| \sin \theta
で与えられ，その向きはa,bを含む平面に垂直で，aからbへの右ねじの進む向きになります。外積はsinを使って定義されている点と，外積の結果はベクトルになる点に注意してください。なお，定義より，外積の大きさはベクトルa,bが成す平行四辺形の面積になります。

また，内積は上式の他，次のベクトル演算によって求めることもできます：
a \cdot b = a_x b_x + a_y b_y + a_z b_z

外積も行列式 (determinant) を用いて次のように求めることができます（eは単位ベクトル）：
a \times b = \begin{vmatrix} e_x & e_y & e_z \\ a_x & a_y & a_z \\ b_x & b_y & b_z \end{vmatrix} = (a_y b_z - a_z b_y, a_z b_x - a_x b_z, a_x b_y - a_y b_x)

これらの式はcos/sinを用いた前式と等価ですが，計算は簡単になります。さらに，これらは3次元空間に定義されているので，z 値に 0 を代入すれば2次元平面上の内積・外積の式を導くことができます：
a \cdot b = |a||b| \cos \theta = a_x b_x + a_y b_y
a \times b = |a||b| \sin \theta = \begin{vmatrix} a_x & a_y \\ b_x & b_y \end{vmatrix} = a_x b_y - a_y b_x

cos/sinを用いた定義はほぼ上式と同じですが，外積の結果は（ベクトルではなく）スカラーになっていることに注意してください*1。

ここではプログラム上では最後の2式を採用し，次のように関数を定義することにします（内積・外積はその演算子の形からドット積，クロス積とも呼ばれており，それがこれらの関数名の由来です*2）：

// 内積 (dot product) : a・b = |a||b|cosΘ
double dot(P a, P b) {
  return (a.real() * b.real() + a.imag() * b.imag());
}
// 外積 (cross product) : a×b = |a||b|sinΘ
double cross(P a, P b) {
  return (a.real() * b.imag() - a.imag() * b.real());
}

ベクトルの内積・外積は様々な用途に使えるので，その性質と使い方を順に見ていきましょう。

2直線の直交判定・平行判定

2次元平面上のベクトルの内積・外積の幾何学的な定義をもう一度見てみましょう：
a \cdot b = |a||b| \cos \theta
a \times b = |a||b| \sin \theta

内積にはcos()が，外積にはsin()が使われています。cos()は90度と-90度のときに値が0となるので，2ベクトルが直交しているときに内積は0になることが分かります。ですから，ベクトルの直交判定に使うことができます。同様に，sin()は0度と180度のときに値が0となるので，外積はベクトルの平行判定に用いることができます。数式でこのことを宣言すれば，
a \cdot b = 0 \Leftrightarrow a \perp b
a \times b = 0 \Leftrightarrow a // b
ということです。

2点a1,a2を通る直線または線分はベクトルa=a1-a2で表すことができるので，直交判定と平行判定をするプログラムは次のようになります：

// 2直線の直交判定 : a⊥b <=> dot(a, b) = 0
int is_orthogonal(P a1, P a2, P b1, P b2) {
  return EQ( dot(a1-a2, b1-b2), 0.0 );
}
// 2直線の平行判定 : a//b <=> cross(a, b) = 0
int is_parallel(P a1, P a2, P b1, P b2) {
  return EQ( cross(a1-a2, b1-b2), 0.0 );
}

なお，ここで関数 EQ は，

#define EPS (1e-10)
#define EQ(a,b) (fabs((a)-(b)) < EPS)

と定義しています。これは，小数は計算誤差（丸め誤差）を持つために == では比較できないため，非常に小さい数 EPS 以下の差であれば同じとみなす，という処理をしています。ここでは EPS は 1 e^{-10} としていますが，この値は問題に合わして変更する必要があるかもしれません。とにかく，コンピュータで小数を扱う場合は常に計算誤差を意識する必要があります。

点が線上にあるかないかの判定

2点a,bを通る直線と点cが与えられたときに，点cが直線ab上にあるかないかを判定することを考えてみましょう（下図の左）。

ここで直線を表すベクトル x=b-a と，aからcに伸びるベクトル y=c-a を考えます（上図の右）。上図から明らかなように，もし点cが直線状にあればベクトルxとベクトルyは重なり，点cが直線状になければ重なりません。ベクトルが重なっているということは，それらのベクトルは必ず並行になっているといえます。従って，点cが直線ab上にあるかどうかは，前述の平行判定（外積が0かどうか）と同様に，

// 点cが直線a,b上にあるかないか
int is_point_on_line(P a, P b, P c) {
  return EQ( cross(b-a, c-a), 0.0 );
}

と書くことができます。

次に，直線ではなく線分の場合を考えてみましょう。つまり，2点a,bを端点とする線分と点cが与えられたときに，点cが線分ab上にあるかないかを判定するということです。

まず，上述の直線の場合は必ず満たしていなければなりません。しかし，線分には長さがあるので，直線ab上に乗っていても線分abの範囲からは外れている可能性があります。

単純に思い浮かぶのは，点a,bと点cの座標を比較する方法です。点cが線分上にあるためには，点a,bを対角点とする長方形の中になければならないので（下図の左），長方形の内か外を調べるコードを次のように付け足せばよいことになります（最初の行は，点cが点aや点bに重なる場合，後の比較が（簡潔な表現をするためにEPSを用いていないので）計算誤差によって誤って判定されることを防いでいます）：

// 点cが線分a,b上にあるかないか(1)
int is_point_on_line(P a, P b, P c) {
  if ( EQ(abs(a-c), 0.0) || EQ(abs(b-c), 0.0) ) return 1;

  if (a.real() < b.real()) {
    if (c.real() < a.real() || b.real() < c.real) return 0;
  } else {
    if (c.real() < b.real() || a.real() < c.real) return 0;
  }
  
  if (a.imag() < b.imag()) {
    if (c.imag() < a.imag() || b.imag() < c.imag) return 0;
  } else {
    if (c.imag() < b.imag() || a.imag() < c.imag) return 0;
  }

  return EQ( cross(b-a, c-a), 0.0 );
}

また，内積を使う（ちょっとややこしい）方法もあります。cos()の結果は-90度〜90度（右半円）であれば正，90度〜270度（左半円）であれば負となります（上図の右）。ですから，ベクトル x=b-a と ベクトル y=c-a の内積が正となれば，点aよりも線分側であると判定できます。それを点bと点cについても行えばいいので，次のようなプログラムが書けます：

// 点cが線分a,b上にあるかないか(2)
int is_point_on_line(P a, P b, P c) {
  return EQ( cross(b-a, c-a), 0.0 ) &&
         (dot(b-a, c-a) > -EPS) &&
         (dot(a-b, c-b) > -EPS);
}

計算誤差に対処するために，0.0 ではなく -EPS と比較していることに注意してください。

あるいは，三角不等式を使う簡潔な方法もあります。もし点cが線分上にあれば a->b の道のりと a->c->b の道のりは等しく，もし点cが線分上になければ　a->b の道のりよりも a->c->b の道のりの方が長くなります。

このことを素直にコーディングすれば，

// 点cが線分a,b上にあるかないか(3)
int is_point_on_line(P a, P b, P c) {
  // |a-c| + |c-b| <= |a-b| なら線分上
  return (abs(a-c) + abs(c-b) < abs(a-b) + EPS);
}

となります。

さて，次回は最後の仕上げとして，直線（線分）と点の距離，線分の交差判定と交点計算を取り上げます。

*1:2次元平面上の外積の結果は，必ず（3次元で言えば）zの値のみを持つベクトルとなるのでわざわざベクトルとはせずにスカラーとして扱う，という解釈でも良いでしょう。

*2:外積はクロス積のcrossではなく行列式のdetにしようかとも思いましたが，こちらの方が記憶しやすそうで混同しなさそうなのでcrossを採用しました。