python 三维凸包_hdu 4266 三维凸包（增量法）

The Worm in the Apple

Time Limit: 50000/20000 MS (Java/Others)    Memory Limit: 65536/65536 K (Java/Others) Total Submission(s): 7    Accepted Submission(s): 5

Problem Description

Willy the Worm was living happily in an apple – until some vile human picked the apple, and started to eat it! Now, Willy must escape! Given a description of the apple (defined as a convex shape in 3D space), and a list of possible positions in the apple for Willy (defined as 3D points), determine the minimum distance Willy must travel to get to the surface of the apple from each point.

Input

There will be several test cases in the input. Each test case will begin with a line with a single integer n (4≤n≤1,000), which tells the number of points describing the apple. On the next n lines will be three integers x, y and z (-10,000≤x,y,z≤10,000), where each point (x,y,z) is either on the surface of, or within, the apple. The apple is the convex hull of these points. No four points will be coplanar. Following the description of the apple, there will be a line with a single integer q (1≤q≤100,000), which is the number of queries – that is, the number of points where Willy might be inside the apple. Each of the following q lines will contain three integers x, y and z (-10,000≤x,y,z≤10,000), representing a point (x,y,z) where Willy might be. All of Willy’s points are guaranteed to be inside the apple. The input will end with a line with a single 0.

Output

For each query, output a single floating point number, indicating the minimum distance Willy must travel to exit the apple. Output this number with exactly 4 decimal places of accuracy, using standard 5 up / 4 down rounding (e.g. 2.12344 rounds to 2.1234, 2.12345 rounds to 2.1235). Output each number on its own line, with no spaces, and do not print any blank lines between answers.

Sample Input

6

0 0 0

100 0 0

0 100 0

0 0 100

20 20 20

30 20 10

4

1 1 1

30 30 35

7 8 9

90 2 2

0

Sample Output

1.0000 2.8868 7.0000 2.0000

题意：求三维凸包中的点到凸包的最短距离

思路：利用增量算法求出三维凸包的每个面，再用点到面的距离，暴力找出最小的距离。

三维凸包的增量法：

初始时需要一个四面体。可以先找两个不同点P1，P2，寻找和它们不共线的第三个点P3，再找不共面的第四个点P4,。如果找不到，则调用二维凸包算法。

接下来计算剩下点的随机排列。每次加一个点，有两种情况：

情况1：新点在当前凸包内部，只需简单地忽略该点，如图1所示。

情况2：新点在当前凸包外部，需要计算并的凸包，在这种情况下，首先需要计算原凸包相对于Pr的水平面，即Pr可以看到的封闭区域，如图2所示。

将P点能看到的所有平面删去，同时求出相对于Pr的水平面，将水平面上的边与Pr相连组成面，如图3所示。

判断一个点P在凸包内还是凸包外，可利用有向体积的方法。在存储面时，保证面的法线方向朝向凸包外部，若存在某一平面和点P所组成的四面体的有向体积为正，则P点在凸包外部，并且此面可被P点看见。此时，只要将此面删去并用同样的方法判断与它相邻的其他平面是否可被P点看见，可以使用深度优先搜索。

此算法对于每个点需要求出它所能看到的所有平面，最简单的方法便是枚举所有当前凸包上的平面，一一判断，此时算法时间复杂度是O(n^2)。此外还有更高效的方法来寻找P点的可见区域，算法的时间复杂度为O(nlog(n))。

View Code

1 #include

2 #include

3 #include

4 #include

5 #include

6 using namespacestd;7 #define PR 1e-9

8 #define N 1100

9 structTPoint10 {11 doublex,y,z;12 TPoint(){}13 TPoint(double _x,double _y,double_z):x(_x),y(_y),z(_z){}14 TPoint operator-(const TPoint p) {return TPoint(x-p.x,y-p.y,z-p.z);}15 TPoint operator*(const TPoint p) {return TPoint(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);}//叉积

16 double operator^(const TPoint p) {return x*p.x+y*p.y+z*p.z;}//点积

17 };18 TPoint dd;19 struct fac//20 {21 int a,b,c;//凸包一个面上的三个点的编号

22 bool ok;//该面是否是最终凸包中的面

23 };24 TPoint xmult(TPoint u,TPoint v)25 {26 return TPoint(u.y*v.z-v.y*u.z,u.z*v.x-u.x*v.z,u.x*v.y-u.y*v.x);27 }28 doubledmult(TPoint u,TPoint v)29 {30 return u.x*v.x+u.y*v.y+u.z*v.z;31 }32 TPoint subt(TPoint u,TPoint v)33 {34 return TPoint(u.x-v.x,u.y-v.y,u.z-v.z);35 }36 doublevlen(TPoint u)37 {38 return sqrt(u.x*u.x+u.y*u.y+u.z*u.z);39 }40 TPoint pvec(TPoint a,TPoint b,TPoint c)41 {42 returnxmult(subt(a,b),subt(b,c));43 }44 doubleDis(TPoint a,TPoint b,TPoint c,TPoint d)45 {46 return fabs(dmult(pvec(a,b,c),subt(d,a)))/vlen(pvec(a,b,c));47 }48 structT3dhull49 {50 int n;//初始点数

51 TPoint ply[N];//初始点

52 int trianglecnt;//凸包上三角形数

53 fac tri[N];//凸包三角形

54 int vis[N][N];//点i到点j是属于哪个面

55 double dist(TPoint a){return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);}//两点长度

