三维凸包模板题……只是听了听计算几何的课之后心血来潮想写的……
我的做法很无脑是吧……暴力枚举三个点组成的三角形,然后枚举剩下的点,判断其余点是否都在这个三角形的同一侧,是的话则说明这个三角形是凸包的一个面。
理论复杂度应该是$O(n^4)$,不过看上去跑得飞快?人帅自带小常数哈哈
这个故事告诉我们:大力出奇迹……
1 #include<cstdio> 2 #include<cstring> 3 #include<cmath> 4 #include<algorithm> 5 #include<iostream> 6 #include<fstream> 7 #include<iomanip> 8 using namespace std; 9 const int maxn=110; 10 const long double eps=1e-15; 11 struct Point{ 12 long double x,y,z; 13 Point(long double x=0.0,long double y=0.0,long double z=0.0):x(x),y(y),z(z){} 14 Point operator-(const Point &a)const{return Point(a.x-x,a.y-y,a.z-z);}//A - B = B到A的位移 15 Point operator/(const long double &a)const{return Point(x/a,y/a,z/a);} 16 }a[maxn]; 17 typedef Point Vector; 18 long double noise(); 19 long double Dot(const Vector&,const Vector&); 20 Vector Cross(const Vector&,const Vector&); 21 long double Length(const Vector&); 22 long double Area2(const Point&,const Point&,const Point&); 23 Vector LawVector(const Vector&,const Vector&); 24 long double Distance(const Point&,const Point&,const Vector&); 25 ifstream fin("enwrap.in"); 26 ofstream fout("enwrap.out"); 27 long double ans=0.0; 28 int n,t; 29 int main(){ 30 fin>>n; 31 for(int i=1;i<=n;i++){ 32 fin>>a[i].x>>a[i].y>>a[i].z; 33 a[i].x+=noise(); 34 a[i].y+=noise(); 35 a[i].z+=noise(); 36 } 37 Vector A; 38 long double d; 39 bool bad; 40 for(int i=1;i<=n;i++)for(int j=1;j<i;j++)for(int k=1;k<j;k++){ 41 bad=false; 42 t=0; 43 A=LawVector(a[j]-a[i],a[k]-a[i]); 44 for(int l=1;l<=n;l++)if(i!=l&&j!=l&&k!=l){ 45 d=Distance(a[l],a[i],A); 46 if(!t)t=(d<0?-1:1); 47 else if(t*d<0){ 48 bad=true; 49 break; 50 } 51 } 52 if(!bad)ans+=Area2(a[i],a[j],a[k]); 53 } 54 ans/=2.0; 55 fout<<fixed<<ans; 56 return 0; 57 } 58 long double noise(){ 59 static int a=1213,b=678424137,p=998244353,x=753497501; 60 x=a*x+b;x%=p; 61 if(x<0)x+=p; 62 return (long double)(x-p)/p*5e-10; 63 } 64 inline long double Dot(const Vector &A,const Vector &B){return A.x*B.x+A.y*B.y+A.z*B.z;} 65 inline Vector Cross(const Vector &A,const Vector &B){return Vector(A.y*B.z-B.y*A.z,A.z*B.x-B.z*A.x,A.x*B.y-B.x*A.y);} 66 inline long double Length(const Vector &A){return sqrt(Dot(A,A));} 67 inline long double Area2(const Point &A,const Point &B,const Point &C){return Length(Cross(B-A,C-A));} 68 inline Vector LawVector(const Vector &A,const Vector &B){ 69 Vector n=Cross(A,B); 70 return n/Length(n); 71 }//两个向量的叉积一定同时垂直于这两个向量 72 inline long double Distance(const Point &A,const Point &P,const Vector &n){return Dot(A-P,n);}//(P,n)是平面的点法式,n是单位向量 73 /* 74 三维凸包——暴力法 75 暴力枚举三个点组成的三角形, 76 判断其他点是否都在这个三角形的同侧, 77 是则说明这个三角形一定是凸包的一个面。 78 */
话说一开始忘了怎么求平面的法向量,手推了个公式然后发现搞出NaN了……因为我推的公式默认法向量在z维的长度不为0,但其实有很多平面的法向量是垂直于z轴的,然后就除零爆炸了……又脑补了很久,后来才想起来两个向量的叉积必定同时垂直于这两个向量……这人没救了
还有一件事,这题卡精度,随机扰乱搞得太大会炸精度……计算几何毁我青春