計算幾何

AaW

Who Am I

吳亞倫
鹽亞倫 / 火車倫 / AaW
建電一三學術長
因為不會計算幾何所以來當計算幾何講師
蔡嘉晉欠我鐘點費

實作懶人工具包

今天所有的實作函數～

#include <bits/stdc++.h>
using namespace std;

const double eps = 1e-7;

struct pt{
    double x, y;
    double len2() {
        return x*x + y*y;
    }

    double len() {
        return sqrt(x*x + y*y);
    }
};

// 向量相加
inline pt operator + (const pt &a, const pt &b) {
    return {a.x + b.x, a.y + b.y};
}

// 向量相減
inline pt operator - (const pt &a, const pt &b) {
    return {a.x - b.x, a.y - b.y};
} 

// 向量內積dot
inline double operator * (const pt &a, const pt &b) {
    return a.x * b.x + a.y * b.y;
}

// 向量外積cross
inline double operator ^ (const pt &a, const pt &b) {
    return a.x * b.y - a.y * b.x;
} 

// 係數積 (含交換率)
inline pt operator * (double k, const pt &a) {
    return {a.x * k, a.y * k};
} 
inline pt operator*(const pt &a, double k) {
    return {a.x * k, a.y * k};
}

// 向量除法 (係數除法)
inline pt operator / (const pt &a, double k) {
    return {a.x / k, a.y / k};
} 
 
// 誤差處理
inline bool operator == (const pt &a, const pt &b) {
    if ((abs(a.x - b.x) <= eps) && (abs(a.y - b.y) <= eps)) {
        return true;
    }
    else {
        return false;
    }
}

// 方向函數 ori，回傳OA向量轉到OB向量的方向（利用外積）
int ori (const pt &a , const pt &b, const pt &o) {
    double cross = ( a - o ) ^ ( b - o ); // OA x OB
    if( fabs( cross ) < eps ) return 0;
    return cross > 0 ? 1 : -1;
}

// 角度函數 angleType，判斷角 AOB 為銳角、直角或鈍角（利用內積）
int angleType (const pt &a, const pt &o, const pt &b) {
    double dot = ( a - o ) * ( b - o ); // OA * OB
    if( fabs( dot ) < eps ) return 0;
    return dot > 0 ? 1 : -1;
}

// k 是否在 ab 線段上
bool onseg(pt a, pt b, pt k){         
    int cross = (a - k) ^ (b - k);    // 是否平行
    int dot = (a - k) * (b - k);      // 是否在線段中
    return (cross == 0)&&(dot <= 0);
}

// 線段香蕉，線段ab是否與cd相交
bool banana(pt a, pt b, pt c, pt d){      
    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);
}

// 交點座標
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));
}


// 點到線段距離
double pointToSegmentDist (const pt &a, const pt &b, const pt &k) {
    if (angleType(a, b, k) > 0 && angleType(b, a, k) > 0) {
        // 計算垂足
        return fabs((a-k)^(b-k) / (b-a).len());
    }
    return sqrt(min( (a-k).len2(), (b-k).len2() ));
}

座標與向量

直角坐標系

$(x,\, y)$

極坐標系

$(r,\, \theta)$

向量

有向線段（向：方向、量：長度）
純量：只有大小、沒有方向的量
$\overrightarrow{AB} $ ：由 A 點起始指向 B 點的向量
可任意移動
向量長度（加絕對值表示）：兩點的歐基里德距離

$\left|\overrightarrow{AB}\right| = \sqrt{{(x_A - x_B)}^2 + {(y_A - y_B)}^2}$

向量的運算

在數學上，我們一樣會替向量定義常見的加減乘乘符號。

向量加法 $\vec a + \vec b$
向量減法 $\vec a - \vec b$
向量外積 $\vec a \times \vec b$
向量內積 $\vec a \cdot \vec b$

向量的運算

向量加法

頭尾相連：$ \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$
- $\Rightarrow \,a_2 + b_2 = c_2$

向量減法

終點減起點

加上反向向量

Ex. $\overrightarrow{AB} - \overrightarrow{CD} = \overrightarrow{AB} + \overrightarrow{DC}$

向量乘法

分為三種

向量係數積：純量 * 向量 -> 向量
向量內積：　向量 $\cdot$ 向量 -> 純量
向量外積：　向量 $\times$ 向量 -> 向量

向量係數積

若 $\vec A = (x, y)$

