Speed up min, max, minratio, maxratio, and intersection calculations …

…by increasing fuzz as depth decreases.

drawsurface.cc

@@ -143,17 +143,17 @@ void drawBezierPatch::bounds(const double* t, bbox3& b)
     }
 
     double c0=cx[0];
-    double fuzz=sqrtFuzz*run::norm(cx,16);
+    double fuzz=Fuzz*run::norm(cx,16);
     x=bound(cx,min,b.empty ? c0 : min(c0,b.left),fuzz,maxdepth);
     X=bound(cx,max,b.empty ? c0 : max(c0,b.right),fuzz,maxdepth);
 
     c0=cy[0];
-    fuzz=sqrtFuzz*run::norm(cy,16);
+    fuzz=Fuzz*run::norm(cy,16);
     y=bound(cy,min,b.empty ? c0 : min(c0,b.bottom),fuzz,maxdepth);
     Y=boundtri(cy,max,b.empty ? c0 : max(c0,b.top),fuzz,maxdepth);
 
     c0=cz[0];
-    fuzz=sqrtFuzz*run::norm(cz,16);
+    fuzz=Fuzz*run::norm(cz,16);
     z=bound(cz,min,b.empty ? c0 : min(c0,b.lower),fuzz,maxdepth);
     Z=bound(cz,max,b.empty ? c0 : max(c0,b.upper),fuzz,maxdepth);
   }  
@@ -357,17 +357,17 @@ void drawBezierTriangle::bounds(const double* t, bbox3& b)
     }
 
     double c0=cx[0];
-    double fuzz=sqrtFuzz*run::norm(cx,10);
+    double fuzz=Fuzz*run::norm(cx,10);
     x=boundtri(cx,min,b.empty ? c0 : min(c0,b.left),fuzz,maxdepth);
     X=boundtri(cx,max,b.empty ? c0 : max(c0,b.right),fuzz,maxdepth);
 
     c0=cy[0];
-    fuzz=sqrtFuzz*run::norm(cy,10);
+    fuzz=Fuzz*run::norm(cy,10);
     y=boundtri(cy,min,b.empty ? c0 : min(c0,b.bottom),fuzz,maxdepth);
     Y=boundtri(cy,max,b.empty ? c0 : max(c0,b.top),fuzz,maxdepth);
 
     c0=cz[0];
-    fuzz=sqrtFuzz*run::norm(cz,10);
+    fuzz=Fuzz*run::norm(cz,10);
     z=boundtri(cz,min,b.empty ? c0 : min(c0,b.lower),fuzz,maxdepth);
     Z=boundtri(cz,max,b.empty ? c0 : max(c0,b.upper),fuzz,maxdepth);
   }

path.cc

