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<stdio.h> 2 #include<string.h> 3 #include<math.h> 4 #include<algorithm> 5 #include <iostream> 6 using namespace std; 7 #define PR 1e-9 8 #define N 1100 9 struct TPoint 10 { 11 double x,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 double dmult(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 double vlen(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 return xmult(subt(a,b),subt(b,c)); 43 } 44 double Dis(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 struct T3dhull 49 { 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,int b) 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]))<PR 91 &&fabs(volume(a,b,c,ply[tri[e].b]))<PR 92 &&fabs(volume(a,b,c,ply[tri[e].c]))<PR; 93 } 94 void construct()//構建凸包 95 { 96 int i,j; 97 trianglecnt=0; 98 if(n<4) return ; 99 bool tmp=true; 100 for(i=1;i<n;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<n;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<n;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<n;i++)//構建更新凸包 135 { 136 for(j=0;j<trianglecnt;j++)//對每個點判斷是否在當前3維凸包內或外(i表示當前點,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;i<cnt;i++) 147 { 148 if(tri[i].ok) 149 tri[trianglecnt++]=tri[i]; 150 } 151 } 152 double res() 153 { 154 double _min=1e300; 155 for(int i=0;i<trianglecnt;i++) 156 { 157 double now=Dis(ply[tri[i].a],ply[tri[i].b],ply[tri[i].c],dd); 158 if(_min>now) _min=now; 159 } 160 return _min; 161 } 162 }hull; 163 int main() 164 { 165 double now,_min; 166 while(scanf("%d",&hull.n)!=EOF) 167 { 168 if(hull.n==0) break; 169 int i,j,q; 170 171 for(i=0;i<hull.n;i++) 172 scanf("%lf%lf%lf",&hull.ply[i].x,&hull.ply[i].y,&hull.ply[i].z); 173 hull.construct(); 174 scanf("%d",&q); 175 for(j=0;j<q;j++) 176 { 177 scanf("%lf%lf%lf",&dd.x,&dd.y,&dd.z); 178 printf("%.4lf\n",hull.res()); 179 } 180 } 181 return 0; 182 }