diff options
author | Karl Berry <karl@freefriends.org> | 2009-05-16 00:19:13 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2009-05-16 00:19:13 +0000 |
commit | bab45528d65eaafe68a705dbb2a57075c7b7cbd8 (patch) | |
tree | 10b4ae2b5195c8dede153ab89359ec00f55f325f /Build/source/utils/asymptote/path3.cc | |
parent | 8643d90372e9c31e0f461c15c596b60a545bd7d3 (diff) |
asymptote 1.72 sources (not integrated into build yet)
git-svn-id: svn://tug.org/texlive/trunk@13110 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/utils/asymptote/path3.cc')
-rw-r--r-- | Build/source/utils/asymptote/path3.cc | 914 |
1 files changed, 914 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..7cfa16647d5 --- /dev/null +++ b/Build/source/utils/asymptote/path3.cc @@ -0,0 +1,914 @@ +/***** + * 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"); +} + +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); + + 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 > n-1) + b = n-1; + } + + 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 > n-1) + a = n-1; + } + 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)); + + 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; +} + +pair path3::bounds(double (*m)(double, double), + double (*x)(const triple&, double*), + double (*y)(const triple&, double*), double *t) const +{ + checkEmpty3(n); + + triple v=point((Int) 0); + pair B=pair(x(v,t),y(v,t)); + + Int n=length(); + for(Int i=0; i <= n; ++i) { + if(straight(i)) { + triple v=point(i); + B=pair(m(B.getx(),x(v,t)),m(B.gety(),y(v,t))); + } 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,x,t,B.getx()), + bound(z0,c0,c1,z1,m,y,t,B.gety())); + } + } + 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 using adaptive simpson integration. +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; + } + const triple& c0=postcontrol(i); + const triple& c1=precontrol(i+1); + + 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"); + 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 >= -Fuzz && t <= 1.0+Fuzz) { + double s=i+t; + if((g.point(s)-v).abs2() <= fuzz2) { + if(cycles && s >= n-Fuzz) 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); + } +} + +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 fuzz, 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); + } +} + +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) +{ + size_t n=S1.size(); + if(single) { + if(n > 0) { + s=S1[0]; + 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); + } +} + +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(); + + 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); + 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,fuzz,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; + if((maxp-minp).length()+(maxq-minq).length() <= fuzz || depth == 0) { + if(single) { + s=0; + t=0; + } else { + S.push_back(0.0); + T.push_back(0.0); + } + return true; + } + + path3 p1,p2; + double pscale,poffset; + + 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); + 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) { + 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); + 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; + + 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); + 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,fuzz,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,fuzz,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,fuzz,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()); +} + +double xproject(const triple& v, double *t) +{ + double x=v.getx(); + double y=v.gety(); + double z=v.getz(); + double f=t[12]*x+t[13]*y+t[14]*z+t[15]; + if(f == 0.0) run::dividebyzero(); + return (t[0]*x+t[1]*y+t[2]*z+t[3])/f; +} + +double yproject(const triple& v, double *t) +{ + double x=v.getx(); + double y=v.gety(); + double z=v.getz(); + double f=t[12]*x+t[13]*y+t[14]*z+t[15]; + if(f == 0.0) run::dividebyzero(); + return (t[4]*x+t[5]*y+t[6]*z+t[7])/f; +} + +double xratio(const triple& v, double *) +{ + double z=v.getz(); + return v.getx()/z; +} + +double yratio(const triple& v, double *) +{ + double z=v.getz(); + return v.gety()/z; +} + +struct Split { + double m0,m1,m2,m3,m4,m5; + Split(double z0, double c0, double c1, double 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 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 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 bound(double *p, double (*m)(double, double), double b, int depth) +{ + b=m(b,cornerbound(p,m)); + if(m(-1.0,1.0)*(b-controlbound(p,m)) >= -sqrtFuzz || 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 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]); + + // 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.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, + c7.m4,c8.m4,c9.m4,c10.m4,c7.m5,c8.m5,c9.m5,c10.m5}; + b=bound(s2,m,b,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); +} + +double bound(triple *p, double (*m)(double, double), + double (*f)(const triple&, double*), double *t, + double b, 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) + 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]); + + 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]); + + // 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.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, + c7.m4,c8.m4,c9.m4,c10.m4,c7.m5,c8.m5,c9.m5,c10.m5}; + b=bound(s2,m,f,t,b,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); +} + +} //namespace camp + |