[筆記] 計算幾何

IrisT發表於2022-05-02

筆記

其他的咕咕咕了，不喜歡計算幾何 qaq

基本定義

向量

\(\overrightarrow{PQ}=Q-P\)

叉積

二維 `a.x * b.y-a.y * b.x`

三維

`sum += a.x * b.y * c.z + a.y * b.z * c.x + a.z * b.x * c.y`

`sum -= a.z * b.y * c.x + a.y * b.x * c.z + a.x * b.z * c.y`

基本操作

將一個向量旋轉一定角度

Vec rotate(Vec a, LD b){
	LD s = sin(b), c = cos(b);
	return {a.x * c - a.y * s, a.x * s + a.y * c};
}

判斷一個點在直線的哪邊

有直線上的一點 \(P\)，直線的方向向量 \(v\)，想知道 \(Q\) 在直線哪邊：

利用叉積的性質，若 \(\overrightarrow{PQ}\times v >0\)，則 \(Q\) 在直線逆時針方向，否則在順時針方向。

bool Left(Vec a, Line b){ return cross(a - b.p , b.v) > 0; }

判斷兩圓之間的關係

返回的是公切線個數

int Pos(Circle a, Circle b){
	LD dis = Dis(a.O, b.O); if(a.r < b.r) swap(a, b);
	if(dis > a.r + b.r) return 4;
	if(dis == a.r + b.r) return 3;
	if(dis > a.r - b.r) return 2;
	if(dis == a.r - b.r) return 1;
	return 0;
}

求兩條直線的交點

有直線 \(AB,CD\)，求交點 \(E\)：

首先確定是否只有一個交點，然後因為記錄的是直線上的一個點和直線的方向向量，所以只需要知道這個點與交點的距離 \(l\) ，再將這個點沿方向向量平移 \(l\) 個單位長度即可。利用正弦定理求解

由上圖可知，\(|\mathbf{a\times b|=|a||b|\sin\beta,|u\times b|=|u||b|\sin\theta}\)。

作商得：

\[T=\frac{\mathbf{|u\times b|}}{\mathbf{|a\times b|}}=\frac{|\mathbf{u}|\sin\theta}{|\mathbf{a}|\sin\beta}=l \]

交點即點 \(B+T\mathbf{a}\)。

Vec Inter(Line a, Line b){ return a.p + cross(b.v, (b.p - a.p)) / cross(b.v, a.v) * a.v;}

求兩圓的交點

注意要先判斷有沒有交點

pair <Vec, Vec> Inter(Circle a, Circle b){
	LD x = a.r, y = b.r, z = Dis(a.O, b.O);
	LD tar1 = acos((x * x + z * z - y * y) / (2 * x * z));
	Vec i1 = a.O + rotate(Line(a.O, b.O), tar1) / z * x;
	LD tar2 = acos((y * y + z * z - x * x) / (2 * y * z));
	Vec i2 = b.O + rotate(Line(b.O, a.O), tar2) / z * y;
	return {i1, i2};
}

求多邊形周長

double PloygonDis(Vec *a, int n){
	double sum = 0;
	lfor(i, 1, n) sum += Dis(a[i], a[i % n + 1]);
	return sum;
}

求多邊形面積

double PloygonArea(Vec *a, int n){
	double sum = 0;
	lfor(i, 3, n) sum += cross(a[i - 1] - a[1], a[i] - a[1]);
	return sum / 2;
}

極角排序

實現

使用 atan2(y, x) 函式，返回值的範圍是 \([-\pi,\pi]\)；
使用叉積大於 \(0\) 的性質，注意因為叉積無法判 180 度以上，所以可能要結合象限排序；

凸包

二維凸包

特殊結論

二維平面四個點求凸包面積 \(\rightarrow\) 任選三個點面積之和 / 2

二維平面三個點面積 \(\rightarrow\) 二個二維向量行列式值的絕對值 / 2
三維空間五個點求凸包體積 \(\rightarrow\) 任選四個點體積之和 / 2

