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(Geom::Point const &iS, Geom::Point const &isD,
  74                           Geom::Point const &iE, Geom::Point const &ieD,
  75                           Geom::Point &pt, double current, int lev, double st, double et)
  76 {
  77     if ( lev <= 0 ) {
  78         return current;
  79     }
  80
  81     Geom::Point const m = 0.5 * (iS + iE) + 0.125 * (isD - ieD);
  82     Geom::Point const md = 0.75 * (iE - iS) - 0.125 * (isD + ieD);
  83     double const mt = (st + et) / 2;
  84
  85     Geom::Point const hisD = 0.5 * isD;
  86     Geom::Point const hieD = 0.5 * ieD;
  87
  88     Geom::Point const mp = pt - m;
  89     double nle = Geom::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(Geom::Point const &start, PathDescrCubicTo res, Geom::Point &pt)
 120 {
 121     Geom::Point const sp = pt - start;
 122     Geom::Point const ep = pt - res.p;
 123     double nle = Geom::dot(sp, sp);
 124     double nnle = Geom::dot(ep, ep);
 125     if ( nnle < nle ) {
 126         nle = nnle;
 127     }
 128
 129     Geom::Point seg = res.p - start;
 130     nnle = NR::cross(seg, sp);
 131     nnle *= nnle;
 132     nnle /= Geom::dot(seg, seg);
 133     if ( nnle < nle ) {
 134         if ( Geom::dot(sp,seg) >= 0 ) {
 135             seg = start - res.p;
 136             if ( Geom::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     Geom::Point const moveToPt = pts[off].p;
 167     MoveTo(moveToPt);
 168     Geom::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(Geom::Point(0, 0), Geom::Point(0, 0), Geom::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 (Geom::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(Geom::Point const &start, PathDescrCubicTo &res,
 250                     double *Xk, double *Yk, double *Qk, double *tk, int nbPt)
 251 {
 252     Geom::Point const end = res.p;
 253
 254     // la matrice tNN
 255     Geom::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     Geom::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     Geom::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     Geom::Point P = Q * M;
 292     Geom::Point cp1;
 293     Geom::Point cp2;
 294     cp1[Geom::X] = P[Geom::X];
 295     cp2[Geom::X] = P[Geom::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 = Geom::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[Geom::Y] = P[Geom::X];
 311     cp2[Geom::Y] = P[Geom::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[Geom::X];
 335             data.Yk[i] = pts[off + i].p[Geom::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             Geom::Point diff;
 349             diff[Geom::X] = data.Xk[i] - data.Xk[i - 1];
 350             diff[Geom::Y] = data.Yk[i] - data.Yk[i - 1];
 351             data.lk[i] = Geom::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         // We've gone too far; we'll have to recalulate the .tk.
 369         data.totLen = 0;
 370         data.tk[0] = 0;
 371         data.lk[0] = 0;
 372         for (int i = 1; i < N; i++) {
 373             data.totLen += data.lk[i];
 374             data.tk[i] = data.totLen;
 375         }
 376
 377         for (int i = 1; i < N; i++) {
 378             data.tk[i] /= data.totLen;
 379         }
 380     }
 381
 382     data.nbPt = N;
 383
 384     if ( data.nbPt <= 0 ) {
 385         return false;
 386     }
 387
 388     res.p[0] = data.Xk[data.nbPt - 1];
 389     res.p[1] = data.Yk[data.nbPt - 1];
 390     res.start[0] = res.start[1] = 0;
 391     res.end[0] = res.end[1] = 0;
 392     worstP = 1;
 393     if ( N <= 2 ) {
 394         return true;
 395     }
 396
 397     if ( data.totLen < 0.0001 ) {
 398         double worstD = 0;
 399         Geom::Point start;
 400         worstP = -1;
 401         start[0] = data.Xk[0];
 402         start[1] = data.Yk[0];
 403         for (int i = 1; i < N; i++) {
 404             Geom::Point nPt;
 405             bool isForced = data.fk[i];
 406             nPt[0] = data.Xk[i];
 407             nPt[1] = data.Yk[i];
 408
 409             double nle = DistanceToCubic(start, res, nPt);
 410             if ( isForced ) {
 411                 // forced points are favored for splitting the recursion; we do this by increasing their distance
 412                 if ( worstP < 0 || 2*nle > worstD ) {
 413                     worstP = i;
 414                     worstD = 2*nle;
 415                 }
 416             } else {
 417                 if ( worstP < 0 || nle > worstD ) {
 418                     worstP = i;
 419                     worstD = nle;
 420                 }
 421             }
 422         }
 423
 424         return true;
 425     }
 426
 427     return AttemptSimplify(data, treshhold, res, worstP);
 428 }
 429
 430
 431 // fit a polyline to a bezier patch, return true is treshhold not exceeded (ie: you can continue)
 432 // version that uses tables from the previous iteration, to minimize amount of work done
 433 bool Path::AttemptSimplify (fitting_tables &data,double treshhold, PathDescrCubicTo & res,int &worstP)
 434 {
 435     Geom::Point start,end;
 436     // pour une coordonnee
 437     Geom::Point cp1, cp2;
 438
 439     worstP = 1;
 440     if (pts.size() == 2) {
 441         return true;
 442     }
 443
 444     start[0] = data.Xk[0];
 445     start[1] = data.Yk[0];
 446     cp1[0] = data.Xk[1];
 447     cp1[1] = data.Yk[1];
 448     end[0] = data.Xk[data.nbPt - 1];
 449     end[1] = data.Yk[data.nbPt - 1];
 450     cp2 = cp1;
 451
 452     if (pts.size()  == 3) {
 453         // start -> cp1 -> end
 454         res.start = cp1 - start;
 455         res.end = end - cp1;
 456         worstP = 1;
 457         return true;
 458     }
 459
 460     if ( FitCubic(start, res, data.Xk, data.Yk, data.Qk, data.tk, data.nbPt) ) {
 461         cp1 = start + res.start / 3;
 462         cp2 = end - res.end / 3;
 463     } else {
 464         // aie, non-inversible
 465         double worstD = 0;
 466         worstP = -1;
 467         for (int i = 1; i < data.nbPt; i++) {
 468             Geom::Point nPt;
 469             nPt[Geom::X] = data.Xk[i];
 470             nPt[Geom::Y] = data.Yk[i];
 471             double nle = DistanceToCubic(start, res, nPt);
 472             if ( data.fk[i] ) {
 473                 // forced points are favored for splitting the recursion; we do this by increasing their distance
 474                 if ( worstP < 0 || 2 * nle > worstD ) {
 475                     worstP = i;
 476                     worstD = 2 * nle;
 477                 }
 478             } else {
 479                 if ( worstP < 0 || nle > worstD ) {
 480                     worstP = i;
 481                     worstD = nle;
 482                 }
 483             }
 484         }
 485         return false;
 486     }
 487
 488     // calcul du delta= pondere par les longueurs des segments
 489     double delta = 0;
 490     {
 491         double worstD = 0;
 492         worstP = -1;
 493         Geom::Point prevAppP;
 494         Geom::Point prevP;
 495         double prevDist;
 496         prevP[Geom::X] = data.Xk[0];
 497         prevP[Geom::Y] = data.Yk[0];
 498         prevAppP = prevP; // le premier seulement
 499         prevDist = 0;
 500 #ifdef with_splotch_killer
 501         if ( data.nbPt <= 20 ) {
 502             for (int i = 1; i < data.nbPt - 1; i++) {
 503                 Geom::Point curAppP;
 504                 Geom::Point curP;
 505                 double curDist;
 506                 Geom::Point midAppP;
 507                 Geom::Point midP;
 508                 double midDist;
 509
 510                 curAppP[Geom::X] = N13(data.tk[i]) * cp1[Geom::X] +
 511                     N23(data.tk[i]) * cp2[Geom::X] +
 512                     N03(data.tk[i]) * data.Xk[0] +
 513                     N33(data.tk[i]) * data.Xk[data.nbPt - 1];
 514
 515                 curAppP[Geom::Y] = N13(data.tk[i]) * cp1[Geom::Y] +
 516                     N23(data.tk[i]) * cp2[Geom::Y] +
 517                     N03(data.tk[i]) * data.Yk[0] +
 518                     N33(data.tk[i]) * data.Yk[data.nbPt - 1];
 519
 520                 curP[Geom::X] = data.Xk[i];
 521                 curP[Geom::Y] = data.Yk[i];
 522                 double mtk = 0.5 * (data.tk[i] + data.tk[i - 1]);
 523
 524                 midAppP[Geom::X] = N13(mtk) * cp1[Geom::X] +
 525                     N23(mtk) * cp2[Geom::X] +
 526                     N03(mtk) * data.Xk[0] +
 527                     N33(mtk) * data.Xk[data.nbPt - 1];
 528
 529                 midAppP[Geom::Y] = N13(mtk) * cp1[Geom::Y] +
 530                     N23(mtk) * cp2[Geom::Y] +
 531                     N03(mtk) * data.Yk[0] +
 532                     N33(mtk) * data.Yk[data.nbPt - 1];
 533
 534                 midP = 0.5 * (curP + prevP);
 535
 536                 Geom::Point diff = curAppP - curP;
 537                 curDist = dot(diff, diff);
 538                 diff = midAppP - midP;
 539                 midDist = dot(diff, diff);
 540
 541                 delta += 0.3333 * (curDist + prevDist + midDist) * data.lk[i];
 542                 if ( curDist > worstD ) {
 543                     worstD = curDist;
 544                     worstP = i;
 545                 } else if ( data.fk[i] && 2 * curDist > worstD ) {
 546                     worstD = 2*curDist;
 547                     worstP = i;
 548                 }
 549                 prevP = curP;
 550                 prevAppP = curAppP;
 551                 prevDist = curDist;
 552             }
 553             delta /= data.totLen;
 554
 555         } else {
 556 #endif
 557             for (int i = 1; i < data.nbPt - 1; i++) {
 558                 Geom::Point curAppP;
 559                 Geom::Point curP;
 560                 double    curDist;
 561
 562                 curAppP[Geom::X] = N13(data.tk[i]) * cp1[Geom::X] +
 563                     N23(data.tk[i]) * cp2[Geom::X] +
 564                     N03(data.tk[i]) * data.Xk[0] +
 565                     N33(data.tk[i]) * data.Xk[data.nbPt - 1];
 566
 567                 curAppP[Geom::Y] = N13(data.tk[i]) * cp1[Geom::Y] +
 568                     N23(data.tk[i]) * cp2[Geom::Y] +
 569                     N03(data.tk[i]) * data.Yk[0] +
 570                     N33(data.tk[i]) * data.Yk[data.nbPt - 1];
 571
 572                 curP[Geom::X] = data.Xk[i];
 573                 curP[Geom::Y] = data.Yk[i];
 574
 575                 Geom::Point diff = curAppP-curP;
 576                 curDist = dot(diff, diff);
 577                 delta += curDist;
 578
 579                 if ( curDist > worstD ) {
 580                     worstD = curDist;
 581                     worstP = i;
 582                 } else if ( data.fk[i] && 2 * curDist > worstD ) {
 583                     worstD = 2*curDist;
 584                     worstP = i;
 585                 }
 586                 prevP = curP;
 587                 prevAppP = curAppP;
 588                 prevDist = curDist;
 589             }
 590 #ifdef with_splotch_killer
 591         }
 592 #endif
 593     }
 594
 595     if (delta < treshhold * treshhold) {
 596         // premier jet
 597
 598         // Refine a little.
 599         for (int i = 1; i < data.nbPt - 1; i++) {
 600             Geom::Point pt(data.Xk[i], data.Yk[i]);
 601             data.tk[i] = RaffineTk(pt, start, cp1, cp2, end, data.tk[i]);
 602             if (data.tk[i] < data.tk[i - 1]) {
 603                 // Force tk to be monotonic non-decreasing.
 604                 data.tk[i] = data.tk[i - 1];
 605             }
 606         }
 607
 608         if ( FitCubic(start, res, data.Xk, data.Yk, data.Qk, data.tk, data.nbPt) == false) {
 609             // ca devrait jamais arriver, mais bon
 610             res.start = 3.0 * (cp1 - start);
 611             res.end = 3.0 * (end - cp2 );
 612             return true;
 613         }
 614
 615         double ndelta = 0;
 616         {
 617             double worstD = 0;
 618             worstP = -1;
 619             Geom::Point prevAppP;
 620             Geom::Point prevP(data.Xk[0], data.Yk[0]);
 621             double prevDist = 0;
 622             prevAppP = prevP; // le premier seulement
 623 #ifdef with_splotch_killer
 624             if ( data.nbPt <= 20 ) {
 625                 for (int i = 1; i < data.nbPt - 1; i++) {
 626                     Geom::Point curAppP;
 627                     Geom::Point curP;
 628                     double  curDist;
 629                     Geom::Point midAppP;
 630                     Geom::Point midP;
 631                     double  midDist;
 632
 633                     curAppP[Geom::X] = N13(data.tk[i]) * cp1[Geom::X] +
 634                         N23(data.tk[i]) * cp2[Geom::X] +
 635                         N03(data.tk[i]) * data.Xk[0] +
 636                         N33(data.tk[i]) * data.Xk[data.nbPt - 1];
 637
 638                     curAppP[Geom::Y] = N13(data.tk[i]) * cp1[Geom::Y] +
 639                         N23(data.tk[i]) * cp2[Geom::Y] +
 640                         N03(data.tk[i]) * data.Yk[0] +
 641                         N33(data.tk[i]) * data.Yk[data.nbPt - 1];
 642
 643                     curP[Geom::X] = data.Xk[i];
 644                     curP[Geom::Y] = data.Yk[i];
 645                     double mtk = 0.5 * (data.tk[i] + data.tk[i - 1]);
 646
 647                     midAppP[Geom::X] = N13(mtk) * cp1[Geom::X] +
 648                         N23(mtk) * cp2[Geom::X] +
 649                         N03(mtk) * data.Xk[0] +
 650                         N33(mtk) * data.Xk[data.nbPt - 1];
 651
 652                     midAppP[Geom::Y] = N13(mtk) * cp1[Geom::Y] +
 653                         N23(mtk) * cp2[Geom::Y] +
 654                         N03(mtk) * data.Yk[0] +
 655                         N33(mtk) * data.Yk[data.nbPt - 1];
 656
 657                     midP = 0.5 * (curP + prevP);
 658
 659                     Geom::Point diff = curAppP - curP;
 660                     curDist = dot(diff, diff);
 661
 662                     diff = midAppP - midP;
 663                     midDist = dot(diff, diff);
 664
 665                     ndelta += 0.3333 * (curDist + prevDist + midDist) * data.lk[i];
 666
 667                     if ( curDist > worstD ) {
 668                         worstD = curDist;
 669                         worstP = i;
 670                     } else if ( data.fk[i] && 2 * curDist > worstD ) {
 671                         worstD = 2*curDist;
 672                         worstP = i;
 673                     }
 674
 675                     prevP = curP;
 676                     prevAppP = curAppP;
 677                     prevDist = curDist;
 678                 }
 679                 ndelta /= data.totLen;
 680             } else {
 681 #endif
 682                 for (int i = 1; i < data.nbPt - 1; i++) {
 683                     Geom::Point curAppP;
 684                     Geom::Point curP;
 685                     double    curDist;
 686
 687                     curAppP[Geom::X] = N13(data.tk[i]) * cp1[Geom::X] +
 688                         N23(data.tk[i]) * cp2[Geom::X] +
 689                         N03(data.tk[i]) * data.Xk[0] +
 690                         N33(data.tk[i]) * data.Xk[data.nbPt - 1];
 691
 692                     curAppP[Geom::Y] = N13(data.tk[i]) * cp1[Geom::Y] +
 693                         N23(data.tk[i]) * cp2[1] +
 694                         N03(data.tk[i]) * data.Yk[0] +
 695                         N33(data.tk[i]) * data.Yk[data.nbPt - 1];
 696
 697                     curP[Geom::X] = data.Xk[i];
 698                     curP[Geom::Y] = data.Yk[i];
 699
 700                     Geom::Point diff = curAppP - curP;
 701                     curDist = dot(diff, diff);
 702
 703                     ndelta += curDist;
 704
 705                     if ( curDist > worstD ) {
 706                         worstD = curDist;
 707                         worstP = i;
 708                     } else if ( data.fk[i] && 2 * curDist > worstD ) {
 709                         worstD = 2 * curDist;
 710                         worstP = i;
 711                     }
 712                     prevP = curP;
 713                     prevAppP = curAppP;
 714                     prevDist = curDist;
 715                 }
 716 #ifdef with_splotch_killer
 717             }
 718 #endif
 719         }
 720
 721         if (ndelta < delta + 0.00001) {
 722             return true;
 723         } else {
 724             // nothing better to do
 725             res.start = 3.0 * (cp1 - start);
 726             res.end = 3.0 * (end - cp2 );
 727         }
 728
 729         return true;
 730     }
 731
 732   return false;
 733 }
 734
 735
 736 bool Path::AttemptSimplify(int off, int N, double treshhold, PathDescrCubicTo &res,int &worstP)
 737 {
 738     Geom::Point start;
 739     Geom::Point end;
 740
 741     // pour une coordonnee
 742     double *Xk;                         // la coordonnee traitee (x puis y)
 743     double *Yk;                         // la coordonnee traitee (x puis y)
 744     double *lk;                         // les longueurs de chaque segment
 745     double *tk;                         // les tk
 746     double *Qk;                         // les Qk
 747     char *fk;       // si point force
 748
 749     Geom::Point cp1;
 750     Geom::Point cp2;
 751
 752     if (N == 2) {
 753         worstP = 1;
 754         return true;
 755     }
 756
 757     start = pts[off].p;
 758     cp1 = pts[off + 1].p;
 759     end = pts[off + N - 1].p;
 760
 761     res.p = end;
 762     res.start[0] = res.start[1] = 0;
 763     res.end[0] = res.end[1] = 0;
 764     if (N == 3) {
 765         // start -> cp1 -> end
 766         res.start = cp1 - start;
 767         res.end = end - cp1;
 768         worstP = 1;
 769         return true;
 770     }
 771
 772     // Totally inefficient, allocates & deallocates all the time.
 773     tk = (double *) g_malloc(N * sizeof(double));
 774     Qk = (double *) g_malloc(N * sizeof(double));
 775     Xk = (double *) g_malloc(N * sizeof(double));
 776     Yk = (double *) g_malloc(N * sizeof(double));
 777     lk = (double *) g_malloc(N * sizeof(double));
 778     fk = (char *) g_malloc(N * sizeof(char));
 779
 780     // chord length method
 781     tk[0] = 0.0;
 782     lk[0] = 0.0;
 783     {
 784         Geom::Point prevP = start;
 785         for (int i = 1; i < N; i++) {
 786             Xk[i] = pts[off + i].p[Geom::X];
 787             Yk[i] = pts[off + i].p[Geom::Y];
 788
 789             if ( pts[off + i].isMoveTo == polyline_forced ) {
 790                 fk[i] = 0x01;
 791             } else {
 792                 fk[i] = 0;
 793             }
 794
 795             Geom::Point diff(Xk[i] - prevP[Geom::X], Yk[i] - prevP[1]);
 796             prevP[0] = Xk[i];
 797             prevP[1] = Yk[i];
 798             lk[i] = Geom::L2(diff);
 799             tk[i] = tk[i - 1] + lk[i];
 800         }
 801     }
 802
 803     if (tk[N - 1] < 0.00001) {
 804         // longueur nulle
 805         res.start[0] = res.start[1] = 0;
 806         res.end[0] = res.end[1] = 0;
 807         double worstD = 0;
 808         worstP = -1;
 809         for (int i = 1; i < N; i++) {
 810             Geom::Point nPt;
 811             bool isForced = fk[i];
 812             nPt[0] = Xk[i];
 813             nPt[1] = Yk[i];
 814
 815             double nle = DistanceToCubic(start, res, nPt);
 816             if ( isForced ) {
 817                 // forced points are favored for splitting the recursion; we do this by increasing their distance
 818                 if ( worstP < 0 || 2 * nle > worstD ) {
 819                     worstP = i;
 820                     worstD = 2 * nle;
 821                 }
 822             } else {
 823                 if ( worstP < 0 || nle > worstD ) {
 824                     worstP = i;
 825                     worstD = nle;
 826                 }
 827             }
 828         }
 829
 830         g_free(tk);
 831         g_free(Qk);
 832         g_free(Xk);
 833         g_free(Yk);
 834         g_free(fk);
 835         g_free(lk);
 836
 837         return false;
 838     }
 839
 840     double totLen = tk[N - 1];
 841     for (int i = 1; i < N - 1; i++) {
 842         tk[i] /= totLen;
 843     }
 844
 845     res.p = end;
 846     if ( FitCubic(start, res, Xk, Yk, Qk, tk, N) ) {
 847         cp1 = start + res.start / 3;
 848         cp2 = end + res.end / 3;
 849     } else {
 850         // aie, non-inversible
 851         res.start[0] = res.start[1] = 0;
 852         res.end[0] = res.end[1] = 0;
 853         double worstD = 0;
 854         worstP = -1;
 855         for (int i = 1; i < N; i++) {
 856             Geom::Point nPt(Xk[i], Yk[i]);
 857             bool isForced = fk[i];
 858             double nle = DistanceToCubic(start, res, nPt);
 859             if ( isForced ) {
 860                 // forced points are favored for splitting the recursion; we do this by increasing their distance
 861                 if ( worstP < 0 || 2 * nle > worstD ) {
 862                     worstP = i;
 863                     worstD = 2 * nle;
 864                 }
 865             } else {
 866                 if ( worstP < 0 || nle > worstD ) {
 867                     worstP = i;
 868                     worstD = nle;
 869                 }
 870             }
 871         }
 872
 873         g_free(tk);
 874         g_free(Qk);
 875         g_free(Xk);
 876         g_free(Yk);
 877         g_free(fk);
 878         g_free(lk);
 879         return false;
 880     }
 881
 882     // calcul du delta= pondere par les longueurs des segments
 883     double delta = 0;
 884     {
 885         double worstD = 0;
 886         worstP = -1;
 887         Geom::Point prevAppP;
 888     Geom::Point   prevP;
 889     double      prevDist;
 890     prevP[0] = Xk[0];
 891     prevP[1] = Yk[0];
 892     prevAppP = prevP; // le premier seulement
 893     prevDist = 0;
 894 #ifdef with_splotch_killer
 895     if ( N <= 20 ) {
 896       for (int i = 1; i < N - 1; i++)
 897       {
 898         Geom::Point curAppP;
 899         Geom::Point curP;
 900         double    curDist;
 901         Geom::Point midAppP;
 902         Geom::Point midP;
 903         double    midDist;
 904
 905         curAppP[0] = N13 (tk[i]) * cp1[0] + N23 (tk[i]) * cp2[0] + N03 (tk[i]) * Xk[0] + N33 (tk[i]) * Xk[N - 1];
 906         curAppP[1] = N13 (tk[i]) * cp1[1] + N23 (tk[i]) * cp2[1] + N03 (tk[i]) * Yk[0] + N33 (tk[i]) * Yk[N - 1];
 907         curP[0] = Xk[i];
 908         curP[1] = Yk[i];
 909         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];
 910         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];
 911         midP=0.5*(curP+prevP);
 912
 913         Geom::Point diff;
 914         diff = curAppP-curP;
 915         curDist = dot(diff,diff);
 916
 917         diff = midAppP-midP;
 918         midDist = dot(diff,diff);
 919
 920         delta+=0.3333*(curDist+prevDist+midDist)/**lk[i]*/;
 921
 922         if ( curDist > worstD ) {
 923           worstD = curDist;
 924           worstP = i;
 925         } else if ( fk[i] && 2*curDist > worstD ) {
 926           worstD = 2*curDist;
 927           worstP = i;
 928         }
 929         prevP = curP;
 930         prevAppP = curAppP;
 931         prevDist = curDist;
 932       }
 933       delta/=totLen;
 934     } else {
 935 #endif
 936       for (int i = 1; i < N - 1; i++)
 937       {
 938         Geom::Point curAppP;
 939         Geom::Point curP;
 940         double    curDist;
 941
 942         curAppP[0] = N13 (tk[i]) * cp1[0] + N23 (tk[i]) * cp2[0] + N03 (tk[i]) * Xk[0] + N33 (tk[i]) * Xk[N - 1];
 943         curAppP[1] = N13 (tk[i]) * cp1[1] + N23 (tk[i]) * cp2[1] + N03 (tk[i]) * Yk[0] + N33 (tk[i]) * Yk[N - 1];
 944         curP[0] = Xk[i];
 945         curP[1] = Yk[i];
 946
 947         Geom::Point diff;
 948         diff = curAppP-curP;
 949         curDist = dot(diff,diff);
 950         delta += curDist;
 951         if ( curDist > worstD ) {
 952           worstD = curDist;
 953           worstP = i;
 954         } else if ( fk[i] && 2*curDist > worstD ) {
 955           worstD = 2*curDist;
 956           worstP = i;
 957         }
 958         prevP = curP;
 959         prevAppP = curAppP;
 960         prevDist = curDist;
 961       }
 962 #ifdef with_splotch_killer
 963     }
 964 #endif
 965   }
 966
 967   if (delta < treshhold * treshhold)
 968   {
 969     // premier jet
 970     res.start = 3.0 * (cp1 - start);
 971     res.end = -3.0 * (cp2 - end);
 972     res.p = end;
 973
 974     // Refine a little.
 975     for (int i = 1; i < N - 1; i++)
 976     {
 977       Geom::Point
 978             pt;
 979       pt[0] = Xk[i];
 980       pt[1] = Yk[i];
 981       tk[i] = RaffineTk (pt, start, cp1, cp2, end, tk[i]);
 982       if (tk[i] < tk[i - 1])
 983             {
 984               // Force tk to be monotonic non-decreasing.
 985               tk[i] = tk[i - 1];
 986             }
 987     }
 988
 989     if ( FitCubic(start,res,Xk,Yk,Qk,tk,N) ) {
 990     } else {
 991       // ca devrait jamais arriver, mais bon
 992       res.start = 3.0 * (cp1 - start);
 993       res.end = -3.0 * (cp2 - end);
 994       g_free(tk);
 995       g_free(Qk);
 996       g_free(Xk);
 997       g_free(Yk);
 998       g_free(fk);
 999       g_free(lk);
1000       return true;
1001     }
1002     double ndelta = 0;
1003     {
1004       double worstD = 0;
1005       worstP = -1;
1006       Geom::Point   prevAppP;
1007       Geom::Point   prevP;
1008       double      prevDist;
1009       prevP[0] = Xk[0];
1010       prevP[1] = Yk[0];
1011       prevAppP = prevP; // le premier seulement
1012       prevDist = 0;
1013 #ifdef with_splotch_killer
1014       if ( N <= 20 ) {
1015         for (int i = 1; i < N - 1; i++)
1016         {
1017           Geom::Point curAppP;
1018           Geom::Point curP;
1019           double    curDist;
1020           Geom::Point midAppP;
1021           Geom::Point midP;
1022           double    midDist;
1023
1024           curAppP[0] = N13 (tk[i]) * cp1[0] + N23 (tk[i]) * cp2[0] + N03 (tk[i]) * Xk[0] + N33 (tk[i]) * Xk[N - 1];
1025           curAppP[1] = N13 (tk[i]) * cp1[1] + N23 (tk[i]) * cp2[1] + N03 (tk[i]) * Yk[0] + N33 (tk[i]) * Yk[N - 1];
1026           curP[0] = Xk[i];
1027           curP[1] = Yk[i];
1028           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];
1029           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];
1030           midP = 0.5*(curP+prevP);
1031
1032           Geom::Point diff;
1033           diff = curAppP-curP;
1034           curDist = dot(diff,diff);
1035           diff = midAppP-midP;
1036           midDist = dot(diff,diff);
1037
1038           ndelta+=0.3333*(curDist+prevDist+midDist)/**lk[i]*/;
1039
1040           if ( curDist > worstD ) {
1041             worstD = curDist;
1042             worstP = i;
1043           } else if ( fk[i] && 2*curDist > worstD ) {
1044             worstD = 2*curDist;
1045             worstP = i;
1046           }
1047           prevP = curP;
1048           prevAppP = curAppP;
1049           prevDist = curDist;
1050         }
1051         ndelta /= totLen;
1052       } else {
1053 #endif
1054         for (int i = 1; i < N - 1; i++)
1055         {
1056           Geom::Point curAppP;
1057           Geom::Point curP;
1058           double    curDist;
1059
1060           curAppP[0] = N13 (tk[i]) * cp1[0] + N23 (tk[i]) * cp2[0] + N03 (tk[i]) * Xk[0] + N33 (tk[i]) * Xk[N - 1];
1061           curAppP[1] = N13 (tk[i]) * cp1[1] + N23 (tk[i]) * cp2[1] + N03 (tk[i]) * Yk[0] + N33 (tk[i]) * Yk[N - 1];
1062           curP[0]=Xk[i];
1063           curP[1]=Yk[i];
1064
1065           Geom::Point diff;
1066           diff=curAppP-curP;
1067           curDist=dot(diff,diff);
1068           ndelta+=curDist;
1069
1070           if ( curDist > worstD ) {
1071             worstD=curDist;
1072             worstP=i;
1073           } else if ( fk[i] && 2*curDist > worstD ) {
1074             worstD=2*curDist;
1075             worstP=i;
1076           }
1077           prevP=curP;
1078           prevAppP=curAppP;
1079           prevDist=curDist;
1080         }
1081 #ifdef with_splotch_killer
1082       }
1083 #endif
1084     }
1085
1086     g_free(tk);
1087     g_free(Qk);
1088     g_free(Xk);
1089     g_free(Yk);
1090     g_free(fk);
1091     g_free(lk);
1092
1093     if (ndelta < delta + 0.00001)
1094     {
1095       return true;
1096     } else {
1097       // nothing better to do
1098       res.start = 3.0 * (cp1 - start);
1099       res.end = -3.0 * (cp2 - end);
1100     }
1101     return true;
1102   } else {
1103     // nothing better to do
1104   }
1105
1106   g_free(tk);
1107   g_free(Qk);
1108   g_free(Xk);
1109   g_free(Yk);
1110   g_free(fk);
1111   g_free(lk);
1112   return false;
1113 }
1114
1115 double Path::RaffineTk (Geom::Point pt, Geom::Point p0, Geom::Point p1, Geom::Point p2, Geom::Point p3, double it)
1116 {
1117     // Refinement of the tk values.
1118     // Just one iteration of Newtow Raphson, given that we're approaching the curve anyway.
1119     // [fr: vu que de toute facon la courbe est approchC)e]
1120     double const Ax = pt[Geom::X] -
1121         p0[Geom::X] * N03(it) -
1122         p1[Geom::X] * N13(it) -
1123         p2[Geom::X] * N23(it) -
1124         p3[Geom::X] * N33(it);
1125
1126     double const Bx = (p1[Geom::X] - p0[Geom::X]) * N02(it) +
1127         (p2[Geom::X] - p1[Geom::X]) * N12(it) +
1128         (p3[Geom::X] - p2[Geom::X]) * N22(it);
1129
1130     double const Cx = (p0[Geom::X] - 2 * p1[Geom::X] + p2[Geom::X]) * N01(it) +
1131         (p3[Geom::X] - 2 * p2[Geom::X] + p1[Geom::X]) * N11(it);
1132
1133     double const Ay =  pt[Geom::Y] -
1134         p0[Geom::Y] * N03(it) -
1135         p1[Geom::Y] * N13(it) -
1136         p2[Geom::Y] * N23(it) -
1137         p3[Geom::Y] * N33(it);
1138
1139     double const By = (p1[Geom::Y] - p0[Geom::Y]) * N02(it) +
1140         (p2[Geom::Y] - p1[Geom::Y]) * N12(it) +
1141         (p3[Geom::Y] - p2[Geom::Y]) * N22(it);
1142
1143     double const Cy = (p0[Geom::Y] - 2 * p1[Geom::Y] + p2[Geom::Y]) * N01(it) +
1144         (p3[Geom::Y] - 2 * p2[Geom::Y] + p1[Geom::Y]) * N11(it);
1145
1146     double const dF = -6 * (Ax * Bx + Ay * By);
1147     double const ddF = 18 * (Bx * Bx + By * By) - 12 * (Ax * Cx + Ay * Cy);
1148     if (fabs (ddF) > 0.0000001) {
1149         return it - dF / ddF;
1150     }
1151
1152     return it;
1153 }
1154
1155 // Variation on the fitting theme: try to merge path commands into cubic bezier patches.
1156 // The goal is to reduce the number of path commands, especially when operations on path produce
1157 // lots of small path elements; ideally you could get rid of very small segments at reduced visual cost.
1158 void Path::Coalesce(double tresh)
1159 {
1160     if ( descr_flags & descr_adding_bezier ) {
1161         CancelBezier();
1162     }
1163
1164     if ( descr_flags & descr_doing_subpath ) {
1165         CloseSubpath();
1166     }
1167
1168     if (descr_cmd.size() <= 2) {
1169         return;
1170     }
1171
1172     SetBackData(false);
1173     Path* tempDest = new Path();
1174     tempDest->SetBackData(false);
1175
1176     ConvertEvenLines(0.25*tresh);
1177
1178     int lastP = 0;
1179     int lastAP = -1;
1180     // As the elements are stored in a separate array, it's no longer worth optimizing
1181     // the rewriting in the same array.
1182     // [[comme les elements sont stockes dans un tableau a part, plus la peine d'optimiser
1183     // la réécriture dans la meme tableau]]
1184
1185     int lastA = descr_cmd[0]->associated;
1186     int prevA = lastA;
1187     Geom::Point firstP;
1188
1189     /* FIXME: the use of this variable probably causes a leak or two.
1190     ** It's a hack anyway, and probably only needs to be a type rather than
1191     ** a full PathDescr.
1192     */
1193     PathDescr *lastAddition = new PathDescrMoveTo(Geom::Point(0, 0));
1194     bool containsForced = false;
1195     PathDescrCubicTo pending_cubic(Geom::Point(0, 0), Geom::Point(0, 0), Geom::Point(0, 0));
1196
1197     for (int curP = 0; curP < int(descr_cmd.size()); curP++) {
1198         int typ = descr_cmd[curP]->getType();
1199         int nextA = lastA;
1200
1201         if (typ == descr_moveto) {
1202
1203             if (lastAddition->flags != descr_moveto) {
1204                 FlushPendingAddition(tempDest,lastAddition,pending_cubic,lastAP);
1205             }
1206             lastAddition = descr_cmd[curP];
1207             lastAP = curP;
1208             FlushPendingAddition(tempDest, lastAddition, pending_cubic, lastAP);
1209             // Added automatically (too bad about multiple moveto's).
1210             // [fr: (tant pis pour les moveto multiples)]
1211             containsForced = false;
1212
1213             PathDescrMoveTo *nData = dynamic_cast<PathDescrMoveTo *>(descr_cmd[curP]);
1214             firstP = nData->p;
1215             lastA = descr_cmd[curP]->associated;
1216             prevA = lastA;
1217             lastP = curP;
1218
1219         } else if (typ == descr_close) {
1220             nextA = descr_cmd[curP]->associated;
1221             if (lastAddition->flags != descr_moveto) {
1222
1223                 PathDescrCubicTo res(Geom::Point(0, 0), Geom::Point(0, 0), Geom::Point(0, 0));
1224                 int worstP = -1;
1225                 if (AttemptSimplify(lastA, nextA - lastA + 1, (containsForced) ? 0.05 * tresh : tresh, res, worstP)) {
1226                     lastAddition = new PathDescrCubicTo(Geom::Point(0, 0),
1227                                                           Geom::Point(0, 0),
1228                                                           Geom::Point(0, 0));
1229                     pending_cubic = res;
1230                     lastAP = -1;
1231                 }
1232
1233                 FlushPendingAddition(tempDest, lastAddition, pending_cubic, lastAP);
1234                 FlushPendingAddition(tempDest, descr_cmd[curP], pending_cubic, curP);
1235
1236             } else {
1237                 FlushPendingAddition(tempDest,descr_cmd[curP],pending_cubic,curP);
1238             }
1239
1240             containsForced = false;
1241             lastAddition = new PathDescrMoveTo(Geom::Point(0, 0));
1242             prevA = lastA = nextA;
1243             lastP = curP;
1244             lastAP = curP;
1245
1246         } else if (typ == descr_forced) {
1247
1248             nextA = descr_cmd[curP]->associated;
1249             if (lastAddition->flags != descr_moveto) {
1250
1251                 PathDescrCubicTo res(Geom::Point(0, 0), Geom::Point(0, 0), Geom::Point(0, 0));
1252                 int worstP = -1;
1253                 if (AttemptSimplify(lastA, nextA - lastA + 1, 0.05 * tresh, res, worstP)) {
1254                     // plus sensible parce que point force
1255                     // ca passe
1256                     /* (Possible translation: More sensitive because contains a forced point.) */
1257                     containsForced = true;
1258                 } else  {
1259                     // Force the addition.
1260                     FlushPendingAddition(tempDest, lastAddition, pending_cubic, lastAP);
1261                     lastAddition = new PathDescrMoveTo(Geom::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(Geom::Point(0, 0), Geom::Point(0, 0), Geom::Point(0, 0));
1275                 int worstP = -1;
1276                 if (AttemptSimplify(lastA, nextA - lastA + 1, tresh, res, worstP)) {
1277                     lastAddition = new PathDescrCubicTo(Geom::Point(0, 0),
1278                                                           Geom::Point(0, 0),
1279                                                           Geom::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                         /* (possible translation: Could be overwritten by the next line.) */
1287                     FlushPendingAddition(tempDest, lastAddition, pending_cubic, lastAP);
1288                     lastAddition = descr_cmd[curP];
1289                     if ( typ == descr_cubicto ) {
1290                         pending_cubic = *(dynamic_cast<PathDescrCubicTo*>(descr_cmd[curP]));
1291                     }
1292                     lastAP = curP;
1293                     containsForced = false;
1294                 }
1295
1296             } else {
1297                 lastA = prevA /*descr_cmd[curP-1]->associated */ ;
1298                 lastAddition = descr_cmd[curP];
1299                 if ( typ == descr_cubicto ) {
1300                     pending_cubic = *(dynamic_cast<PathDescrCubicTo*>(descr_cmd[curP]));
1301                 }
1302                 lastAP = curP;
1303                 containsForced = false;
1304             }
1305             prevA = nextA;
1306
1307         } else if (typ == descr_bezierto) {
1308
1309             if (lastAddition->flags != descr_moveto) {
1310                 FlushPendingAddition(tempDest, lastAddition, pending_cubic, lastAP);
1311                 lastAddition = new PathDescrMoveTo(Geom::Point(0, 0));
1312             }
1313             lastAP = -1;
1314             lastA = descr_cmd[curP]->associated;
1315             lastP = curP;
1316             PathDescrBezierTo *nBData = dynamic_cast<PathDescrBezierTo*>(descr_cmd[curP]);
1317             for (int i = 1; i <= nBData->nb; i++) {
1318                 FlushPendingAddition(tempDest, descr_cmd[curP + i], pending_cubic, curP + i);
1319             }
1320             curP += nBData->nb;
1321             prevA = nextA;
1322
1323         } else if (typ == descr_interm_bezier) {
1324             continue;
1325         } else {
1326             continue;
1327         }
1328     }
1329
1330     if (lastAddition->flags != descr_moveto) {
1331         FlushPendingAddition(tempDest, lastAddition, pending_cubic, lastAP);
1332     }
1333
1334     Copy(tempDest);
1335     delete tempDest;
1336 }
1337
1338
1339 void Path::FlushPendingAddition(Path *dest, PathDescr *lastAddition,
1340                                 PathDescrCubicTo &lastCubic, int lastAP)
1341 {
1342     switch (lastAddition->getType()) {
1343
1344     case descr_moveto:
1345         if ( lastAP >= 0 ) {
1346             PathDescrMoveTo* nData = dynamic_cast<PathDescrMoveTo *>(descr_cmd[lastAP]);
1347             dest->MoveTo(nData->p);
1348         }
1349         break;
1350
1351     case descr_close:
1352         dest->Close();
1353         break;
1354
1355     case descr_cubicto:
1356         dest->CubicTo(lastCubic.p, lastCubic.start, lastCubic.end);
1357         break;
1358
1359     case descr_lineto:
1360         if ( lastAP >= 0 ) {
1361             PathDescrLineTo *nData = dynamic_cast<PathDescrLineTo *>(descr_cmd[lastAP]);
1362             dest->LineTo(nData->p);
1363         }
1364         break;
1365
1366     case descr_arcto:
1367         if ( lastAP >= 0 ) {
1368             PathDescrArcTo *nData = dynamic_cast<PathDescrArcTo *>(descr_cmd[lastAP]);
1369             dest->ArcTo(nData->p, nData->rx, nData->ry, nData->angle, nData->large, nData->clockwise);
1370         }
1371         break;
1372
1373     case descr_bezierto:
1374         if ( lastAP >= 0 ) {
1375             PathDescrBezierTo *nData = dynamic_cast<PathDescrBezierTo *>(descr_cmd[lastAP]);
1376             dest->BezierTo(nData->p);
1377         }
1378         break;
1379
1380     case descr_interm_bezier:
1381         if ( lastAP >= 0 ) {
1382             PathDescrIntermBezierTo *nData = dynamic_cast<PathDescrIntermBezierTo*>(descr_cmd[lastAP]);
1383             dest->IntermBezierTo(nData->p);
1384         }
1385         break;
1386     }
1387 }
1388
1389 /*
1390   Local Variables:
1391   mode:c++
1392   c-file-style:"stroustrup"
1393   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
1394   indent-tabs-mode:nil
1395   fill-column:99
1396   End:
1397 */
1398 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :