Computational geometry (2)

Computational geometry （2）

2022 year 4 month 18 Japan 21:31:04

Find convex hull （ Surround all points The smallest perimeter , The area is not necessarily the smallest The polygon of ） The algorithm of .

Graham Algorithm

Andrew Algorithm $O(nlogn)$

1. Sort the points ,x For the first keyword ,y For the second keyword . bottleneck
2. Maintain the upper half of the convex polygon from left to right , Then maintain the lower part from right to left .（ Use Stack maintain ）
3. For a clockwise convex polygon , First add the first two points to the stack , Judge from the third point , Judge the relationship between the extension line of the straight line composed of the current point and the midpoint of the stack , If the point to be added is to the left of the extension line of the straight line （ Anti-clockwise ）, Push stack out of stack , Then judge , Until the point to be added is on the right of the straight extension , Then add the point to the stack , Finally, update the stack with the starting point .

Try rolling one by yourself .

2022 year 4 month 19 Japan 14:28:13

Hand rolling failed , Study y Master code .

#include<iostream>
#include<algorithm>
#include<cstring>
#include<vector>
#include<cstdlib>
#include<stack>
#include<map>
#include<queue>
#include<cmath>
#include<set>
#include<cassert>
#define fo(i,a,n) for(int i=(a);i<=(n);i++)
#define fo_(i,a,n) for(int i=(a);i<(n);i++)
#define debug(x) cout<<#x<<":"<<x<<endl
#define all(x) (x).begin(),(x).end()
#define pb push_back
#define endl '\n'
using namespace std;
const int N = 1e5+10;
const int M = 510*510;
const double inf = 1e8;
typedef long long ll;
#define xx first
#define yy second

typedef pair<double,double> PDD;
typedef pair<double,double> Vector;
PDD a[N];
int n,stk[N],top;
bool used[N];

double cross(double x1,double y1,double x2,double y2){
    // The first dot crosses the second dot 
    return x1*y2-x2*y1;
}
double dis(Vector A,Vector B){
    
    return sqrt((A.xx-B.xx)*(A.xx-B.xx)+(A.yy-B.yy)*(A.yy-B.yy));
}
int main(){
    
    cin>>n;
    fo(i,1,n){
    
        cin>>a[i].xx>>a[i].yy;
    }

    sort(a+1,a+n+1,[&](PDD a,PDD b){
    
        if(a.xx==b.xx)return a.yy>b.yy;
        return a.xx<b.xx;
    });
    
    stk[++top]={
    1};
    stk[++top]={
    2};
    used[1]=used[2]=1;
    fo(i,3,n){
    
        Vector C = a[i];    
        while(top>1){
    //  There must be at least one point in the convex hull set .
            Vector B = a[stk[top]];
            Vector A = a[stk[top-1]]; 
            if(cross(C.xx-A.xx,C.yy-A.yy,B.xx-A.xx,B.yy-A.yy) < 0){
    
                used[top]=false;
                top--;
            }
        }
        stk[++top]=i;
        used[i]=true;
        break;
    }
    
    for(int i=n;i>=2;i--){
    
        if(used[i])continue;
        Vector C = a[i];    
        while(top>1){
    //  There must be at least one point in the convex hull set .
            Vector B = a[stk[top]];
            Vector A = a[stk[top-1]]; 
            if(cross(C.xx-A.xx,C.yy-A.yy,B.xx-A.xx,B.yy-A.yy) < 0){
    
                used[top]=false;
                top--;
            }
        }
        stk[++top]=i;
        used[i]=true;
        break;
    }
    //  The first point of special judgment 
    Vector C = a[1];   
    while(1){
    
        if(top>1){
    //  There must be at least one point in the convex hull set .
            Vector B = a[stk[top]];
            Vector A = a[stk[top-1]]; 
            if(cross(C.xx-A.xx,C.yy-A.yy,B.xx-A.xx,B.yy-A.yy) < 0){
    
                used[top]=false;
                top--;
            }
        }
        stk[++top]=1;
        used[1]=true;
        break;
    }
    
    
    // fo(i,1,n){
    
    // if(used[i]){
    
    // cout<<i<<" ";
    // }
    // }
    // cout<<endl;
    
    bool flag=0;
    int last,first;
    double res=0;
    fo(i,1,n){
    
        if(!used[i])continue;
        if(!flag){
    
            flag=1;
            last = i;
            first = i;
            continue;
        }
        res+=dis(a[i],a[last]);
        last = i;
    }
    res+=dis(a[last],a[first]);
    printf("%.2lf",res);
    return 0;
}

