Geometry

計算幾何

Introduction to Competitive Programming

2022/06/01

浮點數誤差
向量與內外積
面積
各種計幾問題
凸包
旋轉卡尺
- 最遠點對
掃描線
- 矩形覆蓋面積

浮點數誤差

可以嘗試跑以下程式碼




double x = 0.0;
while(x != 1.0){
    x += 0.1;
}

理論上在迴圈裡跑加 10 次就會結束了
但實際上執行會因為精確度誤差造成程式會停不下來

為了解決浮點數誤差的問題，通常會設一個很小的數字例如\(10^{-9}\)
在比較大小關係時，如果在誤差範圍內則視為相等

改成當前的值與目標若差值小於 EPS 則視為相同







#define EPS 1e-9
double x=0.0;
while(abs(x-1.0) > EPS){
    x += 0.1;
}
// 控制輸出到小數點後 10 位
cout << fixed << setprecision(10) << x << endl;

而 EPS 大小則視不同題目而定

不同題目都有不同的精度限制

由於浮點數誤差是個很麻煩的問題
因此在寫計算幾何時，如果能避免使用浮點數則會盡量使用整數型態

像是判斷距離，當前點到點 u, v 哪個比較近
距離公式 \(\sqrt{x^2+y^2}\) 會開根號
但實際上只需求出 \(x^2+y^2\) 即可判斷兩者距離大小
只需使用 long long 型態即可

歐基里德距離

\(d(p_1, p_2) = \sqrt{(x_1-x_2)^2+(y_1-y_2)^2)}\)

曼哈頓距離

\(d(p_1, p_2) = |x_1-x_2|+|y_1-y_2|\)

向量

\(\vec{v} = (a, b)\)，代表往 +x 軸走單位 \(a\)，往 +y 軸走單位 \(b\)

長度 \(|\vec{v}|\) 為 \(\sqrt{a^2+b^2}\)

向量四則運算

假設有兩向量 \(v_{1} = (x_1, y_1), v_{2} = (x_2, y_2)\)

\(v_1 + v_2 = (x_1 + x_2, y_1 + y_2)\)
\(v_1 - v_2 = (x_1 - x_2, y_1 - y_2)\)
\(v_1 × k = (x_1 × k, y_1 × k)\) (向量伸長成 \(k\) 倍)
\(v_1 ÷ k = (x_1 ÷ k, y_1 ÷ k)\) (向量縮短成 \(k\) 倍)

內積

\(\vec{v_{1}} =(x_1, y_1), \vec{v_{2}}=(x_2, y_2)\)
\(\vec{v_{1}} \cdot \vec{v_{2}} = x_1x_2 + y_1y_2 = |\vec{v_1}||\vec{v_2}| cos \theta\) \((\theta\) 為兩向量夾角)

兩向量同向，則內積為正
兩向量反向，則內積為負
兩向量垂直，則內積為 0

外積

\(\vec{v_{1}} =(x_1, y_1), \vec{v_{2}}=(x_2, y_2)\)
\(\vec{v_{1}} \times \vec{v_{2}} = x_1y_2 -x_2y_1=|\vec{v_1}||\vec{v_2}| sin \theta\) \((\theta\) 為兩向量夾角)

\(\vec{v_{1}} \times \vec{v_{2}} = - \vec{v_{2}} \times \vec{v_{1}}\) 具有負交換律

如果 \(\vec{v_{1}}\) 到 \(\vec{v_{2}}\) 為逆時鐘方向，則外積為正
如果 \(\vec{v_{1}}\) 到 \(\vec{v_{2}}\) 為順時鐘方向，則外積為負
兩向量平行則外積為 0

同時外積取絕對值(| \(\vec{v_{1}} \times \vec{v_{2}}\) |)等同於兩向量夾起來的平行四邊形面積

儲存方法

點 (儲存座標 x, y 表示)

struct Pt{ 
    ld x, y; 
};

線段 (儲存線段兩端點儲存)

struct Line{ 
    Pt st, ed; 
};

圓 (儲存圓心與半徑)




struct Circle{ 
    Pt o; // 圓心
    ld r; // 半徑
};

多邊形 ( 用 vector 儲存所有點)




struct poly{
    int n; // n 邊形
    vector<Pt> pts; 
};

default code



























#define ld long double
struct Pt {
  ld x, y;
  Pt(ld _x=0, ld _y=0):x(_x), y(_y) {}
  
  Pt operator+(const Pt &a){
    return Pt(x+a.x, y+a.y);
  }
  Pt operator-(const Pt &a){
    return Pt(x-a.x, y-a.y);
  }
  Pt operator*(const ld &a){
    return Pt(x*a, y*a);
  }
  Pt operator/(const ld &a){
    return Pt(x/a, y/a);
  }
  ld operator*(const Pt &a){ //程式碼中內積通常用*表示
    return x*a.x + y*a.y;
  }
  ld operator^(const Pt &a){ //程式碼中外積通常用^表示
    return x*a.y - y*a.x;
  }
  bool operator<(const Pt &a) const {
    return x < a.x || (x == a.x && y < a.y);
  }
};

