UVA1438 Asteroids（增量法求三维凸包，加权所有三棱锥质量求多面体重心）

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <cmath>
  4 #include <cstdlib>
  5 #include <vector>
  6 #include <algorithm>
  7  
  8 using namespace std;
  9 const double eps = 1e-9;
 10  
 11 inline int dcmp (double x) { if (fabs(x) < eps) return 0; else return x < 0 ? -1 : 1; }
 12  
 13 struct Point3 {
 14     double x, y, z;
 15  
 16     Point3 (double x = 0, double y = 0, double z = 0): x(x), y(y), z(z) {}
 17     bool operator < (const Point3& u) const { return dcmp(x-u.x)<0 || (dcmp(x-u.x)==0 && dcmp(y-u.y)<0) || (dcmp(x-u.x)==0 && dcmp(y-u.y)==0 && dcmp(z-u.z) < 0); }
 18     bool operator > (const Point3& u) const { return u < (*this); }
 19     bool operator == (const Point3& u) const { return !(u < (*this) || (*this) < u); }
 20     bool operator != (const Point3& u) const { return !((*this) == u); }
 21     Point3 operator + (const Point3& u) const { return Point3(x+u.x, y+u.y, z+u.z); }
 22     Point3 operator - (const Point3& u) const { return Point3(x-u.x, y-u.y, z-u.z); }
 23     Point3 operator * (const double u) const { return Point3(x*u, y*u, z*u); }
 24     Point3 operator / (const double u) const { return Point3(x/u, y/u, z/u); }
 25  
 26     void read () { scanf("%lf%lf%lf", &x, &y, &z); }
 27 };
 28  
 29 typedef Point3 Vector3;
 30  
 31 namespace Vectorial {
 32     double getDot(Vector3 a, Vector3 b) { return a.x*b.x + a.y*b.y + a.z*b.z; }
 33     double getLength(Vector3 a) { return sqrt(getDot(a, a)); }
 34     double getAngle(Vector3 a, Vector3 b) { return acos(getDot(a, b) / getLength(a) / getLength(b)); }
 35     Vector3 getCross (Vector3 a, Vector3 b) { 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); }
 36     Vector3 getNormal(Point3 a, Point3 b, Point3 c) {  
 37         Vector3 u = a-b, v = b-c;
 38         Vector3 k = getCross(u, v);
 39         return k / getLength(k);
 40     }
 41     double getDistancePointToPlane (Point3 p, Point3 p0, Vector3 v) { return fabs(getDot(p-p0, v)); }
 42     Point3 getPlaneProjection (Point3 p, Point3 p0, Vector3 v) { return p - v * getDot(p-p0, v); }
 43 };
 44  
 45 namespace Linear {
 46     using namespace Vectorial;
 47     double getDistancePointToLine(Point3 p, Point3 a, Point3 b) {
 48         Vector3 v1 = b-a, v2 = p-a;
 49         return getLength(getCross(v1,v2)) / getLength(v1);
 50     }
 51     double getDistancePointToSegment(Point3 p, Point3 a, Point3 b) {
 52         if (a == b) return getLength(p-a);
 53         Vector3 v1 = b-a, v2 = p-a, v3 = p-b;
 54         if (dcmp(getDot(v1, v2)) < 0) return getLength(v2);
 55         else if (dcmp(getDot(v1, v3)) > 0) return getLength(v3);
 56         else return getLength(getCross(v1, v2)) / getLength(v1);
 57     }
 58  
 59     bool getPointLineToLine (Point3 a, Vector3 u, Point3 b, Vector3 v, double& s) {
 60         double p = getDot(u, u) * getDot(v, v) - getDot(u, v) * getDot(u, v);
 61         if (dcmp(p) == 0) return false;
 62         double q = getDot(u, v) * getDot(v, a-b) - getDot(v, v) * getDot(u, a-b);
 63         s = p/q;
 64         return true;
 65     }
 66  
 67     double getDistanceLineToLine (Point3 a, Vector3 u, Point3 b, Vector3 v) {
 68         double s, t;
 69         bool flag1 = getPointLineToLine(a, u, b, v, s);
 70         bool flag2 = getPointLineToLine(b, v, a, u, t);
 71         if (flag1 && flag2) {
 72             Point3 p = a + u * s, q = b + v * t;
 73             return getLength(p-q);
 74         }
 75         return 0;
 76     }
 77  
 78     double getDistanceSegmentToSegment(Point3 a, Point3 b, Point3 c, Point3 d) {
 79         double s, t;
 80         bool flag1 = getPointLineToLine(a, b-a, c, d-c, s);
 81         bool flag2 = getPointLineToLine(c, d-c, a, b-a, t);
 82         if (flag1 && flag2 && dcmp(s) > 0 && dcmp(s - 1) < 0 && dcmp(t) > 0 && dcmp(t-1) < 0) {
 83             Vector3 u = b-a, v = d-c;
 84             Point3 p = a + u * s, q = b + v * t;
 85             return getLength(p-q);
 86         } else {
 87             double ans = 1e20;
 88             ans = min(ans, getDistancePointToSegment(a, c, d));
 89             ans = min(ans, getDistancePointToSegment(b, c, d));
 90             ans = min(ans, getDistancePointToSegment(c, a, b));
 91             ans = min(ans, getDistancePointToSegment(d, a, b));
 92             return ans;
 93         }
 94     }
 95 };
 96  
 97 namespace Triangular {
 98     using namespace Vectorial;
 99     double getArea (Point3 a, Point3 b, Point3 c) { return getLength(getCross(b-a, c-a)); }
100     bool onTriangle (Point3 p, Point3 a, Point3 b, Point3 c) {
101         double area1 = getArea(p, a, b);
102         double area2 = getArea(p, b, c);
103         double area3 = getArea(p, c, a);
104         return dcmp(area1 + area2 + area3 - getArea(a, b, c)) == 0;
105     }
106     bool haveIntersectionTriSeg (Point3 p0, Point3 p1, Point3 p2, Point3 a, Point3 b, Point3& p) {
107         Vector3 v = getCross(p1-p0, p2-p0);
108         if (dcmp(getDot(v, b-a)) == 0) return false;
109         else {
110             double t = getDot(v, p0 - a) / getDot(v, b-a);
111             if (dcmp(t) < 0 || dcmp(t-2) > 0) return false;
112             p = a + (b-a) * t;
113             return onTriangle(p, p0, p1, p2);
114         }
115     }
116 };
117  
118 struct Face {
119     int v[3];
120     Face (int a = 0, int b = 0, int c = 0) { v[0] = a, v[1] = b, v[2] = c;}
121     Vector3 normal (Point3 *p) const { return Vectorial::getCross(p[v[1]] - p[v[0]], p[v[2]]-p[v[0]]); }
122     int cansee (Point3 *p, int i) const {
123         return Vectorial::getDot(p[i]-p[v[0]], normal(p)) > 0 ? 1 : 0;
124     }
125 };
126  
127 namespace Polygonal {
128     using namespace Vectorial;
129  
130     double getVolume (Point3 a, Point3 b, Point3 c, Point3 d) { return getDot(d-a, getCross(b-a, c-a)) / 6; }
131  
132     int vis[1005][1005];
133     double rand01() { return rand() / (double) RAND_MAX; }
134     double randeps() { return (rand01() - 0.5) * eps; }
135     Point3 addNoise(Point3 p) { return Point3(p.x+randeps(), p.y+randeps(), p.z+randeps()); }
136     vector<Face> CH3D (Point3 *o, int n, Point3* p) {
137         for (int i = 0; i < n; i++) p[i] = addNoise(o[i]);
138  
139         memset(vis, -1, sizeof(vis));
140         vector<Face> cur;
141         cur.push_back(Face(0, 1, 2));
142         cur.push_back(Face(2, 1, 0));
143  
144         for (int i = 3; i < n; i++) {
145             vector<Face> net;
146             for (int j = 0; j < cur.size(); j++) {
147                 Face& f = cur[j];
148                 int res = f.cansee(p, i);
149                 if (!res) net.push_back(f);
150                 for (int k = 0; k < 3; k++) vis[f.v[k]][f.v[(k+1)%3]] = res;
151             }
152  
153             for (int j = 0; j < cur.size(); j++) {
154                 for (int k = 0; k < 3; k++) {
155                     int a = cur[j].v[k], b = cur[j].v[(k+1)%3];
156                     if (vis[a][b] != vis[b][a] && vis[a][b])
157                         net.push_back(Face(a, b, i));
158                 }
159             }
160             cur = net;
161         }
162         return cur;
163     }
164  
165     Point3 getCenter (const vector<Face>& f, Point3* p) {
166         int n = f.size();
167         double sv = 0, sx = 0, sy = 0, sz = 0;
168         for (int i = 0; i < n; i++) {
169             double v = getVolume(Point3(0, 0, 0), p[f[i].v[0]], p[f[i].v[1]], p[f[i].v[2]]);
170             sv += v;
171             sx += (p[f[i].v[0]].x + p[f[i].v[1]].x + p[f[i].v[2]].x) * v;
172             sy += (p[f[i].v[0]].y + p[f[i].v[1]].y + p[f[i].v[2]].y) * v;
173             sz += (p[f[i].v[0]].z + p[f[i].v[1]].z + p[f[i].v[2]].z) * v;
174         }
175         return Point3(sx/sv/4, sy/sv/4, sz/sv/4);
176     }
177 };
178  
179  
180 using namespace Vectorial;
181 using namespace Polygonal;
182 const int maxn = 105;
183 const double inf = 1e20;
184  
185 int N, M;
186 vector<Face> Poly1, Poly2;
187 Point3 P[maxn], Q[maxn], T[maxn];
188  
189 int main () {
190     while (scanf("%d", &N) == 1) {
191         for (int i = 0; i < N; i++) P[i].read();
192         scanf("%d", &M);
193         for (int i = 0; i < M; i++) Q[i].read();
194         sort(P, P + N);
195         sort(Q, Q + M);
196         N = unique(P, P + N) - P;
197         M = unique(Q, Q + M) - Q;
198  
199         Poly1 = CH3D(P, N, T);
200         Poly2 = CH3D(Q, M, T);
201         Point3 center1 = getCenter(Poly1, P), center2 = getCenter(Poly2, Q);
202         double dis1 = inf, dis2 = inf;
203         for (int i = 0; i < Poly1.size(); i++) {
204             int a = Poly1[i].v[0], b = Poly1[i].v[1], c = Poly1[i].v[2];
205             dis1 = min(dis1, getDistancePointToPlane(center1, P[a], getNormal(P[a], P[b], P[c])));
206         }
207         for (int i = 0; i < Poly2.size(); i++) {
208             int a = Poly2[i].v[0], b = Poly2[i].v[1], c = Poly2[i].v[2];
209             dis2 = min(dis2, getDistancePointToPlane(center2, Q[a], getNormal(Q[a], Q[b], Q[c])));
210         }
211         printf("%.6lf\n", dis1 + dis2);
212     }
213     return 0;
214 }