"Fossies" - the Fresh Open Source Software Archive

Member "asymptote-2.60/path.cc" (6 Nov 2019, 33988 Bytes) of package /linux/misc/asymptote-2.60.src.tgz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "path.cc" see the Fossies "Dox" file reference documentation.

    1 /*****
    2  * path.cc
    3  * Andy Hammerlindl 2002/06/06
    4  *
    5  * Stores and returns information on a predefined path.
    6  *
    7  * When changing the path algorithms, also update the corresponding 
    8  * three-dimensional algorithms in path3.cc.
    9  *****/
   10 
   11 #include "path.h"
   12 #include "util.h"
   13 #include "angle.h"
   14 #include "camperror.h"
   15 #include "mathop.h"
   16 #include "predicates.h"
   17 #include "rounding.h"
   18 
   19 namespace camp {
   20 
   21 const double Fuzz2=1000.0*DBL_EPSILON;
   22 const double Fuzz=sqrt(Fuzz2);
   23 const double Fuzz4=Fuzz2*Fuzz2;
   24 const double BigFuzz=10.0*Fuzz2;
   25 const double fuzzFactor=100.0;
   26 
   27 const double third=1.0/3.0;
   28 
   29 path nullpath;
   30   
   31 const char *nopoints="nullpath has no points";
   32 
   33 void checkEmpty(Int n) {
   34   if(n == 0)
   35     reportError(nopoints);
   36 }
   37 
   38 // Accurate computation of sqrt(1+x)-1.
   39 inline double sqrt1pxm1(double x)
   40 {
   41   return x/(sqrt(1.0+x)+1.0);
   42 }
   43 inline pair sqrt1pxm1(pair x)
   44 {
   45   return x/(Sqrt(1.0+x)+1.0);
   46 }
   47   
   48 // Solve for the real roots of the quadratic equation ax^2+bx+c=0.
   49 quadraticroots::quadraticroots(double a, double b, double c)
   50 {
   51   // Remove roots at numerical infinity.
   52   if(fabs(a) <= Fuzz2*fabs(b)+Fuzz4*fabs(c)) {
   53     if(fabs(b) > Fuzz2*fabs(c)) {
   54       distinct=quadraticroots::ONE;
   55       roots=1;
   56       t1=-c/b;
   57     } else if(c == 0.0) {
   58       distinct=quadraticroots::MANY;
   59       roots=1;
   60       t1=0.0;
   61     } else {
   62       distinct=quadraticroots::NONE;
   63       roots=0;
   64     }
   65   } else {
   66     double factor=0.5*b/a;
   67     double denom=b*factor;
   68     if(fabs(denom) <= Fuzz2*fabs(c)) {
   69       double x=-c/a;
   70       if(x >= 0.0) {
   71         distinct=quadraticroots::TWO;
   72         roots=2;
   73         t2=sqrt(x);
   74         t1=-t2;
   75       } else {
   76         distinct=quadraticroots::NONE;
   77         roots=0;
   78       }
   79     } else {
   80       double x=-2.0*c/denom;
   81       if(x > -1.0) {
   82         distinct=quadraticroots::TWO;
   83         roots=2;
   84         double r2=factor*sqrt1pxm1(x);
   85         double r1=-r2-2.0*factor;
   86         if(r1 <= r2) {
   87           t1=r1;
   88           t2=r2;
   89         } else {
   90           t1=r2;
   91           t2=r1;
   92         }
   93       } else if(x == -1.0) {
   94         distinct=quadraticroots::ONE;
   95         roots=2;
   96         t1=t2=-factor;
   97       } else {
   98         distinct=quadraticroots::NONE;
   99         roots=0;
  100       }
  101     }
  102   }
  103 }
  104 
  105 // Solve for the complex roots of the quadratic equation ax^2+bx+c=0.
  106 Quadraticroots::Quadraticroots(pair a, pair b, pair c)
  107 {
  108   if(a == 0.0) {
  109     if(b != 0.0) {
  110       roots=1;
  111       z1=-c/b;
  112     } else if(c == 0.0) {
  113       roots=1;
  114       z1=0.0;
  115     } else
  116       roots=0;
  117   } else {
  118     roots=2;
  119     pair factor=0.5*b/a;
  120     pair denom=b*factor;
  121     if(denom == 0.0) {
  122       z1=Sqrt(-c/a);
  123       z2=-z1;
  124     } else {
  125       z1=factor*sqrt1pxm1(-2.0*c/denom);
  126       z2=-z1-2.0*factor;
  127     }
  128   }
  129 }
  130 
  131 inline bool goodroot(double a, double b, double c, double t)
  132 {
  133   return goodroot(t) && quadratic(a,b,c,t) >= 0.0;
  134 }
  135 
  136 // Accurate computation of cbrt(sqrt(1+x)+1)-cbrt(sqrt(1+x)-1).
  137 inline double cbrtsqrt1pxm(double x)
  138 {
  139   double s=sqrt1pxm1(x);
  140   return 2.0/(cbrt(x+2.0*(sqrt(1.0+x)+1.0))+cbrt(x)+cbrt(s*s));
  141 }
  142   
  143 // Taylor series of cos((atan(1.0/w)+pi)/3.0).
  144 static inline double costhetapi3(double w)
  145 {
  146   static const double c1=1.0/3.0;
  147   static const double c3=-19.0/162.0;
  148   static const double c5=425.0/5832.0;
  149   static const double c7=-16829.0/314928.0;
  150   double w2=w*w;
  151   double w3=w2*w;
  152   double w5=w3*w2;
  153   return c1*w+c3*w3+c5*w5+c7*w5*w2;
  154 }
  155       
  156 // Solve for the real roots of the cubic equation ax^3+bx^2+cx+d=0.
  157 cubicroots::cubicroots(double a, double b, double c, double d) 
  158 {
  159   static const double ninth=1.0/9.0;
  160   static const double fiftyfourth=1.0/54.0;
  161   
  162   // Remove roots at numerical infinity.
  163   if(fabs(a) <= Fuzz2*(fabs(b)+fabs(c)*Fuzz2+fabs(d)*Fuzz4)) {
  164     quadraticroots q(b,c,d);
  165     roots=q.roots;
  166     if(q.roots >= 1) t1=q.t1;
  167     if(q.roots == 2) t2=q.t2;
  168     return;
  169   }
  170   
  171   // Detect roots at numerical zero.
  172   if(fabs(d) <= Fuzz2*(fabs(c)+fabs(b)*Fuzz2+fabs(a)*Fuzz4)) {
  173     quadraticroots q(a,b,c);
  174     roots=q.roots+1;
  175     t1=0;
  176     if(q.roots >= 1) t2=q.t1;
  177     if(q.roots == 2) t3=q.t2;
  178     return;
  179   }
  180   
  181   b /= a;
  182   c /= a;
  183   d /= a;
  184   
  185   double b2=b*b;
  186   double Q=3.0*c-b2;
  187   if(fabs(Q) < Fuzz2*(3.0*fabs(c)+fabs(b2)))
  188     Q=0.0;
  189   
  190   double R=(3.0*Q+b2)*b-27.0*d;
  191   if(fabs(R) < Fuzz2*((3.0*fabs(Q)+fabs(b2))*fabs(b)+27.0*fabs(d)))
  192     R=0.0;
  193   
  194   Q *= ninth;
  195   R *= fiftyfourth;
  196   
  197   double Q3=Q*Q*Q;
  198   double R2=R*R;
  199   double D=Q3+R2;
  200   double mthirdb=-b*third;
  201   
  202   if(D > 0.0) {
  203     roots=1;
  204     t1=mthirdb;
  205     if(R2 != 0.0) t1 += cbrt(R)*cbrtsqrt1pxm(Q3/R2);
  206   } else {
  207     roots=3;
  208     double v=0.0,theta;
  209     if(R2 > 0.0) {
  210       v=sqrt(-D/R2);
  211       theta=atan(v);
  212     } else theta=0.5*PI;
  213     double factor=2.0*sqrt(-Q)*(R >= 0 ? 1 : -1);
  214       
  215     t1=mthirdb+factor*cos(third*theta);
  216     t2=mthirdb-factor*cos(third*(theta-PI));
  217     t3=mthirdb;
  218     if(R2 > 0.0)
  219       t3 -= factor*((v < 100.0) ? cos(third*(theta+PI)) : costhetapi3(1.0/v)); 
  220   }
  221 }
  222   
  223 pair path::point(double t) const
  224 {
  225   checkEmpty(n);
  226     
  227   Int i = Floor(t);
  228   Int iplus;
  229   t = fmod(t,1);
  230   if (t < 0) t += 1;
  231 
  232   if (cycles) {
  233     i = imod(i,n);
  234     iplus = imod(i+1,n);
  235   }
  236   else if (i < 0)
  237     return nodes[0].point;
  238   else if (i >= n-1)
  239     return nodes[n-1].point;
  240   else
  241     iplus = i+1;
  242 
  243   double one_t = 1.0-t;
  244 
  245   pair a = nodes[i].point,
  246     b = nodes[i].post,
  247     c = nodes[iplus].pre,
  248     d = nodes[iplus].point,
  249     ab   = one_t*a   + t*b,
  250     bc   = one_t*b   + t*c,
  251     cd   = one_t*c   + t*d,
  252     abc  = one_t*ab  + t*bc,
  253     bcd  = one_t*bc  + t*cd,
  254     abcd = one_t*abc + t*bcd;
  255 
  256   return abcd;
  257 }
  258 
  259 pair path::precontrol(double t) const
  260 {
  261   checkEmpty(n);
  262                      
  263   Int i = Floor(t);
  264   Int iplus;
  265   t = fmod(t,1);
  266   if (t < 0) t += 1;
  267 
  268   if (cycles) {
  269     i = imod(i,n);
  270     iplus = imod(i+1,n);
  271   }
  272   else if (i < 0)
  273     return nodes[0].pre;
  274   else if (i >= n-1)
  275     return nodes[n-1].pre;
  276   else
  277     iplus = i+1;
  278 
  279   double one_t = 1.0-t;
  280 
  281   pair a = nodes[i].point,
  282     b = nodes[i].post,
  283     c = nodes[iplus].pre,
  284     ab   = one_t*a   + t*b,
  285     bc   = one_t*b   + t*c,
  286     abc  = one_t*ab  + t*bc;
  287 
  288   return (abc == a) ? nodes[i].pre : abc;
  289 }
  290         
  291  
  292 pair path::postcontrol(double t) const
  293 {
  294   checkEmpty(n);
  295   
  296   Int i = Floor(t);
  297   Int iplus;
  298   t = fmod(t,1);
  299   if (t < 0) t += 1;
  300 
  301   if (cycles) {
  302     i = imod(i,n);
  303     iplus = imod(i+1,n);
  304   }
  305   else if (i < 0)
  306     return nodes[0].post;
  307   else if (i >= n-1)
  308     return nodes[n-1].post;
  309   else
  310     iplus = i+1;
  311 
  312   double one_t = 1.0-t;
  313   
  314   pair b = nodes[i].post,
  315     c = nodes[iplus].pre,
  316     d = nodes[iplus].point,
  317     bc   = one_t*b   + t*c,
  318     cd   = one_t*c   + t*d,
  319     bcd  = one_t*bc  + t*cd;
  320 
  321   return (bcd == d) ? nodes[iplus].post : bcd;
  322 }
  323 
  324 path path::reverse() const
  325 {
  326   mem::vector<solvedKnot> nodes(n);
  327   Int len=length();
  328   for (Int i = 0, j = len; i < n; i++, j--) {
  329     nodes[i].pre = postcontrol(j);
  330     nodes[i].point = point(j);
  331     nodes[i].post = precontrol(j);
  332     nodes[i].straight = straight(j-1);
  333   }
  334   return path(nodes, n, cycles);
  335 }
  336 
  337 path path::subpath(Int a, Int b) const
  338 {
  339   if(empty()) return path();
  340 
  341   if (a > b) {
  342     const path &rp = reverse();
  343     Int len=length();
  344     path result = rp.subpath(len-a, len-b);
  345     return result;
  346   }
  347 
  348   if (!cycles) {
  349     if (a < 0) {
  350       a = 0;
  351       if(b < 0)
  352         b = 0;
  353     }
  354     if (b > n-1) {
  355       b = n-1;
  356       if(a > b)
  357         a = b;
  358     }
  359   }
  360 
  361   Int sn = b-a+1;
  362   mem::vector<solvedKnot> nodes(sn);
  363 
  364   for (Int i = 0, j = a; j <= b; i++, j++) {
  365     nodes[i].pre = precontrol(j);
  366     nodes[i].point = point(j);
  367     nodes[i].post = postcontrol(j);
  368     nodes[i].straight = straight(j);
  369   }
  370   nodes[0].pre = nodes[0].point;
  371   nodes[sn-1].post = nodes[sn-1].point;
  372 
  373   return path(nodes, sn);
  374 }
  375 
  376 inline pair split(double t, const pair& x, const pair& y) { return x+(y-x)*t; }
  377 
  378 inline void splitCubic(solvedKnot sn[], double t, const solvedKnot& left_,
  379                        const solvedKnot& right_)
  380 {
  381   solvedKnot &left=(sn[0]=left_), &mid=sn[1], &right=(sn[2]=right_);
  382   if(left.straight) {
  383     mid.point=split(t,left.point,right.point);
  384     pair deltaL=third*(mid.point-left.point);
  385     left.post=left.point+deltaL;
  386     mid.pre=mid.point-deltaL;
  387     pair deltaR=third*(right.point-mid.point);
  388     mid.post=mid.point+deltaR;
  389     right.pre=right.point-deltaR;
  390     mid.straight=true;
  391   } else {
  392     pair x=split(t,left.post,right.pre); // m1
  393     left.post=split(t,left.point,left.post); // m0
  394     right.pre=split(t,right.pre,right.point); // m2
  395     mid.pre=split(t,left.post,x); // m3
  396     mid.post=split(t,x,right.pre); // m4 
  397     mid.point=split(t,mid.pre,mid.post); // m5
  398   }
  399 }
  400 
  401 path path::subpath(double a, double b) const
  402 {
  403   if(empty()) return path();
  404   
  405   if (a > b) {
  406     const path &rp = reverse();
  407     Int len=length();
  408     return rp.subpath(len-a, len-b);
  409   }
  410 
  411   solvedKnot aL, aR, bL, bR;
  412   if (!cycles) {
  413     if (a < 0) {
  414       a = 0;
  415       if (b < 0)
  416         b = 0;
  417     }   
  418     if (b > n-1) {
  419       b = n-1;
  420       if (a > b)
  421         a = b;
  422     }
  423     aL = nodes[(Int)floor(a)];
  424     aR = nodes[(Int)ceil(a)];
  425     bL = nodes[(Int)floor(b)];
  426     bR = nodes[(Int)ceil(b)];
  427   } else {
  428     if(run::validInt(a) && run::validInt(b)) {
  429       aL = nodes[imod((Int) floor(a),n)];
  430       aR = nodes[imod((Int) ceil(a),n)];
  431       bL = nodes[imod((Int) floor(b),n)];
  432       bR = nodes[imod((Int) ceil(b),n)];
  433     } else reportError("invalid path index");
  434   }
  435 
  436   if (a == b) return path(point(a));
  437 
  438   solvedKnot sn[3];
  439   path p = subpath(Ceil(a), Floor(b));
  440   if (a > floor(a)) {
  441     if (b < ceil(a)) {
  442       splitCubic(sn,a-floor(a),aL,aR);
  443       splitCubic(sn,(b-a)/(ceil(b)-a),sn[1],sn[2]);
  444       return path(sn[0],sn[1]);
  445     }
  446     splitCubic(sn,a-floor(a),aL,aR);
  447     p=concat(path(sn[1],sn[2]),p);
  448   }
  449   if (ceil(b) > b) {
  450     splitCubic(sn,b-floor(b),bL,bR);
  451     p=concat(p,path(sn[0],sn[1]));
  452   }
  453   return p;
  454 }
  455 
  456 // Special case of subpath for paths of length 1 used by intersect.
  457 void path::halve(path &first, path &second) const
  458 {
  459   solvedKnot sn[3];
  460   splitCubic(sn,0.5,nodes[0],nodes[1]);
  461   first=path(sn[0],sn[1]);
  462   second=path(sn[1],sn[2]);
  463 }
  464   
  465 // Calculate the coefficients of a Bezier derivative divided by 3.
  466 static inline void derivative(pair& a, pair& b, pair& c,
  467                               const pair& z0, const pair& c0,
  468                               const pair& c1, const pair& z1)
  469 {
  470   a=z1-z0+3.0*(c0-c1);
  471   b=2.0*(z0+c1)-4.0*c0;
  472   c=c0-z0;
  473 }
  474 
  475 bbox path::bounds() const
  476 {
  477   if(!box.empty) return box;
  478   
  479   if (empty()) {
  480     // No bounds
  481     return bbox();
  482   }
  483   
  484   Int len=length();
  485   box.add(point(len));
  486   times=bbox(len,len,len,len);
  487 
  488   for (Int i = 0; i < len; i++) {
  489     addpoint(box,i);
  490     if(straight(i)) continue;
  491     
  492     pair a,b,c;
  493     derivative(a,b,c,point(i),postcontrol(i),precontrol(i+1),point(i+1));
  494     
  495     // Check x coordinate
  496     quadraticroots x(a.getx(),b.getx(),c.getx());
  497     if(x.distinct != quadraticroots::NONE && goodroot(x.t1))
  498       addpoint(box,i+x.t1);
  499     if(x.distinct == quadraticroots::TWO && goodroot(x.t2))
  500       addpoint(box,i+x.t2);
  501     
  502     // Check y coordinate
  503     quadraticroots y(a.gety(),b.gety(),c.gety());
  504     if(y.distinct != quadraticroots::NONE && goodroot(y.t1))
  505       addpoint(box,i+y.t1);
  506     if(y.distinct == quadraticroots::TWO && goodroot(y.t2))
  507       addpoint(box,i+y.t2);
  508   }
  509   return box;
  510 }
  511 
  512 bbox path::bounds(double min, double max) const
  513 {
  514   bbox box;
  515   
  516   Int len=length();
  517   for (Int i = 0; i < len; i++) {
  518     addpoint(box,i,min,max);
  519     if(straight(i)) continue;
  520     
  521     pair a,b,c;
  522     derivative(a,b,c,point(i),postcontrol(i),precontrol(i+1),point(i+1));
  523     
  524     // Check x coordinate
  525     quadraticroots x(a.getx(),b.getx(),c.getx());
  526     if(x.distinct != quadraticroots::NONE && goodroot(x.t1))
  527       addpoint(box,i+x.t1,min,max);
  528 
  529     if(x.distinct == quadraticroots::TWO && goodroot(x.t2))
  530       addpoint(box,i+x.t2,min,max);
  531     
  532     // Check y coordinate
  533     quadraticroots y(a.gety(),b.gety(),c.gety());
  534     if(y.distinct != quadraticroots::NONE && goodroot(y.t1))
  535       addpoint(box,i+y.t1,min,max);
  536     if(y.distinct == quadraticroots::TWO && goodroot(y.t2))
  537       addpoint(box,i+y.t2,min,max);
  538   }
  539   addpoint(box,len,min,max);
  540   return box;
  541 }
  542   
  543 inline void add(bbox& box, const pair& z, const pair& min, const pair& max)
  544 {
  545   box += z+min;
  546   box += z+max;
  547 }
  548 
  549 bbox path::internalbounds(const bbox& padding) const
  550 {
  551   bbox box;
  552   
  553   // Check interior nodes.
  554   Int len=length();
  555   for (Int i = 1; i < len; i++) {
  556     pair pre=point(i)-precontrol(i);
  557     pair post=postcontrol(i)-point(i);
  558     
  559     // Check node x coordinate
  560     if((pre.getx() >= 0.0) ^ (post.getx() >= 0))
  561       add(box,point(i),padding.left,padding.right);
  562                               
  563     // Check node y coordinate
  564     if((pre.gety() >= 0.0) ^ (post.gety() >= 0))
  565       add(box,point(i),pair(0,padding.bottom),pair(0,padding.top));
  566   }
  567                               
  568   // Check interior segments.
  569   for (Int i = 0; i < len; i++) {
  570     if(straight(i)) continue;
  571     
  572     pair a,b,c;
  573     derivative(a,b,c,point(i),postcontrol(i),precontrol(i+1),point(i+1));
  574     
  575     // Check x coordinate
  576     quadraticroots x(a.getx(),b.getx(),c.getx());
  577     if(x.distinct != quadraticroots::NONE && goodroot(x.t1))
  578       add(box,point(i+x.t1),padding.left,padding.right);
  579     if(x.distinct == quadraticroots::TWO && goodroot(x.t2))
  580       add(box,point(i+x.t2),padding.left,padding.right);
  581     
  582     // Check y coordinate
  583     quadraticroots y(a.gety(),b.gety(),c.gety());
  584     if(y.distinct != quadraticroots::NONE && goodroot(y.t1))
  585       add(box,point(i+y.t1),pair(0,padding.bottom),pair(0,padding.top));
  586     if(y.distinct == quadraticroots::TWO && goodroot(y.t2))
  587       add(box,point(i+y.t2),pair(0,padding.bottom),pair(0,padding.top));
  588   }
  589   return box;
  590 }
  591 
  592 // {{{ Arclength Calculations
  593 
  594 static pair a,b,c;
  595 
  596 static double ds(double t)
  597 {
  598   double dx=quadratic(a.getx(),b.getx(),c.getx(),t);
  599   double dy=quadratic(a.gety(),b.gety(),c.gety(),t);
  600   return sqrt(dx*dx+dy*dy);
  601 }
  602 
  603 // Calculates arclength of a cubic using adaptive simpson integration.
  604 double path::cubiclength(Int i, double goal) const
  605 {
  606   const pair& z0=point(i);
  607   const pair& z1=point(i+1);
  608   double L;
  609   if(straight(i)) {
  610     L=(z1-z0).length();
  611     return (goal < 0 || goal >= L) ? L : -goal/L;
  612   }
  613   
  614   const pair& c0=postcontrol(i);
  615   const pair& c1=precontrol(i+1);
  616   
  617   double integral;
  618   derivative(a,b,c,z0,c0,c1,z1);
  619   
  620   if(!simpson(integral,ds,0.0,1.0,DBL_EPSILON,1.0))
  621     reportError("nesting capacity exceeded in computing arclength");
  622   L=3.0*integral;
  623   if(goal < 0 || goal >= L) return L;
  624   
  625   double t=goal/L;
  626   goal *= third;
  627   static double dxmin=sqrt(DBL_EPSILON);
  628   if(!unsimpson(goal,ds,0.0,t,100.0*DBL_EPSILON,integral,1.0,dxmin))
  629     reportError("nesting capacity exceeded in computing arctime");
  630   return -t;
  631 }
  632 
  633 double path::arclength() const 
  634 {
  635   if (cached_length != -1) return cached_length;
  636 
  637   double L=0.0;
  638   for (Int i = 0; i < n-1; i++) {
  639     L += cubiclength(i);
  640   }
  641   if(cycles) L += cubiclength(n-1);
  642   cached_length = L;
  643   return cached_length;
  644 }
  645 
  646 double path::arctime(double goal) const
  647 {
  648   if (cycles) {
  649     if (goal == 0 || cached_length == 0) return 0;
  650     if (goal < 0)  {
  651       const path &rp = this->reverse();
  652       double result = -rp.arctime(-goal);
  653       return result;
  654     }
  655     if (cached_length > 0 && goal >= cached_length) {
  656       Int loops = (Int)(goal / cached_length);
  657       goal -= loops*cached_length;
  658       return loops*n+arctime(goal);
  659     }      
  660   } else {
  661     if (goal <= 0)
  662       return 0;
  663     if (cached_length > 0 && goal >= cached_length)
  664       return n-1;
  665   }
  666     
  667   double l,L=0;
  668   for (Int i = 0; i < n-1; i++) {
  669     l = cubiclength(i,goal);
  670     if (l < 0)
  671       return (-l+i);
  672     else {
  673       L += l;
  674       goal -= l;
  675       if (goal <= 0)
  676         return i+1;
  677     }
  678   }
  679   if (cycles) {
  680     l = cubiclength(n-1,goal);
  681     if (l < 0)
  682       return -l+n-1;
  683     if (cached_length > 0 && cached_length != L+l) {
  684       reportError("arclength != length.\n"
  685                   "path::arclength(double) must have broken semantics.\n"
  686                   "Please report this error.");
  687     }
  688     cached_length = L += l;
  689     goal -= l;
  690     return arctime(goal)+n;
  691   }
  692   else {
  693     cached_length = L;
  694     return length();
  695   }
  696 }
  697 
  698 // }}}
  699 
  700 // {{{ Direction Time Calulation
  701 // Algorithm Stolen from Knuth's MetaFont
  702 inline double cubicDir(const solvedKnot& left, const solvedKnot& right,
  703                        const pair& rot)
  704 {
  705   pair a,b,c;
  706   derivative(a,b,c,left.point,left.post,right.pre,right.point);
  707   a *= rot; b *= rot; c *= rot;
  708   
  709   quadraticroots ret(a.gety(),b.gety(),c.gety());
  710   switch(ret.distinct) {
  711     case quadraticroots::MANY:
  712     case quadraticroots::ONE:
  713     {
  714       if(goodroot(a.getx(),b.getx(),c.getx(),ret.t1)) return ret.t1;
  715     } break;
  716 
  717     case quadraticroots::TWO:
  718     {
  719       if(goodroot(a.getx(),b.getx(),c.getx(),ret.t1)) return ret.t1;
  720       if(goodroot(a.getx(),b.getx(),c.getx(),ret.t2)) return ret.t2;
  721     } break;
  722 
  723     case quadraticroots::NONE:
  724       break;
  725   }
  726 
  727   return -1;
  728 }
  729 
  730 // TODO: Check that we handle corner cases.
  731 // Velocity(t) == (0,0)
  732 double path::directiontime(const pair& dir) const {
  733   if (dir == pair(0,0)) return 0;
  734   pair rot = pair(1,0)/unit(dir);
  735     
  736   double t; double pre,post;
  737   for (Int i = 0; i < n-1+cycles; ) {
  738     t = cubicDir(this->nodes[i],(cycles && i==n-1) ? nodes[0]:nodes[i+1],rot);
  739     if (t >= 0) return i+t;
  740     i++;
  741     if (cycles || i != n-1) {
  742       pair Pre = (point(i)-precontrol(i))*rot;
  743       pair Post = (postcontrol(i)-point(i))*rot;
  744       static pair zero(0.0,0.0);
  745       if(Pre != zero && Post != zero) {
  746         pre = angle(Pre);
  747         post = angle(Post);
  748         if ((pre <= 0 && post >= 0 && pre >= post - PI) ||
  749             (pre >= 0 && post <= 0 && pre <= post + PI))
  750           return i;
  751       }
  752     }
  753   }
  754   
  755   return -1;
  756 }
  757 // }}}
  758 
  759 // {{{ Path Intersection Calculations
  760 
  761 const unsigned maxdepth=DBL_MANT_DIG;
  762 const unsigned mindepth=maxdepth-16;
  763 
  764 void roots(std::vector<double> &roots, double a, double b, double c, double d)
  765 {
  766   cubicroots r(a,b,c,d);
  767   if(r.roots >= 1) roots.push_back(r.t1);
  768   if(r.roots >= 2) roots.push_back(r.t2);
  769   if(r.roots == 3) roots.push_back(r.t3);
  770 }
  771   
  772 void roots(std::vector<double> &r, double x0, double c0, double c1, double x1,
  773            double x)
  774 {
  775   double a=x1-x0+3.0*(c0-c1);
  776   double b=3.0*(x0+c1)-6.0*c0;
  777   double c=3.0*(c0-x0);
  778   double d=x0-x;
  779   roots(r,a,b,c,d);
  780 }
  781 
  782 // Return all intersection times of path g with the pair z.
  783 void intersections(std::vector<double>& T, const path& g, const pair& z,
  784                    double fuzz)
  785 {
  786   double fuzz2=fuzz*fuzz;
  787   Int n=g.length();
  788   bool cycles=g.cyclic();
  789   for(Int i=0; i < n; ++i) {
  790     // Check both directions to circumvent degeneracy.
  791     std::vector<double> r;
  792     roots(r,g.point(i).getx(),g.postcontrol(i).getx(),
  793           g.precontrol(i+1).getx(),g.point(i+1).getx(),z.getx());
  794     roots(r,g.point(i).gety(),g.postcontrol(i).gety(),
  795           g.precontrol(i+1).gety(),g.point(i+1).gety(),z.gety());
  796     
  797     size_t m=r.size();
  798     for(size_t j=0 ; j < m; ++j) {
  799       double t=r[j];
  800       if(t >= -Fuzz2 && t <= 1.0+Fuzz2) {
  801         double s=i+t;
  802         if((g.point(s)-z).abs2() <= fuzz2) {
  803           if(cycles && s >= n-Fuzz2) s=0;
  804           T.push_back(s);
  805         }
  806       }
  807     }
  808   }
  809 }
  810 
  811 inline bool online(const pair&p, const pair& q, const pair& z, double fuzz)
  812 {
  813   if(p == q) return (z-p).abs2() <= fuzz*fuzz;
  814   return (z.getx()-p.getx())*(q.gety()-p.gety()) ==
  815     (q.getx()-p.getx())*(z.gety()-p.gety());
  816 }
  817 
  818 // Return all intersection times of path g with the (infinite)
  819 // line through p and q; if there are an infinite number of intersection points,
  820 // the returned list is guaranteed to include the endpoint times of
  821 // the intersection if endpoints=true.
  822 void lineintersections(std::vector<double>& T, const path& g,
  823                        const pair& p, const pair& q, double fuzz,
  824                        bool endpoints=false)
  825 {
  826   Int n=g.length();
  827   if(n == 0) {
  828     if(online(p,q,g.point((Int) 0),fuzz)) T.push_back(0.0);
  829     return;
  830   }
  831   bool cycles=g.cyclic();
  832   double dx=q.getx()-p.getx();
  833   double dy=q.gety()-p.gety();
  834   double det=p.gety()*q.getx()-p.getx()*q.gety();
  835   for(Int i=0; i < n; ++i) {
  836     pair z0=g.point(i);
  837     pair c0=g.postcontrol(i);
  838     pair c1=g.precontrol(i+1);
  839     pair z1=g.point(i+1);
  840     pair t3=z1-z0+3.0*(c0-c1);
  841     pair t2=3.0*(z0+c1)-6.0*c0;
  842     pair t1=3.0*(c0-z0);
  843     double a=dy*t3.getx()-dx*t3.gety();
  844     double b=dy*t2.getx()-dx*t2.gety();
  845     double c=dy*t1.getx()-dx*t1.gety();
  846     double d=dy*z0.getx()-dx*z0.gety()+det;
  847     std::vector<double> r;
  848     if(max(max(max(a*a,b*b),c*c),d*d) >
  849        Fuzz4*max(max(max(z0.abs2(),z1.abs2()),c0.abs2()),c1.abs2()))
  850       roots(r,a,b,c,d);
  851     else r.push_back(0.0);
  852     if(endpoints) {
  853       path h=g.subpath(i,i+1);
  854       intersections(r,h,p,fuzz);
  855       intersections(r,h,q,fuzz);
  856       if(online(p,q,z0,fuzz)) r.push_back(0.0);
  857       if(online(p,q,z1,fuzz)) r.push_back(1.0);
  858     }
  859     size_t m=r.size();
  860     for(size_t j=0 ; j < m; ++j) {
  861       double t=r[j];
  862       if(t >= -Fuzz2 && t <= 1.0+Fuzz2) {
  863         double s=i+t;
  864         if(cycles && s >= n-Fuzz2) s=0;
  865         T.push_back(s);
  866       }
  867     }
  868   }
  869 }
  870 
  871 // An optimized implementation of intersections(g,p--q);
  872 // if there are an infinite number of intersection points, the returned list is
  873 // only guaranteed to include the endpoint times of the intersection.
  874 void intersections(std::vector<double>& S, std::vector<double>& T,
  875                    const path& g, const pair& p, const pair& q, double fuzz)
  876 {
  877   double length2=(q-p).abs2();
  878   if(length2 == 0.0) {
  879     std::vector<double> S1;
  880     intersections(S1,g,p,fuzz);
  881     size_t n=S1.size();
  882     for(size_t i=0; i < n; ++i) {
  883       S.push_back(S1[i]);
  884       T.push_back(0.0);
  885     }
  886   } else {
  887     pair factor=(q-p)/length2;
  888     std::vector<double> S1;
  889     lineintersections(S1,g,p,q,fuzz,true);
  890     size_t n=S1.size();
  891     for(size_t i=0; i < n; ++i) {
  892       double s=S1[i];
  893       double t=dot(g.point(s)-p,factor);
  894       if(t >= -Fuzz2 && t <= 1.0+Fuzz2) {
  895         S.push_back(s);
  896         T.push_back(t);
  897       }
  898     }
  899   }
  900 }
  901 
  902 void add(std::vector<double>& S, double s, const path& p, double fuzz2)
  903 {
  904   pair z=p.point(s);
  905   size_t n=S.size();
  906   for(size_t i=0; i < n; ++i)
  907     if((p.point(S[i])-z).abs2() <= fuzz2) return;
  908   S.push_back(s);
  909 }
  910   
  911 void add(std::vector<double>& S, std::vector<double>& T, double s, double t,
  912          const path& p, double fuzz2)
  913 {
  914   pair z=p.point(s);
  915   size_t n=S.size();
  916   for(size_t i=0; i < n; ++i)
  917     if((p.point(S[i])-z).abs2() <= fuzz2) return;
  918   S.push_back(s);
  919   T.push_back(t);
  920 }
  921   
  922 void add(double& s, double& t, std::vector<double>& S, std::vector<double>& T,
  923          std::vector<double>& S1, std::vector<double>& T1,
  924          double pscale, double qscale, double poffset, double qoffset,
  925          const path& p, double fuzz2, bool single)
  926 {
  927   if(single) {
  928     s=s*pscale+poffset;
  929     t=t*qscale+qoffset;
  930   } else {
  931     size_t n=S1.size();
  932     for(size_t i=0; i < n; ++i)
  933       add(S,T,pscale*S1[i]+poffset,qscale*T1[i]+qoffset,p,fuzz2);
  934   }
  935 }
  936 
  937 void add(double& s, double& t, std::vector<double>& S, std::vector<double>& T,
  938          std::vector<double>& S1, std::vector<double>& T1,
  939          const path& p, double fuzz2, bool single)
  940 {
  941   size_t n=S1.size();
  942   if(single) {
  943     if(n > 0) {
  944       s=S1[0];
  945       t=T1[0];
  946     }
  947   } else {
  948     for(size_t i=0; i < n; ++i)
  949       add(S,T,S1[i],T1[i],p,fuzz2);
  950   }
  951 }
  952 
  953 void intersections(std::vector<double>& S, path& g,
  954                    const pair& p, const pair& q, double fuzz)
  955 {       
  956   double fuzz2=max(fuzzFactor*fuzz*fuzz,Fuzz2);
  957   std::vector<double> S1;
  958   lineintersections(S1,g,p,q,fuzz);
  959   size_t n=S1.size();
  960   for(size_t i=0; i < n; ++i)
  961     add(S,S1[i],g,fuzz2);
  962 }
  963 
  964 bool intersections(double &s, double &t, std::vector<double>& S,
  965                    std::vector<double>& T, path& p, path& q,
  966                    double fuzz, bool single, bool exact, unsigned depth)
  967 {
  968   if(errorstream::interrupt) throw interrupted();
  969   
  970   double fuzz2=max(fuzzFactor*fuzz*fuzz,Fuzz2);
  971   
  972   Int lp=p.length();
  973   if(((lp == 1 && p.straight(0)) || lp == 0) && exact) {
  974     std::vector<double> T1,S1;
  975     intersections(T1,S1,q,p.point((Int) 0),p.point(lp),fuzz);
  976     add(s,t,S,T,S1,T1,p,fuzz2,single);
  977     return S1.size() > 0;
  978   }
  979 
  980   Int lq=q.length();
  981   if(((lq == 1 && q.straight(0)) || lq == 0) && exact) {
  982     std::vector<double> S1,T1;
  983     intersections(S1,T1,p,q.point((Int) 0),q.point(lq),fuzz);
  984     add(s,t,S,T,S1,T1,p,fuzz2,single);
  985     return S1.size() > 0;
  986   }
  987   
  988   pair maxp=p.max();
  989   pair minp=p.min();
  990   pair maxq=q.max();
  991   pair minq=q.min();
  992   
  993   if(maxp.getx()+fuzz >= minq.getx() &&
  994      maxp.gety()+fuzz >= minq.gety() && 
  995      maxq.getx()+fuzz >= minp.getx() &&
  996      maxq.gety()+fuzz >= minp.gety()) {
  997     // Overlapping bounding boxes
  998 
  999     --depth;
 1000     fuzz *= 2.0;
 1001     fuzz2=max(fuzzFactor*fuzz*fuzz,Fuzz2);
 1002 
 1003     if((maxp-minp).length()+(maxq-minq).length() <= fuzz || depth == 0) {
 1004       if(single) {
 1005         s=0.5;
 1006         t=0.5;
 1007       } else {
 1008         S.push_back(0.5);
 1009         T.push_back(0.5);
 1010       }
 1011       return true;
 1012     }
 1013     
 1014     path p1,p2;
 1015     double pscale,poffset;
 1016     std::vector<double> S1,T1;
 1017     
 1018     if(lp <= 1) {
 1019       if(lp == 1) p.halve(p1,p2);
 1020       if(lp == 0 || p1 == p || p2 == p) {
 1021         intersections(T1,S1,q,p.point((Int) 0),p.point((Int) 0),fuzz);
 1022         add(s,t,S,T,S1,T1,p,fuzz2,single);
 1023         return S1.size() > 0;
 1024       }
 1025       pscale=poffset=0.5;
 1026     } else {
 1027       Int tp=lp/2;
 1028       p1=p.subpath(0,tp);
 1029       p2=p.subpath(tp,lp);
 1030       poffset=tp;
 1031       pscale=1.0;
 1032     }
 1033       
 1034     path q1,q2;
 1035     double qscale,qoffset;
 1036     
 1037     if(lq <= 1) {
 1038       if(lq == 1) q.halve(q1,q2);
 1039       if(lq == 0 || q1 == q || q2 == q) {
 1040         intersections(S1,T1,p,q.point((Int) 0),q.point((Int) 0),fuzz);
 1041         add(s,t,S,T,S1,T1,p,fuzz2,single);
 1042         return S1.size() > 0;
 1043       }
 1044       qscale=qoffset=0.5;
 1045     } else {
 1046       Int tq=lq/2;
 1047       q1=q.subpath(0,tq);
 1048       q2=q.subpath(tq,lq);
 1049       qoffset=tq;
 1050       qscale=1.0;
 1051     }
 1052       
 1053     bool Short=lp == 1 && lq == 1;
 1054     
 1055     static size_t maxcount=9;
 1056     size_t count=0;
 1057     
 1058     if(intersections(s,t,S1,T1,p1,q1,fuzz,single,exact,depth)) {
 1059       add(s,t,S,T,S1,T1,pscale,qscale,0.0,0.0,p,fuzz2,single);
 1060       if(single || depth <= mindepth)
 1061         return true;
 1062       count += S1.size();
 1063       if(Short && count > maxcount) return true;
 1064     }
 1065      
 1066     S1.clear();
 1067     T1.clear();
 1068     if(intersections(s,t,S1,T1,p1,q2,fuzz,single,exact,depth)) {
 1069       add(s,t,S,T,S1,T1,pscale,qscale,0.0,qoffset,p,fuzz2,single);
 1070       if(single || depth <= mindepth)
 1071         return true;
 1072       count += S1.size();
 1073       if(Short && count > maxcount) return true;
 1074     }
 1075     
 1076     S1.clear();
 1077     T1.clear();
 1078     if(intersections(s,t,S1,T1,p2,q1,fuzz,single,exact,depth)) {
 1079       add(s,t,S,T,S1,T1,pscale,qscale,poffset,0.0,p,fuzz2,single);
 1080       if(single || depth <= mindepth)
 1081         return true;
 1082       count += S1.size();
 1083       if(Short && count > maxcount) return true;
 1084     }
 1085     
 1086     S1.clear();
 1087     T1.clear();
 1088     if(intersections(s,t,S1,T1,p2,q2,fuzz,single,exact,depth)) {
 1089       add(s,t,S,T,S1,T1,pscale,qscale,poffset,qoffset,p,fuzz2,single);
 1090       if(single || depth <= mindepth)
 1091         return true;
 1092       count += S1.size();
 1093       if(Short && count > maxcount) return true;
 1094     }
 1095     
 1096     return S.size() > 0;
 1097   }
 1098   return false;
 1099 }
 1100 
 1101 // }}}
 1102 
 1103 ostream& operator<< (ostream& out, const path& p)
 1104 {
 1105   Int n = p.length();
 1106   if(n < 0)
 1107     out << "<nullpath>";
 1108   else {
 1109     for(Int i = 0; i < n; i++) {
 1110       out << p.point(i);
 1111       if(p.straight(i)) out << "--";
 1112       else
 1113         out << ".. controls " << p.postcontrol(i) << " and "
 1114             << p.precontrol(i+1) << newl << " ..";
 1115     }
 1116     if(p.cycles) 
 1117       out << "cycle";
 1118     else
 1119       out << p.point(n);
 1120   }
 1121   return out;
 1122 }
 1123 
 1124 path concat(const path& p1, const path& p2)
 1125 {
 1126   Int n1 = p1.length(), n2 = p2.length();
 1127 
 1128   if (n1 == -1) return p2;
 1129   if (n2 == -1) return p1;
 1130 
 1131   mem::vector<solvedKnot> nodes(n1+n2+1);
 1132 
 1133   Int i = 0;
 1134   nodes[0].pre = p1.point((Int) 0);
 1135   for (Int j = 0; j < n1; j++) {
 1136     nodes[i].point = p1.point(j);
 1137     nodes[i].straight = p1.straight(j);
 1138     nodes[i].post = p1.postcontrol(j);
 1139     nodes[i+1].pre = p1.precontrol(j+1);
 1140     i++;
 1141   }
 1142   for (Int j = 0; j < n2; j++) {
 1143     nodes[i].point = p2.point(j);
 1144     nodes[i].straight = p2.straight(j);
 1145     nodes[i].post = p2.postcontrol(j);
 1146     nodes[i+1].pre = p2.precontrol(j+1);
 1147     i++;
 1148   }
 1149   nodes[i].point = nodes[i].post = p2.point(n2);
 1150 
 1151   return path(nodes, i+1);
 1152 }
 1153 
 1154 // Interface to orient2d predicate optimized for pairs.
 1155 double orient2d(const pair& a, const pair& b, const pair& c)
 1156 {
 1157   double detleft, detright, det;
 1158   double detsum, errbound;
 1159   double orient;
 1160 
 1161   FPU_ROUND_DOUBLE;
 1162 
 1163   detleft = (a.getx() - c.getx()) * (b.gety() - c.gety());
 1164   detright = (a.gety() - c.gety()) * (b.getx() - c.getx());
 1165   det = detleft - detright;
 1166 
 1167   if (detleft > 0.0) {
 1168     if (detright <= 0.0) {
 1169       FPU_RESTORE;
 1170       return det;
 1171     } else {
 1172       detsum = detleft + detright;
 1173     }
 1174   } else if (detleft < 0.0) {
 1175     if (detright >= 0.0) {
 1176       FPU_RESTORE;
 1177       return det;
 1178     } else {
 1179       detsum = -detleft - detright;
 1180     }
 1181   } else {
 1182     FPU_RESTORE;
 1183     return det;
 1184   }
 1185 
 1186   errbound = ccwerrboundA * detsum;
 1187   if ((det >= errbound) || (-det >= errbound)) {
 1188     FPU_RESTORE;
 1189     return det;
 1190   }
 1191 
 1192   double pa[]={a.getx(),a.gety()};
 1193   double pb[]={b.getx(),b.gety()};
 1194   double pc[]={c.getx(),c.gety()};
 1195   
 1196   orient = orient2dadapt(pa, pb, pc, detsum);
 1197   FPU_RESTORE;
 1198   return orient;
 1199 }
 1200 
 1201 // Returns true iff the point z lies in or on the bounding box
 1202 // of a,b,c, and d.
 1203 bool insidebbox(const pair& a, const pair& b, const pair& c, const pair& d,
 1204                 const pair& z)
 1205 {
 1206   bbox B(a);
 1207   B.addnonempty(b);
 1208   B.addnonempty(c);
 1209   B.addnonempty(d);
 1210   return B.left <= z.getx() && z.getx() <= B.right && B.bottom <= z.gety() 
 1211     && z.gety() <= B.top;
 1212 }
 1213 
 1214 inline bool inrange(double x0, double x1, double x)
 1215 {
 1216   return (x0 <= x && x <= x1) || (x1 <= x && x <= x0);
 1217 }
 1218 
 1219 // Return true if point z is on z0--z1; otherwise compute contribution to 
 1220 // winding number.
 1221 bool checkstraight(const pair& z0, const pair& z1, const pair& z, Int& count)
 1222 {
 1223   if(z0.gety() <= z.gety() && z.gety() <= z1.gety()) {
 1224     double side=orient2d(z0,z1,z);
 1225     if(side == 0.0 && inrange(z0.getx(),z1.getx(),z.getx()))
 1226       return true;
 1227     if(z.gety() < z1.gety() && side > 0) ++count;
 1228   } else if(z1.gety() <= z.gety() && z.gety() <= z0.gety()) {
 1229     double side=orient2d(z0,z1,z);
 1230     if(side == 0.0 && inrange(z0.getx(),z1.getx(),z.getx()))
 1231       return true;
 1232     if(z.gety() < z0.gety() && side < 0) --count;
 1233   }
 1234   return false;
 1235 }
 1236 
 1237 // returns true if point is on curve; otherwise compute contribution to 
 1238 // winding number.
 1239 bool checkcurve(const pair& z0, const pair& c0, const pair& c1,
 1240                const pair& z1, const pair& z, Int& count, unsigned depth) 
 1241 {
 1242   if(depth == 0) return true;
 1243   --depth;
 1244 
 1245   if(insidebbox(z0,c0,c1,z1,z)) {
 1246     const pair m0=0.5*(z0+c0);
 1247     const pair m1=0.5*(c0+c1);
 1248     const pair m2=0.5*(c1+z1);
 1249     const pair m3=0.5*(m0+m1);
 1250     const pair m4=0.5*(m1+m2);
 1251     const pair m5=0.5*(m3+m4);
 1252     if(checkcurve(z0,m0,m3,m5,z,count,depth) || 
 1253        checkcurve(m5,m4,m2,z1,z,count,depth)) return true;
 1254   } else
 1255     if(checkstraight(z0,z1,z,count)) return true;
 1256   return false;
 1257 }
 1258 
 1259 // Return the winding number of the region bounded by the (cyclic) path
 1260 // relative to the point z, or the largest odd integer if the point lies on
 1261 // the path.
 1262 Int path::windingnumber(const pair& z) const
 1263 {
 1264   static const Int undefined=Int_MAX % 2 ? Int_MAX : Int_MAX-1;
 1265   
 1266   if(!cycles)
 1267     reportError("path is not cyclic");
 1268   
 1269   bbox b=bounds();
 1270   
 1271   if(z.getx() < b.left || z.getx() > b.right ||
 1272      z.gety() < b.bottom || z.gety() > b.top) return 0;
 1273   
 1274   Int count=0;
 1275   for(Int i=0; i < n; ++i)
 1276     if(straight(i)) {
 1277       if(checkstraight(point(i),point(i+1),z,count))
 1278         return undefined;
 1279     } else
 1280       if(checkcurve(point(i),postcontrol(i),precontrol(i+1),point(i+1),z,count,
 1281                     maxdepth)) return undefined;
 1282   return count;
 1283 }
 1284 
 1285 path path::transformed(const transform& t) const
 1286 {
 1287   mem::vector<solvedKnot> nodes(n);
 1288 
 1289   for (Int i = 0; i < n; ++i) {
 1290     nodes[i].pre = t * this->nodes[i].pre;
 1291     nodes[i].point = t * this->nodes[i].point;
 1292     nodes[i].post = t * this->nodes[i].post;
 1293     nodes[i].straight = this->nodes[i].straight;
 1294   }
 1295 
 1296   path p(nodes, n, cyclic());
 1297   return p;
 1298 }
 1299 
 1300 path transformed(const transform& t, const path& p)
 1301 {
 1302   Int n = p.size();
 1303   mem::vector<solvedKnot> nodes(n);
 1304 
 1305   for (Int i = 0; i < n; ++i) {
 1306     nodes[i].pre = t * p.precontrol(i);
 1307     nodes[i].point = t * p.point(i);
 1308     nodes[i].post = t * p.postcontrol(i);
 1309     nodes[i].straight = p.straight(i);
 1310   }
 1311 
 1312   return path(nodes, n, p.cyclic());
 1313 }
 1314 
 1315 path nurb(pair z0, pair z1, pair z2, pair z3,
 1316           double w0, double w1, double w2, double w3, Int m)
 1317 {
 1318   mem::vector<solvedKnot> nodes(m+1);
 1319 
 1320   if(m < 1) reportError("invalid sampling interval");
 1321 
 1322   double step=1.0/m;
 1323   for(Int i=0; i <= m; ++i) { 
 1324     double t=i*step;
 1325     double t2=t*t;
 1326     double onemt=1.0-t;
 1327     double onemt2=onemt*onemt;
 1328     double W0=w0*onemt2*onemt;
 1329     double W1=w1*3.0*t*onemt2;
 1330     double W2=w2*3.0*t2*onemt;
 1331     double W3=w3*t2*t;
 1332     nodes[i].point=(W0*z0+W1*z1+W2*z2+W3*z3)/(W0+W1+W2+W3);
 1333   }
 1334   
 1335   static const double twothirds=2.0/3.0;
 1336   pair z=nodes[0].point;
 1337   nodes[0].pre=z;
 1338   nodes[0].post=twothirds*z+third*nodes[1].point;
 1339   for(int i=1; i < m; ++i) {
 1340     pair z0=nodes[i].point;
 1341     pair zm=nodes[i-1].point;
 1342     pair zp=nodes[i+1].point;
 1343     pair pre=twothirds*z0+third*zm;
 1344     pair pos=twothirds*z0+third*zp;
 1345     pair dir=unit(pos-pre);
 1346     nodes[i].pre=z0-length(z0-pre)*dir;
 1347     nodes[i].post=z0+length(pos-z0)*dir;
 1348   }
 1349   z=nodes[m].point;
 1350   nodes[m].pre=twothirds*z+third*nodes[m-1].point;
 1351   nodes[m].post=z;
 1352   return path(nodes,m+1);
 1353 }
 1354 
 1355 } //namespace camp