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));
}