src/livarot/PathSimplify.cpp

   1 /*
   2  *  PathSimplify.cpp
   3  *  nlivarot
   4  *
   5  *  Created by fred on Fri Dec 12 2003.
   6  *
   7  */
   8
   9 #include <glib/gmem.h>
  10 #include <libnr/nr-point-matrix-ops.h>
  11 #include "livarot/Path.h"
  12 #include "livarot/path-description.h"
  13
  14 /*
  15  * Reassembling polyline segments into cubic bezier patches
  16  * thes functions do not need the back data. but they are slower than recomposing
  17  * path descriptions when you have said back data (it's always easier with a model)
  18  * there's a bezier fitter in bezier-utils.cpp too. the main difference is the way bezier patch are split
  19  * here: walk on the polyline, trying to extend the portion you can fit by respecting the treshhold, split when
  20  * treshhold is exceeded. when encountering a "forced" point, lower the treshhold to favor splitting at that point
  21  * in bezier-utils: fit the whole polyline, get the position with the higher deviation to the fitted curve, split
  22  * there and recurse
  23  */
  24
  25
  26 // algo d'origine: http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/INT-APP/CURVE-APP-global.html
  27
  28 // need the b-spline basis for cubic splines
  29 // pas oublier que c'est une b-spline clampee
  30 // et que ca correspond a une courbe de bezier normale
  31 #define N03(t) ((1.0-t)*(1.0-t)*(1.0-t))
  32 #define N13(t) (3*t*(1.0-t)*(1.0-t))
  33 #define N23(t) (3*t*t*(1.0-t))
  34 #define N33(t) (t*t*t)
  35 // quadratic b-splines (jsut in case)
  36 #define N02(t) ((1.0-t)*(1.0-t))
  37 #define N12(t) (2*t*(1.0-t))
  38 #define N22(t) (t*t)
  39 // linear interpolation b-splines
  40 #define N01(t) ((1.0-t))
  41 #define N11(t) (t)
  42
  43
  44
  45 void Path::Simplify(double treshhold)
  46 {
  47     if (pts.size() <= 1) {
  48         return;
  49     }
  50
  51     Reset();
  52
  53     int lastM = 0;
  54     while (lastM < int(pts.size())) {
  55         int lastP = lastM + 1;
  56         while (lastP < int(pts.size())
  57                && (pts[lastP].isMoveTo == polyline_lineto
  58                    || pts[lastP].isMoveTo == polyline_forced))
  59         {
  60             lastP++;
  61         }
  62
  63         DoSimplify(lastM, lastP - lastM, treshhold);
  64
  65         lastM = lastP;
  66     }
  67 }
  68
  69
  70 // dichomtomic method to get distance to curve approximation
  71 // a real polynomial solver would get the minimum more efficiently, but since the polynom
  72 // would likely be of degree >= 5, that would imply using some generic solver, liek using the sturm metod
  73 double RecDistanceToCubic(NR::Point const &iS, NR::Point const &isD,
  74                           NR::Point const &iE, NR::Point const &ieD,
  75                           NR::Point &pt, double current, int lev, double st, double et)
  76 {
  77     if ( lev <= 0 ) {
  78         return current;
  79     }
  80
  81     NR::Point const m = 0.5 * (iS + iE) + 0.125 * (isD - ieD);
  82     NR::Point const md = 0.75 * (iE - iS) - 0.125 * (isD + ieD);
  83     double const mt = (st + et) / 2;
  84
  85     NR::Point const hisD = 0.5 * isD;
  86     NR::Point const hieD = 0.5 * ieD;
  87
  88     NR::Point const mp = pt - m;
  89     double nle = NR::dot(mp, mp);
  90
  91     if ( nle < current ) {
  92
  93         current = nle;
  94         nle = RecDistanceToCubic(iS, hisD, m, md, pt, current, lev - 1, st, mt);
  95         if ( nle < current ) {
  96             current = nle;
  97         }
  98         nle = RecDistanceToCubic(m, md, iE, hieD, pt, current, lev - 1, mt, et);
  99         if ( nle < current ) {
 100             current = nle;
 101         }
 102
 103     } else if ( nle < 2 * current ) {
 104
 105         nle = RecDistanceToCubic(iS, hisD, m, md, pt, current, lev - 1, st, mt);
 106         if ( nle < current ) {
 107             current = nle;
 108         }
 109         nle = RecDistanceToCubic(m, md, iE, hieD, pt, current, lev - 1, mt, et);
 110         if ( nle < current ) {
 111             current = nle;
 112         }
 113     }
 114
 115     return current;
 116 }
 117
 118
 119 double DistanceToCubic(NR::Point const &start, PathDescrCubicTo res, NR::Point &pt)
 120 {
 121     NR::Point const sp = pt - start;
 122     NR::Point const ep = pt - res.p;
 123     double nle = NR::dot(sp, sp);
 124     double nnle = NR::dot(ep, ep);
 125     if ( nnle < nle ) {
 126         nle = nnle;
 127     }
 128
 129     NR::Point seg = res.p - start;
 130     nnle = NR::cross(seg, sp);
 131     nnle *= nnle;
 132     nnle /= NR::dot(seg, seg);
 133     if ( nnle < nle ) {
 134         if ( NR::dot(sp,seg) >= 0 ) {
 135             seg = start - res.p;
 136             if ( NR::dot(ep,seg) >= 0 ) {
 137                 nle = nnle;
 138             }
 139         }
 140     }
 141
 142     return nle;
 143 }
 144
 145
 146 /**
 147  *    Simplification on a subpath.
 148  */
 149
 150 void Path::DoSimplify(int off, int N, double treshhold)
 151 {
 152   // non-dichotomic method: grow an interval of points approximated by a curve, until you reach the treshhold, and repeat
 153     if (N <= 1) {
 154         return;
 155     }
 156
 157     int curP = 0;
 158
 159     fitting_tables data;
 160     data.Xk = data.Yk = data.Qk = NULL;
 161     data.tk = data.lk = NULL;
 162     data.fk = NULL;
 163     data.totLen = 0;
 164     data.nbPt = data.maxPt = data.inPt = 0;
 165
 166     NR::Point const moveToPt = pts[off].p;
 167     MoveTo(moveToPt);
 168     NR::Point endToPt = moveToPt;
 169
 170     while (curP < N - 1) {
 171
 172         int lastP = curP + 1;
 173         int M = 2;
 174
 175         // remettre a zero
 176         data.inPt = data.nbPt = 0;
 177
 178         PathDescrCubicTo res(NR::Point(0, 0), NR::Point(0, 0), NR::Point(0, 0));
 179         bool contains_forced = false;
 180         int step = 64;
 181
 182         while ( step > 0 ) {
 183             int forced_pt = -1;
 184             int worstP = -1;
 185
 186             do {
 187                 if (pts[off + lastP].isMoveTo == polyline_forced) {
 188                     contains_forced = true;
 189                 }
 190                 forced_pt = lastP;
 191                 lastP += step;
 192                 M += step;
 193             } while (lastP < N && ExtendFit(off + curP, M, data,
 194                                             (contains_forced) ? 0.05 * treshhold : treshhold,
 195                                             res, worstP) );
 196             if (lastP >= N) {
 197
 198                 lastP -= step;
 199                 M -= step;
 200
 201             } else {
 202                 // le dernier a echoue
 203                 lastP -= step;
 204                 M -= step;
 205
 206                 if ( contains_forced ) {
 207                     lastP = forced_pt;
 208                     M = lastP - curP + 1;
 209                 }
 210
 211                 AttemptSimplify(off + curP, M, treshhold, res, worstP);       // ca passe forcement
 212             }
 213             step /= 2;
 214         }
 215
 216         endToPt = pts[off + lastP].p;
 217         if (M <= 2) {
 218             LineTo(endToPt);
 219         } else {
 220             CubicTo(endToPt, res.start, res.end);
 221         }
 222
 223         curP = lastP;
 224     }
 225
 226     if (NR::LInfty(endToPt - moveToPt) < 0.00001) {
 227         Close();
 228     }
 229
 230     g_free(data.Xk);
 231     g_free(data.Yk);
 232     g_free(data.Qk);
 233     g_free(data.tk);
 234     g_free(data.lk);
 235     g_free(data.fk);
 236 }
 237
 238
 239 // warning: slow
 240 // idea behing this feature: splotches appear when trying to fit a small number of points: you can
 241 // get a cubic bezier that fits the points very well but doesn't fit the polyline itself
 242 // so we add a bit of the error at the middle of each segment of the polyline
 243 // also we restrict this to <=20 points, to avoid unnecessary computations
 244 #define with_splotch_killer
 245
 246 // primitive= calc the cubic bezier patche that fits Xk and Yk best
 247 // Qk est deja alloue
 248 // retourne false si probleme (matrice non-inversible)
 249 bool Path::FitCubic(NR::Point const &start, PathDescrCubicTo &res,
 250                     double *Xk, double *Yk, double *Qk, double *tk, int nbPt)
 251 {
 252     NR::Point const end = res.p;
 253
 254     // la matrice tNN
 255     NR::Matrix M(0, 0, 0, 0, 0, 0);
 256     for (int i = 1; i < nbPt - 1; i++) {
 257         M[0] += N13(tk[i]) * N13(tk[i]);
 258         M[1] += N23(tk[i]) * N13(tk[i]);
 259         M[2] += N13(tk[i]) * N23(tk[i]);
 260         M[3] += N23(tk[i]) * N23(tk[i]);
 261     }
 262
 263     double const det = M.det();
 264     if (fabs(det) < 0.000001) {
 265         res.start[0]=res.start[1]=0.0;
 266         res.end[0]=res.end[1]=0.0;
 267         return false;
 268     }
 269
 270     NR::Matrix const iM = M.inverse();
 271     M = iM;
 272
 273     // phase 1: abcisses
 274     // calcul des Qk
 275     Xk[0] = start[0];
 276     Yk[0] = start[1];
 277     Xk[nbPt - 1] = end[0];
 278     Yk[nbPt - 1] = end[1];
 279
 280     for (int i = 1; i < nbPt - 1; i++) {
 281         Qk[i] = Xk[i] - N03 (tk[i]) * Xk[0] - N33 (tk[i]) * Xk[nbPt - 1];
 282     }
 283
 284     // le vecteur Q
 285     NR::Point Q(0, 0);
 286     for (int i = 1; i < nbPt - 1; i++) {
 287         Q[0] += N13 (tk[i]) * Qk[i];
 288         Q[1] += N23 (tk[i]) * Qk[i];
 289     }
 290
 291     NR::Point P = Q * M;
 292     NR::Point cp1;
 293     NR::Point cp2;
 294     cp1[NR::X] = P[NR::X];
 295     cp2[NR::X] = P[NR::Y];
 296
 297     // phase 2: les ordonnees
 298     for (int i = 1; i < nbPt - 1; i++) {
 299         Qk[i] = Yk[i] - N03 (tk[i]) * Yk[0] - N33 (tk[i]) * Yk[nbPt - 1];
 300     }
 301
 302     // le vecteur Q
 303     Q = NR::Point(0, 0);
 304     for (int i = 1; i < nbPt - 1; i++) {
 305         Q[0] += N13 (tk[i]) * Qk[i];
 306         Q[1] += N23 (tk[i]) * Qk[i];
 307     }
 308
 309     P = Q * M;
 310     cp1[NR::Y] = P[NR::X];
 311     cp2[NR::Y] = P[NR::Y];
 312
 313     res.start = 3.0 * (cp1 - start);
 314     res.end = 3.0 * (end - cp2 );
 315
 316     return true;
 317 }
 318
 319
 320 bool Path::ExtendFit(int off, int N, fitting_tables &data, double treshhold, PathDescrCubicTo &res, int &worstP)
 321 {
 322     if ( N >= data.maxPt ) {
 323         data.maxPt = 2 * N + 1;
 324         data.Xk = (double *) g_realloc(data.Xk, data.maxPt * sizeof(double));
 325         data.Yk = (double *) g_realloc(data.Yk, data.maxPt * sizeof(double));
 326         data.Qk = (double *) g_realloc(data.Qk, data.maxPt * sizeof(double));
 327         data.tk = (double *) g_realloc(data.tk, data.maxPt * sizeof(double));
 328         data.lk = (double *) g_realloc(data.lk, data.maxPt * sizeof(double));
 329         data.fk = (char *) g_realloc(data.fk, data.maxPt * sizeof(char));
 330     }
 331
 332     if ( N > data.inPt ) {
 333         for (int i = data.inPt; i < N; i++) {
 334             data.Xk[i] = pts[off + i].p[NR::X];
 335             data.Yk[i] = pts[off + i].p[NR::Y];
 336             data.fk[i] = ( pts[off + i].isMoveTo == polyline_forced ) ? 0x01 : 0x00;
 337         }
 338         data.lk[0] = 0;
 339         data.tk[0] = 0;
 340
 341         double prevLen = 0;
 342         for (int i = 0; i < data.inPt; i++) {
 343             prevLen += data.lk[i];
 344         }
 345         data.totLen = prevLen;
 346
 347         for (int i = ( (data.inPt > 0) ? data.inPt : 1); i < N; i++) {
 348             NR::Point diff;
 349             diff[NR::X] = data.Xk[i] - data.Xk[i - 1];
 350             diff[NR::Y] = data.Yk[i] - data.Yk[i - 1];
 351             data.lk[i] = NR::L2(diff);
 352             data.totLen += data.lk[i];
 353             data.tk[i] = data.totLen;
 354         }
 355
 356         for (int i = 0; i < data.inPt; i++) {
 357             data.tk[i] *= prevLen;
 358             data.tk[i] /= data.totLen;
 359         }
 360
 361         for (int i = data.inPt; i < N; i++) {
 362             data.tk[i] /= data.totLen;
 363         }
 364         data.inPt = N;
 365     }
 366
 367     if ( N < data.nbPt ) {
 368         // on est all\8e trop loin
 369         // faut recalculer les tk
 370         data.totLen = 0;
 371         data.tk[0] = 0;
 372         data.lk[0] = 0;
 373         for (int i = 1; i < N; i++) {
 374             data.totLen += data.lk[i];
 375             data.tk[i] = data.totLen;
 376         }
 377
 378         for (int i = 1; i < N; i++) {
 379             data.tk[i] /= data.totLen;
 380         }
 381     }
 382
 383     data.nbPt = N;
 384
 385     if ( data.nbPt <= 0 ) {
 386         return false;
 387     }
 388
 389     res.p[0] = data.Xk[data.nbPt - 1];
 390     res.p[1] = data.Yk[data.nbPt - 1];
 391     res.start[0] = res.start[1] = 0;
 392     res.end[0] = res.end[1] = 0;
 393     worstP = 1;
 394     if ( N <= 2 ) {
 395         return true;
 396     }
 397
 398     if ( data.totLen < 0.0001 ) {
 399         double worstD = 0;
 400         NR::Point start;
 401         worstP = -1;
 402         start[0] = data.Xk[0];
 403         start[1] = data.Yk[0];
 404         for (int i = 1; i < N; i++) {
 405             NR::Point nPt;
 406             bool isForced = data.fk[i];
 407             nPt[0] = data.Xk[i];
 408             nPt[1] = data.Yk[i];
 409
 410             double nle = DistanceToCubic(start, res, nPt);
 411             if ( isForced ) {
 412                 // forced points are favored for splitting the recursion; we do this by increasing their distance
 413                 if ( worstP < 0 || 2*nle > worstD ) {
 414                     worstP = i;
 415                     worstD = 2*nle;
 416                 }
 417             } else {
 418                 if ( worstP < 0 || nle > worstD ) {
 419                     worstP = i;
 420                     worstD = nle;
 421                 }
 422             }
 423         }
 424
 425         return true;
 426     }
 427
 428     return AttemptSimplify(data, treshhold, res, worstP);
 429 }
 430
 431
 432 // fit a polyline to a bezier patch, return true is treshhold not exceeded (ie: you can continue)
 433 // version that uses tables from the previous iteration, to minimize amount of work done
 434 bool Path::AttemptSimplify (fitting_tables &data,double treshhold, PathDescrCubicTo & res,int &worstP)
 435 {
 436     NR::Point start,end;
 437     // pour une coordonnee
 438     NR::Point cp1, cp2;
 439
 440     worstP = 1;
 441     if (pts.size() == 2) {
 442         return true;
 443     }
 444
 445     start[0] = data.Xk[0];
 446     start[1] = data.Yk[0];
 447     cp1[0] = data.Xk[1];
 448     cp1[1] = data.Yk[1];
 449     end[0] = data.Xk[data.nbPt - 1];
 450     end[1] = data.Yk[data.nbPt - 1];
 451     cp2 = cp1;
 452
 453     if (pts.size()  == 3) {
 454         // start -> cp1 -> end
 455         res.start = cp1 - start;
 456         res.end = end - cp1;
 457         worstP = 1;
 458         return true;
 459     }
 460
 461     if ( FitCubic(start, res, data.Xk, data.Yk, data.Qk, data.tk, data.nbPt) ) {
 462         cp1 = start + res.start / 3;
 463         cp2 = end - res.end / 3;
 464     } else {
 465         // aie, non-inversible
 466         double worstD = 0;
 467         worstP = -1;
 468         for (int i = 1; i < data.nbPt; i++) {
 469             NR::Point nPt;
 470             nPt[NR::X] = data.Xk[i];
 471             nPt[NR::Y] = data.Yk[i];
 472             double nle = DistanceToCubic(start, res, nPt);
 473             if ( data.fk[i] ) {
 474                 // forced points are favored for splitting the recursion; we do this by increasing their distance
 475                 if ( worstP < 0 || 2 * nle > worstD ) {
 476                     worstP = i;
 477                     worstD = 2 * nle;
 478                 }
 479             } else {
 480                 if ( worstP < 0 || nle > worstD ) {
 481                     worstP = i;
 482                     worstD = nle;
 483                 }
 484             }
 485         }
 486         return false;
 487     }
 488
 489     // calcul du delta= pondere par les longueurs des segments
 490     double delta = 0;
 491     {
 492         double worstD = 0;
 493         worstP = -1;
 494         NR::Point prevAppP;
 495         NR::Point prevP;
 496         double prevDist;
 497         prevP[NR::X] = data.Xk[0];
 498         prevP[NR::Y] = data.Yk[0];
 499         prevAppP = prevP; // le premier seulement
 500         prevDist = 0;
 501 #ifdef with_splotch_killer
 502         if ( data.nbPt <= 20 ) {
 503             for (int i = 1; i < data.nbPt - 1; i++) {
 504                 NR::Point curAppP;
 505                 NR::Point curP;
 506                 double curDist;
 507                 NR::Point midAppP;
 508                 NR::Point midP;
 509                 double midDist;
 510
 511                 curAppP[NR::X] = N13(data.tk[i]) * cp1[NR::X] +
 512                     N23(data.tk[i]) * cp2[NR::X] +
 513                     N03(data.tk[i]) * data.Xk[0] +
 514                     N33(data.tk[i]) * data.Xk[data.nbPt - 1];
 515
 516                 curAppP[NR::Y] = N13(data.tk[i]) * cp1[NR::Y] +
 517                     N23(data.tk[i]) * cp2[NR::Y] +
 518                     N03(data.tk[i]) * data.Yk[0] +
 519                     N33(data.tk[i]) * data.Yk[data.nbPt - 1];
 520
 521                 curP[NR::X] = data.Xk[i];
 522                 curP[NR::Y] = data.Yk[i];
 523                 double mtk = 0.5 * (data.tk[i] + data.tk[i - 1]);
 524
 525                 midAppP[NR::X] = N13(mtk) * cp1[NR::X] +
 526                     N23(mtk) * cp2[NR::X] +
 527                     N03(mtk) * data.Xk[0] +
 528                     N33(mtk) * data.Xk[data.nbPt - 1];
 529
 530                 midAppP[NR::Y] = N13(mtk) * cp1[NR::Y] +
 531                     N23(mtk) * cp2[NR::Y] +
 532                     N03(mtk) * data.Yk[0] +
 533                     N33(mtk) * data.Yk[data.nbPt - 1];
 534
 535                 midP = 0.5 * (curP + prevP);
 536
 537                 NR::Point diff = curAppP - curP;
 538                 curDist = dot(diff, diff);
 539                 diff = midAppP - midP;
 540                 midDist = dot(diff, diff);
 541
 542                 delta += 0.3333 * (curDist + prevDist + midDist) * data.lk[i];
 543                 if ( curDist > worstD ) {
 544                     worstD = curDist;
 545                     worstP = i;
 546                 } else if ( data.fk[i] && 2 * curDist > worstD ) {
 547                     worstD = 2*curDist;
 548                     worstP = i;
 549                 }
 550                 prevP = curP;
 551                 prevAppP = curAppP;
 552                 prevDist = curDist;
 553             }
 554             delta /= data.totLen;
 555
 556         } else {
 557 #endif
 558             for (int i = 1; i < data.nbPt - 1; i++) {
 559                 NR::Point curAppP;
 560                 NR::Point curP;
 561                 double    curDist;
 562
 563                 curAppP[NR::X] = N13(data.tk[i]) * cp1[NR::X] +
 564                     N23(data.tk[i]) * cp2[NR::X] +
 565                     N03(data.tk[i]) * data.Xk[0] +
 566                     N33(data.tk[i]) * data.Xk[data.nbPt - 1];
 567
 568                 curAppP[NR::Y] = N13(data.tk[i]) * cp1[NR::Y] +
 569                     N23(data.tk[i]) * cp2[NR::Y] +
 570                     N03(data.tk[i]) * data.Yk[0] +
 571                     N33(data.tk[i]) * data.Yk[data.nbPt - 1];
 572
 573                 curP[NR::X] = data.Xk[i];
 574                 curP[NR::Y] = data.Yk[i];
 575
 576                 NR::Point diff = curAppP-curP;
 577                 curDist = dot(diff, diff);
 578                 delta += curDist;
 579
 580                 if ( curDist > worstD ) {
 581                     worstD = curDist;
 582                     worstP = i;
 583                 } else if ( data.fk[i] && 2 * curDist > worstD ) {
 584                     worstD = 2*curDist;
 585                     worstP = i;
 586                 }
 587                 prevP = curP;
 588                 prevAppP = curAppP;
 589                 prevDist = curDist;
 590             }
 591 #ifdef with_splotch_killer
 592         }
 593 #endif
 594     }
 595
 596     if (delta < treshhold * treshhold) {
 597         // premier jet
 598
 599         // Refine a little.
 600         for (int i = 1; i < data.nbPt - 1; i++) {
 601             NR::Point pt(data.Xk[i], data.Yk[i]);
 602             data.tk[i] = RaffineTk(pt, start, cp1, cp2, end, data.tk[i]);
 603             if (data.tk[i] < data.tk[i - 1]) {
 604                 // Force tk to be monotonic non-decreasing.
 605                 data.tk[i] = data.tk[i - 1];
 606             }
 607         }
 608
 609         if ( FitCubic(start, res, data.Xk, data.Yk, data.Qk, data.tk, data.nbPt) == false) {
 610             // ca devrait jamais arriver, mais bon
 611             res.start = 3.0 * (cp1 - start);
 612             res.end = 3.0 * (end - cp2 );
 613             return true;
 614         }
 615
 616         double ndelta = 0;
 617         {
 618             double worstD = 0;
 619             worstP = -1;
 620             NR::Point prevAppP;
 621             NR::Point prevP(data.Xk[0], data.Yk[0]);
 622             double prevDist = 0;
 623             prevAppP = prevP; // le premier seulement
 624 #ifdef with_splotch_killer
 625             if ( data.nbPt <= 20 ) {
 626                 for (int i = 1; i < data.nbPt - 1; i++) {
 627                     NR::Point curAppP;
 628                     NR::Point curP;
 629                     double  curDist;
 630                     NR::Point midAppP;
 631                     NR::Point midP;
 632                     double  midDist;
 633
 634                     curAppP[NR::X] = N13(data.tk[i]) * cp1[NR::X] +
 635                         N23(data.tk[i]) * cp2[NR::X] +
 636                         N03(data.tk[i]) * data.Xk[0] +
 637                         N33(data.tk[i]) * data.Xk[data.nbPt - 1];
 638
 639                     curAppP[NR::Y] = N13(data.tk[i]) * cp1[NR::Y] +
 640                         N23(data.tk[i]) * cp2[NR::Y] +
 641                         N03(data.tk[i]) * data.Yk[0] +
 642                         N33(data.tk[i]) * data.Yk[data.nbPt - 1];
 643
 644                     curP[NR::X] = data.Xk[i];
 645                     curP[NR::Y] = data.Yk[i];
 646                     double mtk = 0.5 * (data.tk[i] + data.tk[i - 1]);
 647
 648                     midAppP[NR::X] = N13(mtk) * cp1[NR::X] +
 649                         N23(mtk) * cp2[NR::X] +
 650                         N03(mtk) * data.Xk[0] +
 651                         N33(mtk) * data.Xk[data.nbPt - 1];
 652
 653                     midAppP[NR::Y] = N13(mtk) * cp1[NR::Y] +
 654                         N23(mtk) * cp2[NR::Y] +
 655                         N03(mtk) * data.Yk[0] +
 656                         N33(mtk) * data.Yk[data.nbPt - 1];
 657
 658                     midP = 0.5 * (curP + prevP);
 659
 660                     NR::Point diff = curAppP - curP;
 661                     curDist = dot(diff, diff);
 662
 663                     diff = midAppP - midP;
 664                     midDist = dot(diff, diff);
 665
 666                     ndelta += 0.3333 * (curDist + prevDist + midDist) * data.lk[i];
 667
 668                     if ( curDist > worstD ) {
 669                         worstD = curDist;
 670                         worstP = i;
 671                     } else if ( data.fk[i] && 2 * curDist > worstD ) {
 672                         worstD = 2*curDist;
 673                         worstP = i;
 674                     }
 675
 676                     prevP = curP;
 677                     prevAppP = curAppP;
 678                     prevDist = curDist;
 679                 }
 680                 ndelta /= data.totLen;
 681             } else {
 682 #endif
 683                 for (int i = 1; i < data.nbPt - 1; i++) {
 684                     NR::Point curAppP;
 685                     NR::Point curP;
 686                     double    curDist;
 687
 688                     curAppP[NR::X] = N13(data.tk[i]) * cp1[NR::X] +
 689                         N23(data.tk[i]) * cp2[NR::X] +
 690                         N03(data.tk[i]) * data.Xk[0] +
 691                         N33(data.tk[i]) * data.Xk[data.nbPt - 1];
 692
 693                     curAppP[NR::Y] = N13(data.tk[i]) * cp1[NR::Y] +
 694                         N23(data.tk[i]) * cp2[1] +
 695                         N03(data.tk[i]) * data.Yk[0] +
 696                         N33(data.tk[i]) * data.Yk[data.nbPt - 1];
 697
 698                     curP[NR::X] = data.Xk[i];
 699                     curP[NR::Y] = data.Yk[i];
 700
 701                     NR::Point diff = curAppP - curP;
 702                     curDist = dot(diff, diff);
 703
 704                     ndelta += curDist;
 705
 706                     if ( curDist > worstD ) {
 707                         worstD = curDist;
 708                         worstP = i;
 709                     } else if ( data.fk[i] && 2 * curDist > worstD ) {
 710                         worstD = 2 * curDist;
 711                         worstP = i;
 712                     }
 713                     prevP = curP;
 714                     prevAppP = curAppP;
 715                     prevDist = curDist;
 716                 }
 717 #ifdef with_splotch_killer
 718             }
 719 #endif
 720         }
 721
 722         if (ndelta < delta + 0.00001) {
 723             return true;
 724         } else {
 725             // nothing better to do
 726             res.start = 3.0 * (cp1 - start);
 727             res.end = 3.0 * (end - cp2 );
 728         }
 729
 730         return true;
 731     }
 732
 733   return false;
 734 }
 735
 736
 737 bool Path::AttemptSimplify(int off, int N, double treshhold, PathDescrCubicTo &res,int &worstP)
 738 {
 739     NR::Point start;
 740     NR::Point end;
 741
 742     // pour une coordonnee
 743     double *Xk;                         // la coordonnee traitee (x puis y)
 744     double *Yk;                         // la coordonnee traitee (x puis y)
 745     double *lk;                         // les longueurs de chaque segment
 746     double *tk;                         // les tk
 747     double *Qk;                         // les Qk
 748     char *fk;       // si point force
 749
 750     NR::Point cp1;
 751     NR::Point cp2;
 752
 753     if (N == 2) {
 754         worstP = 1;
 755         return true;
 756     }
 757
 758     start = pts[off].p;
 759     cp1 = pts[off + 1].p;
 760     end = pts[off + N - 1].p;
 761
 762     res.p = end;
 763     res.start[0] = res.start[1] = 0;
 764     res.end[0] = res.end[1] = 0;
 765     if (N == 3) {
 766         // start -> cp1 -> end
 767         res.start = cp1 - start;
 768         res.end = end - cp1;
 769         worstP = 1;
 770         return true;
 771     }
 772
 773     // Totally inefficient, allocates & deallocates all the time.
 774     tk = (double *) g_malloc(N * sizeof(double));
 775     Qk = (double *) g_malloc(N * sizeof(double));
 776     Xk = (double *) g_malloc(N * sizeof(double));
 777     Yk = (double *) g_malloc(N * sizeof(double));
 778     lk = (double *) g_malloc(N * sizeof(double));
 779     fk = (char *) g_malloc(N * sizeof(char));
 780
 781     // chord length method
 782     tk[0] = 0.0;
 783     lk[0] = 0.0;
 784     {
 785         NR::Point prevP = start;
 786         for (int i = 1; i < N; i++) {
 787             Xk[i] = pts[off + i].p[NR::X];
 788             Yk[i] = pts[off + i].p[NR::Y];
 789
 790             if ( pts[off + i].isMoveTo == polyline_forced ) {
 791                 fk[i] = 0x01;
 792             } else {
 793                 fk[i] = 0;
 794             }
 795
 796             NR::Point diff(Xk[i] - prevP[NR::X], Yk[i] - prevP[1]);
 797             prevP[0] = Xk[i];
 798             prevP[1] = Yk[i];
 799             lk[i] = NR::L2(diff);
 800             tk[i] = tk[i - 1] + lk[i];
 801         }
 802     }
 803
 804     if (tk[N - 1] < 0.00001) {
 805         // longueur nulle
 806         res.start[0] = res.start[1] = 0;
 807         res.end[0] = res.end[1] = 0;
 808         double worstD = 0;
 809         worstP = -1;
 810         for (int i = 1; i < N; i++) {
 811             NR::Point nPt;
 812             bool isForced = fk[i];
 813             nPt[0] = Xk[i];
 814             nPt[1] = Yk[i];
 815
 816             double nle = DistanceToCubic(start, res, nPt);
 817             if ( isForced ) {
 818                 // forced points are favored for splitting the recursion; we do this by increasing their distance
 819                 if ( worstP < 0 || 2 * nle > worstD ) {
 820                     worstP = i;
 821                     worstD = 2 * nle;
 822                 }
 823             } else {
 824                 if ( worstP < 0 || nle > worstD ) {
 825                     worstP = i;
 826                     worstD = nle;
 827                 }
 828             }
 829         }
 830
 831         g_free(tk);
 832         g_free(Qk);
 833         g_free(Xk);
 834         g_free(Yk);
 835         g_free(fk);
 836         g_free(lk);
 837
 838         return false;
 839     }
 840
 841     double totLen = tk[N - 1];
 842     for (int i = 1; i < N - 1; i++) {
 843         tk[i] /= totLen;
 844     }
 845
 846     res.p = end;
 847     if ( FitCubic(start, res, Xk, Yk, Qk, tk, N) ) {
 848         cp1 = start + res.start / 3;
 849         cp2 = end + res.end / 3;
 850     } else {
 851         // aie, non-inversible
 852         res.start[0] = res.start[1] = 0;
 853         res.end[0] = res.end[1] = 0;
 854         double worstD = 0;
 855         worstP = -1;
 856         for (int i = 1; i < N; i++) {
 857             NR::Point nPt(Xk[i], Yk[i]);
 858             bool isForced = fk[i];
 859             double nle = DistanceToCubic(start, res, nPt);
 860             if ( isForced ) {
 861                 // forced points are favored for splitting the recursion; we do this by increasing their distance
 862                 if ( worstP < 0 || 2 * nle > worstD ) {
 863                     worstP = i;
 864                     worstD = 2 * nle;
 865                 }
 866             } else {
 867                 if ( worstP < 0 || nle > worstD ) {
 868                     worstP = i;
 869                     worstD = nle;
 870                 }
 871             }
 872         }
 873
 874         g_free(tk);
 875         g_free(Qk);
 876         g_free(Xk);
 877         g_free(Yk);
 878         g_free(fk);
 879         g_free(lk);
 880         return false;
 881     }
 882
 883     // calcul du delta= pondere par les longueurs des segments
 884     double delta = 0;
 885     {
 886         double worstD = 0;
 887         worstP = -1;
 888         NR::Point prevAppP;
 889     NR::Point   prevP;
 890     double      prevDist;
 891     prevP[0] = Xk[0];
 892     prevP[1] = Yk[0];
 893     prevAppP = prevP; // le premier seulement
 894     prevDist = 0;
 895 #ifdef with_splotch_killer
 896     if ( N <= 20 ) {
 897       for (int i = 1; i < N - 1; i++)
 898       {
 899         NR::Point curAppP;
 900         NR::Point curP;
 901         double    curDist;
 902         NR::Point midAppP;
 903         NR::Point midP;
 904         double    midDist;
 905
 906         curAppP[0] = N13 (tk[i]) * cp1[0] + N23 (tk[i]) * cp2[0] + N03 (tk[i]) * Xk[0] + N33 (tk[i]) * Xk[N - 1];
 907         curAppP[1] = N13 (tk[i]) * cp1[1] + N23 (tk[i]) * cp2[1] + N03 (tk[i]) * Yk[0] + N33 (tk[i]) * Yk[N - 1];
 908         curP[0] = Xk[i];
 909         curP[1] = Yk[i];
 910         midAppP[0] = N13 (0.5*(tk[i]+tk[i-1])) * cp1[0] + N23 (0.5*(tk[i]+tk[i-1])) * cp2[0] + N03 (0.5*(tk[i]+tk[i-1])) * Xk[0] + N33 (0.5*(tk[i]+tk[i-1])) * Xk[N - 1];
 911         midAppP[1] = N13 (0.5*(tk[i]+tk[i-1])) * cp1[1] + N23 (0.5*(tk[i]+tk[i-1])) * cp2[1] + N03 (0.5*(tk[i]+tk[i-1])) * Yk[0] + N33 (0.5*(tk[i]+tk[i-1])) * Yk[N - 1];
 912         midP=0.5*(curP+prevP);
 913
 914         NR::Point diff;
 915         diff = curAppP-curP;
 916         curDist = dot(diff,diff);
 917
 918         diff = midAppP-midP;
 919         midDist = dot(diff,diff);
 920
 921         delta+=0.3333*(curDist+prevDist+midDist)/**lk[i]*/;
 922
 923         if ( curDist > worstD ) {
 924           worstD = curDist;
 925           worstP = i;
 926         } else if ( fk[i] && 2*curDist > worstD ) {
 927           worstD = 2*curDist;
 928           worstP = i;
 929         }
 930         prevP = curP;
 931         prevAppP = curAppP;
 932         prevDist = curDist;
 933       }
 934       delta/=totLen;
 935     } else {
 936 #endif
 937       for (int i = 1; i < N - 1; i++)
 938       {
 939         NR::Point curAppP;
 940         NR::Point curP;
 941         double    curDist;
 942
 943         curAppP[0] = N13 (tk[i]) * cp1[0] + N23 (tk[i]) * cp2[0] + N03 (tk[i]) * Xk[0] + N33 (tk[i]) * Xk[N - 1];
 944         curAppP[1] = N13 (tk[i]) * cp1[1] + N23 (tk[i]) * cp2[1] + N03 (tk[i]) * Yk[0] + N33 (tk[i]) * Yk[N - 1];
 945         curP[0] = Xk[i];
 946         curP[1] = Yk[i];
 947
 948         NR::Point diff;
 949         diff = curAppP-curP;
 950         curDist = dot(diff,diff);
 951         delta += curDist;
 952         if ( curDist > worstD ) {
 953           worstD = curDist;
 954           worstP = i;
 955         } else if ( fk[i] && 2*curDist > worstD ) {
 956           worstD = 2*curDist;
 957           worstP = i;
 958         }
 959         prevP = curP;
 960         prevAppP = curAppP;
 961         prevDist = curDist;
 962       }
 963 #ifdef with_splotch_killer
 964     }
 965 #endif
 966   }
 967
 968   if (delta < treshhold * treshhold)
 969   {
 970     // premier jet
 971     res.start = 3.0 * (cp1 - start);
 972     res.end = -3.0 * (cp2 - end);
 973     res.p = end;
 974
 975     // Refine a little.
 976     for (int i = 1; i < N - 1; i++)
 977     {
 978       NR::Point
 979             pt;
 980       pt[0] = Xk[i];
 981       pt[1] = Yk[i];
 982       tk[i] = RaffineTk (pt, start, cp1, cp2, end, tk[i]);
 983       if (tk[i] < tk[i - 1])
 984             {
 985               // Force tk to be monotonic non-decreasing.
 986               tk[i] = tk[i - 1];
 987             }
 988     }
 989
 990     if ( FitCubic(start,res,Xk,Yk,Qk,tk,N) ) {
 991     } else {
 992       // ca devrait jamais arriver, mais bon
 993       res.start = 3.0 * (cp1 - start);
 994       res.end = -3.0 * (cp2 - end);
 995       g_free(tk);
 996       g_free(Qk);
 997       g_free(Xk);
 998       g_free(Yk);
 999       g_free(fk);
1000       g_free(lk);
1001       return true;
1002     }
1003     double ndelta = 0;
1004     {
1005       double worstD = 0;
1006       worstP = -1;
1007       NR::Point   prevAppP;
1008       NR::Point   prevP;
1009       double      prevDist;
1010       prevP[0] = Xk[0];
1011       prevP[1] = Yk[0];
1012       prevAppP = prevP; // le premier seulement
1013       prevDist = 0;
1014 #ifdef with_splotch_killer
1015       if ( N <= 20 ) {
1016         for (int i = 1; i < N - 1; i++)
1017         {
1018           NR::Point curAppP;
1019           NR::Point curP;
1020           double    curDist;
1021           NR::Point midAppP;
1022           NR::Point midP;
1023           double    midDist;
1024
1025           curAppP[0] = N13 (tk[i]) * cp1[0] + N23 (tk[i]) * cp2[0] + N03 (tk[i]) * Xk[0] + N33 (tk[i]) * Xk[N - 1];
1026           curAppP[1] = N13 (tk[i]) * cp1[1] + N23 (tk[i]) * cp2[1] + N03 (tk[i]) * Yk[0] + N33 (tk[i]) * Yk[N - 1];
1027           curP[0] = Xk[i];
1028           curP[1] = Yk[i];
1029           midAppP[0] = N13 (0.5*(tk[i]+tk[i-1])) * cp1[0] + N23 (0.5*(tk[i]+tk[i-1])) * cp2[0] + N03 (0.5*(tk[i]+tk[i-1])) * Xk[0] + N33 (0.5*(tk[i]+tk[i-1])) * Xk[N - 1];
1030           midAppP[1] = N13 (0.5*(tk[i]+tk[i-1])) * cp1[1] + N23 (0.5*(tk[i]+tk[i-1])) * cp2[1] + N03 (0.5*(tk[i]+tk[i-1])) * Yk[0] + N33 (0.5*(tk[i]+tk[i-1])) * Yk[N - 1];
1031           midP = 0.5*(curP+prevP);
1032
1033           NR::Point diff;
1034           diff = curAppP-curP;
1035           curDist = dot(diff,diff);
1036           diff = midAppP-midP;
1037           midDist = dot(diff,diff);
1038
1039           ndelta+=0.3333*(curDist+prevDist+midDist)/**lk[i]*/;
1040
1041           if ( curDist > worstD ) {
1042             worstD = curDist;
1043             worstP = i;
1044           } else if ( fk[i] && 2*curDist > worstD ) {
1045             worstD = 2*curDist;
1046             worstP = i;
1047           }
1048           prevP = curP;
1049           prevAppP = curAppP;
1050           prevDist = curDist;
1051         }
1052         ndelta /= totLen;
1053       } else {
1054 #endif
1055         for (int i = 1; i < N - 1; i++)
1056         {
1057           NR::Point curAppP;
1058           NR::Point curP;
1059           double    curDist;
1060
1061           curAppP[0] = N13 (tk[i]) * cp1[0] + N23 (tk[i]) * cp2[0] + N03 (tk[i]) * Xk[0] + N33 (tk[i]) * Xk[N - 1];
1062           curAppP[1] = N13 (tk[i]) * cp1[1] + N23 (tk[i]) * cp2[1] + N03 (tk[i]) * Yk[0] + N33 (tk[i]) * Yk[N - 1];
1063           curP[0]=Xk[i];
1064           curP[1]=Yk[i];
1065
1066           NR::Point diff;
1067           diff=curAppP-curP;
1068           curDist=dot(diff,diff);
1069           ndelta+=curDist;
1070
1071           if ( curDist > worstD ) {
1072             worstD=curDist;
1073             worstP=i;
1074           } else if ( fk[i] && 2*curDist > worstD ) {
1075             worstD=2*curDist;
1076             worstP=i;
1077           }
1078           prevP=curP;
1079           prevAppP=curAppP;
1080           prevDist=curDist;
1081         }
1082 #ifdef with_splotch_killer
1083       }
1084 #endif
1085     }
1086
1087     g_free(tk);
1088     g_free(Qk);
1089     g_free(Xk);
1090     g_free(Yk);
1091     g_free(fk);
1092     g_free(lk);
1093
1094     if (ndelta < delta + 0.00001)
1095     {
1096       return true;
1097     } else {
1098       // nothing better to do
1099       res.start = 3.0 * (cp1 - start);
1100       res.end = -3.0 * (cp2 - end);
1101     }
1102     return true;
1103   } else {
1104     // nothing better to do
1105   }
1106
1107   g_free(tk);
1108   g_free(Qk);
1109   g_free(Xk);
1110   g_free(Yk);
1111   g_free(fk);
1112   g_free(lk);
1113   return false;
1114 }
1115
1116 double Path::RaffineTk (NR::Point pt, NR::Point p0, NR::Point p1, NR::Point p2, NR::Point p3, double it)
1117 {
1118     // Refinement of the tk values.
1119     // Just one iteration of Newtow Raphson, given that we're approaching the curve anyway.
1120     // [fr: vu que de toute facon la courbe est approchC)e]
1121     double const Ax = pt[NR::X] -
1122         p0[NR::X] * N03(it) -
1123         p1[NR::X] * N13(it) -
1124         p2[NR::X] * N23(it) -
1125         p3[NR::X] * N33(it);
1126
1127     double const Bx = (p1[NR::X] - p0[NR::X]) * N02(it) +
1128         (p2[NR::X] - p1[NR::X]) * N12(it) +
1129         (p3[NR::X] - p2[NR::X]) * N22(it);
1130
1131     double const Cx = (p0[NR::X] - 2 * p1[NR::X] + p2[NR::X]) * N01(it) +
1132         (p3[NR::X] - 2 * p2[NR::X] + p1[NR::X]) * N11(it);
1133
1134     double const Ay =  pt[NR::Y] -
1135         p0[NR::Y] * N03(it) -
1136         p1[NR::Y] * N13(it) -
1137         p2[NR::Y] * N23(it) -
1138         p3[NR::Y] * N33(it);
1139
1140     double const By = (p1[NR::Y] - p0[NR::Y]) * N02(it) +
1141         (p2[NR::Y] - p1[NR::Y]) * N12(it) +
1142         (p3[NR::Y] - p2[NR::Y]) * N22(it);
1143
1144     double const Cy = (p0[NR::Y] - 2 * p1[NR::Y] + p2[NR::Y]) * N01(it) +
1145         (p3[NR::Y] - 2 * p2[NR::Y] + p1[NR::Y]) * N11(it);
1146
1147     double const dF = -6 * (Ax * Bx + Ay * By);
1148     double const ddF = 18 * (Bx * Bx + By * By) - 12 * (Ax * Cx + Ay * Cy);
1149     if (fabs (ddF) > 0.0000001) {
1150         return it - dF / ddF;
1151     }
1152
1153     return it;
1154 }
1155
1156 // variation on the fitting theme: try to merge path commands into cubic bezier patches
1157 // the goal was to reduce the number of path commands, especially when ooperations on path produce
1158 // lots of small path elements; ideally you could get rid of very small segments at reduced visual cost
1159 void Path::Coalesce(double tresh)
1160 {
1161     if ( descr_flags & descr_adding_bezier ) {
1162         CancelBezier();
1163     }
1164
1165     if ( descr_flags & descr_doing_subpath ) {
1166         CloseSubpath();
1167     }
1168
1169     if (descr_cmd.size() <= 2) {
1170         return;
1171     }
1172
1173     SetBackData(false);
1174     Path* tempDest = new Path();
1175     tempDest->SetBackData(false);
1176
1177     ConvertEvenLines(0.25*tresh);
1178
1179     int lastP = 0;
1180     int lastAP = -1;
1181     // As the elements are stored in a separate tableau, it's no longer worth optimizing
1182     // the rewriting in the same tableau.
1183     // [[comme les elements sont stockes dans un tableau a part, plus la peine d'optimiser
1184     // la r\8e\8ecriture dans la meme tableau]]
1185
1186     int lastA = descr_cmd[0]->associated;
1187     int prevA = lastA;
1188     NR::Point firstP;
1189
1190     /* FIXME: the use of this variable probably causes a leak or two.
1191     ** It's a hack anyway, and probably only needs to be a type rather than
1192     ** a full PathDescr.
1193     */
1194     PathDescr *lastAddition = new PathDescrMoveTo(NR::Point(0, 0));
1195     bool containsForced = false;
1196     PathDescrCubicTo pending_cubic(NR::Point(0, 0), NR::Point(0, 0), NR::Point(0, 0));
1197
1198     for (int curP = 0; curP < int(descr_cmd.size()); curP++) {
1199         int typ = descr_cmd[curP]->getType();
1200         int nextA = lastA;
1201
1202         if (typ == descr_moveto) {
1203
1204             if (lastAddition->flags != descr_moveto) {
1205                 FlushPendingAddition(tempDest,lastAddition,pending_cubic,lastAP);
1206             }
1207             lastAddition = descr_cmd[curP];
1208             lastAP = curP;
1209             FlushPendingAddition(tempDest, lastAddition, pending_cubic, lastAP);
1210             // Added automatically (too bad about multiple moveto's).
1211             // [fr: (tant pis pour les moveto multiples)]
1212             containsForced = false;
1213
1214             PathDescrMoveTo *nData = dynamic_cast<PathDescrMoveTo *>(descr_cmd[curP]);
1215             firstP = nData->p;
1216             lastA = descr_cmd[curP]->associated;
1217             prevA = lastA;
1218             lastP = curP;
1219
1220         } else if (typ == descr_close) {
1221             nextA = descr_cmd[curP]->associated;
1222             if (lastAddition->flags != descr_moveto) {
1223
1224                 PathDescrCubicTo res(NR::Point(0, 0), NR::Point(0, 0), NR::Point(0, 0));
1225                 int worstP = -1;
1226                 if (AttemptSimplify(lastA, nextA - lastA + 1, (containsForced) ? 0.05 * tresh : tresh, res, worstP)) {
1227                     lastAddition = new PathDescrCubicTo(NR::Point(0, 0),
1228                                                           NR::Point(0, 0),
1229                                                           NR::Point(0, 0));
1230                     pending_cubic = res;
1231                     lastAP = -1;
1232                 }
1233
1234                 FlushPendingAddition(tempDest, lastAddition, pending_cubic, lastAP);
1235                 FlushPendingAddition(tempDest, descr_cmd[curP], pending_cubic, curP);
1236
1237             } else {
1238                 FlushPendingAddition(tempDest,descr_cmd[curP],pending_cubic,curP);
1239             }
1240
1241             containsForced = false;
1242             lastAddition = new PathDescrMoveTo(NR::Point(0, 0));
1243             prevA = lastA = nextA;
1244             lastP = curP;
1245             lastAP = curP;
1246
1247         } else if (typ == descr_forced) {
1248
1249             nextA = descr_cmd[curP]->associated;
1250             if (lastAddition->flags != descr_moveto) {
1251
1252                 PathDescrCubicTo res(NR::Point(0, 0), NR::Point(0, 0), NR::Point(0, 0));
1253                 int worstP = -1;
1254                 if (AttemptSimplify(lastA, nextA - lastA + 1, 0.05 * tresh, res, worstP)) {
1255                     // plus sensible parce que point force
1256                     // ca passe
1257                     containsForced = true;
1258                 } else  {
1259                     // on force l'addition
1260                     FlushPendingAddition(tempDest, lastAddition, pending_cubic, lastAP);
1261                     lastAddition = new PathDescrMoveTo(NR::Point(0, 0));
1262                     prevA = lastA = nextA;
1263                     lastP = curP;
1264                     lastAP = curP;
1265                     containsForced = false;
1266                 }
1267             }
1268
1269         } else if (typ == descr_lineto || typ == descr_cubicto || typ == descr_arcto) {
1270
1271             nextA = descr_cmd[curP]->associated;
1272             if (lastAddition->flags != descr_moveto) {
1273
1274                 PathDescrCubicTo res(NR::Point(0, 0), NR::Point(0, 0), NR::Point(0, 0));
1275                 int worstP = -1;
1276                 if (AttemptSimplify(lastA, nextA - lastA + 1, tresh, res, worstP)) {
1277                     lastAddition = new PathDescrCubicTo(NR::Point(0, 0),
1278                                                           NR::Point(0, 0),
1279                                                           NR::Point(0, 0));
1280                     pending_cubic = res;
1281                     lastAddition->associated = lastA;
1282                     lastP = curP;
1283                     lastAP = -1;
1284                 }  else {
1285                     lastA = descr_cmd[lastP]->associated;       // pourrait etre surecrit par la ligne suivante
1286                     FlushPendingAddition(tempDest, lastAddition, pending_cubic, lastAP);
1287                     lastAddition = descr_cmd[curP];
1288                     if ( typ == descr_cubicto ) {
1289                         pending_cubic = *(dynamic_cast<PathDescrCubicTo*>(descr_cmd[curP]));
1290                     }
1291                     lastAP = curP;
1292                     containsForced = false;
1293                 }
1294
1295             } else {
1296                 lastA = prevA /*descr_cmd[curP-1]->associated */ ;
1297                 lastAddition = descr_cmd[curP];
1298                 if ( typ == descr_cubicto ) {
1299                     pending_cubic = *(dynamic_cast<PathDescrCubicTo*>(descr_cmd[curP]));
1300                 }
1301                 lastAP = curP;
1302                 containsForced = false;
1303             }
1304             prevA = nextA;
1305
1306         } else if (typ == descr_bezierto) {
1307
1308             if (lastAddition->flags != descr_moveto) {
1309                 FlushPendingAddition(tempDest, lastAddition, pending_cubic, lastAP);
1310                 lastAddition = new PathDescrMoveTo(NR::Point(0, 0));
1311             }
1312             lastAP = -1;
1313             lastA = descr_cmd[curP]->associated;
1314             lastP = curP;
1315             PathDescrBezierTo *nBData = dynamic_cast<PathDescrBezierTo*>(descr_cmd[curP]);
1316             for (int i = 1; i <= nBData->nb; i++) {
1317                 FlushPendingAddition(tempDest, descr_cmd[curP + i], pending_cubic, curP + i);
1318             }
1319             curP += nBData->nb;
1320             prevA = nextA;
1321
1322         } else if (typ == descr_interm_bezier) {
1323             continue;
1324         } else {
1325             continue;
1326         }
1327     }
1328
1329     if (lastAddition->flags != descr_moveto) {
1330         FlushPendingAddition(tempDest, lastAddition, pending_cubic, lastAP);
1331     }
1332
1333     Copy(tempDest);
1334     delete tempDest;
1335 }
1336
1337
1338 void Path::FlushPendingAddition(Path *dest, PathDescr *lastAddition,
1339                                 PathDescrCubicTo &lastCubic, int lastAP)
1340 {
1341     switch (lastAddition->getType()) {
1342
1343     case descr_moveto:
1344         if ( lastAP >= 0 ) {
1345             PathDescrMoveTo* nData = dynamic_cast<PathDescrMoveTo *>(descr_cmd[lastAP]);
1346             dest->MoveTo(nData->p);
1347         }
1348         break;
1349
1350     case descr_close:
1351         dest->Close();
1352         break;
1353
1354     case descr_cubicto:
1355         dest->CubicTo(lastCubic.p, lastCubic.start, lastCubic.end);
1356         break;
1357
1358     case descr_lineto:
1359         if ( lastAP >= 0 ) {
1360             PathDescrLineTo *nData = dynamic_cast<PathDescrLineTo *>(descr_cmd[lastAP]);
1361             dest->LineTo(nData->p);
1362         }
1363         break;
1364
1365     case descr_arcto:
1366         if ( lastAP >= 0 ) {
1367             PathDescrArcTo *nData = dynamic_cast<PathDescrArcTo *>(descr_cmd[lastAP]);
1368             dest->ArcTo(nData->p, nData->rx, nData->ry, nData->angle, nData->large, nData->clockwise);
1369         }
1370         break;
1371
1372     case descr_bezierto:
1373         if ( lastAP >= 0 ) {
1374             PathDescrBezierTo *nData = dynamic_cast<PathDescrBezierTo *>(descr_cmd[lastAP]);
1375             dest->BezierTo(nData->p);
1376         }
1377         break;
1378
1379     case descr_interm_bezier:
1380         if ( lastAP >= 0 ) {
1381             PathDescrIntermBezierTo *nData = dynamic_cast<PathDescrIntermBezierTo*>(descr_cmd[lastAP]);
1382             dest->IntermBezierTo(nData->p);
1383         }
1384         break;
1385     }
1386 }
1387
1388 /*
1389   Local Variables:
1390   mode:c++
1391   c-file-style:"stroustrup"
1392   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
1393   indent-tabs-mode:nil
1394   fill-column:99
1395   End:
1396 */
1397 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :