３．直線、線分の計算幾何 シリーズ：「基礎的計算幾何ライブラリの作り方」

シリーズの構成

１．計算幾何の基本
２．計算幾何のインフラ整備
３．直線、線分の計算幾何　　　←いまここ
４．円の計算幾何
５．多角形の計算幾何
６．AOJにある問題の解説

今回の内容

　今回は、直線や線分についての計算幾何を説明します。計算幾何の基本的な単位であるまっすぐな線について、交点、距離などを求めますが、これらは今後の円や多角形の計算幾何の基礎になるものが多いです。ぜひ身につけましょう。
　
　また、今回の特徴として、できるだけベクトルの内積や外積でアプローチする、というものです。誤差を産みづらく、コーディング量も減ってlibrary使用できないコンテストでも簡単に復元できるからです。

全体的な構造

　今回の幾何ライブラリは、すべてgeometry2dネームスペース下で書いていましたが、今回の内容、すなわちとりわけ直線や線分にかかわる部分はgeometry2d直下にline2dというネームスペースを作り、そこで直線の構造体やさまざまな機能を定義します。　

　ちなみに、ライブラリの全容の詳細はここにあります。全体図把握のために一読することを勧めます。
幾何ライブラリの構造.md - Google ドライブ

今回実装する機能

　ネームスペースgeometry2d直下で次の機能を実装します。

• ネームスペースline2d

　このline2dのネームスペース内で次の機能を実装します。

• 直線を表すLine構造体。

• lineIntersection()　与えられた直線の交点を返す　存在しないならerror_val
• segmentIntersection()　与えられた線分の交点を返す　複雑な仕様あり
• distanceBetweenPointAndLine()　与えられた点と直線の距離を返す
• distanceBetweenPointAndRay()　与えられた点と半直線の距離を返す
• distanceBetweenSegmentAndSegment()　与えられた線分と線分の距離を返す
• projection()　与えられた点と直線について、正射影を返す
• reflection()　与えられた点と直線について、線対称の点(鏡映)を返す

Line構造体の具体的な実装

　直線は、通る2点を指定すれば一意に定まることから、直線の通る2点(Point型)を構造体に持たせます。普通の直線には関係のないことですが、向きをつけるために片方の点の名前をbegin、もう片方の点の名前をendとしましょう。もちろん普通の直線を扱うのには、この2点の順番に関係はありません。

コンストラクタ

　コンストラクタとしては、

• 何も入力されなかったら、直線の2点に適当な値を入れる
• 2点を入力されて、それに該当する直線を作る
• という式から直線を作る

を用意します。第一話でも第二話でも言いましたが、計算幾何は基本的にはベクトルや平面幾何を用いたアプローチをします。なぜ代数的なアプローチは避けるのかというと、

• 例外の処理が面倒(直線の場合、y=ax+bという形で表せないタイプには例外処理が必要)
• 連立させて方程式解くことに帰着できるが、コードがごちゃごちゃになってわかりづらくなる
• 平面幾何やベクトルの方が視覚的にも理解しやすく、コーディング量が減り覚えやすいことが多い

などが挙げられる。というわけで、代数的な式から、通る2点を抽出するいわばベクトルの形に直すのである。

　具体的な直し方はコードを見てください。y軸と平行な直線(x=定数)は、傾きを定義できないため場合分けすればよい。抽出する二点は、ある程度(1.0ぐらい)離れてる2点を選ぼう。

クリックで展開

struct Line {
//直線の通る二点。有向直線でないのならば、beginやendに違いはない。
	Point begin, end;
	Line() {
		begin = Point(), end = Point();
	}

	Line(const Point& b, const Point& e) {
		begin = b, end = e;
	}

	//ax+by+c=0
	Line(const double a, const double b, const double c) {
		if (sgn(a) == 0 && sgn(b) == 0) {
			assert(-1);//警告を出してプログラムを終了させる。
		}
		if (sgn(b) == 0) {
			//ax+c=0になる。
			begin = Point(-c / a, 0.0);
			end = Point(-c / a, 1.0);
		}
		else {
			//y=-(ax+c)/b　傾きは-a/bで、y切片が-c/b
			begin = Point(0, -c / b);
			end = Point(1.0, -(a + c) / b);
		}
	}
};

vec(), countervec()　直線をベクトル化する

　平面幾何は、直線をベクトル化して、内積や外積などのベクトル演算の恩恵に預かるように実装されてる場合が多いです。それゆえ、直線を代表するベクトルを要求される場合が多いですが、毎回(line1.end - line1.begin)と書くのは冗長です。