三維空間四個點體積 \(\rightarrow\) 三個三維向量行列式值的絕對值 / 6

半平面交

定義

有若干條有向直線，要求保留每條直線其中一側，求最後保留的範圍。

離線 \(O(n\log n)\)

用有向直線（一個點和一個方向向量）表示半平面，以下預設半平面在有向直線的左側。
對有向直線按方向向量的極角排序，維護一個雙端佇列，儲存當前構成半平面的直線以及相鄰兩直線的交點。
每次加入一條有向直線，如果隊首 / 隊尾的交點在直線右側（用叉積判）則彈掉隊首 / 隊尾的直線。
需要注意的細節：

加入直線時，先彈隊尾，再彈隊首。
特判平行直線，在右側的要彈掉。
最後還要檢查隊尾交點是否在隊首直線的右側，如果是也要彈掉。
如果題目給出的半平面不一定有限制邊界，則應該手動加入一個 INF 邊界。

double HPI(Line *a, int n){
	static deque <Vec> I; while(!I.empty()) I.pop_back(); 
	static deque <Line> Q; while(!Q.empty()) Q.pop_back();
	sort(a + 1, a + n + 1);
	Q.push_back(a[1]);
	lfor(i, 2, n){
		while(!I.empty() && Cross(a[i].v, I.back() - a[i].p) <= 0) I.pop_back(), Q.pop_back();
		while(!I.empty() && Cross(a[i].v, I.front() - a[i].p) <= 0) I.pop_front(), Q.pop_front();
		if(a[i].at2 != Q.back().at2) Q.push_back(a[i]);
		else if(Cross(a[i].v, Q.back().p - a[i].p) <= 0){
			Q.back() = a[i]; if(!I.empty()) I.pop_back();
		}else continue;
		auto qwq = Q.rbegin();
		if(Q.size() > 1) I.push_back(Inter(*(++qwq), Q.back()));
	}
	while(!I.empty() && Cross(Q.front().v, I.back() - Q.front().p) <= 0) I.pop_back(), Q.pop_back(); 
	if(Q.size() > 1) I.push_back(Inter(Q.front(), Q.back()));
	int cnt = 0; static Vec *b = new Vec[I.size() + 1];
	for(auto x : I) b[++cnt] = x;
	return PloygonDis(b, cnt);
}

積分

對積分的感性理解

\[\int_a^bf(x)\mathrm{d}x \]

有一個函式 \(f(x)\)，求其在區間 \([a,b]\) 與 \(x\) 軸圍成的面積，\(x\) 軸上為正，\(x\) 軸下為負。
那麼 \(\int\) 類比與 \(\sum\) 符號，同時用 \(\mathrm{d}x\) 表示將 \(x\) 分成很多很多很小的份。
同時也不難意識到，這個空間是封閉的才可以求面積。

如何積分

大部分的函式都是無法精確積分的，於是採用一些公式來逼近。

自適應辛普森積分

公式：

\[\int_a^b f(x) \mathrm{d}x = \frac{(b - a)(f(a) + f(b)+ 4f(\frac{a+b}{2}))}{6} \]

二次函式的積分可以精確計算，辛普森積分即是一種拿二次函式來擬合的方式。
在二次函式的情況下，該式求出的即準確積分，推導過程。

自適應：

因為 \(f(x)\) 不是二次函式，那麼當然不能直接積分，於是就有了根據誤差調整的自適應做法。

LD Ars(LD l, LD r, LD eps, LD val){
	LD mid = (l + r) / 2;
	LD L = simpson(l, mid), R = simpson(mid, r);
	if(fabs(L + R - val) <= eps) return L + R;
	return Ars(l, mid, eps / 2, L) + Ars(mid, r, eps / 2, R);
}

非常直接的想法，如果不滿足 eps，就繼續細分割槽間。
注意因為要合併兩個區間，所以對誤差的要求在提高。
因為是自適應的，所以複雜度玄學。
對於擬合易錯的函式，可強制迭代一定層數。

