HDU 4573 Throw the Stones 三维凸包模版应用

点击打开链接

本题的意思：n孩子轮流扔石头构成三维凸包，求出第几个孩子扔的石子使三维凸包体积增量最大并输出体积增量

  体积增量＝Vi－Vi-1 (i代表第i个孩子扔石子，Vi-1 代表前i－1个孩子构成三维凸包的体积)

1<= N <= 10^4, 1 <= i <= N, -10^4 <= x i , y i, z i <= 10^4.  (<=x i , y i, z i 代表第i个孩子扔石子的坐标 )

解题思路：如果直接暴力求增量的话复杂度n^2(1e8) 会超时。

那么我们可以根据三维凸包的增量法进行边输入点 边算出体积增量 同时维护凸包。我么知道石子 要么在现有的凸包里 要么在凸包外面，在里面的话 我们可以知道体积增量为 0，在外面我们可以发现算法维护凸包删面的过程中 这个面和这个点构成的四面体就是体积增量的一部分，所以我们只需要在删每个面的过程中求出这个四面体的体积并求和 就是第i个石子的增量最后比较出最大的就可以了

代码：

#include <iostream>
#include <cmath>
#include <cstdio>
using namespace std;
const int MAXV=1e4+1;
const double eps = 1e-9;
//三维点
struct pt {
    double x, y, z;
    pt() {}
    pt(double _x, double _y, double _z): x(_x), y(_y), z(_z) {}
    pt operator - (const pt p1) {
        return pt(x - p1.x, y - p1.y, z - p1.z);
    }
    pt operator + (const pt p1) {
        return pt(x + p1.x, y + p1.y, z + p1.z);
    }
    pt operator *(const pt p){
        return pt(y*p.z-z*p.y,z*p.x-x*p.z, x*p.y-y*p.x);
    }
    pt operator *(double d)
    {
        return pt(x*d,y*d,z*d);
    }
    
    pt operator / (double d)
    {
        return pt(x/d,y/d,z/d);
    }
    double operator ^ (pt p) {
        return x*p.x+y*p.y+z*p.z;    //点乘
    }
    double len(){
        return sqrt(x*x+y*y+z*z);
    }
}P[MAXV];;
struct pp{
    double x,y;
    pp(){}
    pp(double x,double y):x(x),y(y){}
};
int to[MAXV][MAXV];
struct _3DCH {
    struct fac {
        int a, b, c;
        bool ok;
    };
    
    int n;
    int cnt;
    fac F[MAXV*8];
    
    
    double max_va=0,diff_va=0;
    pt Cross3(pt a,pt p){
        return pt(a.y*p.z-a.z*p.y, a.z*p.x-a.x*p.z, a.x*p.y-a.y*p.x);  }
    double vlen(pt a) {
        return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
    }
    double volume(pt a, pt b, pt c, pt d) {
        return Cross3((b-a),(c-a))^(d-a);    //四面体有向体积*6
    }
    double Dot3( pt u, pt v )
    {
        return u.x * v.x + u.y * v.y + u.z * v.z;
    }
    
    //平面的法向量
    pt pvec(pt a,pt b,pt c)
    {
        return (Cross3((a-b),(b-c)));
    }
    //点到面的距离
    double Dis(pt a,pt b,pt c,pt d)
    {
        return fabs(pvec(a,b,c)^(d-a))/vlen(pvec(a,b,c));
    }
    //正：点在面同向
    double ptof(pt &p, fac &f) {
        pt m = P[f.b]-P[f.a], n = P[f.c]-P[f.a], t = p-P[f.a];
        return Cross3(m , n) ^ t;
    }
    
    void deal(int p, int a, int b) {
        int f = to[a][b];
        fac add;
        if (F[f].ok) {
            if (ptof(P[p], F[f]) > eps)
                dfs(p, f);
            else {
                add.a = b, add.b = a, add.c = p, add.ok = 1;
                to[p][b] = to[a][p] = to[b][a] = cnt;
                F[cnt++] = add;
            }
        }
    }
    
    void dfs(int p, int cur) {
        F[cur].ok = 0;
        diff_va+=fabs(volume(P[p],P[F[cur].a],P[F[cur].b],P[F[cur].c]));
        deal(p, F[cur].b, F[cur].a);
        deal(p, F[cur].c, F[cur].b);
        deal(p, F[cur].a, F[cur].c);
    }
    
    //体积
    double volume() {
        pt O(0, 0, 0);
        double ret = 0.0;
        for (int i = 0; i < cnt; i++) {
            if(F[i].ok)
            ret += volume(O, P[F[i].a], P[F[i].b], P[F[i].c]);
        }
        return fabs(ret / 6.0);
    }
    //构建三维凸包
    void construct() {
        int num=1;max_va=0;
        cnt = 0;int sb=0;
        for(int i=0;i<n;++i){
            scanf("%lf%lf%lf",&P[i].x,&P[i].y,&P[i].z);
            if(sb==0){
                if (vlen(P[0] - P[i]) > eps) {
                    swap(P[1], P[i]);
                    sb++;
                }
                continue;
            }
            if(sb==1){
                if (vlen(Cross3((P[0] - P[1]) , (P[1] - P[i]))) > eps) {
                    swap(P[2], P[i]);
                    sb++;
                }
                continue;
            }
            if(sb==2){
                if (fabs(Cross3((P[0] - P[1]) , (P[1] - P[2])) ^ (P[0] - P[i])) > eps) {
                    swap(P[3], P[i]);
                    sb++;
                }
                if(sb==3){
                    fac add;
                    num=i+1;
                    for (int k = 0;k<4;k++) {
                        add.a = (k+1)%4, add.b = (k+2)%4, add.c = (k+3)%4, add.ok = 1;
                        if(ptof(P[k], add)>0)
                            swap(add.b,add.c);
                        to[add.a][add.b]=to[add.b][add.c] = to[add.c][add.a]=cnt;
                        F[cnt++] = add;
                    }
                    for (int k = 4; k <=i; k++) {
                        for (int j = 0; j < cnt; j++) {
                            if (F[j].ok && ptof(P[k], F[j]) > eps) {
                                dfs(k, j);
                                break;
                            }
                        }
                    }
                  //  cout<<diff_va/6<<endl;
                    max_va=volume();
                   // cout<<max_va<<endl;
                    ++sb;
                }
                continue;
            }
            if(sb>3){
                diff_va=0;
                for (int j = 0; j < cnt; j++) {
                    if (F[j].ok && ptof(P[i], F[j]) > eps) {
                        dfs(i, j);
                        break;
                    }
                }
                diff_va/=6;
               // cout<<diff_va<<endl;
                if(diff_va>max_va){
                    max_va=diff_va;
                    num=i+1;
                }
            }
        }
        printf("%d %.2f\n",num,max_va);
    }
}lou;
int main(int argc, const char * argv[]) {
    int num=0;
    while(scanf("%d",&lou.n)==1){
        cout<<"Case #"<<++num<<":"<<endl;
        lou.construct();
    }
    return 0;
}