计算几何

导引问题

求三边形面积：

• 在解析几何里，\(\triangle\text{ABC}\)的面积可以通过如下方法求得：

• 点坐标 \(\Rightarrow\) 边长 \(\Rightarrow\) 海伦公式 \(\Rightarrow\) 面积

• 缺点？

• 计算量大、精度损失，更好的方法？

边的属性

题目：从三个问题谈起

1、给定具有共同端点的两个有向线段\(P_0P_1\)和\(P_0P_2\) ，请问——\(P_0P_1\)是否在\(P_0P_2\)的顺时针方向上？

2、给定两个有向线段\(P_0P_1\)和\(P_0P_2\)，如果我们首先沿着\(P_0P_1\)的方向走，然后沿着\(P_1P_2\)走，请问——在\(P_1\)点，是左转还是右转？

3、给定两个向量\(P_1P_2\)和\(P_3P_4\)，如何判断它们是否相交？

叉乘

问题：是否在顺时针方向上：

转向

交点

方法一 使用解析几何：解方程组+高斯消元，精度低

方法二利用叉乘：

x,y 都需要判断，这样可以保证垂直，水平下的特殊情况

交点代码

太疯狂了，总共写了3个多小时 gxx’s Problem

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define ll long long
#define inf 0x3f3f3f3f

int sign(ll x){
    if(x>0){
        return 1;
    }
    else if(x<0){
        return -1;
    }
    return 0;
}

struct frac{
    ll z;
    ll m;
    ll gcd(ll a,ll b){
        while(b){
            ll r=a%b;
            a=b;
            b=r;
        }
        return a;
    }
    void create(ll a,ll b=1){
        ll t=gcd(a,b);
        z=a/t;
        m=b/t;
        if(m<0){
            m=-m; z=-z;
        }
    }
    void print(){
        if(m==1){
            printf("%lld",z);
        }   
        else{
            printf("%lld/%lld",z,m);
        }
    }
};

struct point{
    ll x;
    ll y;
    bool operator ==(const point &p){
        return x==p.x && y==p.y;
    }
    ll cross(const point &p1,const point &p2){
        return (p1.x-x)*(p2.y-y)-(p2.x-x)*(p1.y-y);
    }
};

struct segment{
    point s;
    point t;
    bool ispoint(){
        return s==t;
    }
    ll cross(point &p){
        return p.cross(s,t);
    }
    bool has(point &p){
        if(cross(p)!=0){
            return false;
        }
        if(s.x==t.x){
            return p.y<=max(s.y, t.y) 
                && p.y >= min(s.y, t.y);
        }
        else{
            return p.x <= max(s.x, t.x) 
                && p.x >= min(s.x, t.x);
        }
    }
    bool operator ==(const segment &seg){
        return (s==seg.s && t==seg.t)  (s==seg.t && t==seg.s);
    }
};

