計算幾何_三維凸包
阿新 • • 發佈：2018-11-11
1.hdoj3662 3D Convex Hull
傳送：http://acm.hdu.edu.cn/showproblem.php?pid=3662
題意：給出空間n個點，問凸包表面的多邊形個數。
分析：rt。
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=305; 4 const double eps=1e-8; 5 int sgn(double x){ 6 if (fabs(x)<eps) return 0; 7 ifhdoj3662(x<0) return -1; 8 return 1; 9 } 10 struct point3{ 11 double x,y,z; 12 point3(double _x=0,double _y=0,double _z=0){x=_x;y=_y;z=_z;} 13 bool operator ==(const point3 &b)const{ 14 return sgn(x-b.x)==0 && sgn(y-b.y)==0 && sgn(z-b.z)==0; 15 } 16 point3 operator-(const point3 &b)const{ 17 return point3(x-b.x,y-b.y,z-b.z); 18 } 19 //點乘 20 double operator *(const point3 &b)const{ 21 return x*b.x+y*b.y+z*b.z; 22 } 23 //叉乘 24 point3 operator ^(const point3 &b)const{ 25 return point3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);26 } 27 double len(){ 28 return sqrt(x*x+y*y+z*z); 29 } 30 }; 31 32 struct CH3D{ 33 struct face{ 34 int a,b,c; //表示凸包一個面的三個點的編號 35 bool ok; //表示該面是否屬於最終的凸包上的面 36 }; 37 int n; //初始頂點數 38 point3 p[maxn]; 39 int num; //凸包表面三角形數 40 face f[8*maxn]; //凸包表面的三角形 41 int g[maxn][maxn]; 42 point3 cross(const point3 &a,const point3 &b,const point3 &c){ 43 return (b-a)^(c-a); 44 } 45 //四面體有向面積 *6 46 double volume(point3 a,point3 b,point3 c,point3 d){ 47 return ((b-a)^(c-a))*(d-a); 48 } 49 double dblcmp(point3 &P,face &F){ 50 point3 p1=p[F.b]-p[F.a],p2=p[F.c]-p[F.a],p3=P-p[F.a]; 51 return (p1^p2)*p3; 52 } 53 void deal(int x,int a,int b){ 54 int k=g[a][b]; 55 face add; 56 if (f[k].ok){ 57 if (dblcmp(p[x],f[k])>eps) dfs(x,k); 58 else{ 59 add.a=b; add.b=a; add.c=x; add.ok=true; 60 g[x][b]=g[a][x]=g[b][a]=num; 61 f[num++]=add; 62 } 63 } 64 } 65 //遞迴搜尋應該從凸包中刪除的面 66 void dfs(int p,int now){ 67 f[now].ok=false; 68 deal(p,f[now].b,f[now].a); 69 deal(p,f[now].c,f[now].b); 70 deal(p,f[now].a,f[now].c); 71 } 72 //構建三維凸包 73 void create(){ 74 num=0; 75 face add; 76 77 //保證前四個點不共面 78 bool flag=true; 79 for (int i=1;i<n;i++) 80 if (!(p[0]==p[i])){ 81 swap(p[1],p[i]); 82 flag=false; break; 83 } 84 if (flag) return ; 85 flag=true; 86 for (int i=2;i<n;i++) 87 if (((p[1]-p[0])^(p[i]-p[0])).len()>eps){ 88 swap(p[2],p[i]); 89 flag=false; break; 90 } 91 if (flag) return ; 92 flag=true; 93 for (int i=3;i<n;i++){ 94 if (fabs(((p[1]-p[0])^(p[2]-p[0]))*(p[i]-p[0]))>eps){ 95 swap(p[3],p[i]); 96 flag=false; break; 97 } 98 } 99 if (flag) return ; 100 101 for (int i=0;i<4;i++){ 102 add.a=(i+1)%4; add.b=(i+2)%4; add.c=(i+3)%4; add.ok=true; 103 if (dblcmp(p[i],add)>0) swap(add.b,add.c); 104 g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num; 105 f[num++]=add; 106 } 107 for (int i=4;i<n;i++){ 108 for (int j=0;j<num;j++) 109 if (f[j].ok && dblcmp(p[i],f[j])>eps){ 110 dfs(i,j); 111 break; 112 } 113 } 114 int tmp=num; num=0; 115 for (int i=0;i<tmp;i++) if (f[i].ok) f[num++]=f[i]; 116 } 117 bool same(int s,int t){ 118 point3 &a=p[f[s].a],&b=p[f[s].b],&c=p[f[s].c]; 119 return fabs(volume(a,b,c,p[f[t].a]))<eps && fabs(volume(a,b,c,p[f[t].b]))<eps 120 && fabs(volume(a,b,c,p[f[t].c]))<eps; 121 } 122 //凸包上多邊形面積 123 int polygon(){ 124 int res=0; 125 for (int i=0;i<num;i++){ 126 int flag=1; 127 for (int j=0;j<i;j++) 128 if (same(i,j)){flag=0;break;} 129 res+=flag; 130 } 131 return res; 132 } 133 }; 134 CH3D C; 135 int main(){ 136 int n; double x,y,z; 137 while (cin >> n){ 138 C.n=n; 139 for (int i=0;i<n;i++) cin >> C.p[i].x >> C.p[i].y >> C.p[i].z; 140 C.create(); 141 cout << C.polygon() << endl; 142 } 143 return 0; 144 }
2.hdoj4266 The Worm in the Apple
傳送：http://acm.hdu.edu.cn/showproblem.php?pid=4266
題意：對於每個詢問，問到凸包表面的最小距離。
分析：rt。
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=1005; 4 const double eps=1e-8; 5 const double inf=1e20; 6 int sgn(double x){ 7 if (fabs(x)<eps) return 0; 8 if (x<0) return -1; 9 return 1; 10 } 11 struct point3{ 12 double x,y,z; 13 point3(double _x=0,double _y=0,double _z=0){x=_x;y=_y;z=_z;} 14 bool operator ==(const point3 &b)const{ 15 return sgn(x-b.x)==0 && sgn(y-b.y)==0 && sgn(z-b.z)==0; 16 } 17 point3 operator -(const point3 &b)const{ 18 return point3(x-b.x,y-b.y,z-b.z); 19 } 20 //點乘 21 double operator *(const point3 &b)const{ 22 return x*b.x+y*b.y+z*b.z; 23 } 24 //叉乘 25 point3 operator ^(const point3 &b)const{ 26 return point3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x); 27 } 28 double len(){ 29 return sqrt(x*x+y*y+z*z); 30 } 31 }; 32 33 struct CH3D{ 34 struct face{ 35 int a,b,c; //表示凸包一個面的三個點的編號 36 bool ok; //表示該面是否屬於最終的凸包上的面 37 }; 38 int n; //初始頂點數 39 point3 p[maxn]; 40 int num; //凸包表面三角形數 41 face f[8*maxn]; //凸包表面的三角形 42 int g[maxn][maxn]; 43 point3 cross(const point3 &a,const point3 &b,const point3 &c){ 44 return (b-a)^(c-a); 45 } 46 //四面體有向面積 *6 47 double volume(point3 a,point3 b,point3 c,point3 d){ 48 return ((b-a)^(c-a))*(d-a); 49 } 50 double dblcmp(point3 &P,face &F){ 51 point3 p1=p[F.b]-p[F.a],p2=p[F.c]-p[F.a],p3=P-p[F.a]; 52 return (p1^p2)*p3; 53 } 54 void deal(int x,int a,int b){ 55 int k=g[a][b]; 56 face add; 57 if (f[k].ok){ 58 if (dblcmp(p[x],f[k])>eps) dfs(x,k); 59 else{ 60 add.a=b; add.b=a; add.c=x; add.ok=true; 61 g[x][b]=g[a][x]=g[b][a]=num; 62 f[num++]=add; 63 } 64 } 65 } 66 //遞迴搜尋應該從凸包中刪除的面 67 void dfs(int p,int now){ 68 f[now].ok=false; 69 deal(p,f[now].b,f[now].a); 70 deal(p,f[now].c,f[now].b); 71 deal(p,f[now].a,f[now].c); 72 } 73 //構建三維凸包 74 void create(){ 75 num=0; 76 face add; 77 78 //保證前四個點不共面 79 bool flag=true; 80 for (int i=1;i<n;i++) 81 if (!(p[0]==p[i])){ 82 swap(p[1],p[i]); 83 flag=false; break; 84 } 85 if (flag) return ; 86 flag=true; 87 for (int i=2;i<n;i++) 88 if (((p[1]-p[0])^(p[i]-p[0])).len()>eps){ 89 swap(p[2],p[i]); 90 flag=false; break; 91 } 92 if (flag) return ; 93 flag=true; 94 for (int i=3;i<n;i++){ 95 if (fabs(((p[1]-p[0])^(p[2]-p[0]))*(p[i]-p[0]))>eps){ 96 swap(p[3],p[i]); 97 flag=false; break; 98 } 99 } 100 if (flag) return ; 101 102 for (int i=0;i<4;i++){ 103 add.a=(i+1)%4; add.b=(i+2)%4; add.c=(i+3)%4; add.ok=true; 104 if (dblcmp(p[i],add)>0) swap(add.b,add.c); 105 g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num; 106 f[num++]=add; 107 } 108 for (int i=4;i<n;i++){ 109 for (int j=0;j<num;j++) 110 if (f[j].ok && dblcmp(p[i],f[j])>eps){ 111 dfs(i,j); 112 break; 113 } 114 } 115 int tmp=num; num=0; 116 for (int i=0;i<tmp;i++) if (f[i].ok) f[num++]=f[i]; 117 } 118 double ptoface(point3 pp,int i){ 119 double tmp1=fabs(volume(p[f[i].a],p[f[i].b],p[f[i].c],pp)); 120 double tmp2=((p[f[i].b]-p[f[i].a])^(p[f[i].c]-p[f[i].a])).len(); 121 return tmp1/tmp2; 122 } 123 double slove(point3 pp){ 124 double mi=inf; 125 for (int i=0;i<num;i++) mi=min(mi,ptoface(pp,i)); 126 return mi; 127 } 128 }; 129 CH3D C; 130 int main(){ 131 int n,q; point3 pp; 132 while (cin >> n && n){ 133 C.n=n; 134 for (int i=0;i<n;i++) cin >> C.p[i].x >> C.p[i].y >> C.p[i].z; 135 C.create(); 136 cin >> q; 137 while (q--){ 138 cin >> pp.x >> pp.y >> pp.z; 139 double ans=inf; 140 double tmp=C.slove(pp); 141 ans=min(ans,tmp); 142 printf("%.4f\n",(ans)); 143 } 144 } 145 return 0; 146 }hdoj4266
3.hdoj4273 Rescue
傳送：http://acm.hdu.edu.cn/showproblem.php?pid=4273
題意：給出一個三維凸包，求凸包重心到凸包表面的最短距離 。
分析：rt。
1 #include<bits/stdc++.h> 2 using namespace std; 3 const int maxn=305; 4 const double eps=1e-8; 5 const double inf=1e20; 6 int sgn(double x){ 7 if (fabs(x)<eps) return 0; 8 if (x<0) return -1; 9 return 1; 10 } 11 struct point3{ 12 double x,y,z; 13 point3(double _x=0,double _y=0,double _z=0){x=_x;y=_y;z=_z;} 14 bool operator ==(const point3 &b)const{ 15 return sgn(x-b.x)==0 && sgn(y-b.y)==0 && sgn(z-b.z)==0; 16 } 17 point3 operator +(const point3 &b)const{ 18 return point3(x+b.x,y+b.y,z+b.z); 19 } 20 point3 operator -(const point3 &b)const{ 21 return point3(x-b.x,y-b.y,z-b.z); 22 } 23 point3 operator *(const double k)const{ 24 return point3(x*k,y*k,z*k); 25 } 26 point3 operator /(const double k)const{ 27 return point3(x/k,y/k,z/k); 28 } 29 //點乘 30 double operator *(const point3 &b)const{ 31 return x*b.x+y*b.y+z*b.z; 32 } 33 //叉乘 34 point3 operator ^(const point3 &b)const{ 35 return point3(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x); 36 } 37 double len(){ 38 return sqrt(x*x+y*y+z*z); 39 } 40 }; 41 42 struct CH3D{ 43 struct face{ 44 int a,b,c; //表示凸包一個面的三個點的編號 45 bool ok; //表示該面是否屬於最終的凸包上的面 46 }; 47 int n; //初始頂點數 48 point3 p[maxn]; 49 int num; //凸包表面三角形數 50 face f[8*maxn]; //凸包表面的三角形 51 int g[maxn][maxn]; 52 point3 cross(const point3 &a,const point3 &b,const point3 &c){ 53 return (b-a)^(c-a); 54 } 55 //四面體有向面積 *6 56 double volume(point3 a,point3 b,point3 c,point3 d){ 57 return ((b-a)^(c-a))*(d-a); 58 } 59 double dblcmp(point3 &P,face &F){ 60 point3 p1=p[F.b]-p[F.a],p2=p[F.c]-p[F.a],p3=P-p[F.a]; 61 return (p1^p2)*p3; 62 } 63 void deal(int x,int a,int b){ 64 int k=g[a][b]; 65 face add; 66 if (f[k].ok){ 67 if (dblcmp(p[x],f[k])>eps) dfs(x,k); 68 else{ 69 add.a=b; add.b=a; add.c=x; add.ok=true; 70 g[x][b]=g[a][x]=g[b][a]=num; 71 f[num++]=add; 72 } 73 } 74 } 75 //遞迴搜尋應該從凸包中刪除的面 76 void dfs(int p,int now){ 77 f[now].ok=false; 78 deal(p,f[now].b,f[now].a); 79 deal(p,f[now].c,f[now].b); 80 deal(p,f[now].a,f[now].c); 81 } 82 //構建三維凸包 83 void create(){ 84 num=0; 85 face add; 86 87 //保證前四個點不共面 88 bool flag=true; 89 for (int i=1;i<n;i++) 90 if (!(p[0]==p[i])){ 91 swap(p[1],p[i]); 92 flag=false; break; 93 } 94 if (flag) return ; 95 flag=true; 96 for (int i=2;i<n;i++) 97 if (((p[1]-p[0])^(p[i]-p[0])).len()>eps){ 98 swap(p[2],p[i]); 99 flag=false; break; 100 } 101 if (flag) return ; 102 flag=true; 103 for (int i=3;i<n;i++){ 104 if (fabs(((p[1]-p[0])^(p[2]-p[0]))*(p[i]-p[0]))>eps){ 105 swap(p[3],p[i]); 106 flag=false; break; 107 } 108 } 109 if (flag) return ; 110 111 for (int i=0;i<4;i++){ 112 add.a=(i+1)%4; add.b=(i+2)%4; add.c=(i+3)%4; add.ok=true; 113 if (dblcmp(p[i],add)>0) swap(add.b,add.c); 114 g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num; 115 f[num++]=add; 116 } 117 for (int i=4;i<n;i++){ 118 for (int j=0;j<num;j++) 119 if (f[j].ok && dblcmp(p[i],f[j])>eps){ 120 dfs(i,j); 121 break; 122 } 123 } 124 int tmp=num; num=0; 125 for (int i=0;i<tmp;i++) if (f[i].ok) f[num++]=f[i]; 126 } 127 //求凸包重心 128 point3 barycenter(){ 129 point3 ans=point3(0,0,0),o=point3(0,0,0); 130 double all=0; 131 for (int i=0;i<num;i++){ 132 double vol=volume(o,p[f[i].a],p[f[i].b],p[f[i].c]); 133 ans=ans+(((o+p[f[i].a]+p[f[i].b]+p[f[i].c])/4.0)*vol); 134 all+=vol; 135 } 136 ans=ans/all; 137 return ans; 138 } 139 //點到平面距離 140 double ptoface(point3 pp,int i){ 141 double tmp1=fabs(volume(p[f[i].a],p[f[i].b],p[f[i].c],pp)); 142 double tmp2=((p[f[i].b]-p[f[i].a])^(p[f[i].c]-p[f[i].a])).len(); 143 return tmp1/tmp2; 144 } 145 }; 146 CH3D C; 147 int main(){ 148 int n; double x,y,z; 149 while (cin >> n){ 150 C.n=n; 151 for (int i=0;i<n;i++) cin >> C.p[i].x >> C.p[i].y >> C.p[i].z; 152 C.create(); 153 point3 pp=C.barycenter(); 154 double ans=inf; 155 for (int i=0;i<C.num;i++) 156 ans=min(ans,C.ptoface(pp,i)); 157 printf("%.3f\n",ans); 158 } 159 return 0; 160 }hdoj4273