[筆記] 計算幾何

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;
}

極角排序

實現

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

凸包

二維凸包

特殊結論

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

    二維平面三個點面積 \(\rightarrow\) 二個二維向量行列式值的絕對值 / 2

  • 三維空間五個點求凸包體積 \(\rightarrow\) 任選四個點體積之和 / 2

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

半平面交

定義

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

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

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

  1. 加入直線時,先彈隊尾,再彈隊首。
  2. 特判平行直線,在右側的要彈掉。
  3. 最後還要檢查隊尾交點是否在隊首直線的右側,如果是也要彈掉。
  4. 如果題目給出的半平面不一定有限制邊界,則應該手動加入一個 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);
}

相關文章