diff --git a/src/mesh_fix/CMakeLists.txt b/src/mesh_fix/CMakeLists.txt index ecc43b415..d0a195131 100644 --- a/src/mesh_fix/CMakeLists.txt +++ b/src/mesh_fix/CMakeLists.txt @@ -31,6 +31,4 @@ target_include_directories(meshfix PUBLIC include/Kernel ) -target_compile_options(meshfix PRIVATE -fpermissive) - ################################################################################ diff --git a/src/mesh_fix/include/Kernel/basics.h b/src/mesh_fix/include/Kernel/basics.h index 406076e00..809f81b39 100644 --- a/src/mesh_fix/include/Kernel/basics.h +++ b/src/mesh_fix/include/Kernel/basics.h @@ -116,12 +116,6 @@ typedef unsigned short UINT16; typedef signed short INT16; #endif -#ifdef IS64BITPLATFORM -typedef long int j_voidint; -#else -typedef int j_voidint; -#endif - #define UBYTE_MAX 255 #ifndef UINT16_MAX @@ -138,18 +132,169 @@ typedef int j_voidint; #define MAX(a,b)(((a)>(b))?(a):(b)) #endif - //////// Swaps two pointers. /////////////////////////////// inline void p_swap(void **a, void **b) {void *t = *a; *a = *b; *b = t;} ///////////////////////////////////////////////////////////////////////////////////////////// + +template +class Primitive { +protected: + T value; + +public: + + // we must type cast to child to so + // a += 3 += 5 ... and etc.. work the same way + // as on primitives + Child &childRef(){ + return *((Child*)this); + } + + // you can overload to give a default value if you want + Primitive(){} + explicit Primitive(T v):value(v){} + + T get(){ + return value; + } + + #define OP(op) Child &operator op(Child const &v){\ + value op v.value; \ + return childRef(); \ + }\ + Child &operator op(T const &v){\ + value op v; \ + return childRef(); \ + } + + // all with equals + OP(+=) + OP(-=) + OP(*=) + OP(/=) + OP(<<=) + OP(>>=) + OP(|=) + OP(^=) + OP(&=) + OP(%=) + + #undef OP + + #define OP(p) Child operator p(Child const &v){\ + Child other = childRef();\ + other p ## = v;\ + return other;\ + }\ + Child operator p(T const &v){\ + Child other = childRef();\ + other p ## = v;\ + return other;\ + } + + OP(+) + OP(-) + OP(*) + OP(/) + OP(<<) + OP(>>) + OP(|) + OP(^) + OP(&) + OP(%) + + #undef OP + + + #define OP(p) bool operator p(Child const &v){\ + return value p v.value;\ + }\ + bool operator p(T const &v){\ + return value p v;\ + } + + OP(&&) + OP(||) + OP(<) + OP(<=) + OP(>) + OP(>=) + OP(==) + OP(!=) + + #undef OP + + Child operator +(){return Child(value);} + Child operator -(){return Child(-value);} + Child &operator ++(){++value; return childRef();} + Child operator ++(int){ + Child ret(value); + ++value; + return childRef(); + } + Child operator --(int){ + Child ret(value); + --value; + return childRef(); + } + + bool operator!(){return !value;} + Child operator~(){return Child(~value);} + +}; + + class Data { public: virtual ~Data() = default; }; +class intWrapper: public Data { +private: + int val; + public: + intWrapper(int val = 0) : + val(val) { + } + operator int &() { + return val; + } + int* operator &() { + return &val; + } + operator int() const { + return val; + } + operator int*() { + return &val; + } +}; + +class doubleWrapper: public Data, public Primitive { +public: + doubleWrapper(double val = 0) { + this->value = val; + } + operator double &() { + return value; + } + double* operator &() { + return &value; + } + operator double() const { + return value; + } + operator double*() { + return &value; + } +}; + +inline int to_int(Data *d) { + return static_cast(d)->operator int(); +} } //namespace T_MESH diff --git a/src/mesh_fix/include/Kernel/coordinates.h b/src/mesh_fix/include/Kernel/coordinates.h index 0d42228df..d58d279c4 100644 --- a/src/mesh_fix/include/Kernel/coordinates.h +++ b/src/mesh_fix/include/Kernel/coordinates.h @@ -149,10 +149,10 @@ PM_Rational ceil(const PM_Rational& a); PM_Rational floor(const PM_Rational& a); PM_Rational round(const PM_Rational& a); -/**************** I/O operators ****************/ - -inline std::ostream & operator<<(std::ostream &o, const PM_Rational& c) { return o << c.toRational(); } -inline std::istream & operator>>(std::istream &i, PM_Rational& c) { EXACT_NT a; i >> a; c.setFromRational(a); return i; } +/**************** I/O operators ****************/ + +inline std::ostream & operator<<(std::ostream &o, const PM_Rational& c) { return o << c.toRational(); } +inline std::istream & operator>>(std::istream &i, PM_Rational& c) { EXACT_NT a; i >> a; c.setFromRational(a); return i; } #define TMESH_TO_DOUBLE(x) ((x).toDouble()) #define TMESH_TO_FLOAT(x) ((x).toFloat()) diff --git a/src/mesh_fix/include/Kernel/heap.h b/src/mesh_fix/include/Kernel/heap.h index bfe7164e2..4768089bf 100644 --- a/src/mesh_fix/include/Kernel/heap.h +++ b/src/mesh_fix/include/Kernel/heap.h @@ -54,7 +54,7 @@ namespace T_MESH class abstractHeap { protected: - void **heap; //!< Heap data is stored here + Data **heap; //!< Heap data is stored here int numels; //!< Current number of elements int maxels; //!< Maximum number of elements int *positions; //!< Optional pointer to an array of positions @@ -67,7 +67,7 @@ class abstractHeap //! This function must be implemented in the extended class. //! The return value must be <0 if a0 if a>b or 0 if a=b. - virtual int compare(const void *a, const void *b) = 0; + virtual int compare(const Data *a, const Data *b) = 0; public : @@ -81,11 +81,11 @@ class abstractHeap //! returned, otherwise the index position of the newly inserted element is //! returned. - int insert(void *e); + int insert(Data *e); int isEmpty() const {return (numels == 0);} //!< Returns TRUE if the heap is empty - void *getHead() const {return heap[1];} //!< Returns the first element of the heap - void *removeHead(); //!< Removes and returns the first element after rearranging the heap + Data *getHead() const {return heap[1];} //!< Returns the first element of the heap + Data *removeHead(); //!< Removes and returns the first element after rearranging the heap void flush() {numels=0;} //!< Removes all the elements }; diff --git a/src/mesh_fix/include/Kernel/point.h b/src/mesh_fix/include/Kernel/point.h index 049b6b72d..c35dbcb38 100644 --- a/src/mesh_fix/include/Kernel/point.h +++ b/src/mesh_fix/include/Kernel/point.h @@ -72,7 +72,7 @@ class Point : public Data { public : coord x,y,z; //!< Coordinates - void *info; //!< Further information + Data *info; //!< Further information //! Creates a new point with coordinates (0,0,0). Point() {x = y = z = 0; info = NULL;} diff --git a/src/mesh_fix/include/TMesh/edge.h b/src/mesh_fix/include/TMesh/edge.h index dba8b6316..89b8bf5ed 100644 --- a/src/mesh_fix/include/TMesh/edge.h +++ b/src/mesh_fix/include/TMesh/edge.h @@ -55,7 +55,7 @@ class Edge : public Data Vertex *v1,*v2; //!< End-points class Triangle *t1,*t2; //!< Incident triangles unsigned char mask; //!< bit-mask for marking purposes - void *info; //!< Further information + Data *info; //!< Further information Edge(); //!< AMF_ADD 1.1-2 > diff --git a/src/mesh_fix/include/TMesh/marchIntersections.h b/src/mesh_fix/include/TMesh/marchIntersections.h index b1c822145..37a575fb7 100644 --- a/src/mesh_fix/include/TMesh/marchIntersections.h +++ b/src/mesh_fix/include/TMesh/marchIntersections.h @@ -77,7 +77,7 @@ class mc_cell : public Data } void polygonize(Basic_TMesh *tin); - static int compare(const void *e1, const void *e2); + static int compare(const Data *e1, const Data *e2); void merge(mc_cell *m); diff --git a/src/mesh_fix/include/TMesh/triangle.h b/src/mesh_fix/include/TMesh/triangle.h index 710d5e3ef..832679fee 100644 --- a/src/mesh_fix/include/TMesh/triangle.h +++ b/src/mesh_fix/include/TMesh/triangle.h @@ -51,7 +51,7 @@ class Triangle : public Data public : Edge *e1, *e2, *e3; //!< Edges of the triangle - void *info; //!< Further information + Data* info; //!< Further information unsigned char mask; //!< bit-mask for marking purposes Triangle(); diff --git a/src/mesh_fix/src/Algorithms/checkAndRepair.cpp b/src/mesh_fix/src/Algorithms/checkAndRepair.cpp index d524bb071..263046a1d 100644 --- a/src/mesh_fix/src/Algorithms/checkAndRepair.cpp +++ b/src/mesh_fix/src/Algorithms/checkAndRepair.cpp @@ -282,7 +282,7 @@ int Basic_TMesh::mergeCoincidentEdges() if (e->v2->info != e->v2) e->v2 = (Vertex *)e->v2->info; e->v1->e0 = e->v2->e0 = e; } - int rv = removeVertices(); + removeVertices(); // At this point the mesh should no longer have duplicated vertices, but may have duplicated edges E.sort(&vtxEdgeCompare); @@ -296,7 +296,6 @@ int Basic_TMesh::mergeCoincidentEdges() { Triangle *t1 = e->getBoundaryTriangle(); Edge *f = ((Edge *)e->info); - Triangle *t2 = f->getBoundaryTriangle(); t1->replaceEdge(e, f); ((f->t1 == NULL) ? (f->t1) : (f->t2)) = t1; e->v1 = e->v2 = NULL; @@ -355,14 +354,14 @@ bool Basic_TMesh::rebuildConnectivity(bool fixconnectivity) //!< AMF_CHANGE 1.1> Triangle *t; ExtVertex **var = new ExtVertex *[V.numels()]; int i=0; - FOREACHVERTEX(v, n) { v->e0 = NULL; var[i] = new ExtVertex(v); v->info = (void *)(intptr_t)i; i++; } + FOREACHVERTEX(v, n) { v->e0 = NULL; var[i] = new ExtVertex(v); v->info = new intWrapper(i); i++; } int nt = T.numels(); int *triangles = new int[nt*3]; i = 0; FOREACHTRIANGLE(t, n) { - triangles[i * 3] = (j_voidint)t->v1()->info; - triangles[i*3+1] = (j_voidint)t->v2()->info; - triangles[i*3+2] = (j_voidint)t->v3()->info; + triangles[i * 3] = ((intWrapper*)t->v1()->info)->operator int(); + triangles[i*3+1] = ((intWrapper*)t->v2()->info)->operator int(); + triangles[i*3+2] = ((intWrapper*)t->v3()->info)->operator int(); i++; } T.freeNodes(); @@ -450,7 +449,7 @@ int multiSplitEdge(Basic_TMesh *tin, Edge *e) while (splitVertices.numels()) { coord ad, mind = DBL_MAX; - Vertex *gv; + Vertex *gv = nullptr; FOREACHVVVERTEX((&splitVertices), v, n) if ((ad = v->squaredDistance(e->v2)) < mind) { gv = v; mind = ad; } splitVertices.removeNode(gv); tin->splitEdge(e, gv); diff --git a/src/mesh_fix/src/Algorithms/holeFilling.cpp b/src/mesh_fix/src/Algorithms/holeFilling.cpp index 793d45c95..c19557ec6 100644 --- a/src/mesh_fix/src/Algorithms/holeFilling.cpp +++ b/src/mesh_fix/src/Algorithms/holeFilling.cpp @@ -580,7 +580,7 @@ int Basic_TMesh::refineSelectedHolePatches(Triangle *t0) Edge *e, *f; Vertex *v; List *ve, toswap, reg, all_edges, interior_edges, boundary_edges, boundary_vertices, interior_vertices; - coord sigma, l, sv1, sv2, sv3, dv1, dv2, dv3; + doubleWrapper sigma, l, sv1, sv2, sv3, dv1, dv2, dv3; int swaps, totits, nee, ntb, nnt=-1, pnnt, gits=0; const double alpha = sqrt(2.0); Point vc; @@ -630,7 +630,7 @@ int Basic_TMesh::refineSelectedHolePatches(Triangle *t0) { ve = v->VE(); sigma=0; nee=0; FOREACHVEEDGE(ve, e, m) if (!IS_BIT(e, 5)) {nee++; sigma += e->length();} - sigma /= nee; v->info = new coord(sigma); + sigma /= double(nee); v->info = new doubleWrapper(sigma); delete(ve); } @@ -644,9 +644,9 @@ int Basic_TMesh::refineSelectedHolePatches(Triangle *t0) FOREACHVTTRIANGLE((®), t, n) { vc = t->getCenter(); - sv1 = (*(coord *)t->v1()->info); - sv2 = (*(coord *)t->v2()->info); - sv3 = (*(coord *)t->v3()->info); + sv1 = (*(doubleWrapper *)t->v1()->info); + sv2 = (*(doubleWrapper *)t->v2()->info); + sv3 = (*(doubleWrapper *)t->v3()->info); sigma = (sv1+sv2+sv3)/3.0; dv1 = alpha*(t->v1()->distance(&vc)); dv2 = alpha*(t->v2()->distance(&vc)); @@ -658,7 +658,7 @@ int Basic_TMesh::refineSelectedHolePatches(Triangle *t0) nnt += (T.numels()-ntb); if (T.numels() == ntb+2) { - v->info = new coord(sigma); + v->info = new doubleWrapper(sigma); interior_vertices.appendHead(v); interior_edges.appendHead(v->e0); interior_edges.appendHead(v->e0->leftTriangle(v)->prevEdge(v->e0)); @@ -700,8 +700,8 @@ int Basic_TMesh::refineSelectedHolePatches(Triangle *t0) } while (nnt && gits<10); //FOREACHVEEDGE((&boundary_edges), e, n) UNMARK_BIT(e, 6); - FOREACHVVVERTEX((&boundary_vertices), v, n) { delete((coord *)v->info); v->info = NULL; MARK_BIT(v, 5);} - FOREACHVVVERTEX((&interior_vertices), v, n) { delete((coord *)v->info); v->info = NULL; MARK_BIT(v, 6);} + FOREACHVVVERTEX((&boundary_vertices), v, n) { delete((Data *)v->info); v->info = NULL; MARK_BIT(v, 5);} + FOREACHVVVERTEX((&interior_vertices), v, n) { delete((Data *)v->info); v->info = NULL; MARK_BIT(v, 6);} if (gits>=10) {TMesh::warning("Fill holes: Refinement stage failed to converge. Breaking.\n"); return 1;} diff --git a/src/mesh_fix/src/Algorithms/marchIntersections.cpp b/src/mesh_fix/src/Algorithms/marchIntersections.cpp index 7d9a83351..8da0f440c 100644 --- a/src/mesh_fix/src/Algorithms/marchIntersections.cpp +++ b/src/mesh_fix/src/Algorithms/marchIntersections.cpp @@ -134,7 +134,7 @@ UBYTE mc_cell::lookdown() return i; } -int mc_cell::compare(const void *e1, const void *e2) +int mc_cell::compare(const Data *e1, const Data *e2) { mc_cell *a = (mc_cell *)e1; mc_cell *b = (mc_cell *)e2; @@ -680,7 +680,7 @@ void mc_cell::merge(mc_cell *m) List *mc_grid::createCells() { - int i,j,k,numcells=numrays+1; + int i,j,k=numrays+1; mc_ints *m; Node *n; List *ac = new List; diff --git a/src/mesh_fix/src/Kernel/heap.cpp b/src/mesh_fix/src/Kernel/heap.cpp index c488251a6..59fc5a418 100644 --- a/src/mesh_fix/src/Kernel/heap.cpp +++ b/src/mesh_fix/src/Kernel/heap.cpp @@ -37,7 +37,7 @@ namespace T_MESH abstractHeap::abstractHeap(int size) { - heap = new void *[size+1]; + heap = new Data *[size+1]; numels = 0; maxels = size; positions = NULL; @@ -52,9 +52,9 @@ int abstractHeap::upheap(int k) { if (k < 2) return k; - void *t = heap[k]; + Data *t = heap[k]; int fk = (k%2)?((k-1)/2):(k/2); - void *f = heap[fk]; + Data *f = heap[fk]; if (compare(t, f) <= 0) { @@ -62,8 +62,8 @@ int abstractHeap::upheap(int k) heap[fk] = t; if (positions != NULL) { - positions[(j_voidint)f] = k; - positions[(j_voidint)t] = fk; + positions[to_int(f)] = k; + positions[to_int(t)] = fk; } return upheap(fk); } @@ -74,21 +74,21 @@ int abstractHeap::downheap(int k) { int j; - void *t = heap[k]; + Data *t = heap[k]; int fk = (numels%2)?((numels-1)/2):(numels/2); if (k > fk) return k; j = k+k; if (j < numels && compare(heap[j], heap[j+1]) >= 0) j++; - void *f = heap[j]; + Data *f = heap[j]; if (compare(t, f) >= 0) { heap[k] = f; heap[j] = t; if (positions != NULL) { - positions[(j_voidint)f] = k; - positions[(j_voidint)t] = j; + positions[to_int(f)] = k; + positions[to_int(t)] = j; } return downheap(j); } @@ -96,23 +96,23 @@ int abstractHeap::downheap(int k) return k; } -int abstractHeap::insert(void *t) +int abstractHeap::insert(Data *t) { if (numels == maxels) return -1; heap[++numels] = t; - if (positions != NULL) positions[(j_voidint)t] = numels; + if (positions != NULL) positions[to_int(t)] = numels; return upheap(numels); } -void *abstractHeap::removeHead() +Data *abstractHeap::removeHead() { - void *t = heap[1]; - if (positions != NULL) positions[(j_voidint)t] = 0; + Data *t = heap[1]; + if (positions != NULL) positions[to_int(t)] = 0; heap[1] = heap[numels--]; if (numels) { - if (positions != NULL) positions[(j_voidint)heap[1]] = 1; + if (positions != NULL) positions[to_int(heap[1])] = 1; downheap(1); } diff --git a/src/mesh_fix/src/Kernel/orientation.c b/src/mesh_fix/src/Kernel/orientation.c index ab33dfd00..97a0591d2 100644 --- a/src/mesh_fix/src/Kernel/orientation.c +++ b/src/mesh_fix/src/Kernel/orientation.c @@ -28,673 +28,673 @@ * * ****************************************************************************/ -// This code is inspired on ideas first published in the following paper: -// Jonathan Richard Shewchuk. Adaptive Precision Floating-Point Arithmetic -// and Fast Robust Geometric Predicates, Discrete & Computational Geometry -// 18(3):305–363, October 1997. -// - -#include - -#ifdef SPECIFY_FP_PRECISION -#include -#endif - -/*****************************************************************************/ -/* */ -/* This section contains private macros and functions. */ -/* These are not part of the public interface. */ -/* */ -/*****************************************************************************/ - -#define FABS(a) (((a)>=0.0)?(a):(-(a))) -#define FTST(a,b,x,y) _bvr=x-a; y=b-_bvr -#define FTS(a,b,x,y) x=(double)(a+b); FTST(a,b,x,y) -#define TST(a,b,x,y) _bvr=(double)(x-a); _avr=x-_bvr; _brn=b-_bvr; _arn=a-_avr; y=_arn+_brn -#define TWS(a,b,x,y) x=(double)(a+b); TST(a,b,x,y) -#define TDT(a,b,x,y) _bvr=(double)(a-x); _avr=x+_bvr; _brn=_bvr-b; _arn=a-_avr; y=_arn+_brn -#define TWD(a,b,x,y) x=(double)(a-b); TDT(a,b,x,y) -#define SPLT(a,ahi,alo) c=(double)(_spl*a); abig=(double)(c-a); ahi=c-abig; alo=a-ahi -#define TPT(a,b,x,y) SPLT(a,ahi,alo); SPLT(b,bhi,blo); err1=x-(ahi*bhi); err2=err1-(alo*bhi); err3=err2-(ahi*blo); y=(alo*blo)-err3 -#define TWP(a,b,x,y) x=(double)(a*b); TPT(a,b,x,y) -#define TPP(a,b,bhi,blo,x,y) x=(double)(a*b); SPLT(a,ahi,alo); err1=x-(ahi*bhi); err2=err1-(alo*bhi); err3=err2-(ahi*blo); y=(alo*blo)-err3 -#define TOD(a1,a0,b,x2,x1,x0) TWD(a0,b,_i,x0); TWS(a1,_i,x2,x1) -#define TTD(a1,a0,b1,b0,x3,x2,x1,x0) TOD(a1,a0,b0,_j,_0,x0); TOD(_j,_0,b1,x3,x2,x1) -#define TOP(a1,a0,b,x3,x2,x1,x0) SPLT(b,bhi,blo); TPP(a0,b,bhi,blo,_i,x0); TPP(a1,b,bhi,blo,_j,_0); TWS(_i,_0,_k,x1); FTS(_j,_k,x3,x2) - -double _spl, _eps, _reb, _ccwebA, _ccwebB, _ccwebC, _o3ebA, _o3ebB, _o3ebC; -double _iccebA, _iccebB, _iccebC, _ispebA, _ispebB, _ispebC; - -int _fesze(int elen, double *e, int flen, double *f, double *h) -{ - double Q, Qnew, hh, _bvr, _avr, _brn, _arn, enow, fnow; - int eindex, findex, hindex; - - enow = e[0]; - fnow = f[0]; - eindex = findex = 0; - if ((fnow > enow) == (fnow > -enow)) { - Q = enow; - enow = e[++eindex]; - } else { - Q = fnow; - fnow = f[++findex]; - } - hindex = 0; - if ((eindex < elen) && (findex < flen)) { - if ((fnow > enow) == (fnow > -enow)) { - FTS(enow, Q, Qnew, hh); - enow = e[++eindex]; - } else { - FTS(fnow, Q, Qnew, hh); - fnow = f[++findex]; - } - Q = Qnew; - if (hh != 0.0) { - h[hindex++] = hh; - } - while ((eindex < elen) && (findex < flen)) { - if ((fnow > enow) == (fnow > -enow)) { - TWS(Q, enow, Qnew, hh); - enow = e[++eindex]; - } else { - TWS(Q, fnow, Qnew, hh); - fnow = f[++findex]; - } - Q = Qnew; - if (hh != 0.0) { - h[hindex++] = hh; - } - } - } - while (eindex < elen) { - TWS(Q, enow, Qnew, hh); - enow = e[++eindex]; - Q = Qnew; - if (hh != 0.0) { - h[hindex++] = hh; - } - } - while (findex < flen) { - TWS(Q, fnow, Qnew, hh); - fnow = f[++findex]; - Q = Qnew; - if (hh != 0.0) { - h[hindex++] = hh; - } - } - if ((Q != 0.0) || (hindex == 0)) { - h[hindex++] = Q; - } - return hindex; -} - -int _seze(int elen, double *e, double b, double *h) -{ - double Q, sum, hh, product1, product0, enow, _bvr, _avr, _brn, _arn, c; - double abig, ahi, alo, bhi, blo, err1, err2, err3; - int eindex, hindex; - - SPLT(b, bhi, blo); - TPP(e[0], b, bhi, blo, Q, hh); - hindex = 0; - if (hh != 0) { - h[hindex++] = hh; - } - for (eindex = 1; eindex < elen; eindex++) { - enow = e[eindex]; - TPP(enow, b, bhi, blo, product1, product0); - TWS(Q, product0, sum, hh); - if (hh != 0) { - h[hindex++] = hh; - } - FTS(product1, sum, Q, hh); - if (hh != 0) { - h[hindex++] = hh; - } - } - if ((Q != 0.0) || (hindex == 0)) { - h[hindex++] = Q; - } - return hindex; -} - -double _estm(int elen, double *e) -{ - int eindex; - double Q = e[0]; - for (eindex = 1; eindex < elen; eindex++) Q += e[eindex]; - return Q; -} - - -double _adaptive2dorientation(double *pa, double *pb, double *pc, double detsum) -{ - double acx, acy, bcx, bcy,acxtail, acytail, bcxtail, bcytail, detleft, detright; - double detlefttail, detrighttail, det, errbound, B[4], C1[8], C2[12], D[16]; - double B3, u[4], u3, s1, t1, s0, t0, _bvr, _avr, _brn, _arn, c; - double abig, ahi, alo, bhi, blo, err1, err2, err3, _i, _j, _0; - int C1length, C2length, Dlength; - - acx = (double) (pa[0] - pc[0]); - bcx = (double) (pb[0] - pc[0]); - acy = (double) (pa[1] - pc[1]); - bcy = (double) (pb[1] - pc[1]); - - TWP(acx, bcy, detleft, detlefttail); - TWP(acy, bcx, detright, detrighttail); - - TTD(detleft, detlefttail, detright, detrighttail, - B3, B[2], B[1], B[0]); - B[3] = B3; - - det = _estm(4, B); - errbound = _ccwebB * detsum; - if ((det >= errbound) || (-det >= errbound)) { - return det; - } - - TDT(pa[0], pc[0], acx, acxtail); - TDT(pb[0], pc[0], bcx, bcxtail); - TDT(pa[1], pc[1], acy, acytail); - TDT(pb[1], pc[1], bcy, bcytail); - - if ((acxtail == 0.0) && (acytail == 0.0) - && (bcxtail == 0.0) && (bcytail == 0.0)) { - return det; - } - - errbound = _ccwebC * detsum + _reb * FABS(det); - det += (acx * bcytail + bcy * acxtail) - - (acy * bcxtail + bcx * acytail); - if ((det >= errbound) || (-det >= errbound)) { - return det; - } - - TWP(acxtail, bcy, s1, s0); - TWP(acytail, bcx, t1, t0); - TTD(s1, s0, t1, t0, u3, u[2], u[1], u[0]); - u[3] = u3; - C1length = _fesze(4, B, 4, u, C1); - - TWP(acx, bcytail, s1, s0); - TWP(acy, bcxtail, t1, t0); - TTD(s1, s0, t1, t0, u3, u[2], u[1], u[0]); - u[3] = u3; - C2length = _fesze(C1length, C1, 4, u, C2); - - TWP(acxtail, bcytail, s1, s0); - TWP(acytail, bcxtail, t1, t0); - TTD(s1, s0, t1, t0, u3, u[2], u[1], u[0]); - u[3] = u3; - Dlength = _fesze(C2length, C2, 4, u, D); - - return(D[Dlength - 1]); -} - -double _adaptive3dorientation(double *pa, double *pb, double *pc, double *pd, double permanent) -{ - double adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz, det, errbound; - double bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1; - double bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0, bc[4], ca[4], ab[4]; - double bc3, ca3, ab3, adet[8], bdet[8], cdet[8]; - double abdet[16], *finnow, *finother, *finswap, fin1[192], fin2[192]; - double adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail, adztail, bdztail, cdztail; - double at_blarge, at_clarge, bt_clarge, bt_alarge, ct_alarge, ct_blarge; - double at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4]; - double bdxt_cdy1, cdxt_bdy1, cdxt_ady1, adxt_cdy1, adxt_bdy1, bdxt_ady1; - double bdxt_cdy0, cdxt_bdy0, cdxt_ady0, adxt_cdy0, adxt_bdy0, bdxt_ady0; - double bdyt_cdx1, cdyt_bdx1, cdyt_adx1, adyt_cdx1, adyt_bdx1, bdyt_adx1; - double bdyt_cdx0, cdyt_bdx0, cdyt_adx0, adyt_cdx0, adyt_bdx0, bdyt_adx0; - double bct[8], cat[8], abt[8], bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1; - double adxt_cdyt1, adxt_bdyt1, bdxt_adyt1, bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0; - double adxt_cdyt0, adxt_bdyt0, bdxt_adyt0, u[4], v[12], w[16], u3, negate; - double _bvr, _avr, _brn, _arn, c, abig, ahi, alo, bhi, blo; - double err1, err2, err3, _i, _j, _k, _0; - int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen, finlength; - int vlength, wlength, alen, blen, clen, ablen, bctlen, catlen, abtlen; - - adx = (double) (pa[0] - pd[0]); - bdx = (double) (pb[0] - pd[0]); - cdx = (double) (pc[0] - pd[0]); - ady = (double) (pa[1] - pd[1]); - bdy = (double) (pb[1] - pd[1]); - cdy = (double) (pc[1] - pd[1]); - adz = (double) (pa[2] - pd[2]); - bdz = (double) (pb[2] - pd[2]); - cdz = (double) (pc[2] - pd[2]); - - TWP(bdx, cdy, bdxcdy1, bdxcdy0); - TWP(cdx, bdy, cdxbdy1, cdxbdy0); - TTD(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]); - bc[3] = bc3; - alen = _seze(4, bc, adz, adet); - - TWP(cdx, ady, cdxady1, cdxady0); - TWP(adx, cdy, adxcdy1, adxcdy0); - TTD(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]); - ca[3] = ca3; - blen = _seze(4, ca, bdz, bdet); - - TWP(adx, bdy, adxbdy1, adxbdy0); - TWP(bdx, ady, bdxady1, bdxady0); - TTD(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]); - ab[3] = ab3; - clen = _seze(4, ab, cdz, cdet); - - ablen = _fesze(alen, adet, blen, bdet, abdet); - finlength = _fesze(ablen, abdet, clen, cdet, fin1); - - det = _estm(finlength, fin1); - errbound = _o3ebB * permanent; - if ((det >= errbound) || (-det >= errbound)) { - return det; - } - - TDT(pa[0], pd[0], adx, adxtail); - TDT(pb[0], pd[0], bdx, bdxtail); - TDT(pc[0], pd[0], cdx, cdxtail); - TDT(pa[1], pd[1], ady, adytail); - TDT(pb[1], pd[1], bdy, bdytail); - TDT(pc[1], pd[1], cdy, cdytail); - TDT(pa[2], pd[2], adz, adztail); - TDT(pb[2], pd[2], bdz, bdztail); - TDT(pc[2], pd[2], cdz, cdztail); - - if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) - && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0) - && (adztail == 0.0) && (bdztail == 0.0) && (cdztail == 0.0)) { - return det; - } - - errbound = _o3ebC * permanent + _reb * FABS(det); - det += (adz * ((bdx * cdytail + cdy * bdxtail) - - (bdy * cdxtail + cdx * bdytail)) - + adztail * (bdx * cdy - bdy * cdx)) - + (bdz * ((cdx * adytail + ady * cdxtail) - - (cdy * adxtail + adx * cdytail)) - + bdztail * (cdx * ady - cdy * adx)) - + (cdz * ((adx * bdytail + bdy * adxtail) - - (ady * bdxtail + bdx * adytail)) - + cdztail * (adx * bdy - ady * bdx)); - if ((det >= errbound) || (-det >= errbound)) { - return det; - } - - finnow = fin1; - finother = fin2; - - if (adxtail == 0.0) { - if (adytail == 0.0) { - at_b[0] = 0.0; - at_blen = 1; - at_c[0] = 0.0; - at_clen = 1; - } else { - negate = -adytail; - TWP(negate, bdx, at_blarge, at_b[0]); - at_b[1] = at_blarge; - at_blen = 2; - TWP(adytail, cdx, at_clarge, at_c[0]); - at_c[1] = at_clarge; - at_clen = 2; - } - } else { - if (adytail == 0.0) { - TWP(adxtail, bdy, at_blarge, at_b[0]); - at_b[1] = at_blarge; - at_blen = 2; - negate = -adxtail; - TWP(negate, cdy, at_clarge, at_c[0]); - at_c[1] = at_clarge; - at_clen = 2; - } else { - TWP(adxtail, bdy, adxt_bdy1, adxt_bdy0); - TWP(adytail, bdx, adyt_bdx1, adyt_bdx0); - TTD(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0, - at_blarge, at_b[2], at_b[1], at_b[0]); - at_b[3] = at_blarge; - at_blen = 4; - TWP(adytail, cdx, adyt_cdx1, adyt_cdx0); - TWP(adxtail, cdy, adxt_cdy1, adxt_cdy0); - TTD(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0, - at_clarge, at_c[2], at_c[1], at_c[0]); - at_c[3] = at_clarge; - at_clen = 4; - } - } - if (bdxtail == 0.0) { - if (bdytail == 0.0) { - bt_c[0] = 0.0; - bt_clen = 1; - bt_a[0] = 0.0; - bt_alen = 1; - } else { - negate = -bdytail; - TWP(negate, cdx, bt_clarge, bt_c[0]); - bt_c[1] = bt_clarge; - bt_clen = 2; - TWP(bdytail, adx, bt_alarge, bt_a[0]); - bt_a[1] = bt_alarge; - bt_alen = 2; - } - } else { - if (bdytail == 0.0) { - TWP(bdxtail, cdy, bt_clarge, bt_c[0]); - bt_c[1] = bt_clarge; - bt_clen = 2; - negate = -bdxtail; - TWP(negate, ady, bt_alarge, bt_a[0]); - bt_a[1] = bt_alarge; - bt_alen = 2; - } else { - TWP(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0); - TWP(bdytail, cdx, bdyt_cdx1, bdyt_cdx0); - TTD(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0, - bt_clarge, bt_c[2], bt_c[1], bt_c[0]); - bt_c[3] = bt_clarge; - bt_clen = 4; - TWP(bdytail, adx, bdyt_adx1, bdyt_adx0); - TWP(bdxtail, ady, bdxt_ady1, bdxt_ady0); - TTD(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0, - bt_alarge, bt_a[2], bt_a[1], bt_a[0]); - bt_a[3] = bt_alarge; - bt_alen = 4; - } - } - if (cdxtail == 0.0) { - if (cdytail == 0.0) { - ct_a[0] = 0.0; - ct_alen = 1; - ct_b[0] = 0.0; - ct_blen = 1; - } else { - negate = -cdytail; - TWP(negate, adx, ct_alarge, ct_a[0]); - ct_a[1] = ct_alarge; - ct_alen = 2; - TWP(cdytail, bdx, ct_blarge, ct_b[0]); - ct_b[1] = ct_blarge; - ct_blen = 2; - } - } else { - if (cdytail == 0.0) { - TWP(cdxtail, ady, ct_alarge, ct_a[0]); - ct_a[1] = ct_alarge; - ct_alen = 2; - negate = -cdxtail; - TWP(negate, bdy, ct_blarge, ct_b[0]); - ct_b[1] = ct_blarge; - ct_blen = 2; - } else { - TWP(cdxtail, ady, cdxt_ady1, cdxt_ady0); - TWP(cdytail, adx, cdyt_adx1, cdyt_adx0); - TTD(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0, - ct_alarge, ct_a[2], ct_a[1], ct_a[0]); - ct_a[3] = ct_alarge; - ct_alen = 4; - TWP(cdytail, bdx, cdyt_bdx1, cdyt_bdx0); - TWP(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0); - TTD(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0, - ct_blarge, ct_b[2], ct_b[1], ct_b[0]); - ct_b[3] = ct_blarge; - ct_blen = 4; - } - } - - bctlen = _fesze(bt_clen, bt_c, ct_blen, ct_b, bct); - wlength = _seze(bctlen, bct, adz, w); - finlength = _fesze(finlength, finnow, wlength, w, - finother); - finswap = finnow; finnow = finother; finother = finswap; - - catlen = _fesze(ct_alen, ct_a, at_clen, at_c, cat); - wlength = _seze(catlen, cat, bdz, w); - finlength = _fesze(finlength, finnow, wlength, w, - finother); - finswap = finnow; finnow = finother; finother = finswap; - - abtlen = _fesze(at_blen, at_b, bt_alen, bt_a, abt); - wlength = _seze(abtlen, abt, cdz, w); - finlength = _fesze(finlength, finnow, wlength, w, - finother); - finswap = finnow; finnow = finother; finother = finswap; - - if (adztail != 0.0) { - vlength = _seze(4, bc, adztail, v); - finlength = _fesze(finlength, finnow, vlength, v, - finother); - finswap = finnow; finnow = finother; finother = finswap; - } - if (bdztail != 0.0) { - vlength = _seze(4, ca, bdztail, v); - finlength = _fesze(finlength, finnow, vlength, v, - finother); - finswap = finnow; finnow = finother; finother = finswap; - } - if (cdztail != 0.0) { - vlength = _seze(4, ab, cdztail, v); - finlength = _fesze(finlength, finnow, vlength, v, - finother); - finswap = finnow; finnow = finother; finother = finswap; - } - - if (adxtail != 0.0) { - if (bdytail != 0.0) { - TWP(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0); - TOP(adxt_bdyt1, adxt_bdyt0, cdz, u3, u[2], u[1], u[0]); - u[3] = u3; - finlength = _fesze(finlength, finnow, 4, u, - finother); - finswap = finnow; finnow = finother; finother = finswap; - if (cdztail != 0.0) { - TOP(adxt_bdyt1, adxt_bdyt0, cdztail, u3, u[2], u[1], u[0]); - u[3] = u3; - finlength = _fesze(finlength, finnow, 4, u, - finother); - finswap = finnow; finnow = finother; finother = finswap; - } - } - if (cdytail != 0.0) { - negate = -adxtail; - TWP(negate, cdytail, adxt_cdyt1, adxt_cdyt0); - TOP(adxt_cdyt1, adxt_cdyt0, bdz, u3, u[2], u[1], u[0]); - u[3] = u3; - finlength = _fesze(finlength, finnow, 4, u, - finother); - finswap = finnow; finnow = finother; finother = finswap; - if (bdztail != 0.0) { - TOP(adxt_cdyt1, adxt_cdyt0, bdztail, u3, u[2], u[1], u[0]); - u[3] = u3; - finlength = _fesze(finlength, finnow, 4, u, - finother); - finswap = finnow; finnow = finother; finother = finswap; - } - } - } - if (bdxtail != 0.0) { - if (cdytail != 0.0) { - TWP(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0); - TOP(bdxt_cdyt1, bdxt_cdyt0, adz, u3, u[2], u[1], u[0]); - u[3] = u3; - finlength = _fesze(finlength, finnow, 4, u, - finother); - finswap = finnow; finnow = finother; finother = finswap; - if (adztail != 0.0) { - TOP(bdxt_cdyt1, bdxt_cdyt0, adztail, u3, u[2], u[1], u[0]); - u[3] = u3; - finlength = _fesze(finlength, finnow, 4, u, - finother); - finswap = finnow; finnow = finother; finother = finswap; - } - } - if (adytail != 0.0) { - negate = -bdxtail; - TWP(negate, adytail, bdxt_adyt1, bdxt_adyt0); - TOP(bdxt_adyt1, bdxt_adyt0, cdz, u3, u[2], u[1], u[0]); - u[3] = u3; - finlength = _fesze(finlength, finnow, 4, u, - finother); - finswap = finnow; finnow = finother; finother = finswap; - if (cdztail != 0.0) { - TOP(bdxt_adyt1, bdxt_adyt0, cdztail, u3, u[2], u[1], u[0]); - u[3] = u3; - finlength = _fesze(finlength, finnow, 4, u, - finother); - finswap = finnow; finnow = finother; finother = finswap; - } - } - } - if (cdxtail != 0.0) { - if (adytail != 0.0) { - TWP(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0); - TOP(cdxt_adyt1, cdxt_adyt0, bdz, u3, u[2], u[1], u[0]); - u[3] = u3; - finlength = _fesze(finlength, finnow, 4, u, - finother); - finswap = finnow; finnow = finother; finother = finswap; - if (bdztail != 0.0) { - TOP(cdxt_adyt1, cdxt_adyt0, bdztail, u3, u[2], u[1], u[0]); - u[3] = u3; - finlength = _fesze(finlength, finnow, 4, u, - finother); - finswap = finnow; finnow = finother; finother = finswap; - } - } - if (bdytail != 0.0) { - negate = -cdxtail; - TWP(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0); - TOP(cdxt_bdyt1, cdxt_bdyt0, adz, u3, u[2], u[1], u[0]); - u[3] = u3; - finlength = _fesze(finlength, finnow, 4, u, - finother); - finswap = finnow; finnow = finother; finother = finswap; - if (adztail != 0.0) { - TOP(cdxt_bdyt1, cdxt_bdyt0, adztail, u3, u[2], u[1], u[0]); - u[3] = u3; - finlength = _fesze(finlength, finnow, 4, u, - finother); - finswap = finnow; finnow = finother; finother = finswap; - } - } - } - - if (adztail != 0.0) { - wlength = _seze(bctlen, bct, adztail, w); - finlength = _fesze(finlength, finnow, wlength, w, - finother); - finswap = finnow; finnow = finother; finother = finswap; - } - if (bdztail != 0.0) { - wlength = _seze(catlen, cat, bdztail, w); - finlength = _fesze(finlength, finnow, wlength, w, - finother); - finswap = finnow; finnow = finother; finother = finswap; - } - if (cdztail != 0.0) { - wlength = _seze(abtlen, abt, cdztail, w); - finlength = _fesze(finlength, finnow, wlength, w, - finother); - finswap = finnow; finnow = finother; finother = finswap; - } - - return finnow[finlength - 1]; -} - - - - - -/*****************************************************************************/ -/* */ -/* PUBLIC FUNCTIONS */ -/* initPredicates() Sets the variables used for exact arithmetic. This */ -/* must be called once before using the other two functions. */ -/* orient2d() Computes the orientation of three 2D points. */ -/* orient3d() Computes the orientation of four 3D points. */ -/* */ -/*****************************************************************************/ - -void initPredicates() -{ - static char a_c=0; - double hf, ck, lc; - int e_o; - - if (a_c) return; else a_c = 1; - -#ifdef SPECIFY_FP_PRECISION - unsigned int old_cfp; - _controlfp_s(&old_cfp, _PC_53, MCW_PC); -#endif - - e_o = 1; - _eps = _spl = ck = 1.0; - hf = 0.5; - - do - { - lc=ck; - _eps *= hf; - if (e_o) _spl *= 2.0; - e_o = !e_o; - ck = 1.0 + _eps; - } while ((ck != 1.0) && (ck != lc)); - _spl += 1.0; - - _reb = (3.0 + 8.0 * _eps) * _eps; - _ccwebA = (3.0 + 16.0 * _eps) * _eps; - _ccwebB = (2.0 + 12.0 * _eps) * _eps; - _ccwebC = (9.0 + 64.0 * _eps) * _eps * _eps; - _o3ebA = (7.0 + 56.0 * _eps) * _eps; - _o3ebB = (3.0 + 28.0 * _eps) * _eps; - _o3ebC = (26.0 + 288.0 * _eps) * _eps * _eps; - _iccebA = (10.0 + 96.0 * _eps) * _eps; - _iccebB = (4.0 + 48.0 * _eps) * _eps; - _iccebC = (44.0 + 576.0 * _eps) * _eps * _eps; - _ispebA = (16.0 + 224.0 * _eps) * _eps; - _ispebB = (5.0 + 72.0 * _eps) * _eps; - _ispebC = (71.0 + 1408.0 * _eps) * _eps * _eps; - -#ifdef SPECIFY_FP_PRECISION - _controlfp_s(&old_cfp, _CW_DEFAULT, MCW_PC); -#endif -} - -double orient2d(double *pa, double *pb, double *pc) -{ - double dlf, drg, det, dsm, eb; - - dlf = (pa[0]-pc[0])*(pb[1]-pc[1]); - drg = (pa[1]-pc[1])*(pb[0]-pc[0]); - det = dlf - drg; - - if (dlf > 0.0) {if (drg <= 0.0) return det; else dsm = dlf + drg;} - else if (dlf < 0.0) {if (drg >= 0.0) return det; else dsm = -dlf - drg;} - else return det; - - eb = _ccwebA*dsm; - if ((det>=eb) || (-det>=eb)) return det; - - return _adaptive2dorientation(pa, pb, pc, dsm); -} - -double orient3d(double *pa, double *pb, double *pc, double *pd) -{ - double adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz, pm, eb; - double bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady, det; - - adx = pa[0]-pd[0]; bdx = pb[0]-pd[0]; cdx = pc[0]-pd[0]; - ady = pa[1]-pd[1]; bdy = pb[1]-pd[1]; cdy = pc[1]-pd[1]; - adz = pa[2]-pd[2]; bdz = pb[2]-pd[2]; cdz = pc[2]-pd[2]; - - bdxcdy = bdx*cdy; cdxbdy = cdx*bdy; - cdxady = cdx*ady; adxcdy = adx*cdy; - adxbdy = adx*bdy; bdxady = bdx*ady; - - det = adz*(bdxcdy-cdxbdy)+bdz*(cdxady-adxcdy)+cdz*(adxbdy-bdxady); - pm=(FABS(bdxcdy)+FABS(cdxbdy))*FABS(adz)+(FABS(cdxady)+FABS(adxcdy))*FABS(bdz)+(FABS(adxbdy)+FABS(bdxady))*FABS(cdz); - eb = _o3ebA*pm; - if ((det>eb) || (-det>eb)) return det; - return _adaptive3dorientation(pa, pb, pc, pd, pm); -} +// This code is inspired on ideas first published in the following paper: +// Jonathan Richard Shewchuk. Adaptive Precision Floating-Point Arithmetic +// and Fast Robust Geometric Predicates, Discrete & Computational Geometry +// 18(3):305–363, October 1997. +// + +#include + +#ifdef SPECIFY_FP_PRECISION +#include +#endif + +/*****************************************************************************/ +/* */ +/* This section contains private macros and functions. */ +/* These are not part of the public interface. */ +/* */ +/*****************************************************************************/ + +#define FABS(a) (((a)>=0.0)?(a):(-(a))) +#define FTST(a,b,x,y) _bvr=x-a; y=b-_bvr +#define FTS(a,b,x,y) x=(double)(a+b); FTST(a,b,x,y) +#define TST(a,b,x,y) _bvr=(double)(x-a); _avr=x-_bvr; _brn=b-_bvr; _arn=a-_avr; y=_arn+_brn +#define TWS(a,b,x,y) x=(double)(a+b); TST(a,b,x,y) +#define TDT(a,b,x,y) _bvr=(double)(a-x); _avr=x+_bvr; _brn=_bvr-b; _arn=a-_avr; y=_arn+_brn +#define TWD(a,b,x,y) x=(double)(a-b); TDT(a,b,x,y) +#define SPLT(a,ahi,alo) c=(double)(_spl*a); abig=(double)(c-a); ahi=c-abig; alo=a-ahi +#define TPT(a,b,x,y) SPLT(a,ahi,alo); SPLT(b,bhi,blo); err1=x-(ahi*bhi); err2=err1-(alo*bhi); err3=err2-(ahi*blo); y=(alo*blo)-err3 +#define TWP(a,b,x,y) x=(double)(a*b); TPT(a,b,x,y) +#define TPP(a,b,bhi,blo,x,y) x=(double)(a*b); SPLT(a,ahi,alo); err1=x-(ahi*bhi); err2=err1-(alo*bhi); err3=err2-(ahi*blo); y=(alo*blo)-err3 +#define TOD(a1,a0,b,x2,x1,x0) TWD(a0,b,_i,x0); TWS(a1,_i,x2,x1) +#define TTD(a1,a0,b1,b0,x3,x2,x1,x0) TOD(a1,a0,b0,_j,_0,x0); TOD(_j,_0,b1,x3,x2,x1) +#define TOP(a1,a0,b,x3,x2,x1,x0) SPLT(b,bhi,blo); TPP(a0,b,bhi,blo,_i,x0); TPP(a1,b,bhi,blo,_j,_0); TWS(_i,_0,_k,x1); FTS(_j,_k,x3,x2) + +double _spl, _eps, _reb, _ccwebA, _ccwebB, _ccwebC, _o3ebA, _o3ebB, _o3ebC; +double _iccebA, _iccebB, _iccebC, _ispebA, _ispebB, _ispebC; + +int _fesze(int elen, double *e, int flen, double *f, double *h) +{ + double Q, Qnew, hh, _bvr, _avr, _brn, _arn, enow, fnow; + int eindex, findex, hindex; + + enow = e[0]; + fnow = f[0]; + eindex = findex = 0; + if ((fnow > enow) == (fnow > -enow)) { + Q = enow; + enow = e[++eindex]; + } else { + Q = fnow; + fnow = f[++findex]; + } + hindex = 0; + if ((eindex < elen) && (findex < flen)) { + if ((fnow > enow) == (fnow > -enow)) { + FTS(enow, Q, Qnew, hh); + enow = e[++eindex]; + } else { + FTS(fnow, Q, Qnew, hh); + fnow = f[++findex]; + } + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + while ((eindex < elen) && (findex < flen)) { + if ((fnow > enow) == (fnow > -enow)) { + TWS(Q, enow, Qnew, hh); + enow = e[++eindex]; + } else { + TWS(Q, fnow, Qnew, hh); + fnow = f[++findex]; + } + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + } + while (eindex < elen) { + TWS(Q, enow, Qnew, hh); + enow = e[++eindex]; + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + while (findex < flen) { + TWS(Q, fnow, Qnew, hh); + fnow = f[++findex]; + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + if ((Q != 0.0) || (hindex == 0)) { + h[hindex++] = Q; + } + return hindex; +} + +int _seze(int elen, double *e, double b, double *h) +{ + double Q, sum, hh, product1, product0, enow, _bvr, _avr, _brn, _arn, c; + double abig, ahi, alo, bhi, blo, err1, err2, err3; + int eindex, hindex; + + SPLT(b, bhi, blo); + TPP(e[0], b, bhi, blo, Q, hh); + hindex = 0; + if (hh != 0) { + h[hindex++] = hh; + } + for (eindex = 1; eindex < elen; eindex++) { + enow = e[eindex]; + TPP(enow, b, bhi, blo, product1, product0); + TWS(Q, product0, sum, hh); + if (hh != 0) { + h[hindex++] = hh; + } + FTS(product1, sum, Q, hh); + if (hh != 0) { + h[hindex++] = hh; + } + } + if ((Q != 0.0) || (hindex == 0)) { + h[hindex++] = Q; + } + return hindex; +} + +double _estm(int elen, double *e) +{ + int eindex; + double Q = e[0]; + for (eindex = 1; eindex < elen; eindex++) Q += e[eindex]; + return Q; +} + + +double _adaptive2dorientation(double *pa, double *pb, double *pc, double detsum) +{ + double acx, acy, bcx, bcy,acxtail, acytail, bcxtail, bcytail, detleft, detright; + double detlefttail, detrighttail, det, errbound, B[5], C1[8], C2[12], D[16]; + double B3, u[5], u3, s1, t1, s0, t0, _bvr, _avr, _brn, _arn, c; + double abig, ahi, alo, bhi, blo, err1, err2, err3, _i, _j, _0; + int C1length, C2length, Dlength; + + acx = (double) (pa[0] - pc[0]); + bcx = (double) (pb[0] - pc[0]); + acy = (double) (pa[1] - pc[1]); + bcy = (double) (pb[1] - pc[1]); + + TWP(acx, bcy, detleft, detlefttail); + TWP(acy, bcx, detright, detrighttail); + + TTD(detleft, detlefttail, detright, detrighttail, + B3, B[2], B[1], B[0]); + B[3] = B3; + + det = _estm(4, B); + errbound = _ccwebB * detsum; + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + TDT(pa[0], pc[0], acx, acxtail); + TDT(pb[0], pc[0], bcx, bcxtail); + TDT(pa[1], pc[1], acy, acytail); + TDT(pb[1], pc[1], bcy, bcytail); + + if ((acxtail == 0.0) && (acytail == 0.0) + && (bcxtail == 0.0) && (bcytail == 0.0)) { + return det; + } + + errbound = _ccwebC * detsum + _reb * FABS(det); + det += (acx * bcytail + bcy * acxtail) + - (acy * bcxtail + bcx * acytail); + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + TWP(acxtail, bcy, s1, s0); + TWP(acytail, bcx, t1, t0); + TTD(s1, s0, t1, t0, u3, u[2], u[1], u[0]); + u[3] = u3; + C1length = _fesze(4, B, 4, u, C1); + + TWP(acx, bcytail, s1, s0); + TWP(acy, bcxtail, t1, t0); + TTD(s1, s0, t1, t0, u3, u[2], u[1], u[0]); + u[3] = u3; + C2length = _fesze(C1length, C1, 4, u, C2); + + TWP(acxtail, bcytail, s1, s0); + TWP(acytail, bcxtail, t1, t0); + TTD(s1, s0, t1, t0, u3, u[2], u[1], u[0]); + u[3] = u3; + Dlength = _fesze(C2length, C2, 4, u, D); + + return(D[Dlength - 1]); +} + +double _adaptive3dorientation(double *pa, double *pb, double *pc, double *pd, double permanent) +{ + double adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz, det, errbound; + double bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1; + double bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0, bc[4], ca[4], ab[4]; + double bc3, ca3, ab3, adet[8], bdet[8], cdet[8]; + double abdet[16], *finnow, *finother, *finswap, fin1[192], fin2[192]; + double adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail, adztail, bdztail, cdztail; + double at_blarge, at_clarge, bt_clarge, bt_alarge, ct_alarge, ct_blarge; + double at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4]; + double bdxt_cdy1, cdxt_bdy1, cdxt_ady1, adxt_cdy1, adxt_bdy1, bdxt_ady1; + double bdxt_cdy0, cdxt_bdy0, cdxt_ady0, adxt_cdy0, adxt_bdy0, bdxt_ady0; + double bdyt_cdx1, cdyt_bdx1, cdyt_adx1, adyt_cdx1, adyt_bdx1, bdyt_adx1; + double bdyt_cdx0, cdyt_bdx0, cdyt_adx0, adyt_cdx0, adyt_bdx0, bdyt_adx0; + double bct[8], cat[8], abt[8], bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1; + double adxt_cdyt1, adxt_bdyt1, bdxt_adyt1, bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0; + double adxt_cdyt0, adxt_bdyt0, bdxt_adyt0, u[4], v[12], w[16], u3, negate; + double _bvr, _avr, _brn, _arn, c, abig, ahi, alo, bhi, blo; + double err1, err2, err3, _i, _j, _k, _0; + int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen, finlength; + int vlength, wlength, alen, blen, clen, ablen, bctlen, catlen, abtlen; + + adx = (double) (pa[0] - pd[0]); + bdx = (double) (pb[0] - pd[0]); + cdx = (double) (pc[0] - pd[0]); + ady = (double) (pa[1] - pd[1]); + bdy = (double) (pb[1] - pd[1]); + cdy = (double) (pc[1] - pd[1]); + adz = (double) (pa[2] - pd[2]); + bdz = (double) (pb[2] - pd[2]); + cdz = (double) (pc[2] - pd[2]); + + TWP(bdx, cdy, bdxcdy1, bdxcdy0); + TWP(cdx, bdy, cdxbdy1, cdxbdy0); + TTD(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]); + bc[3] = bc3; + alen = _seze(4, bc, adz, adet); + + TWP(cdx, ady, cdxady1, cdxady0); + TWP(adx, cdy, adxcdy1, adxcdy0); + TTD(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]); + ca[3] = ca3; + blen = _seze(4, ca, bdz, bdet); + + TWP(adx, bdy, adxbdy1, adxbdy0); + TWP(bdx, ady, bdxady1, bdxady0); + TTD(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]); + ab[3] = ab3; + clen = _seze(4, ab, cdz, cdet); + + ablen = _fesze(alen, adet, blen, bdet, abdet); + finlength = _fesze(ablen, abdet, clen, cdet, fin1); + + det = _estm(finlength, fin1); + errbound = _o3ebB * permanent; + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + TDT(pa[0], pd[0], adx, adxtail); + TDT(pb[0], pd[0], bdx, bdxtail); + TDT(pc[0], pd[0], cdx, cdxtail); + TDT(pa[1], pd[1], ady, adytail); + TDT(pb[1], pd[1], bdy, bdytail); + TDT(pc[1], pd[1], cdy, cdytail); + TDT(pa[2], pd[2], adz, adztail); + TDT(pb[2], pd[2], bdz, bdztail); + TDT(pc[2], pd[2], cdz, cdztail); + + if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) + && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0) + && (adztail == 0.0) && (bdztail == 0.0) && (cdztail == 0.0)) { + return det; + } + + errbound = _o3ebC * permanent + _reb * FABS(det); + det += (adz * ((bdx * cdytail + cdy * bdxtail) + - (bdy * cdxtail + cdx * bdytail)) + + adztail * (bdx * cdy - bdy * cdx)) + + (bdz * ((cdx * adytail + ady * cdxtail) + - (cdy * adxtail + adx * cdytail)) + + bdztail * (cdx * ady - cdy * adx)) + + (cdz * ((adx * bdytail + bdy * adxtail) + - (ady * bdxtail + bdx * adytail)) + + cdztail * (adx * bdy - ady * bdx)); + if ((det >= errbound) || (-det >= errbound)) { + return det; + } + + finnow = fin1; + finother = fin2; + + if (adxtail == 0.0) { + if (adytail == 0.0) { + at_b[0] = 0.0; + at_blen = 1; + at_c[0] = 0.0; + at_clen = 1; + } else { + negate = -adytail; + TWP(negate, bdx, at_blarge, at_b[0]); + at_b[1] = at_blarge; + at_blen = 2; + TWP(adytail, cdx, at_clarge, at_c[0]); + at_c[1] = at_clarge; + at_clen = 2; + } + } else { + if (adytail == 0.0) { + TWP(adxtail, bdy, at_blarge, at_b[0]); + at_b[1] = at_blarge; + at_blen = 2; + negate = -adxtail; + TWP(negate, cdy, at_clarge, at_c[0]); + at_c[1] = at_clarge; + at_clen = 2; + } else { + TWP(adxtail, bdy, adxt_bdy1, adxt_bdy0); + TWP(adytail, bdx, adyt_bdx1, adyt_bdx0); + TTD(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0, + at_blarge, at_b[2], at_b[1], at_b[0]); + at_b[3] = at_blarge; + at_blen = 4; + TWP(adytail, cdx, adyt_cdx1, adyt_cdx0); + TWP(adxtail, cdy, adxt_cdy1, adxt_cdy0); + TTD(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0, + at_clarge, at_c[2], at_c[1], at_c[0]); + at_c[3] = at_clarge; + at_clen = 4; + } + } + if (bdxtail == 0.0) { + if (bdytail == 0.0) { + bt_c[0] = 0.0; + bt_clen = 1; + bt_a[0] = 0.0; + bt_alen = 1; + } else { + negate = -bdytail; + TWP(negate, cdx, bt_clarge, bt_c[0]); + bt_c[1] = bt_clarge; + bt_clen = 2; + TWP(bdytail, adx, bt_alarge, bt_a[0]); + bt_a[1] = bt_alarge; + bt_alen = 2; + } + } else { + if (bdytail == 0.0) { + TWP(bdxtail, cdy, bt_clarge, bt_c[0]); + bt_c[1] = bt_clarge; + bt_clen = 2; + negate = -bdxtail; + TWP(negate, ady, bt_alarge, bt_a[0]); + bt_a[1] = bt_alarge; + bt_alen = 2; + } else { + TWP(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0); + TWP(bdytail, cdx, bdyt_cdx1, bdyt_cdx0); + TTD(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0, + bt_clarge, bt_c[2], bt_c[1], bt_c[0]); + bt_c[3] = bt_clarge; + bt_clen = 4; + TWP(bdytail, adx, bdyt_adx1, bdyt_adx0); + TWP(bdxtail, ady, bdxt_ady1, bdxt_ady0); + TTD(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0, + bt_alarge, bt_a[2], bt_a[1], bt_a[0]); + bt_a[3] = bt_alarge; + bt_alen = 4; + } + } + if (cdxtail == 0.0) { + if (cdytail == 0.0) { + ct_a[0] = 0.0; + ct_alen = 1; + ct_b[0] = 0.0; + ct_blen = 1; + } else { + negate = -cdytail; + TWP(negate, adx, ct_alarge, ct_a[0]); + ct_a[1] = ct_alarge; + ct_alen = 2; + TWP(cdytail, bdx, ct_blarge, ct_b[0]); + ct_b[1] = ct_blarge; + ct_blen = 2; + } + } else { + if (cdytail == 0.0) { + TWP(cdxtail, ady, ct_alarge, ct_a[0]); + ct_a[1] = ct_alarge; + ct_alen = 2; + negate = -cdxtail; + TWP(negate, bdy, ct_blarge, ct_b[0]); + ct_b[1] = ct_blarge; + ct_blen = 2; + } else { + TWP(cdxtail, ady, cdxt_ady1, cdxt_ady0); + TWP(cdytail, adx, cdyt_adx1, cdyt_adx0); + TTD(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0, + ct_alarge, ct_a[2], ct_a[1], ct_a[0]); + ct_a[3] = ct_alarge; + ct_alen = 4; + TWP(cdytail, bdx, cdyt_bdx1, cdyt_bdx0); + TWP(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0); + TTD(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0, + ct_blarge, ct_b[2], ct_b[1], ct_b[0]); + ct_b[3] = ct_blarge; + ct_blen = 4; + } + } + + bctlen = _fesze(bt_clen, bt_c, ct_blen, ct_b, bct); + wlength = _seze(bctlen, bct, adz, w); + finlength = _fesze(finlength, finnow, wlength, w, + finother); + finswap = finnow; finnow = finother; finother = finswap; + + catlen = _fesze(ct_alen, ct_a, at_clen, at_c, cat); + wlength = _seze(catlen, cat, bdz, w); + finlength = _fesze(finlength, finnow, wlength, w, + finother); + finswap = finnow; finnow = finother; finother = finswap; + + abtlen = _fesze(at_blen, at_b, bt_alen, bt_a, abt); + wlength = _seze(abtlen, abt, cdz, w); + finlength = _fesze(finlength, finnow, wlength, w, + finother); + finswap = finnow; finnow = finother; finother = finswap; + + if (adztail != 0.0) { + vlength = _seze(4, bc, adztail, v); + finlength = _fesze(finlength, finnow, vlength, v, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdztail != 0.0) { + vlength = _seze(4, ca, bdztail, v); + finlength = _fesze(finlength, finnow, vlength, v, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdztail != 0.0) { + vlength = _seze(4, ab, cdztail, v); + finlength = _fesze(finlength, finnow, vlength, v, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + + if (adxtail != 0.0) { + if (bdytail != 0.0) { + TWP(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0); + TOP(adxt_bdyt1, adxt_bdyt0, cdz, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = _fesze(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + if (cdztail != 0.0) { + TOP(adxt_bdyt1, adxt_bdyt0, cdztail, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = _fesze(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + if (cdytail != 0.0) { + negate = -adxtail; + TWP(negate, cdytail, adxt_cdyt1, adxt_cdyt0); + TOP(adxt_cdyt1, adxt_cdyt0, bdz, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = _fesze(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + if (bdztail != 0.0) { + TOP(adxt_cdyt1, adxt_cdyt0, bdztail, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = _fesze(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + } + if (bdxtail != 0.0) { + if (cdytail != 0.0) { + TWP(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0); + TOP(bdxt_cdyt1, bdxt_cdyt0, adz, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = _fesze(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + if (adztail != 0.0) { + TOP(bdxt_cdyt1, bdxt_cdyt0, adztail, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = _fesze(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + if (adytail != 0.0) { + negate = -bdxtail; + TWP(negate, adytail, bdxt_adyt1, bdxt_adyt0); + TOP(bdxt_adyt1, bdxt_adyt0, cdz, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = _fesze(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + if (cdztail != 0.0) { + TOP(bdxt_adyt1, bdxt_adyt0, cdztail, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = _fesze(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + } + if (cdxtail != 0.0) { + if (adytail != 0.0) { + TWP(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0); + TOP(cdxt_adyt1, cdxt_adyt0, bdz, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = _fesze(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + if (bdztail != 0.0) { + TOP(cdxt_adyt1, cdxt_adyt0, bdztail, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = _fesze(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + if (bdytail != 0.0) { + negate = -cdxtail; + TWP(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0); + TOP(cdxt_bdyt1, cdxt_bdyt0, adz, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = _fesze(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + if (adztail != 0.0) { + TOP(cdxt_bdyt1, cdxt_bdyt0, adztail, u3, u[2], u[1], u[0]); + u[3] = u3; + finlength = _fesze(finlength, finnow, 4, u, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + } + + if (adztail != 0.0) { + wlength = _seze(bctlen, bct, adztail, w); + finlength = _fesze(finlength, finnow, wlength, w, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdztail != 0.0) { + wlength = _seze(catlen, cat, bdztail, w); + finlength = _fesze(finlength, finnow, wlength, w, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdztail != 0.0) { + wlength = _seze(abtlen, abt, cdztail, w); + finlength = _fesze(finlength, finnow, wlength, w, + finother); + finswap = finnow; finnow = finother; finother = finswap; + } + + return finnow[finlength - 1]; +} + + + + + +/*****************************************************************************/ +/* */ +/* PUBLIC FUNCTIONS */ +/* initPredicates() Sets the variables used for exact arithmetic. This */ +/* must be called once before using the other two functions. */ +/* orient2d() Computes the orientation of three 2D points. */ +/* orient3d() Computes the orientation of four 3D points. */ +/* */ +/*****************************************************************************/ + +void initPredicates() +{ + static char a_c=0; + double hf, ck, lc; + int e_o; + + if (a_c) return; else a_c = 1; + +#ifdef SPECIFY_FP_PRECISION + unsigned int old_cfp; + _controlfp_s(&old_cfp, _PC_53, MCW_PC); +#endif + + e_o = 1; + _eps = _spl = ck = 1.0; + hf = 0.5; + + do + { + lc=ck; + _eps *= hf; + if (e_o) _spl *= 2.0; + e_o = !e_o; + ck = 1.0 + _eps; + } while ((ck != 1.0) && (ck != lc)); + _spl += 1.0; + + _reb = (3.0 + 8.0 * _eps) * _eps; + _ccwebA = (3.0 + 16.0 * _eps) * _eps; + _ccwebB = (2.0 + 12.0 * _eps) * _eps; + _ccwebC = (9.0 + 64.0 * _eps) * _eps * _eps; + _o3ebA = (7.0 + 56.0 * _eps) * _eps; + _o3ebB = (3.0 + 28.0 * _eps) * _eps; + _o3ebC = (26.0 + 288.0 * _eps) * _eps * _eps; + _iccebA = (10.0 + 96.0 * _eps) * _eps; + _iccebB = (4.0 + 48.0 * _eps) * _eps; + _iccebC = (44.0 + 576.0 * _eps) * _eps * _eps; + _ispebA = (16.0 + 224.0 * _eps) * _eps; + _ispebB = (5.0 + 72.0 * _eps) * _eps; + _ispebC = (71.0 + 1408.0 * _eps) * _eps * _eps; + +#ifdef SPECIFY_FP_PRECISION + _controlfp_s(&old_cfp, _CW_DEFAULT, MCW_PC); +#endif +} + +double orient2d(double *pa, double *pb, double *pc) +{ + double dlf, drg, det, dsm, eb; + + dlf = (pa[0]-pc[0])*(pb[1]-pc[1]); + drg = (pa[1]-pc[1])*(pb[0]-pc[0]); + det = dlf - drg; + + if (dlf > 0.0) {if (drg <= 0.0) return det; else dsm = dlf + drg;} + else if (dlf < 0.0) {if (drg >= 0.0) return det; else dsm = -dlf - drg;} + else return det; + + eb = _ccwebA*dsm; + if ((det>=eb) || (-det>=eb)) return det; + + return _adaptive2dorientation(pa, pb, pc, dsm); +} + +double orient3d(double *pa, double *pb, double *pc, double *pd) +{ + double adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz, pm, eb; + double bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady, det; + + adx = pa[0]-pd[0]; bdx = pb[0]-pd[0]; cdx = pc[0]-pd[0]; + ady = pa[1]-pd[1]; bdy = pb[1]-pd[1]; cdy = pc[1]-pd[1]; + adz = pa[2]-pd[2]; bdz = pb[2]-pd[2]; cdz = pc[2]-pd[2]; + + bdxcdy = bdx*cdy; cdxbdy = cdx*bdy; + cdxady = cdx*ady; adxcdy = adx*cdy; + adxbdy = adx*bdy; bdxady = bdx*ady; + + det = adz*(bdxcdy-cdxbdy)+bdz*(cdxady-adxcdy)+cdz*(adxbdy-bdxady); + pm=(FABS(bdxcdy)+FABS(cdxbdy))*FABS(adz)+(FABS(cdxady)+FABS(adxcdy))*FABS(bdz)+(FABS(adxbdy)+FABS(bdxady))*FABS(cdz); + eb = _o3ebA*pm; + if ((det>eb) || (-det>eb)) return det; + return _adaptive3dorientation(pa, pb, pc, pd, pm); +} diff --git a/src/mesh_fix/src/Kernel/tmesh.cpp b/src/mesh_fix/src/Kernel/tmesh.cpp index 97ab75a06..c301cebeb 100644 --- a/src/mesh_fix/src/Kernel/tmesh.cpp +++ b/src/mesh_fix/src/Kernel/tmesh.cpp @@ -253,7 +253,7 @@ void TMesh::printElapsedTime(bool reset) { static clock_t beginning_time; if (reset) beginning_time = clock(); - else printf("\n\n********** PARTIAL ELAPSED: %d msecs\n\n", (clock() - beginning_time)); + else printf("\n\n********** PARTIAL ELAPSED: %d msecs\n\n", int((clock() - beginning_time))); } } //namespace T_MESH diff --git a/src/mesh_fix/src/MeshFix/meshfix.cpp b/src/mesh_fix/src/MeshFix/meshfix.cpp index d05a7cbaa..d864c7964 100644 --- a/src/mesh_fix/src/MeshFix/meshfix.cpp +++ b/src/mesh_fix/src/MeshFix/meshfix.cpp @@ -39,14 +39,14 @@ bool joinClosestComponents(Basic_TMesh *tin) { i++; triList.appendHead(t); - t->info = (void *)(intptr_t)i; + t->info = new intWrapper(i); while (triList.numels()) { t = (Triangle *)triList.popHead(); - if ((s = t->t1()) != NULL && s->info == NULL) { triList.appendHead(s); s->info = (void *)(intptr_t)i; } - if ((s = t->t2()) != NULL && s->info == NULL) { triList.appendHead(s); s->info = (void *)(intptr_t)i; } - if ((s = t->t3()) != NULL && s->info == NULL) { triList.appendHead(s); s->info = (void *)(intptr_t)i; } + if ((s = t->t1()) != NULL && s->info == NULL) { triList.appendHead(s); s->info = new intWrapper(i); } + if ((s = t->t2()) != NULL && s->info == NULL) { triList.appendHead(s); s->info = new intWrapper(i); } + if ((s = t->t3()) != NULL && s->info == NULL) { triList.appendHead(s); s->info = new intWrapper(i); } } } diff --git a/src/mesh_fix/src/TMesh/io.cpp b/src/mesh_fix/src/TMesh/io.cpp index 7efae0131..1dda5ff41 100644 --- a/src/mesh_fix/src/TMesh/io.cpp +++ b/src/mesh_fix/src/TMesh/io.cpp @@ -722,7 +722,7 @@ int Basic_TMesh::saveVRML1(const char *fname, const int mode) fprintf(fp,"Material {\n diffuseColor [\n"); FOREACHTRIANGLE(t, n) { - pkc = (unsigned int)((j_voidint)t->info); + pkc = ((intWrapper*)(t->info))->operator int(); fprintf(fp," %f %f %f,\n",((pkc>>24)&0x000000ff)/255.0,((pkc>>16)&0x000000ff)/255.0,((pkc>>8)&0x000000ff)/255.0); } fprintf(fp," ]\n}\nMaterialBinding {\n value PER_FACE_INDEXED\n}\n"); @@ -731,7 +731,7 @@ int Basic_TMesh::saveVRML1(const char *fname, const int mode) fprintf(fp,"Material {\n diffuseColor [\n"); FOREACHVERTEX(v, n) { - pkc = (unsigned int)((j_voidint)v->info); + pkc = ((intWrapper*)(v->info))->operator int(); fprintf(fp," %f %f %f,\n",((pkc>>24)&0x000000ff)/255.0,((pkc>>16)&0x000000ff)/255.0,((pkc>>8)&0x000000ff)/255.0); } fprintf(fp," ]\n}\nMaterialBinding {\n value PER_VERTEX_INDEXED\n}\n"); @@ -893,7 +893,7 @@ int Basic_TMesh::saveVerTri(const char *fname) ocds = new coord[V.numels()]; i=0; FOREACHVERTEX(v, n) ocds[i++] = v->x; i=0; FOREACHVERTEX(v, n) v->x = ++i; - i=0; FOREACHTRIANGLE(t, n) {i++; t->info = (void *)(intptr_t)i;} + i=0; FOREACHTRIANGLE(t, n) {i++; t->info = new intWrapper(i);} fprintf(fpt,"%d\n",T.numels()); FOREACHTRIANGLE(t, n) @@ -1460,7 +1460,7 @@ int Basic_TMesh::loadSTL(const char *fname) { FILE *fp; int nt=0, i=0; - char kw[64]="", kw2[64]="", *line, facet[50]; + char kw[65]="", kw2[65]="", *line, facet[50]; float x,y,z; bool binary=0; Vertex *v, *v1=NULL, *v2=NULL, *v3=NULL; diff --git a/src/mesh_fix/src/TMesh/tin.cpp b/src/mesh_fix/src/TMesh/tin.cpp index 3240e17d9..7cd3d8913 100644 --- a/src/mesh_fix/src/TMesh/tin.cpp +++ b/src/mesh_fix/src/TMesh/tin.cpp @@ -187,11 +187,11 @@ void Basic_TMesh::init(const Basic_TMesh *tin, const bool clone_info) Triangle *t, *nt; int i; - void **t_info = new void *[tin->T.numels()]; + Data **t_info = new Data *[tin->T.numels()]; i=0; FOREACHVTTRIANGLE((&(tin->T)), t, n) t_info[i++]=t->info; - void **e_info = new void *[tin->E.numels()]; + Data **e_info = new Data *[tin->E.numels()]; i=0; FOREACHVEEDGE((&(tin->E)), e, n) e_info[i++]=e->info; - void **v_info = new void *[tin->V.numels()]; + Data **v_info = new Data *[tin->V.numels()]; i=0; FOREACHVVVERTEX((&(tin->V)), v, n) v_info[i++]=v->info; FOREACHVVVERTEX((&(tin->V)), v, n) @@ -838,14 +838,14 @@ Basic_TMesh *Basic_TMesh::createSubMeshFromSelection(Triangle *t0, bool keep_ref FOREACHVEEDGE((&sE), e, n) e->v1->e0 = e->v2->e0 = e; int i; - void **v_info = NULL, **e_info = NULL, **t_info = NULL; + Data **v_info = NULL, **e_info = NULL, **t_info = NULL; if (!keep_ref) { - v_info = new void *[sV.numels()]; + v_info = new Data *[sV.numels()]; i=0; FOREACHVVVERTEX((&sV), v, n) v_info[i++] = v->info; - e_info = new void *[sE.numels()]; + e_info = new Data *[sE.numels()]; i=0; FOREACHVEEDGE((&sE), e, n) e_info[i++] = e->info; - t_info = new void *[sT.numels()]; + t_info = new Data *[sT.numels()]; i=0; FOREACHVTTRIANGLE((&sT), t, n) t_info[i++] = t->info; } diff --git a/src/slic3r/Utils/FixModelByMeshFix.cpp b/src/slic3r/Utils/FixModelByMeshFix.cpp index ba39c9e83..150b7c103 100644 --- a/src/slic3r/Utils/FixModelByMeshFix.cpp +++ b/src/slic3r/Utils/FixModelByMeshFix.cpp @@ -50,6 +50,8 @@ namespace detail { using namespace T_MESH; +static int ensure_tmesh_initialization = (TMesh::init(), 0); + double closestPair(List *bl1, List *bl2, Vertex **closest_on_bl1, Vertex **closest_on_bl2) { Node *n, *m; @@ -84,14 +86,14 @@ bool joinClosestComponents(Basic_TMesh *tin) { i++; triList.appendHead(t); - t->info = (void *)(intptr_t)i; + t->info = new intWrapper(i); while (triList.numels()) { t = (Triangle *)triList.popHead(); - if ((s = t->t1()) != NULL && s->info == NULL) {triList.appendHead(s); s->info = (void *)(intptr_t)i;} - if ((s = t->t2()) != NULL && s->info == NULL) {triList.appendHead(s); s->info = (void *)(intptr_t)i;} - if ((s = t->t3()) != NULL && s->info == NULL) {triList.appendHead(s); s->info = (void *)(intptr_t)i;} + if ((s = t->t1()) != NULL && s->info == NULL) {triList.appendHead(s); s->info = new intWrapper(i);} + if ((s = t->t2()) != NULL && s->info == NULL) {triList.appendHead(s); s->info = new intWrapper(i);} + if ((s = t->t3()) != NULL && s->info == NULL) {triList.appendHead(s); s->info = new intWrapper(i);} } } @@ -298,7 +300,8 @@ bool fix_model_by_meshfix(ModelObject &model_object, int volume_idx, wxProgressD tin.deselectTriangles(); tin.boundaries(); // Keep only the largest component (i.e. with most triangles) - on_progress(L("Remove smallest components"), unsigned(progress_part_base + 0.2 * percent_per_part)); + on_progress(L("Remove smallest components"), + unsigned(progress_part_base + 0.2 * percent_per_part)); if (canceled) throw RepairCanceledException(); tin.removeSmallestComponents(); @@ -308,7 +311,8 @@ bool fix_model_by_meshfix(ModelObject &model_object, int volume_idx, wxProgressD if (canceled) throw RepairCanceledException(); if (tin.boundaries()) { - on_progress(L("Patch small holes"), unsigned(progress_part_base + 0.4 * percent_per_part)); + on_progress(L("Patch small holes"), + unsigned(progress_part_base + 0.4 * percent_per_part)); if (canceled) throw RepairCanceledException(); tin.fillSmallBoundaries(0, true); @@ -320,22 +324,25 @@ bool fix_model_by_meshfix(ModelObject &model_object, int volume_idx, wxProgressD // Run geometry correction if (!tin.boundaries()) { int iteration = 0; - on_progress(L("Start iterative correction"), unsigned(progress_part_base + 0.55 * percent_per_part)); + on_progress(L("Start iterative correction"), + unsigned(progress_part_base + 0.55 * percent_per_part)); tin.deselectTriangles(); tin.invertSelection(); bool fixed = false; while (iteration < 10 && !fixed) { //default constants taken from TMesh library fixed = tin.meshclean_single_iteration(3, canceled); - on_progress(L("Fixing geometry"), progress_part_base + percent_per_part * std::min(0.9, 0.6 + iteration*0.08)); // majority of objects should finish in 4 iterations + on_progress(L("Fixing geometry"), + progress_part_base + + percent_per_part * std::min(0.9, 0.6 + iteration * 0.08)); // majority of objects will finish in 4 iterations if (canceled) throw RepairCanceledException(); iteration++; } } - if (tin.boundaries() || tin.T.numels() == 0) { - throw RepairFailedException(); - } +// if (tin.boundaries() || tin.T.numels() == 0) { +// throw RepairFailedException(); +// } parts[part_idx] = tin.to_indexed_triangle_set(); }