2線分の交点座標(2次元)

Aizu Online JudgeのCoursesを埋めていたところ、2線分の交点を求める問題に出会った。
そこで2線分の交点導出方法を考える。

2線分の交点導出アルゴリズム

ここでは同一平面上に存在し、並行でない線分 $AB, CD$ について考える。
4点 $A, B, C, D$ の2次元座標が与えられたときの交点 $X$ の座標を求めたい。

点 $X$ は線分 $AB, CD$ 上に存在するため媒介変数 $s, t$ を用いて

X = A + s\vec{AB} = C + t \vec{CD}

と表現できる。

$\vec{AB} = B - A, \vec{CD} = D - C$ であるため、各点に関して $x, y$ 座標の関係式が求まる。

\begin{equation}
\left \{
\begin{array}{l}
A_x + s(B_x - A_x) = C_x + t(D_x - C_x) \\
A_y + s(B_y - A_y) = C_y + t(D_y - C_y)
\end{array}
\right.
\end{equation}

これを $s, t$ に関して解く。

s = \frac{(C_x - A_x)(D_y - C_y) - (C_y - A_y)(D_x - C_x)}{(B_x - A_x)(D_y - C_y) - (B_y - A_y)(D_x - C_x)} \\
t = \frac{(A_x - C_x)(B_y - A_y) - (A_y - C_y)(B_x - A_x)}{(D_x - C_x)(B_y - A_y) - (D_y - C_y)(B_x - A_x)}

線分 $AB, CD$ が並行でないとき、この式より2線分の端点座標から $s, t$ の値が求まり、交点の座標を計算できる。

2線分の交差判定

$s, t$ の値から2線分の交差判定を行うことができる。
直線 $AB, CD$ の交点が線分 $AB, CD$ 上に存在するとき線分 $AB, CD$ の交点が存在するといえる。
したがって $s, t$ が共に $[0, 1]$ の範囲に収まっているとき線分 $AB, CD$ は交差する。

0 \leq s, t \leq 1 を満たすとき交差する

パラメータ s, t の幾何的意味

2線分の交点導出アルゴリズム で得られた $s, t$ の式を見ると分母分子が共に外積の形をしていることに気づく。

s = \frac{(C_x - A_x)(D_y - C_y) - (C_y - A_y)(D_x - C_x)}{(B_x - A_x)(D_y - C_y) - (B_y - A_y)(D_x - C_x)}
  = \frac{|\vec{AC} \times \vec{CD}|}{|\vec{AB} \times \vec{CD}|}

2次元ベクトルにおける外積の大きさは2ベクトルで成す並行四辺形の面積に等しい。
ここで $\vec{AC}, \vec{AB}, \vec{CD}$ の関係を図示する。

$|\vec{AC} \times \vec{CD}|$ と等しい面積を持つ並行四辺形を緑、
$|\vec{AB} \times \vec{CD}|$ と等しい面積を持つ並行四辺形を青で表した。

図に表現してみると $|\vec{AC} \times \vec{CD}| = |\vec{AX} \times \vec{CD}|$ であることがよく理解できる。
またパラメータ $s$ が青と緑の平行四辺形の面積比を表していることがわかる。

C++ ソースコード

引数に4点 $A, B, C, D$ の座標を渡すと2線分 $AB, CD$ の交点の座標を返す関数IntersectionをC++で書く。

#define INF 1e60
typedef struct Point_Coordinates {
	double x, y;
	friend Point_Coordinates operator-(const Point_Coordinates& l, const Point_Coordinates& r) {
		return { l.x - r.x, l.y - r.y };
	}
} point;

double Cross(const point& a, const point& b) {
	return a.x * b.y - a.y * b.x;
}

point Intersection(const point& a, const point& b, const point& c, const point& d) {
	double s, t, deno = Cross(b - a, d - c);
	point error = { INF, INF };
	if (deno == 0.0) {
		// 線分が平行
		return error;
	}
	s = Cross(c - a, d - c) / deno;
	t = Cross(b - a, a - c) / deno;
	if (s < 0.0 || 1.0 < s || t < 0.0 || 1.0 < t) {
		// 線分が交差していない
		return error;
	}
	return { a.x + s * (b - a).x, a.y + s * (b - a).y };
}

まとめ

図で表すとその式が出てくる理由が透けることがあって面白いです。
AOJ埋めのんびり頑張ります。