幾何ライブラリ解説2 - プログラミング雑記帳

前回

距離
交点

距離

点と直線の距離

高校数学なのですが、射影を使うよりも、外積を使って値を出す方が楽なので一応説明します。
まず、点を $P$ 、直線上の異なる二点を $L_1$ , $L_2$ とし、また求める解を $d$ とします。
$\vec{a} = \overrightarrow{L_1L_2}$ , $\vec{b} = \overrightarrow{L_1P}$ とし、 $\vec{a}$ , $\vec{b}$ のなす角を $\theta$ とすると、
$\vec{a} \times \vec{b} = \|\vec{a}\|\|\vec{b}\|sin\theta$ より、
$\displaystyle \frac{\vec{a} \times \vec{b}}{\|\vec{a}\|} = \|\vec{b}\|sin\theta = \|\vec{b}\|\frac{d}{\|\vec{b}\|} = d$ となります。

double getDistanceLP(Line l, Point p) {
	return abs(cross(l.p2 - l.p1, p - l.p1) / abs(l.p2 - l.p1));
}

点と線分の距離

まず、点を $P$ 、線分の二端点を $S_1$ , $S_2$ とします。
射影を求めた時にその点が線分上にあれば、点と直線の距離と同じようにできます。
しかし、その点が線分上にない場合、点と線分の二端点の近い方の距離が解となります。
この判定は内積を使い正負を見ることで分かります。
$\overrightarrow{S_1S_2} \cdot \overrightarrow{S_1P}$ の値が負ならば $P$ と $S_1$ との距離を求めればよいです。
同様に、 $\overrightarrow{S_2S_1} \cdot \overrightarrow{S_2P}$ の値が負ならば $P$ と $S_2$ との距離を求めればよいです。

double getDistanceSP(Segment s, Point p) {
	if (dot(s.p2 - s.p1, p - s.p1) < 0.0) return abs(p - s.p1);
	if (dot(s.p1 - s.p2, p - s.p2) < 0.0) return abs(p - s.p2);
	return getDistanceLP(s, p);
}

線分と線分の距離

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/2/CGL_2_D
線分が交差していれば勿論0です。
そうでないときは線分と片側の線分の端点との距離をそれぞれ4通り求め、それの最小値を出せばよいです。

double getDistanceSS(Segment s1, Segment s2) {
	if (intersectSS(s1, s2)) return 0.0;
	return min({ getDistanceSP(s1, s2.p1), getDistanceSP(s1, s2.p2),
			   getDistanceSP(s2, s1.p1), getDistanceSP(s2, s1.p2) });
}

交点

直線と直線の交点

下の図を用いて説明しますが、直角や平行などを書き忘れたのでなんとなく分かってもらえると助かります。
一つ目の直線の任意の二点を $L_1$ , $L_2$ とし、二つ目の直線の任意の二点を $L_3$ , $L_4$ としています。
$\displaystyle \vec{L_1} + \frac{d_2}{d_1}\overrightarrow{L_1L_2}$ が求めたい答えになることが図より分かります。
これは相似から簡単に証明することができます。
したがって、 $\displaystyle \frac{d_2}{d_1}$ を求めます。
$\displaystyle \frac{\overrightarrow{L_1L_3} \times \overrightarrow{L_3L_4}}{\overrightarrow{L_1L_2} \times \overrightarrow{L_3L_4}} = \frac{\|\overrightarrow{L_1L_3}\|\|\overrightarrow{L_3L_4}\|sin\varphi}{\|\overrightarrow{L_1L_2}\|\|\overrightarrow{L_3L_4}\|sin\theta} = \frac{\|\overrightarrow{L_1L_3}\|sin\varphi}{\|\overrightarrow{L_1L_2}\|sin\theta} = \frac{d_2}{d_1}$
と求めることができます。
直線が一致する場合の処理は問題に依存するので適宜変えますが、別の問題で使うためその時にまた説明します。

上記の説明とコードがほんの少しだけ違いますが問題ないのでこのままにします。

Point getCrossPointLL(Line l1, Line l2) {
	double a = cross(l1.p2 - l1.p1, l2.p2 - l2.p1);
	double b = cross(l1.p2 - l1.p1, l1.p2 - l2.p1);
	if (abs(a) < EPS && abs(b) < EPS) return l2.p1;
	return l2.p1 + (l2.p2 - l2.p1) * (b / a);
}

