n维空间的多面体的有向测度和重心
//#define LOCAL#pragma GCC optimize(2)
#pragma G++ optimize(2)
#pragma warning(disable:4996)
#pragma comment(linker, "/STACK:1024000000,1024000000")
#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <string>
#include <ctype.h>
#include <string.h>
#include <fstream>
#include <sstream>
#include <math.h>
#include <map>
//#include <unordered采用map>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <set>
#include <time.h>
#include <stdlib.h>
#include <bitset>
using namespace std;
//#define int unsigned long long
//#define int long long
#define re register int
#define ci const int
#define ui unsigned int
typedef pair<int, int> P;
#define FE(cur) for(re h = head, to; ~h; h = g.nxt)
#define ilv inline void
#define ili inline int
#define ilc inline char
#define ild inline double
#define ilp inline P
#define LEN(cur) (hjt.r - hjt.l)
#define MID(cur) (hjt.l + hjt.r >> 1)
#define SQUARE(x) ((x) * (x))
typedef vector<int>::iterator vit;
typedef set<int>::iterator sit;
typedef map<int, int>::iterator mit;
const int inf = ~0u >> 1;
//const int inf = 0x3f3f3f3f;
const double PI = acos(-1.0), eps = 1e-8;
namespace fastio
{
const int BUF = 1 << 21;
char fr, fw, * pr1 = fr, * pr2 = fr; int pw;
ilc gc() { return pr1 == pr2 && (pr2 = (pr1 = fr) + fread(fr, 1, BUF, stdin), *pr2 = 0, pr1 == pr2) ? EOF : *pr1++; }
ilv flush() { fwrite(fw, 1, pw, stdout); pw = 0; }
ilv pc(char c) { if (pw >= BUF) flush(); fw = c; }
ili read(int& x)
{
x = 0; int f = 1; char c = gc(); if (!~c) return EOF;
while (!isdigit(c)) { if (c == '-') f = -1; c = gc(); }
while (isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = gc();
x *= f; return 1;
}
ili read(double& x)
{
int xx = 0; double f = 1.0, fraction = 1.0; char c = gc(); if (!~c) return EOF;
while (!isdigit(c)) { if (c == '-') f = -1.0; c = gc(); }
while (isdigit(c)) { xx = (xx << 3) + (xx << 1) + (c ^ 48), c = gc(); }
x = xx;
if (c ^ '.') { x = f * xx; return 1; }
c = gc();
while (isdigit(c)) x += (c ^ 48) * (fraction /= 10), c = gc();
x *= f; return 1;
}
ilv write(int x) { if (x < 0) pc('-'), x = -x; if (x > 9) write(x / 10); pc(x % 10 + 48); }
ilv writeln(int x) { write(x); pc(10); }
ili read(char* x)
{
char c = gc(); if (!~c) return EOF;
while (!isalpha(c) && !isdigit(c)) c = gc();
while (isalpha(c) || isdigit(c)) *x++ = c, c = gc();
*x = 0; return 1;
}
ili readln(char* x)
{
char c = gc(); if (!~c) return EOF;
while (c == 10) c = gc();
while (c >= 32 && c <= 126) *x++ = c, c = gc();
*x = 0; return 1;
}
ilv write(char* x) { while (*x) pc(*x++); }
ilv write(const char* x) { while (*x) pc(*x++); }
ilv writeln(char* x) { write(x); pc(10); }
ilv writeln(const char* x) { write(x); pc(10); }
ilv write(char c) { pc(c); }
ilv writeln(char c) { write(c); pc(10); }
} using namespace fastio;
#define cp const Point
#define cs const Surface
ci maxn = 105;
ili dbcmp(double x)
{
if (fabs(x) < eps)
{
return 0;
}
return x > eps ? 1 : -1;
}
ilv mmin(double& a, double b)
{
if (a > b)
{
a = b;
}
}
struct CH3D
{
struct Point
{
double x, y, z;
Point(double x = 0, double y = 0, double z = 0) :x(x), y(y), z(z) {}
Point operator + (cp& o) const
{
return Point(x + o.x, y + o.y, z + o.z);
}
Point operator - (cp& o) const
{
return Point(x - o.x, y - o.y, z - o.z);
}
double operator * (cp& o) const
{
return x * o.x + y * o.y + z * o.z;
}
Point operator / (cp& o) const
{
return Point(y * o.z - z * o.y, z * o.x - x * o.z, x * o.y - y * o.x);
}
double magnitude() const
{
return sqrt(SQUARE(x) + SQUARE(y) + SQUARE(z));
}
Point scalar(double lambda) const
{
return Point(x * lambda, y * lambda, z * lambda);
}
} ps;
int n;
struct Surface
{
int a, b, c, flag;
Surface(int a = 0, int b = 0, int c = 0) :a(a), b(b), c(c) { flag = 1; }
} surfaces;
int num;
int g;
ild area(cp& a, cp& b, cp& c)
{
return ((b - a) / (c - a)).magnitude() / 2.0;
}
ild area(cs& sur)
{
return area(ps, ps, ps);
}
ild area()
{
double ans = 0;
if (n == 3)
{
return area(ps, ps, ps);
}
for (re i = 0; i < num; i++)
{
ans += area(ps.a], ps.b], ps.c]);
}
return ans;
}
ild volume(cp& a, cp& b, cp& c, cp& d)
{
return (((b - a) / (c - a)) * (d - a)) / 6.0;
}
ild volume(cp& p, cs& sur)
{
return volume(ps, ps, ps, p);
}
ild volume()
{
double ans = 0;
Point o;
for (re i = 0; i < num; i++)
{
ans += volume(o, ps.a], ps.b], ps.c]);
}
return ans;
}
ili ck(cp& p, cs& sur)
{
return dbcmp(volume(ps, ps, ps, p)) == 1;
}
ili onsurface(cp& p, cs& sur)
{
return !dbcmp(volume(ps, ps, ps, p));
}
ili same(cs& s, cs& t)
{
return onsurface(ps, t) && onsurface(ps, t) && onsurface(ps, t);
}
ili prework()
{
if (n < 4)
{
return 0;
}
int ok = 1;
for (re i = 1; i < n; i++)
{
if (dbcmp((ps - ps).magnitude()))
{
swap(ps, ps);
ok = 0;
break;
}
}
if (ok)
{
return 0;
}
ok = 1;
for (re i = 2; i < n; i++)
{
if (dbcmp(area(ps, ps, ps)))
{
swap(ps, ps);
ok = 0;
break;
}
}
if (ok)
{
return 0;
}
ok = 1;
for (re i = 3; i < n; i++)
{
if (dbcmp(volume(ps, ps, ps, ps)))
{
swap(ps, ps);
ok = 0;
break;
}
}
return !ok;
}
ilv construct()
{
num = 0;
if (!prework())
{
return;
}
Surface add;
for (re i = 0; i < 4; i++)
{
Surface add((i + 1) % 4, (i + 2) % 4, (i + 3) % 4);
if (ck(ps, add))
{
swap(add.b, add.c);
}
g = g = g = num;
surfaces = add;
}
for (re i = 4; i < n; i++)
{
for (re j = 0; j < num; j++)
{
if (surfaces.flag && ck(ps, surfaces))
{
dfs(i, j);
break;
}
}
}
int tmp = num; num = 0;
for (re i = 0; i < tmp; i++)
{
if (surfaces.flag)
{
surfaces = surfaces;
}
}
}
void dfs(int p, int j)
{
surfaces.flag = 0;
rmv(p, surfaces.b, surfaces.a);
rmv(p, surfaces.c, surfaces.b);
rmv(p, surfaces.a, surfaces.c);
}
void rmv(int p, int a, int b)
{
int f = g;
if (surfaces.flag)
{
if (ck(ps, surfaces))
{
dfs(p, f);
}
else
{
Surface add(b, a, p);
g = g = g = num;
surfaces = add;
}
}
}
ild dis(cp& p, cs& sur)
{
return fabs(volume(p, sur)) / area(sur) * 3.0;
}
ild dis(cp& p)
{
double ans = inf;
for (int i = 0; i < num; i++)
{
mmin(ans, dis(p, surfaces));
}
return ans;
}
Point gg()
{
Point ans, o = ps;
double v = 0, t;
for (re i = 0; i < num; i++)
{
cp& a = ps.a];
cp& b = ps.b];
cp& c = ps.c];
t = volume(o, surfaces);
ans = ans + (a + b + c + o).scalar(t / 4);
v += t;
}
return ans.scalar(1 / v);
}
}ch;
signed main()
{
#ifdef LOCAL
FILE* ALLIN = freopen("d:\\data.in", "r", stdin);
freopen("d:\\my.out", "w", stdout);
#endif
while (~scanf("%d", &ch.n))
//while (~read(ch.n))
{
// for (re i = 0; i < ch.n; i++) read(ch.ps.x), read(ch.ps.y), read(ch.ps.z);
for (re i = 0; i < ch.n; i++) scanf("%lf%lf%lf", &ch.ps.x, &ch.ps.y, &ch.ps.z);
ch.construct();
CH3D::Point centroid = ch.gg();
printf("%.3lf\n", ch.dis(centroid));
}
flush();
#ifdef LOCAL
fclose(ALLIN);
#endif
return 0;
}
https://cloud.tencent.com/developer/article/1695459?areaId=106001
页:
[1]