$k$ 為純量

則 $k\vec A = (kx, ky)$

幾何意義：

將向量長度放大 k 倍

向量內積（DOT)

定義：若 $\vec a$ 與 $\vec b$ 為二度或三度空間的向量，且 $\theta$ 為$\vec a$ 與 $\vec b$ 的夾角（滿足$0\leq\theta\leq\pi$），

則點積 $\vec a\cdot\vec b$ 定義為

a\cdot b= \begin{cases} |a||b|\cos\theta& ,a ≠ 0\land b ≠ 0\\ 0& a = 0\lor b = 0 \end{cases}

內積（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$ 為一向量定義為：

a\times b = (\left |\begin{matrix}a_2 &a_3\\b_2 &b_3 \\\end{matrix} \right | ,\left |\begin{matrix}a_3&a_1\\b_3 &b_1 \\\end{matrix} \right | ,\left |\begin{matrix}a_1&a_2\\b_1 &b_2 \\\end{matrix} \right | )

向量的外積只有在三度空間才有定義！

外積結果的向量由右手定則決定

外積（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

外積（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 {
    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}$ 來處理

向量運算要怎麼實作？

用運算子多載

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)$

三角函數

（補充）

廣義三角函數

\sin(\theta) = \frac{y}{r} \\\cos(\theta) = \frac{x}{r} \\\tan(\theta) = \frac{y}{x} \\atan2(y,x) = \theta

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

有向面積

性質

利用外積計算多邊形面積
有向角：逆時針為正、順時針為負
外積有正有負，稱為「有向面積」

可以幹嘛？

計算面積
評估角度

三角形面積

兩向量夾出的三角形有向面積為：

\triangle ABC = \frac{1}{2} \overrightarrow{AB}\times \overrightarrow{AC}

WHY?

三角形面積

幾何證明：簡單明瞭

\begin{aligned}\triangle ABC &=\frac{1}{2}|\overrightarrow{AB}\times \overrightarrow{AC}| \\&=\frac{1}{2}|\overrightarrow{AB}||\overrightarrow{AC}|\sin\theta \end{aligned}

多邊形面積

在平面上任選一點（通常是原點）
以此點為一端點，將多邊形切割成許多三角形。
依照逆時鐘的方向將所有三角形的有向面積相加即是多邊形的有向面積。

Source: https://reurl.cc/NZ5ZNq

多邊形面積

假設多邊形的頂點依序為（要特別注意結束的點必須是起始點）：

P_0,P_1,...,P_{n-1},P_n = P_0

則多邊形的面積為（aka測量師公式）：

AREA = \frac{1}{2}|\sum_{i = 0}^{n-1}{\overrightarrow{P_i}\times \overrightarrow{P_{i+1}}}|

多邊形面積

Pick's Theorem

給定頂點座標均是格子點的簡單多邊形，皮克定理說明了其面積$A$和內部格點數目$i$、邊上格點數目$b$的關係：

A = i + \frac{b}{2} - 1

證明在這裡

（補充）

Pick's Theorem

A = i + \frac{b}{2} - 1

Source: Wiki

\textcolor{#318CE7}{A = 10} \\\textcolor{red}{i = 7} \\\textcolor{Aquamarine}{B = 8}

題目

ZJ a871: Museum Area
- n個點圍成的多邊形，求面積。
CSES: Polygon Area
- 一樣是裸題
Polygon Lattice Points
- 裸皮克定理

線段相交

~~segment banana~~

二元一次方程式求焦點

如何判斷兩個直線是否香蕉？

看有沒有解

更快的做法？

直線法向量

若一直線 $ax+by = c$

則其法向量為 (a, b)

法向量方向與直線垂直

proof ? 畫看看

已知兩直線法向量，

可透過法向量是否平行得知兩向量有沒有焦點

直接求出直線交點，再判斷交點是否在線段內

如何判斷兩個線段是否香蕉？

這樣好嗎？
聽起來很麻煩誒

求交點的問題

沒有交點？
無限多交點？
交點不好求 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$異側

ori(A,B,C)\times ori(A,B,D) <= 0

＝0代表C或D有一點在直線 AB 上面

點與線段的關係

判斷線段香蕉之前，先處理點在線段上的情況
onseg() 判斷點 O 是否在線段 AB 之上
當點在線段上時，OA向量與OB向量之外積為零、內積為負（夾角為180度）

// k 是否在 ab 線段上
bool onseg(pt a, pt b, pt k){         
    int cross = (a - k) ^ (b - k);    // 是否平行
    int dot = (a - k) * (b - k);      // 是否在線段中
    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);
}

