diff options
Diffstat (limited to 'Build/source/utils/asymptote/path3.cc')
-rw-r--r-- | Build/source/utils/asymptote/path3.cc | 186 |
1 files changed, 174 insertions, 12 deletions
diff --git a/Build/source/utils/asymptote/path3.cc b/Build/source/utils/asymptote/path3.cc index 86029b6e60e..07f31ad934d 100644 --- a/Build/source/utils/asymptote/path3.cc +++ b/Build/source/utils/asymptote/path3.cc @@ -853,16 +853,16 @@ double cornerbound(triple *P, double (*m)(double, double), return m(b,f(P[15])); } -// 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. +// 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&)) + double (*f)(const triple&), int n) { double MX=-P[0].getx(); double MY=-P[0].gety(); double Z=P[0].getz(); - double MZ=-P[0].getz(); - for(int i=1; i < 16; ++i) { + double MZ=-Z; + for(int i=1; i < n; ++i) { triple v=P[i]; MX=m(MX,-v.getx()); MY=m(MY,-v.gety()); @@ -892,7 +892,7 @@ 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)) >= -fuzz || depth == 0) + if(m(-1.0,1.0)*(b-ratiobound(P,m,f,16)) >= -fuzz || depth == 0) return b; --depth; @@ -924,6 +924,168 @@ double bound(triple *P, double (*m)(double, double), 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; + + 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; + + 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) @@ -956,24 +1118,24 @@ void bounds(triple& Pmin, triple& Pmax, triple *P, double fuzz) 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 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); - double ymax=bound(Py,max,by,fuzz); + 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); - double zmax=bound(Pz,max,bz,fuzz); + 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); } |