計算幾何（凸包與旋轉卡殼）
一、二維向量
1、定義，
- 二維向量是二維平面中的有向線段，
- \(\vec{a}=(x,y)\) 表示從 \((0,0)\) 運動到 \((x,y)\)，
- 向量的模：向量長度，記為 \(|\vec{a}|\)，
2、向量運算，
- 單位向量：三維空間中，\(\vec{i}=(1,0,0)，\vec{j}=(0,1,0)，\vec{k}=(0,0,1)\)，則三維空間中任一向量 \(\vec{a}\) 均唯一可被表示為 \(\vec{i}，\vec{j}，\vec{k}\) 的線性組合，即 \(\vec{a}=x\vec{i} +y\vec{j} +z\vec{k}\)，此時 \(\vec{i}，\vec{j}，\vec{k}\) 張成 \(R^3\)，為 \(R^3\) 的單位向量（標準基），
（1）向量加法，
- 對于向量 \(\vec{a}=(x_a,y_a)，\vec{b}=(x_b,y_b)\)，有 \(\vec{a}+\vec{b}=(x_a+x_b,y_a+y_b)\)，
- 對于減法而言，將上述式子中加號換成減號即可，
（2) 向量積，
- 點積（標量積，內積）
- 定義：兩向量水平分量的乘積，
- 幾何定義式：\(\vec{a} \cdot \vec{b}=|\vec{a}||\vec{b}| \cos \theta\)，
- 代數定義式：\(\vec{a} \cdot \vec{b}=x_ax_b +y_ay_b\)，
- 線性代數定義式：\(\vec{a}^T \vec{b}=\begin{bmatrix}x_1 & x_2\end{bmatrix} \begin{bmatrix}y_1\\y_2\end{bmatrix}=x_1y_1 +x_2y_2\)，
- 更一般地，有 \(n\) 維向量點積：\(\vec{a} \cdot \vec{b}=\sum_{i=1}^{n}a_ib_i\)，
- 點積有分配律，但因為它的結果是標量，所以點積沒有方向，
- 叉積（矢量積，外積）
- 定義：兩向量垂直分量的乘積，叉積的結果是向量，所以有方向，
- 幾何定義式：\(\vec{a} \times \vec{b}=|\vec{a}||\vec{b}| \sin \theta\)，注意：該式子只能求出叉積的模，不能得到方向，
- 代數定義式：三維空間中，設 \(\vec{a}=(x_a,y_a,z_a)，\vec{b}=(x_b,y_b,z_b)\)，\(則\vec{a} \times \vec{b}=det\begin{vmatrix}\vec{i} & \vec{j} & \vec{k} \\ x_a & y_a & z_a \\ x_b & y_b & z_b \end{vmatrix}=(y_az_b-z_ay_b)\vec{i} + (x_az_b-z_ax_b)\vec{j} +(x_ay_b-x_by_a)\vec{k}\)
- 代數定義式二維情況下，\(\vec{a} \times \vec{b}=(x_ay_b-x_by_a)\vec{k}\)，
- 線性代數定義式（外積展開）:\(\vec{a} \times \vec{b}=\begin{bmatrix} x_a \\ y_a \\ z_a \end{bmatrix} \begin{bmatrix}x_b & y_b & z_b \end{bmatrix}=\begin{bmatrix} x_ax_b & x_ay_b & x_az_b \\ y_ax_b & y_ay_b & y_az_b \\ z_ax_b & z_ay_b & z_az_b \end{bmatrix}\)，
- 右手法則：\(\vec{a} \times \vec{b}\) 時，右手由 \(\vec{a}\) 握向 \(\vec{b}\)，大拇指方向為叉積方向，
- 叉積可快速計算由相鄰向量形成的平行四邊形面積，
- 應用：
- 判斷向量位置關系（見凸包），
二、凸包
1、二維凸包，
（1）定義，
- 凸多邊形：所有內角大小都在 \([0,\pi]\) 范圍內的簡單多邊形，
- 凸包：平面上能包含所有給定點的最小凸多邊形，
（2）構建方法，
- 找到一個必定在凸包上的點，然后按順時針或逆時針順序依次加入所有點
- 1、極角序，Graham掃描法，
- 演算法流程：
- 找到縱坐標最小的點，該點必在凸包上，
- 將該點設為原點，剩下的點用叉積進行極角排序，
- 用單調堆疊維護，若兩向量叉積 \(\geq 0\)，則另一個向量不在凸包上，
- 主要函式實作：
- 向量結構體：
struct Vector{ double x, y; Vector(const double &x, const double &y){ this->x = x; this->y = y; } Vector(){} void read(){ scanf("%lf%lf", &x, &y); } double length(){ return sqrt(x * x + y * y); } } Vector a[N], s[N];//N由題目確定,a陣列存盤給定點，s陣列存盤凸包上的點， Vector operator - (const Vector &a, const Vector &b){ return Vector(a.x - b.x, a.y - b.y); }- 叉積函式：
double cross(Vector a, Vector b){ return a.x * b.y - a.y * b.x; }- 極角序計算函式：（這里使用相對距離計算極角序，精確度更高）
double dis(Vector a, Vector b){ Vector c = a - b; return c.length(); }- 極角序比較函式：
int cmp(Vector p, Vector q){ double x = cross(p - a[1], q - a[1]); if(x > 0) return 1; if(x == 0 && dis(p, a[1]) < dis(q, a[1])) return 1; return 0; }- 叉積判斷向量位置關系：（左正右負）
double multi(Vector a, Vector b, Vector c){ return cross(b - a, c - a); }- 主函式核心部分：
int main(){//n為給定點的個數 ....... int k = 1; for(int i = 2; i <= n; i++) if(a[i].y < a[k].y || (a[k].y == a[i].y && a[i].x < a[k].x)) k = i; swap(a[1], a[k]);//得出原點， sort(a + 2, a + 1 + n, cmp); s[1] = a[1]; s[2] = a[2]; int t = 2; for(int i = 3; i <= n; i++){ while(t >= 2 && multi(s[t - 1], s[t], a[i]) <= 0) t--;//題目要求凸包上的點可共線時取等號，否則不取， s[++t] = a[i]; } ....... }
- 演算法流程：
- 2、水平序，Andrew演算法
- 演算法流程：
- 以 \(x\) 為第一關鍵字，\(y\) 為第二關鍵字排序，排序后，最小與最大元素必在凸包上，
- 最大與最小元素連線處切開分上下凸殼掃描，
- 演算法流程：
- 判斷點是否在凸包內，
- 以凸包最下面的點為原點，
- 變換坐標，
- 找到與該點極角最接近的兩點形成的線段，用叉積來判定，
（3）例題
- 例題1：城墻
一個貪婪的國王要求他的建筑師建一堵墻圍繞他的城堡，且墻與城堡之間的距離總不小于一個數 \(L\) ，輸入城堡各個節點（圖中實線與實線的交點）的坐標和 \(L\) ，要求最小的墻壁周長，
輸入格式：
輸入第一行 \(N(3\leq N \leq 1000)\)和 \(L(1\leq L \leq 1000)\)，其中 \(N\) 為節點個數，以下 \(N\) 行每行是各個節點的橫坐標 \(x_i\) 和縱坐標 \(y_i\) ，其中 \(-10^4\leq x_i,y_i \leq 10^4\)，
輸出格式：
輸出僅一個數，為最小的墻壁周長(四舍五入至整數)，
輸入樣例：
9 100
200 400
300 400
300 300
400 300
400 400
500 400
500 200
350 200
200 200
輸出樣例:
1628
這道題的核心是墻與城堡之間的距離總不小于一個數 \(L\) ，但是幾何這玩意有難想象，所以這時候就要發揮 \(OIer\) 的傳統手藝（畫圖）了！
隨便畫了個圖出來，圓弧部分距離每個點的距離均為 \(L\)，仔細觀察一下紅線圈出來的區域，有什么發現？
這些圓弧拼起來恰巧是一個圓，不信的同志可以自己畫來算，很簡單的，
再仔（隨）細（便）觀（一）察（猜），發現每一條直線就是兩個點間的距離，所以問題就轉化成了：求凸包周長，求出來后再加上一個圓的周長就OK了，妥妥的模板題了，
//AC代碼，壓行嚴重，勿噴，
#include<bits/stdc++.h>
using namespace std;
const double PI = acos(-1.0);
int n, L;
struct point{
double x, y;
point(const double &x, const double &y){ this->x = x; this->y = y; }
point(){}
void read(){ scanf("%lf%lf", &x, &y); }
double length(){ return sqrt(x * x + y * y); }
}a[1005], s[1005];
point operator - (const point &a, const point &b){return point(a.x - b.x, a.y - b.y);}
double X(point a, point b){ return a.x * b.y - a.y * b.x; }
double dis(point a, point b){ point c = a - b; return c.length(); }
int cmp(point p, point q){
double x = X( p - a[1], q - a[1] );
if(x > 0) return 1;
if(x == 0 && dis(p, a[1]) < dis(q, a[1])) return 1;
return 0;
}
double multi(point a, point b, point c){ return X(b - a, c - a); }
int Round(double x){
int r;
if(x - (int)x > 0.5)r = (int)x + 1;
else r = (int)x;
return r;
}
int main(){
scanf("%d%d", &n, &L);
for(int i = 1; i <= n; i++) a[i].read();
int k = 1;
for(int i = 2; i <= n; i++) if(a[i].y < a[k].y || (a[k].y == a[i].y && a[i].x < a[k].x)) k = i;
swap(a[1], a[k]);
sort(a + 2, a + 1 + n, cmp);
s[1] = a[1];
s[2] = a[2];
int t = 2;
for(int i = 3; i <= n; i++){
while(t >= 2 && multi(s[t - 1], s[t], a[i]) <= 0) t--;
s[++t] = a[i];
}
// cout<<endl;
// for(int i = 1; i <= t; i++)cout<<s[i].x<<" "<<s[i].y<<endl;
double sum = 0;
for(int i = 1; i < t; i++) {
sum += dis(s[i], s[i + 1]);
}
sum += dis(s[1], s[t]) + PI * L * 2;
cout<<Round(sum);
return 0;
}
- 例題2：信用卡凸包（SHOI2012）
信用卡是一個矩形，唯四個角作了圓滑處理，使它們都是與矩形的兩邊相切的 1/4 圓，如下圖所示，現在平面上有一些規格相同的信用卡，試求其凸包的周長，注意凸包未必是多邊形，因為它可能包含若干段圓弧，
輸入格式：
輸入的第一行是一個正整數 \(n\)，表示信用卡的張數，第二行包含三個實數 \(a,b,r\)，分別表示信用卡（圓滑處理前）豎直方向的長度、水平方向的長度，以及 \(\frac{1}{4}\) 圓的半徑，
之后 \(n\) 行，每行包含三個實數 \(x,y,\theta\)，分別表示一張信用卡中心（即對角線交點）的橫、縱坐標，以及繞中心 逆時針旋轉的 弧度，
輸出格式：
輸出只有一行，包含一個實數，表示凸包的周長，四舍五入精確到小數點后2位，
輸入樣例1：
2
6.0 2.0 0.0
0.0 0.0 0.0
2.0 -2.0 1.5707963268
輸出樣例1：
21.66輸入樣例2:
3
6.0 6.0 1.0
4.0 4.0 0.0
0.0 8.0 0.0
0.0 0.0 0.0
輸出樣例2：
41.60輸入樣例3：
3
6.0 6.0 1.0
4.0 4.0 0.1745329252
0.0 8.0 0.3490658504
0.0 0.0 0.5235987756
輸出樣例3：
41.63樣例3說明：
分析：
- 考慮 \(r = 0\) 的情況，發現此時信用卡就是矩形，所以跑一遍Graham演算法就可以了，
- 考慮 \(r \geq 1\) 的情況，經過多次傳統手藝（畫圖）分析，發現在多張信用卡形成的帶圓弧的凸包上，圓弧所對應的圓心角之和總為 \(2\pi\)，適用我們只需要跑一遍Graham再加上一個整圓的周長就可以了，這個規律與例題1相同，
- 設當前點為 \(A(x,y)\)，則它繞原點逆時針旋轉 \(\theta\) 后的坐標為 \(A′(xcos?\theta ?ysin?\theta ,xsin?\theta +ycos\theta )\)，
有了這些，這個題目就成了一個裸的凸包模板，
//AC代碼
#include<bits/stdc++.h>
#define MAXN 10005
using namespace std;
const double PI = 3.141592653589793;
const double eps = 1e-23;
int n, t1;
double A, B, R;
struct point{
double x, y;
point(const double &x, const double &y){ this->x = x; this->y = y; }
point(){}
double length(){ return sqrt(x * x + y * y); }
}a[MAXN<<2], s[MAXN<<2];
bool operator == (const point &a, const point &b){
if(a.x == b.x && a.y == b.y) return true;
return false;
}
point operator + (const point &a, const point &b){ return point(a.x + b.x, a.y + b.y); }
point operator - (const point &a, const point &b){ return point(a.x - b.x, a.y - b.y); }
double X(point a, point b){ return a.x * b.y - a.y * b.x; }
double dis(point a, point b){ point c = a - b; return c.length(); }
double multi(point a, point b, point c){ return X(b - a, c - a); }
int cmp(point p, point q){
double x = X(p - a[1], q - a[1]);
if(x > 0) return 1;
if(x == 0 && dis(p, a[1]) < dis(q, a[1])) return 1;
return 0;
}
point rotate(double x, double y, double t){
double c = cos(t), s = sin(t);
return point( x * c - y * s, x * s + y * c );
}
double Round(double x){
double r;
if(x - floor(x) - 0.99 > 0.005)r = x - floor(x) + 0.01 + (int)x;
else r = x;
return x;
}
int main(){
scanf("%d", &n);
scanf("%lf%lf%lf", &A, &B, &R);
A/=2, B/=2;
double cir = PI * 2.0 * R;
for(int i = 1; i <= n; i++){
double x, y, rou;
scanf("%lf%lf%lf", &x, &y, &rou);
x+=eps; y += eps; rou += eps;
point pp = point(x, y);
a[++t1] = rotate(B - R, A - R, rou) + pp;
a[++t1] = rotate(-B + R, A - R, rou) + pp;
a[++t1] = rotate(B - R, -A + R, rou) + pp;
a[++t1] = rotate(-B + R, -A + R, rou) + pp;
}
int k = 1;
for(int i = 1; i <= t1; i++) if(a[i].y < a[k].y || (a[i].x < a[k].x && a[i].y == a[k].y)) k = i;
swap(a[k], a[1]);
sort(a + 2, a + 1 + t1, cmp);
s[1] = a[1];
int cnt = 1;
for(int i = 2; i <= t1; i++){
while(cnt > 1 && multi(s[cnt - 1], s[cnt], a[i]) <= 0) cnt--;
s[++cnt] = a[i];
}
s[++cnt] = a[1];
double sum = 0;
for(int i = 1; i < cnt; i++)sum += dis(s[i], s[i + 1]);
printf("%.2lf", Round(sum + cir));
return 0;
}
2、動態凸包，
- 使用平衡樹進行動態插入與排序
- 主要思路：快速定位凸包需要調整的部分，然后進行調整，
- 為統一排序標準，將凸包分為上下凸殼，
- 上凸殼關于 \(x\) 軸對稱后就變為下凸殼，方便操作，
- 水平序，極角序均可，
- 插入一個點時，判斷是否在凸包下方，
- 如果是，則插入且不斷往兩邊刪點，直到所有點都滿足凸包定義為止，
- 注意邊界處的點不要同時加入上下凸包，
- 為統一排序標準，將凸包分為上下凸殼，
- 每個點只會加入和彈出平衡樹一次，所以時間復雜度為 \(O(N\log^2N)\)
- 動態凸包不支持洗掉，
- 離線后倒序處理，將洗掉轉為加點，
- 主要思路：快速定位凸包需要調整的部分，然后進行調整，
- 線段樹分治：建立凸包時，按照歸并思想合并左右子樹凸包，避免重新排序，
- \(CDQ\) 分治
- 二進制分組 \(\&\) 定期重構
- 缺點：不支持洗掉操作，
三、旋轉卡殼
1、前置知識，
（1）凸多邊形的直徑，
- 定義：一個凸多邊形上任意兩點距離的最大值為凸多邊形的直徑，
- 確定直徑的點可能多于一對
- 一個凸 \(N\) 邊形可能存在 \(n\) 對直徑，
（2）支撐線，
- 直線 \(L\) 過凸多邊形 \(P\) 一個頂點，使 \(P\) 全在 \(L\) 一側，此時稱直線 \(L\) 為凸多邊形 \(P\) 的支撐線，
（3）對踵點，
- 兩條平行直線過凸多邊形的兩個頂點且將凸多邊形夾在其之間，此時這對頂點稱作對踵點，
- 兩條平行支撐線所過兩點是一對對踵點，
- 一個凸 \(N\) 邊形的對踵點最多有 \(\left \lceil \frac{3N}{2} \right \rceil\) 對，
(4）凸多邊形直徑定理，
- 凸多邊形 \(P\) 的直徑 \(=\) 凸多邊形 \(P\) 的所有平行支撐線之間距離的最大值，
- 凸多邊形的直徑可在對踵點對中尋找，
2、旋轉卡殼演算法，
（1）旋轉卡殼：尋找凸多邊形對踵點的演算法，
- 主要思路：利用凸包單調性，列舉點 -> 更新點對 -> 更新答案，
- 演算法流程：
- 求出凸多邊形凸包，此時最遠點對必在凸包上，時間復雜度為 \(O(N\log N)\)，
- 找到一對初始對踵點，作水平切線，
- 選 \(y_{max}\) 與 \(y_{min}\)，
- 更新答案，
- 沿一個方向旋轉這兩條切線，知道其中一條與凸包的邊重疊，
- 產生新對踵點對，更新答案，
- 不斷重復這一程序，直到所有點對都已生成對踵點，
- 該演算法本質上是對于每一條凸包的邊，計算最遠的過凸包上點的平行直線，
- 兩直線距離可看作三角形面積除以底邊，用叉積計算，
- 代碼中心：
int rotating_calipers(){ int re = 0; if(t == 2) return dis(s[1], s[2]); int j = 2; for(int i = 1; i <=t; i++){ while(multi(s[i], s[i % t + 1], s[j]) <= multi(s[i], s[i % t + 1], s[j % t + 1])) j = j % t+1; re = max(re, max(dis(s[i], s[j]), dis(s[i % t+ 1], s[j]))); } return re;//回傳的是凸多邊形直徑， }
（2）例題
給定平面上 \(n\) 個點，求凸包直徑，
輸入格式：
第一行一個正整數 \(n\)，
接下來 \(n\) 行，每行兩個整數 \(x,y\)，表示一個點的坐標，
輸出格式：
輸出一行一個整數，表示答案的平方，
輸入樣例：
4
0 0
0 1
1 1
1 0
輸出樣例：
2
說明/提示：
【資料范圍】
對于 \(100\%\) 的資料，\(2\leq n \leq 50000，|x|,|y|\leq 10^4\)，
簡單的模板題，直接上程式，
//AC代碼
#include<bits/stdc++.h>
#define MAXN 50005
using namespace std;
int n, t;
struct point{
int x, y;
point(int x, int y):x(x), y(y){}
point(){}
void read(){ scanf("%d%d", &x, &y); }
int length(){ return x * x + y * y; }
}a[MAXN], s[MAXN];
point operator - (const point &a, const point &b){
return point(a.x - b.x, a.y - b.y);
}
int X(point a, point b){
return a.x * b.y - a.y * b.x;
}
int multi(point a, point b, point c){ return X(b - a, c - a); }
int dis(point a, point b){
point c = a - b;
return c.length();
}
int cmp(point p, point q){
int x = X(p - a[1], q - a[1]);
if(x > 0) return 1;
if(x == 0 && dis(p, a[1]) < dis(q, a[1])) return 1;
return 0;
}
void convex_hull(){
int k = 1;
for(int i = 1; i <= n; i++)
if(a[i].y < a[k].y || (a[i].y == a[k].y && a[i].x < a[k].x))
k = i;
swap(a[1], a[k]);
sort(a + 2, a + 1 + n, cmp);
s[1] = a[1];
t = 1;
for(int i = 2; i <= n; i++){
while(t > 1 && multi(s[t - 1], s[t], a[i]) <= 0)
t--;
s[++t] = a[i];
}
}
int rotating_calipers(){
int re = 0;
if(t == 2) return dis(s[1], s[2]);
int j = 2;
for(int i = 1; i <=t; i++){
while(multi(s[i], s[i % t + 1], s[j]) <= multi(s[i], s[i % t + 1], s[j % t + 1]))
j = j % t+1;
re = max(re, max(dis(s[i], s[j]), dis(s[i % t+ 1], s[j])));
}
return re;
}
int main(){
scanf("%d", &n);
for(int i = 1; i <= n; i++) a[i].read();
convex_hull();
printf("%d", rotating_calipers());
return 0;
}
（3）應用，
- 凸多邊形直徑與寬，
- 凸多邊形間距離，
- 凸多邊形最小面積/周長外接矩形，
- 定理：關于凸多邊形 \(P\) 的一個外接矩形與 \(P\) 的一邊重合，
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/297725.html
標籤：其他