@@ -18,10 +18,11 @@
 
 namespace camp {
 
-const double Fuzz=1000.0*DBL_EPSILON;
+const double AFuzz=1000.0*DBL_EPSILON;
+const double Fuzz=sqrt(AFuzz);
+
 const double BigFuzz=10.0*Fuzz;
 const double Fuzz2=Fuzz*Fuzz;
-const double sqrtFuzz=sqrt(Fuzz);
 const double fuzzFactor=10.0;
 
 const double third=1.0/3.0;
@@ -999,6 +1000,8 @@ bool intersections(double &s, double &t, std::vector<double>& S,
     // Overlapping bounding boxes
 
     --depth;
+    fuzz *= 2;
+
     if((maxp-minp).length()+(maxq-minq).length() <= fuzz || depth == 0) {
       if(single) {
         s=0.5;
@@ -1240,6 +1243,7 @@ bool checkcurve(const pair& z0, const pair& c0, const pair& c1,
 {
   if(depth == 0) return true;
   --depth;
+
   if(insidebbox(z0,c0,c1,z1,z)) {
     const pair m0=0.5*(z0+c0);
     const pair m1=0.5*(c0+c1);

path3.cc

@@ -336,7 +336,9 @@ double bound(triple z0, triple c0, triple c1, triple z1,
   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);
@@ -352,7 +354,7 @@ double bound(triple z0, triple c0, triple c1, triple z1,
 
 pair path3::ratio(double (*m)(double, double)) const
 {
-  double fuzz=sqrtFuzz*(max()-min()).length();
+  double fuzz=Fuzz*(max()-min()).length();
   checkEmpty3(n);
 
   triple v=point((Int) 0);
@@ -611,6 +613,8 @@ bool intersections(double &s, double &t, std::vector<double>& S,
     // Overlapping bounding boxes
 
     --depth;
+    fuzz *= 2;
+
     if((maxp-minp).length()+(maxq-minq).length() <= fuzz || depth == 0) {
       if(single) {
         s=0.5;
@@ -815,7 +819,9 @@ double bound(double *P, double (*m)(double, double), double b,
   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]);
@@ -894,8 +900,10 @@ double bound(triple *P, double (*m)(double, double),
   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;
 
+  --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]);
@@ -1038,7 +1046,9 @@ double boundtri(double *P, double (*m)(double, double), double b,
   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);
 
@@ -1065,7 +1075,9 @@ double boundtri(triple *P, double (*m)(double, double),
   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);
 
@@ -1178,6 +1190,8 @@ bool intersections(double& U, double& V, const triple& v, triple *P,
      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;
@@ -1271,6 +1285,7 @@ bool intersections(std::vector<double>& T, std::vector<double>& U,
      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);

picture.cc

@@ -323,7 +323,7 @@ pair picture::ratio(double (*m)(double, double))
   bool first=true;
   pair b;
   bounds3();
-  double fuzz=sqrtFuzz*(b3.Max()-b3.Min()).length();
+  double fuzz=Fuzz*(b3.Max()-b3.Min()).length();
   matrixstack ms;
   for(nodelist::const_iterator p=nodes.begin(); p != nodes.end(); ++p) {
     assert(*p);

runarray.in

@@ -1915,9 +1915,9 @@ triple minbezier(triplearray2 *P, triple b)
   size_t N;
   real *A=copyTripleArray2Components(P,N);
   bound_double *B=bounddouble(N);
-  b=triple(B(A,::min,b.getx(),sqrtFuzz*norm(A,N),maxdepth),
-           B(A+N,::min,b.gety(),sqrtFuzz*norm(A+N,N),maxdepth),
-           B(A+2*N,::min,b.getz(),sqrtFuzz*norm(A+2*N,N),maxdepth));
+  b=triple(B(A,::min,b.getx(),Fuzz*norm(A,N),maxdepth),
+           B(A+N,::min,b.gety(),Fuzz*norm(A+N,N),maxdepth),
+           B(A+2*N,::min,b.getz(),Fuzz*norm(A+2*N,N),maxdepth));
   delete[] A;
   return b;
 }
@@ -1927,9 +1927,9 @@ triple maxbezier(triplearray2 *P, triple b)
   size_t N;
   real *A=copyTripleArray2Components(P,N);
   bound_double *B=bounddouble(N);
-  b=triple(B(A,::max,b.getx(),sqrtFuzz*norm(A,N),maxdepth),
-           B(A+N,::max,b.gety(),sqrtFuzz*norm(A+N,N),maxdepth),
-           B(A+2*N,::max,b.getz(),sqrtFuzz*norm(A+2*N,N),maxdepth));
+  b=triple(B(A,::max,b.getx(),Fuzz*norm(A,N),maxdepth),
+           B(A+N,::max,b.gety(),Fuzz*norm(A+N,N),maxdepth),
+           B(A+2*N,::max,b.getz(),Fuzz*norm(A+2*N,N),maxdepth));
   delete[] A;
   return b;
 }
@@ -1938,7 +1938,7 @@ pair minratio(triplearray2 *P, pair b)
 {
   size_t N;
   triple *A=copyTripleArray2C(P,N);
-  real fuzz=sqrtFuzz*norm(A,N);
+  real fuzz=Fuzz*norm(A,N);
   bound_triple *B=boundtriple(N);
   b=pair(B(A,::min,xratio,b.getx(),fuzz,maxdepth),
          B(A,::min,yratio,b.gety(),fuzz,maxdepth));
@@ -1951,7 +1951,7 @@ pair maxratio(triplearray2 *P, pair b)
   size_t N;
   triple *A=copyTripleArray2C(P,N);
   bound_triple *B=boundtriple(N);
-  real fuzz=sqrtFuzz*norm(A,N);
+  real fuzz=Fuzz*norm(A,N);
   b=pair(B(A,::max,xratio,b.getx(),fuzz,maxdepth),
          B(A,::max,yratio,b.gety(),fuzz,maxdepth));
   delete[] A;