Back to index

texmacs  1.0.7.15
curve.cpp
Go to the documentation of this file.
00001 
00002 /******************************************************************************
00003 * MODULE     : curve.cpp
00004 * DESCRIPTION: mathematical curves
00005 * COPYRIGHT  : (C) 2003  Joris van der Hoeven and Henri Lesourd
00006 *******************************************************************************
00007 * This software falls under the GNU general public license version 3 or later.
00008 * It comes WITHOUT ANY WARRANTY WHATSOEVER. For details, see the file LICENSE
00009 * in the root directory or <http://www.gnu.org/licenses/gpl-3.0.html>.
00010 ******************************************************************************/
00011 
00012 #include "path.hpp"
00013 #include "curve.hpp"
00014 #include "frame.hpp"
00015 #include "equations.hpp"
00016 #include "math_util.hpp"
00017 #include "polynomial.hpp"
00018 #include "merge_sort.hpp"
00019 
00020 /******************************************************************************
00021 * General routines
00022 ******************************************************************************/
00023 
00024 array<point>
00025 curve_rep::rectify (double eps) {
00026   array<point> a (1);
00027   a[0]= evaluate (0.0);
00028   rectify_cumul (a, eps);
00029   return a;
00030 }
00031 
00032 double
00033 curve_rep::bound (double t, double eps) {
00034   //TODO: Improve this, as soon as the curvature ()
00035   //      for transformed_curves will be implemented
00036   bool b;
00037   double ng= norm (grad (t, b));
00038   if (ng <= 1.0e-12)
00039     return tm_infinity;
00040   double delta= eps / ng;
00041   while (delta >= 1.0e-6) {
00042     point val1 = evaluate (t);
00043     point val2 = evaluate (max (t-delta, 0.0));
00044     point val3 = evaluate (min (t+delta, 1.0));
00045     if (norm (val2 - val1) <= eps+1.0e-6 &&
00046        norm (val3 - val1) <= eps+1.0e-6)
00047       return delta;
00048     delta /= 2.0;
00049   }
00050   return delta;
00051 }
00052 
00053 int
00054 curve_rep::get_control_points (
00055   array<double>&abs, array<point>& pts, array<path>& cip)
00056 {
00057   abs= array<double>();
00058   pts= array<point>();
00059   cip= array<path>();
00060   return 0;
00061 }
00062 
00063 struct curvet {
00064   double t, dist;
00065 };
00066 
00067 struct less_eq_curvet {
00068   static inline bool leq (curvet& a, curvet &b) {
00069     return a.dist <= b.dist; }
00070 };
00071 
00072 static array<curvet>
00073 curvet_closest_points (
00074   curve c, double t1, double t2, point p, double eps)
00075 {
00076   array<curvet> res;
00077   if (t1>t2) {
00078     array<curvet> r1, r2;
00079     r1= curvet_closest_points (c, 0.0, t2, p, eps);
00080     r2= curvet_closest_points (c, t1, 1.0, p, eps);
00081     res= r1 << r2;
00082   }
00083   else {
00084     double t;
00085     double closest= -1;
00086     point pclosest;
00087     double n0= tm_infinity;
00088     bool stored= true;
00089     double nprec= n0;
00090     bool decreasing= false;
00091     for (t=t1; t<=t2;) {
00092       point pt= c->evaluate(t);
00093       double n= norm (pt - p);
00094       if (n < n0) {
00095        n0= n;
00096        closest= t;
00097        pclosest= pt;
00098        stored= false;
00099       }
00100       decreasing= n<nprec;
00101       if (!stored && !decreasing) {
00102         curvet ct;
00103         ct.dist= norm (pclosest - p);
00104         ct.t= closest;
00105        res << ct;
00106        stored= true;
00107       }
00108       if (stored && decreasing)
00109         n0= tm_infinity;
00110       double delta= (n - eps) / 2;
00111       t += max (0.00001, c->bound (t, max (eps, delta)));
00112       nprec= n;
00113     }
00114     if (!stored && decreasing) {
00115       curvet ct;
00116       ct.dist= norm (pclosest - p);
00117       ct.t= closest;
00118       res << ct;
00119     }
00120   }
00121   return res;
00122 }
00123 
00124 array<double>
00125 curve_rep::find_closest_points (
00126   double t1, double t2, point p, double eps)
00127 {
00128   array<curvet> res= curvet_closest_points (curve (this), t1, t2, p, eps);
00129   merge_sort_leq <curvet, less_eq_curvet> (res);
00130   array<double> rest= array<double> (N(res));
00131   int i;
00132   for (i=0; i<N(res); i++)
00133     rest[i]= res[i].t;
00134   return rest;
00135 }
00136 
00137 double
00138 curve_rep::find_closest_point (
00139   double t1, double t2, point p, double eps, bool& found)
00140 {
00141   array<double> res= find_closest_points (t1, t2, p, eps);
00142   found= N(res)>0;
00143   if (found)
00144     return res[0];
00145   else
00146     return -1;
00147 }
00148 
00149 point
00150 closest (curve f, point p) {
00151   array<double> abs;
00152   array<point> pts;
00153   array<path> rcip;
00154   f->get_control_points (abs, pts, rcip);
00155   if (N(abs) == 0) return f(0);
00156   double t1= abs[0];
00157   double t2= abs[N(abs)-1];
00158   double best= 0;
00159   double eps = norm (f(0) - p);
00160   for (int i=0; i<10; i++) {
00161     bool found= false;
00162     double t= f->find_closest_point (t1, t2, p, eps, found);
00163     if (found) best= t;
00164     else break;
00165     double eps2= norm (f(t) - p);
00166     if (eps2 >= 0.9 * eps) break;
00167     eps= eps2;
00168   }
00169   return f(best);
00170 }
00171 
00172 point
00173 grad (curve f, double t) {
00174   bool error= false;
00175   return f->grad (t, error);
00176 }
00177 
00178 bool
00179 intersection (curve f, curve g, double& t, double& u) {
00180   // for two dimensional curves only
00181   double d= norm (f (t) - g (u));
00182   while (!fnull (d, 1.0e-9)) {
00183     point  ft = f (t);
00184     point  gu = g (u);
00185     point  gf = grad (f, t);
00186     point  gg = grad (g, u);
00187     double det= gf[0] * gg[1] - gf[1] * gg[0];
00188     if (fnull (det, 1.0e-12)) break;
00189     double dt = ((gu[0] - ft[0]) * gg[1] - (gu[1] - ft[1]) * gg[0]) / det;
00190     double du = ((gu[0] - ft[0]) * gf[1] - (gu[1] - ft[1]) * gf[0]) / det;
00191     double T  = t + dt;
00192     double U  = u + du;
00193     if (T < 0.0 || T > 1.0 || U < 0.0 || U > 1.0) break;
00194     double D  = norm (f (T) - g (U));
00195     if (D > 0.9 * d) break;
00196     t= T; u= U; d= D;
00197   }
00198   return fnull (d, 1.0e-9);
00199 }
00200 
00201 array<point>
00202 intersection (curve f, curve g, point p0, double eps) {
00203   // For local intersections only
00204   array<point> res;
00205   if (is_nil (f) || is_nil (g)) return res;
00206   bool found= false;
00207   double d1, d2;
00208   if (f==g) {
00209     array<double> pts= f->find_closest_points (0.0, 1.0, p0, eps);
00210     if (N(pts)>=2) {
00211       d1= pts[0];
00212       d2= pts[1];
00213       found= intersection (f, f, d1, d2);
00214     }
00215   }
00216   else {
00217     bool res1, res2;
00218     d1= f->find_closest_point (0.0, 1.0, p0, eps, res1);
00219     d2= g->find_closest_point (0.0, 1.0, p0, eps, res2);
00220     if (res1 && res2)
00221       found= intersection (f, g, d1, d2);
00222   }
00223   if (found)
00224     res << f (d1);
00225   return res;
00226 }
00227 
00228 /******************************************************************************
00229 * Segments
00230 ******************************************************************************/
00231 
00232 struct segment_rep: public curve_rep {
00233   point p1, p2;
00234   path cip1, cip2;
00235   segment_rep (point p1b, point p2b): p1 (p1b), p2 (p2b) {}
00236   point evaluate (double t) { return (1.0-t)*p1 + t*p2; }
00237   void rectify_cumul (array<point>& a, double eps) { (void) eps; a << p2; }
00238   double bound (double t, double eps) {
00239     return curve_rep::bound (t, eps);
00240   }
00241   point grad (double t, bool& error) {
00242     (void) t;
00243     error= false;
00244     return p2 - p1;
00245   }
00246   double curvature (double t1, double t2) {
00247     (void) t1; (void) t2;
00248     return tm_infinity;
00249   }
00250   int get_control_points (
00251     array<double>&abs, array<point>& pts, array<path>& cip);
00252 };
00253 
00254 int
00255 segment_rep::get_control_points (
00256   array<double>&abs, array<point>& pts, array<path>& cip)
00257 {
00258   array<double> u;
00259   u << 0.0;
00260   u << 1.0;
00261   abs= u;
00262   array<point> a;
00263   a << p1;
00264   a << p2;
00265   pts= a;
00266   cip= array<path> ();
00267   cip << cip1;
00268   cip << cip2;
00269   return 2;
00270 }
00271 
00272 curve
00273 segment (point p1, point p2) {
00274   return tm_new<segment_rep> (p1, p2);
00275 }
00276 
00277 /******************************************************************************
00278 * Poly-segments
00279 ******************************************************************************/
00280 
00281 struct poly_segment_rep: public curve_rep {
00282   array<point> a;
00283   array<path> cip;
00284   int n;
00285   poly_segment_rep (array<point> a2, array<path> cip2):
00286     a (a2), cip (cip2), n(N(a)-1) {}
00287   int nr_components () { return n; }
00288   point evaluate (double t) {
00289     int i= min ((int) (n*t), n-1);
00290     return (i+1 - n*t)*a[i] + (n*t - i)*a[i+1];
00291   }
00292   void rectify_cumul (array<point>& cum, double eps) {
00293     (void) eps;
00294     for (int i=0; i<n; i++)
00295       cum << a[i+1];
00296   }
00297   double bound (double t, double eps) {
00298     return curve_rep::bound (t, eps);
00299   }
00300   double curvature (double t1, double t2) {
00301     (void) t1; (void) t2;
00302     return tm_infinity;
00303   }
00304   point grad (double t, bool& error) {
00305     error= false;
00306     int i= min ((int) (n*t), n-1);
00307     return n * (a[i+1] - a[i]);
00308   }
00309   int get_control_points (
00310     array<double>&abs, array<point>& pts, array<path>& cip);
00311 };
00312 
00313 int
00314 poly_segment_rep::get_control_points (
00315   array<double>&abs, array<point>& pts, array<path>& rcip)
00316 {
00317   array<double> u(n+1);
00318   int i;
00319   for (i=0; i<n+1; i++)
00320     u[i]= ((double)i) / (n==0?1:n);
00321   abs = u;
00322   pts = a;
00323   rcip= cip;
00324   return n+1;
00325 }
00326 
00327 curve
00328 poly_segment (array<point> a, array<path> cip) {
00329   return tm_new<poly_segment_rep> (a, cip);
00330 }
00331 
00332 /******************************************************************************
00333 * Splines
00334 ******************************************************************************/
00335 
00336 static const double epsilon=0.01;//0.00005;
00337 typedef polynomial<double> dpol;
00338 typedef vector<polynomial<double> > dpols;
00339 
00340 struct spline_rep: public curve_rep {
00341   array<point> a;
00342   array<path> cip;
00343   int n;
00344   array<double> U;
00345   array<dpols> p;
00346   bool close, interpol;
00347 
00348   spline_rep (
00349     array<point> a, array<path> cip, bool close=false, bool interpol=true);
00350 
00351   inline double d (int i,int k) { return U[i]-U[i-k]; }
00352   inline double m (int i) { return (U[i]+U[i+1])/2; }
00353 
00354   double convert (double u) { return U[2]+u*(U[n+1]-U[2]); }
00355   double unconvert (double u) { return (u-U[2])/(U[n+1]-U[2]); }
00356   int interval_no (double u);
00357 
00358   point spline (int i,double u,int o=0);
00359   inline double S (
00360     array<dpol> p1, array<dpol> p2, array<dpol> p3,
00361     int i, double u);
00362   point evaluate (double t,int o);
00363   point evaluate (double t);
00364   double bound (double t, double eps);
00365   point grad (double t, bool& error);
00366 
00367   double curvature (int i, double t1, double t2);
00368   double curvature (double t1, double t2);
00369 
00370   bool approx (int i, double u1, double u2, double eps);
00371   void rectify_cumul (array<point>& cum, int i,
00372                     double u1, double u2, double eps);
00373   void rectify_cumul (array<point>& cum, double eps);
00374   int get_control_points (
00375     array<double>&abs, array<point>& pts, array<path>& cip);
00376 };
00377 
00378 // Creation
00379 spline_rep::spline_rep (
00380   array<point> a2, array<path> cip2, bool close2,bool interpol2):
00381   a(a2), cip(cip2), n(N(a)-1), close(close2), interpol(interpol2)
00382 {
00383   array<dpol> p1,p2,p3;
00384   p1= array<dpol> (n+3);
00385   p2= array<dpol> (n+3);
00386   p3= array<dpol> (n+3);
00387   p = array<dpols> (n+3);
00388   U = array<double> (n+6);
00389   if (close) {
00390     int i;
00391     a->resize(N(a)+2);
00392     for (i=0;i<=1;i++)
00393       a[n+1+i]=a[i];
00394     n+=2;
00395   }
00396   int i;
00397   double x=0;
00398   if (!close) {
00399     for (i=0;i<3;i++) {
00400       U[i]=x;
00401       x+=epsilon; 
00402     }
00403     x+=1-epsilon;
00404     for (i=3;i<=n;i++) {
00405       U[i]=x;
00406       x+=1;
00407     }
00408     for (i=n+1;i<=n+3;i++) {
00409       U[i]=x;
00410       x+=epsilon;
00411     }
00412   }
00413   else {
00414     for (i=0;i<=n+3;i++) {
00415       U[i]=x;
00416       x+=1;
00417     }
00418   }
00419   for (i=0;i<=n;i++) {
00420     double di22= d(i+2,2);
00421     double di11= d(i+1,1);
00422     double di21= d(i+2,1);
00423     double di32= d(i+3,2);
00424     double di31= d(i+3,1);
00425     p1[i]= dpol (square(U[i])/di22/di11,
00426                -2*U[i]/di22/di11,
00427                1/di22/di11);
00428     p2[i]= dpol (-U[i+2]*U[i]/di22/di21-U[i+3]*U[i+1]/di32/di21,
00429                (U[i+2]+U[i])/di22/di21+(U[i+3]+U[i+1])/di32/di21,
00430                -1/di22/di21-1/di32/di21);
00431     p3[i]= dpol (square(U[i+3])/di32/di31,
00432                -2*U[i+3]/di32/di31,
00433                1/di32/di31);
00434     /*
00435     p1[i][2]= 1/di22/di11;
00436     p1[i][1]= -2*U[i]/di22/di11;
00437     p1[i][0]= square(U[i])/di22/di11;
00438     p2[i][2]= -1/di22/di21-1/di32/di21;
00439     p2[i][1]= (U[i+2]+U[i])/di22/di21+(U[i+3]+U[i+1])/di32/di21;
00440     p2[i][0]= -U[i+2]*U[i]/di22/di21-U[i+3]*U[i+1]/di32/di21;
00441     p3[i][2]= 1/di32/di31;
00442     p3[i][1]= -2*U[i+3]/di32/di31;
00443     p3[i][0]= square(U[i+3])/di32/di31;
00444     */
00445   }
00446   if (interpol) {
00447     array<point> x(n+1), y(n+1);
00448     y= a;
00449     {
00450       array<double> a(n+1), b(n+1), c(n+1);
00451       if (close) {
00452         a[n-2]= S (p1, p2, p3, n-3, m(n-1));
00453         b[0]= S (p1, p2, p3, 1, m(2));
00454         b[n-2]= S (p1, p2, p3, n-2, m(n-1));
00455         c[0]=S (p1, p2, p3, 2, m(2));
00456         for (i=1; i<n-2; i++) {
00457            a[i]= S (p1, p2, p3, i, m(i+2));
00458            b[i]= S (p1, p2, p3, i+1, m(i+2));
00459            c[i]= S (p1, p2, p3, i+2, m(i+2));
00460         }
00461         xtridiag_solve (a, b, c, S(p1, p2, p3, n-1, m(n-1)),
00462                                  S(p1, p2, p3, 0, m(2)),
00463                                  x, y, n-1);
00464      /* Rotate the result in order to have an appropriate
00465         parametrization (to keep in sync the ordering of
00466         the points & the paths in cip, in particular) */
00467         point p=x[n-2];
00468         for (i=n-2;i>=1;i--)
00469           x[i]=x[i-1];
00470         x[0]=p;
00471      // Splice the spline
00472         x[n-1]= x[0];
00473         x[n]= x[1];
00474       }
00475       else {
00476         a[n]= S (p1, p2, p3, n-1, U[n+1]);
00477         b[0]= S (p1, p2, p3, 0, U[2]);
00478         b[n]= S (p1, p2, p3, n, U[n+1]);
00479         c[0]= S (p1, p2, p3, 1, U[2]);
00480         for (i=1; i<n; i++) {
00481            a[i]= S (p1, p2, p3, i-1, m(i+1));
00482            b[i]= S (p1, p2, p3, i, m(i+1));
00483            c[i]= S (p1, p2, p3, i+1, m(i+1));
00484         }
00485         tridiag_solve (a, b, c, x, y, n+1);
00486       }
00487     }
00488     a= x;
00489   }
00490   for (i=2; i<=n; i++)
00491     p[i]= a[i]*p1[i] + a[i-1]*p2[i-1] + a[i-2]*p3[i-2];
00492 }
00493 
00494 // Evaluation
00495 double
00496 spline_rep::S (
00497   array<dpol> p1, array<dpol> p2, array<dpol> p3,
00498   int i, double u)
00499 {
00500   if (i<0 || i>n) return 0;
00501   else if (u<U[i] || u>=U[i+3]) return 0;
00502   else if (u<U[i+1]) return p1[i](u);
00503   else if (u<U[i+2]) return p2[i](u);
00504   else if (u<U[i+3]) return p3[i](u);
00505   else FAILED ("we should **never** go here");
00506   return 0.0;
00507 }
00508 
00509 point
00510 spline_rep::spline (int i, double u, int o) {
00511   int j, n= N(p[i]);
00512   point res (n);
00513   if (o<0 || o>2) o=0;
00514   for (j=0; j<n; j++)
00515     res[j]= p[i][j] (u, o);
00516   return res;
00517 }
00518 
00519 int
00520 spline_rep::interval_no (double u) {
00521   int i;
00522   for (i=0;i<N(U);i++)
00523     if (u>=U[i] && u<U[i+1]) return i;
00524   return -1;
00525 }
00526 
00527 static double
00528 prod (double x, int n) {
00529   double r;
00530   if (n<=0) return 1;
00531   r=x; n--;
00532   while (n--)
00533     r*=(x-n);
00534   return r;
00535 }
00536 
00537 point
00538 spline_rep::evaluate (double t, int o) {
00539   point res;
00540   int no;
00541   t=convert(t);
00542   double k=U[n+1]-U[2];
00543   no=interval_no(t);
00544   if (no<2) res=spline(2,U[2],o);
00545   else if (no>n) res=spline(n,U[n+1],o);
00546   else res=spline(no,t,o);
00547   return prod(k,o)*res;
00548 }
00549 
00550 point
00551 spline_rep::evaluate (double t) {
00552   return evaluate(t,0);
00553 }
00554 
00555 double
00556 spline_rep::bound (double t, double eps) {
00557   return eps/norm(evaluate(t,1));
00558 }
00559 
00560 point
00561 spline_rep::grad (double t, bool& error) {
00562   error= false;
00563   return evaluate(t,1);
00564 }
00565 
00566 // Rectification
00567 bool
00568 spline_rep::approx (int i, double u1, double u2, double eps) {
00569   double l,R;
00570   point p1,p2;
00571   p1=spline(i,u1);
00572   p2=spline(i,u2);
00573   l=norm(p1-p2);
00574   // When l and R are very small, the test l<=R
00575   // can fail forever. So we set l to exactly 0
00576   if (l!=0 && fnull (l,1.0e-6)) l=0;
00577   R=curvature(i,u1,u2);
00578   return l<=2*sqrt(2*R*eps);
00579 }
00580 
00581 void
00582 spline_rep::rectify_cumul (array<point>& cum, int i,
00583                            double u1, double u2, double eps) {
00584   if (approx(i,u1,u2,eps))
00585     cum << spline(i,u2);
00586   else {
00587     double u=(u1+u2)/2;
00588     rectify_cumul(cum,i,u1,u,eps);
00589     rectify_cumul(cum,i,u,u2,eps);
00590   }
00591 }
00592 
00593 void
00594 spline_rep::rectify_cumul (array<point>& cum, double eps) {
00595   int i;
00596   for (i=2;i<=n;i++)
00597     rectify_cumul(cum,i,U[i],U[i+1],eps);
00598 }
00599 
00600 // Curvature
00601 double
00602 spline_rep::curvature (int i, double t1, double t2) {
00603   point a, b;
00604   a= extract (p[i], 2);
00605   b= extract (p[i], 1);
00606   double t,R;
00607   point pp,ps;
00608   if (norm(a)==0) return tm_infinity;
00609   t=-(a*b)/(2*a*a);
00610   if (t1>t) t=t1;
00611   else if (t2<t) t=t2;
00612   pp=spline(i,t,1);
00613   ps=spline(i,t,2);
00614   if (norm(ps)==0) return tm_infinity;
00615   R=square(norm(pp))/norm(ps);
00616   return R;
00617 }
00618 
00619 double
00620 spline_rep::curvature (double t1, double t2) {
00621   double res;
00622   int no,no1,no2;
00623   t1=convert(t1);
00624   t2=convert(t2);
00625   no1=interval_no(t1);
00626   no2=interval_no(t2);
00627   res=tm_infinity;
00628   for (no=no1;no<=no2;no++) {
00629     double r;
00630     if (no<2) r=curvature(2,t1,t2);
00631     else if (no>n) r=curvature(n,t1,t2);
00632     else r=curvature(no,t1,t2);
00633     if (r<=res) res=r;
00634   }
00635   return res;
00636 }
00637 
00638 // Control points
00639 int
00640 spline_rep::get_control_points (
00641   array<double>&abs, array<point>& pts, array<path>& rcip)
00642 {
00643   array<double> u(n+1);
00644   array<point> p(n+1);
00645   int i;
00646   if (interpol) {
00647     if (close) {
00648       u->resize (n-1);
00649       p->resize (n-1);
00650       for (i=0; i<n-1; i++)
00651        u[i]= unconvert ((U[i+2]+U[i+3])/2);
00652     }
00653     else {
00654       u[0]= unconvert (U[2]);
00655       u[n]= unconvert (U[n+1]);
00656       for (i=2; i<=n; i++)
00657        u[i-1]= unconvert ((U[i]+U[i+1])/2);
00658     }
00659     for (i=0; i<=(close ? n-2 : n); i++)
00660       p[i]= evaluate (u[i]);
00661   }
00662   else FAILED ("not yet implemented");
00663   abs = u;
00664   pts = p;
00665   rcip= cip;
00666   return close ? n-1 : n+1;
00667 }
00668 
00669 curve
00670 spline (array<point> a, array<path> cip, bool close, bool interpol) {
00671   return tm_new<spline_rep> (a, cip, close, interpol);
00672 }
00673 
00674 /******************************************************************************
00675 * Arcs
00676 ******************************************************************************/
00677 
00678 struct arc_rep: public curve_rep {
00679   array<point> a;
00680   array<path> cip;
00681   array<double> u;
00682   point center;  // Center
00683   point i, j;    // The two base vectors of the ellipsis's 2D plane
00684   double r1, r2; // The two radiuses of the ellipsis
00685   double e1, e2; // Coordinates of the two extremal points of the arc
00686   arc_rep (array<point> a, array<path> cip, bool close);
00687   point evaluate (double t);
00688   void rectify_cumul (array<point>& cum, double eps);
00689   double bound (double t, double eps);
00690   point grad (double t, bool& error);
00691   double curvature (double t1, double t2);
00692   int get_control_points (
00693     array<double>&abs, array<point>& pts, array<path>& cip);
00694 };
00695 
00696 arc_rep::arc_rep (array<point> a2, array<path> cip2, bool close):
00697   a(a2), cip(cip2)
00698 {
00699   int n= N(a);
00700   point o1= (n>0 ? a[0] : point (0,0));
00701   point o2= (n>1 ? a[1] : point (0,0));
00702   point o3= (n>2 ? a[2] : point (0,0));
00703   if (n!=3 || linearly_dependent (o1, o2, o3) ||
00704      (N (center= intersection (midperp (o1, o2, o3),
00705                             midperp (o2, o3, o1))) == 0))
00706   {
00707     center= i= j= 0*o1;
00708     r1= r2= 1;
00709     e1= 0;
00710     e2= 1;
00711     a= array<point> (1);
00712     a[0]= o1;
00713     if (N(cip)) {
00714       path p= cip[0];
00715       cip= array<path> (1);
00716       cip[0]= p;
00717     }
00718     u= array<double> (1);
00719     u[0]= 0;
00720     return;
00721   }
00722   r1= r2= norm (o1-center);
00723   if (orthogonalize (i, j, center, o1, o2)) ;
00724   else
00725     orthogonalize (i, j, center, o1, o3);
00726   e1= 0;
00727   point o2b (2), o3b (2);
00728   o2b[0]= (o2-center) * i;
00729   o2b[1]= (o2-center) * j;
00730   o3b[0]= (o3-center) * i;
00731   o3b[1]= (o3-center) * j;
00732   e2= arg (o3b) / (2*tm_PI);
00733   double e3= arg (o2b) / (2*tm_PI);
00734   if (e2<e3) {
00735     j= -j;
00736     o3b[1]= -o3b[1];
00737     e2= arg (o3b) / (2*tm_PI);
00738     o2b[1]= -o2b[1];
00739     e3= arg (o2b) / (2*tm_PI);
00740   }
00741   u= array<double> (3);
00742   u[0]= 0;
00743   u[1]= e3;
00744   u[2]= e2;
00745   if (close) e2= 1.0;
00746 }
00747 
00748 point
00749 arc_rep::evaluate (double t) {
00750   t= e1 + t*(e2 - e1);
00751   return center + r1*cos(2*tm_PI*t)*i
00752                 + r2*sin(2*tm_PI*t)*j;
00753 }
00754 
00755 void
00756 arc_rep::rectify_cumul (array<point>& cum, double eps) {
00757   double t, step;
00758   step= sqrt (2*eps / max (r1, r2) ) / tm_PI;
00759   for (t=step; t<=1.0; t+=step)
00760     cum << evaluate (t);
00761   if (t-step != 1.0)
00762     cum << evaluate (1.0);
00763 }
00764 
00765 double
00766 arc_rep::bound (double t, double eps) {
00767   return curve_rep::bound (t, eps);
00768 }
00769 
00770 point
00771 arc_rep::grad (double t, bool& error) {
00772   error= false;
00773   t= e1 + t*(e2 - e1);
00774   return -2*tm_PI*r1*sin(2*tm_PI*t)*i
00775          +2*tm_PI*r2*cos(2*tm_PI*t)*j;
00776 }
00777 
00778 double
00779 arc_rep::curvature (double t1, double t2) {
00780   (void) t1; (void) t2;
00781   if (r1 >= r2)
00782     return fnull (r1,1.0e-6) ? tm_infinity : square(r2)/r1;
00783   else
00784     return fnull (r2,1.0e-6) ? tm_infinity : square(r1)/r2;
00785 }
00786 
00787 curve
00788 arc (array<point> a, array<path> cip, bool close) {
00789   return tm_new<arc_rep> (a, cip, close);
00790 }
00791 
00792 int
00793 arc_rep::get_control_points (
00794   array<double>&abs, array<point>& pts, array<path>& rcip)
00795 {
00796   abs = u;
00797   pts = a;
00798   rcip= cip;
00799   return N(a);
00800 }
00801 
00802 /******************************************************************************
00803 * Compound curves
00804 ******************************************************************************/
00805 
00806 struct compound_curve_rep: public curve_rep {
00807   curve c1, c2;
00808   int n1, n2;
00809   compound_curve_rep (curve c1b, curve c2b):
00810     c1 (c1b), c2 (c2b), n1 (c1->nr_components()), n2 (c2->nr_components()) {}
00811   int nr_components () { return n1 + n2; }
00812   point evaluate (double t) {
00813     double n= n1+n2;
00814     if (t <= n1/n) return c1 (t*n/n1);
00815     else return c2 (t*n/n2 - n1); }
00816   void rectify_cumul (array<point>& a, double eps) {
00817     c1->rectify_cumul (a, eps);
00818     c2->rectify_cumul (a, eps);
00819   }
00820   double bound (double t, double eps) {
00821     return curve_rep::bound (t, eps);
00822   }
00823   point grad (double t, bool& error) {
00824     double n= n1+n2;
00825     if (t <= n1/n) return (n * c1->grad ((t*n)/n1, error)) / n1;
00826     else return (n * c2->grad ((t*n)/n2 - n1, error)) / n2;
00827   }
00828   double curvature (double t1, double t2) {
00829     return max (c1->curvature (t1, t2), c2->curvature (t1, t2));
00830   }
00831   /*int get_control_points (
00832     array<double>&abs, array<point>& pts, array<path>& cip);
00833   */
00834 };
00835 
00836 curve
00837 operator * (curve c1, curve c2) {
00838   // FIXME: we might want to test whether c1(1.0) is approx equal to c2(0.0)
00839   return tm_new<compound_curve_rep> (c1, c2);
00840 }
00841 
00842 /******************************************************************************
00843 * Inverted curves
00844 ******************************************************************************/
00845 
00846 struct inverted_curve_rep: public curve_rep {
00847   curve c;
00848   int n;
00849   inverted_curve_rep (curve c2): c (c2), n (c->nr_components()) {}
00850   int nr_components () { return n; }
00851   point evaluate (double t) { return c (1.0 - t); }
00852   void rectify_cumul (array<point>& a, double eps) {
00853     array<point> b= c->rectify (eps);
00854     int i, k= N(b);
00855     for (i=k-1; i>=0; i--) a << b[i];
00856   }
00857   double bound (double t, double eps) {
00858     return curve_rep::bound (t, eps);
00859   }
00860   point grad (double t, bool& error) {
00861     return - c->grad (1.0 - t, error);
00862   }
00863   double curvature (double t1, double t2) {
00864     return c->curvature (1-t2, 1-t1);
00865   }
00866   int get_control_points (
00867     array<double>&abs, array<point>& pts, array<path>& cip);
00868 };
00869 
00870 int
00871 inverted_curve_rep::get_control_points (
00872   array<double>&abs, array<point>& pts, array<path>& cip)
00873 {
00874   int res= c->get_control_points (abs, pts, cip);
00875   int i;
00876   abs= copy (abs);
00877   for (i=0; i<res; i++)
00878     abs[i]= 1 - abs[i];
00879   return res;
00880 }
00881 
00882 curve
00883 invert (curve c) {
00884   return tm_new<inverted_curve_rep> (c);
00885 }
00886 
00887 /******************************************************************************
00888 * Transformed curves
00889 ******************************************************************************/
00890 
00891 struct transformed_curve_rep: public curve_rep {
00892   frame f;
00893   curve c;
00894   int n;
00895   transformed_curve_rep (frame f2, curve c2):
00896     f (f2), c (c2), n (c->nr_components()) {}
00897   int nr_components () { return n; }
00898   point evaluate (double t) { return f (c (t)); }
00899   void rectify_cumul (array<point>& a, double eps);
00900   double bound (double t, double eps) {
00901     return curve_rep::bound (t, eps);
00902   }
00903   point grad (double t, bool& error);
00904   double curvature (double t1, double t2) {
00905     (void) t1; (void) t2;
00906     FAILED ("not yet implemented");
00907     return 0.0;
00908   }
00909   int get_control_points (
00910     array<double>&abs, array<point>& pts, array<path>& cip);
00911 };
00912 
00913 void
00914 transformed_curve_rep::rectify_cumul (array<point>& a, double eps) {
00915   if (f->linear) {
00916     double delta= f->direct_bound (c(0.0), eps);
00917     array<point> b= c->rectify (delta);
00918     int i, k= N(b);
00919     for (i=0; i<k; i++) a << f(b[i]);
00920   }
00921   else FAILED ("not yet implemented");
00922 }
00923 
00924 point
00925 transformed_curve_rep::grad (double t, bool& error) {
00926   bool error2;
00927   point w2= c->grad (t, error2);
00928   point w1= f->jacobian (c(t), w2, error);
00929   error |= error2;
00930   return w1;
00931 }
00932 
00933 int
00934 transformed_curve_rep::get_control_points (
00935   array<double>&abs, array<point>& pts, array<path>& cip)
00936 {
00937   int res= c->get_control_points (abs, pts, cip);
00938   int i;
00939   pts= copy (pts);
00940   for (i=0; i<N(pts); i++)
00941     pts[i]= f (pts[i]);
00942   return res;
00943 }
00944 
00945 curve
00946 frame::operator () (curve c) {
00947   return tm_new<transformed_curve_rep> (*this, c);
00948 }
00949 
00950 curve
00951 frame::operator [] (curve c) {
00952   return tm_new<transformed_curve_rep> (invert (*this), c);
00953 }