三角形面積

海龍公式

\(s = \frac{a+b+c}{2}\) (a,b,c 為三角形三邊長)

\(\triangle abc = \sqrt{s(s-a)(s-b)(s-c)}\)
外積求法

|\(\vec{v_1} \times \vec{v_2}\)| 為所夾的平行四邊形面積
而平行四邊形的一半為三角形面積

因此 \(\frac{|\vec{v_1} \times \vec{v_2}|}{2}\) 為三角形面積

多邊形面積

一個 \(n\) 個點的多邊形(\(P_1, P_2,... P_n\))
可以拆成很多三角形，再分別套用外積公式

\(\frac{1}{2}|\sum\limits_{i=1}^{n} \overrightarrow {OP_{i}} \times \overrightarrow {OP_{i+1}}|\)

複雜度 \(O(N)\)

使用外積公式求面積會發現

\(\frac{1}{2}|\sum\limits_{i=1}^{n} \overrightarrow {OP_{i}} \times \overrightarrow {OP_{i+1}}|\)

如果所有點都為整數點，則除了最後乘 \(\frac{1}{2}\) 使得面積有可能為小數，否則在用外積公式時都為整數運算

因此整數點的面積不是整數就是 .5 結尾
過程中可以完全不用用到小數可以避免浮點數誤差

是否三點共線

如果 3 個點連成一點線，則圍起來的三角形面積為 0
則直接判斷外積是否為 0 即可



bool collinearity(const Pt& a, const Pt& b, const Pt& c){
    return (b-a)^(c-a) < EPS;
}

判斷點是否在線段上

必要條件為點與線段共線(外積為 0)
且點往兩端點的方向會相反(內積為負)



bool inLine(const Pt& p, const Line& li){
    return collinearity(li.st, li.ed, p) && (li.st-p)*(li.ed-p) < EPS;
}

判斷兩線段是否相交

如果相交則分成兩種 case

交點在其中一個端點上
交點不在端點上

第一種 case 直接判斷以下四種情況即可

\(\texttt{inLine(line1.st, line2)}\)
\(\texttt{inLine(line1.ed, line2)}\)
\(\texttt{inLine(line2.st, line1)}\)
\(\texttt{inLine(line2.ed, line1)}\)

第二種 case :交點不在端點上

圖會長成以下

呈現交叉情況

則點 A、B 會在線段 \(\overline{CD}\)相異側
則點 C、D 會在線段 \(\overline{AB}\)相異側

如果在相異側則外積相乘為負

\((\overrightarrow{AB}\times\overrightarrow{AC})(\overrightarrow{AB}\times\overrightarrow{AD}) <0\)
\((\overrightarrow{CD}\times\overrightarrow{CA})(\overrightarrow{CD}\times\overrightarrow{CB}) <0\)

點在多邊形內部

兩種判法

若點在多邊形內，則選一個方向的射線出現會碰到奇數次邊
(需要特判點是否在多邊形的邊或頂點上)
並且選的射線隨機選，避免剛好交到多邊形的頂點上
若點在凸多邊形內部，則點與多邊形上所有頂點 \(P_1,P_2,...P_n\) 依序對於所有 \(\overrightarrow{AP_i} \times \overrightarrow{AP_{i+1}}\) 皆為同號(全部都是正或者全部都是負的)
若有異號的情況則代表點在凸包外
若出現一個 0 則代表點在邊上，出現兩個 0 則代表點在頂點上

凸包

你有 \(n\) 個點，然後你要找出一個凸多邊形可以圍住這n個點且面積最小，
這個凸多邊形稱為凸包

你也可以想成凸包是用一條橡皮筋把牆壁上的圖釘圈起來的範圍

Andrew's Monotone Chain

(競賽常用的做法)

作法為分別求出上下凸包，求完之後再和再一起，
上半凸包+下半凸包=完整的凸包

先求下凸包

先將全部點照 \(x\) 大小排序,再排 \(y\) 大小
依序嘗試將點加入凸包判斷是否合法
不合法刪掉前面最後一個點，合法則將點加進凸包

如何判斷合法

首先先觀察凸包，若我們從左到右依序看點，會發現下一個向量都是在前一個向量的逆時鐘方向
像是 \(\overrightarrow{P_2 P_3}\) 會在 \(\overrightarrow{P_2 P_4}\) 的逆時鐘方向

而 \(\overrightarrow{P_{i} P_{i+1}}\) 若是在\(\overrightarrow{P_{i} P_{i+2}}\) 的逆時鐘方向，則\(\overrightarrow{P_{i} P_{i+1}} \times \overrightarrow{P_{i} P_{i+2}} > 0\)

