UVA1438 Asteroids(增量法求三维凸包,加权所有三棱锥质量求多面体重心)
https://www.luogu.com.cn/problem/UVA1438
题解建议参考lrj白书。
代码来自牛逼网友。
我只是存个板子。
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 }