src/trace/potrace/trace.cpp

   1 /* Copyright (C) 2001-2005 Peter Selinger.
   2    This file is part of potrace. It is free software and it is covered
   3    by the GNU General Public License. See the file COPYING for details. */
   4
   5 /* $Id$ */
   6 /* transform jaggy paths into smooth curves */
   7
   8 #include <math.h>
   9 #include <stdlib.h>
  10
  11 #include "curve.h"
  12 #include "lists.h"
  13 #include "progress.h"
  14
  15 #define INFTY 10000000  /* it suffices that this is longer than any
  16                            path; it need not be really infinite */
  17 #define COS179 -0.999847695156   /* the cosine of 179 degrees */
  18
  19 /* ---------------------------------------------------------------------- */
  20 /* some useful macros. Note: the "mod" macro works correctly for
  21    negative a. Also note that the test for a>=n, while redundant,
  22    speeds up the mod function by 70% in the average case (significant
  23    since the program spends about 16% of its time here - or 40%
  24    without the test). The "floordiv" macro returns the largest integer
  25    <= a/n, and again this works correctly for negative a, as long as
  26    a,n are integers and n>0. */
  27
  28 #define SAFE_MALLOC(var, n, typ) \
  29   if ((var = (typ *)malloc((n)*sizeof(typ))) == NULL) goto malloc_error
  30
  31 /* ---------------------------------------------------------------------- */
  32 /* auxiliary functions */
  33
  34 /* return a direction that is 90 degrees counterclockwise from p2-p0,
  35    but then restricted to one of the major wind directions (n, nw, w, etc) */
  36 static inline point_t dorth_infty(dpoint_t p0, dpoint_t p2) {
  37   point_t r;
  38
  39   r.y = sign(p2.x-p0.x);
  40   r.x = -sign(p2.y-p0.y);
  41
  42   return r;
  43 }
  44
  45 /* return (p1-p0)x(p2-p0), the area of the parallelogram */
  46 static inline double dpara(dpoint_t p0, dpoint_t p1, dpoint_t p2) {
  47   double x1, y1, x2, y2;
  48
  49   x1 = p1.x-p0.x;
  50   y1 = p1.y-p0.y;
  51   x2 = p2.x-p0.x;
  52   y2 = p2.y-p0.y;
  53
  54   return x1*y2 - x2*y1;
  55 }
  56
  57 /* ddenom/dpara have the property that the square of radius 1 centered
  58    at p1 intersects the line p0p2 iff |dpara(p0,p1,p2)| <= ddenom(p0,p2) */
  59 static inline double ddenom(dpoint_t p0, dpoint_t p2) {
  60   point_t r = dorth_infty(p0, p2);
  61
  62   return r.y*(p2.x-p0.x) - r.x*(p2.y-p0.y);
  63 }
  64
  65 /* return 1 if a <= b < c < a, in a cyclic sense (mod n) */
  66 static inline int cyclic(int a, int b, int c) {
  67   if (a<=c) {
  68     return (a<=b && b<c);
  69   } else {
  70     return (a<=b || b<c);
  71   }
  72 }
  73
  74 /* determine the center and slope of the line i..j. Assume i<j. Needs
  75    "sum" components of p to be set. */
  76 static void pointslope(privpath_t *pp, int i, int j, dpoint_t *ctr, dpoint_t *dir) {
  77   /* assume i<j */
  78
  79   int n = pp->len;
  80   sums_t *sums = pp->sums;
  81
  82   double x, y, x2, xy, y2;
  83   double k;
  84   double a, b, c, lambda2, l;
  85   int r=0; /* rotations from i to j */
  86
  87   while (j>=n) {
  88     j-=n;
  89     r+=1;
  90   }
  91   while (i>=n) {
  92     i-=n;
  93     r-=1;
  94   }
  95   while (j<0) {
  96     j+=n;
  97     r-=1;
  98   }
  99   while (i<0) {
 100     i+=n;
 101     r+=1;
 102   }
 103
 104   x = sums[j+1].x-sums[i].x+r*sums[n].x;
 105   y = sums[j+1].y-sums[i].y+r*sums[n].y;
 106   x2 = sums[j+1].x2-sums[i].x2+r*sums[n].x2;
 107   xy = sums[j+1].xy-sums[i].xy+r*sums[n].xy;
 108   y2 = sums[j+1].y2-sums[i].y2+r*sums[n].y2;
 109   k = j+1-i+r*n;
 110
 111   ctr->x = x/k;
 112   ctr->y = y/k;
 113
 114   a = (x2-(double)x*x/k)/k;
 115   b = (xy-(double)x*y/k)/k;
 116   c = (y2-(double)y*y/k)/k;
 117
 118   lambda2 = (a+c+sqrt((a-c)*(a-c)+4*b*b))/2; /* larger e.value */
 119
 120   /* now find e.vector for lambda2 */
 121   a -= lambda2;
 122   c -= lambda2;
 123
 124   if (fabs(a) >= fabs(c)) {
 125     l = sqrt(a*a+b*b);
 126     if (l!=0) {
 127       dir->x = -b/l;
 128       dir->y = a/l;
 129     }
 130   } else {
 131     l = sqrt(c*c+b*b);
 132     if (l!=0) {
 133       dir->x = -c/l;
 134       dir->y = b/l;
 135     }
 136   }
 137   if (l==0) {
 138     dir->x = dir->y = 0;   /* sometimes this can happen when k=4:
 139                               the two eigenvalues coincide */
 140   }
 141 }
 142
 143 /* the type of (affine) quadratic forms, represented as symmetric 3x3
 144    matrices.  The value of the quadratic form at a vector (x,y) is v^t
 145    Q v, where v = (x,y,1)^t. */
 146 typedef double quadform_t[3][3];
 147
 148 /* Apply quadratic form Q to vector w = (w.x,w.y) */
 149 static inline double quadform(quadform_t Q, dpoint_t w) {
 150   double v[3];
 151   int i, j;
 152   double sum;
 153
 154   v[0] = w.x;
 155   v[1] = w.y;
 156   v[2] = 1;
 157   sum = 0.0;
 158
 159   for (i=0; i<3; i++) {
 160     for (j=0; j<3; j++) {
 161       sum += v[i] * Q[i][j] * v[j];
 162     }
 163   }
 164   return sum;
 165 }
 166
 167 /* calculate p1 x p2 */
 168 static inline int xprod(point_t p1, point_t p2) {
 169   return p1.x*p2.y - p1.y*p2.x;
 170 }
 171
 172 /* calculate (p1-p0)x(p3-p2) */
 173 static inline double cprod(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3) {
 174   double x1, y1, x2, y2;
 175
 176   x1 = p1.x - p0.x;
 177   y1 = p1.y - p0.y;
 178   x2 = p3.x - p2.x;
 179   y2 = p3.y - p2.y;
 180
 181   return x1*y2 - x2*y1;
 182 }
 183
 184 /* calculate (p1-p0)*(p2-p0) */
 185 static inline double iprod(dpoint_t p0, dpoint_t p1, dpoint_t p2) {
 186   double x1, y1, x2, y2;
 187
 188   x1 = p1.x - p0.x;
 189   y1 = p1.y - p0.y;
 190   x2 = p2.x - p0.x;
 191   y2 = p2.y - p0.y;
 192
 193   return x1*x2 + y1*y2;
 194 }
 195
 196 /* calculate (p1-p0)*(p3-p2) */
 197 static inline double iprod1(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3) {
 198   double x1, y1, x2, y2;
 199
 200   x1 = p1.x - p0.x;
 201   y1 = p1.y - p0.y;
 202   x2 = p3.x - p2.x;
 203   y2 = p3.y - p2.y;
 204
 205   return x1*x2 + y1*y2;
 206 }
 207
 208 /* calculate distance between two points */
 209 static inline double ddist(dpoint_t p, dpoint_t q) {
 210   return sqrt(sq(p.x-q.x)+sq(p.y-q.y));
 211 }
 212
 213 /* calculate point of a bezier curve */
 214 static inline dpoint_t bezier(double t, dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3) {
 215   double s = 1-t;
 216   dpoint_t res;
 217
 218   /* Note: a good optimizing compiler (such as gcc-3) reduces the
 219      following to 16 multiplications, using common subexpression
 220      elimination. */
 221
 222   res.x = s*s*s*p0.x + 3*(s*s*t)*p1.x + 3*(t*t*s)*p2.x + t*t*t*p3.x;
 223   res.y = s*s*s*p0.y + 3*(s*s*t)*p1.y + 3*(t*t*s)*p2.y + t*t*t*p3.y;
 224
 225   return res;
 226 }
 227
 228 /* calculate the point t in [0..1] on the (convex) bezier curve
 229    (p0,p1,p2,p3) which is tangent to q1-q0. Return -1.0 if there is no
 230    solution in [0..1]. */
 231 static double tangent(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3, dpoint_t q0, dpoint_t q1) {
 232   double A, B, C;   /* (1-t)^2 A + 2(1-t)t B + t^2 C = 0 */
 233   double a, b, c;   /* a t^2 + b t + c = 0 */
 234   double d, s, r1, r2;
 235
 236   A = cprod(p0, p1, q0, q1);
 237   B = cprod(p1, p2, q0, q1);
 238   C = cprod(p2, p3, q0, q1);
 239
 240   a = A - 2*B + C;
 241   b = -2*A + 2*B;
 242   c = A;
 243
 244   d = b*b - 4*a*c;
 245
 246   if (a==0 || d<0) {
 247     return -1.0;
 248   }
 249
 250   s = sqrt(d);
 251
 252   r1 = (-b + s) / (2 * a);
 253   r2 = (-b - s) / (2 * a);
 254
 255   if (r1 >= 0 && r1 <= 1) {
 256     return r1;
 257   } else if (r2 >= 0 && r2 <= 1) {
 258     return r2;
 259   } else {
 260     return -1.0;
 261   }
 262 }
 263
 264 /* ---------------------------------------------------------------------- */
 265 /* Preparation: fill in the sum* fields of a path (used for later
 266    rapid summing). Return 0 on success, 1 with errno set on
 267    failure. */
 268 static int calc_sums(privpath_t *pp) {
 269   int i, x, y;
 270   int n = pp->len;
 271
 272   SAFE_MALLOC(pp->sums, pp->len+1, sums_t);
 273
 274   /* origin */
 275   pp->x0 = pp->pt[0].x;
 276   pp->y0 = pp->pt[0].y;
 277
 278   /* preparatory computation for later fast summing */
 279   pp->sums[0].x2 = pp->sums[0].xy = pp->sums[0].y2 = pp->sums[0].x = pp->sums[0].y = 0;
 280   for (i=0; i<n; i++) {
 281     x = pp->pt[i].x - pp->x0;
 282     y = pp->pt[i].y - pp->y0;
 283     pp->sums[i+1].x = pp->sums[i].x + x;
 284     pp->sums[i+1].y = pp->sums[i].y + y;
 285     pp->sums[i+1].x2 = pp->sums[i].x2 + x*x;
 286     pp->sums[i+1].xy = pp->sums[i].xy + x*y;
 287     pp->sums[i+1].y2 = pp->sums[i].y2 + y*y;
 288   }
 289   return 0;
 290
 291  malloc_error:
 292   return 1;
 293 }
 294
 295 /* ---------------------------------------------------------------------- */
 296 /* Stage 1: determine the straight subpaths (Sec. 2.2.1). Fill in the
 297    "lon" component of a path object (based on pt/len).  For each i,
 298    lon[i] is the furthest index such that a straight line can be drawn
 299    from i to lon[i]. Return 1 on error with errno set, else 0. */
 300
 301 /* this algorithm depends on the fact that the existence of straight
 302    subpaths is a triplewise property. I.e., there exists a straight
 303    line through squares i0,...,in iff there exists a straight line
 304    through i,j,k, for all i0<=i<j<k<=in. (Proof?) */
 305
 306 /* this implementation of calc_lon is O(n^2). It replaces an older
 307    O(n^3) version. A "constraint" means that future points must
 308    satisfy xprod(constraint[0], cur) >= 0 and xprod(constraint[1],
 309    cur) <= 0. */
 310
 311 /* Remark for potrace 1.1: the current implementation of calc_lon is
 312    more complex than the implementation found in potrace 1.0, but it
 313    is considerably faster. The introduction of the "nc" data structure
 314    means that we only have to test the constraints for "corner"
 315    points. On a typical input file, this speeds up the calc_lon
 316    function by a factor of 31.2, thereby decreasing its time share
 317    within the overall potrace algorithm from 72.6% to 7.82%, and
 318    speeding up the overall algorithm by a factor of 3.36. On another
 319    input file, calc_lon was sped up by a factor of 6.7, decreasing its
 320    time share from 51.4% to 13.61%, and speeding up the overall
 321    algorithm by a factor of 1.78. In any case, the savings are
 322    substantial. */
 323
 324 /* returns 0 on success, 1 on error with errno set */
 325 static int calc_lon(privpath_t *pp) {
 326   point_t *pt = pp->pt;
 327   int n = pp->len;
 328   int i, j, k, k1;
 329   int ct[4], dir;
 330   point_t constraint[2];
 331   point_t cur;
 332   point_t off;
 333   int *pivk = NULL;  /* pivk[n] */
 334   int *nc = NULL;    /* nc[n]: next corner */
 335   point_t dk;  /* direction of k-k1 */
 336   int a, b, c, d;
 337
 338   SAFE_MALLOC(pivk, n, int);
 339   SAFE_MALLOC(nc, n, int);
 340
 341   /* initialize the nc data structure. Point from each point to the
 342      furthest future point to which it is connected by a vertical or
 343      horizontal segment. We take advantage of the fact that there is
 344      always a direction change at 0 (due to the path decomposition
 345      algorithm). But even if this were not so, there is no harm, as
 346      in practice, correctness does not depend on the word "furthest"
 347      above.  */
 348   k = 0;
 349   for (i=n-1; i>=0; i--) {
 350     if (pt[i].x != pt[k].x && pt[i].y != pt[k].y) {
 351       k = i+1;  /* necessarily i<n-1 in this case */
 352     }
 353     nc[i] = k;
 354   }
 355
 356   SAFE_MALLOC(pp->lon, n, int);
 357
 358   /* determine pivot points: for each i, let pivk[i] be the furthest k
 359      such that all j with i<j<k lie on a line connecting i,k. */
 360
 361   for (i=n-1; i>=0; i--) {
 362     ct[0] = ct[1] = ct[2] = ct[3] = 0;
 363
 364     /* keep track of "directions" that have occurred */
 365     dir = (3+3*(pt[mod(i+1,n)].x-pt[i].x)+(pt[mod(i+1,n)].y-pt[i].y))/2;
 366     ct[dir]++;
 367
 368     constraint[0].x = 0;
 369     constraint[0].y = 0;
 370     constraint[1].x = 0;
 371     constraint[1].y = 0;
 372
 373     /* find the next k such that no straight line from i to k */
 374     k = nc[i];
 375     k1 = i;
 376     while (1) {
 377
 378       dir = (3+3*sign(pt[k].x-pt[k1].x)+sign(pt[k].y-pt[k1].y))/2;
 379       ct[dir]++;
 380
 381       /* if all four "directions" have occurred, cut this path */
 382       if (ct[0] && ct[1] && ct[2] && ct[3]) {
 383         pivk[i] = k1;
 384         goto foundk;
 385       }
 386
 387       cur.x = pt[k].x - pt[i].x;
 388       cur.y = pt[k].y - pt[i].y;
 389
 390       /* see if current constraint is violated */
 391       if (xprod(constraint[0], cur) < 0 || xprod(constraint[1], cur) > 0) {
 392         goto constraint_viol;
 393       }
 394
 395       /* else, update constraint */
 396       if (abs(cur.x) <= 1 && abs(cur.y) <= 1) {
 397         /* no constraint */
 398       } else {
 399         off.x = cur.x + ((cur.y>=0 && (cur.y>0 || cur.x<0)) ? 1 : -1);
 400         off.y = cur.y + ((cur.x<=0 && (cur.x<0 || cur.y<0)) ? 1 : -1);
 401         if (xprod(constraint[0], off) >= 0) {
 402           constraint[0] = off;
 403         }
 404         off.x = cur.x + ((cur.y<=0 && (cur.y<0 || cur.x<0)) ? 1 : -1);
 405         off.y = cur.y + ((cur.x>=0 && (cur.x>0 || cur.y<0)) ? 1 : -1);
 406         if (xprod(constraint[1], off) <= 0) {
 407           constraint[1] = off;
 408         }
 409       }
 410       k1 = k;
 411       k = nc[k1];
 412       if (!cyclic(k,i,k1)) {
 413         break;
 414       }
 415     }
 416   constraint_viol:
 417     /* k1 was the last "corner" satisfying the current constraint, and
 418        k is the first one violating it. We now need to find the last
 419        point along k1..k which satisfied the constraint. */
 420     dk.x = sign(pt[k].x-pt[k1].x);
 421     dk.y = sign(pt[k].y-pt[k1].y);
 422     cur.x = pt[k1].x - pt[i].x;
 423     cur.y = pt[k1].y - pt[i].y;
 424     /* find largest integer j such that xprod(constraint[0], cur+j*dk)
 425        >= 0 and xprod(constraint[1], cur+j*dk) <= 0. Use bilinearity
 426        of xprod. */
 427     a = xprod(constraint[0], cur);
 428     b = xprod(constraint[0], dk);
 429     c = xprod(constraint[1], cur);
 430     d = xprod(constraint[1], dk);
 431     /* find largest integer j such that a+j*b>=0 and c+j*d<=0. This
 432        can be solved with integer arithmetic. */
 433     j = INFTY;
 434     if (b<0) {
 435       j = floordiv(a,-b);
 436     }
 437     if (d>0) {
 438       j = min(j, floordiv(-c,d));
 439     }
 440     pivk[i] = mod(k1+j,n);
 441   foundk:
 442     ;
 443   } /* for i */
 444
 445   /* clean up: for each i, let lon[i] be the largest k such that for
 446      all i' with i<=i'<k, i'<k<=pivk[i']. */
 447
 448   j=pivk[n-1];
 449   pp->lon[n-1]=j;
 450   for (i=n-2; i>=0; i--) {
 451     if (cyclic(i+1,pivk[i],j)) {
 452       j=pivk[i];
 453     }
 454     pp->lon[i]=j;
 455   }
 456
 457   for (i=n-1; cyclic(mod(i+1,n),j,pp->lon[i]); i--) {
 458     pp->lon[i] = j;
 459   }
 460
 461   free(pivk);
 462   free(nc);
 463   return 0;
 464
 465  malloc_error:
 466   free(pivk);
 467   free(nc);
 468   return 1;
 469 }
 470
 471
 472 /* ---------------------------------------------------------------------- */
 473 /* Stage 2: calculate the optimal polygon (Sec. 2.2.2-2.2.4). */
 474
 475 /* Auxiliary function: calculate the penalty of an edge from i to j in
 476    the given path. This needs the "lon" and "sum*" data. */
 477
 478 static double penalty3(privpath_t *pp, int i, int j) {
 479   int n = pp->len;
 480   point_t *pt = pp->pt;
 481   sums_t *sums = pp->sums;
 482
 483   /* assume 0<=i<j<=n  */
 484   double x, y, x2, xy, y2;
 485   double k;
 486   double a, b, c, s;
 487   double px, py, ex, ey;
 488
 489   int r=0; /* rotations from i to j */
 490
 491   if (j>=n) {
 492     j-=n;
 493     r+=1;
 494   }
 495
 496   x = sums[j+1].x-sums[i].x+r*sums[n].x;
 497   y = sums[j+1].y-sums[i].y+r*sums[n].y;
 498   x2 = sums[j+1].x2-sums[i].x2+r*sums[n].x2;
 499   xy = sums[j+1].xy-sums[i].xy+r*sums[n].xy;
 500   y2 = sums[j+1].y2-sums[i].y2+r*sums[n].y2;
 501   k = j+1-i+r*n;
 502
 503   px = (pt[i].x+pt[j].x)/2.0-pt[0].x;
 504   py = (pt[i].y+pt[j].y)/2.0-pt[0].y;
 505   ey = (pt[j].x-pt[i].x);
 506   ex = -(pt[j].y-pt[i].y);
 507
 508   a = ((x2-2*x*px)/k+px*px);
 509   b = ((xy-x*py-y*px)/k+px*py);
 510   c = ((y2-2*y*py)/k+py*py);
 511
 512   s = ex*ex*a + 2*ex*ey*b + ey*ey*c;
 513
 514   return sqrt(s);
 515 }
 516
 517 /* find the optimal polygon. Fill in the m and po components. Return 1
 518    on failure with errno set, else 0. Non-cyclic version: assumes i=0
 519    is in the polygon. Fixme: ### implement cyclic version. */
 520 static int bestpolygon(privpath_t *pp)
 521 {
 522   int i, j, m, k;
 523   int n = pp->len;
 524   double *pen = NULL; /* pen[n+1]: penalty vector */
 525   int *prev = NULL;   /* prev[n+1]: best path pointer vector */
 526   int *clip0 = NULL;  /* clip0[n]: longest segment pointer, non-cyclic */
 527   int *clip1 = NULL;  /* clip1[n+1]: backwards segment pointer, non-cyclic */
 528   int *seg0 = NULL;    /* seg0[m+1]: forward segment bounds, m<=n */
 529   int *seg1 = NULL;   /* seg1[m+1]: backward segment bounds, m<=n */
 530   double thispen;
 531   double best;
 532   int c;
 533
 534   SAFE_MALLOC(pen, n+1, double);
 535   SAFE_MALLOC(prev, n+1, int);
 536   SAFE_MALLOC(clip0, n, int);
 537   SAFE_MALLOC(clip1, n+1, int);
 538   SAFE_MALLOC(seg0, n+1, int);
 539   SAFE_MALLOC(seg1, n+1, int);
 540
 541   /* calculate clipped paths */
 542   for (i=0; i<n; i++) {
 543     c = mod(pp->lon[mod(i-1,n)]-1,n);
 544     if (c == i) {
 545       c = mod(i+1,n);
 546     }
 547     if (c < i) {
 548       clip0[i] = n;
 549     } else {
 550       clip0[i] = c;
 551     }
 552   }
 553
 554   /* calculate backwards path clipping, non-cyclic. j <= clip0[i] iff
 555      clip1[j] <= i, for i,j=0..n. */
 556   j = 1;
 557   for (i=0; i<n; i++) {
 558     while (j <= clip0[i]) {
 559       clip1[j] = i;
 560       j++;
 561     }
 562   }
 563
 564   /* calculate seg0[j] = longest path from 0 with j segments */
 565   i = 0;
 566   for (j=0; i<n; j++) {
 567     seg0[j] = i;
 568     i = clip0[i];
 569   }
 570   seg0[j] = n;
 571   m = j;
 572
 573   /* calculate seg1[j] = longest path to n with m-j segments */
 574   i = n;
 575   for (j=m; j>0; j--) {
 576     seg1[j] = i;
 577     i = clip1[i];
 578   }
 579   seg1[0] = 0;
 580
 581   /* now find the shortest path with m segments, based on penalty3 */
 582   /* note: the outer 2 loops jointly have at most n interations, thus
 583      the worst-case behavior here is quadratic. In practice, it is
 584      close to linear since the inner loop tends to be short. */
 585   pen[0]=0;
 586   for (j=1; j<=m; j++) {
 587     for (i=seg1[j]; i<=seg0[j]; i++) {
 588       best = -1;
 589       for (k=seg0[j-1]; k>=clip1[i]; k--) {
 590         thispen = penalty3(pp, k, i) + pen[k];
 591         if (best < 0 || thispen < best) {
 592           prev[i] = k;
 593           best = thispen;
 594         }
 595       }
 596       pen[i] = best;
 597     }
 598   }
 599
 600   pp->m = m;
 601   SAFE_MALLOC(pp->po, m, int);
 602
 603   /* read off shortest path */
 604   for (i=n, j=m-1; i>0; j--) {
 605     i = prev[i];
 606     pp->po[j] = i;
 607   }
 608
 609   free(pen);
 610   free(prev);
 611   free(clip0);
 612   free(clip1);
 613   free(seg0);
 614   free(seg1);
 615   return 0;
 616
 617  malloc_error:
 618   free(pen);
 619   free(prev);
 620   free(clip0);
 621   free(clip1);
 622   free(seg0);
 623   free(seg1);
 624   return 1;
 625 }
 626
 627 /* ---------------------------------------------------------------------- */
 628 /* Stage 3: vertex adjustment (Sec. 2.3.1). */
 629
 630 /* Adjust vertices of optimal polygon: calculate the intersection of
 631    the two "optimal" line segments, then move it into the unit square
 632    if it lies outside. Return 1 with errno set on error; 0 on
 633    success. */
 634
 635 static int adjust_vertices(privpath_t *pp) {
 636   int m = pp->m;
 637   int *po = pp->po;
 638   int n = pp->len;
 639   point_t *pt = pp->pt;
 640   int x0 = pp->x0;
 641   int y0 = pp->y0;
 642
 643   dpoint_t *ctr = NULL;      /* ctr[m] */
 644   dpoint_t *dir = NULL;      /* dir[m] */
 645   quadform_t *q = NULL;      /* q[m] */
 646   double v[3];
 647   double d;
 648   int i, j, k, l;
 649   dpoint_t s;
 650   int r;
 651
 652   SAFE_MALLOC(ctr, m, dpoint_t);
 653   SAFE_MALLOC(dir, m, dpoint_t);
 654   SAFE_MALLOC(q, m, quadform_t);
 655
 656   r = privcurve_init(&pp->curve, m);
 657   if (r) {
 658     goto malloc_error;
 659   }
 660
 661   /* calculate "optimal" point-slope representation for each line
 662      segment */
 663   for (i=0; i<m; i++) {
 664     j = po[mod(i+1,m)];
 665     j = mod(j-po[i],n)+po[i];
 666     pointslope(pp, po[i], j, &ctr[i], &dir[i]);
 667   }
 668
 669   /* represent each line segment as a singular quadratic form; the
 670      distance of a point (x,y) from the line segment will be
 671      (x,y,1)Q(x,y,1)^t, where Q=q[i]. */
 672   for (i=0; i<m; i++) {
 673     d = sq(dir[i].x) + sq(dir[i].y);
 674     if (d == 0.0) {
 675       for (j=0; j<3; j++) {
 676         for (k=0; k<3; k++) {
 677           q[i][j][k] = 0;
 678         }
 679       }
 680     } else {
 681       v[0] = dir[i].y;
 682       v[1] = -dir[i].x;
 683       v[2] = - v[1] * ctr[i].y - v[0] * ctr[i].x;
 684       for (l=0; l<3; l++) {
 685         for (k=0; k<3; k++) {
 686           q[i][l][k] = v[l] * v[k] / d;
 687         }
 688       }
 689     }
 690   }
 691
 692   /* now calculate the "intersections" of consecutive segments.
 693      Instead of using the actual intersection, we find the point
 694      within a given unit square which minimizes the square distance to
 695      the two lines. */
 696   for (i=0; i<m; i++) {
 697     quadform_t Q;
 698     dpoint_t w;
 699     double dx, dy;
 700     double det;
 701     double min, cand; /* minimum and candidate for minimum of quad. form */
 702     double xmin, ymin;  /* coordinates of minimum */
 703     int z;
 704
 705     /* let s be the vertex, in coordinates relative to x0/y0 */
 706     s.x = pt[po[i]].x-x0;
 707     s.y = pt[po[i]].y-y0;
 708
 709     /* intersect segments i-1 and i */
 710
 711     j = mod(i-1,m);
 712
 713     /* add quadratic forms */
 714     for (l=0; l<3; l++) {
 715       for (k=0; k<3; k++) {
 716         Q[l][k] = q[j][l][k] + q[i][l][k];
 717       }
 718     }
 719
 720     while(1) {
 721       /* minimize the quadratic form Q on the unit square */
 722       /* find intersection */
 723
 724 #ifdef HAVE_GCC_LOOP_BUG
 725       /* work around gcc bug #12243 */
 726       free(NULL);
 727 #endif
 728
 729       det = Q[0][0]*Q[1][1] - Q[0][1]*Q[1][0];
 730       if (det != 0.0) {
 731         w.x = (-Q[0][2]*Q[1][1] + Q[1][2]*Q[0][1]) / det;
 732         w.y = ( Q[0][2]*Q[1][0] - Q[1][2]*Q[0][0]) / det;
 733         break;
 734       }
 735
 736       /* matrix is singular - lines are parallel. Add another,
 737          orthogonal axis, through the center of the unit square */
 738       if (Q[0][0]>Q[1][1]) {
 739         v[0] = -Q[0][1];
 740         v[1] = Q[0][0];
 741       } else if (Q[1][1]) {
 742         v[0] = -Q[1][1];
 743         v[1] = Q[1][0];
 744       } else {
 745         v[0] = 1;
 746         v[1] = 0;
 747       }
 748       d = sq(v[0]) + sq(v[1]);
 749       v[2] = - v[1] * s.y - v[0] * s.x;
 750       for (l=0; l<3; l++) {
 751         for (k=0; k<3; k++) {
 752           Q[l][k] += v[l] * v[k] / d;
 753         }
 754       }
 755     }
 756     dx = fabs(w.x-s.x);
 757     dy = fabs(w.y-s.y);
 758     if (dx <= .5 && dy <= .5) {
 759       pp->curve.vertex[i].x = w.x+x0;
 760       pp->curve.vertex[i].y = w.y+y0;
 761       continue;
 762     }
 763
 764     /* the minimum was not in the unit square; now minimize quadratic
 765        on boundary of square */
 766     min = quadform(Q, s);
 767     xmin = s.x;
 768     ymin = s.y;
 769
 770     if (Q[0][0] == 0.0) {
 771       goto fixx;
 772     }
 773     for (z=0; z<2; z++) {   /* value of the y-coordinate */
 774       w.y = s.y-0.5+z;
 775       w.x = - (Q[0][1] * w.y + Q[0][2]) / Q[0][0];
 776       dx = fabs(w.x-s.x);
 777       cand = quadform(Q, w);
 778       if (dx <= .5 && cand < min) {
 779         min = cand;
 780         xmin = w.x;
 781         ymin = w.y;
 782       }
 783     }
 784   fixx:
 785     if (Q[1][1] == 0.0) {
 786       goto corners;
 787     }
 788     for (z=0; z<2; z++) {   /* value of the x-coordinate */
 789       w.x = s.x-0.5+z;
 790       w.y = - (Q[1][0] * w.x + Q[1][2]) / Q[1][1];
 791       dy = fabs(w.y-s.y);
 792       cand = quadform(Q, w);
 793       if (dy <= .5 && cand < min) {
 794         min = cand;
 795         xmin = w.x;
 796         ymin = w.y;
 797       }
 798     }
 799   corners:
 800     /* check four corners */
 801     for (l=0; l<2; l++) {
 802       for (k=0; k<2; k++) {
 803         w.x = s.x-0.5+l;
 804         w.y = s.y-0.5+k;
 805         cand = quadform(Q, w);
 806         if (cand < min) {
 807           min = cand;
 808           xmin = w.x;
 809           ymin = w.y;
 810         }
 811       }
 812     }
 813
 814     pp->curve.vertex[i].x = xmin + x0;
 815     pp->curve.vertex[i].y = ymin + y0;
 816     continue;
 817   }
 818
 819   free(ctr);
 820   free(dir);
 821   free(q);
 822   return 0;
 823
 824  malloc_error:
 825   free(ctr);
 826   free(dir);
 827   free(q);
 828   return 1;
 829 }
 830
 831 /* ---------------------------------------------------------------------- */
 832 /* Stage 4: smoothing and corner analysis (Sec. 2.3.3) */
 833
 834 /* Always succeeds and returns 0 */
 835 static int smooth(privcurve_t *curve, int sign, double alphamax) {
 836   int m = curve->n;
 837
 838   int i, j, k;
 839   double dd, denom, alpha;
 840   dpoint_t p2, p3, p4;
 841
 842   if (sign == '-') {
 843     /* reverse orientation of negative paths */
 844     for (i=0, j=m-1; i<j; i++, j--) {
 845       dpoint_t tmp;
 846       tmp = curve->vertex[i];
 847       curve->vertex[i] = curve->vertex[j];
 848       curve->vertex[j] = tmp;
 849     }
 850   }
 851
 852   /* examine each vertex and find its best fit */
 853   for (i=0; i<m; i++) {
 854     j = mod(i+1, m);
 855     k = mod(i+2, m);
 856     p4 = interval(1/2.0, curve->vertex[k], curve->vertex[j]);
 857
 858     denom = ddenom(curve->vertex[i], curve->vertex[k]);
 859     if (denom != 0.0) {
 860       dd = dpara(curve->vertex[i], curve->vertex[j], curve->vertex[k]) / denom;
 861       dd = fabs(dd);
 862       alpha = dd>1 ? (1 - 1.0/dd) : 0;
 863       alpha = alpha / 0.75;
 864     } else {
 865       alpha = 4/3.0;
 866     }
 867     curve->alpha0[j] = alpha;    /* remember "original" value of alpha */
 868
 869     if (alpha > alphamax) {  /* pointed corner */
 870       curve->tag[j] = POTRACE_CORNER;
 871       curve->c[j][1] = curve->vertex[j];
 872       curve->c[j][2] = p4;
 873     } else {
 874       if (alpha < 0.55) {
 875         alpha = 0.55;
 876       } else if (alpha > 1) {
 877         alpha = 1;
 878       }
 879       p2 = interval(.5+.5*alpha, curve->vertex[i], curve->vertex[j]);
 880       p3 = interval(.5+.5*alpha, curve->vertex[k], curve->vertex[j]);
 881       curve->tag[j] = POTRACE_CURVETO;
 882       curve->c[j][0] = p2;
 883       curve->c[j][1] = p3;
 884       curve->c[j][2] = p4;
 885     }
 886     curve->alpha[j] = alpha;    /* store the "cropped" value of alpha */
 887     curve->beta[j] = 0.5;
 888   }
 889   curve->alphacurve = 1;
 890   return 0;
 891 }
 892
 893 /* ---------------------------------------------------------------------- */
 894 /* Stage 5: Curve optimization (Sec. 2.4) */
 895
 896 /* a private type for the result of opti_penalty */
 897 struct opti_s {
 898   double pen;      /* penalty */
 899   dpoint_t c[2];   /* curve parameters */
 900   double t, s;     /* curve parameters */
 901   double alpha;    /* curve parameter */
 902 };
 903 typedef struct opti_s opti_t;
 904
 905 /* calculate best fit from i+.5 to j+.5.  Assume i<j (cyclically).
 906    Return 0 and set badness and parameters (alpha, beta), if
 907    possible. Return 1 if impossible. */
 908 static int opti_penalty(privpath_t *pp, int i, int j, opti_t *res, double opttolerance, int *convc, double *areac) {
 909   int m = pp->curve.n;
 910   int k, k1, k2, conv, i1;
 911   double area, alpha, d, d1, d2;
 912   dpoint_t p0, p1, p2, p3, pt;
 913   double A, R, A1, A2, A3, A4;
 914   double s, t;
 915
 916   /* check convexity, corner-freeness, and maximum bend < 179 degrees */
 917
 918   if (i==j) {  /* sanity - a full loop can never be an opticurve */
 919     return 1;
 920   }
 921
 922   k = i;
 923   i1 = mod(i+1, m);
 924   k1 = mod(k+1, m);
 925   conv = convc[k1];
 926   if (conv == 0) {
 927     return 1;
 928   }
 929   d = ddist(pp->curve.vertex[i], pp->curve.vertex[i1]);
 930   for (k=k1; k!=j; k=k1) {
 931     k1 = mod(k+1, m);
 932     k2 = mod(k+2, m);
 933     if (convc[k1] != conv) {
 934       return 1;
 935     }
 936     if (sign(cprod(pp->curve.vertex[i], pp->curve.vertex[i1], pp->curve.vertex[k1], pp->curve.vertex[k2])) != conv) {
 937       return 1;
 938     }
 939     if (iprod1(pp->curve.vertex[i], pp->curve.vertex[i1], pp->curve.vertex[k1], pp->curve.vertex[k2]) < d * ddist(pp->curve.vertex[k1], pp->curve.vertex[k2]) * COS179) {
 940       return 1;
 941     }
 942   }
 943
 944   /* the curve we're working in: */
 945   p0 = pp->curve.c[mod(i,m)][2];
 946   p1 = pp->curve.vertex[mod(i+1,m)];
 947   p2 = pp->curve.vertex[mod(j,m)];
 948   p3 = pp->curve.c[mod(j,m)][2];
 949
 950   /* determine its area */
 951   area = areac[j] - areac[i];
 952   area -= dpara(pp->curve.vertex[0], pp->curve.c[i][2], pp->curve.c[j][2])/2;
 953   if (i>=j) {
 954     area += areac[m];
 955   }
 956
 957   /* find intersection o of p0p1 and p2p3. Let t,s such that o =
 958      interval(t,p0,p1) = interval(s,p3,p2). Let A be the area of the
 959      triangle (p0,o,p3). */
 960
 961   A1 = dpara(p0, p1, p2);
 962   A2 = dpara(p0, p1, p3);
 963   A3 = dpara(p0, p2, p3);
 964   /* A4 = dpara(p1, p2, p3); */
 965   A4 = A1+A3-A2;
 966
 967   if (A2 == A1) {  /* this should never happen */
 968     return 1;
 969   }
 970
 971   t = A3/(A3-A4);
 972   s = A2/(A2-A1);
 973   A = A2 * t / 2.0;
 974
 975   if (A == 0.0) {  /* this should never happen */
 976     return 1;
 977   }
 978
 979   R = area / A;  /* relative area */
 980   alpha = 2 - sqrt(4 - R / 0.3);  /* overall alpha for p0-o-p3 curve */
 981
 982   res->c[0] = interval(t * alpha, p0, p1);
 983   res->c[1] = interval(s * alpha, p3, p2);
 984   res->alpha = alpha;
 985   res->t = t;
 986   res->s = s;
 987
 988   p1 = res->c[0];
 989   p2 = res->c[1];  /* the proposed curve is now (p0,p1,p2,p3) */
 990
 991   res->pen = 0;
 992
 993   /* calculate penalty */
 994   /* check tangency with edges */
 995   for (k=mod(i+1,m); k!=j; k=k1) {
 996     k1 = mod(k+1,m);
 997     t = tangent(p0, p1, p2, p3, pp->curve.vertex[k], pp->curve.vertex[k1]);
 998     if (t<-.5) {
 999       return 1;
1000     }
1001     pt = bezier(t, p0, p1, p2, p3);
1002     d = ddist(pp->curve.vertex[k], pp->curve.vertex[k1]);
1003     if (d == 0.0) {  /* this should never happen */
1004       return 1;
1005     }
1006     d1 = dpara(pp->curve.vertex[k], pp->curve.vertex[k1], pt) / d;
1007     if (fabs(d1) > opttolerance) {
1008       return 1;
1009     }
1010     if (iprod(pp->curve.vertex[k], pp->curve.vertex[k1], pt) < 0 || iprod(pp->curve.vertex[k1], pp->curve.vertex[k], pt) < 0) {
1011       return 1;
1012     }
1013     res->pen += sq(d1);
1014   }
1015
1016   /* check corners */
1017   for (k=i; k!=j; k=k1) {
1018     k1 = mod(k+1,m);
1019     t = tangent(p0, p1, p2, p3, pp->curve.c[k][2], pp->curve.c[k1][2]);
1020     if (t<-.5) {
1021       return 1;
1022     }
1023     pt = bezier(t, p0, p1, p2, p3);
1024     d = ddist(pp->curve.c[k][2], pp->curve.c[k1][2]);
1025     if (d == 0.0) {  /* this should never happen */
1026       return 1;
1027     }
1028     d1 = dpara(pp->curve.c[k][2], pp->curve.c[k1][2], pt) / d;
1029     d2 = dpara(pp->curve.c[k][2], pp->curve.c[k1][2], pp->curve.vertex[k1]) / d;
1030     d2 *= 0.75 * pp->curve.alpha[k1];
1031     if (d2 < 0) {
1032       d1 = -d1;
1033       d2 = -d2;
1034     }
1035     if (d1 < d2 - opttolerance) {
1036       return 1;
1037     }
1038     if (d1 < d2) {
1039       res->pen += sq(d1 - d2);
1040     }
1041   }
1042
1043   return 0;
1044 }
1045
1046 /* optimize the path p, replacing sequences of Bezier segments by a
1047    single segment when possible. Return 0 on success, 1 with errno set
1048    on failure. */
1049 static int opticurve(privpath_t *pp, double opttolerance) {
1050   int m = pp->curve.n;
1051   int *pt = NULL;     /* pt[m+1] */
1052   double *pen = NULL; /* pen[m+1] */
1053   int *len = NULL;    /* len[m+1] */
1054   opti_t *opt = NULL; /* opt[m+1] */
1055   int om;
1056   int i,j,r;
1057   opti_t o;
1058   dpoint_t p0;
1059   int i1;
1060   double area;
1061   double alpha;
1062   double *s = NULL;
1063   double *t = NULL;
1064
1065   int *convc = NULL; /* conv[m]: pre-computed convexities */
1066   double *areac = NULL; /* cumarea[m+1]: cache for fast area computation */
1067
1068   SAFE_MALLOC(pt, m+1, int);
1069   SAFE_MALLOC(pen, m+1, double);
1070   SAFE_MALLOC(len, m+1, int);
1071   SAFE_MALLOC(opt, m+1, opti_t);
1072   SAFE_MALLOC(convc, m, int);
1073   SAFE_MALLOC(areac, m+1, double);
1074
1075   /* pre-calculate convexity: +1 = right turn, -1 = left turn, 0 = corner */
1076   for (i=0; i<m; i++) {
1077     if (pp->curve.tag[i] == POTRACE_CURVETO) {
1078       convc[i] = sign(dpara(pp->curve.vertex[mod(i-1,m)], pp->curve.vertex[i], pp->curve.vertex[mod(i+1,m)]));
1079     } else {
1080       convc[i] = 0;
1081     }
1082   }
1083
1084   /* pre-calculate areas */
1085   area = 0.0;
1086   areac[0] = 0.0;
1087   p0 = pp->curve.vertex[0];
1088   for (i=0; i<m; i++) {
1089     i1 = mod(i+1, m);
1090     if (pp->curve.tag[i1] == POTRACE_CURVETO) {
1091       alpha = pp->curve.alpha[i1];
1092       area += 0.3*alpha*(4-alpha)*dpara(pp->curve.c[i][2], pp->curve.vertex[i1], pp->curve.c[i1][2])/2;
1093       area += dpara(p0, pp->curve.c[i][2], pp->curve.c[i1][2])/2;
1094     }
1095     areac[i+1] = area;
1096   }
1097
1098   pt[0] = -1;
1099   pen[0] = 0;
1100   len[0] = 0;
1101
1102   /* Fixme: we always start from a fixed point -- should find the best
1103      curve cyclically ### */
1104
1105   for (j=1; j<=m; j++) {
1106     /* calculate best path from 0 to j */
1107     pt[j] = j-1;
1108     pen[j] = pen[j-1];
1109     len[j] = len[j-1]+1;
1110
1111     for (i=j-2; i>=0; i--) {
1112       r = opti_penalty(pp, i, mod(j,m), &o, opttolerance, convc, areac);
1113       if (r) {
1114         break;
1115       }
1116       if (len[j] > len[i]+1 || (len[j] == len[i]+1 && pen[j] > pen[i] + o.pen)) {
1117         pt[j] = i;
1118         pen[j] = pen[i] + o.pen;
1119         len[j] = len[i] + 1;
1120         opt[j] = o;
1121       }
1122     }
1123   }
1124   om = len[m];
1125   r = privcurve_init(&pp->ocurve, om);
1126   if (r) {
1127     goto malloc_error;
1128   }
1129   SAFE_MALLOC(s, om, double);
1130   SAFE_MALLOC(t, om, double);
1131
1132   j = m;
1133   for (i=om-1; i>=0; i--) {
1134     if (pt[j]==j-1) {
1135       pp->ocurve.tag[i]    = pp->curve.tag[mod(j,m)];
1136       pp->ocurve.c[i][0]   = pp->curve.c[mod(j,m)][0];
1137       pp->ocurve.c[i][1]   = pp->curve.c[mod(j,m)][1];
1138       pp->ocurve.c[i][2]   = pp->curve.c[mod(j,m)][2];
1139       pp->ocurve.vertex[i] = pp->curve.vertex[mod(j,m)];
1140       pp->ocurve.alpha[i]  = pp->curve.alpha[mod(j,m)];
1141       pp->ocurve.alpha0[i] = pp->curve.alpha0[mod(j,m)];
1142       pp->ocurve.beta[i]   = pp->curve.beta[mod(j,m)];
1143       s[i] = t[i] = 1.0;
1144     } else {
1145       pp->ocurve.tag[i] = POTRACE_CURVETO;
1146       pp->ocurve.c[i][0] = opt[j].c[0];
1147       pp->ocurve.c[i][1] = opt[j].c[1];
1148       pp->ocurve.c[i][2] = pp->curve.c[mod(j,m)][2];
1149       pp->ocurve.vertex[i] = interval(opt[j].s, pp->curve.c[mod(j,m)][2], pp->curve.vertex[mod(j,m)]);
1150       pp->ocurve.alpha[i] = opt[j].alpha;
1151       pp->ocurve.alpha0[i] = opt[j].alpha;
1152       s[i] = opt[j].s;
1153       t[i] = opt[j].t;
1154     }
1155     j = pt[j];
1156   }
1157
1158   /* calculate beta parameters */
1159   for (i=0; i<om; i++) {
1160     i1 = mod(i+1,om);
1161     pp->ocurve.beta[i] = s[i] / (s[i] + t[i1]);
1162   }
1163   pp->ocurve.alphacurve = 1;
1164
1165   free(pt);
1166   free(pen);
1167   free(len);
1168   free(opt);
1169   free(s);
1170   free(t);
1171   free(convc);
1172   free(areac);
1173   return 0;
1174
1175  malloc_error:
1176   free(pt);
1177   free(pen);
1178   free(len);
1179   free(opt);
1180   free(s);
1181   free(t);
1182   free(convc);
1183   free(areac);
1184   return 1;
1185 }
1186
1187 /* ---------------------------------------------------------------------- */
1188
1189 #define TRY(x) if (x) goto try_error
1190
1191 /* return 0 on success, 1 on error with errno set. */
1192 int process_path(path_t *plist, const potrace_param_t *param, progress_t *progress) {
1193   path_t *p;
1194   double nn = 0, cn = 0;
1195
1196   if (progress->callback) {
1197     /* precompute task size for progress estimates */
1198     nn = 0;
1199     list_forall (p, plist) {
1200       nn += p->priv->len;
1201     }
1202     cn = 0;
1203   }
1204
1205   /* call downstream function with each path */
1206   list_forall (p, plist) {
1207     TRY(calc_sums(p->priv));
1208     TRY(calc_lon(p->priv));
1209     TRY(bestpolygon(p->priv));
1210     TRY(adjust_vertices(p->priv));
1211     TRY(smooth(&p->priv->curve, p->sign, param->alphamax));
1212     if (param->opticurve) {
1213       TRY(opticurve(p->priv, param->opttolerance));
1214       p->priv->fcurve = &p->priv->ocurve;
1215     } else {
1216       p->priv->fcurve = &p->priv->curve;
1217     }
1218     privcurve_to_curve(p->priv->fcurve, &p->curve);
1219
1220     if (progress->callback) {
1221       cn += p->priv->len;
1222       progress_update(cn/nn, progress);
1223     }
1224   }
1225
1226   progress_update(1.0, progress);
1227
1228   return 0;
1229
1230  try_error:
1231   return 1;
1232 }