int intersection(segment &p1, segment &p2, frac &x, frac &y){
    //0 none
    //1 one point
    //-1 inf
    bool isp1=p1.ispoint();
    bool isp2=p2.ispoint();
    if(isp1 && isp2){
        if(p1==p2){
            //When 2 segments are the same point.
            x.create(p1.s.x);
            y.create(p1.s.y);
            return 1;
        }
        else{
            return 0;
        }
    }
    else if(isp1  isp2){
        if(isp1 && p2.has(p1.s)){
            x.create(p1.s.x);
            y.create(p1.s.y);
            return 1;
        }
        else if(isp2 && p1.has(p2.s)){
            x.create(p2.s.x);
            y.create(p2.s.y);
            return 1;
        }
        else
            return 0;
    }
    else{
        ll a,b,c,d;
        int s1,s2,s3,s4;
        s1 = sign(a = p2.cross(p1.s));
        s2 = sign(b = p2.cross(p1.t));
        s3 = sign(c = p1.cross(p2.s));
        s4 = sign(d = p1.cross(p2.t));
        if((s1^s2) == -2 && (s3 ^ s4) == -2){
            x.create(b*p1.s.x - a*p1.t.x, b-a);
            y.create(b*p1.s.y - a*p1.t.y, b-a);
            return 1;
        }
        else if(a==0 && b==0){ //overlay
            if(p2.has(p1.s) && p2.has(p1.t)){
                return -1;
            }
            else if(p2.has(p1.s)){//p1's on p2. Here are several cases
                if((p1.s==p2.s && !p1.has(p2.t)) 
                     (p1.s==p2.t && !p1.has(p2.s))){
                    //When p1's == p2.s or t but p2 is not on p1
                    x.create(p1.s.x);
                    y.create(p1.s.y);
                    return 1;
                }
                else
                    return -1;
            }
            else if(p2.has(p1.t)){
                if((p1.t==p2.t && !p1.has(p2.s) ) 
                     (p1.t == p2.s && !p1.has(p2.t) )){
                    // When p1't == p2.s or t but p2 is not on p1
                    x.create(p1.t.x);
                    y.create(p1.t.y);
                    return 1;
                }
                else{
                    return -1;
                }
            }
            else{
                if(p1.has(p2.s)){
                    return -1;
                }
                else{
                    return 0;
                }
            }
        }
        else if(a==0){//By default, b!=0
            if(p2.has(p1.s)){
                x.create(p1.s.x);
                y.create(p1.s.y);
                return 1;
            }
            else{
                return 0;
            }
        }
        else if(b==0){
            if(p2.has(p1.t)){
                x.create(p1.t.x);
                y.create(p1.t.y);
                return 1;
            }
            else{
                return 0;
            }
        }
        else if(c==0){
            if(p1.has(p2.s)){
                x.create(p2.s.x);
                y.create(p2.s.y);
                return 1;
            }
            else{
                return 0;
            }
        }
        else if(d==0){
            if(p1.has(p2.t)){
                x.create(p2.t.x);
                y.create(p2.t.y);
                return 1;
            }
            else {
                return 0;
            }
        }
        else{
            return 0;
        }
    }
}

int main(){
    int t;
    cin>>t;
    frac x,y;
    int flag=0;
    segment ab,cd;
    while(t--){
        scanf("%lld %lld %lld %lld",&ab.s.x,&ab.s.y,&ab.t.x,&ab.t.y);
        scanf("%lld %lld %lld %lld",&cd.s.x,&cd.s.y,&cd.t.x,&cd.t.y);
        flag= intersection(ab,cd,x,y);
        if(flag==1){
            printf("1\n");
            x.print(),printf(" "),y.print();
            printf("\n");
        }
        else if(flag==0){
            printf("0\n");
        }
        else {
            printf("INF\n");
        }
    }
    return 0;
}

叉乘：没有除法和开根号，避免传统方法的精度问题

多边形面积

三角形面积：叉乘二分之一

abs绝对值

凸多边形分割为三角形

凹多边形面积

先判断凹凸性：依次相叉乘

公式也成立（忽略符号的情况下）

依然成立！！！

多边形面积公式：\(A = \sum (A_i) \ (i=1 \cdots N-2)\)

结论：“有向面积”\(A\)比“面积”\(S\)其实更本质！

中间点

外部点

公式依然成立

原点剖分（荐

方便直接可用

简化的公式

\[ A = \frac{1}{2} \sum_{i=1}^{N} \begin{vmatrix} X_i & Y_i \\ X_{i+1} & Y_{i+1} \end{vmatrix} \quad (i = 1 \dots N) \]

多边形重心

求多边形重心

• 给定一个简单多边形，求其重心。
• 输入：多边形（顶点按逆时针顺序排列）
• 输出：重心点C的坐标
• 三角形的重心公式是？
• \(x_0=(x_1+x_2+x_3)/3\)
• \(y_0=(y_1+y_2+y_3)/3\)