Where you can learn .

stk[++top]={
    1};
stk[++top]={
    2};
used[1]=used[2]=1;
fo(i,3,n){
    
    Vector C = a[i];    
    while(top>1){
    //  There must be at least one point in the convex hull set .
        Vector B = a[stk[top]];
        Vector A = a[stk[top-1]]; 
        if(cross(C.xx-A.xx,C.yy-A.yy,B.xx-A.xx,B.yy-A.yy) < 0){
    
            used[top]=false;
            top--;
        }
    }
    stk[++top]=i;
    used[i]=true;
    break;
}
----------------------------------------------
y The general way of writing ：
    for (int i = 0; i < n; i ++ )
    {
    
        while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) <= 0)
        {
    
            //  Even if the point on the convex hull boundary is deleted from the stack , Nor can it be deleted used The mark on the wall 
            if (area(q[stk[top - 1]], q[stk[top]], q[i]) < 0)
                used[stk[top -- ]] = false;
            else top -- ;
        }
        stk[ ++ top] = i;
        used[i] = true;
    }

What you learned ： At present, several can be directly added , But it needs to be used while When circulating special judgment , You can refer to y The total , Put all the code in for Inside the loop ,while Add another judgment .

    used[0] = false;
    for (int i = n - 1; i >= 0; i -- )
    {
    
        if (used[i]) continue;
        while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) <= 0)
            top -- ;
        stk[ ++ top] = i;
    }

You need to use the 1 The added points update the entire convex hull , Very clever use used[0] = 0;

such ,stk Will appear twice 0 Number point , Last stk[top] It must be equal to 0 , It's just convenient to calculate the perimeter .

    double res = 0;
    for (int i = 2; i <= top; i ++ )
        res += get_dist(q[stk[i - 1]], q[stk[i]]);
    return res;

 Computational geometry templates 
double get_dist(PDD a, PDD b)
{
    
    double dx = a.x - b.x;
    double dy = a.y - b.y;
    return sqrt(dx * dx + dy * dy);
}
double cross(PDD a, PDD b)
{
    
    return a.x * b.y - a.y * b.x;
}
// a->b  Vector and  a->c  Find the cross product of vectors .
double area(PDD a, PDD b, PDD c)
{
    
    return cross(b - a, c - a);//  First reload two points and subtract 
}

Refer to the big guy's code

 Can pass acwing and luogu（ No visual inspection bug）, There is no need to judge the problem of collinearity 
 Examples ：
    2
    0 0 
    0 1
 Output ：
    2

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>

#define x first
#define y second

using namespace std;

typedef pair<double, double> PDD;
const int N = 10010;

int n;
PDD q[N];
int stk[N];

double get_dist(PDD a, PDD b)
{
    
    double dx = a.x - b.x;
    double dy = a.y - b.y;
    return sqrt(dx * dx + dy * dy);
}

PDD operator-(PDD a, PDD b)
{
    
    return {
    a.x - b.x, a.y - b.y};
}

double cross(PDD a, PDD b)
{
    
    return a.x * b.y - a.y * b.x;
}

double area(PDD a, PDD b, PDD c)
{
    
    return cross(b - a, c - a);
}

double andrew()
{
    
    sort(q, q + n);
    int top = 0;
    for (int i = 0; i < n; i ++ )
    {
    
        while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) <= 0)
            top -- ;
        stk[ ++ top] = i;
    }
    int k = top;
    for (int i = n - 1; i >= 0; i -- )
    {
    
        //  When you don't want to delete the convex hull , The point on the far right  ,  The first point on the left can still be added again when finding the convex hull stk in .
        while (top > k && area(q[stk[top - 1]], q[stk[top]], q[i]) <= 0)
            top -- ;
        stk[ ++ top] = i;
    }

    double res = 0;
    for (int i = 2; i <= top; i ++ )
        res += get_dist(q[stk[i - 1]], q[stk[i]]);
    return res;
}

