bzoj 1964: hull 三维凸包 计算几何

1964: hull 三维凸包

Time Limit: 1 Sec  Memory Limit: 64 MB
Submit: 54  Solved: 39
[Submit][Status][Discuss]

Description

三维凸包问题是一个基础的三维计算几何问题，不过这次你只需要做一个简单版的三维凸包问题就行了。

Input

输入数据一共有若干行，每行三个整数，表示一个点的坐标。点的个数为五个。

Output

输出一个实数，保留两位小数，表示三维凸包的体积。

Sample Input

0 0 0
2 0 0
0 2 0
2 2 0
1 1 1

Sample Output

1.33

HINT

对于所有数据，坐标范围在[0,100]*[0,100]*[0,100]。

 

　　直接水过，不过网上的三维凸包的标程一个比一个不靠谱，由一组数据拍出来，加上我的程序，三个程序三个答案。整个人都不好了。

　　枚举一个面，使得剩下两个点与这条面组成四面体的有向体积符号不同，则两点在面的两侧，计算即可。

　　四面体的计算公式是三个向量组成的矩阵行列式的1/6.

　　这是居然我的第一个三维算几题。

　　另外，这个方法的正确性我还是比较怀疑，因为存在数据是他和网上的三维算几模板矛盾，虽然大部分可以看出是网上的程序的问题，但是居然他与我自己写的版有矛盾，就姑且把它当成骗分暴力吧。（最后有三维凸包程序）

#include<iostream>
#include<cstring>
#include<algorithm>
#include<cstdio>
#include<cmath>
using namespace std;
#define MAXN 10
typedef double real;
struct point
{
        real x,y,z;
        point(real x,real y,real z):x(x),y(y),z(z){}
        point(){}
        void read()
        {
                scanf("%lf %lf %lf",&x,&y,&z);
        }
}pl[MAXN];;
real area(point p1,point p2,point p3)
{
        return (p1.x*p2.y*p3.z + p1.z*p2.x*p3.y + p1.y*p2.z*p3.x
                -p1.x*p2.z*p3.y - p1.z*p2.y*p3.x - p1.y*p2.x*p3.z)/6;
}
point operator -(point p1,point p2)
{
        return point(p2.x-p1.x,p2.y-p1.y,p2.z-p1.z);
}
typedef point line;