处理方式将处理成三角形重心，然后加权

凸包

围绳子

凸包上的点一定是原来就有的点，凸多边形

Graham Scan扫描法

选定一个y最低的点，按与X轴的角度sort排序叉乘比较（p1,p2）叉乘比较相对大小

判断左右转

凸包代码

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <queue>
#include <stack>
using namespace std;
#define ll long long
#define inf 0x3f3f3f3f

struct node{
    ll x;
    ll y;
}tree[10010];

void swapnode(node &p1,node &p2){
    node p3=p1;
    p1=p2;
    p2=p3;
    return ;
}

ll cross(node tree,node p1,node p2){
    return (p1.x-tree.x)*(p2.y-tree.y)-
        (p2.x-tree.x)*(p1.y-tree.y);
}

bool cmp(node p1,node p2){
    ll sum=cross(tree[0],p1,p2);
    if(sum>0 ){
        return true;
    }
    return false;
}

double dis(node p1,node p2){
    return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}

int main(){
    int n;
    while(cin>>n,n){
        double ans=0;
        memset(tree,0,sizeof(tree));
        int pos=0;
        ll miny=inf;
        for(int i=0;i<n;i++){
            ll x,y;
            scanf("%lld%lld",&x,&y);
            tree[i].x=x;
            tree[i].y=y;
            if(y<miny){
                miny=y;
                pos=i;
            }
        }
        if(n==1){
            cout<<0<<endl;
            continue;
        }
        else if(n==2){
            printf("%.2lf\n",dis(tree[0],tree[1]));
            continue;
        }
        swapnode(tree[0],tree[pos]);
        tree[n]=tree[0];
        sort(tree+1,tree+n,cmp);
        // //testsort:::
        // for(int i=0;i<n;i++){
        //  cout<<tree[i].x<<" "<<tree[i].y<<endl;
        // }
        // cout<<"END"<<endl;

        stack<node> s;
        s.push(tree[0]);
        s.push(tree[1]);
        for(int i=2;i<=n;i++){
            node p3=tree[i];
            node p2=s.top();
            s.pop();
            node p1=s.top();
            s.pop();
            bool flag=false;
            // printf("shuzi::%lld %lld %lld %lld %lld %lld\n",p1.x,p1.y,p2.x,p2.y,p3.x,p3.y);
            // cout<<"cross::"<<cross(p2,p3,p1)<<endl;
            while(cross(p2,p3,p1)<0){
                if(s.empty()){
                    s.push(p1);
                    s.push(p3);
                    flag=true;
                    break;
                }
                p2=p1;
                p1=s.top();
                s.pop();
            }
            if(flag){
                continue;
            }
            s.push(p1);
            s.push(p2);
            s.push(p3);
        }
        //testSTACK:::
        // while(s.size()!=0){
        //  cout<<s.top().x<<" "<<s.top().y<<endl;
        //  s.pop();
        // }

        while(s.size()!=1){
            // cout<<s.top().x<<" "<<s.top().y<<endl;
            node p1=s.top();
            s.pop();
            node p2=s.top();
            ans+=dis(p1,p2);
        }
        printf("%.2lf\n",ans);
    }
    return 0;
}

伪代码

procedure GrahamScan(set Q)
    let p₀ be the point with the minimum y-coordinate
    let <p₁, ..., pₘ> be the remaining points in Q, sorted by the angle in counterclockwise order around p₀
    Top(S) ← 0
    Push(p₀, S); Push(p₁, S); Push(p₂, S)
    for i ← 3 to m do
        while the angle formed by points NextToTop(S), Top(S) and pᵢ makes a non-left turn do
            Pop(S)
        Push(pᵢ, S)
    return S

Jarvis March步进法

如同拉一根长的绳子，每次看角度最小的

计算几何卡精度，需要选一个人进行操作，考精度和运气，一般放在最后