　というわけで、vec()では、end - beginを返しcountervec()ではその逆のベクトルのbegin - endを返しましょう。

クリックで展開

Point vec() const {
	return end - begin;
}

Point countervec() const {
	return begin - end;
}

半直線、線分

　半直線や線分は直線とよく似ています。というか直線の一部です。
　
　ですので、次のように定めて、この直線の構造体を流用します。

• 半直線は、beginから始まり、end方向へと無限に通る。
• 線分は、beginから始まり、endで止まる。

　同じLineという構造体でも、名前は、半直線はRay、線分はSegmentと別名をつけて区別します。これは、ヒューマンエラーを減らすのに役に立ちます。

クリックで展開

//半直線なら、Lineと同じだけど、はっきりと区別する。
typedef Line Ray;
//線分の場合、Lineと同じものになるが、はっきりと区別する
typedef Line Segment;

具体的な機能の実装

lineIntersection()　直線の交点

　一番シンプルの、直線同士の交点です。直線をベクトル表示して、変数同士の方程式を解くかと思いきや、これは平面幾何的に次のように解けます。

　外積で求めた、底辺が同じ平行四辺形の面積比は、それらの高さの比、という事実を利用して求めています。この場合は、煩雑なコードになることなく簡潔に実装できます。[出典1]

クリックで展開

//直線の交点を返す。交わってなければ、error_valを返す。
Point lineIntersection(const Line& l1, const Line& l2) {
	if (sgn(l1.vec().cross(l2.vec())) == 0)return error_val;
		Point ret;
	ret = l1.begin + l1.vec() *
		abs((l2.end - l1.begin).cross(l2.vec()) / l1.vec().cross(l2.vec()));
	return ret;
}

segmentIntersection()　線分の交点

　返す要素としては、交わってるor交わってないと、交点となります。(一部がピッタリ重なってる場合は、交点が存在しないけど交わってるので、交わってるけど交点はerror_val、というのを返しておく。)

クリックで展開

//線分が共通部分を持つかどうか？と線分の交点を返す。共通部分がない、もしくは交点が一意ではないなら、error_valを返す。
//trueなら、共通部分を持つ。falseなら、共通部分を持たない。
pair<bool, Point> segmentIntersection(const Segment& s1, const Segment& s2) {
	if (iSP(s1.begin, s1.end, s2.begin) * iSP(s1.begin, s1.end, s2.end) <= 0 &&
	iSP(s2.begin, s2.end, s1.begin) * iSP(s2.begin, s2.end, s1.end) <= 0) {
	//平行ならば、交点は定まらない。(完全に重なってるので)
	if (s1.vec().cross(s2.vec()) == 0)
		return make_pair(true, error_val);
	else //そうでないのなら、lineIntersection()で交点を返しておく。
		return make_pair(true, lineIntersection(s1, s2));
	}
	return make_pair(false, error_val);
}

試し打ち問題
onlinejudge.u-aizu.ac.jp

distanceBetweenPointAndLine()　点と直線の距離

　点と直線の距離と言えば、ヘッセの公式を思い浮かぶ人は多いと思いますが、ここも代数学的ではなく、幾何学的にアプローチしましょう。

　外積で平行四辺形の面積を求めて、それの高さとして距離を求めるものです。[出典1]

クリックで展開

//点と直線の距離。引数は、点、直線上の2点
double distanceBetweenPointAndLine(const Point& p, const Line& l) {
	return abs(l.vec().cross(p - l.begin) / l.vec().length());
}

distanceBetweenPointAndRay()　点と半直線の距離

　これは直線バージョンと非常に似ています。イメージとしては、無限に伸びる方が最短距離ならば、点と直線の距離と同じ扱いです。違うのは、始点方向であり、直線と違って伸びてません。下図を参照してください。[出典1]

　図のタイトルが間違ってるので注意！

　コードでは直接内積で書いてますが、もちろん第二話で実装した、鋭角、直角、鈍角を判定するangletype()を利用しても構いません。

クリックで展開

//点と半直線の距離。引数は、点、半直線(始点から終点方向に延びる)
double distanceBetweenPointAndRay(const Point& p, const Ray& r) {
	//始点との距離のパターン
	if (sgn((p - r.begin).dot(r.vec())) < 0)
		return r.begin.distanceFrom(p);
	return abs(r.vec().cross(p - r.begin) / r.vec().length());
}

distanceBetweenPointAndSegment()　点と線分の距離

