２．計算幾何のインフラ整備 シリーズ：「基礎的計算幾何ライブラリの作り方」

シリーズの構成

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

今回の内容

　今回は、計算幾何をするにあたって、必要なインフラ整備です。僕のメイン言語のC++で実装していきます。

　前回では、外積、内積、ベクトルなどといろいろな道具の説明をしましたが、残念ながらC++には2次元の直交座標やベクトルをそのまま扱えるフレームワークは存在しません。(終)
　
　というわけで、自分で作りましょう。

既存のフレームワークからの派生

　さきほど、そのまま扱えるフレームワークは存在しない！と言いましたけど、少し加工すれば使えるフレームワークは存在します。std::complex、複素数です。

　もともと、数3を履修した人ならばわかる通り、複素数(と複素平面)は2次元直交座標において、回転なども複素数の積で表せたりと、平面幾何と親和性が高いフレームワークです。そして、STLのすでに用意されてるstd::complexについては、後述する演算のほとんどが、すでにオーバーロードされてます。つまり、特殊なコーディングなしに、直感的な記述ができる、ということです。

　今回の講座では、一応一から作り直しますが、実際のコンテストでコピペができない場合は、std::complexを継承した、もしくはそのまま使っていくつかの機能(内積、外積)を定義するのが時間短縮になります。

　また、今回実装するPointという構造体は、座標(x, y)というものも、ベクトル(x, y)というのも意味します。(座標自体、(0,0)からそこへ伸びてるベクトルとみなせる)これを位置ベクトルと言います。

全体的な構造

　今回はライブラリという形なので、コピペして使用できる(競プロ以外の分野にも)ということで、ある程度の厳密性などのために、競プロでは主に使わないような構造になっています。主に

• ライブラリ全体はgeometry2dというネームスペース内にある。
• その中に、直線関連の機能はline2d、円関連の機能はcircle2d、多角形関連の機能はpolygon2dというネームスペースをそれぞれ設け、その中に格納する。
• 2次元座標(ベクトル)のクラスのほかにも、直線を表すLine、円を表すCircleなどを実装。

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

今回実装する機能

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

• EPS(浮動小数点精度)　2つの数の差がこれ以下ならば、一致するとみなす
• sgn()　与えられた数がEPSを加味して、正 or 負 or 0　を返す。
• 構造体　Point(これは自作する場合の話です)
• 入出力ストリームのオーバーロード
• 定数倍乗算の演算子(どの順番でもかけられるようにするため)
• '>'演算子　x座標を優先的に比較し、その次にy座標を比較する
• error_val　解なしなどになったときに返される値。
• ==演算子　(どの順番でも比較できるようにするため)
• iSP()　3点の位置関係を返す関数。
• angletype()　角度が鋭角なのか直角なのか鈍角なのかを返す(角度には符号はつかない)

具体的な実装(前編)

EPS, sgn()

　前の記事でも紹介したのと同じように、2つの数の差がある値以下ならば、それら2つを等しいものとみなす、というものです。その値をEPSとして、僕は定数として、としています。

　また、そのEPSを加味した判定関数として、sgn()も実装します。与えられた数が0より大きいか、小さいか、0と等しいかを返すので、a - b > 0という条件式は、sgn(a - b) > 0とそのまま書けます。

クリックで展開

	//eps　浮動小数点誤差
	const double EPS = 1e-10;

	//a > 0ならば+1, a == 0ならば0, a < 0ならば-1　を返す。　基本的にEPS込みの評価はこれで行う。
	//不等式は、加減算に直してこれに適用する。
	int sgn(const double a) {
		return (a < -EPS ? -1 : (a > EPS ? +1 : 0));
	}

Point型の構造体

　作り方はこの記事を大いに参考にさせていただきました。
qiita.com
　余談ですが、作者はSiv3DというWindowsアプリのフレームワークも制作してます。僕もそれなり使用したことがありますが、非常におすすめです！

　この記事の通りに実装しましょう(丸投げ)。ここでは、記事に存在しない昨日の実装や、変更があった実装にのみ説明します。ちなみに、内積はdot()で、外積はcross()です。(内積の記号はドット、外積は×なのでクロス)ちなみに、ほかの人の実装を見てる限り、
外積はdet()としてる人もいます。(detは行列式(determinant)の略。2次元版外積は2次行列式と値が同じなので)

これを実装することで、次のようなことができます。

• ベクトル同士の加減算
• ベクトルの要素に定数倍を掛け、割る
• ベクトルの長さ
• ある点から今の点までの距離
• 今の点の(ベクトルの)偏角
• 今の点の偏角を指定の度数だけ回す

.dot()や.cross()の試し打ち問題
onlinejudge.u-aizu.ac.jp

normalUnitVector()　正規法線ベクトル

　正規化されたベクトルというのは、長さが1のベクトルと言います。向きはしっかり持ってるので、伸ばしたい長さぶんだけ、x成分やy成分にかけてあげれば簡単にその向きの指定の長さを持つベクトルが作れます。