資芽版的香蕉

線段交點座標

知道兩線段是否相交之後，進一步求出交點座標

\vec{P} = \vec{A} + t\overrightarrow{AB}

其中，$t$為ACD和

線段交點座標

知道兩線段是否相交之後，進一步求出交點座標

線段交點座標

幾個算式列一列得到：

t = \frac{\overrightarrow{CA}\times\overrightarrow{CD}}{\overrightarrow{CD}\times\overrightarrow{AB}}

帶回 $\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

https://cses.fi/problemset/task/2192

給定平面上$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";
    }
}

題目

NEOJ 400：向左轉向右轉
- 給你平面上n個點，依序走訪每一個點，試問走訪過程中共執行幾次的左轉、右轉以及迴轉。
CSES 2190：Line Segment Intersection
- 裸的線段相交
NEOJ 401：線段相交
- 裸的線段相交

點到線段距離

問 C 點到 $\overline{\rm AB}$ 的距離平方

當ABC皆為格子點時，此平方必為有理數

case 1

垂足在外面

$\min(\rm |C-A|, |C-B|)$

case 2

垂足在線段上

用面積以及底邊長度求高！

所以到底怎麼判垂足在哪？

$\angle CAB$ 和 $\angle CBA$ 都是銳角時，長度為高！

用內積！

實作

// 角度函數 angleType，判斷角 AOB 為銳角、直角或鈍角（利用內積）
int angleType (const pt &a, const pt &o, const pt &b) {
    double dot = ( a - o ) * ( b - o ); // OA * OB
    if( fabs( dot ) < eps ) return 0;
    return dot > 0 ? 1 : -1;
}

// 點到線段距離
double pointToSegmentDist (const pt &a, const pt &b, const pt &k) {
    if (angleType(a, b, k) > 0 && angleType(b, a, k) > 0) {
        // 計算垂足
        return fabs((a-k)^(b-k) / (b-a).len());
    }
    return sqrt(min( (a-k).len2(), (b-k).len2() ));
}

資芽給的實作

極角排序

逆時鐘一個一個排

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 operator < (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());
//以id為原點進行極角排序

TIOJ 1205 直角三角形

https://tioj.ck.tp.edu.tw/problems/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

https://leetcode.com/problems/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

凸包是什麼

凸
- 凸包是一個凸多邊形
- 每個內角 $< 180^{\circ}$
包
- 凸包要包住所有平面上的點

如何找出凸包？

Gift wrapping (a.k.a. Jarvis march)
Graham scan
Quickhull
Divide and conquer
Monotone chain (a.k.a. Andrew's algorithm)
Incremental convex hull algorithm
Kirkpatrick–Seidel algorithm
Chan's algorithm

一堆凸包演算法

偷

下課

學弟妹都來考學術

演算法社課計算幾何（Ver. 一五）

By Aaron Wu

計算幾何

Who Am I

目錄

實作懶人工具包

座標與向量

向量

向量的運算

向量係數積

向量內積 （DOT)

內積（DOT）

內積（DOT）

外積（CROSS）

外積（CROSS）

外積（CROSS）

比較

向量怎麼實作？

向量運算要怎麼實作？

Eps是什麼？

三角函數

廣義三角函數

弧度

三角函數運算

有向面積

性質

可以幹嘛？

三角形面積

三角形面積

多邊形面積

多邊形面積

多邊形面積

Pick's Theorem

Pick's Theorem

題目

線段相交

方向函式

方向函式

更精準的求線段相交

點與線段的關係

線段相交完整程式碼

線段交點座標

線段交點座標

線段交點座標

動腦一下

Point in Polygon

來看看凸多邊形的情況

射線作法

程式碼

題目

點到線段距離

點到線段距離

實作

極角排序

Polar coordinate system

C++ 三角函數

應用問題

Sort by Angle

Sort by Cross

TIOJ 1205 直角三角形

想法

程式碼

Counting Triangles

想法

程式碼

Max Points on a Line

想法

題目

凸包演算法

凸包是什麼

如何找出凸包？

下課

演算法社課 計算幾何（Ver. 一五）

More from Aaron Wu

向量內積（DOT)

演算法社課計算幾何（Ver. 一五）