计算几何_三维凸包

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     if (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 }

hdoj3662

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
",(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
",ans);
158     }
159     return 0;
160 }

hdoj4273