Three-dimensional computational geometry template

On the net to find a three-dimensional computational geometry template, perfect a bit, so that it can use ...

#include <cstdio> #include <cstring> #include <algorithm>using namespace std;/*********** Foundation ********    /const double eps=0.000001;typedef struct POINT_3D {Double x, y, Z; Point_3d (Double xx = 0, double yy = 0, double zz = 0): X (xx), Y (yy), Z (ZZ) {} bool operator = = (Const point_3d& A)    const {return x==a.x && y==a.y && z==a.z; }}vector_3d;    Point_3d read_point_3d () {Double x, y, Z;    scanf ("%lf%lf%lf", &x,&y,&z); return point_3d (x, y, z);} VECTOR_3D operator + (const VECTOR_3D & A, const VECTOR_3D & B) {return vector_3d (a.x + b.x, A.y + b.y, A.z + B.Z);} VECTOR_3D operator-(const POINT_3D & A, const POINT_3D & B) {return vector_3d (a.x-b.x, A.y-b.y, A.z-b. z);} VECTOR_3D operator * (const VECTOR_3D & A, double p) {return vector_3d (a.x * p, A.Y * p, A.Z * p);} VECTOR_3D operator/(const VECTOR_3D & A, double p) {return vector_3d (a.x/p, a.y/p, a.z/p);} Double Dot (const VECTOR_3D & A, const VECTOR_3D & B) {return a.x * b.x + a.y * b.y + a.z * B.Z;} Double Length (const VECTOR_3D & a) {return sqrt (Dot (a));} Double Angle (const VECTOR_3D & A, const VECTOR_3D & B) {return ACOs (Dot (A, b)/Length (a)/length (B));}  VECTOR_3D Cross (const VECTOR_3D & A, const VECTOR_3D & B) {return vector_3d (A.Y * b.z-a.z * b.y, A.Z * b.x- a.x * b.z, a.x * b.y-a.y * b.x);} Double Area2 (const POINT_3D & A, const POINT_3D & B, const POINT_3D & C) {return Length (b-a, c-a ));}  Double Volume6 (const POINT_3D & A, const POINT_3D & B, const POINT_3D & C, const POINT_3D & D) {return Dot (D-a, Cross (b-a, c-a));}  Tetrahedron center of Gravity point_3d centroid (const POINT_3D & A, const POINT_3D & B, const POINT_3D & C, const POINT_3D & D) {return (A + B + C + D)/4.0;} /************ Point line *************///the distance from the point p to the plane p0-n. n must be a unit vector double Distancetoplane (const POINT_3D & P, const POINT_3D & P0, conSt Vector_3d & N) {return fabs (Dot (P-P0, N)), or//if the absolute value is not obtained, a projection of the direction distance}//point p on the plane p0-n is given.  n must be a unit vector point_3d getplaneprojection (const POINT_3D & P, const POINT_3D & p0, const VECTOR_3D & N) {return P -N * Dot (P-P0, n);}  Line P1-P2 to the intersection of the plane p0-n point_3d lineplaneintersection (point_3d p1, point_3d p2, point_3d p0, vector_3d N) {vector_3d v= p2    -P1; Double t = (dot (n, p0-p1)/dot (n, p2-p1)); The denominator is 0, the line is parallel to the plane or return p1 + V * t on the plane; If the line segment determines if T is between 0~1}//point P to line ab distance double distancetoline (const POINT_3D & P, const POINT_3D & A, const POINT_3D &amp ;    B) {Vector_3d V1 = b-a, v2 = p-a; Return Length (Cross (v1, v2))/Length (v1);}    Distance from point to segment Double distancetoseg (point_3d p, point_3d A, point_3d b) {if (a = = b) {return Length (P-A);    } vector_3d v1 = b-a, v2 = p-a, V3 = p-b;    if (Dot (v1, v2) + EPS < 0) {return Length (v2);      } else {if (Dot (v1, v3)-EPS > 0) {return Length (V3);  } else {return Length (cross (v1, v2))/Length (v1); }}}//the p1+s*u of the opposite plane corresponds to P2+t*v's crossover if parallel | coincident, return falsebool Linedistance3d (point_3d p1, vector_3d u, point_3d p2, vector_3d v    , double & s) {double b = dot (u, u) * dot (V, v)-dot (U, v) * DOT (U, v);    if (ABS (b) <= EPS) {return false;    Double A = dot (u, v) * DOT (V, p1-p2)-dot (V, v) * DOT (U, P1-P2);    s = A/b; return true;} P1 and P2 are the same side bool Sameside (const POINT_3D & p1, const POINT_3D & p2, const POINT_3D & A, const point_3d on line B) & B) {return Dot (cross (b-a, p1-a), Cross (B-a, p2-a))-EPS >= 0;} Point P in Triangle P0, P1, P2 bool Pointintri (const POINT_3D & P, const POINT_3D & P0, const POINT_3D & P1, const POINT_ 3D & P2) {return sameside (p, P0, P1, P2) && sameside (p, P1, P0, P2) && sameside (p, P2, P0, P1);} Whether triangle P0P1P2 and line segments AB intersect bool Trisegintersection (const POINT_3D & P0, const POINT_3D & P1, const POINT_3D & P2, const POINT_3D & A, const POINT_3D & B, Point_3d & P) {vector_3d n = Cross (p1-p0, p2-p0);    if (ABS (Dot (n, b-a)) <= EPS) {return false;        Segment A-B and planar p0p1p2 parallel or coplanar} else//plane A and line p1-p2 have unique intersection {double T = dot (n, p0-a)/dot (n, b-a);    if (t + EPS < 0 | | t-1-EPS > 0) {return false; Not on line ab} P = A + (b-a) * t;    Intersection return Pointintri (P, P0, P1, P2);    }}//space two or three whether the angle intersects BOOL Tritriintersection (POINT_3D * T1, point_3d * T2) {point_3d P;            for (int i = 0; i < 3; i++) {if (Trisegintersection (t1[0], t1[1], t1[2], t2[i], t2[(i + 1)% 3], P)) {        return true;        } if (Trisegintersection (t2[0], t2[1], t2[2], t1[i], t1[(i + 1)% 3], P)) {return true; }} return false;}  Space two straight lines on the nearest point pair return nearest distance point pair saved in ans1 Ans2 double segsegdistance (point_3d A1, Point_3d B1, point_3d A2, point_3d B2, point_3d& ANS1, point_3d& ans2) {vector_3d V1 = (A1-B1), V2 = (A2-B2);    Vector_3d N = Cross (v1, v2);    Vector_3d ab = (A1-A2);    Double ans = Dot (n, AB)/Length (n);    point_3d p1 = a1, p2 = A2;    vector_3d D1 = b1-a1, d2 = B2-A2;    Double T1, T2;    T1 = Dot (Cross (P2-P1, D2), Cross (D1, D2));    T2 = Dot (Cross (P2-P1, D1), Cross (D1, D2));    Double DD = Length (Cross (d1, D2));    T1/= DD * DD;    T2/= DD * DD;    ans1 = (A1 + (B1-A1) * T1);    Ans2 = (A2 + (B2-A2) * T2); Return fabs (ANS);} Determine if p is in triangles A, B, C, and the distance to three edges is at least mindist. Guaranteed P, A, B, C coplanar bool Insidewithmindistance (const POINT_3D & P, const POINT_3D & A, const POINT_3D & B, const POI Nt_3d & C, double mindist) {if (!    Pointintri (P, A, B, C)) {return false;    } if (Distancetoline (P, A, B) < Mindist) {return false;    } if (Distancetoline (P, B, C) < Mindist) {return false;    } if (Distancetoline (P, C, A) < Mindist) {return false;    }return true;} Determine if p is in convex quadrilateral abcd (clockwise or counterclockwise), and the distance to four edges is at least mindist. Guaranteed P, A, B, C, D coplanar bool Insidewithmindistance (const POINT_3D & P, const POINT_3D & A, const POINT_3D & B, const POINT_3D & C, const POINT_3D & D, double mindist) {if (!    Pointintri (P, A, B, C)) {return false; } if (!    Pointintri (P, C, D, A)) {return false;    } if (Distancetoline (P, A, B) < Mindist) {return false;    } if (Distancetoline (P, B, C) < Mindist) {return false;    } if (Distancetoline (P, C, D) < Mindist) {return false;    } if (Distancetoline (P, D, A) < Mindist) {return false; } return true; /************* convex hull related problems *******************///plus interference double rand01 () {return rand ()/(double) Rand_max;} Double Randeps () {return (RAND01 ()-0.5) * EPS;} Point_3d add_noise (const POINT_3D & P) {return point_3d (p.x + randeps (), P.y + randeps (), P.z + randeps ());}    struct face{int v[3]; Face (int A, int b, int c) {v[0] = A;        V[1] = b;    V[2] = c; } vector_3d Normal (const vector<point_3d> & P) const {return cross (p[v[1])-p[v[0] [p[v[2]]-p[    V[0]]);  }//F can see p[i] int cansee (const vector<point_3d> & P, int i) const {return Dot (P[i]-p[v[0]],    Normal (P)) > 0; }};//incremental method for three-dimensional convex hull//NOTE: No special cases (such as four-point coplanar) are considered.    In practice, vector<face> ch3d (const vector<point_3d> & P) {int n = p.size () when the input point is slightly disturbed before the call;    vector<vector<int> > Vis (n);    for (int i = 0; i < n; i++) {vis[i].resize (n);    } vector<face> cur; Cur.push_back (Face (0, 1, 2));    Due to disturbances, the first three points are not collinear cur.push_back (Face (2, 1, 0));        for (int i = 3; i < n; i++) {vector<face> next;            Calculates the visibility of the "left" of each edge for (int j = 0; J < Cur.size (), j + +) {face & f = cur[j];            int res = F.cansee (P, i);      if (!res) {next.push_back (f);      } for (int k = 0; k < 3; k++) {vis[f.v[k]][f.v[(k + 1)% 3]] = res;                }} for (int j = 0; J < Cur.size (); j + +) for (int k = 0; k < 3; k++) {                int a = Cur[j].v[k], B = cur[j].v[(k + 1)% 3]; if (vis[a][b] = Vis[b][a] && vis[a][b])//(A, B) is the dividing line, the left side of P[i] is visible {next.push_back (                Face (A, B, i));    }} cur = next; } return cur;}    struct convexpolyhedron{int n;    Vector<point_3d> P, P2;    Vector<face> faces;        BOOL Read () {if (scanf ("%d", &n)! = 1) {return false;        } p.resize (n);        P2.resize (n);            for (int i = 0; i < n; i++) {P[i] = read_point_3d ();        P2[i] = add_noise (P[i]);        } faces = Ch3d (P2);    return true; }//three-dimensional convex hull center of gravity point_3d centroid () {point_3d C = p[0];       Double Totv = 0;        Point_3d tot (0, 0, 0); for (int i = 0; i < faces.size (); i++) {point_3d P1 = p[faces[i].v[0]], p2 = p[faces[i].v[1]], p3 = P            [Faces[i].v[2]];            Double v =-volume6 (P1, p2, p3, C);            Totv + = V;        tot = tot + centroid (p1, p2, p3, C) * v;    } return TOT/TOTV;        }//Bump center of gravity to surface nearest distance double mindist (point_3d C) {double ans = 1e30; for (int i = 0; i < faces.size (); i++) {point_3d P1 = p[faces[i].v[0]], p2 = p[faces[i].v[1]], p3 = P            [Faces[i].v[2]];        ans = min (ans, fabs (-volume6 (P1, p2, p3, C)/Area2 (P1, p2, p3)));    } return ans; }};int Main () {return 0;}

Three-dimensional computational geometry template