56 double area(TPoint a,TPoint b,TPoint c){return dist((b-a)*(c-a));}//三角形面积*2

57 double volume(TPoint a,TPoint b,TPoint c,TPoint d){return (b-a)*(c-a)^(d-a);}//四面体有向体积*6

58 double ptoplane(TPoint &p,fac &f)//正：点在面同向

59 {60 TPoint m=ply[f.b]-ply[f.a],n=ply[f.c]-ply[f.a],t=p-ply[f.a];61 return (m*n)^t;62 }63 void deal(int p,int a,intb)64 {65 int f=vis[a][b];//与当前面(cnt)共边(ab)的那个面

66 fac add;67 if(tri[f].ok)68 {69 if((ptoplane(ply[p],tri[f]))>PR) dfs(p,f);//如果p点能看到该面f，则继续深度探索f的3条边，以便更新新的凸包面

70 else//否则因为p点只看到cnt面，看不到f面，则p点和a、b点组成一个三角形。

71 {72 add.a=b,add.b=a,add.c=p,add.ok=1;73 vis[p][b]=vis[a][p]=vis[b][a]=trianglecnt;74 tri[trianglecnt++]=add;75 }76 }77 }78 void dfs(int p,int cnt)//维护凸包，如果点p在凸包外更新凸包

79 {80 tri[cnt].ok=0;//当前面需要删除，因为它在更大的凸包里面81 //下面把边反过来(先b,后a)，以便在deal()中判断与当前面(cnt)共边(ab)的那个面。即判断与当头面(cnt)相邻的3个面(它们与当前面的共边是反向的，如下图中(1)的法线朝外(即逆时针)的面130和312,它们共边13，但一个方向是13,另一个方向是31)

82

83 deal(p,tri[cnt].b,tri[cnt].a);84 deal(p,tri[cnt].c,tri[cnt].b);85 deal(p,tri[cnt].a,tri[cnt].c);86 }87 bool same(int s,int e)//判断两个面是否为同一面

88 {89 TPoint a=ply[tri[s].a],b=ply[tri[s].b],c=ply[tri[s].c];90 return fabs(volume(a,b,c,ply[tri[e].a]))

95 {96 inti,j;97 trianglecnt=0;98 if(n<4) return;99 bool tmp=true;100 for(i=1;i

101 {102 if((dist(ply[0]-ply[i]))>PR)103 {104 swap(ply[1],ply[i]); tmp=false; break;105 }106 }107 if(tmp) return;108 tmp=true;109 for(i=2;i

110 {111 if((dist((ply[0]-ply[1])*(ply[1]-ply[i])))>PR)112 {113 swap(ply[2],ply[i]); tmp=false; break;114 }115 }116 if(tmp) return;117 tmp=true;118 for(i=3;i

119 {120 if(fabs((ply[0]-ply[1])*(ply[1]-ply[2])^(ply[0]-ply[i]))>PR)121 {122 swap(ply[3],ply[i]); tmp=false; break;123 }124 }125 if(tmp) return;126 fac add;127 for(i=0;i<4;i++)//构建初始四面体(4个点为ply[0],ply[1],ply[2],ply[3])

128 {129 add.a=(i+1)%4,add.b=(i+2)%4,add.c=(i+3)%4,add.ok=1;130 if((ptoplane(ply[i],add))>0) swap(add.b,add.c);//保证逆时针，即法向量朝外，这样新点才可看到。

131 vis[add.a][add.b]=vis[add.b][add.c]=vis[add.c][add.a]=trianglecnt;//逆向的有向边保存

132 tri[trianglecnt++]=add;133 }134 for(i=4;i

135 {136 for(j=0;j

137 {138 if(tri[j].ok&&(ptoplane(ply[i],tri[j]))>PR)//对当前凸包面进行判断，看是否点能否看到这个面

139 {140 dfs(i,j); break;//点能看到当前面，更新凸包的面(递归，可能不止更新一个面)。当前点更新完成后break跳出循环

141 }142 }143 }144 int cnt=trianglecnt;//这些面中有一些tri[i].ok=0，它们属于开始建立但后来因为在更大凸包内故需删除的，所以下面几行代码的作用是只保存最外层的凸包

145 trianglecnt=0;146 for(i=0;inow) _min=now;159 }160 return_min;161 }162 }hull;163 intmain()164 {165 doublenow,_min;166 while(scanf("%d",&hull.n)!=EOF)167 {168 if(hull.n==0) break;169 inti,j,q;170

171 for(i=0;i