diff options
Diffstat (limited to 'Build/source/utils/asymptote/path3.cc')
-rw-r--r-- | Build/source/utils/asymptote/path3.cc | 1460 |
1 files changed, 1460 insertions, 0 deletions
diff --git a/Build/source/utils/asymptote/path3.cc b/Build/source/utils/asymptote/path3.cc new file mode 100644 index 00000000000..417fed23845 --- /dev/null +++ b/Build/source/utils/asymptote/path3.cc @@ -0,0 +1,1460 @@ +/***** + * path3.cc + * John Bowman + * + * Compute information for a three-dimensional path. + *****/ + +#include <cfloat> + +#include "path3.h" +#include "util.h" +#include "camperror.h" +#include "mathop.h" + +namespace camp { + +using run::operator *; +using vm::array; + +path3 nullpath3; + +void checkEmpty3(Int n) { + if(n == 0) + reportError("nullpath3 has no points"); +} + +triple path3::point(double t) const +{ + checkEmpty3(n); + + Int i = Floor(t); + Int iplus; + t = fmod(t,1); + if (t < 0) t += 1; + + if (cycles) { + i = imod(i,n); + iplus = imod(i+1,n); + } + else if (i < 0) + return nodes[0].point; + else if (i >= n-1) + return nodes[n-1].point; + else + iplus = i+1; + + double one_t = 1.0-t; + + triple a = nodes[i].point, + b = nodes[i].post, + c = nodes[iplus].pre, + d = nodes[iplus].point, + ab = one_t*a + t*b, + bc = one_t*b + t*c, + cd = one_t*c + t*d, + abc = one_t*ab + t*bc, + bcd = one_t*bc + t*cd, + abcd = one_t*abc + t*bcd; + + return abcd; +} + +triple path3::precontrol(double t) const +{ + checkEmpty3(n); + + Int i = Floor(t); + Int iplus; + t = fmod(t,1); + if (t < 0) t += 1; + + if (cycles) { + i = imod(i,n); + iplus = imod(i+1,n); + } + else if (i < 0) + return nodes[0].pre; + else if (i >= n-1) + return nodes[n-1].pre; + else + iplus = i+1; + + double one_t = 1.0-t; + + triple a = nodes[i].point, + b = nodes[i].post, + c = nodes[iplus].pre, + ab = one_t*a + t*b, + bc = one_t*b + t*c, + abc = one_t*ab + t*bc; + + return (abc == a) ? nodes[i].pre : abc; +} + + +triple path3::postcontrol(double t) const +{ + checkEmpty3(n); + + Int i = Floor(t); + Int iplus; + t = fmod(t,1); + if (t < 0) t += 1; + + if (cycles) { + i = imod(i,n); + iplus = imod(i+1,n); + } + else if (i < 0) + return nodes[0].post; + else if (i >= n-1) + return nodes[n-1].post; + else + iplus = i+1; + + double one_t = 1.0-t; + + triple b = nodes[i].post, + c = nodes[iplus].pre, + d = nodes[iplus].point, + bc = one_t*b + t*c, + cd = one_t*c + t*d, + bcd = one_t*bc + t*cd; + + return (bcd == d) ? nodes[iplus].post : bcd; +} + +path3 path3::reverse() const +{ + mem::vector<solvedKnot3> nodes(n); + Int len=length(); + for (Int i = 0, j = len; i < n; i++, j--) { + nodes[i].pre = postcontrol(j); + nodes[i].point = point(j); + nodes[i].post = precontrol(j); + nodes[i].straight = straight(j-1); + } + return path3(nodes, n, cycles); +} + +path3 path3::subpath(Int a, Int b) const +{ + if(empty()) return path3(); + + if (a > b) { + const path3 &rp = reverse(); + Int len=length(); + path3 result = rp.subpath(len-a, len-b); + return result; + } + + if (!cycles) { + if (a < 0) { + a = 0; + if(b < 0) + b = 0; + } + if (b > n-1) { + b = n-1; + if(a > b) + a = b; + } + } + + Int sn = b-a+1; + mem::vector<solvedKnot3> nodes(sn); + + for (Int i = 0, j = a; j <= b; i++, j++) { + nodes[i].pre = precontrol(j); + nodes[i].point = point(j); + nodes[i].post = postcontrol(j); + nodes[i].straight = straight(j); + } + nodes[0].pre = nodes[0].point; + nodes[sn-1].post = nodes[sn-1].point; + + return path3(nodes, sn); +} + +inline triple split(double t, const triple& x, const triple& y) { + return x+(y-x)*t; +} + +inline void splitCubic(solvedKnot3 sn[], double t, const solvedKnot3& left_, + const solvedKnot3& right_) +{ + solvedKnot3 &left=(sn[0]=left_), &mid=sn[1], &right=(sn[2]=right_); + if(left.straight) { + mid.point=split(t,left.point,right.point); + triple deltaL=third*(mid.point-left.point); + left.post=left.point+deltaL; + mid.pre=mid.point-deltaL; + triple deltaR=third*(right.point-mid.point); + mid.post=mid.point+deltaR; + right.pre=right.point-deltaR; + mid.straight=true; + } else { + triple x=split(t,left.post,right.pre); // m1 + left.post=split(t,left.point,left.post); // m0 + right.pre=split(t,right.pre,right.point); // m2 + mid.pre=split(t,left.post,x); // m3 + mid.post=split(t,x,right.pre); // m4 + mid.point=split(t,mid.pre,mid.post); // m5 + } +} + +path3 path3::subpath(double a, double b) const +{ + if(empty()) return path3(); + + if (a > b) { + const path3 &rp = reverse(); + Int len=length(); + return rp.subpath(len-a, len-b); + } + + solvedKnot3 aL, aR, bL, bR; + if (!cycles) { + if (a < 0) { + a = 0; + if (b < 0) + b = 0; + } + if (b > n-1) { + b = n-1; + if (a > b) + a = b; + } + aL = nodes[(Int)floor(a)]; + aR = nodes[(Int)ceil(a)]; + bL = nodes[(Int)floor(b)]; + bR = nodes[(Int)ceil(b)]; + } else { + if(run::validInt(a) && run::validInt(b)) { + aL = nodes[imod((Int) floor(a),n)]; + aR = nodes[imod((Int) ceil(a),n)]; + bL = nodes[imod((Int) floor(b),n)]; + bR = nodes[imod((Int) ceil(b),n)]; + } else reportError("invalid path3 index"); + } + + if (a == b) return path3(point(a)); + + solvedKnot3 sn[3]; + path3 p = subpath(Ceil(a), Floor(b)); + if (a > floor(a)) { + if (b < ceil(a)) { + splitCubic(sn,a-floor(a),aL,aR); + splitCubic(sn,(b-a)/(ceil(b)-a),sn[1],sn[2]); + return path3(sn[0],sn[1]); + } + splitCubic(sn,a-floor(a),aL,aR); + p=concat(path3(sn[1],sn[2]),p); + } + if (ceil(b) > b) { + splitCubic(sn,b-floor(b),bL,bR); + p=concat(p,path3(sn[0],sn[1])); + } + return p; +} + +// Special case of subpath for paths of length 1 used by intersect. +void path3::halve(path3 &first, path3 &second) const +{ + solvedKnot3 sn[3]; + splitCubic(sn,0.5,nodes[0],nodes[1]); + first=path3(sn[0],sn[1]); + second=path3(sn[1],sn[2]); +} + +// Calculate the coefficients of a Bezier derivative divided by 3. +static inline void derivative(triple& a, triple& b, triple& c, + const triple& z0, const triple& c0, + const triple& c1, const triple& z1) +{ + a=z1-z0+3.0*(c0-c1); + b=2.0*(z0+c1)-4.0*c0; + c=c0-z0; +} + +bbox3 path3::bounds() const +{ + if(!box.empty) return box; + + if (empty()) { + // No bounds + return bbox3(); + } + + Int len=length(); + box.add(point(len)); + times=bbox3(len,len,len,len,len,len); + + for (Int i = 0; i < len; i++) { + addpoint(box,i); + if(straight(i)) continue; + + triple a,b,c; + derivative(a,b,c,point(i),postcontrol(i),precontrol(i+1),point(i+1)); + + // Check x coordinate + quadraticroots x(a.getx(),b.getx(),c.getx()); + if(x.distinct != quadraticroots::NONE && goodroot(x.t1)) + addpoint(box,i+x.t1); + if(x.distinct == quadraticroots::TWO && goodroot(x.t2)) + addpoint(box,i+x.t2); + + // Check y coordinate + quadraticroots y(a.gety(),b.gety(),c.gety()); + if(y.distinct != quadraticroots::NONE && goodroot(y.t1)) + addpoint(box,i+y.t1); + if(y.distinct == quadraticroots::TWO && goodroot(y.t2)) + addpoint(box,i+y.t2); + + // Check z coordinate + quadraticroots z(a.getz(),b.getz(),c.getz()); + if(z.distinct != quadraticroots::NONE && goodroot(z.t1)) + addpoint(box,i+z.t1); + if(z.distinct == quadraticroots::TWO && goodroot(z.t2)) + addpoint(box,i+z.t2); + } + return box; +} + +// Return f evaluated at controlling vertex of bounding box of convex hull for +// similiar-triangle transform x'=x/z, y'=y/z, where z < 0. +double ratiobound(triple z0, triple c0, triple c1, triple z1, + double (*m)(double, double), + double (*f)(const triple&)) +{ + double MX=m(m(m(-z0.getx(),-c0.getx()),-c1.getx()),-z1.getx()); + double MY=m(m(m(-z0.gety(),-c0.gety()),-c1.gety()),-z1.gety()); + double Z=m(m(m(z0.getz(),c0.getz()),c1.getz()),z1.getz()); + double MZ=m(m(m(-z0.getz(),-c0.getz()),-c1.getz()),-z1.getz()); + return m(f(triple(-MX,-MY,Z)),f(triple(-MX,-MY,-MZ))); +} + +double bound(triple z0, triple c0, triple c1, triple z1, + double (*m)(double, double), + double (*f)(const triple&), double b, double fuzz, int depth) +{ + b=m(b,m(f(z0),f(z1))); + if(m(-1.0,1.0)*(b-ratiobound(z0,c0,c1,z1,m,f)) >= -fuzz || depth == 0) + return b; + + --depth; + fuzz *= 2; + + triple m0=0.5*(z0+c0); + triple m1=0.5*(c0+c1); + triple m2=0.5*(c1+z1); + triple m3=0.5*(m0+m1); + triple m4=0.5*(m1+m2); + triple m5=0.5*(m3+m4); + + // Check both Bezier subpaths. + b=bound(z0,m0,m3,m5,m,f,b,fuzz,depth); + return bound(m5,m4,m2,z1,m,f,b,fuzz,depth); +} + +pair path3::ratio(double (*m)(double, double)) const +{ + double fuzz=Fuzz*(max()-min()).length(); + checkEmpty3(n); + + triple v=point((Int) 0); + pair B=pair(xratio(v),yratio(v)); + + Int n=length(); + for(Int i=0; i <= n; ++i) { + if(straight(i)) { + triple v=point(i); + B=pair(m(B.getx(),xratio(v)),m(B.gety(),yratio(v))); + } else { + triple z0=point(i); + triple c0=postcontrol(i); + triple c1=precontrol(i+1); + triple z1=point(i+1); + B=pair(bound(z0,c0,c1,z1,m,xratio,B.getx(),fuzz), + bound(z0,c0,c1,z1,m,yratio,B.gety(),fuzz)); + } + } + return B; +} + +// {{{ Arclength Calculations + +static triple a,b,c; + +static double ds(double t) +{ + double dx=quadratic(a.getx(),b.getx(),c.getx(),t); + double dy=quadratic(a.gety(),b.gety(),c.gety(),t); + double dz=quadratic(a.getz(),b.getz(),c.getz(),t); + return sqrt(dx*dx+dy*dy+dz*dz); +} + +// Calculates arclength of a cubic Bezier curve using adaptive Simpson +// integration. +double arcLength(const triple& z0, const triple& c0, const triple& c1, + const triple& z1) +{ + double integral; + derivative(a,b,c,z0,c0,c1,z1); + + if(!simpson(integral,ds,0.0,1.0,DBL_EPSILON,1.0)) + reportError("nesting capacity exceeded in computing arclength"); + return integral; +} + +double path3::cubiclength(Int i, double goal) const +{ + const triple& z0=point(i); + const triple& z1=point(i+1); + double L; + if(straight(i)) { + L=(z1-z0).length(); + return (goal < 0 || goal >= L) ? L : -goal/L; + } + + double integral=arcLength(z0,postcontrol(i),precontrol(i+1),z1); + + L=3.0*integral; + if(goal < 0 || goal >= L) return L; + + double t=goal/L; + goal *= third; + static double dxmin=sqrt(DBL_EPSILON); + if(!unsimpson(goal,ds,0.0,t,100.0*DBL_EPSILON,integral,1.0,dxmin)) + reportError("nesting capacity exceeded in computing arctime"); + return -t; +} + +double path3::arclength() const +{ + if (cached_length != -1) return cached_length; + + double L=0.0; + for (Int i = 0; i < n-1; i++) { + L += cubiclength(i); + } + if(cycles) L += cubiclength(n-1); + cached_length = L; + return cached_length; +} + +double path3::arctime(double goal) const +{ + if (cycles) { + if (goal == 0 || cached_length == 0) return 0; + if (goal < 0) { + const path3 &rp = this->reverse(); + double result = -rp.arctime(-goal); + return result; + } + if (cached_length > 0 && goal >= cached_length) { + Int loops = (Int)(goal / cached_length); + goal -= loops*cached_length; + return loops*n+arctime(goal); + } + } else { + if (goal <= 0) + return 0; + if (cached_length > 0 && goal >= cached_length) + return n-1; + } + + double l,L=0; + for (Int i = 0; i < n-1; i++) { + l = cubiclength(i,goal); + if (l < 0) + return (-l+i); + else { + L += l; + goal -= l; + if (goal <= 0) + return i+1; + } + } + if (cycles) { + l = cubiclength(n-1,goal); + if (l < 0) + return -l+n-1; + if (cached_length > 0 && cached_length != L+l) { + reportError("arclength != length.\n" + "path3::arclength(double) must have broken semantics.\n" + "Please report this error."); + } + cached_length = L += l; + goal -= l; + return arctime(goal)+n; + } + else { + cached_length = L; + return length(); + } +} + +// }}} + +// {{{ Path3 Intersection Calculations + +// Return all intersection times of path3 g with the triple v. +void intersections(std::vector<double>& T, const path3& g, const triple& v, + double fuzz) +{ + double fuzz2=fuzz*fuzz; + Int n=g.length(); + bool cycles=g.cyclic(); + for(Int i=0; i < n; ++i) { + // Check all directions to circumvent degeneracy. + std::vector<double> r; + roots(r,g.point(i).getx(),g.postcontrol(i).getx(), + g.precontrol(i+1).getx(),g.point(i+1).getx(),v.getx()); + roots(r,g.point(i).gety(),g.postcontrol(i).gety(), + g.precontrol(i+1).gety(),g.point(i+1).gety(),v.gety()); + roots(r,g.point(i).getz(),g.postcontrol(i).getz(), + g.precontrol(i+1).getz(),g.point(i+1).getz(),v.getz()); + + size_t m=r.size(); + for(size_t j=0 ; j < m; ++j) { + double t=r[j]; + if(t >= -Fuzz2 && t <= 1.0+Fuzz2) { + double s=i+t; + if((g.point(s)-v).abs2() <= fuzz2) { + if(cycles && s >= n-Fuzz2) s=0; + T.push_back(s); + } + } + } + } +} + +// An optimized implementation of intersections(g,p--q); +// if there are an infinite number of intersection points, the returned list is +// only guaranteed to include the endpoint times of the intersection. +void intersections(std::vector<double>& S, std::vector<double>& T, + const path3& g, const triple& p, double fuzz) +{ + std::vector<double> S1; + intersections(S1,g,p,fuzz); + size_t n=S1.size(); + for(size_t i=0; i < n; ++i) { + S.push_back(S1[i]); + T.push_back(0.0); + } +} + +void add(std::vector<double>& S, std::vector<double>& T, double s, double t, + const path3& p, const path3& q, double fuzz2) +{ + triple P=p.point(s); + for(size_t i=0; i < S.size(); ++i) + if((p.point(S[i])-P).abs2() <= fuzz2) return; + S.push_back(s); + T.push_back(t); +} + +void add(double& s, double& t, std::vector<double>& S, std::vector<double>& T, + std::vector<double>& S1, std::vector<double>& T1, + double pscale, double qscale, double poffset, double qoffset, + const path3& p, const path3& q, double fuzz2, bool single) +{ + if(single) { + s=s*pscale+poffset; + t=t*qscale+qoffset; + } else { + size_t n=S1.size(); + for(size_t i=0; i < n; ++i) + add(S,T,pscale*S1[i]+poffset,qscale*T1[i]+qoffset,p,q,fuzz2); + } +} + +void add(double& s, double& t, std::vector<double>& S, std::vector<double>& T, + std::vector<double>& S1, std::vector<double>& T1, + const path3& p, const path3& q, double fuzz2, bool single) +{ + size_t n=S1.size(); + if(single) { + if(n > 0) { + s=S1[0]; + t=T1[0]; + } + } else { + for(size_t i=0; i < n; ++i) + add(S,T,S1[i],T1[i],p,q,fuzz2); + } +} + +bool intersections(double &s, double &t, std::vector<double>& S, + std::vector<double>& T, path3& p, path3& q, + double fuzz, bool single, bool exact, unsigned depth) +{ + if(errorstream::interrupt) throw interrupted(); + + double fuzz2=max(fuzzFactor*fuzz*fuzz,Fuzz2); + + Int lp=p.length(); + if(lp == 0 && exact) { + std::vector<double> T1,S1; + intersections(T1,S1,q,p.point(lp),fuzz); + add(s,t,S,T,S1,T1,p,q,fuzz2,single); + return S1.size() > 0; + } + + Int lq=q.length(); + if(lq == 0 && exact) { + std::vector<double> S1,T1; + intersections(S1,T1,p,q.point(lq),fuzz); + add(s,t,S,T,S1,T1,p,q,fuzz2,single); + return S1.size() > 0; + } + + triple maxp=p.max(); + triple minp=p.min(); + triple maxq=q.max(); + triple minq=q.min(); + + if(maxp.getx()+fuzz >= minq.getx() && + maxp.gety()+fuzz >= minq.gety() && + maxp.getz()+fuzz >= minq.getz() && + maxq.getx()+fuzz >= minp.getx() && + maxq.gety()+fuzz >= minp.gety() && + maxq.getz()+fuzz >= minp.getz()) { + // Overlapping bounding boxes + + --depth; +// fuzz *= 2; + + if((maxp-minp).length()+(maxq-minq).length() <= fuzz || depth == 0) { + if(single) { + s=0.5; + t=0.5; + } else { + S.push_back(0.5); + T.push_back(0.5); + } + return true; + } + + path3 p1,p2; + double pscale,poffset; + + std::vector<double> S1,T1; + +// fuzz2=max(fuzzFactor*fuzz*fuzz,Fuzz2); + + if(lp <= 1) { + if(lp == 1) p.halve(p1,p2); + if(lp == 0 || p1 == p || p2 == p) { + intersections(T1,S1,q,p.point((Int) 0),fuzz); + add(s,t,S,T,S1,T1,p,q,fuzz2,single); + return S1.size() > 0; + } + pscale=poffset=0.5; + } else { + Int tp=lp/2; + p1=p.subpath(0,tp); + p2=p.subpath(tp,lp); + poffset=tp; + pscale=1.0; + } + + path3 q1,q2; + double qscale,qoffset; + + if(lq <= 1) { + if(lq == 1) q.halve(q1,q2); + if(lq == 0 || q1 == q || q2 == q) { + intersections(S1,T1,p,q.point((Int) 0),fuzz); + add(s,t,S,T,S1,T1,p,q,fuzz2,single); + return S1.size() > 0; + } + qscale=qoffset=0.5; + } else { + Int tq=lq/2; + q1=q.subpath(0,tq); + q2=q.subpath(tq,lq); + qoffset=tq; + qscale=1.0; + } + + bool Short=lp == 1 && lq == 1; + + static size_t maxcount=9; + size_t count=0; + + if(intersections(s,t,S1,T1,p1,q1,fuzz,single,exact,depth)) { + add(s,t,S,T,S1,T1,pscale,qscale,0.0,0.0,p,q,fuzz2,single); + if(single || depth <= mindepth) + return true; + count += S1.size(); + if(Short && count > maxcount) return true; + } + + S1.clear(); + T1.clear(); + if(intersections(s,t,S1,T1,p1,q2,fuzz,single,exact,depth)) { + add(s,t,S,T,S1,T1,pscale,qscale,0.0,qoffset,p,q,fuzz2,single); + if(single || depth <= mindepth) + return true; + count += S1.size(); + if(Short && count > maxcount) return true; + } + + S1.clear(); + T1.clear(); + if(intersections(s,t,S1,T1,p2,q1,fuzz,single,exact,depth)) { + add(s,t,S,T,S1,T1,pscale,qscale,poffset,0.0,p,q,fuzz2,single); + if(single || depth <= mindepth) + return true; + count += S1.size(); + if(Short && count > maxcount) return true; + } + + S1.clear(); + T1.clear(); + if(intersections(s,t,S1,T1,p2,q2,fuzz,single,exact,depth)) { + add(s,t,S,T,S1,T1,pscale,qscale,poffset,qoffset,p,q,fuzz2,single); + if(single || depth <= mindepth) + return true; + count += S1.size(); + if(Short && count > maxcount) return true; + } + + return S.size() > 0; + } + return false; +} + +// }}} + +path3 concat(const path3& p1, const path3& p2) +{ + Int n1 = p1.length(), n2 = p2.length(); + + if (n1 == -1) return p2; + if (n2 == -1) return p1; + triple a=p1.point(n1); + triple b=p2.point((Int) 0); + + mem::vector<solvedKnot3> nodes(n1+n2+1); + + Int i = 0; + nodes[0].pre = p1.point((Int) 0); + for (Int j = 0; j < n1; j++) { + nodes[i].point = p1.point(j); + nodes[i].straight = p1.straight(j); + nodes[i].post = p1.postcontrol(j); + nodes[i+1].pre = p1.precontrol(j+1); + i++; + } + for (Int j = 0; j < n2; j++) { + nodes[i].point = p2.point(j); + nodes[i].straight = p2.straight(j); + nodes[i].post = p2.postcontrol(j); + nodes[i+1].pre = p2.precontrol(j+1); + i++; + } + nodes[i].point = nodes[i].post = p2.point(n2); + + return path3(nodes, i+1); +} + +path3 transformed(const array& t, const path3& p) +{ + Int n = p.size(); + mem::vector<solvedKnot3> nodes(n); + + for (Int i = 0; i < n; ++i) { + nodes[i].pre = t * p.precontrol(i); + nodes[i].point = t * p.point(i); + nodes[i].post = t * p.postcontrol(i); + nodes[i].straight = p.straight(i); + } + + return path3(nodes, n, p.cyclic()); +} + +path3 transformed(const double* t, const path3& p) +{ + Int n = p.size(); + mem::vector<solvedKnot3> nodes(n); + + for(Int i=0; i < n; ++i) { + nodes[i].pre=t*p.precontrol(i); + nodes[i].point=t*p.point(i); + nodes[i].post=t*p.postcontrol(i); + nodes[i].straight=p.straight(i); + } + + return path3(nodes, n, p.cyclic()); +} + +template<class T> +struct Split { + T m0,m1,m2,m3,m4,m5; + Split(T z0, T c0, T c1, T z1) { + m0=0.5*(z0+c0); + m1=0.5*(c0+c1); + m2=0.5*(c1+z1); + m3=0.5*(m0+m1); + m4=0.5*(m1+m2); + m5=0.5*(m3+m4); + } +}; + +double cornerbound(double *P, double (*m)(double, double)) +{ + double b=m(P[0],P[3]); + b=m(b,P[12]); + return m(b,P[15]); +} + +double controlbound(double *P, double (*m)(double, double)) +{ + double b=m(P[1],P[2]); + b=m(b,P[4]); + b=m(b,P[5]); + b=m(b,P[6]); + b=m(b,P[7]); + b=m(b,P[8]); + b=m(b,P[9]); + b=m(b,P[10]); + b=m(b,P[11]); + b=m(b,P[13]); + return m(b,P[14]); +} + +double bound(double *P, double (*m)(double, double), double b, + double fuzz, int depth) +{ + b=m(b,cornerbound(P,m)); + if(m(-1.0,1.0)*(b-controlbound(P,m)) >= -fuzz || depth == 0) + return b; + + --depth; + fuzz *= 2; + + Split<double> c0(P[0],P[1],P[2],P[3]); + Split<double> c1(P[4],P[5],P[6],P[7]); + Split<double> c2(P[8],P[9],P[10],P[11]); + Split<double> c3(P[12],P[13],P[14],P[15]); + + Split<double> c4(P[12],P[8],P[4],P[0]); + Split<double> c5(c3.m0,c2.m0,c1.m0,c0.m0); + Split<double> c6(c3.m3,c2.m3,c1.m3,c0.m3); + Split<double> c7(c3.m5,c2.m5,c1.m5,c0.m5); + Split<double> c8(c3.m4,c2.m4,c1.m4,c0.m4); + Split<double> c9(c3.m2,c2.m2,c1.m2,c0.m2); + Split<double> c10(P[15],P[11],P[7],P[3]); + + // Check all 4 Bezier subpatches. + double s0[]={c4.m5,c5.m5,c6.m5,c7.m5,c4.m3,c5.m3,c6.m3,c7.m3, + c4.m0,c5.m0,c6.m0,c7.m0,P[12],c3.m0,c3.m3,c3.m5}; + b=bound(s0,m,b,fuzz,depth); + double s1[]={P[0],c0.m0,c0.m3,c0.m5,c4.m2,c5.m2,c6.m2,c7.m2, + c4.m4,c5.m4,c6.m4,c7.m4,c4.m5,c5.m5,c6.m5,c7.m5}; + b=bound(s1,m,b,fuzz,depth); + double s2[]={c0.m5,c0.m4,c0.m2,P[3],c7.m2,c8.m2,c9.m2,c10.m2, + c7.m4,c8.m4,c9.m4,c10.m4,c7.m5,c8.m5,c9.m5,c10.m5}; + b=bound(s2,m,b,fuzz,depth); + double s3[]={c7.m5,c8.m5,c9.m5,c10.m5,c7.m3,c8.m3,c9.m3,c10.m3, + c7.m0,c8.m0,c9.m0,c10.m0,c3.m5,c3.m4,c3.m2,P[15]}; + return bound(s3,m,b,fuzz,depth); +} + +double cornerbound(triple *P, double (*m)(double, double), + double (*f)(const triple&)) +{ + double b=m(f(P[0]),f(P[3])); + b=m(b,f(P[12])); + return m(b,f(P[15])); +} + +// Return f evaluated at controlling vertex of bounding box of n control +// net points for similiar-triangle transform x'=x/z, y'=y/z, where z < 0. +double ratiobound(triple *P, double (*m)(double, double), + double (*f)(const triple&), int n) +{ + double MX=-P[0].getx(); + double MY=-P[0].gety(); + double Z=P[0].getz(); + double MZ=-Z; + for(int i=1; i < n; ++i) { + triple v=P[i]; + MX=m(MX,-v.getx()); + MY=m(MY,-v.gety()); + Z=m(Z,v.getz()); + MZ=m(MZ,-v.getz()); + } + return m(f(triple(-MX,-MY,Z)),f(triple(-MX,-MY,-MZ))); +} + +double controlbound(triple *P, double (*m)(double, double), + double (*f)(const triple&)) +{ + double b=m(f(P[1]),f(P[2])); + b=m(b,f(P[4])); + b=m(b,f(P[5])); + b=m(b,f(P[6])); + b=m(b,f(P[7])); + b=m(b,f(P[8])); + b=m(b,f(P[9])); + b=m(b,f(P[10])); + b=m(b,f(P[11])); + b=m(b,f(P[13])); + return m(b,f(P[14])); +} + +double bound(triple *P, double (*m)(double, double), + double (*f)(const triple&), double b, double fuzz, int depth) +{ + b=m(b,cornerbound(P,m,f)); + if(m(-1.0,1.0)*(b-ratiobound(P,m,f,16)) >= -fuzz || depth == 0) + return b; + + --depth; + fuzz *= 2; + + Split<triple> c0(P[0],P[1],P[2],P[3]); + Split<triple> c1(P[4],P[5],P[6],P[7]); + Split<triple> c2(P[8],P[9],P[10],P[11]); + Split<triple> c3(P[12],P[13],P[14],P[15]); + + Split<triple> c4(P[12],P[8],P[4],P[0]); + Split<triple> c5(c3.m0,c2.m0,c1.m0,c0.m0); + Split<triple> c6(c3.m3,c2.m3,c1.m3,c0.m3); + Split<triple> c7(c3.m5,c2.m5,c1.m5,c0.m5); + Split<triple> c8(c3.m4,c2.m4,c1.m4,c0.m4); + Split<triple> c9(c3.m2,c2.m2,c1.m2,c0.m2); + Split<triple> c10(P[15],P[11],P[7],P[3]); + + // Check all 4 Bezier subpatches. + triple s0[]={c4.m5,c5.m5,c6.m5,c7.m5,c4.m3,c5.m3,c6.m3,c7.m3, + c4.m0,c5.m0,c6.m0,c7.m0,P[12],c3.m0,c3.m3,c3.m5}; + b=bound(s0,m,f,b,fuzz,depth); + triple s1[]={P[0],c0.m0,c0.m3,c0.m5,c4.m2,c5.m2,c6.m2,c7.m2, + c4.m4,c5.m4,c6.m4,c7.m4,c4.m5,c5.m5,c6.m5,c7.m5}; + b=bound(s1,m,f,b,fuzz,depth); + triple s2[]={c0.m5,c0.m4,c0.m2,P[3],c7.m2,c8.m2,c9.m2,c10.m2, + c7.m4,c8.m4,c9.m4,c10.m4,c7.m5,c8.m5,c9.m5,c10.m5}; + b=bound(s2,m,f,b,fuzz,depth); + triple s3[]={c7.m5,c8.m5,c9.m5,c10.m5,c7.m3,c8.m3,c9.m3,c10.m3, + c7.m0,c8.m0,c9.m0,c10.m0,c3.m5,c3.m4,c3.m2,P[15]}; + return bound(s3,m,f,b,fuzz,depth); +} + +template<class T> +struct Splittri { + T l003,p102,p012,p201,p111,p021,r300,p210,p120,u030; + T u021,u120; + T p033,p231,p330; + T p123; + T l012,p312,r210,l102,p303,r201; + T u012,u210,l021,p4xx,r120,px4x,pxx4,l201,r102; + T l210,r012,l300; + T r021,u201,r030; + T u102,l120,l030; + T l111,r111,u111,c111; + + Splittri(const T *p) { + l003=p[0]; + p102=p[1]; + p012=p[2]; + p201=p[3]; + p111=p[4]; + p021=p[5]; + r300=p[6]; + p210=p[7]; + p120=p[8]; + u030=p[9]; + + u021=0.5*(u030+p021); + u120=0.5*(u030+p120); + + p033=0.5*(p021+p012); + p231=0.5*(p120+p111); + p330=0.5*(p120+p210); + + p123=0.5*(p012+p111); + + l012=0.5*(p012+l003); + p312=0.5*(p111+p201); + r210=0.5*(p210+r300); + + l102=0.5*(l003+p102); + p303=0.5*(p102+p201); + r201=0.5*(p201+r300); + + u012=0.5*(u021+p033); + u210=0.5*(u120+p330); + l021=0.5*(p033+l012); + p4xx=0.5*p231+0.25*(p111+p102); + r120=0.5*(p330+r210); + px4x=0.5*p123+0.25*(p111+p210); + pxx4=0.25*(p021+p111)+0.5*p312; + l201=0.5*(l102+p303); + r102=0.5*(p303+r201); + + l210=0.5*(px4x+l201); // = m120 + r012=0.5*(px4x+r102); // = m021 + l300=0.5*(l201+r102); // = r003 = m030 + + r021=0.5*(pxx4+r120); // = m012 + u201=0.5*(u210+pxx4); // = m102 + r030=0.5*(u210+r120); // = u300 = m003 + + u102=0.5*(u012+p4xx); // = m201 + l120=0.5*(l021+p4xx); // = m210 + l030=0.5*(u012+l021); // = u003 = m300 + + l111=0.5*(p123+l102); + r111=0.5*(p312+r210); + u111=0.5*(u021+p231); + c111=0.25*(p033+p330+p303+p111); + } +}; + +// Return the extremum of the vertices of a Bezier triangle. +double cornerboundtri(double *P, double (*m)(double, double)) +{ + double b=m(P[0],P[6]); + return m(b,P[9]); +} + +double cornerboundtri(triple *P, double (*m)(double, double), + double (*f)(const triple&)) +{ + double b=m(f(P[0]),f(P[6])); + return m(b,f(P[9])); +} + +// Return the extremum of the non-vertex control points of a Bezier triangle. +double controlboundtri(double *P, double (*m)(double, double)) +{ + double b=m(P[1],P[2]); + b=m(b,P[3]); + b=m(b,P[4]); + b=m(b,P[5]); + b=m(b,P[7]); + return m(b,P[8]); +} + +double controlboundtri(triple *P, double (*m)(double, double), + double (*f)(const triple&)) +{ + double b=m(f(P[1]),f(P[2])); + b=m(b,f(P[3])); + b=m(b,f(P[4])); + b=m(b,f(P[5])); + b=m(b,f(P[7])); + return m(b,f(P[8])); +} + +// Return the global bound of a Bezier triangle. +double boundtri(double *P, double (*m)(double, double), double b, + double fuzz, int depth) +{ + b=m(b,cornerboundtri(P,m)); + if(m(-1.0,1.0)*(b-controlboundtri(P,m)) >= -fuzz || depth == 0) + return b; + + --depth; + fuzz *= 2; + + Splittri<double> s(P); + + double l[]={s.l003,s.l102,s.l012,s.l201,s.l111, + s.l021,s.l300,s.l210,s.l120,s.l030}; // left + b=boundtri(l,m,b,fuzz,depth); + + double r[]={s.l300,s.r102,s.r012,s.r201,s.r111, + s.r021,s.r300,s.r210,s.r120,s.r030}; // right + b=boundtri(r,m,b,fuzz,depth); + + double u[]={s.l030,s.u102,s.u012,s.u201,s.u111, + s.u021,s.r030,s.u210,s.u120,s.u030}; // up + b=boundtri(u,m,b,fuzz,depth); + + double c[]={s.r030,s.u201,s.r021,s.u102,s.c111, + s.r012,s.l030,s.l120,s.l210,s.l300}; // center + return boundtri(c,m,b,fuzz,depth); +} + +double boundtri(triple *P, double (*m)(double, double), + double (*f)(const triple&), double b, double fuzz, int depth) +{ + b=m(b,cornerboundtri(P,m,f)); + if(m(-1.0,1.0)*(b-ratiobound(P,m,f,10)) >= -fuzz || depth == 0) + return b; + + --depth; + fuzz *= 2; + + Splittri<triple> s(P); + + triple l[]={s.l003,s.l102,s.l012,s.l201,s.l111, + s.l021,s.l300,s.l210,s.l120,s.l030}; // left + b=boundtri(l,m,f,b,fuzz,depth); + + triple r[]={s.l300,s.r102,s.r012,s.r201,s.r111, + s.r021,s.r300,s.r210,s.r120,s.r030}; // right + b=boundtri(r,m,f,b,fuzz,depth); + + triple u[]={s.l030,s.u102,s.u012,s.u201,s.u111, + s.u021,s.r030,s.u210,s.u120,s.u030}; // up + b=boundtri(u,m,f,b,fuzz,depth); + + triple c[]={s.r030,s.u201,s.r021,s.u102,s.c111, + s.r012,s.l030,s.l120,s.l210,s.l300}; // center + return boundtri(c,m,f,b,fuzz,depth); +} + +inline void add(std::vector<double>& T, std::vector<double>& U, + std::vector<double>& V, double t, double u, double v, + const path3& p, double fuzz2) +{ + triple z=p.point(t); + size_t n=T.size(); + for(size_t i=0; i < n; ++i) + if((p.point(T[i])-z).abs2() <= fuzz2) return; + T.push_back(t); + U.push_back(u); + V.push_back(v); +} + +void add(std::vector<double>& T, std::vector<double>& U, + std::vector<double>& V, std::vector<double>& T1, + std::vector<double>& U1, std::vector<double>& V1, + const path3& p, double tscale, double toffset, + double uoffset, double voffset, double fuzz2) +{ + size_t n=T1.size(); + for(size_t i=0; i < n; ++i) + add(T,U,V,tscale*T1[i]+toffset,0.5*U1[i]+uoffset,0.5*V1[i]+voffset,p, + fuzz2); +} + +void bounds(triple& Pmin, triple& Pmax, triple *P, double fuzz) +{ + double Px[]={P[0].getx(),P[1].getx(),P[2].getx(),P[3].getx(), + P[4].getx(),P[5].getx(),P[6].getx(),P[7].getx(), + P[8].getx(),P[9].getx(),P[10].getx(),P[11].getx(), + P[12].getx(),P[13].getx(),P[14].getx(),P[15].getx()}; + double bx=Px[0]; + double xmin=bound(Px,min,bx,fuzz,maxdepth); + double xmax=bound(Px,max,bx,fuzz,maxdepth); + + double Py[]={P[0].gety(),P[1].gety(),P[2].gety(),P[3].gety(), + P[4].gety(),P[5].gety(),P[6].gety(),P[7].gety(), + P[8].gety(),P[9].gety(),P[10].gety(),P[11].gety(), + P[12].gety(),P[13].gety(),P[14].gety(),P[15].gety()}; + double by=Py[0]; + double ymin=bound(Py,min,by,fuzz,maxdepth); + double ymax=bound(Py,max,by,fuzz,maxdepth); + + double Pz[]={P[0].getz(),P[1].getz(),P[2].getz(),P[3].getz(), + P[4].getz(),P[5].getz(),P[6].getz(),P[7].getz(), + P[8].getz(),P[9].getz(),P[10].getz(),P[11].getz(), + P[12].getz(),P[13].getz(),P[14].getz(),P[15].getz()}; + double bz=Pz[0]; + double zmin=bound(Pz,min,bz,fuzz,maxdepth); + double zmax=bound(Pz,max,bz,fuzz,maxdepth); + Pmin=triple(xmin,ymin,zmin); + Pmax=triple(xmax,ymax,zmax); +} + +inline double abs2(double x, double y, double z) +{ + return x*x+y*y+z*z; +} + +bool intersections(double& U, double& V, const triple& v, triple *P, + double fuzz, unsigned depth) +{ + if(errorstream::interrupt) throw interrupted(); + + triple Pmin,Pmax; + bounds(Pmin,Pmax,P,fuzz); + + double x=P[0].getx(); + double y=P[0].gety(); + double z=P[0].getz(); + double X=x, Y=y, Z=z; + for(int i=1; i < 16; ++i) { + triple v=P[i]; + double vx=v.getx(); + x=min(x,vx); + X=max(X,vx); + double vy=v.gety(); + y=min(y,vy); + Y=max(Y,vy); + double vz=v.getz(); + z=min(z,vz); + Z=max(Z,vz); + } + + if(X+fuzz >= v.getx() && + Y+fuzz >= v.gety() && + Z+fuzz >= v.getz() && + v.getx()+fuzz >= x && + v.gety()+fuzz >= y && + v.getz()+fuzz >= z) { // Overlapping bounding boxes + + --depth; +// fuzz *= 2; + + if(abs2(X-x,Y-y,Z-z) <= fuzz*fuzz || depth == 0) { + U=0.5; + V=0.5; + return true; + } + +// Compute the control points of the four subpatches obtained by splitting +// the patch with control points P at u=v=1/2. + Split<triple> c0(P[0],P[1],P[2],P[3]); + Split<triple> c1(P[4],P[5],P[6],P[7]); + Split<triple> c2(P[8],P[9],P[10],P[11]); + Split<triple> c3(P[12],P[13],P[14],P[15]); + + Split<triple> c4(P[12],P[8],P[4],P[0]); + Split<triple> c5(c3.m0,c2.m0,c1.m0,c0.m0); + Split<triple> c6(c3.m3,c2.m3,c1.m3,c0.m3); + Split<triple> c7(c3.m5,c2.m5,c1.m5,c0.m5); + Split<triple> c8(c3.m4,c2.m4,c1.m4,c0.m4); + Split<triple> c9(c3.m2,c2.m2,c1.m2,c0.m2); + Split<triple> c10(P[15],P[11],P[7],P[3]); + + // Check all 4 Bezier subpatches. + + double U1,V1; + triple Q0[]={P[0],c0.m0,c0.m3,c0.m5,c4.m2,c5.m2,c6.m2,c7.m2, + c4.m4,c5.m4,c6.m4,c7.m4,c4.m5,c5.m5,c6.m5,c7.m5}; + if(intersections(U1,V1,v,Q0,fuzz,depth)) { + U=0.5*U1; + V=0.5*V1; + return true; + } + + triple Q1[]={c0.m5,c0.m4,c0.m2,P[3],c7.m2,c8.m2,c9.m2,c10.m2, + c7.m4,c8.m4,c9.m4,c10.m4,c7.m5,c8.m5,c9.m5,c10.m5}; + if(intersections(U1,V1,v,Q1,fuzz,depth)) { + U=0.5*U1; + V=0.5*V1+0.5; + return true; + } + + triple Q2[]={c7.m5,c8.m5,c9.m5,c10.m5,c7.m3,c8.m3,c9.m3,c10.m3, + c7.m0,c8.m0,c9.m0,c10.m0,c3.m5,c3.m4,c3.m2,P[15]}; + if(intersections(U1,V1,v,Q2,fuzz,depth)) { + U=0.5*U1+0.5; + V=0.5*V1+0.5; + return true; + } + + triple Q3[]={c4.m5,c5.m5,c6.m5,c7.m5,c4.m3,c5.m3,c6.m3,c7.m3, + c4.m0,c5.m0,c6.m0,c7.m0,P[12],c3.m0,c3.m3,c3.m5}; + if(intersections(U1,V1,v,Q3,fuzz,depth)) { + U=0.5*U1+0.5; + V=0.5*V1; + return true; + } + } + return false; +} + +bool intersections(std::vector<double>& T, std::vector<double>& U, + std::vector<double>& V, path3& p, triple *P, + double fuzz, bool single, unsigned depth) +{ + if(errorstream::interrupt) throw interrupted(); + + double fuzz2=max(fuzzFactor*fuzz*fuzz,Fuzz2); + + triple pmin=p.min(); + triple pmax=p.max(); + + double x=P[0].getx(); + double y=P[0].gety(); + double z=P[0].getz(); + double X=x, Y=y, Z=z; + for(int i=1; i < 16; ++i) { + triple v=P[i]; + double vx=v.getx(); + x=min(x,vx); + X=max(X,vx); + double vy=v.gety(); + y=min(y,vy); + Y=max(Y,vy); + double vz=v.getz(); + z=min(z,vz); + Z=max(Z,vz); + } + + if(X+fuzz >= pmin.getx() && + Y+fuzz >= pmin.gety() && + Z+fuzz >= pmin.getz() && + pmax.getx()+fuzz >= x && + pmax.gety()+fuzz >= y && + pmax.getz()+fuzz >= z) { // Overlapping bounding boxes + + --depth; +// fuzz *= 2; + + if(((pmax-pmin).length()+sqrt(abs2(X-x,Y-y,Z-z)) <= fuzz) || depth == 0) { + T.push_back(0.5); + U.push_back(0.5); + V.push_back(0.5); + return true; + } + + Int lp=p.length(); + + path3 p0,p1; + p.halve(p0,p1); + + std::vector<double> T1,U1,V1; + double tscale,toffset; + + if(lp <= 1) { + if(lp == 1) p.halve(p0,p1); + if(lp == 0 || p0 == p || p1 == p) { + double u,v; + if(intersections(u,v,p.point((Int) 0),P,fuzz,depth)) { + T1.push_back(0.0); + U1.push_back(u); + V1.push_back(v); + add(T,U,V,T1,U1,V1,p,1.0,0.0,0.0,0.0,fuzz2); + } + return T1.size() > 0; + } + tscale=toffset=0.5; + } else { + Int tp=lp/2; + p0=p.subpath(0,tp); + p1=p.subpath(tp,lp); + toffset=tp; + tscale=1.0; + } + + Split<triple> c0(P[0],P[1],P[2],P[3]); + Split<triple> c1(P[4],P[5],P[6],P[7]); + Split<triple> c2(P[8],P[9],P[10],P[11]); + Split<triple> c3(P[12],P[13],P[14],P[15]); + + Split<triple> c4(P[12],P[8],P[4],P[0]); + Split<triple> c5(c3.m0,c2.m0,c1.m0,c0.m0); + Split<triple> c6(c3.m3,c2.m3,c1.m3,c0.m3); + Split<triple> c7(c3.m5,c2.m5,c1.m5,c0.m5); + Split<triple> c8(c3.m4,c2.m4,c1.m4,c0.m4); + Split<triple> c9(c3.m2,c2.m2,c1.m2,c0.m2); + Split<triple> c10(P[15],P[11],P[7],P[3]); + + static size_t maxcount=9; + size_t count=0; + + bool Short=lp == 1; + + // Check all 4 Bezier subpatches against p0. + triple Q0[]={P[0],c0.m0,c0.m3,c0.m5,c4.m2,c5.m2,c6.m2,c7.m2, + c4.m4,c5.m4,c6.m4,c7.m4,c4.m5,c5.m5,c6.m5,c7.m5}; + if(intersections(T1,U1,V1,p0,Q0,fuzz,single,depth)) { + add(T,U,V,T1,U1,V1,p,tscale,0.0,0.0,0.0,fuzz2); + if(single || depth <= mindepth) + return true; + count += T1.size(); + if(Short && count > maxcount) return true; + } + + T1.clear(); + U1.clear(); + V1.clear(); + triple Q1[]={c0.m5,c0.m4,c0.m2,P[3],c7.m2,c8.m2,c9.m2,c10.m2, + c7.m4,c8.m4,c9.m4,c10.m4,c7.m5,c8.m5,c9.m5,c10.m5}; + if(intersections(T1,U1,V1,p0,Q1,fuzz,single,depth)) { + add(T,U,V,T1,U1,V1,p,tscale,0.0,0.0,0.5,fuzz2); + if(single || depth <= mindepth) + return true; + count += T1.size(); + if(Short && count > maxcount) return true; + } + + T1.clear(); + U1.clear(); + V1.clear(); + triple Q2[]={c7.m5,c8.m5,c9.m5,c10.m5,c7.m3,c8.m3,c9.m3,c10.m3, + c7.m0,c8.m0,c9.m0,c10.m0,c3.m5,c3.m4,c3.m2,P[15]}; + if(intersections(T1,U1,V1,p0,Q2,fuzz,single,depth)) { + add(T,U,V,T1,U1,V1,p,tscale,0.0,0.5,0.5,fuzz2); + if(single || depth <= mindepth) + return true; + count += T1.size(); + if(Short && count > maxcount) return true; + } + + T1.clear(); + U1.clear(); + V1.clear(); + triple Q3[]={c4.m5,c5.m5,c6.m5,c7.m5,c4.m3,c5.m3,c6.m3,c7.m3, + c4.m0,c5.m0,c6.m0,c7.m0,P[12],c3.m0,c3.m3,c3.m5}; + if(intersections(T1,U1,V1,p0,Q3,fuzz,single,depth)) { + add(T,U,V,T1,U1,V1,p,tscale,0.0,0.5,0.0,fuzz2); + if(single || depth <= mindepth) + return true; + count += T1.size(); + if(Short && count > maxcount) return true; + } + + // Check all 4 Bezier subpatches against p1. + T1.clear(); + U1.clear(); + V1.clear(); + if(intersections(T1,U1,V1,p1,Q0,fuzz,single,depth)) { + add(T,U,V,T1,U1,V1,p,tscale,toffset,0.0,0.0,fuzz2); + if(single || depth <= mindepth) + return true; + count += T1.size(); + if(Short && count > maxcount) return true; + } + + T1.clear(); + U1.clear(); + V1.clear(); + if(intersections(T1,U1,V1,p1,Q1,fuzz,single,depth)) { + add(T,U,V,T1,U1,V1,p,tscale,toffset,0.0,0.5,fuzz2); + if(single || depth <= mindepth) + return true; + count += T1.size(); + if(Short && count > maxcount) return true; + } + + T1.clear(); + U1.clear(); + V1.clear(); + if(intersections(T1,U1,V1,p1,Q2,fuzz,single,depth)) { + add(T,U,V,T1,U1,V1,p,tscale,toffset,0.5,0.5,fuzz2); + if(single || depth <= mindepth) + return true; + count += T1.size(); + if(Short && count > maxcount) return true; + } + + T1.clear(); + U1.clear(); + V1.clear(); + if(intersections(T1,U1,V1,p1,Q3,fuzz,single,depth)) { + add(T,U,V,T1,U1,V1,p,tscale,toffset,0.5,0.0,fuzz2); + if(single || depth <= mindepth) + return true; + count += T1.size(); + if(Short && count > maxcount) return true; + } + + return T.size() > 0; + } + return false; +} + +} //namespace camp |