n维空间的多面体的有向测度和重心

admin

1. //#define LOCAL
2. #pragma GCC optimize(2)
3. #pragma G++ optimize(2)
4. #pragma warning(disable:4996)
5. #pragma comment(linker, "/STACK:1024000000,1024000000")
6. #include <stdio.h>
7. #include <iostream>
8. #include <iomanip>
9. #include <string>
10. #include <ctype.h>
11. #include <string.h>
12. #include <fstream>
13. #include <sstream>
14. #include <math.h>
15. #include <map>
16. //#include <unordered采用map>
17. #include <algorithm>
18. #include <vector>
19. #include <stack>
20. #include <queue>
21. #include <set>
22. #include <time.h>
23. #include <stdlib.h>
24. #include <bitset>
25. using namespace std;
26. //#define int unsigned long long
27. //#define int long long
28. #define re register int
29. #define ci const int
30. #define ui unsigned int
31. typedef pair<int, int> P;
32. #define FE(cur) for(re h = head[cur], to; ~h; h = g[h].nxt)
33. #define ilv inline void
34. #define ili inline int
35. #define ilc inline char
36. #define ild inline double
37. #define ilp inline P
38. #define LEN(cur) (hjt[cur].r - hjt[cur].l)
39. #define MID(cur) (hjt[cur].l + hjt[cur].r >> 1)
40. #define SQUARE(x) ((x) * (x))
41. typedef vector<int>::iterator vit;
42. typedef set<int>::iterator sit;
43. typedef map<int, int>::iterator mit;
44. const int inf = ~0u >> 1;
45. //const int inf = 0x3f3f3f3f;
46. const double PI = acos(-1.0), eps = 1e-8;
47. namespace fastio
48. {
49.     const int BUF = 1 << 21;
50.     char fr[BUF], fw[BUF], * pr1 = fr, * pr2 = fr; int pw;
51.     ilc gc() { return pr1 == pr2 && (pr2 = (pr1 = fr) + fread(fr, 1, BUF, stdin), *pr2 = 0, pr1 == pr2) ? EOF : *pr1++; }
52.     ilv flush() { fwrite(fw, 1, pw, stdout); pw = 0; }
53.     ilv pc(char c) { if (pw >= BUF) flush(); fw[pw++] = c; }
54.     ili read(int& x)
55.     {
56.         x = 0; int f = 1; char c = gc(); if (!~c) return EOF;
57.         while (!isdigit(c)) { if (c == '-') f = -1; c = gc(); }
58.         while (isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = gc();
59.         x *= f; return 1;
60.     }
61.     ili read(double& x)
62.     {
63.         int xx = 0; double f = 1.0, fraction = 1.0; char c = gc(); if (!~c) return EOF;
64.         while (!isdigit(c)) { if (c == '-') f = -1.0; c = gc(); }
65.         while (isdigit(c)) { xx = (xx << 3) + (xx << 1) + (c ^ 48), c = gc(); }
66.         x = xx;
67.         if (c ^ '.') { x = f * xx; return 1; }
68.         c = gc();
69.         while (isdigit(c)) x += (c ^ 48) * (fraction /= 10), c = gc();
70.         x *= f; return 1;
71.     }
72.     ilv write(int x) { if (x < 0) pc('-'), x = -x; if (x > 9) write(x / 10); pc(x % 10 + 48); }
73.     ilv writeln(int x) { write(x); pc(10); }
74.     ili read(char* x)
75.     {
76.         char c = gc(); if (!~c) return EOF;
77.         while (!isalpha(c) && !isdigit(c)) c = gc();
78.         while (isalpha(c) || isdigit(c)) *x++ = c, c = gc();
79.         *x = 0; return 1;
80.     }
81.     ili readln(char* x)
82.     {
83.         char c = gc(); if (!~c) return EOF;
84.         while (c == 10) c = gc();
85.         while (c >= 32 && c <= 126) *x++ = c, c = gc();
86.         *x = 0; return 1;
87.     }
88.     ilv write(char* x) { while (*x) pc(*x++); }
89.     ilv write(const char* x) { while (*x) pc(*x++); }
90.     ilv writeln(char* x) { write(x); pc(10); }
91.     ilv writeln(const char* x) { write(x); pc(10); }
92.     ilv write(char c) { pc(c); }
93.     ilv writeln(char c) { write(c); pc(10); }
94. } using namespace fastio;
95. #define cp const Point
96. #define cs const Surface
97. ci maxn = 105;
99. ili dbcmp(double x)
100. {
101.     if (fabs(x) < eps)
102.     {
103.         return 0;
104.     }
105.     return x > eps ? 1 : -1;
106. }
108. ilv mmin(double& a, double b)
109. {
110.     if (a > b)
111.     {
112.         a = b;
113.     }
114. }
116. struct CH3D
117. {
118.     struct Point
119.     {
120.         double x, y, z;
121.         Point(double x = 0, double y = 0, double z = 0) :x(x), y(y), z(z) {}
123.         Point operator + (cp& o) const
124.         {
125.             return Point(x + o.x, y + o.y, z + o.z);
126.         }
128.         Point operator - (cp& o) const
129.         {
130.             return Point(x - o.x, y - o.y, z - o.z);
131.         }
133.         double operator * (cp& o) const
134.         {
135.             return x * o.x + y * o.y + z * o.z;
136.         }
138.         Point operator / (cp& o) const
139.         {
140.             return Point(y * o.z - z * o.y, z * o.x - x * o.z, x * o.y - y * o.x);
141.         }
143.         double magnitude() const
144.         {
145.             return sqrt(SQUARE(x) + SQUARE(y) + SQUARE(z));
146.         }
148.         Point scalar(double lambda) const
149.         {
150.             return Point(x * lambda, y * lambda, z * lambda);
151.         }
153.     } ps[maxn];
154.     int n;
156.     struct Surface
157.     {
158.         int a, b, c, flag;
159.         Surface(int a = 0, int b = 0, int c = 0) :a(a), b(b), c(c) { flag = 1; }
160.     } surfaces[maxn << 3];
161.     int num;
163.     int g[maxn][maxn];
165.     ild area(cp& a, cp& b, cp& c)
166.     {
167.         return ((b - a) / (c - a)).magnitude() / 2.0;
168.     }
170.     ild area(cs& sur)
171.     {
172.         return area(ps[sur.a], ps[sur.b], ps[sur.c]);
173.     }
175.     ild area()
176.     {
177.         double ans = 0;
178.         if (n == 3)
179.         {
180.             return area(ps[0], ps[1], ps[2]);
181.         }
182.         for (re i = 0; i < num; i++)
183.         {
184.             ans += area(ps[surfaces[i].a], ps[surfaces[i].b], ps[surfaces[i].c]);
185.         }
186.         return ans;
187.     }
189.     ild volume(cp& a, cp& b, cp& c, cp& d)
190.     {
191.         return (((b - a) / (c - a)) * (d - a)) / 6.0;
192.     }
194.     ild volume(cp& p, cs& sur)
195.     {
196.         return volume(ps[sur.a], ps[sur.b], ps[sur.c], p);
197.     }
199.     ild volume()
200.     {
201.         double ans = 0;
202.         Point o;
203.         for (re i = 0; i < num; i++)
204.         {
205.             ans += volume(o, ps[surfaces[i].a], ps[surfaces[i].b], ps[surfaces[i].c]);
206.         }
207.         return ans;
208.     }
210.     ili ck(cp& p, cs& sur)
211.     {
212.         return dbcmp(volume(ps[sur.a], ps[sur.b], ps[sur.c], p)) == 1;
213.     }
215.     ili onsurface(cp& p, cs& sur)
216.     {
217.         return !dbcmp(volume(ps[sur.a], ps[sur.b], ps[sur.c], p));
218.     }
220.     ili same(cs& s, cs& t)
221.     {
222.         return onsurface(ps[s.a], t) && onsurface(ps[s.b], t) && onsurface(ps[s.c], t);
223.     }
225.     ili prework()
226.     {
227.         if (n < 4)
228.         {
229.             return 0;
230.         }
231.         int ok = 1;
232.         for (re i = 1; i < n; i++)
233.         {
234.             if (dbcmp((ps[0] - ps[i]).magnitude()))
235.             {
236.                 swap(ps[1], ps[i]);
237.                 ok = 0;
238.                 break;
239.             }
240.         }
241.         if (ok)
242.         {
243.             return 0;
244.         }
245.         ok = 1;
246.         for (re i = 2; i < n; i++)
247.         {
248.             if (dbcmp(area(ps[0], ps[1], ps[i])))
249.             {
250.                 swap(ps[2], ps[i]);
251.                 ok = 0;
252.                 break;
253.             }
254.         }
255.         if (ok)
256.         {
257.             return 0;
258.         }
259.         ok = 1;
260.         for (re i = 3; i < n; i++)
261.         {
262.             if (dbcmp(volume(ps[0], ps[1], ps[2], ps[i])))
263.             {
264.                 swap(ps[3], ps[i]);
265.                 ok = 0;
266.                 break;
267.             }
268.         }
269.         return !ok;
270.     }
272.     ilv construct()
273.     {
274.         num = 0;
275.         if (!prework())
276.         {
277.             return;
278.         }
279.         Surface add;
280.         for (re i = 0; i < 4; i++)
281.         {
282.             Surface add((i + 1) % 4, (i + 2) % 4, (i + 3) % 4);
283.             if (ck(ps[i], add))
284.             {
285.                 swap(add.b, add.c);
286.             }
287.             g[add.a][add.b] = g[add.b][add.c] = g[add.c][add.a] = num;
288.             surfaces[num++] = add;
289.         }
290.         for (re i = 4; i < n; i++)
291.         {
292.             for (re j = 0; j < num; j++)
293.             {
294.                 if (surfaces[j].flag && ck(ps[i], surfaces[j]))
295.                 {
296.                     dfs(i, j);
297.                     break;
298.                 }
299.             }
300.         }
301.         int tmp = num; num = 0;
302.         for (re i = 0; i < tmp; i++)
303.         {
304.             if (surfaces[i].flag)
305.             {
306.                 surfaces[num++] = surfaces[i];
307.             }
308.         }
309.     }
311.     void dfs(int p, int j)
312.     {
313.         surfaces[j].flag = 0;
314.         rmv(p, surfaces[j].b, surfaces[j].a);
315.         rmv(p, surfaces[j].c, surfaces[j].b);
316.         rmv(p, surfaces[j].a, surfaces[j].c);
317.     }
319.     void rmv(int p, int a, int b)
320.     {
321.         int f = g[a][b];
322.         if (surfaces[f].flag)
323.         {
324.             if (ck(ps[p], surfaces[f]))
325.             {
326.                 dfs(p, f);
327.             }
328.             else
329.             {
330.                 Surface add(b, a, p);
331.                 g[b][a] = g[a][p] = g[p][b] = num;
332.                 surfaces[num++] = add;
333.             }
334.         }
335.     }
337.     ild dis(cp& p, cs& sur)
338.     {
339.         return fabs(volume(p, sur)) / area(sur) * 3.0;
340.     }
342.     ild dis(cp& p)
343.     {
344.         double ans = inf;
345.         for (int i = 0; i < num; i++)
346.         {
347.             mmin(ans, dis(p, surfaces[i]));
348.         }
349.         return ans;
350.     }
352.     Point gg()
353.     {
354.         Point ans, o = ps[0];
355.         double v = 0, t;
356.         for (re i = 0; i < num; i++)
357.         {
358.             cp& a = ps[surfaces[i].a];
359.             cp& b = ps[surfaces[i].b];
360.             cp& c = ps[surfaces[i].c];
361.             t = volume(o, surfaces[i]);
362.             ans = ans + (a + b + c + o).scalar(t / 4);
363.             v += t;
364.         }
365.         return ans.scalar(1 / v);
366.     }
368. }ch;
370. signed main()
371. {
372. #ifdef LOCAL
373.     FILE* ALLIN = freopen("d:\\data.in", "r", stdin);
374.     freopen("d:\\my.out", "w", stdout);
375. #endif
376.     while (~scanf("%d", &ch.n))
377. //  while (~read(ch.n))
378.     {
379. //      for (re i = 0; i < ch.n; i++) read(ch.ps[i].x), read(ch.ps[i].y), read(ch.ps[i].z);
380.         for (re i = 0; i < ch.n; i++) scanf("%lf%lf%lf", &ch.ps[i].x, &ch.ps[i].y, &ch.ps[i].z);
381.         ch.construct();
382.         CH3D::Point centroid = ch.gg();
383.         printf("%.3lf\n", ch.dis(centroid));
384.     }
385.     flush();
386. #ifdef LOCAL
387.     fclose(ALLIN);
388. #endif
389.     return 0;
390. }
392. https://cloud.tencent.com/developer/article/1695459?areaId=106001

复制代码