　ベクトルAの法線ベクトルというのは、Aと垂直なベクトルすべてです。(は？めちゃくちゃ一杯ありすぎるだろ・・・)単に垂直と言っても、2つの向きがありますが、ここではベクトルAの向きを反時計回りに90度回したものと定めます。逆の向きがほしいのならば、全部マイナスの符号をつければいいです。

　Aの正規法線ベクトルは、この2つを組み合わせたもので、「長さが1でベクトルAから反時計回りに90度回したベクトル」と僕は定めています。

　ちなみに法線ベクトルは英語でUnitVectorと言います。normalizeは一般化という意味ですが、めちゃくちゃ混ざりませんか？

クリックで展開

//単位法線ベクトル
Point normalUnitVector() const {
	return{ -normalized().y, normalized().x };
}

rotaion(), angle()　座標(ベクトル)の偏角、回転

クリックで展開

//(0, 0)中心にarg(弧度法)回転した座標
	Point rotation(double arg) const {
		double cs = cos(arg), sn = sin(arg);
		return Point(x * cs - y * sn, x * sn + y * cs);
	}
	//(0, 0)中心の円上に乗ってるとしたときの、偏角
	double angle() const {
		return atan2(y, x);
	}

具体的な実装(後編)

比較演算子

　ソートする際に、わざわざ第三引数に判定関数を渡さなくて済むようにあらかじめ定義しておきましょう。判定方法は、「まずx成分の大小で並び替え、x成分が同じものはy成分の大小で並び替え」というものです。

　elseのあとのreturnは書き忘れないようにしましょう！(2敗)青い予約語が並ぶと、書いたと錯覚して後々つらい目にあることがあるので、念には念を入れてチェックしよう！

クリックで展開

//第一にx, 第二にyを小さい順にソートするためのオーバーロード
	inline bool operator<(const Point& a, const Point& b) {
		if (sgn(a.x - b.x) != 0)return sgn(a.x - b.x) < 0;
		else return sgn(a.y - b.y) < 0;
	}

iSP()　3点の位置関係

　3点(順にA, B, Cとして)を入れて、その位置関係は、ABから見て

• BCは左に曲がってる
• BCは右に曲がってる
• ABC(CBA)の順に並んでいる
• ACB(BCA)の順に並んでいる
• BAC(BCA)の順に並んでいる

を返す関数です。

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

クリックで展開

/*
	3点A, B, Cの位置関係を返す関数 A, Bがすべて異なった点であるのが前提
	//ABから見てBCは左に曲がるのなら +1
	//ABから見てBCは右に曲がるのなら -1
	//ABC(CBA)の順番で一直線上に並ぶなら +2
	//ACB(BCA)の順番で一直線上に並ぶなら 0
	BAC(CAB)の順番で一直線上に並ぶなら -2
	*/
	int iSP(const Point& a, const Point& b, const Point& c) {
		int flg = sgn((b - a).cross(c - a));
		if (flg == 1) {
			return +1;
		}
		else if (flg == -1) {
			return -1;
		}
		else {
			//ABC(CBA)
			if (sgn((b - a).dot(c - b)) > 0)
				return +2;
			//BAC(CAB)
			else if (sgn((a - b).dot(c - a)) > 0)
				return -2;
			//ACB(BCA)　CがA or Bと一致しても、こっちに含まれる。
			else
				return 0;
		}
	}

angletype()　鋭角か？鈍角か？直角か？

　三点A, B, Cを与えて、角ABC(正の角度とみなす。もちろん180度より小さい方)が鋭角、鈍角、直角のいずれを返す関数。

クリックで展開

//角ABCが鋭角なら0、直角なら1、鈍角なら2を返す。
int angletype(const Point& a, const Point& b, const Point& c) {
	auto v = (a - b).dot(c - b);
	if (sgn(v) > 0)return 0;
	else if (sgn(v) == 0)return 1;
	else return 2;
}

std::complexでの流用(応用編)

　std::complexには、[出典1]のように、高速フーリエ変換で使われるときに、少し遅いと指摘されていますが、幾何の演算に関してはそんなギリギリを攻めることはないと思います。コーディング時間の短縮というメリットの方が大きいと思います。　

　そして、ここまで実装した機能の一部は、std::complexの以下の機能で流用可能です。

std::complex直下のメンバ関数

メンバ関数ですよ！

• x成分の代わりに、.real(), y成分の代わりに.imag()
• <演算子を除くすべての演算子は、オーバーロード済である。

非メンバ関数

• std::abs()　引数の絶対値を返す。これは実装した.length()や.distanceFrom()に相当します。
• std::arg()　偏角を返します。.angle()に相当します。std::norm()という似た機能も存在する。

現時点での完成品

クリックで展開

namespace geometry2d {

	//eps　浮動小数点誤差
	const double EPS = 1e-10;

	//a > 0ならば+1, a == 0ならば0, a < 0ならば-1　を返す。　基本的にEPS込みの評価はこれで行う。
	//不等式は、加減算に直してこれに適用する。
	int sgn(const double a) {
		return (a < -EPS ? -1 : (a > EPS ? +1 : 0));
	}

