src/trace/potrace/trace.cpp

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