計算幾何
AaW
目錄
-
座標與向量
-
有向面積
-
線段相交
-
極角排序
-
凸包演算法
今天講的幾何都是在二維平面的拉
「計算幾何是以數學方式精確與高效的計算出幾何結構與其性質的嚴謹應用。」
座標與向量
直角坐標系
\((x,\, y)\)
極坐標系
\((r,\, \theta)\)
向量
- 有向線段(向:方向、量:長度)
- 純量:只有大小、沒有方向的量
- \(\overrightarrow{AB} \) :由 A 點起始指向 B 點的向量
- 可任意移動
- 向量長度(加絕對值表示):兩點的歐基里德距離
\(\left|\overrightarrow{AB}\right| = \sqrt{{(x_A - x_B)}^2 + {(y_A - y_B)}^2}\)
向量怎麼實作?
struct pt {
T X, Y; // T 可以為 int 或是 double
}
#define pt pair<int, int>
#define X first
#define Y second
pt v1;
cin >> v1.X >> v1.Y;
或是
座標也可以用一樣的方式實作,例如:
\(A\)點的座標可以視為\(\overrightarrow{OA}\) 來處理
向量的運算
向量加法
- 頭尾相連:\( \overrightarrow{AB} + \overrightarrow{BC} = \overrightarrow{AC}\)
- 在圖上可以用三角形法或平行四邊形法
- 代數運算:
- \( \overrightarrow{a} = (a_1, a_2) , \overrightarrow{b} = (b_1, b_2), \overrightarrow{c} = (c_1, c_2)\)
-
\( \overrightarrow{a} + \overrightarrow{b} = \overrightarrow{c}\)
- \(\Rightarrow \, a_1+b_1= c_1, \,c_2 = a_2 + b_2\)
向量減法
終點減起點
加上反向向量
Ex. \(\overrightarrow{AB} - \overrightarrow{CD} = \overrightarrow{AB} + \overrightarrow{DC}\)
向量內積 (DOT)
定義:若 \(\vec a\) 與 \(\vec b\) 為二度或三度空間的向量,且 \(\theta\) 為\(\vec a\) 與 \(\vec b\) 的夾角(滿足\(0\leq\theta\leq\pi\)),
則點積 \(\vec a\cdot\vec b\) 定義為
內積(DOT)
Dot Product
- 直觀看法:兩個向量在「同一方向」的程度
- 兩個向量的內積(點積)以「\(\cdot\)」表示
- 運算結果為一個純量
- 幾何關係:\(\vec A\cdot \vec B = |\vec A||\vec B|\cos\theta\)
- 正定性:\(|\vec A|^2 = \vec A\cdot\vec A\)
- 座標公式:假設\(A(x_A,y_A),B(x_B,y_B)\),則內積滿足:\(\vec A\cdot\vec B = x_Ax_B+y_Ay_B\)
- 怎麼證明這兩個公式相等?
- 將cos用餘弦定理代換!
Source: Wiki
內積(DOT)
Dot Product
- \(\vec A\cdot\vec B > 0\):兩向量之間的夾角為銳角
- \(\vec A\cdot\vec B = 0\):兩向量垂直
- \(\vec A\cdot\vec B < 0\):兩向量夾角為鈍角
等等會一直用到
外積(CROSS)
Cross product
定義:若\(a = (a_1,a_2,a_3)\)與\(b = (b_1,b_2,b_3)\)為三度空間裡的向量,則外積(叉積)\(a\times b\) 為一向量定義為:
向量的外積只有在三度空間才有定義!
外積結果的向量由右手定則決定
外積(CROSS)
Cross product
- 直觀看法:兩向量在逆時針方向的「垂直」的程度
- 兩個向量的外積(叉積)以「\(\times\)」表示
- 幾何關係:\(\vec A\times \vec B = \hat n|\vec A||\vec B|\sin\theta\)(\(\theta\)為有向角)
- 運算結果為一個向量(三維空間),方向由右手定則決定
- 外積結果同時也是兩向量所張的平行四邊形面積,我們在二維平面上將它拿來當作純量來使用
- 座標公式:假設\(A(x_A,y_A),B(x_B,y_B)\),則外積滿足:\(\vec A\times\vec B = x_Ay_B-x_By_A\)
Source: Wiki
Source: Wiki
外積(CROSS)
Cross product
- \(\vec A\times\vec B > 0\):A向量到B向量逆時針旋轉
- \(\vec A\times\vec B = 0\):兩向量平行
- \(\vec A\times\vec B < 0\):A向量到B向量順時針旋轉
比較
- 向量內積(Dot Product)
- 判斷向量是否垂直
- 結果為純量
- 向量外積(Cross Product)
- 判斷向量是否平行
- 結果為向量(今天都當純量來用)
- 判斷向量夾角方向 (順時針或是逆時針)
- 可用來計算面積
向量運算要怎麼實作?
用運算子多載
struct pt{
int x,y;
};
pt operator+(const pt &a, const pt &b) {
return {a.x + b.x, a.y + b.y};
} //向量相加
pt operator-(const pt &a, const pt &b) {
return {a.x - b.x, a.y - b.y};
} //向量相減
pt operator*(int k, const pt &a) {
return {a.x * k, a.y * k};
} //向量乘法(係數積)
pt operator*(const pt &a, int k) {
return {a.x * k, a.y * k};
} // 乘法有交換率
pt operator/(const pt &a, int k) {
return {a.x / k, a.y / k};
} //向量除法 (係數除法)
int operator^(const pt &a, const pt &b) {
return a.x * b.y - a.y * b.x;
} //向量外積cross
int operator*(const pt &a, const pt &b) {
return a.x * b.x + a.y * b.y;
} //向量內積dot
inline bool operator == (const pt &a, const pt &b){
if(abs(a.x-b.x) <= eps && abs(a.y - b.y) <= eps) return true;
return false;
}
Eps是什麼?
- 在進行有小數的計算幾何時,常常會有很多的浮點誤差
- eps為誤差容忍值
- 兩數相差eps內視爲相等
\(x ⇒ (x − eps, x + eps)\)
廣義三角函數
Source: https://pse.is/3suhzc
弧度
- \(2\pi = 360^{\circ},\pi = 180^{\circ}\)
- \(x = \frac{x^{\circ}}{180}\times\pi\)
- \(x^{\circ} = \frac{x}{\pi}\times 180^{\circ}\)
- 弧度沒有單位,以「rad」表示!
- C++ 裡面的三角函數都是弧度!
Source: Wiki
三角函數運算
- 平方關係:\(\sin^2\theta +\cos^2\theta = 1\)
- 正弦定理:\(\frac{c}{\sin(\gamma)} = \frac{a}{\sin(\alpha)} = \frac{b}{\sin(\beta)}\)
- 餘弦定理:\(cos(\gamma) = \frac{a^2+b^2-c^2}{2ab}\)
有向面積
性質
- 利用外積計算多邊形面積
- 有向角:逆時針為正、順時針為負
- 外積有正有負,稱為「有向面積」
+
-
可以幹嘛?
- 計算面積
- 評估角度
三角形面積
兩向量夾出的三角形有向面積為:
A
B
C
WHY?
三角形面積
幾何證明:簡單明瞭
多邊形面積
- 在平面上任選一點(通常是原點)
- 以此點為一端點,將多邊形切割成許多三角形。
- 依照逆時鐘的方向將所有三角形的有向面積相加即是多邊形的有向面積。
Source: https://reurl.cc/NZ5ZNq
多邊形面積
假設多邊形的頂點依序為(要特別注意結束的點必須是起始點):
則多邊形的面積為(aka測量師公式):
多邊形面積
Pick's Theorem
給定頂點座標均是格子點的簡單多邊形,皮克定理說明了其面積\(A\)和內部格點數目\(i\)、邊上格點數目\(b\)的關係:
證明在這裡
Pick's Theorem
Source: Wiki
題目
-
ZJ a871: Museum Area
- n個點圍成的多邊形,求面積。
-
CSES: Polygon Area
- 一樣是裸題
-
Polygon Lattice Points
- 裸皮克定理
線段相交
segment banana
直接求出交點,再判斷交點是否在線段內
如何判斷?
這樣好嗎?
聽起來很麻煩誒
求交點的問題
- 沒有交點?
- 無限多交點?
- 交點不好求 e.g. 垂直線
- 交點位置誤差?
方向函式
- 利用外積找兩向量相對的方向,方便判斷
- 共線則為回傳0,逆時針回傳1,順時針回傳-1
//方向函數ori,回傳OA向量與OB向量的方向
int ori( const Pt& a , const Pt& b, const Pt& o){ //方向函數
double cross = ( a - o ) ^ ( b - o );
if( fabs( cross ) < eps ) return 0;
return cross > 0 ? 1 : -1;
}
方向函式
方向函式能做什麼?
判斷點是否在線段兩端點的異側!
更精準的求線段相交
若線段\(AB\) 與\(CD\) 相交
則點\(C\)和\(D\)會在線段\(AB\)異側
=0代表C或D有一點在直線 AB 上面
點與線段的關係
- 判斷線段
香蕉之前,先處理點在線段上的情況 - onseg() 判斷點 O 是否在線段 AB 之上
- 當點在線段上時,OA向量與OB向量之外積為零、內積為負(夾角為180度)
bool onseg(pt a, pt b, pt o){ //o是否在ab線段上
int cross = (a - o) ^ (b - o); //是否平行
int dot = (a - o) * (b - o); //是否在線段中
return (cross == 0)&&(dot <= 0);
}
線段相交完整程式碼
- 首先判斷特例,也就是線段有一端點在另外線段上
- 利用方向函式判斷點是否在另一線段端點的異側
bool banana(pt a, pt b, pt c, pt d){ //線段ab是否與cd相交
if(onseg(a,b,c)||onseg(a,b,d)) return true; //點c、d是否洽在線段ab上
if(onseg(c,d,a)||onseg(c,d,b)) return true; //點a、b是否洽在線段cd上
return (ori(a,b,c)*ori(a,b,d) < 0 && ori(c,d,a)*ori(c,d,b) < 0);
}
資芽版的香蕉
線段交點座標
知道兩線段是否相交之後,進一步求出交點座標
其中,\(t\)為ACD和
線段交點座標
知道兩線段是否相交之後,進一步求出交點座標
線段交點座標
幾個算式列一列得到:
帶回 \(\vec{P} = \vec{A} + t\overrightarrow{AB}\) 即可
pt banana_point(pt a, pt b, pt c, pt d){
assert(banana(a, b, c, d)); //是否相交
return a + ((a - c)^(d - c) * (b - a))/((d - c)^(b - a));
}
Point in Polygon
給定平面上\(n\)個點構成的簡單多邊形,還有\(m\)個詢問的點。對每一個點輸出它是在多邊形內部、外部或邊界上。
$$3\leq n\leq 1000$$
$$1\leq m\leq 1000$$
$$10^9\leq x_i,y_i\leq 10^9$$
來看看凸多邊形的情況
- 由P點指向多邊形上點的向量依序繞一圈做外積,過程中結果皆為正,表示在多邊形內部。
- 若遇到凹多邊形,則這個方法會失敗。
射線作法
- 對每一個點做一條往正x方向的射線(實作上就是和一個x座標很大的點,連成線段)
- 計算每一條邊與這條線段相交的次數。
- 如果是奇數則在多邊形內部,反之則是在外部。
- 小心處理相交於多邊形頂點的情況
- 頂點會被算到兩次
- 相交頂點時只考慮 \(y\) 較小的點
- 把射線改成往右上角斜?
- 右端點設成\((\infty,y_i+1)\)
Source: geeksforgeeks
程式碼
#include <bits/stdc++.h>
#define int long long
#define double long double
#define pii pair<int, int>
#define N 1005
#define INF 1000000005
#define x first
#define y second
#define IOS ios::sync_with_stdio(0),cin.tie(0)
using namespace std;
int n,m;
struct point{
int x,y;
bool operator == (point b){
if(x == b.x && y == b.y)return true;
return false;
}
point operator - (point b){return {x - b.x , y - b.y};}
point operator + (point b){return {x + b.x , y + b.y};}
int operator ^ (point b){return (x * b.y - y * b.x);}
int operator * (point b){return (x * b.x + y * b.y);}
}pt[N];
bool onseg(point a,point b,point o){
int cross = (a - o) ^ (b - o);
int dot = (a - o) * (b - o);
return (cross == 0) && (dot <= 0);
}
int dir(point a,point b,point o){
int cross = (a - o) ^ (b - o);
if(cross == 0)return 0;
else if(cross > 0)return 1;
else return -1;
}
bool inter(point a,point b,point c,point d){
if(onseg(a,b,c) || onseg(a,b,d))return true;
if(onseg(c,d,a) || onseg(c,d,b))return true;
if(dir(a,b,c) * dir(a,b,d) < 0 && dir(c,d,a) * dir(c,d,b) < 0)
return true;
return false;
}
int solve(point cur){
int sum = 0;
for(int i=0;i<n;i++){
if(onseg(pt[i],pt[i+1],cur) == 1)return -1;
if(inter(pt[i],pt[i+1],cur,point{INF,cur.y}))sum++;
point temp = pt[i].y > pt[i+1].y ? pt[i] : pt[i+1];
if(temp.y == cur.y && temp.x > cur.x)sum--;
}
return sum;
}
signed main(){
IOS;
cin>>n>>m;
for(int i=0;i<n;i++)cin>>pt[i].x>>pt[i].y;
pt[n] = pt[0];
for(int i=0;i<m;i++){
point temp;cin>>temp.x>>temp.y;
int ans = solve(temp);
if(ans == -1)cout<<"BOUNDARY"<<"\n";
else if(ans & 1)cout<<"INSIDE"<<"\n";
else cout<<"OUTSIDE"<<"\n";
}
}
點到線段距離
問 C 點到 \(\overline{\rm AB}\) 的距離平方
- 當ABC皆為格子點時,此平方必為有理數
case 1
垂足在外面
\(\min(\rm |C-A|, |C-B|)\)
case 2
垂足在線段上
用面積以及底邊長度求高!
所以到底怎麼判垂足在哪?
\(\angle CAB\) 和 \(\angle CBA\) 都是銳角時,長度為高!
用內積!
資芽給的實作
題目
-
NEOJ 400:向左轉向右轉
- 給你平面上n個點,依序走訪每一個點,試問走訪過程中共執行幾次的左轉、右轉以及迴轉。
-
CSES 2190:Line Segment Intersection
- 裸的線段相交
-
NEOJ 401:線段相交
- 裸的線段相交
極角排序
逆時鐘一個一個排
Polar coordinate system
利用每一個對特定點(通常是原點)的極座標的角度進行排序
C++ 三角函數
- 用 asin(),acos() 常常有值域的問題
- 用 atan2() 很棒!
- 值域 \((-\pi,\pi]\),可以換算成0~360度
Source: Wiki
應用問題
- 點集共線問題
- 點集能構成多少個三角形
- 常常搭配雙指針找範圍
Sort by Angle
- 優點:直觀
- 缺點:慢、浮點數誤差
- 利用\(atan2(y,x)\)求出角度
bool cmp(pt a, pt b){
return atan2(a.y,a.x) < atan2(b.y,b.x);
}
sort(p.begin(),p.end(),cmp);
Sort by Cross
- 利用有向角的正負決定前後順序
- 因為外積的夾角都是較小的角,不能直接Cross,要先判斷點在左或右半平面
bool cmp(pt a, pt b){
bool f1 = a < pt{0,0};
bool f2 = b < pt{0,0};
if(f1 != f2)return f1 < f2;
return (a ^ b) > 0;
//逆時針將點進行極角排序,從270度開始逆時針
}
sort(p.begin(),p.end(),cmp);
//以id為原點進行極角排序
TIOJ 1205 直角三角形
給你 \(N\) 個座標平面上的點\((x_i,y_i)\),請問總共可形成多少個直角三角形呢?
\(3\leq N \leq 1500, -10^9\leq x_i,y_i\leq 10^9\)
註:保證同一座標不會有兩個點
想法
- 枚舉三角形的三個頂點,排除三點共線的情況,檢查兩兩線段夾的角度是不是直角。
- 複雜度:\(O(n^3)\)
- 可能會TLE,但時限7秒應該會過
- 極角排序?
- 枚舉每一個點當作原點,對所有點進行極角排序
- 預處理將所有共線的點縮成一點
- 雙指針走訪每一個點,與原點夾90度則答案加一
- 繞完一圈就結束了!
- 複雜度:\(O(n^2\log n)\)
程式碼
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 100001
#define all(x) x.begin(),x.end()
#define eps 1e-9
#define x first
#define y second
using namespace std;
struct pt{
int x,y;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(x-b.x == 0 && y-b.y == 0)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
int operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
int operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
};
vector<pt> p,temp,pp;
vector<int> cnt;
int n,ans = 0;
bool cmp(pt a, pt b){
bool f1 = a < pt{0,0};
bool f2 = b < pt{0,0};
if(f1 != f2)return f1 < f2;
return (a ^ b) > 0;
//逆時針將點進行極角排序,從270度開始逆時針
}
//O(n)枚舉每個點當直角情況
int solve(pt id){
pp.clear();cnt.clear();temp.clear();
for(pt i : p){
pt cur = i - id;
if(cur == pt{0,0})continue;
temp.push_back(cur);
}
sort(all(temp),cmp); //以id為原點進行極角排序
pp.push_back(temp[0]); //pp每一角度只存至多一個點
cnt.push_back(1); //考慮每個點共線情況
int len = temp.size();
rep(i,1,len-1){
int cross = temp[i]^temp[i-1],dot = temp[i]*temp[i-1];
if(cross == 0 && dot >= 0)cnt[cnt.size()-1] += 1; //共線數量+=1
else {pp.push_back(temp[i]);cnt.push_back(1);} //非共線設定數量為1
}
len = pp.size(); //考慮橫跨一周的情況
rep(i,0,len-1){ //雙指針i,p1可能會超過一圈
pp.push_back(pp[i]); //將點再繞一圈
cnt.push_back(cnt[i]);
}
int ans = 0,p1 = 0;
rep(i, 0, len-1){
while(p1 < i+len && (pp[i]^pp[p1]) >= 0 && (pp[i]*pp[p1]) > 0)p1 += 1;
//夾銳角的情況要p1+=1
if((pp[i]^pp[p1]) > 0 && (pp[i]*pp[p1]) == 0)ans += cnt[i]*cnt[p1];
//正向的直角三角形,若共線則兩者數量相乘
}
return ans;
}
signed main(){
Orz;
while(cin>>n){
if(n == 0)break;
p.assign(n,{0,0});
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
int ans = 0;
rep(i,0,n-1){
ans += solve(p[i]);
}
cout<<ans<<endl;
}
}
Counting Triangles
給你二維平面上\(n\)個點(n≤2400),每一個點座標皆不相同,計算總共可以畫出多少個三角形?
$$-10^9\leq x_i,y_i\leq 10^9$$
想法
- 在任三點沒有共線的情況下,總共可以畫出 \(C_3^n\) 個三角形
- 如果有\(m\)個點共線(在其他點不共線的情況下),則會少掉 \(C_3^m\) 個三角形
- 對於每一個點進行極角排序,雙指針找共線的區間長度,總複雜度\(O(n^2\log n)\)
- 把點全部移到上半平面較好處理!
程式碼
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 100001
#define all(x) x.begin(),x.end()
#define eps 1e-9
#define x first
#define y second
using namespace std;
struct pt{
int x,y;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(x-b.x == 0 && y-b.y == 0)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
int operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
int operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
};
vector<pt> p;
vector<int> cnt;
int n,ans = 0;
bool cmp(pt a, pt b){
bool f1 = a < pt{0,0};
bool f2 = b < pt{0,0};
if(f1 != f2)return f1 < f2;
return (a ^ b) > 0;
//逆時針將點進行極角排序,從270度開始逆時針
}
//用cnt[i]統計區間長度為i的線段數量
void solve(pt id){
vector<pt> pp;
for(auto i : p){ //以id為原點
pt cur = i-id;
if(cur == pt{0,0})continue;
if(cur.y < 0){cur.x = -cur.x;cur.y = -cur.y;}
if(cur.x < 0 && cur.y==0){cur.x = -cur.x;}
pp.push_back(cur);
}
sort(all(pp),cmp); //將id當作原點進行排序
int p1 = 0,p2 = 0,len = pp.size(); //雙指針找共線區間
while(p1 < n-1){ //最大化區間
while(p2+1 < len && (pp[p1]^pp[p2+1]) == 0)p2++;
cnt[p2-p1+2]+=1;
p1 = p2+1;
}
}
signed main(){
Orz;
cin>>n;
p.assign(n,{0,0});
cnt.resize(n+1,0);
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
rep(i,0,n-1)solve(p[i]);
int ans = (n*(n-1)*(n-2))/6;
rep(i,3,n)ans-=(cnt[i]*(i-1)*(i-2))/6;
cout<<ans<<endl;
}
Max Points on a Line
給定二維平面上\(n\)個點,問一條直線最多能穿越幾個點?
$$1\leq n\leq 300$$
$$-10^4\leq x_i,y_i\leq 10^4$$
想法
- 枚舉每兩個點構成的直線,接著看有多少點會落在這條直線上
- 複雜度:\(O(n^3)\)
- 類似上一題的作法,極角排序之後用雙指針找出共線的最大區間
- 複雜度:\(O(n^2\log n)\)
- 有沒有更好的作法?
- 把每一個點當作是原點,把其他的點與原點的斜率丟進hash當中,計算同一斜率上有最多點的個數
- 斜率可以用分數處理,用浮點數有誤差的問題
- 複雜度:\(O(n^2)\)
題目
-
Leetcode 1610:Maximum Number of Visible Points
- 平面上\(n\)個點,以及一個視野的角度,找出一次最大能看到多少個點。
凸包演算法
Convex Hull