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

幾何ライブラリの説明をしていきます。
ここで扱うものはAOJの計算幾何学に入っている問題です。
参照している問題に書いてある変数等の宣言はしないので、URLが書いてあるものについては一度見てください。
おそらく幾つかに分けることになると思います。
高校数学を前提知識としていますのでその部分の説明はあまりしません。
また、include文などを逐次記述していると膨れてしまうので、記述しません。
最低限理解できればいいと思って書いているので、式変形や変数の定義など書かない時が多々あります。

ここでは、一度記述したものについては再度記述することなく使用しますので、標準のヘッダから提供される関数等と混同しないように注意してください。

分かりにくい場合や説明が足りない場合、コメントを残してもらえれば、追記するかもしれません。

高校数学範囲外の話
テンプレート
- マクロ
- 構造体とそれに伴う基本演算
時計・反時計回り
交差判定

高校数学範囲外の話

・ $asin(x)$ , $acos(x)$ , $atan(x)$
下記のサイトの説明がすごく分かりやすいです。
$Arctan x$ についてのみ $atan2(y, x)$ という関数があります。
この関数は対辺と隣辺の比ではなく符号付の長さを引数にすることで定義域を広くしています。
定義域は $[-\pi, \pi]$ となっています。
mathtrain.jp

・外積
二次元幾何以外を扱わない為、二次元幾何における外積の演算内容と簡単な意味のみ説明します。
$\vec{a} = (a_x, a_y)$ , $\vec{b} = (b_x, b_y)$ とします。
この時外積 $\vec{a} \times \vec{b}$ は、
$\vec{a} \times \vec{b} = a_xb_y - a_yb_x = \|\vec{a}\|\|\vec{b}\|sin\theta$
となります。プログラム上では二項目を使いますが、三項目も説明でよく使います。
外積の値は $\vec{a}$ , $\vec{b}$ によって作られる平行四辺形の面積となっています。これは三項目より明らかです。

テンプレート

共通して使うものを幾つか並べます。
コメントアウトで軽い説明をしています。

マクロ

const double EPS = 1e-10;
const double PI = acos(-1);
#define equals(a, b) (fabs((a) - (b)) < EPS) // 誤差を考慮した同値判定

構造体とそれに伴う基本演算

vector(動的配列) と Vector(数学のベクトル) 混乱しそうだな...

struct Point { // 点の構造体
	double x, y; // (x, y)

	Point() {}
	Point(double x, double y) :x(x), y(y) {}

	Point operator+(const Point& p) const { return Point(x + p.x, y + p.y); }
	Point operator-(const Point& p) const { return Point(x - p.x, y - p.y); }
	Point operator*(const double& k) const { return Point(x * k, y * k); }
	Point operator/(const double& k) const { return Point(x / k, y / k); }

	friend istream& operator>>(istream& is, Point& p) {
		is >> p.x >> p.y;
		return is;
	}

	bool operator==(const Point& p) const { return (fabs(x - p.x) < EPS && fabs(y - p.y) < EPS); }
	bool operator<(const Point& p) const { return (x != p.x ? x < p.x : y < p.y); }

	double norm() { return x * x + y * y; }
	double abs() { return sqrt(norm()); }
};

typedef Point Vector; // 点とベクトルは同じように表せる

double norm(Vector a) { return a.x * a.x + a.y * a.y; }
double abs(Vector a) { return sqrt(norm(a)); }
double dot(Vector a, Vector b) { return a.x * b.x + a.y * b.y; } // 内積
double cross(Vector a, Vector b) { return a.x * b.y - a.y * b.x; } // 外積

double arg(Vector p) { return atan2(p.y, p.x); } // 偏角
Point polar(double r, double a) { return Point(cos(a) * r, sin(a) * r); } // 極座標から直交座標への変換

bool isParallel(Vector a, Vector b) { return equals(cross(a, b), 0.0); } // 平行判定
bool isOrthogonal(Vector a, Vector b) { return equals(dot(a, b), 0.0); } // 直交判定

struct EndPoint { // 端点の構造体(マンハッタン幾何)
	Point p;
	int seg, st; // 点のid, 点の種別

	EndPoint() {}
	EndPoint(Point p, int seg, int st) :p(p), seg(seg), st(st) {}

	bool operator<(const EndPoint& ep) const {
		if (p.y == ep.p.y) return st < ep.st;
		return p.y < ep.p.y;
	}
};

struct Segment { // 線分の構造体
	Point p1, p2; // 2端点

	Segment() {}
	Segment(Point p1, Point p2) :p1(p1), p2(p2) {}

	friend istream& operator>>(istream& is, Segment& s) {
		is >> s.p1 >> s.p2;
		return is;
	}
};

typedef Segment Line; // 直線は線分に長さを無くしただけのもの

Point project(Segment s, Point p) { // 点の射影
	Vector base = s.p2 - s.p1;
	double r = dot(p - s.p1, base) / base.norm();
	return s.p1 + base * r;
}

Point reflect(Segment s, Point p) { // 点の反射
	return p + (project(s, p) - p) * 2.0;
}

struct Circle { // 円の構造体
	Point c; // 中心点
	double r; // 半径

	Circle() {}
	Circle(Point c, double r) :c(c), r(r) {}
};

