Three-dimensional convex hull template for surface area and volume of three-dimensional convex hull

#include <stdio.h> #include <string.h> #include <stdlib.h> #include <math.h> #include <iost  ream> #include <queue> #include <algorithm> using namespace std;      #define PR 1e-8 #define N 510 struct tpoint{double x, y, Z;  Tpoint () {} tpoint (double _x, double _y, double _z): X (_x), Y (_y), Z (_z) {} tpoint operator-(const Tpoint p) {return      Tpoint (x-p.x, Y-P.Y, z-p.z);}      Tpoint operator* (const Tpoint p) {return tpoint (y*p.z-z*p.y, z*p.x-x*p.z, x*p.y-y*p.x);}  Double operator^ (const Tpoint p) {return x*p.x+y*p.y+z*p.z;}  };      struct fac{int A, b, C;  BOOL OK;  };      struct t3dhull{int n;      Tpoint Ply[n];      int trianglecnt;      FAC tri[n];      int vis[n][n];      Double Dist (Tpoint a) {return sqrt (A.X*A.X+A.Y*A.Y+A.Z*A.Z);}      Double area (Tpoint A, tpoint B, Tpoint c) {return dist ((b-a) * (c-a));}      Double Volume (Tpoint A, tpoint B, Tpoint C, Tpoint D) {return (b-a) * (c-a) ^ (d-a);} DoubLe Ptoplane (tpoint &p, FAC &f) {tpoint m = ply[f.b]-PLY[F.A], n = ply[f.c]-ply[f.a], t = p-ply[f.          A];      Return (m*n) ^t;          } void Deal (int p, int a, int b) {int f = vis[a][b];          FAC add;              if (Tri[f].ok) {if ((Ptoplane (Ply[p], tri[f]) > PR) DFS (P, f);                  else {add.a = b, add.b = a, ADD.C = p, Add.ok = 1;                  VIS[P][B] = vis[a][p] = Vis[b][a] = trianglecnt;              tri[trianglecnt++] = add;          }}} void Dfs (int p, int cnt) {Tri[cnt].ok = 0;          Deal (P, tri[cnt].b, TRI[CNT].A);          Deal (P, tri[cnt].c, tri[cnt].b);      Deal (P, tri[cnt].a, TRI[CNT].C);          } bool Same (int s, int e) {Tpoint a = ply[tri[s].a], B = ply[tri[s].b], c = ply[tri[s].c]; Return Fabs (Volume (A,B,C,PLY[TRI[E].A)) < PR && Fabs (Volume (a,b,c,ply[tri[e].b]) < PR && Fabs (Volume (A,B,C,PLY[TRI[E].C)) < PR;          } void construct () {int i, J;          trianglecnt = 0;          if (n<4) return;          BOOL TMP = TRUE; for (i = 1; i < n; i++) {if ((Dist (ply[0]-ply[i))) > PR) {Swap (P                  LY[1], ply[i]);                  TMP = FALSE;              Break          }} if (TMP) return;          TMP = TRUE;                  for (i = 2; i < n; i++) {if ((Dist ((ply[0]-ply[1]) * (ply[1]-ply[i))) > PR) {                  Swap (ply[2], ply[i]);                  TMP = FALSE;              Break          }} if (TMP) return;          TMP = TRUE;              for (i = 3; i < n; i++) {if (Fabs ((ply[0]-ply[1)) * (ply[1]-ply[2]) ^ (ply[0]-ply[i]) &GT;PR)                  {Swap (ply[3], ply[i]);               TMP =false;   Break          }} if (TMP) return;          FAC add;              for (i = 0; i < 4; i++) {add.a = (i+1)%4, add.b = (i+2)%4, add.c = (i+3)%4, Add.ok = 1;              if ((Ptoplane (Ply[i], add)) >0) swap (add.b, ADD.C);              VIS[ADD.A][ADD.B] = vis[add.b][add.c] = vis[add.c][add.a] = trianglecnt;          tri[trianglecnt++] = add;                  } for (i = 4; i < n; i++) {for (j = 0; J < trianglecnt; J + +) { if (Tri[j].ok && (Ptoplane (Ply[i], tri[j])) > PR) {dfs (i, j); Brea                  K          }}} int cnt = TRIANGLECNT;          trianglecnt = 0;          for (i = 0; i < cnt; i++) {if (Tri[i].ok) tri[trianglecnt++] = Tri[i];          }} double area () {double ret = 0; for (int i = 0; i < trianglecnt;          i++) ret + = Area (PLY[TRI[I].A], ply[tri[i].b], ply[tri[i].c]);      return ret/2.0;          } double Volume () {Tpoint P (0,0,0);          DOUBLE ret = 0;          for (int i = 0; i < trianglecnt; i++) ret + = Volume (P, ply[tri[i].a], ply[tri[i].b], ply[tri[i].c]);      return Fabs (RET/6);    }}hull;      int main () {int Cas = 1;          while (scanf ("%d", &AMP;HULL.N), HULL.N) {int i; for (i = 0; i < HULL.N; i++) scanf ("%lf%lf%lf", &hull.ply[i].x, &hull.ply[i].y, &hull.ply[i].z          );          Hull.construct ();      printf ("Case%d:%.2lf\n", cas++, Hull.area ());  } return 0;   }

Three-dimensional convex hull template for surface area and volume of three-dimensional convex hull