積分在 OI 中的應用

一般用來求各種面積。

[CQOI2005]三角形面積並

沒封裝的簡陋玩意

const LD Pi = acos(-1);

struct Vec{ LD x, y; };
void Out(Vec a){ cerr << a.x << ' ' << a.y << endl; }
void In(Vec &a){ scanf("%Lf%Lf", &a.x, &a.y); }
bool operator ==(Vec a, Vec b){ return a.x == b.x && a.y == b.y; }
bool operator !=(Vec a, Vec b){ return a.x != b.x || a.y != b.y; }
bool operator <(Vec a, Vec b){ return atan2(a.y, a.x) < atan2(b.y, b.x); }
LD atan2(Vec a){ return atan2(a.y, a.x); }
LD Cross(Vec a, Vec b){ return a.x * b.y - a.y * b.x; }
LD Dis(Vec a, Vec b){ return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y)); }
Vec operator -(Vec a, Vec b){ return {a.x - b.x, a.y - b.y}; }
Vec operator +(Vec a, Vec b){ return {a.x + b.x, a.y + b.y}; }
Vec operator /(Vec a, LD b){ return {a.x / b, a.y / b}; }
Vec operator *(Vec a, LD b){ return {a.x * b, a.y * b}; }
Vec operator *(LD b, Vec a){ return {a.x * b, a.y * b}; }
Vec rotate(Vec a, LD b){
	LD s = sin(b), c = cos(b);
	return {a.x * c - a.y * s, a.x * s + a.y * c};
}

struct Line{ 
	Vec p, v; LD at2; 
	Line(){}
	Line(Vec a, Vec b, LD c){ p = a, v = b, at2 = c; }
	Line(Vec a, Vec b){ p = a, v = b - a, at2 = atan2(v.y, v.x); } 
};
bool operator <(Line a, Line b){ return a.at2 < b.at2; }
bool Left(Vec a, Line b){ return Cross(a - b.p , b.v) > 0; }
Vec Inter(Line a, Line b){ return a.p + Cross(b.v, (b.p - a.p)) / Cross(b.v, a.v) * a.v;} 
Vec operator +(Vec a, Line b){ return {a.x + b.v.x, a.y + b.v.y}; }
Line rotate(Line a, LD b){ return {a.p, rotate(a.v, b), a.at2}; }
Line operator *(Line a, LD b){ return (Line){a.p, (Vec){a.v.x * b, a.v.y * b}, a.at2}; }
Line operator /(Line a, LD b){ return (Line){a.p, (Vec){a.v.x / b, a.v.y / b}, a.at2}; }

struct Circle{ Vec O; LD r; }; 
void In(Circle &a){ scanf("%Lf", &a.r), In(a.O); }
int Pos(Circle a, Circle b){
	LD dis = Dis(a.O, b.O); if(a.r < b.r) swap(a, b);
	if(dis > a.r + b.r) return 4;
	if(dis == a.r + b.r) return 3;
	if(dis > a.r - b.r) return 2;
	if(dis == a.r - b.r) return 1;
	return 0;
}
pair <Vec, Vec> Inter(Circle a, Circle b){
	LD x = a.r, y = b.r, z = Dis(a.O, b.O);
	LD tar1 = acos((x * x + z * z - y * y) / (2 * x * z));
	Vec i1 = a.O + rotate(Line(a.O, b.O), tar1) / z * x;
	LD tar2 = acos((y * y + z * z - x * x) / (2 * y * z));
	Vec i2 = b.O + rotate(Line(b.O, a.O), tar2) / z * y;
	return {i1, i2};
}

