diff options
Diffstat (limited to 'Build/source/utils/asymptote/path3.cc')
-rw-r--r-- | Build/source/utils/asymptote/path3.cc | 505 |
1 files changed, 383 insertions, 122 deletions
diff --git a/Build/source/utils/asymptote/path3.cc b/Build/source/utils/asymptote/path3.cc index 7cfa16647d5..da973dc36be 100644 --- a/Build/source/utils/asymptote/path3.cc +++ b/Build/source/utils/asymptote/path3.cc @@ -24,28 +24,6 @@ void checkEmpty3(Int n) { reportError("nullpath3 has no points"); } -double bound(triple z0, triple c0, triple c1, triple z1, - double (*m)(double, double), - double (*f)(const triple&, double*), double *t, - double b, int depth) -{ - b=m(b,m(f(z0,t),f(z1,t))); - if(m(-1.0,1.0)*(b-m(f(c0,t),f(c1,t))) >= -sqrtFuzz || depth == 0) - return b; - --depth; - - 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,t,b,depth); - return bound(m5,m4,m2,z1,m,f,t,b,depth); -} - triple path3::point(double t) const { checkEmpty3(n); @@ -337,6 +315,28 @@ bbox3 path3::bounds() const return box; } +double bound(triple z0, triple c0, triple c1, triple z1, + double (*m)(double, double), + double (*f)(const triple&, double*), double *t, + double b, double fuzz, int depth) +{ + b=m(b,m(f(z0,t),f(z1,t))); + if(m(-1.0,1.0)*(b-m(f(c0,t),f(c1,t))) >= -fuzz || depth == 0) + return b; + --depth; + + 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,t,b,fuzz,depth); + return bound(m5,m4,m2,z1,m,f,t,b,fuzz,depth); +} + pair path3::bounds(double (*m)(double, double), double (*x)(const triple&, double*), double (*y)(const triple&, double*), double *t) const @@ -515,7 +515,7 @@ void intersections(std::vector<double>& S, std::vector<double>& T, size_t n=S1.size(); for(size_t i=0; i < n; ++i) { S.push_back(S1[i]); - T.push_back(0); + T.push_back(0.0); } } @@ -532,13 +532,12 @@ void add(std::vector<double>& S, std::vector<double>& T, double s, double 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 fuzz, bool single) + const path3& p, const path3& q, double fuzz2, bool single) { if(single) { s=s*pscale+poffset; t=t*qscale+qoffset; } else { - double fuzz2=4.0*fuzz*fuzz; 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); @@ -547,7 +546,7 @@ void add(double& s, double& t, std::vector<double>& S, std::vector<double>& T, 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 fuzz, bool single) + const path3& p, const path3& q, double fuzz2, bool single) { size_t n=S1.size(); if(single) { @@ -556,7 +555,6 @@ void add(double& s, double& t, std::vector<double>& S, std::vector<double>& T, t=T1[0]; } } else { - double fuzz2=4.0*fuzz*fuzz; for(size_t i=0; i < n; ++i) add(S,T,S1[i],T1[i],p,q,fuzz2); } @@ -568,11 +566,14 @@ bool intersections(double &s, double &t, std::vector<double>& S, { if(errorstream::interrupt) throw interrupted(); + double fuzz2=max(fuzzFactor*fuzz,Fuzz); + fuzz2=fuzz2*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,fuzz,single); + add(s,t,S,T,S1,T1,p,q,fuzz2,single); return S1.size() > 0; } @@ -580,7 +581,7 @@ bool intersections(double &s, double &t, std::vector<double>& S, 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,fuzz,single); + add(s,t,S,T,S1,T1,p,q,fuzz2,single); return S1.size() > 0; } @@ -600,11 +601,11 @@ bool intersections(double &s, double &t, std::vector<double>& S, --depth; if((maxp-minp).length()+(maxq-minq).length() <= fuzz || depth == 0) { if(single) { - s=0; - t=0; + s=0.5; + t=0.5; } else { - S.push_back(0.0); - T.push_back(0.0); + S.push_back(0.5); + T.push_back(0.5); } return true; } @@ -612,12 +613,13 @@ bool intersections(double &s, double &t, std::vector<double>& S, path3 p1,p2; double pscale,poffset; + std::vector<double> S1,T1; + if(lp <= 1) { if(lp == 1) p.halve(p1,p2); if(lp == 0 || p1 == p || p2 == p) { - std::vector<double> T1,S1; intersections(T1,S1,q,p.point((Int) 0),fuzz); - add(s,t,S,T,S1,T1,p,q,fuzz,single); + add(s,t,S,T,S1,T1,p,q,fuzz2,single); return S1.size() > 0; } pscale=poffset=0.5; @@ -635,9 +637,8 @@ bool intersections(double &s, double &t, std::vector<double>& S, if(lq <= 1) { if(lq == 1) q.halve(q1,q2); if(lq == 0 || q1 == q || q2 == q) { - std::vector<double> S1,T1; intersections(S1,T1,p,q.point((Int) 0),fuzz); - add(s,t,S,T,S1,T1,p,q,fuzz,single); + add(s,t,S,T,S1,T1,p,q,fuzz2,single); return S1.size() > 0; } qscale=qoffset=0.5; @@ -654,9 +655,8 @@ bool intersections(double &s, double &t, std::vector<double>& S, static size_t maxcount=9; size_t count=0; - std::vector<double> S1,T1; 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,fuzz,single); + 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(); @@ -666,7 +666,7 @@ bool intersections(double &s, double &t, std::vector<double>& S, 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,fuzz,single); + 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(); @@ -676,7 +676,7 @@ bool intersections(double &s, double &t, std::vector<double>& S, 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,fuzz,single); + 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(); @@ -686,7 +686,7 @@ bool intersections(double &s, double &t, std::vector<double>& S, 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,fuzz,single); + add(s,t,S,T,S1,T1,pscale,qscale,poffset,qoffset,p,q,fuzz2,single); if(single || depth <= mindepth) return true; count += S1.size(); @@ -779,9 +779,10 @@ double yratio(const triple& v, double *) return v.gety()/z; } +template<class T> struct Split { - double m0,m1,m2,m3,m4,m5; - Split(double z0, double c0, double c1, double z1) { + 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); @@ -791,123 +792,383 @@ struct Split { } }; -double cornerbound(double *p, double (*m)(double, double)) +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 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 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 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 cornerbound(triple *p, double (*m)(double, double), +double cornerbound(triple *P, double (*m)(double, double), double (*f)(const triple&, double*), double *t) { - double b=m(f(p[0],t),f(p[3],t)); - b=m(b,f(p[12],t)); - return m(b,f(p[15],t)); + double b=m(f(P[0],t),f(P[3],t)); + b=m(b,f(P[12],t)); + return m(b,f(P[15],t)); } -double controlbound(triple *p, double (*m)(double, double), +double controlbound(triple *P, double (*m)(double, double), double (*f)(const triple&, double*), double *t) { - double b=m(f(p[1],t),f(p[2],t)); - b=m(b,f(p[4],t)); - b=m(b,f(p[5],t)); - b=m(b,f(p[6],t)); - b=m(b,f(p[7],t)); - b=m(b,f(p[8],t)); - b=m(b,f(p[9],t)); - b=m(b,f(p[10],t)); - b=m(b,f(p[11],t)); - b=m(b,f(p[13],t)); - return m(b,f(p[14],t)); + double b=m(f(P[1],t),f(P[2],t)); + b=m(b,f(P[4],t)); + b=m(b,f(P[5],t)); + b=m(b,f(P[6],t)); + b=m(b,f(P[7],t)); + b=m(b,f(P[8],t)); + b=m(b,f(P[9],t)); + b=m(b,f(P[10],t)); + b=m(b,f(P[11],t)); + b=m(b,f(P[13],t)); + return m(b,f(P[14],t)); } -double bound(double *p, double (*m)(double, double), double b, int depth) +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)) >= -sqrtFuzz || depth == 0) + b=m(b,cornerbound(P,m)); + if(m(-1.0,1.0)*(b-controlbound(P,m)) >= -fuzz || depth == 0) return b; --depth; - Split c0(p[0],p[1],p[2],p[3]); - Split c1(p[4],p[5],p[6],p[7]); - Split c2(p[8],p[9],p[10],p[11]); - Split c3(p[12],p[13],p[14],p[15]); + 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 c4(p[12],p[8],p[4],p[0]); - Split c5(c3.m0,c2.m0,c1.m0,c0.m0); - Split c6(c3.m3,c2.m3,c1.m3,c0.m3); - Split c7(c3.m5,c2.m5,c1.m5,c0.m5); - Split c8(c3.m4,c2.m4,c1.m4,c0.m4); - Split c9(c3.m5,c2.m5,c1.m5,c0.m5); - Split c10(p[15],p[11],p[7],p[3]); + 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,depth); - double s1[]={p[0],c0.m0,c0.m3,c0.m5,c4.m2,c5.m2,c6.m2,c7.m2, + 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,depth); - double s2[]={c0.m5,c0.m4,c0.m2,p[3],c7.m2,c8.m2,c9.m2,c10.m2, + 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,depth); + 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,depth); + c7.m0,c8.m0,c9.m0,c10.m0,c3.m5,c3.m4,c3.m2,P[15]}; + return bound(s3,m,b,fuzz,depth); } -double bound(triple *p, double (*m)(double, double), +double bound(triple *P, double (*m)(double, double), double (*f)(const triple&, double*), double *t, - double b, int depth) + double b, double fuzz, int depth) { - b=m(b,cornerbound(p,m,f,t)); - if(m(-1.0,1.0)*(b-controlbound(p,m,f,t)) >= -sqrtFuzz || depth == 0) + b=m(b,cornerbound(P,m,f,t)); + if(m(-1.0,1.0)*(b-controlbound(P,m,f,t)) >= -fuzz || depth == 0) return b; --depth; - Split3 c0(p[0],p[1],p[2],p[3]); - Split3 c1(p[4],p[5],p[6],p[7]); - Split3 c2(p[8],p[9],p[10],p[11]); - Split3 c3(p[12],p[13],p[14],p[15]); + 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]); - Split3 c4(p[12],p[8],p[4],p[0]); - Split3 c5(c3.m0,c2.m0,c1.m0,c0.m0); - Split3 c6(c3.m3,c2.m3,c1.m3,c0.m3); - Split3 c7(c3.m5,c2.m5,c1.m5,c0.m5); - Split3 c8(c3.m4,c2.m4,c1.m4,c0.m4); - Split3 c9(c3.m5,c2.m5,c1.m5,c0.m5); - Split3 c10(p[15],p[11],p[7],p[3]); + 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,t,b,depth); - triple s1[]={p[0],c0.m0,c0.m3,c0.m5,c4.m2,c5.m2,c6.m2,c7.m2, + c4.m0,c5.m0,c6.m0,c7.m0,P[12],c3.m0,c3.m3,c3.m5}; + b=bound(s0,m,f,t,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,t,b,depth); - triple s2[]={c0.m5,c0.m4,c0.m2,p[3],c7.m2,c8.m2,c9.m2,c10.m2, + b=bound(s1,m,f,t,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,t,b,depth); + b=bound(s2,m,f,t,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,t,b,depth); + c7.m0,c8.m0,c9.m0,c10.m0,c3.m5,c3.m4,c3.m2,P[15]}; + return bound(s3,m,f,t,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); + double xmax=bound(Px,max,bx,fuzz); + + 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); + double ymax=bound(Py,max,by,fuzz); + + 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); + double zmax=bound(Pz,max,bz,fuzz); + Pmin=triple(xmin,ymin,zmin); + Pmax=triple(xmax,ymax,zmax); +} + +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); + + if(Pmax.getx()+fuzz >= v.getx() && + Pmax.gety()+fuzz >= v.gety() && + Pmax.getz()+fuzz >= v.getz() && + v.getx()+fuzz >= Pmin.getx() && + v.gety()+fuzz >= Pmin.gety() && + v.getz()+fuzz >= Pmin.getz()) { // Overlapping bounding boxes + + --depth; + if((Pmax-Pmin).abs2() <= 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; +} + +void intersections(std::vector<double>& T, std::vector<double>& U, + std::vector<double>& V, path3& p, triple *P, + double fuzz, unsigned depth) +{ + if(errorstream::interrupt) throw interrupted(); + + triple pmin=p.min(); + triple pmax=p.max(); + + triple Pmin,Pmax; + bounds(Pmin,Pmax,P,fuzz); + + if(Pmax.getx()+fuzz >= pmin.getx() && + Pmax.gety()+fuzz >= pmin.gety() && + Pmax.getz()+fuzz >= pmin.getz() && + pmax.getx()+fuzz >= Pmin.getx() && + pmax.gety()+fuzz >= Pmin.gety() && + pmax.getz()+fuzz >= Pmin.getz()) { // Overlapping bounding boxes + + --depth; + + if(((pmax-pmin).length()+(Pmax-Pmin).length() <= fuzz) || depth == 0) { + T.push_back(0.5); + U.push_back(0.5); + V.push_back(0.5); + return; + } + + Int lp=p.length(); + + path3 p0,p1; + p.halve(p0,p1); + + std::vector<double> T1,U1,V1; + double tscale,toffset; + + double fuzz2=max(fuzzFactor*fuzz,Fuzz); + fuzz2=fuzz2*fuzz2; + + 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; + } + 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]); + + // 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}; + intersections(T1,U1,V1,p0,Q0,fuzz,depth); + add(T,U,V,T1,U1,V1,p,tscale,0.0,0.0,0.0,fuzz2); + + 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}; + intersections(T1,U1,V1,p0,Q1,fuzz,depth); + add(T,U,V,T1,U1,V1,p,tscale,0.0,0.0,0.5,fuzz2); + + 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]}; + intersections(T1,U1,V1,p0,Q2,fuzz,depth); + add(T,U,V,T1,U1,V1,p,tscale,0.0,0.5,0.5,fuzz2); + + 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}; + intersections(T1,U1,V1,p0,Q3,fuzz,depth); + add(T,U,V,T1,U1,V1,p,tscale,0.0,0.5,0.0,fuzz2); + + // Check all 4 Bezier subpatches against p1. + T1.clear(); + U1.clear(); + V1.clear(); + intersections(T1,U1,V1,p1,Q0,fuzz,depth); + add(T,U,V,T1,U1,V1,p,tscale,toffset,0.0,0.0,fuzz2); + + T1.clear(); + U1.clear(); + V1.clear(); + intersections(T1,U1,V1,p1,Q1,fuzz,depth); + add(T,U,V,T1,U1,V1,p,tscale,toffset,0.0,0.5,fuzz2); + + T1.clear(); + U1.clear(); + V1.clear(); + intersections(T1,U1,V1,p1,Q2,fuzz,depth); + add(T,U,V,T1,U1,V1,p,tscale,toffset,0.5,0.5,fuzz2); + + T1.clear(); + U1.clear(); + V1.clear(); + intersections(T1,U1,V1,p1,Q3,fuzz,depth); + add(T,U,V,T1,U1,V1,p,tscale,toffset,0.5,0.0,fuzz2); + } + } } //namespace camp |