【学习笔记】计算几何

常用函数

一些基本定义

为了简化代码，我打出了

# define ConP const Point&

防止精度误差而常见 $EPS=10^{-9}$（注：为了缩短篇幅，部分代码不合语法，比如大括号的省略）。使用符号函数 getSgn，正数、负数、零分别对应 $+1,-1,0$ 。

const double EPS = 1e-9;
inline int getSgn(const double &x)
	return x < -EPS ? -1 : x > EPS;

那么判断大小就直接使用 getSgn 了呗！进而有比较函数 cmp 。比较像 $\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}$ 。

/** @return if x > y return 1 */
inline int cmp(const double &x,const double &y)
	return getSgn(x-y); // easy access to getSgn

然后是 Point，点（或者向量）的类。它是一切的基础了——没有向量，计算几何怎么算啊？

struct Point{
	double x, y;
	Point(double X,double Y):x(X),y(Y)
	Point operator + (const Point &t)
		return Point(x + t.x, y + t.y);
	Point operator - (const Point &t)
		return Point(x - t.x, y - t.y);
	Point operator * (const double &t)
		return Point(x * t, y * t);
	Point operator / (const double &t)
		return Point(x / t, y / t);

可以再重载一下，以 $x$ 为第一关键字、$y$ 为第二关键字的 < 。

bool operator < (const Point &p) const {
	int c = cmp(x,p.x);
	if(c) return c == -1; // x < p.x
	return cmp(y,p.y) == -1; // y < p.y
}
bool operator == (const Point &p)
	return !cmp(x,p.x) && !cmp(y,p.y);

然后是点积与叉乘。不建议重载为运算符。写成 Point 的成员函数足矣。

double dot(const Point &p) return x*p.x+y*p.y;
double det(const Point &p) return x*p.y-y*p.x;

还有长度。

double len2() return x*x+y*y; // or dot(*this)
double len() return sqrt(len2());

还有其他一些常用的，都定义一下。一般的题目里不会用这么多。

小常识：cmath::atan2(y,x) 将返回一个 $\theta \in (-\pi ,\pi ]$，使得以 $x$ 轴正半轴为始边、$\theta$ 为其弧度的角的终边过 $(x,y)$ 。更具体地说：返回的是 极角。

三角形面积

说海伦公式的拖出去喂狗。给出三个顶点，求 有向面积，就是裸叉乘。但是为什么不调用 det 函数呢？我也不知道 😕

还有我们熟知的，叉乘判断左右。就是上面那个面积的正负呗。注意这个 triangle 计算的是平行四边形面积。

# define triangle(p1,p2,p3) ((p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y))
# define leftright(p1,p2,p3) getSgn(triangle(p1,p2,p3))

判断直线平行

说斜率的拖出去喂狗。我们的核心就是：避免除法、避免精度误差。

还是用叉乘。由于 叉乘有分配率，求 $(\overrightarrow{p_2}-\overrightarrow{p_1})\times (\overrightarrow{q_2}-\overrightarrow{q_1})$ 可以把右边的括号拆开，得到

bool isParallel(ConP p1,ConP p2,ConP q1,ConP q2){
	double a = triangle(p1,p2,q1);
	double b = triangle(p1,p2,q2);
	return getSgn(a-b) == 0;
}

其实本质上就是平行线间形成的三角形面积相等。

求直线交点

众所周知，一种表示共线向量的方法：设交点为 $M$，如果 $|\overrightarrow{P_{1}M}|:|\overrightarrow{P_{2}M}|=a:b$，那么 $\overrightarrow{OM}=\frac{b\overrightarrow{O{P}_1}+a\overrightarrow{O{P}_2}}{a+b}$ 。所以现在的问题还是求出二者长度的比值。

有一个巧妙的转化：三角形的底相同时，高的 长度之比等于面积之比，所以有
$$\frac{|\overrightarrow{P_{1}M}|}{|\overrightarrow{P_{2}M}|}=\frac{S_{△Q_1{P}_1{Q}_2}}{S_{△Q_1{P}_2{Q}_2}}$$

就又成为了三角形面积计算的裸题。

Point intersectLL(ConP p1,ConP p2,ConP q1,ConP q2){
	double S1 = triangle(q1, q2, p1);
	double S2 = -triangle(q1, q2, p2);
	return (p1 * S2 + p2 * S1) / (S1 + S2);
}

点是否在线段上

其实有一个非常简单的做法：看 $l_{OP}$ 与 $l_{AB}$ 的交点是不是 $P$ 。但是这样太暴力了……

不如老老实实的，它等价于 $\overrightarrow{AB}$ 和 $\overrightarrow{PB}$ 共线，且 $P$ 的 $x$ 坐标在 $[x_A,x_B]$ 之内。

bool inRange(double l,double x,double r)
	return cmp(x,l) != cmp(x,r);
bool onSegment(ConP p1,ConP p2,ConP q)
	return !leftright(p1,p2,q) && inRange(p1.x,q.x,p2.x);

当然你仔细想想，$cmp$ 不同也就是相乘为负，而相乘的结果就是点积。

计算几何

求多边形面积

这里有一个非常有用的思想，三角剖分。我们只会求三角形的面积，那就把多边形划分成三角形呗！但是考虑到多边形可能非常畸形，三角剖分不容易进行。

事实上，我们只需要给三角形定一个系数 $\pm 1$ 即可。因为 多边形内任意一点引出的射线与多边形相交奇数次（不考虑在端点上相交的情况）。这是拓扑学基本原理了，相交必然改变内外情况。同理，多边形外的点引出射线就会相交偶数次。

那么，如果我们利用叉乘自带的正负号，容易发现，从左边跃到右边会减，从右边跳到左边会加。那么，相交偶数次的，恰好就正负抵消了；相交奇数次的，我们 不妨让第一次是右至左，那么它恰好就是正的。怎么才能让第一次是右至左呢？逆时针遍历每条边 就行了。

于是我们写出了一个非常简单的代码，却可以处理所有长相怪异的多边形！

double area(const vector<Point> &ps){
	double res = 0; int len = ps.size();
	for(int i=0; i!=len; ++i)
		res += ps[i].det(ps[(i+1)%len]);
	return res/2; // det = 2*S_{triangle}
} // 你也可以对 res 求绝对值，就不需要 ps 逆时针顺序了

注意：边不能自交。如果边可以自交，那么偶数条边不一定是一半左到右、一半右到左。

点是否被包含

完全就是上面那玩意儿的应用。取一个朝向 $x$ 轴负半轴的射线，判断交点个数就可以了。这是 $\mathcal{O}(n)$ 的。

/** @return 2:inside  1:onSegment  0:outside */
short contain(const vector<Point> &ps,ConP q){
	int len = ps.size(), res = 0;
	for(int i=0; i!=len; ++i){
		Point u = ps[i], v = ps[(i+1)%len];
		if(onSegment(u,v,q)) return 1;
		if(cmp(u.y,v.y) <= 0) swap(u,v);
		if(cmp(q.y,u.y) > 0 || cmp(q.y,v.y) != 1)
			continue; // not in y range
		res ^= (leftright(q,u,v) > 0);
	}
	return res<<1;
}

如果已知多边形是一个凸包，取射线为垂直于 $x$ 轴的。可直接 lower_bound 找到 $x$ 恰比这个点大的 $P_0$，取它的前驱节点 $Q_0$ 必定满足其 $x$ 小于这个点，我们就找到了一条线段与该射线相交。再看看这个线段是上半部分还是下半部分，当前点是这个线段的上方还是下方，就可以得到结果。这是 $\mathcal{O}(\log n)$ 的。

凸包与最近点对

这都是经典问题了，并不需要计算几何的知识。才怪呢。

前者可以就用斜率去理解它。后者是 $\mathcal{O}(n\log n)$ 经典分治。另：当 $n\le 10$ 时，建议暴力求解最近点对。因为老师的课件里，这是递归出口。由于题目中遇到了，我顺便贴个代码，仅供参考。

Point ddg[MaxN];
bool cmpAngle(const Point &a,const Point &b){
	int sgn = getSgn(a.det(b));
	if(sgn) return sgn == 1;
	return getSgn(a.len2()-b.len2()) < 0;
}
int convexHull(){
	Point o = p[1]; rep(i,2,n)
		if(p[i].y == o.y && p[i].x < o.x)
			o = p[i]; // most left
		else if(p[i].y < o.y) o = p[i];
	rep(i,1,n) ddg[i] = p[i]-o;
	sort(ddg+1,ddg+n+1,cmpAngle);
	int tot = 2; // how many points
	for(int i=3; i<=n; ++i){
		while(tot >= 2)
			if(getSgn(triangle(ddg[tot-1],
				ddg[tot],ddg[i])) != 1)
					-- tot; // pop_back
			else break; // nothing to do
		++ tot; ddg[tot] = ddg[i];
	}
	rep(i,1,tot) ddg[i] = ddg[i]+o;
	return tot;
}

这是 $\mathrm{G}\mathrm{r}\mathrm{a}\mathrm{h}\mathrm{a}\mathrm{m}$ 扫描法。这招不需要选择最左下角的点，只需要确保你选择的点在凸包内。听上去似乎随便选一个点不就行了？然而原点的极角该怎么规定呢？所以最好取三个点的重心。

比如这道动态凸包，代码逻辑就非常简单了。提交记录在此。结果似乎被上下凸包分开算的方法吊打了。

话说 $\mathrm{M}\mathrm{e}\mathrm{l}\mathrm{k}\mathrm{m}\mathrm{a}\mathrm{n}$ 求凸包 一直没学懂，先贴个链接吧，希望以后自己能够研究清楚……

点集直径

显然最远的点只会在凸包上。然后就是旋转卡壳了。这个我就不讲了，请自己找资料吧。

double convexDiameter(const vector<Point> &ps){
	int len = ps.size(); if(len <= 1) return 0;
	int is = 0, js = 0; // init
	for(int i=0; i!=len; ++i){
		if(ps[i] < ps[is]) is = i; // left
		if(ps[js] < ps[i]) js = i; // right
	}
	double res = ps[is].distTo(ps[js]);
	for(int i=is,j=js; true; ){
		Point ii = ps[(i+1)%len]-ps[i];
		Point jj = ps[(j+1)%len]-ps[j];
		if(ii.det(jj) >= 0) j = (j+1)%len;
		else i = (i+1)%len; // move ahead
		res = max(res,ps[i].distTo(ps[j]));
		if(i == is && j == js) break ;
	}
	return res;
}

圆与线的交点

设圆心为 $O$，直线上两点为 $P,Q$，圆的半径为 $r$，求解直线与圆的交点 $A,B$（其中 $A$ 是更靠近 $P$ 的那一个）。

显然我们可以联立方程然后解交点。但是它的精度就令人堪忧了……而且 $\mathrm{s}\mathrm{q}\mathrm{r}\mathrm{t}$ 不是一个非常快的函数……

我们有强大的工具，点积。用它往往可以求出投影向量。记 $O$ 与过 $P,Q$ 的直线 $l$ 的垂足为 $M$，显然
$$\overrightarrow{PM}=\frac{\overrightarrow{PO}\cdot \overrightarrow{PQ}}{|\overrightarrow{PQ}{|}^2}\times \overrightarrow{PQ}$$

投影向量的 “定义” 嘛。那我们只需要求出 $|\overrightarrow{MB}|$ 就可以确定 $A,B$ 的位置了。

有一个巧妙的方法：注意到 $\{\begin{array}{l}{r}^2=|\overrightarrow{MB}{|}^2+|\overrightarrow{OM}{|}^2\\ |\overrightarrow{OP}{|}^2=|\overrightarrow{MP}{|}^2+|\overrightarrow{OM}{|}^2\end{array}$，二者相减可得
$$|\overrightarrow{OP}{|}^2-r^2=|\overrightarrow{MP}{|}^2-|\overrightarrow{MB}{|}^2$$

而前面我们已经知道了
$$|\overrightarrow{MP}{|}^2={\left(\frac{\overrightarrow{PO}\cdot \overrightarrow{PQ}}{|\overrightarrow{PQ}|}\right)}^2$$

所以很轻松地写出
$$|\overrightarrow{MB}{|}^2\times |\overrightarrow{PQ}{|}^2=(\overrightarrow{PO}\cdot \overrightarrow{PQ}{)}^2-|\overrightarrow{PQ}{|}^2\left(|\overrightarrow{OP}{|}^2-r^2\right)$$

为啥要把 $|\overrightarrow{PQ}{|}^2$ 乘过来呢？因为求 $\overrightarrow{PM}$ 的时候就求出了 $\overrightarrow{PO}\cdot \overrightarrow{PQ}$ 的值，现在就可以二次利用啦！同时我们也要避免对于很小的数求 $\mathrm{s}\mathrm{q}\mathrm{r}\mathrm{t}$，否则精度误差大。在这个式子里，我们可以先把右侧的值算出来，然后开根得到 $|\overrightarrow{MB}|\times |\overrightarrow{PQ}|$，除以 $|\overrightarrow{PQ}{|}^2$ 后作为 $\overrightarrow{PQ}$ 的系数就是 $\overrightarrow{MB}$ 啦！只做了一次平方根！

pair<Point,Point> intersectCL(ConP o,double r,ConP p,ConP q){
	double x = (o-p).det(q-p), y = (q-p).len2();
	double d = x*x-y*((p-o).len2()-r*r);
	Point m = p+(q-p)*(x/y);
	if(cmp((m-o).len2(),r*r) > 0){
		// They have no intersection
		// I don't know what to do ...
	}
	Point dr = (q-p)*(sqrt(d)/y);
	return make_pair(m-dr,m+dr);
}

圆与圆的交

如果只求交集的面积，直接用余弦定理计算夹角就行。似乎 $\mathrm{s}\mathrm{q}\mathrm{r}\mathrm{t}$ 比三角函数更慢……

double area(Point A,double ra,Point B,double rb){
	if(getSgn(ra-rb) < 0) swap(A,B), swap(ra,rb);
	double dis = (B-A).len(), ans;
	if(getSgn(dis-rb-ra) >= 0) return 0; // separated
	if(getSgn(dis+rb-ra) <= 0) return M_PI*rb*rb;
	double xez = (B-A).len2(); // more precision
	double w = acos((ra*ra+xez-rb*rb)/2/ra/dis);
	ans = ra*ra*w-dis*ra*sin(w);
	w = acos((rb*rb+xez-ra*ra)/2/rb/dis);
	ans += rb*rb*w; return ans;
}

圆与三角形的交

事实上只能处理一种特殊的三角形：其中一个顶点与圆心重合。此时只需要大讨论一下各种情况。一般三角形见下方 areaC_（圆与多边形的交）。

由于代码很复杂，这里不给出代码实现，只需要知道这是可以实现就行了。

圆与多边形的交

直接三角剖分。因为三角剖分的本质是考虑每一个点，即使和 $\odot O$ 取了交集也是一样的。

然后就转化为了 areaCT，即圆（$circle$）和三角形（$triangle$）的交。

三角形重心

众所周知，$G$ 是 $△ABC$ 的重心，当且仅当
$$\overrightarrow{GA}+\overrightarrow{GB}+\overrightarrow{GC}=\overrightarrow{0}$$

所以 $G_x=\frac{1}{3}(A_x+B_x+C_x),G_y=\frac{1}{3}(A_y+B_y+C_y)$ 。虽然直接看这个式子也挺直观的。

Point gravityCenter(ConP A,ConP B,ConP C)
	return Point((A.x+B.x+C.x)/3,(A.y+B.y+C.y)/3);

三角形外接圆

没想到很好的求法。但是只求半径会简单一些。根据正弦定理，$R=\frac{a}{2\sin A}$，又根据面积计算公式有 $S=\frac{1}{2}bc\sin A$，代入得 $R=\frac{abc}{4S}$ 。

如果一定要求出圆心呢？目前我只会联立直线解方程……

Point Point::rotate90() return Point(-y,x);
Point circumCenter(Point A,Point B,Point C){
	Point M = (A+B)/2, N = (B+C)/2;
	Point E1 = (B-M).rotate90()+M;
	Point E2 = (B-N).rotate90()+N;
	return intersectLL(M,E1,N,E2);
}

而且之前搞错了，把重心当成了外接圆圆心……我只能说：

半平面交

说人话：找一个区域，使得区域内的点在所有给出向量（有向直线）的左侧。

由于此区域很可能是凸包，所以我们会想到，按照极角排序后处理。用一个双端队列维护凸包（有点像 $\text{Melkman’s Algorithm}$ 的感觉），每次从两端删边。

可以只维护直线，直线的交点（凸包的端点）可以临时计算。判断依据为，交点在直线的右侧，则将队列端头的直线删去。对了，平行的直线请提前判断，只保留一条。

若最后的交集是一条线呢？一个简单的处理方法是，上面的 “交点在直线右侧” 需要是严格的，即允许交点在线上，这样会保留多条交于一点的直线；可以避免两条平行（但方向相反）的直线求交点。

深入应用

最小圆覆盖

即，求一个圆心 $O$ 与半径 $r$，使得 $\odot O$（含边界）包含了点集内所有点。最小化 $r$ 的值。

这是一个神奇的题目。一般来说，圆上有三个点——若只有两个点，且不是圆的直径，可以使圆心略微向这条线段偏移使得半径变小。

我们可以首先打出一个 调整法（或称之为 增量法）的代码：顺序枚举点，如果点不在目前得到的最小圆覆盖中，则将其强制设定为圆上一点，嵌套一个循环来顺序枚举点，最多三层。代码见下。

pair<Point,double> min_circle(vector<Point> ps){
	random_shuffle(ps.begin(),ps.end());
	int len = ps.size(); double r = 0;
	Point o = ps[0]; // whatever
	for(int i=0; i!=len; ++i){
		if(cmp(o.distTo(ps[i]),r) <= 0)
			continue; // contained
		o = ps[i], r = 0; // reset
		for(int j=0; j!=i; ++j){
			if(cmp(o.distTo(ps[j]),r) <= 0)
				continue; // involved
			o = (ps[i]+ps[j])/2; // middle
			r = o.distTo(ps[i]);
			for(int k=0; k!=j; ++k){
				if(cmp(o.distTo(ps[k]),r) <= 0)
					continue; // covered
				o = circumCenter(ps[i],ps[j],ps[k]);
				r = o.distTo(ps[i]); // update
			}
		}
	}
}

正确性的证明：原本这里只有两个字，“显然” 。然吾今日复视之，觉得它并没有我看上去那么简单。想了一下午，无果。遂祈神助。神曰：何难之有？今记之如下：

---

命题：求解问题 “找最小圆，使之过 $0/1/2$ 个特殊点，且包含点集内所有点” 时，如果点集为 $S(|S|>1)$ 的答案是 $\odot O$，且点 $P$ 在 $\odot O$ 外，则点集为 S + { P } S+\{P\} S+{P} 的答案 $\odot O^{\prime }$ 满足 $P$ 在 $\odot O^{\prime }$ 上。

证明：反证法。不妨设 $P$ 在 $\odot O^{\prime }$ 严格内部。根据定义，交集 $D=\odot O^{\prime }\cap \odot O\ne \varnothing$ 且 $S\subset D$ 。由于 $\odot O^{\prime }$ 上必然有至少两个 $S$ 中的点（否则可以调整至更小），所以 $\odot O^{\prime }$ 过 $D$，可知 $\odot O^{\prime }$ 与 $\odot O$ 有两个交点。不妨设其为 $A,B$，显然必须经过的特殊点是 { A , B } \{A,B\} {A,B} 的子集。

当 $D$ 的两侧都是劣弧时，考察以 $A,B$ 为直径的圆 $\odot C$ 。由于两侧均劣弧，故 $D\subset \odot C$，又因为过 $A,B$ 的圆中 $C$ 的直径最小，所以 $\odot C$ 比 $\odot O$ 更小。这说明 $\odot C$ 才是 $S$ 的答案，矛盾。

若 $D$ 的其中一侧是优弧，由于 $\odot O^{\prime }$ 比 $\odot O$ 大，再根据 $P\in \odot O^{\prime }$ 而 $P\notin \odot O$，可以粗略画一张图：

$P$ 一定落在 ${\ell }_{AB}$ 的优弧一侧，因为劣弧一侧总是 $\odot O$ 更凸。此时考虑 $△PAB$ 的外接圆 $\odot C$，弧 $APB$ 仍然是优弧，但比 $\odot O^{\prime }$ 上的优弧 $AB$ 更凹，可见 $\odot C$ 比 $\odot O^{\prime }$ 更小。同时显然有 $D\subset \odot C$，所以 $\odot C$ 才是 S + { P } S+\{P\} S+{P} 的答案，归谬。

---

还得分析一下复杂度。分析每一个循环的可能性——别忘了我们有随机打乱点。$i$ 的 if 语句成立，必要条件是 $i$ 是 $[1,i]$ 的最小圆覆盖的圆上点。我们可以认为 $\frac{3}{i}$ 是其概率。

$j$ 的 if 语句成立，必要条件也是 $j$ 为最小覆盖圆上点，概率是 $\frac{2}{j}$，而进入之后运行次数是 $\mathcal{O}(j)$，所以每个 $j$ 提供的期望贡献是 $\frac{2}{j}\cdot j=2$ 。

所以 $i$ 不变时，$j$ 循环一次是 $i\cdot 2=\mathcal{O}(i)$ 的。而 $i$ 的 if 被触发是 $\frac{3}{i}$ 的概率，故期望为 $\frac{3}{i}\cdot \mathcal{O}(i)=\mathcal{O}(1)$ 的。共 $n$ 个 $i$，总复杂度期望为 $\mathcal{O}(n)$ 的！

类似地，我们会发现，当维度为 $k$ 时，时间复杂度仅仅为 $\mathcal{O}(k!\cdot n)$ 而已！不可思议！

多边形重心

定义：重心 $G(x_0,y_0)$ 是满足
$$\underset{D}{\iint }(x-x_0)\mathrm{d}S=\underset{D}{\iint }(y-y_0)\mathrm{d}S=0$$

的点。其中 $S$ 是面积，$D$ 是平面上多边形内的点。

从物理学的意义上讲：重心是整个物体受重力的等效点。利用这一点，我们可以 $\mathcal{O}(n)$ 的求出它。

首先三角剖分——因为它的本质是考虑每个点，那么减去相当于是质量为 $-m$ 的一个点。然后对三角形求重心。然后两个三角形的重心用带权平均合并——可以想象成，找到使杠杆平衡的支点。