線分と線分の交点

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/2/CGL_2_C
交点を持つことを仮定している場合、直線と直線の交点で代用できます。
そうでない場合でも、線分が交差しているという仮定の下だと簡単に考えることができるので、それを書きます。
下の図を用います。
一つ目の線分の二端点を $S_1$ , $S_2$ とし、二つ目の線分の二端点を $S_3$ , $S_4$ としています。
$PS_1 : PS_2 = t : 1 - t$ とすると、求める交点は $\vec{S_1} + t\overrightarrow{S_1S_2}$ となります。
また相似より、 $t : 1 - t = d_1 : d_2$ なので、 $\displaystyle t = \frac{d_1}{d_1 + d_2}$ と $d_1$ , $d_2$ が分かれば求められます。
$d_1$ , $d_2$ は既に説明した点と直線の距離を使うことで求められますが、除算の回数を減らすために関数を呼び出さずに似たことをしています。

Point getCrossPointSS(Segment s1, Segment s2) {
	Vector base = s2.p2 - s2.p1;
	double d1 = abs(cross(base, s1.p1 - s2.p1));
	double d2 = abs(cross(base, s1.p2 - s2.p1));
	return s1.p1 + (s1.p2 - s1.p1) * (d1 / (d1 + d2));
}

円と直線の交点

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/7/CGL_7_D
高校数学なので説明は省略します。円の弦の長さを求めてちょっとやるだけです。

vector<Point> getCrossPointCL(Circle c, Line l) {
	vector<Point> ps;
	Vector pr = project(l, c.c);
	Vector e = (l.p2 - l.p1) / abs(l.p2 - l.p1);
	if (equals(getDistanceLP(l, c.c), c.r)) return vector<Point>{pr, pr};
	double base = sqrt(c.r * c.r - norm(pr - c.c));
	ps.push_back(pr + e * base); ps.push_back(pr - e * base);
	return ps;
}

円と線分の交点

線分と円の交点が2つある場合や接する場合、円と直線の交点を出力すればよいです。
そうでなければ、つまり交点が1つのとき、円と直線の2交点のうちのどちらかを出力します。これは、2交点のうちどちらの交点が線分上に
あるかを見ればよいです。

vector<Point> getCrossPointCS(Circle c, Segment s) {
	Line l(s);
	vector<Point> ps = getCrossPointCL(c, l);
	if (intersectCS(c, s) == 2) return ps;
	if (dot(l.p1 - ps[0], l.p2 - ps[0]) < 0) ps[1] = ps[0];
	else ps[0] = ps[1];
	return ps;
}

円と円の交点

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/7/CGL_7_E
一つ目の円の中心点と半径を $C_1$ , $r_1$ 、二つ目の円の中心点と半径を $C_2$ , $r_2$ とします。
また、求める2交点を $P_1$ , $P_2$ とし、 $a = \angle P_1C_1C_2 = \angle P_2C_1C_2$ （これは合同な三角形であることから導けます)とします。
$C_1$ と $C_2$ の距離を $d$ とすると、
余弦定理より、 $\displaystyle cos\angle P_1C_1C_2 = \frac{d^2 + r_1^2 - r_2^2}{2dr_1}$ とできるので、これを逆三角関数に入れ、角度を出したものを $a$ とします。
それと、 $\overrightarrow{C_1C_2}$ と $x$ 軸のなす角を $t$ とします。
上記で求めた角度を用いることで、
$\displaystyle \vec{P_1} = \vec{C_1} + \begin{bmatrix} r_1cos(t + a) \\ r_1sin(t + a) \\ \end{bmatrix}$ ,
$\displaystyle \vec{P_2} = \vec{C_1} + \begin{bmatrix} r_1cos(t - a) \\ r_1sin(t - a) \\ \end{bmatrix}$
と2交点を求めることができます。

vector<Point> getCrossPointCC(Circle c1, Circle c2) {
	double d = abs(c1.c - c2.c);
	double a = acos((c1.r * c1.r + d * d - c2.r * c2.r) / (2 * c1.r * d));
	double t = arg(c2.c - c1.c);
	vector<Point> ps;
	ps.push_back(c1.c + polar(c1.r, t + a)); ps.push_back(c1.c + polar(c1.r, t - a));
	return ps;
}