　これは、半直線と同じように考えられます。半直線では、始点となった端では角度判定をしましたが、線分では両方の端に対して同じように判定すればよいのです。すなわち

• 左の端より左にあるのなら、左端との距離が最短
• 右の端より右にあるのなら、右端との距離が最短
• それ以外ならば、直線とみなした時の距離が最短

です。[出典1]

クリックで展開

//点と線分の距離。引数は、点、線分の両端
double distanceBetweenPointAndSegment(const Point& p, const Segment& s) {
	if (sgn(s.vec().dot(p - s.begin)) < 0 || sgn(s.countervec().dot(p - s.end)) < 0) {
		//下した垂線は線分の上にはない
		return min(p.distanceFrom(s.begin), p.distanceFrom(s.end));
	}
	return distanceBetweenPointAndLine(p, s);
}

distanceBetweenSegmentAndSegment()　線分と線分の距離

1. • 直線と直線の距離: 平行なら1つの直線とのこりの直線上の任意の1点との距離、それ以外なら0
2. • 直線と線分の距離: 交差していたら0、それ以外なら線分の両端と直線との距離のうち小さいほう
3. • 線分と線分の距離: 交差していたら0、それ以外な

ら4通りの線分と点の距離のmin

https://www.ioi-jp.org/camp/2017/2017-sp_camp-hide.pdf

　すでに何度も参考にさせていただいてますが、このスライドではこう書かれています。証明は難しくないので皆さん考えましょう(丸投げ)。これを実装します。

クリックで展開

//二線分間の距離
double distanceBetweenSegmentAndSegment(const Segment& s1, const Segment& s2) {
	if (segmentIntersection(s1, s2).first)return 0;//交点を持つ
	double ans = distanceBetweenPointAndSegment(s1.begin, s2);
	ans = min(ans, distanceBetweenPointAndSegment(s1.end, s2));
	ans = min(ans, distanceBetweenPointAndSegment(s2.begin, s1));
	ans = min(ans, distanceBetweenPointAndSegment(s2.end, s1));
	return ans;
}

試し打ち問題
onlinejudge.u-aizu.ac.jp

projection()　正射影

　第一話の内積のところで、少し紹介しました、点の正射影についてです。イメージとしては、ある直線に垂直に光が当たるとき、点が直線のどこで影を作るところを求めることです。

　下図の点Aが直線BCへの正射影は点Hです。これは直線を垂線を下ろしてその足ということでもあります。求め方も以下の通りです。

クリックで展開

//正射影
//引数は点A, B, Cで、Aの直線BC上の正射影を求める。
Point projection(const Point& a, const Line& l) {
	Point ret;
	ret = l.begin +
		l.vec().normalized() * (a - l.begin).dot(l.vec()) / l.vec().length();
	return ret;
}

試し打ち問題
onlinejudge.u-aizu.ac.jp

reflection()　線対称(鏡映)

　正射影と非常に似てる概念として、線対称があります。

　正射影を利用してこれを計算することができますが、以下のようにベクトルのたどる経路について気を付ける必要があります。例えば、試し打ち問題として紹介するAOJの問題では、赤い経路をたどって計算すると、誤差によってWAになります。

クリックで展開

//鏡映変換。引数は点A, B, Cで、直線BCにおいて、Aと線対称な点を求める。
Point reflection(const Point& a, const Line& l) {
	Point ret;
	ret = a + 2 * (projection(a, l) - a);
	return ret;
}

試し打ち問題
onlinejudge.u-aizu.ac.jp

今回の全容

クリックで展開

//向きつき直線を想定。Segment
struct Line {

	//直線の通る二点。有向直線でないのならば、beginやendに違いはない。
	Point begin, end;

	Line() {
		begin = Point(), end = Point();
	}

	Line(const Point& b, const Point& e) {
		begin = b, end = e;
	}

	//ax+by+c=0
	Line(const double a, const double b, const double c) {
		if (sgn(a) == 0 && sgn(b) == 0) {
			assert(-1);
		}

		if (sgn(b) == 0) {
			//ax+c=0になる。
			begin = Point(-c / a, 0.0);
			end = Point(-c / a, 1.0);
		}
		else {
			//y=-(ax+c)/b　傾きは-a/bで、y切片が-c/b
			begin = Point(0, -c / b);
			end = Point(1.0, -(a + c) / b);
		}
	}

	Point vec() const {
		return end - begin;
	}

	Point countervec() const {
		return begin - end;
	}

};