bool vis[MAXN];
int main()
{
        freopen("input.txt","r",stdin);
        for (int i=0;i<5;i++)
        {
                pl[i].read();
        }
//        printf("%.2lf\n",area(pl[1]-pl[0],pl[2]-pl[0],pl[3]-pl[0]));
        real ans=0;
        real s1,s2,s3;
        for (int a=0;a<5;a++)
        {
                vis[a]=true;
                for (int b=0;b<5;b++)
                {
                        if (vis[b])continue;
                        vis[b]=true;
                        for (int c=0;c<5;c++)
                        {
                                if (vis[c])continue;
                                vis[c]=true;
                                for (int d=0;d<5;d++)
                                {
                                        if (vis[d])continue;
                                        vis[d]=true;
                                        for (int e=0;e<5;e++)
                                        {
                                                if (vis[e])continue;
                                                vis[e]=true;
                                                s3=abs(s1=area(pl[c]-pl[a],pl[b]-pl[a],pl[d]-pl[a]))+abs(s2=area(pl[c]-pl[a],pl[b]-pl[a],pl[e]-pl[a]));
                                                if (s1*s2>0)
                                                {
                                                        vis[e]=false;continue;
                                                }
                                                if (s1==0 && s2==0)
                                                {
                                                        vis[e]=false;continue;
                                                }
                                                ans=max(ans,s3);
                                                vis[e]=false;
                                        }
                                        vis[d]=false;
                                }
                                vis[c]=false;
                        }
                        vis[b]=false;
                }
                vis[a]=false;
        }
        printf("%.2lf\n",ans);
}

 

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<assert.h>
#include<vector>
using namespace std;
#define MAXN 120
#define eps 1e-8
typedef double real;
inline int sgn(real x)
{
        if (abs(x)<eps)return 0;
        return x<0?-1:1;
}
struct point
{
        real x,y,z;
        point(real x,real y,real z):x(x),y(y),z(z){}
        point(){}
        int read()
        {
                return scanf("%lf %lf %lf\n",&x,&y,&z);
        }
        void noise()
        {
                x+=(real)(rand()%10000-5000)/10000*eps;
                y+=(real)(rand()%10000-5000)/10000*eps;
                z+=(real)(rand()%10000-5000)/10000*eps;
        }
        real len()
        {
                return sqrt(x*x+y*y+z*z);
        }
};
typedef point vect;
struct line
{
        point ps;
        real x,y,z;
        line(){}
        line(point p1,point p2)
        {
                ps=p1;
                x=p2.x-p1.x;
                y=p2.y-p1.y;
                z=p2.z-p1.z;
        }
        line(point ps,real x,real y,real z):ps(ps),x(x),y(y),z(z){}
};
vect xmul(line l1,line l2)
{
        return vect(l1.y*l2.z-l1.z*l2.y,l1.z*l2.x-l1.x*l2.z,l1.x*l2.y-l1.y*l2.x);
}
real volume(line l1,line l2,line l3)
{
        return + l1.x*l2.y*l3.z - l1.x*l2.z*l3.y 
               - l1.y*l2.x*l3.z + l1.y*l2.z*l3.x
               + l1.z*l2.x*l3.y - l1.z*l2.y*l3.x;
}
//+1 2 3
//-1 3 2
//-2 1 3
//+2 3 1
//+3 1 2
//-3 2 1
struct surface
{
        point ps;
        real x1,y1,z1;
        real x2,y2,z2;
        surface(){}
        surface(point p1,point p2,point p3)
        {
                ps=p1;
                x1=p2.x-p1.x,y1=p2.y-p1.y,z1=p2.z-p1.z;
                x2=p3.x-p1.x,y2=p3.y-p1.y,z2=p3.z-p1.z;
        }
        real volume(point pt)
        {
                return ::volume(line(ps,pt),line(ps,x1,y1,z1),line(ps,x2,y2,z2));
        }
        vect nvect()
        {
                return xmul(line(ps,x1,y1,z1),line(ps,x2,y2,z2));
        }
        void reverse()
        {
                swap(x1,x2);
                swap(y1,y2);
                swap(z1,z2);
        }
};
point pl[MAXN];
struct face
{
        int pt[3];
        face(int x,int y,int z)
        {
                pt[0]=x;pt[1]=y;pt[2]=z;
        }
        surface ToSurface()
        {
                return surface(pl[pt[0]],pl[pt[1]],pl[pt[2]]);
        }
        void print()
        {
                printf("Face:%d %d %d\n",pt[0],pt[1],pt[2]);
        }
};
vector<face> cc;
vector<pair<int,int> > chs;
bool status[MAXN][MAXN];
int main()
{
        freopen("input.txt","r",stdin);
        int n=0;
        while (~pl[n++].read());
        n--;
        for (int i=0;i<n;i++)
                pl[i].noise();
    /*    for (int i=0;i<n;i++)
                swap(pl[rand()%n],pl[rand()%n]);*/
        cc.push_back(face(0,1,2));
        cc.push_back(face(2,1,0));
        for (int i=0;i<3;i++)
                status[i][(i+1)%3]=true;
        for (int i=1;i<4;i++)
                status[i%3][i-1]=true;
        for (int i=3;i<n;i++)
        {
                //for (int j=0;j<cc.size();j++)cc[j].print();    printf("\n");
                chs.clear();
                for (int j=0;j<cc.size();j++)
                {
                        if (cc[j].ToSurface().volume(pl[i])<0)
                        {
                                for (int k=0;k<3;k++)
                                {
                                        status[cc[j].pt[k]][cc[j].pt[(k+1)%3]]=false;
                                        chs.push_back(make_pair(cc[j].pt[k],cc[j].pt[(k+1)%3]));
                                }
                                swap(cc[j],cc[cc.size()-1]);
                                cc.pop_back();
                                j--;
                        }
                }
                for (int j=0;j<chs.size();j++)
                {
                        if (!status[chs[j].first][chs[j].second] && status[chs[j].second][chs[j].first])continue;
                        chs[j]=chs[chs.size()-1];
                        j--;
                        chs.pop_back();
                }
                for (int j=0;j<chs.size();j++)
                {
                        cc.push_back(face(i,chs[j].first,chs[j].second));
                        status[i][chs[j].first]=status[chs[j].first][chs[j].second]=status[chs[j].second][i]=true;
                }
                for (int j=0;j<n;j++)
                        for (int k=0;k<n;k++)
                                assert(!(status[j][k]^status[k][j]));
        }
        //for (int j=0;j<cc.size();j++)cc[j].print();    printf("\n");
        real ans=0;
        for (int i=0;i<cc.size();i++)
        {
                ans+=cc[i].ToSurface().volume(point(0,0,0))/6;
        }
        printf("%.2lf\n",abs(ans));
}