double PloygonArea(Vec *a, int n){
	double sum = 0;
	lfor(i, 3, n) sum += Cross(a[i - 1] - a[1], a[i] - a[1]);
	return sum / 2;
}
double PloygonDis(Vec *a, int n){
	double sum = 0;
	lfor(i, 1, n) sum += Dis(a[i], a[i % n + 1]);
	return sum;
}
double HPI(Line *a, int n){
	static deque <Vec> I; while(!I.empty()) I.pop_back(); 
	static deque <Line> Q; while(!Q.empty()) Q.pop_back();
	sort(a + 1, a + n + 1);
	Q.push_back(a[1]);
	lfor(i, 2, n){
		while(!I.empty() && Cross(a[i].v, I.back() - a[i].p) <= 0) I.pop_back(), Q.pop_back();
		while(!I.empty() && Cross(a[i].v, I.front() - a[i].p) <= 0) I.pop_front(), Q.pop_front();
		if(a[i].at2 != Q.back().at2) Q.push_back(a[i]);
		else if(Cross(a[i].v, Q.back().p - a[i].p) <= 0){
			Q.back() = a[i]; if(!I.empty()) I.pop_back();
		}else continue;
		auto qwq = Q.rbegin();
		if(Q.size() > 1) I.push_back(Inter(*(++qwq), Q.back()));
	}
	while(!I.empty() && Cross(Q.front().v, I.back() - Q.front().p) <= 0) I.pop_back(), Q.pop_back(); 
	if(Q.size() > 1) I.push_back(Inter(Q.front(), Q.back()));
	int cnt = 0; static Vec *b = new Vec[I.size() + 1];
	for(auto x : I) b[++cnt] = x;
	return PloygonDis(b, cnt);
}

【學習筆記】計算幾何
2022-06-07
筆記
計算幾何
2024-11-11
計算幾何：模板
2020-04-04
計算幾何模板
2024-08-25
【計算幾何】向量表示
2017-02-24
【總結】計算幾何模板
2017-05-15
二維計算幾何模板
2017-08-13
【計算幾何】線段相交
2013-08-01
三維計算幾何模板
2015-04-14
Something about 計算幾何
2024-07-08
計算幾何 —— 二維幾何基礎 —— 距離度量方法
2018-10-29
邊緣計算、霧計算、雲端計算區別幾何？
2020-11-08
計算幾何_向量的實現
2017-08-08
【計算幾何】多邊形交集
2013-08-02
計算幾何——平面最近點對
2024-05-23
POJ - 1556 【計算幾何 + 最短路】
2020-10-14
An Easy Problem?! POJ 2826 計算幾何
2017-11-13
計算幾何常用的函式/方法
2015-09-02
函式
BNUOJ 12887 isumi(計算幾何+最大流)
2019-02-21
SGU 124 Broken line(計算幾何)
2019-02-27
【計算幾何】Triangles HUST 1607
2017-05-08
【計算幾何】多邊形點集排序
2013-08-01
排序
C++計算幾何演算法大全
2016-11-23
C++演算法
POJ 2991 Crane（線段樹+計算幾何）
2016-01-12
【計算幾何】點在多邊形內部
2013-08-01
POJ 1556 The Doors(Dijkstra+計算幾何)
2011-09-17
二維幾何常用運算
2013-11-18
微分幾何學習筆記 (2)- 微分流形(Manifold)
2017-11-30
筆記
產品幾何技術規範學習筆記
2024-09-29
筆記
SGU 120 SGU 228 Archipelago(計算幾何)
2019-02-25
Go
計算機視覺—圖片幾何變換（2）
2019-02-27
計算機視覺
【計算幾何】求線段相交交點座標
2017-03-04
BZOJ 1027 合金計算幾何，Floyd判環
2017-01-19
CG_Hadoop：基於MapReduce的計算幾何
2014-12-26
Hadoop
計算機圖形學（四）_幾何變換_1_基本的二維幾何變換(一)
2017-03-02
計算機
POJ 1113 Wall(思維計算幾何數學)
2019-02-14
計算幾何（一）：凸包問題（Convex Hull）
2020-09-22
《計算廣告》學習筆記
2017-10-31
筆記