1. 程式人生 > >計算幾何_三維凸包

計算幾何_三維凸包

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\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