//半直線なら、Lineと同じだけど、はっきりと区別する。
typedef Line Ray;
//線分の場合、Lineと同じものになるが、はっきりと区別する
typedef Line Segment;


//直線の交点を返す。交わってなければ、error_valを返す。
Point lineIntersection(const Line& l1, const Line& l2) {
	if (sgn(l1.vec().cross(l2.vec())) == 0)return error_val;

	Point ret;
	ret = l1.begin + l1.vec() *
		abs((l2.end - l1.begin).cross(l2.vec()) / l1.vec().cross(l2.vec()));
	return ret;
}

//線分が共通部分を持つかどうか？と線分の交点を返す。共通部分がない、もしくは交点が一意ではないなら、error_valを返す。
//trueなら、共通部分を持つ。falseなら、共通部分を持たない。
pair<bool, Point> segmentIntersection(const Segment& s1, const Segment& s2) {

	if (iSP(s1.begin, s1.end, s2.begin) * iSP(s1.begin, s1.end, s2.end) <= 0 &&
		iSP(s2.begin, s2.end, s1.begin) * iSP(s2.begin, s2.end, s1.end) <= 0) {

		//平行ならば、交点は定まらない。(完全に重なってるので)
		if (s1.vec().cross(s2.vec()) == 0)
			return make_pair(true, error_val);
		else //そうでないのなら、lineIntersection()で交点を返しておく。
			return make_pair(true, lineIntersection(s1, s2));

	}
	return make_pair(false, error_val);
}


//点と直線の距離。引数は、点、直線上の2点
double distanceBetweenPointAndLine(const Point& p, const Line& l) {
	return abs(l.vec().cross(p - l.begin) / l.vec().length());
}

//点と半直線の距離。引数は、点、半直線(始点から終点方向に延びる)
double distanceBetweenPointAndRay(const Point& p, const Ray& r) {
	//始点との距離のパターン
	if (sgn((p - r.begin).dot(r.vec())) < 0)
		return r.begin.distanceFrom(p);
	return abs(r.vec().cross(p - r.begin) / r.vec().length());
}

//点と線分の距離。引数は、点、線分の両端
double distanceBetweenPointAndSegment(const Point& p, const Segment& s) {
	if (sgn(s.vec().dot(p - s.begin)) < 0 || sgn(s.countervec().dot(p - s.end)) < 0) {
		//下した垂線は線分の上にはない
		return min(p.distanceFrom(s.begin), p.distanceFrom(s.end));
	}
	return distanceBetweenPointAndLine(p, s);
}

//二線分間の距離
double distanceBetweenSegmentAndSegment(const Segment& s1, const Segment& s2) {
	if (segmentIntersection(s1, s2).first)return 0;//交点を持つ

	double ans = distanceBetweenPointAndSegment(s1.begin, s2);
	ans = min(ans, distanceBetweenPointAndSegment(s1.end, s2));
	ans = min(ans, distanceBetweenPointAndSegment(s2.begin, s1));
	ans = min(ans, distanceBetweenPointAndSegment(s2.end, s1));
	return ans;
}

//正射影
//引数は点A, B, Cで、Aの直線BC上の正射影を求める。
Point projection(const Point& a, const Line& l) {
	Point ret;
	ret = l.begin +
		l.vec().normalized() * (a - l.begin).dot(l.vec()) / l.vec().length();
	return ret;
}

//鏡映変換。引数は点A, B, Cで、直線BCにおいて、Aと線対称な点を求める。
Point reflection(const Point& a, const Line& l) {
	Point ret;
	ret = a + 2 * (projection(a, l) - a);
	return ret;
}

最後に

　直線、半直線、線分については以上です。計算するにあたって、幾何的なアプローチを中心にして、誤差やコーディングの煩雑さを減らした実装になっています。直線関連の判定はシンプルながらも、多くの幾何の問題の幾何要素と直結していて、非常に重要だと思います。それぞれの機能の導出する方法を覚えて、ソラでも書けるのを目指しましょう。僕も頑張ります。

　次回は、円についての計算幾何になります。

　重ね重ね、[出典1]のPDFを作ってくれた秀郁未さん、校正に手を貸してくれた早稲田大学競技プログラミング好きの皆さん、ありがとうございました。

参考文献

[出典1]
geometry
https://www.ioi-jp.org/camp/2017/2017-sp_camp-hide.pdf
このようなPDFを作ってくれた方に改めて大いに感謝します！

全体的に参考になった文献

おりたたみ機能の書き方
bluebirdofoz.hatenablog.com