typedef vector<Point> Polygon; // 多角形(順序付きの点の集合)

時計・反時計回り

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/1/CGL_1_C

まず、
$\vec{a} = \vec{p_1} -\vec{p_0}$ , $\vec{b} = \vec{p_2} - \vec{p_0}$ ,
$\vec{a}$ , $\vec{b}$ のなす角を $\theta$ とします。
この時、外積 $\vec{a} \times \vec{b} = \|\vec{a}\|\|\vec{b}\|sin\theta$ より時計・反時計回りの判定ができます。
外積の値が正であれば、 $sin\theta$ が正であるので反時計回りと判定できます。
同様に、外積の値が負であれば、時計回りと判定できます。
これ以降どちらでもない場合、つまり外積が $0$ の時について考えます。
この場合、三点が一直線上にあることが分かるのでもう少し詳しく分類していきます。
まず、内積 $\vec{a} \cdot \vec{b} = \|\vec{a}\|\|\vec{b}\|cos\theta$ が負の時、 $\theta = \pi$ であるので ONLINE_BACK となります。
そうでないときは $\theta = 0$ です。
この時 ONLINE_FRONT となるには $\|\vec{a}\| < \|\vec{b}\|$ である必要があり、そうでなければ ON_SEGMENT となります。

static const int COUNTER_CLOCKWISE = 1;
static const int CLOCKWISE = -1;
static const int ONLINE_BACK = 2;
static const int ONLINE_FRONT = -2;
static const int ON_SEGMENT = 0;

int ccw(Point p0, Point p1, Point p2) {
	Vector a = p1 - p0, b = p2 - p0;
	if (cross(a, b) > EPS) return COUNTER_CLOCKWISE;
	if (cross(a, b) < -EPS) return CLOCKWISE;
	if (dot(a, b) < -EPS) return ONLINE_BACK;
	if (a.norm() < b.norm()) return ONLINE_FRONT;
	return ON_SEGMENT;
}

交差判定

線分同士の交差判定

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/2/CGL_2_B

線分 $s_1$ を直線とみなすと、2つの領域に分けることができます。このとき、点 $p_2$ , $p_3$ がこの2つの領域の別側に属している、つまりその2点が同じ領域に属していないことが必要条件となります。ただし、 $p_2$ , $p_3$ のうちどちらかが直線上にあるものも条件に含まれます。
同様に、線分 $s_2$ と点 $p_0$ , $p_1$ でも同じことが言えます。
この二つの条件を満たしているとき線分が交差していると判定できます。

bool intersectSS(Point p1, Point p2, Point p3, Point p4) {
	return (ccw(p1, p2, p3) * ccw(p1, p2, p4) <= 0 &&
			ccw(p3, p4, p1) * ccw(p3, p4, p2) <= 0);
}
bool intersectSS(Segment s1, Segment s2) {
	return intersectSS(s1.p1, s1.p2, s2.p1, s2.p2);
}

円と線分の交差判定

線分が円に接する場合は交点が2つあると考えることにします。
まず交点を持たない条件について考えます。
線分を直線とみなしたとき、その直線と円の中心点との距離が円の半径より大きい場合確実に交点を持ちません。これは射影より求められます。この射影で求めた点を $C$ します。
また、円の内部に線分が内包されるときも勿論交点を持ちませんこれは線分の両端点と円の中心点との距離が円の半径より小さい場合と言い換えることができます。
上記の条件だけでは十分ではありません、円の外に線分が位置し、最初の条件を満たさない場合が存在するからです。
しかし、これ以上考えるのは難しいので、上記の条件を満たさなければ、交点を1つもつ条件と2つもつ条件を考え、そのどちらにも属さなかった場合に交点を持たないと判断することにします。
したがって、まず交点を1つもつ条件を考えます。
線分のどちらかの端点が円の内部に、もう片方の端点が円の外部にあるとき、交点をただ1つもつと言えます。これは線分の端点と円の中心点との距離を求めることで分かります。
次に交点を2つもつ条件を考えます。
上記で求めた点 $C$ が線分上にある時に交点を2つもつといえます。
よって上記のどれにも当てはまらなかったときには交点をもちません。

int intersectCS(Circle c, Segment s) {
	if (norm(project(s, c.c) - c.c) - c.r * c.r > EPS) return 0;
	double d1 = abs(c.c - s.p1), d2 = abs(c.c - s.p2);
	if (d1 < c.r + EPS && d2 < c.r + EPS) return 0;
	if ((d1 < c.r - EPS && d2 > c.r + EPS) || (d1 > c.r + EPS && d2 < c.r - EPS)) return 1;
	Point h = project(s, c.c);
	if (dot(s.p1 - h, s.p2 - h) < 0) return 2;
	return 0;
}

円同士の交差判定

https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/7/CGL_7_A

まんま高校数学なので省略します。
戻り値は共通接線の数です。

int intersectCC(Circle c1, Circle c2) {
	if (c1.r < c2.r) swap(c1, c2);
	double d = abs(c1.c - c2.c);
	double r = c1.r + c2.r;
	if (equals(d, r)) return 3; // 外接する
	if (d > r) return 4; // 離れている
	if (equals(d + c2.r, c1.r)) return 1; // 内接する
	if (d + c2.r < c1.r) return 0; // 内包する
	return 2; // 交わる
}