	struct Point {
		//2次元ベクトルのクラス

		double x, y;
		Point(double _x, double _y) {
			x = _x, y = _y;
		}
		Point() {
			x = 0, y = 0;
		}

		Point operator+() const {
			return *this;
		}
		Point operator-() const {
			return{ -x, -y };
		}
		Point operator+ (const Point& b) const {
			return{ x + b.x, y + b.y };
		}
		Point operator- (const Point& b) const {
			return{ x - b.x, y - b.y };
		}
		Point operator* (const double b) const {
			return{ x * b, y * b };
		}
		Point operator/ (const double b) const {
			return{ x / b, y / b };
		}
		Point operator+= (const Point& b) {
			x += b.x, y += b.y;
			return *this;
		}
		Point operator-= (const Point& b) {
			x -= b.x, y -= b.y;
			return *this;
		}
		Point operator*= (const double b) {
			x *= b, y *= b;
			return *this;
		}
		Point operator/= (const double b) {
			x /= b, y /= b;
			return *this;
		}
		bool operator== (const Point& b) {
			return b.x == x && b.y == y;
		}

		double lengthSquare() const {
			return (x * x + y * y);
		}
		double length() const {
			return std::sqrt(lengthSquare());
		}
		double dot(const Point& b) const {
			return x * b.x + y * b.y;
		}
		double cross(const Point& b) const {
			//Generally, cross product is vector, but in 2D, cross product is also scalar.
			return x * b.y - y * b.x;
		}
		double distanceFrom(const Point& b) const {
			return std::sqrt((x - b.x) * (x - b.x) + (y - b.y) * (y - b.y));
		}
		Point normalized() const {
			return{ x / length(), y / length() };
		}
		bool isZero() const {
			return sgn(x) == 0 && sgn(y) == 0;
		}
		//単位法線ベクトル
		Point normalUnitVector() const {
			return{ -normalized().y, normalized().x };
		}
		//(0, 0)中心にarg(弧度法)回転した座標
		Point rotation(double arg) const {
			double cs = cos(arg), sn = sin(arg);
			return Point(x * cs - y * sn, x * sn + y * cs);
		}
		//(0, 0)中心の円上に乗ってるとしたときの、偏角
		double angle() const {
			return atan2(y, x);
		}
	};

	inline Point operator*(double a, const Point& b) {
		return{ b.x * a, b.y * a };
	}

	template <class Char>
	inline std::basic_ostream<Char>& operator <<(std::basic_ostream<Char>& os, const Point& v)
	{
		return os << Char('(') << v.x << Char(',') << v.y << Char(')');
	}

	template <class Char>
	inline std::basic_istream<Char>& operator >> (std::basic_istream<Char>& is, Point& v)
	{
		return is >> v.x >> v.y;
	}

	//幾何ライブラリにおいて、解なしのような状態の時に返される点の値
	const Point error_val = { 114514.0, -191981.0 };

	inline bool operator==(const Point& a, const Point& b) {
		return (sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0);
	}

	inline bool operator!=(const Point& a, const Point& b) {
		return !(a == b);
	}

	//第一にx, 第二にyを小さい順にソートするためのオーバーロード
	inline bool operator<(const Point& a, const Point& b) {
		if (sgn(a.x - b.x) != 0)return sgn(a.x - b.x) < 0;
		else return sgn(a.y - b.y) < 0;
	}

	/*
	3点A, B, Cの位置関係を返す関数 A, Bがすべて異なった点であるのが前提
	//ABから見てBCは左に曲がるのなら +1
	//ABから見てBCは右に曲がるのなら -1
	//ABC(CBA)の順番で一直線上に並ぶなら +2
	//ACB(BCA)の順番で一直線上に並ぶなら 0
	BAC(CAB)の順番で一直線上に並ぶなら -2
	*/
	int iSP(const Point& a, const Point& b, const Point& c) {
		int flg = sgn((b - a).cross(c - a));
		if (flg == 1) {
			return +1;
		}
		else if (flg == -1) {
			return -1;
		}
		else {
			//ABC(CBA)
			if (sgn((b - a).dot(c - b)) > 0)
				return +2;
			//BAC(CAB)
			else if (sgn((a - b).dot(c - a)) > 0)
				return -2;
			//ACB(BCA)　CがA or Bと一致しても、こっちに含まれる。
			else
				return 0;
		}
	}

	//角ABCが鋭角なら0、直角なら1、鈍角なら2を返す。
	int angletype(const Point& a, const Point& b, const Point& c) {
		auto v = (a - b).dot(c - b);
		if (sgn(v) > 0)return 0;
		else if (sgn(v) == 0)return 1;
		else return 2;
	}

最後に

　インフラ整備は以上です。次回からは、直線にまつわる計算幾何について紹介し、実装します！

参考文献

[出典1]　

FFT高速化、一瞬で出来て効果が高いテクがあって、complexの掛け算をベタ打ち

— よすぽ (@yosupot) 2017年11月11日

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

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

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