计算几何/平面和空间

总结:

注意几点:

  • 二维向量的操作可以直接使用STL中的complex
  • 问题中的数值往往是浮点数,此时应该注意误差的问题,不考虑误差会WA掉的,这是非常重要的一点

计算几何基础

POJ 1127 Jack Straws

题意:判断给定的若干组线段是否有交点

#include<cstdio>
#include<cmath>
#include<iostream>
using namespace std;
const int maxn = 25;
bool g[maxn][maxn];
double eps = 1e-10;
double add(double a, double b)
{
    if (abs(a + b) < eps * (abs(a) + abs(b)))  return 0;
    else return a + b;
}

struct P {
    double x, y;
    P() {}
    P(double x, double y) : x(x), y(y) {}
    P operator + (P p) {
        return P(add(x, p.x), add(y, p.y));
    }
    P operator -(P p) {
        return P(add(x, -p.x), add(y, -p.y));
    }
    double dot(P p) {
        return add(x * p.x, y *p.y);
    }
    double det(P p) {
        return add(x * p.y, -y * p.x);
    }
    P operator *(double d) {
        return P(x * d, y * d);
    }
};
P p[maxn], q[maxn];
bool onseg(P p1, P p2, P q)
{
    return (p1 - q).det(p2 - q) == 0 && (p1 - q).dot(p2 - q) <= 0;
}
P intersection(P p1, P p2, P q1, P q2)
{
    return p1 + (p2 - p1) * ((q2 - q1).det(q1 - p1) / (q2 - q1).det(p2 - p1));
}
int main(void)
{
    int n;
    while (scanf("%d", &n) && n) {
        for (int i = 0; i < n; i++) {
            scanf("%lf%lf%lf%lf", &p[i].x, &p[i].y, &q[i].x, &q[i].y);
        }
        for (int i = 0; i < n; i++) {
            g[i][i] = true;
            for (int j = 0; j < i; j++) {
                //如果两个直线平行,只需要他们重合且线段有重合部分
                if ((p[i] - q[i]).det(p[j] - q[j]) == 0) {
                    g[i][j] = g[j][i] = onseg(p[i], q[i], p[j]) 
                        || onseg(p[j], q[j], p[i]) 
                        || onseg(p[i], q[i], q[j]) 
                        || onseg(p[j], q[j], q[i]);
                }
                else {
                    //不平行的话先求交点,然后再判断交点是否在两条线段上
                    P tmp = intersection(p[i], q[i], p[j], q[j]);
                    g[i][j] = g[j][i] = onseg(p[i], q[i], tmp) && onseg(p[j], q[j], tmp);
                }
            }
        }
        //最后用Floyd-Warshall算法求任意两个棍子是否通过其他的棍子相连
        for (int k = 0; k < n; k++) {
            for (int i = 0; i < n; i++) {
                for (int j = 0; j < n; j++)
                    g[i][j] |= g[i][k] && g[k][j];
            }
        }
        int a, b;
        while (scanf("%d%d", &a, &b) && a + b) {
            if (g[a - 1][b - 1])
                printf("CONNECTED
");
            else
                printf("NOT CONNECTED
");
        }
    }
}

极限情况

AOJ 2308

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;

const double g = 9.8;   //重力加速度
const double EPS = 1e-10;
int N;
double V, X, Y;
double L[55], B[55], R[55], T[55];

// 计算以 vy 的速度竖直向上射出 t 秒后的高度
double calc(double vy,double t)
{
    return vy*t-g*t*t/2;
}
// a 相对 lb 和 ub 的位置
int cmp(double lb, double ub, double a)
{ 
   return a < lb + EPS ? -1 : a > ub - EPS ? 1 : 0;
}
// 判断当射出路径经过 (qx, qy) 时,蛋能否击中猪
bool check(double qx, double qy)
{
    // 设初速度在 x 和 y 方向上的分量分别为 vx 和 vy,过(qx,qy)的时间为 t
    // 求解联立方程 vx^2 + vy^2 = v^2, vx * t = qx, vy * t - 1/2 * g * t^2 = qy
    double a = g * g / 4, b = g * qy - V * V, c = qx * qx + qy * qy;
    double D = b * b - 4 * a * c;
    if(D < 0 && D > -EPS) D = 0;
    if(D < 0) return false;
    for(int d = -1; d <= 1; d += 2)  // 验证联立方程的两个解
    {
        double t2 = (-b + d * sqrt(D)) / (2 * a);
        if(t2 <= 0) continue;
        double t = sqrt(t2);
        double vx = qx / t, vy = (qy + g * t * t / 2) / t;
        
        // 判断是否通过猪的正上方
        double yt = calc(vy, X / vx);
        if(yt < Y - EPS) continue;
        
        bool ok = true;
        for(int i = 0; i < N; i++)
        {
            if(L[i] >= X) continue;
            // 判断在猪正上方的鸟和猪之间是否有障碍物
            if(R[i] == X && Y <= T[i] && B[i] <= yt)
                ok = false;
            //判断在飞到猪正上方之前是否会撞到障碍物
            int yL = cmp(B[i], T[i], calc(vy, L[i] / vx));  //左侧的相对位置
            int yR = cmp(B[i], T[i], calc(vy, R[i] / vx));  //右侧的相对位置
            int xH = cmp(L[i], R[i], vx * (vy / g));   //最高点的相对位置
            int yH = cmp(B[i], T[i], calc(vy, vy / g));
            if(xH == 0 && yH >= 0 && yL < 0)
                ok = false;
            if(yL * yR  <= 0)
                ok = false;
        }
        if(ok)
            return true;
    }
    return false;
}
void solve(void)
{
    // 截掉猪以右的障碍物
    for(int i = 0; i < N; i++)
        R[i] = min(R[i], X);
    bool ok = check(X, Y);   // 直接撞上猪
    for(int i = 0; i < N; i++)
    {
        ok |= check(L[i], T[i]);
        ok |= check(R[i], T[i]);
    }
    puts(ok ? "Yes" : "No");
}
int main()
{
    scanf("%d %lf %lf %lf", &N, &V, &X, &Y);
    for(int i = 0; i < N; i++)
        scanf("%lf %lf %lf %lf", &L[i], &B[i], &R[i], &T[i]);
    solve();
    return 0;
}

POJ 1981 Circle and Points

题意:平面上有N个点,用单位圆去套,最多能套几个

思路:当然是先想到枚举两个,再跑一边所有点的O(n^3)的优秀算法,// 滑稽一下,不过这数据看起来也不一定会T掉,不管辣

point getmid(point p1, point p2)
{
    point mid,center;
    mid.x = (p1.x + p2.x) / 2.0;
    mid.y = (p1.y + p2.y) / 2.0;
    double angle = atan2(p1.x - p2.x, p2.y - p1.y); // 从平行四边形的角度去看,求经过中点的那条边的角度
    double dcm = sqrt(1-dis(p1, mid) * dis(p1, mid)); // 因为是单位圆,我们要计算出以二者中心为圆心,距离一半为半径的圆的半径和单位圆半径之差
    center.x = mid.x + dcm * cos(angle);//三角形思路。。。
    center.y = mid.y + dcm * sin(angle);
    return center;
}
计算两点之间圆心的函数
// #include<bits/stdc++.h>
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;

const int INF    = 0x3f3f3f3f;
const int MAX_N  = 500 + 50;
const int MOD    = 1000000007;
const double PI  = acos(-1.0);
const double EPS = 1e-10;

struct point
{
    double x,y;
    point(double x=0,double y=0):x(x),y(y) {}
}p[310];

struct node
{
    double ang;
    int in;
}arc[310*300];

double dis(point a,point b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}

int dcmp(double x)
{
    if(fabs(x)<EPS) return 0;
    else return x < 0 ? -1 : 1 ;
}

bool cmp(node a,node b)
{
    if(dcmp(a.ang-b.ang)==0) return a.in>b.in;
    return dcmp(a.ang-b.ang)<0;
}

int main(void)
{
    int i,j,n;
    while(scanf("%d",&n)&&n)
    {
        for(i = 1 ; i <= n; i++)
            scanf("%lf%lf",&p[i].x,&p[i].y);
        int g = 0;
        int ans = 0,maxz = 1;
        for(i = 1; i <= n ; i++)
        {
            ans = 0; g = 0;
            for(j = 1; j <= n ; j++)
            {
                if(dis(p[i],p[j])>2.0) continue;
                double ang1 = atan2(p[j].y-p[i].y,p[j].x-p[i].x);
                double ang2 = acos(dis(p[i],p[j])/2);
                arc[++g].ang = ang1-ang2;//这里角度的算法很巧妙
                arc[g].in = 1;
                arc[++g].ang = ang1+ang2;
                arc[g].in = -1;
            }
            sort(arc+1,arc+g+1,cmp);
            
            for(j = 1 ; j <= g;j++)
            {
                ans+=arc[j].in;
                maxz = max(maxz,ans);
            }
        }
        printf("%d
",maxz);
    }
    return 0;
}

POJ 1418

题意:给定Confetti的尺寸和位置以及它们的叠放次序,计算出有多少Confetti是可以看见的

#include <iostream>
#include <cmath>
#include <algorithm>
#include <cstdio>
using namespace std;

const int  MAX_N = 128;
const double EPS = 5e-13;
const double PI = acos(-1.0);

typedef struct
{
    double x, y;
} point;

double Distance(const point & p1, const point & p2){return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));}

double MainAngle(double a)
{
    while (a > 2 * PI) a -= 2 * PI;
    while (a < 0) a += 2 * PI;
    return a;
}

int  n;
point o[MAX_N];  //圆心
double r[MAX_N];  //圆的弧度

int  pan;   //与这个圆的交点数目
double pa[2 * MAX_N]; //存放与这个圆所有交点对应的弧度
int  visible[MAX_N];
int  ans;

int main()
{
    int i, j, k, t;
    point tp;
    double a, b, d;
    while (scanf("%d", &n)&&n)
    {
        for (i = 0; i < n; ++i)
        {
            scanf("%lf %lf %lf", &o[i].x, &o[i].y, &r[i]);
            visible[i] = 0;
        }

        for (i = 0; i < n; ++i)
        {
            pan = 0;
            pa[pan++] = 0;
            pa[pan++] = 2 * PI;

            for (j = 0; j < n; ++j)
            {
                if (j == i) continue;

                d = Distance(o[i], o[j]); //判断两个圆心距离

                if (r[i] + r[j] < d || r[i] + d < r[j] || r[j] + d < r[i]) //包含或不相交的
                    continue;

                a = atan2(o[j].y - o[i].y, o[j].x - o[i].x);//atan2(),是求这个点和x轴正方形夹角,*PI/180  得到度数
                b = acos((r[i] * r[i] + d * d - r[j] * r[j]) / (2 * r[i] * d));

                pa[pan] = MainAngle(a + b);
                pan++;
                pa[pan] = MainAngle(a - b);
                pan++;
            }

            sort(pa, pa + pan);

            for (j = 0; j < pan - 1; ++j)
            {
                a = (pa[j] + pa[j + 1]) / 2;

                for (t = -1; t <= 1; t += 2) //t = -1 或 1
                {
                    //将每段圆弧中点往内外各移动很小距离
                    tp.x = o[i].x + (r[i] + t * EPS) * cos(a);
                    tp.y = o[i].y + (r[i] + t * EPS) * sin(a);

                    for (k = n - 1; k >= 0; --k)
                        //如果找到第一个cover point i 的 Arc j 的圆,break
                        if (Distance(tp, o[k]) < r[k])
                            break;
                    visible[k] = 1;
                }
            }
        }
        ans = 0;
        for (i = 0; i < n; ++i)
            if (visible[i] == 1)
                ans++;

        printf("%d
", ans);
    }

    return 0;
}

AOJ 2201 Immortal Jewels

 题意:有n个宝石,宝石为圆形。给出现在有根金属棒,靠近宝石距离m之内就会吸附上去。现在给出每颗宝石的信息,求用这个棒子一次能钓出最多多少个宝石?

看起来简单但是感觉着实挺复杂,先空着吧,之后再来

平面扫描

POJ 2932 Coneology

 题意:平面上有N个两两不相交的圆,求全部最外层的,即不被其它圆包括的圆的个数并输出

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

const int INF    = 0x3f3f3f3f;
const int MAX_N  = 8e5 + 50;
const int MOD    = 1000000007;
const double PI  = acos(-1.0);
const double EPS = 1e-10;

double x[MAX_N],y[MAX_N],r[MAX_N];
int n;

bool inside(int i,int j)
{
    double dx=x[i]-x[j];
    double dy=y[i]-y[j];
    return dx*dx + dy*dy <=r[j]*r[j];
}

void solve()
{
    vector<pair<double,int> > e; // 圆的左右坐标
    e.clear();
    for(int i=0;i<n;i++)
    {
        e.push_back(make_pair(x[i]-r[i],i));
        e.push_back(make_pair(x[i]+r[i],i));
    }
    sort(e.begin(),e.end());

    set<pair<double,int> > out;
    vector<int> res;
    res.clear(); out.clear();

    for(int i=0;i<e.size();i++)
    {
        int id=e[i].second%n;
        if(e[i].second<n)
        {
            set<pair<double,int> >::iterator it=out.lower_bound(make_pair(y[id],id));
            if(it != out.end() && inside (id,it->second)) continue;
            if(it != out.begin() && inside (id ,(--it)->second)) continue;
            res.push_back(id);
            out.insert(make_pair(y[id],id));
        }else{
            out.erase(make_pair(y[id],id));
        }
    }
    
    sort(res.begin(),res.end());
    printf("%d
",res.size());
    for(int i=0;i<res.size();i++)
       printf("%d%c",res[i]+1,i+1==res.size() ? '
':' ');
}

int main(void)
{
    while(scanf("%d",&n)==1)
    {
        for(int i=0;i<n;i++)
            scanf("%lf%lf%lf",&r[i],&x[i],&y[i]);
        solve();
    }

    return 0;
}

POJ Barn Expansion

题意:有N块不重叠的矩形地,由左下角(A,B)和右上角(C,D)决定。如果两块地的边或角相交,则两块地都无法扩大。求多少地可以扩大?

#include<bits/stdc++.h>

const int INF    = 0x3f3f3f3f;
const int MAX_N  = 8e5 + 50;
const int MOD    = 1000000007;
const double PI  = acos(-1.0);
const double EPS = 1e-10;

struct data_
{
    int l,r,now,id;
    data_(){}
    data_(int a,int b,int c,int d)
    {
        l=a,r=b,now=c,id=d;
    }
};

bool cmp(data_ a,data_ b)
{
    if(a.now!=b.now) return a.now<b.now;
    else if(a.l!=b.l) return a.l<b.l;
    else return a.r<b.r;
}

vector<data_> sx,sy;
int n;
int vis[25005];

int main(void)
{
    scanf("%d",&n);
    for(int i=0;i<n;i++)
    {
        int a,b,c,d;
        scanf("%d%d%d%d",&a,&b,&c,&d);
        sx.push_back(data_(b,d,c,i));
        sx.push_back(data_(b,d,a,i));
        sy.push_back(data_(a,c,b,i));
        sy.push_back(data_(a,c,d,i));
    }
    sort(sx.begin(),sx.end(),cmp);
    sort(sy.begin(),sx.end(),cmp);

    int up=sx[0].r;
    for(int i=1;i<sx.size();i++)
    {
        if(sx[i].now == sx[i-1].now)
            if(up>=sx[i].l)
                vis[sx[i].id]=vis[sx[i-1].id]=1;
        else up=sx[i].r;

        up=max(up,sx[i].r);
    }

    up=sy[0].r;
    for(int i=1;i<sy.size();i++)
    {
        if(sy[i].now == sy[i-1].now)
            if(up>=sy[i].l)
                vis[sy[i].id]=vis[sy[i-1].id]=1;

        else up=sy[i].r;

        up=max(up,sy[i].r);
    }
    int ans=0;
    for(int i=1;i<=n;i++)
        if(!vis[i]) ans++;
    printf("%d
",ans);

    return 0;
}

POJ 3293 Rectilinear polygon

题意:给定N个点,问是否能组成直角多边形(每个顶点都与另外两个顶点构成直角,每条边都平行于坐标轴),并求出周长

思路:扫面线移动时,只有点数为偶数的情况,然后有一点就是要判断横竖线段是否相交

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

const int INF    = 0x3f3f3f3f;
const int MAX_N  = 1e5 + 50;
const int MOD    = 1000000007;
const double PI  = acos(-1.0);
const double EPS = 1e-10;

struct point
{
    int x,y;
    int id;
}p[MAX_N];

struct line
{
    int x,y1,y2;
    line(){}
    line(int a,int b,int c){x=a,y1=b,y2=c;}
}L[MAX_N];

int con[MAX_N][2];
int n,ln;

bool cmp_x(const point a,const point b)
{
    if(a.x!=b.x) return a.x<b.x;
    else return a.y<b.y;
}

bool cmp_y(const point a,const point b)
{
    if(a.y!=b.y) return a.y<b.y;
    else return a.x<b.x;
}

bool judge(const point a,const point b)
{
    int y=a.y,x1=a.x,x2=b.x;
    for(int i=0;i<ln;i++)
    {
        if(x1<L[i].x&&x2>L[i].x&&L[i].y1<y&&L[i].y2>y)
            return 1;
    }
    return 0;
}

int main(void)
{
    int T;
    scanf("%d",&T);

    while(T--)
    {
        scanf("%d",&n);
        for(int i=0;i<n;i++)
        {
            scanf("%d%d",&p[i].x,&p[i].y); p[i].id=i;
        }
        
        sort(p,p+n,cmp_x);
        int sum=0,cnt=1,flage=0;
        ln=0;
        for(int i=1;i<n&&!flage;i++)
        {
            if(p[i].x!=p[i-1].x){
                if(cnt&1) flage=1;
                cnt=1;
            }else{
                cnt++;
                if((cnt&1)==0)
                {
                    sum += p[i].y-p[i-1].y;
                    con[p[i].id][0]=p[i-1].id;
                    con[p[i-1].id][0]=p[i].id;
                    L[ln++]=line(p[i].x,p[i-1].y,p[i].y);
                }
            }
        }
        if(flage==1) { printf("-1
");continue;}
        
        cnt=1;
        sort(p,p+n,cmp_y);
        for(int i=1;i<n&&!flage;i++)
        {
            if(p[i].y!=p[i-1].y)
            {
                if(cnt&1)
                    flage=1;
                cnt=1;
            }else{
               cnt++;
               if((cnt&1)==0)
                {
                    sum += p[i].x-p[i-1].x;
                    con[p[i].id][1]=p[i-1].id;
                    con[p[i-1].id][1]=p[i].id;
                    if(judge(p[i-1],p[i]))
                        flage=1;
                }
            }
        }
        if(flage==1){ printf("-1
");continue;}
        int s=1,x=0,ans=0;
        do{
            x=con[x][s];
            s=1^s;
            ans++;
        }while(x!=0 && !flage);
        if(ans!=n)  printf("-1
");
        else printf("%d
",sum);
    }

    return 0;
}

POJ 2482 Stars in Your Window

题意:夜空有n个星星,坐标(x,y)亮度c。用长W宽H的窗户去套,问能套住的星星的亮度之和的最大值?

思路:用RMQ的思路

#include<bits/stdc++.h>

using namespace std;
typedef long long ll;

const int INF    = 0x3f3f3f3f;
const int MAX_N  = 1e4 + 50;
const int MOD    = 1000000007;
const double PI  = acos(-1.0);
const double EPS = 1e-10;

ll xs[MAX_N], ys[MAX_N];
int cs[MAX_N];
ll X[MAX_N * 2], Y[MAX_N * 2];    
pair<pair<int,int>, pair<int,int> > event[MAX_N*2];

int node_value[MAX_N*8],node_max[MAX_N*8];

void add(int from,int to,int v,int i,int l,int r)
{
    if(from<=l && r<=to)
    {
        node_value[i] += v;
        node_max[i] +=v;
        return;
    }
    int m=(l+r)>>1;
    if(from<=m) add(from,to,v,i*2,l,m);
    if(m<to) add(from,to,v,i*2+1,m+1,r);
    node_max[i]=max(node_max[i*2],node_max[i*2+1])+node_value[i];
}

int main(void)
{
    int n,W,H;
    while(~scanf("%d%d%d",&n,&W,&H))
    {
        for(int i=0;i<n;i++)
        {
            scanf("%lld%lld%d",xs+i,ys+i,cs+i);
            xs[i]*=2, ys[i]*=2;
        }
        for(int i=0;i<n;i++)
        {
            X[i*2] = xs[i]-W;
            X[i*2+1] = xs[i]+W;
            Y[i*2] = ys[i]-H;
            Y[i*2+1] = ys[i]-1+H;
        }
        sort(X,X + n*2);
        sort(Y,Y + n*2);
        for(int i=0;i<n;i++)
        {
            event[i*2] = make_pair(make_pair(lower_bound(X,X+n*2,xs[i]-W)-X,cs[i]), make_pair(lower_bound(Y,Y+n*2,ys[i]-H)-Y,lower_bound(Y,Y+n*2,ys[i]+H-1)-Y));
            event[i*2+1] = make_pair(make_pair(lower_bound(X,X+n*2,xs[i]+W)-X,-cs[i]), make_pair(lower_bound(Y,Y+n*2,ys[i]-H)-Y,lower_bound(Y,Y+n*2,ys[i]+H-1)-Y));
        }
        sort(event, event + n * 2);
        int ans = 0;
        for(int i=0;i<n*2;i++)
        {
            add(event[i].second.first,event[i].second.second,event[i].first.second,1,0,n*2);
            ans=max(ans,node_max[1]);
        }
        printf("%d
",ans);
    }
    
    return 0;
}

凸包

POJ 2187 Beauty Contest

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

const int INF    = 0x3f3f3f3f;
const int MAX_N  = 5e4 + 50;
const int MOD    = 1000000007;
const double PI  = acos(-1.0);
const double EPS = 1e-10;

double add(double a,double b)
{
    if(abs(a+b)<EPS*(abs(a+b))) return 0;
    return a+b;
}

struct P
{
    double x,y;
    P(){}
    P(double x,double y):x(x),y(y){}

    P operator+(P p){return P(add(x,p.x),add(y,p.y));}
    P operator-(P p){return P(add(x,-p.x),add(y,-p.y));}
    P operator*(double d){return P(x*d,y*d);}
    double dot(P p){return add(x*p.x,y*p.y);}
    double det(P p){return add(x*p.y,-y*p.x);}
};

bool cmp_x(const P& a,const P& b)
{
    if(a.x!=b.x) return a.x<b.x;
    return a.y<b.y;
}

vector<P> convex_hull(P *ps,int n)
{
    sort(ps,ps+n,cmp_x);
    int k=0;
    vector<P>qs(n*2);
    for(int i=0;i<n;i++)
    {
        while(k>1 && (qs[k-1]-qs[k-2]).det(ps[i]-qs[k-1])<=0)
            k--;
        qs[k++]=ps[i];
    }
    for(int i=n-2,t=k;i>=0;i--)
    {
        while(k>t&&(qs[k-1]-qs[k-2]).det(ps[i]-qs[k-1])<=0)
            k--;
        qs[k++]=ps[i];
    }
    qs.resize(k-1);
    return qs;
}

double dist(P a,P b){return (a-b).dot(a-b);}

int N;
P ps[MAX_N];

void solve()
{
    vector<P>qs = convex_hull(ps,N);
    double res=0;
    for(int i=0;i<qs.size();i++)
        for(int j=0;j<i;j++)
            res=max(res,dist(qs[i],qs[j]));
    printf("%.0f
",res);
}

int main(void)
{
    scanf("%d",&N);
    for(int i=0;i<N;i++)
        scanf("%lf%lf",&ps[i].x,&ps[i].y);
    solve();
    
    return 0;
}

POJ 1113 Wall

题意:皇帝要你造墙将城堡围起来,城堡的顶点有N个,墙必须离城堡的边至少L单位远,并且墙的总长度尽量小。求此长度?

思路:将一个“凹”多边形转化为一个“凸”多边形后求周长

#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
using namespace std;
 
#define MAX_N 1000 + 16
#define M_PI 3.14159265358979323846

struct P
{
    int x, y;
    P() {}
    P(int x, int y) : x(x), y(y) {}
    P operator + (P p){ return P(x + p.x, y + p.y); }
    P operator - (P p){ return P(x - p.x, y - p.y); }
    P operator * (int d){ return P(x*d, y*d); }
    bool operator < (const P& a) const
    {
        if (x != a.x) return x < a.x;
        else return y < a.y;
    }
    int dot(P p) { return x*p.x + y*p.y; }
    int det(P p) { return x*p.y - y*p.x; }
};
 
// 字典序比较
bool cmp_x(P a, P b)
{
    if (a.x != b.x) return a.x < b.x;
    return a.y < b.y;
}
// 求凸包
vector<P> convex_hull(P *ps, int n)
{
    sort(ps, ps + n, cmp_x);
    int k = 0;    // 凸包的顶点数
    vector<P> qs(n * 2);    // 构造中的凸包
    // 构造凸包的下侧
    for (int i = 0; i < n; ++i)
    {
        while (k > 1 && (qs[k - 1] - qs[k - 2]).det(ps[i] - qs[k - 1]) <= 0) --k;
        qs[k++] = ps[i];
    }
    // 构造凸包的上侧
    for (int i = n - 2, t = k; i >= 0; --i)
    {
        while (k > t && (qs[k - 1] - qs[k - 2]).det(ps[i] - qs[k - 1]) <= 0) --k;
        qs[k++] = ps[i];
    }
    qs.resize(k - 1);
    return qs;
}
 
// 距离的平方
double dist(P p, P q)
{
    return sqrt((double)(p - q).dot(p - q));
}
 
// 求解凸包对踵点最大距离
double max_distance(P *ps, int N)
{
    vector<P> qs = convex_hull(ps, N);
    int n = qs.size();
    if (n == 2)    return dist(qs[0], qs[1]);    // 特别处理凸包退化的情况
    int i = 0, j = 0;    // 某个方向上的对踵点对
    // 求出x轴方向上的对踵点对
    for (int k = 0; k < n; k++)
    {
        if (!cmp_x(qs[i], qs[k])) i = k;
        if (cmp_x(qs[j], qs[k])) j = k;
    }
    double res = 0;
    int si = i, sj = j;
    while (i != sj || j != si)    // 将方向逐步旋转180度
    {
        res = max(res, dist(qs[i], qs[j]));
        // 判断先转到边i-(i+1)的法线方向还是边j-(j+1)的法线方向
        if ((qs[(i + 1) % n] - qs[i]).det(qs[(j + 1) % n] - qs[j]) < 0)
            i = (i + 1) % n;    //  先转到边i-(i+1)的法线方向
        else
            j = (j + 1) % n;    // 先转到边j-(j+1)的法线方向
    }
    return res;
}
 
// 求解凸包周长
double total_distance(P *ps, int N)
{
    vector<P> qs = convex_hull(ps, N);
    int n = qs.size();
    double res = 0;
    for (int i = 0; i < n; ++i)
        res += dist(qs[i], qs[(i + 1) % n]);
 
    return res;
}
 
P ps[MAX_N];
 
int main(int argc, char *argv[])
{
    int N, L;
    while (~scanf("%d%d", &N, &L))
    {
        for (int i = 0; i < N; ++i)
            scanf("%d%d", &ps[i].x, &ps[i].y);
        double ans = (2 * M_PI * L + total_distance(ps, N));
        printf("%.0lf
", ans);
    }

    return 0;
}

POJ 1912 A highway and the seven dwarfs 

题意:侏儒岛上有N栋房子,组成一个社区。现给定一条高铁,问该高铁是否分割了社区?

即是问n个点的凸包是否在直线同一侧。

#include<bits/stdc++.h>

using namespace std;
typedef long long ll;

const int INF    = 0x3f3f3f3f;
const int MAX_N  = 1e5 + 50;
const int MOD    = 1000000007;
const double PI  = acos(-1.0);
const double EPS = 1e-10;

double add(double a,double b)
{
    if(abs(a+b)<EPS*(abs(a+b))) return 0;
    return a+b;
}

struct P
{
    double x,y;
    P(){}
    P(double x,double y):x(x),y(y){}

    P operator + (P p){return P(add(x,p.x),add(y,p.y));}
    P operator - (P p){return P(add(x,-p.x),add(y,-p.y));}
    P operator * (double d){return P(x*d,y*d);}   
    bool operator < (const P& a) const
    {
        if(x != a.x) return x<a.x;
        else return y<a.y;
    }
    double dot(P p){return add(x*p.x,y*p.y);}
    double det(P p){return add(x*p.y,-y*p.x);}
};

vector<P> convex_hull(P *ps,int n)
{
    sort(ps,ps+n);
    int k=0;
    vector<P> qs(n*2);
    for(int i=0;i<n;i++)
    {
        while(k>1 && (qs[k-1]-qs[k-2]).det(ps[i]-qs[k-1])<=0) k--;
        qs[k++]=ps[i];
    }
    for(int i=n-2,t=k;i>=0;i--)
    {
        while(k>t&&(qs[k-1]-qs[k-2]).det(ps[i]-qs[k-1])<=0) k--;
        qs[k++]=ps[i];
    }
    qs.resize(k-1);
    return qs;
}

P hs[MAX_N];
double as[MAX_N];

inline double normalize(double r) {if(r < -PI/2.0 + EPS) return r+=PI*2; return r;}

inline double atan2(const P&p) {return normalize(atan2(p.y,p.x));}

inline bool double_cmp(double a,double b){return a+EPS<b;}

int main(void)
{
    int N;
    scanf("%d",&N);
    for(int i=0;i<N;i++)
        scanf("%lf%lf",&hs[i].x,&hs[i].y);
    vector<P> chs;
    int n=0;
    if(N>1)
    {
        chs = convex_hull(hs,N);
        n=chs.size();
        chs.push_back(chs[0]);
    }
    for (int i = 0; i < n; ++i) as[i] = atan2(chs[i + 1] - chs[i]);

    sort(as,as+n,double_cmp);
    P p1,p2;
    while(~scanf("%lf%lf%lf%lf",&p1.x, &p1.y, &p2.x, &p2.y))
    {
        if(N<=1){puts("GOOD"); continue;}
        
        int i = upper_bound(as, as + n, atan2(p2 - p1), double_cmp) - as;
        int j = upper_bound(as, as + n, atan2(p1 - p2), double_cmp) - as;
        puts((((p2 - p1).det(chs[i] - p1) * (p2 - p1).det(chs[j] - p1) > -EPS)) ? "GOOD" : "BAD");
    }

    return 0;
}

POJ 3608 Bridge Across Islands

题意:在两个凸包小岛之间造桥,求最小距离?

#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
using namespace std;
 
#define MAX_N 10000 + 16
#define INF 0x3F3F3F3F
#define EPS 1E-10
 
struct Point
{
    double x, y;
    Point() {}
    Point(double x, double y) : x(x), y(y) {}
    Point operator + (const Point& p){ return Point(x + p.x, y + p.y); }
    Point operator - (const Point& p){ return Point(x - p.x, y - p.y); }
    Point operator * (const double& d){ return Point(x * d, y * d); }
    bool operator < (const Point& a) const
    {
        if (x != a.x) return x < a.x;
        else return y < a.y;
    }
    double dot(const Point& p) { return x * p.x + y * p.y; }
    double det(const Point& p) { return x * p.y - y * p.x; }
};
 
Point P[MAX_N], Q[MAX_N];

inline double cross(Point A, Point B, Point C){ return (B - A).det(C - A);}
 
inline double multi(Point A, Point B, Point C){ return (B - A).dot(C - A);}

inline double dist(Point A, Point B){ return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));}

inline void anticlockwise(Point* p, int N)
{
    for (int i = 0; i < N - 2; ++i)
    {
        double tmp = cross(p[i], p[i + 1], p[i + 2]);
        if (tmp > EPS) return;
        else if (tmp < -EPS){ reverse(p, p + N); return;}
    }
}
 
//计算C点到线段AB的最短距离
inline double point_to_line(Point A, Point B, Point C)
{
    if (dist(A, B) < EPS) return dist(B, C);
    if (multi(A, B, C) < -EPS) return dist(A, C);
    if (multi(B, A, C) < -EPS) return dist(B, C);
    return fabs(cross(A, B, C) / dist(A, B));
}
 
inline double line_to_line(Point A, Point B, Point C, Point D)
{
    return min(min(point_to_line(A, B, C), point_to_line(A, B, D)), min(point_to_line(C, D, A), point_to_line(C, D, B)));
}
 
double solve(Point* P, Point* Q, int n, int m)
{
    int yminP = 0, ymaxQ = 0;
    for (int i = 0; i < n; ++i) if (P[i].y < P[yminP].y) yminP = i;
    for (int i = 0; i < m; ++i) if (Q[i].y > Q[ymaxQ].y) ymaxQ = i;
    P[n] = P[0];Q[m] = Q[0];
    double arg, ans = INF;
    for (int i = 0; i < n; ++i)
    {
        while (arg = cross(P[yminP + 1], Q[ymaxQ + 1], P[yminP]) - cross(P[yminP + 1], Q[ymaxQ], P[yminP]) > EPS) ymaxQ = (ymaxQ + 1) % m;
        ans = min(ans, line_to_line(P[yminP], P[yminP + 1], Q[ymaxQ], Q[ymaxQ + 1]));
        yminP = (yminP + 1) % n;
    }
    return ans;
} 

int main(void)
{

    int N, M;
    while (~scanf("%d%d", &N, &M) && N + M)
    {
        for (int i = 0; i < N; ++i)
            scanf("%lf%lf", &P[i].x, &P[i].y);
        for (int i = 0; i < M; ++i)
            scanf("%lf%lf", &Q[i].x, &Q[i].y);
            
        anticlockwise(P, N);
        anticlockwise(Q, M);
        printf("%.5lf
", solve(P, Q, N, M));
    }

    return 0;
}

POJ 2079 Triangle

求N个点组成的三角形的最大面积

#include<bits/stdc++.h>

using namespace std;
typedef long long ll;

const int INF    = 0x3f3f3f3f;
const int MAX_N  = 5e4 + 50;
const int MOD    = 1000000007;
const double PI  = acos(-1.0);
const double EPS = 1e-10;

struct Point
{
    int x,y;
    Point() {}
    Point(int x, int y) : x(x), y(y) {}
    Point operator + (const Point& p){ return Point(x + p.x, y + p.y); }
    Point operator - (const Point& p){ return Point(x - p.x, y - p.y); }
    Point operator * (const int& d){ return Point(x * d, y * d); }
    bool operator < (const Point& a) const
    {
        if (x != a.x) return x < a.x;
        else return y < a.y;
    }
    int dot(const Point& p) { return x * p.x + y * p.y; }
    int det(const Point& p) { return x * p.y - y * p.x; }
}P[MAX_N];

inline int cross(Point A, Point B, Point C){return (B - A).det(C - A);}

inline int compute_area(Point A, Point B, Point C){    return abs(cross(A, B, C));}

vector <Point> convex_hull(Point *ps, int N)
{
    sort(ps,ps+N);
    int k=0;
    vector<Point> qs(N*2);
    for(int i=0;i<N;i++)
    {
        while(k>1 && (qs[k-1]-qs[k-2]).det(ps[i]-qs[k-1])<=0)
            --k;
        qs[k++]=ps[i];
    }
    for(int i=N-2, t=k ;i>=0;i--)
    {
        while(k > t && (qs[k - 1] - qs[k - 2]).det(ps[i] - qs[k - 1]) <= 0)
            k--;
        qs[k++]=ps[i];
    }
    qs.resize(k-1);
    return qs;
}

int main(void)
{
    int N;
    while(~scanf("%d",&N) && N>0)
    {
        for(int i=0;i<N;i++)
            scanf("%d%d",&P[i].x,&P[i].y);
        
        vector<Point> ps=convex_hull(P,N); N=ps.size();
        int ans=0;

        for(int offset=1; offset<(N+1)/2;offset++)
        {
            int first=(offset+1)%N;
            for(int third=0;third<N;third++)
            {
                int second=(third+offset)%N;
                int prev=compute_area(ps[third],ps[second],ps[first]);
                for(++first;first!=second && first!=third;first++)
                {
                    first=(first+N)%N;                    
                    int cur = compute_area(ps[third], ps[second], ps[first]);
                    ans = max(ans, prev);
                    if (cur <= prev) break;
                    prev = cur;
                }
                --first; first=(first+N)%N;
            }
        }
        printf("%d.%s
", ans / 2, ans % 2 == 1 ? "50" : "00");
    }

    return 0;
}

POJ 3246 Game

题意:N个点中去掉一个得到N个点集,求这些点集构成的凸包的最小面积

#include<bits/stdc++.h>

using namespace std;
typedef long long ll;

const int INF    = 0x3f3f3f3f;
const int MAX_N  = 1e5 + 50;
const int MOD    = 1000000007;
const double PI  = acos(-1.0);
const double EPS = 1e-10;

struct Point
{
    int id;
    int x,y;
    Point() {}
    Point(int x, int y) : x(x), y(y) {}
    Point operator + (const Point& p){ return Point(x + p.x, y + p.y); }
    Point operator - (const Point& p){ return Point(x - p.x, y - p.y); }
    Point operator * (const int& d){ return Point(x * d, y * d); }
    bool operator < (const Point& a) const
    {
        if (x != a.x) return x < a.x;
        else return y < a.y;
    }
    int dot(const Point& p) { return x * p.x + y * p.y; }
    int det(const Point& p) { return x * p.y - y * p.x; }
}P[MAX_N],Q[MAX_N];

inline int cross(Point A, Point B, Point C){return (B - A).det(C - A);}

inline int compute_area(Point A, Point B, Point C){    return abs(cross(A, B, C));}

inline int compute_area(const vector<Point>& ps)
{
    int total=0;
    for(int i=2;i<ps.size();i++)
        total += compute_area(ps[0],ps[i-1],ps[i]);
    return total;
}

vector <Point> convex_hull(Point *ps, int N)
{
    sort(ps,ps+N);
    int k=0;
    vector<Point> qs(N*2);
    for(int i=0;i<N;i++)
    {
        while(k>1 && (qs[k-1]-qs[k-2]).det(ps[i]-qs[k-1])<=0) --k;
        qs[k++]=ps[i];
    }
    for(int i=N-2, t=k ;i>=0;i--)
    {
        while(k > t && (qs[k - 1] - qs[k - 2]).det(ps[i] - qs[k - 1]) <= 0)  k--;
        qs[k++]=ps[i];
    }
    qs.resize(k-1);
    return qs;
}

int main(void)
{
    int N;
    while(~scanf("%d",&N) && N>0)
    {
        for(int i=0;i<N;i++)
        {
            scanf("%d%d",&P[i].x,&P[i].y);
            P[i].id=i;
        }
        memcpy(Q,P,N*sizeof(Point));
        vector<Point> ps = convex_hull(P,N);
        int ans=INF;
        for(int i=0;i<ps.size();i++)
        {
            memcpy(P,Q,N*sizeof(Point));
            swap(P[ps[i].id],P[N-1]);
            ans=min(ans,compute_area(convex_hull(P,N-1)));
        }
        printf("%d.%s
",ans/2,ans%2 ==1 ? "50":"00");
    }

    return 0;
}

POJ 3689 Equations

数值积分

 总结:

解决问题:求凸多面体的体积

重点:Simpson公式,本质上就是用二次曲线来逼近函数,精度比矩形切割要好。

我们只需把两个多边形的顶点的X全部存起来,从小到大排序,那么挨着的两个X代表的区间绝对是连续的,可以看出P1的底面在x处的截距y是关于x的线性函数,

所以他们的乘积f(x)=(z2z1)*(y2y1)是个二次函数,

所以有公式:

#include<bits/stdc++.h>

using namespace std;
typedef long long ll;

const int INF    = 0x3f3f3f3f;

int M,N;
int X1[110],Y1[110],X2[110],Z2[110];

double get_width(int* X, int* Y, int n, double x)
{
    double lb=INF, ub=-INF;
    for(int i=0;i<n;i++)
    {
        double x1=X[i], y1=Y[i], x2 =X[(i+1)%n], y2=Y[(i+1)%n];
        if((x1-x)*(x2-x)<=0 && x1 != x2)
        {
            double y=y1+(x-x1)*(y2*1.0-y1)/(x2-x1);
            lb=min(lb,y);
            ub=max(ub,y);
        }
    }
    return max(0.0,ub-lb);
}

double f(double x){ return get_width(X1,Y1,M,x)*get_width(X2,Z2,N,x);}

double simpson(double a,double b)
{
    double c=(a+b)/2;
    return (b-a)/6.0*(f(a)+4*f(c)+f(b));
}

vector<int> sx;

int main(void)
{
    while(~scanf("%d%d",&M,&N))
    {
        if(M==0 && N==0) return 0;
        sx.clear();
        for(int i=0;i<M;i++)
        {
            scanf("%d%d",&X1[i],&Y1[i]);
            sx.push_back(X1[i]);
        }
        for(int i=0;i<N;i++)
        {
            scanf("%d%d",&X2[i],&Z2[i]);
            sx.push_back(X2[i]);
        }
        sort(sx.begin(),sx.end());
        int min1=*min_element(X1,X1+M),max1=*max_element(X1,X1+M);
        int min2=*min_element(X2,X2+N),max2=*max_element(X2,X2+N);
        double res=0;
        for(int i=0;i<sx.size()-1;i++)
        {
            int a=sx[i],b=sx[i+1];
            if(a>=min1 && a<=max1 && a>=min2 && a<=max2 && b>=min1 && b<=max1 && b>=min2 && b<=max2)
                res+=simpson(a*1.0,b*1.0);
        }
        printf("%.10f
",res);
    }

    return 0;
}
原文地址:https://www.cnblogs.com/jaszzz/p/13081649.html