int main()
{
    
    scanf("%d", &n);
    for (int i = 0; i < n; i ++ ) scanf("%lf%lf", &q[i].x, &q[i].y);
    double res = andrew();
    printf("%.2lf\n", res);

    return 0;
}

Credit card bumps

Spinning bug I've been looking for it for a long time qwq

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

typedef pair<double,double>PDD;
#define xx first
#define yy second
#define pi acos(-1)
const int N = 1e5+10;

int n,stk[N],top,idx;
double a,b,r;
PDD A[N];//  Turn a rectangular square into four corners .

double dis(PDD a,PDD b){
    
    double dx = a.xx-b.xx;
    double dy = a.yy-b.yy;
    return sqrt(dx*dx+dy*dy);
}

PDD operator - (PDD a,PDD b){
    
    return {
    a.xx-b.xx,a.yy-b.yy};
}

double cross(double x1,double y1,double x2,double y2){
    
    return x1*y2-x2*y1;
}

double area(PDD a,PDD b,PDD c){
    
    return cross(b.xx-a.xx,b.yy-a.yy,c.xx-a.xx,c.yy-a.yy);
}

double get(){
    
    sort(A+1,A+idx+1);
    double res=0;
    for(int i=1;i<=idx;i++){
    //idx A little bit 
        while(top>=2 && area(A[stk[top-1]], A[stk[top]], A[i])<=0){
    
            top--;
        }
        stk[++top]=i;
    }
    int k = top;
    for(int i=idx;i>=1;i--){
    
        while(top>k && area(A[stk[top-1]], A[stk[top]], A[i]) <=0){
    
            top--;
        }
        stk[++top]=i;
    }
    for(int i=2;i<=top;i++){
    
        res+=dis(A[stk[i]],A[stk[i-1]]);
    }
    return res;
}

PDD rotate(double x,double y,double theta){
    
    double X = x*cos(theta)+y*sin(theta);
    double Y = -x*sin(theta)+y*cos(theta);
    return {
    X,Y};
}

int main(){
    
    cin>>n>>a>>b>>r;
    /*  The wrong way to write , But it's hard to see bug  The original midpoint (x,y), Don't rotate ,  In the upper right corner is  (x+b/2-r,y+a/2-r)  So the rotation vector is (b/2-r,a/2-r) ,  But you can't use this vector instead of turning four times .  Because the angle between the diagonals of a rectangle is not necessarily 45 degree .  Suppose you rotate  theta  horn , for(int i=1;i<=n;i++){ double x,y,theta; cin>>x>>y>>theta; //  Corner formula ？ PDD t = rotate((b/2)-r,(a/2)-r,-theta); A[++idx] = {x + t.xx, y + t.yy}; A[++idx] = {x + t.xx, y - t.yy}; A[++idx] = {x - t.xx, y + t.yy}; A[++idx] = {x - t.xx, y - t.yy}; } */
    int dx[]={
    1,1,-1,-1};
    int dy[]={
    1,-1,1,-1};
    for(int i=1;i<=n;i++){
    
        double x,y,theta;
        cin>>x>>y>>theta;
        //  Corner formula ？
        for(int j=0;j<4;j++){
    
            //  This bug So hard 
            PDD t = rotate(dx[j]*((b/2)-r),dy[j]*((a/2)-r),-theta);
            A[++idx] = {
    x + t.xx, y + t.yy};    
        }
    }
    // cout<<idx<<endl;
    // for(int i=1;i<=idx;i++){
    
    // cout<<i<<" "<<A[i].xx<<" "<<A[i].yy<<endl;
    // }
    // return 0;
    double res= get();
    printf("%.2lf",res + 2*pi*r);
    return 0;
}