summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/path3.cc
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/utils/asymptote/path3.cc')
-rw-r--r--Build/source/utils/asymptote/path3.cc1460
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