因此如果 \(\overrightarrow{P_{i} P_{i+1}} \times \overrightarrow{P_{i} P_{i+2}} \le 0\) 代表前一個點會在凸包裡面，
而不會在凸包上，就可以把前一個點刪掉，再把新的點加進去凸包裡面

而實作可以用一個stack去維護，
一開始先把前兩個點直接加進去stack之後每一個點\(P_i\)都去判斷，
向量是不是往逆時針方向 ( \(\overrightarrow{P_{i-2} P_{i-1}} \times \overrightarrow{P_{i-2} P_{i}} > 0\) )
若非法則把stack的top刪掉，一直刪到合法為止，
再把新的點加進去凸包





































struct Pt{
    int x,y;
    Pt(){}
    Pt(int _x,int _y){
        x=_x,y=_y;
    }
    Pt operator-(const Pt &a){
        return Pt(x-a.x, y-a.y);
    }
    bool operator<(const Pt &a) const {
        return x < a.x || (x == a.x && y < a.y);
    }
    friend int cross(const Pt& o,const Pt& a,const Pt& b){
        Pt lhs = o-a, rhs = o-b;
        return lhs.x*rhs.y - lhs.y*rhs.x;
    }
};
vector<Pt> convex_hull(vector<Pt> hull){
    sort(hull.begin(),hull.end());
    int top=0;
    vector<Pt> stk;
    for(int i=0;i<hull.size();i++){
        while(top>=2&&cross(stk[top-2],stk[top-1],hull[i])<=0)
            stk.pop_back(),top--;
        stk.push_back(hull[i]);
        top++;
    }
    for(int i=hull.size()-2,t=top+1;i>=0;i--){
        while(top>=t&&cross(stk[top-2],stk[top-1],hull[i])<=0)
            stk.pop_back(),top--;
        stk.push_back(hull[i]);
        top++;
    }
    stk.pop_back();
    return stk;
}

排序 \(O(NlgN)\) + 找凸包 \(O(N)\)
複雜度 \(O(NlgN)\)

旋轉卡尺

最遠點對
最大三角形
最大四角形

最遠點對

一群點當中，距離最遠的兩個點叫作「最遠點對」

觀察一下可以發現，因為凸包上所有的點是包圍所有點的多邊形，因此最遠的兩點一定是凸包的其中兩點。

因此我們在做最遠點對的時候要先做凸包，接著可以用旋轉卡尺的概念去找最遠點對

類似雙指針的概念
每次移動把一個指針往下移一格，另一個則不斷往下移動，直到最大為止










double FarthestPair(vector<Pt> arr){
    double ret=0;
    for(int i=0,j=i+1;i<arr.size();i++){
        while(distance(i,j)<distance(i,(j+1)%arr.size())){
            j=(j+1)%arr.size();
        }
        ret=max(ret,distance(i,j));
    }
    return ret;
}

複雜度 \(O(NlgN)\)
求出凸包 \(O(NlgN)\) + 旋轉卡尺 \(O(N)\)

最大三角形

給一堆點，求其中三個點圍成的三角形面積最大

一樣用旋轉卡尺，窮舉每個點當定點，第二個點每次向下轉一格，
而第三個點就跟著轉到找到最大為止
複雜度 \(O(N^2)\) 窮舉點對

最大四角形

給一堆點，求其中四個點圍成的四角形面積最大

作法：把四角形拆成兩個三角形，兩邊分別旋轉卡尺

複雜度 \(O(N^2)\)

掃描線

解決圖形的周長，面積總和等問題

在二維座標中，給定n個矩形，求出這些矩形覆蓋的面積總和

作法

先將矩形的y座標離散化並且把矩形拆成左界右界
掃描線由左至右掃過去
用線段樹維護區間
每次掃到一個邊界就計算答案並加減值

假設我們掃描線從x軸由左至右掃過去

在區間10~40加值

ans += 30 * (30-10)
(y軸覆蓋總和) * (掃描線現在位置-前一個位置)
在區間30-70加值

ans += 60 * (50-30)
(y軸覆蓋總和) * (掃描線現在位置-前一個位置)
在區間10-40減值

ans += 40 * (80-50)
(y軸覆蓋總和) * (掃描線現在位置-前一個位置)
在區間30-70減值

https://oi-wiki.org/geometry/scanning/

計算幾何是個噁心的單元
除了實作複雜還須調整浮點數，
常常程式寫對但因為浮點數誤差可能就要處理個半小時到幾個小時

但卻又是 ICPC 台灣區常見的題型一定要好好跟他混熟
並且有很多很噁心的內容沒交到有興趣的可以在自己去研究
https://cp-algorithms.com/#geometry

Question Time and Practice

Homework Link

提交期限到下星期三 6/8 23:59 截止

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`	在筆記中貼入程式碼
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.