admin 发表于 2024-3-16 09:04:56

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]
查看完整版本: n维空间的多面体的有向测度和重心