uvalive 4589 Asteroids

  1 #include<cstdio>
  2 #include<cmath>
  3 #include<cstdlib>
  4 #include<cstring>
  5 #include<queue>
  6 using namespace std;
  7 
  8 const double eps = 1e-8;
  9 int dcmp(double x)
 10 {
 11     if(fabs(x) < eps) return 0;
 12     else return x < 0 ? -1 : 1;
 13 }
 14 
 15 struct Point3
 16 {
 17     double x, y, z;
 18     Point3(double x=0, double y=0, double z=0):x(x),y(y),z(z) { }
 19 };
 20 
 21 typedef Point3 Vector3;
 22 
 23 Vector3 operator + (const Vector3& A, const Vector3& B)
 24 {
 25     return Vector3(A.x+B.x, A.y+B.y, A.z+B.z);
 26 }
 27 Vector3 operator - (const Point3& A, const Point3& B)
 28 {
 29     return Vector3(A.x-B.x, A.y-B.y, A.z-B.z);
 30 }
 31 Vector3 operator * (const Vector3& A, double p)
 32 {
 33     return Vector3(A.x*p, A.y*p, A.z*p);
 34 }
 35 Vector3 operator / (const Vector3& A, double p)
 36 {
 37     return Vector3(A.x/p, A.y/p, A.z/p);
 38 }
 39 
 40 bool operator == (const Point3& a, const Point3& b)
 41 {
 42     return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0 && dcmp(a.z-b.z) == 0;
 43 }
 44 
 45 Point3 read_point3()
 46 {
 47     Point3 p;
 48     scanf("%lf%lf%lf", &p.x, &p.y, &p.z);
 49     return p;
 50 }
 51 
 52 double Dot(const Vector3& A, const Vector3& B)
 53 {
 54     return A.x*B.x + A.y*B.y + A.z*B.z;
 55 }
 56 double Length(const Vector3& A)
 57 {
 58     return sqrt(Dot(A, A));
 59 }
 60 double Angle(const Vector3& A, const Vector3& B)
 61 {
 62     return acos(Dot(A, B) / Length(A) / Length(B));
 63 }
 64 Vector3 Cross(const Vector3& A, const Vector3& B)
 65 {
 66     return Vector3(A.y*B.z - A.z*B.y, A.z*B.x - A.x*B.z, A.x*B.y - A.y*B.x);
 67 }
 68 double Area2(const Point3& A, const Point3& B, const Point3& C)
 69 {
 70     return Length(Cross(B-A, C-A));
 71 }
 72 double Volume6(const Point3& A, const Point3& B, const Point3& C, const Point3& D)
 73 {
 74     return Dot(D-A, Cross(B-A, C-A));
 75 }
 76 Point3 Centroid(const Point3& A, const Point3& B, const Point3& C, const Point3& D)
 77 {
 78     return (A + B + C + D)/4.0;
 79 }
 80 
 81 double rand01()
 82 {
 83     return rand() / (double)RAND_MAX;
 84 }
 85 double randeps()
 86 {
 87     return (rand01() - 0.5) * eps;
 88 }
 89 
 90 Point3 add_noise(const Point3& p)
 91 {
 92     return Point3(p.x + randeps(), p.y + randeps(), p.z + randeps());
 93 }
 94 
 95 struct Face
 96 {
 97     int v[3];
 98     Face(int a, int b, int c)
 99     {
100         v[0] = a;
101         v[1] = b;
102         v[2] = c;
103     }
104     Vector3 Normal(const vector<Point3>& P) const
105     {
106         return Cross(P[v[1]]-P[v[0]], P[v[2]]-P[v[0]]);
107     }
108     // f是否能看见P[i]
109     int CanSee(const vector<Point3>& P, int i) const
110     {
111         return Dot(P[i]-P[v[0]], Normal(P)) > 0;
112     }
113 };
114 
115 // 增量法求三维凸包
116 // 注意：没有考虑各种特殊情况（如四点共面）。实践中，请在调用前对输入点进行微小扰动
117 vector<Face> CH3D(const vector<Point3>& P)
118 {
119     int n = P.size();
120     vector<vector<int> > vis(n);
121     for(int i = 0; i < n; i++) vis[i].resize(n);
122 
123     vector<Face> cur;
124     cur.push_back(Face(0, 1, 2)); // 由于已经进行扰动，前三个点不共线
125     cur.push_back(Face(2, 1, 0));
126     for(int i = 3; i < n; i++)
127     {
128         vector<Face> next;
129         // 计算每条边的“左面”的可见性
130         for(int j = 0; j < cur.size(); j++)
131         {
132             Face& f = cur[j];
133             int res = f.CanSee(P, i);
134             if(!res) next.push_back(f);
135             for(int k = 0; k < 3; k++) vis[f.v[k]][f.v[(k+1)%3]] = res;
136         }
137         for(int j = 0; j < cur.size(); j++)
138             for(int k = 0; k < 3; k++)
139             {
140                 int a = cur[j].v[k], b = cur[j].v[(k+1)%3];
141                 if(vis[a][b] != vis[b][a] && vis[a][b]) // (a,b)是分界线，左边对P[i]可见
142                     next.push_back(Face(a, b, i));
143             }
144         cur = next;
145     }
146     return cur;
147 }
148 
149 struct ConvexPolyhedron
150 {
151     int n;
152     vector<Point3> P, P2;
153     vector<Face> faces;
154 
155     bool read()
156     {
157         if(scanf("%d", &n) != 1) return false;
158         P.resize(n);
159         P2.resize(n);
160         for(int i = 0; i < n; i++)
161         {
162             P[i] = read_point3();
163             P2[i] = add_noise(P[i]);
164         }
165         faces = CH3D(P2);
166         return true;
167     }
168 
169     Point3 centroid()
170     {
171         Point3 C = P[0];
172         double totv = 0;
173         Point3 tot(0,0,0);
174         for(int i = 0; i < faces.size(); i++)
175         {
176             Point3 p1 = P[faces[i].v[0]], p2 = P[faces[i].v[1]], p3 = P[faces[i].v[2]];
177             double v = -Volume6(p1, p2, p3, C);
178             totv += v;
179             tot = tot + Centroid(p1, p2, p3, C)*v;
180         }
181         return tot / totv;
182     }
183 
184     double mindist(Point3 C)
185     {
186         double ans = 1e30;
187         for(int i = 0; i < faces.size(); i++)
188         {
189             Point3 p1 = P[faces[i].v[0]], p2 = P[faces[i].v[1]], p3 = P[faces[i].v[2]];
190             ans = min(ans, fabs(-Volume6(p1, p2, p3, C) / Area2(p1, p2, p3)));
191         }
192         return ans;
193     }
194 };
195 
196 int main()
197 {
198     int n, m;
199     ConvexPolyhedron P1, P2;
200     while(P1.read() && P2.read())
201     {
202         Point3 C1 = P1.centroid();
203         double d1 = P1.mindist(C1);
204         Point3 C2 = P2.centroid();
205         double d2 = P2.mindist(C2);
206         printf("%.8lf\n", d1+d2);
207     }
208     return 0;
209 }