Three-dimensional convex hull + Point-to-plane distance + known 3-point plane equation

/* ===================================================== ===============* \ | 3D convex hull | call: build convex hull = construct (); \ * =================================================== ===============*/# define TPN 1010 struct tpoint {Double X, y, Z; void get () {scanf ("% lf", & 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;} // dot product}; struct FAC {int A, B, C; // bool OK of the three vertices on a convex hull surface; // whether the surface is in the final convex hull}; struct t3dhull {int N; // The initial point tpoint ply [TPN]; // The initial point int trianglecnt; // Number of triangles on the convex hull FAC tri [TPN]; // Number of triangles in the convex hull int vis [TPN] [TPN]; // void add () on which the point I-to-point j belongs () {ply [n ++]. get ();} double dist (tpoint A) {return SQRT (. x *. X +. y *. Y +. z *. Z);} // double area (tpoint A, tpoint B, tpoint c) {return dist (B-a) * (c-));} // Triangle Area x 2 double volume (tpoint A, tpoint B, tpoint C, tpoint d) {return (B-a) * (c-) ^ (D-A);} // The directed volume of the triangle x 6 double ptoplane (tpoint & P, FAC & F) {// positive: point in the same direction as face 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]; // The FAC add of the surface with the co-edge (AB) of the current surface (CNT); If (TRI [f]. OK) {If (pTop Lane (ply [p], Tri [f])> EPS) DFS (p, f); // if you can see this F, then, continue to explore the three sides of F in depth, so as to update the new convex hull surface else // otherwise, because point P only sees the CNT surface, not the f surface, then, points A and B form a triangle. {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) {// maintain the convex hull. If point P is outside the convex hull, update the convex hull tri [CNT]. OK = 0; // The current surface needs to be deleted because it places the opposite side in the larger convex hull // (first B, then a), so that in deal () to determine the surface of the co-edge (AB) with the current surface (CNT. That is, to judge the three sides adjacent to the Head Face (CNT) (they are opposite to the common edge of the current surface, such as the normal of the center (1) outward (that is, counterclockwise) surface 130 and 312, they have 13 sides, but one direction is 13 and the other direction is 31) 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) {// determine whether the two sides are the same. 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]) <EPS & FABS (volume (a, B, c, ply [TRI [e]. b]) <EPS & FABS (volume (a, B, c, ply [TRI [e]. c]) <EPS ;} Void construct () {// build a convex packet int I, j; trianglecnt = 0; If (n <4) return; bool TMP = true; for (I = 1; I <n; I ++) // The first two vertices do not share {If (Dist (ply [0]-ply [I])> EPS) {swap (ply [1], ply [I]); TMP = false; break;} If (TMP) return; TMP = true; for (I = 2; I <n; I ++) {// The first three points are not collocated if (Dist (ply [0]-ply [1]) * (ply [1]-ply [I])> EPS) {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])> EPS) {swap (ply [3], ply [I]); TMP = false; break;} If (TMP) return; FAC add; for (I = 0; I <4; I ++) {// construct the initial triangle (four vertices are ply [0], ply [1], ply [2], ply [3]) 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); // counter-clockwise, that is, the forward and outward vector of the normal vector, so that new points can be seen. Vis [add. a] [add. b] = vis [add. b] [add. c] = vis [add. c] [add. a] = trianglecnt; // reverse directed edge storage tri [trianglecnt ++] = add ;}for (I = 4; I <n; I ++) {// build and update a convex hull for (j = 0; j <trianglecnt; j ++) {// determine whether each vertex is inside or outside the current 3-dimensional convex hull (I indicates the current vertex, J indicates the current surface) if (TRI [J]. OK & (ptoplane (ply [I], Tri [J])> EPS) {// judge the current convex surface, check whether you can see the DFS (I, j); break; // click it to view the current surface. Update the convex surface (recursive, more than updating one surface ). After the current vertex is updated, the break jumps out of the loop.} int CNT = trianglecnt; // There are some Tri [I]. OK = 0. They are created but need to be deleted later because they are in a larger convex hull. Therefore, the following lines of code only store the outermost convex hull trianglecnt = 0; for (I = 0; I <CNT; I ++) {If (TRI [I]. OK) tri [trianglecnt ++] = tri [I] ;}} double area () {// Surface 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;} double volume () {// volume tpoint P (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);} int facetri () {return trianglecnt;} // number of surface triangles int facepolygon () {// number of surface polygon int ans = 0, I, J, K; for (I = 0; I <trianglecnt; I ++) {for (j = 0, k = 1; j <I; j ++) {If (same (I, j) {k = 0; break ;}} ans + = K ;} return ans ;}};

// Obtain the coordinate of three vertices. Obtain the plane AX + by + cz + D = 0. Void get_panel (tpoint P1, tpoint P2, tpoint P3, double & A, double & B, double & C, double & D) {A = (p2.y-p1.y) * (p3.z-p1.z)-(p2.z-p1.z) * (p3.y-p1.y )); B = (p2.z-p1.z) * (p3.x-p1.x)-(p2.x-p1.x) * (p3.z-p1.z); C = (p2.x-p1.x) * (p3.y-p1.y)-(p2.y-p1.y) * (p3.x-p1.x )); D = (0-(A * p1.x + B * p1.y + C * p1.z);} // distance from a point to a plane to a double dis_pt2panel (tpoint PT, double A, double B, double C, double D) {return f_abs (A * PT. X + B * PT. Y + C * PT. Z + d)/SQRT (A * A + B * B + C * C );}