trace.c
Go to the documentation of this file.
00001 /* Copyright (C) 2001-2007 Peter Selinger.
00002    This file is part of Potrace. It is free software and it is covered
00003    by the GNU General Public License. See the file COPYING for details. */
00004 
00005 /* $Id: trace.c 147 2007-04-09 00:44:09Z selinger $ */
00006 /* transform jaggy paths into smooth curves */
00007 
00008 #include <stdio.h>
00009 #include <math.h>
00010 #include <stdlib.h>
00011 #include <string.h>
00012 
00013 #include "potracelib.h"
00014 #include "curve.h"
00015 #include "lists.h"
00016 #include "auxiliary.h"
00017 #include "trace.h"
00018 #include "progress.h"
00019 
00020 #define INFTY 10000000  /* it suffices that this is longer than any
00021                            path; it need not be really infinite */
00022 #define COS179 -0.999847695156   /* the cosine of 179 degrees */
00023 
00024 /* ---------------------------------------------------------------------- */
00025 #define SAFE_MALLOC(var, n, typ) \
00026   if ((var = (typ *)malloc((n)*sizeof(typ))) == NULL) goto malloc_error 
00027 
00028 /* ---------------------------------------------------------------------- */
00029 /* auxiliary functions */
00030 
00031 /* return a direction that is 90 degrees counterclockwise from p2-p0,
00032    but then restricted to one of the major wind directions (n, nw, w, etc) */
00033 static inline point_t dorth_infty(dpoint_t p0, dpoint_t p2) {
00034   point_t r;
00035   
00036   r.y = sign(p2.x-p0.x);
00037   r.x = -sign(p2.y-p0.y);
00038 
00039   return r;
00040 }
00041 
00042 /* return (p1-p0)x(p2-p0), the area of the parallelogram */
00043 static inline double dpara(dpoint_t p0, dpoint_t p1, dpoint_t p2) {
00044   double x1, y1, x2, y2;
00045 
00046   x1 = p1.x-p0.x;
00047   y1 = p1.y-p0.y;
00048   x2 = p2.x-p0.x;
00049   y2 = p2.y-p0.y;
00050 
00051   return x1*y2 - x2*y1;
00052 }
00053 
00054 /* ddenom/dpara have the property that the square of radius 1 centered
00055    at p1 intersects the line p0p2 iff |dpara(p0,p1,p2)| <= ddenom(p0,p2) */
00056 static inline double ddenom(dpoint_t p0, dpoint_t p2) {
00057   point_t r = dorth_infty(p0, p2);
00058 
00059   return r.y*(p2.x-p0.x) - r.x*(p2.y-p0.y);
00060 }
00061 
00062 /* return 1 if a <= b < c < a, in a cyclic sense (mod n) */
00063 static inline int cyclic(int a, int b, int c) {
00064   if (a<=c) {
00065     return (a<=b && b<c);
00066   } else {
00067     return (a<=b || b<c);
00068   }
00069 }
00070 
00071 /* determine the center and slope of the line i..j. Assume i<j. Needs
00072    "sum" components of p to be set. */
00073 static void pointslope(privpath_t *pp, int i, int j, dpoint_t *ctr, dpoint_t *dir) {
00074   /* assume i<j */
00075 
00076   int n = pp->len;
00077   sums_t *sums = pp->sums;
00078 
00079   double x, y, x2, xy, y2;
00080   double k;
00081   double a, b, c, lambda2, l;
00082   int r=0; /* rotations from i to j */
00083 
00084   while (j>=n) {
00085     j-=n;
00086     r+=1;
00087   }
00088   while (i>=n) {
00089     i-=n;
00090     r-=1;
00091   }
00092   while (j<0) {
00093     j+=n;
00094     r-=1;
00095   }
00096   while (i<0) {
00097     i+=n;
00098     r+=1;
00099   }
00100   
00101   x = sums[j+1].x-sums[i].x+r*sums[n].x;
00102   y = sums[j+1].y-sums[i].y+r*sums[n].y;
00103   x2 = sums[j+1].x2-sums[i].x2+r*sums[n].x2;
00104   xy = sums[j+1].xy-sums[i].xy+r*sums[n].xy;
00105   y2 = sums[j+1].y2-sums[i].y2+r*sums[n].y2;
00106   k = j+1-i+r*n;
00107   
00108   ctr->x = x/k;
00109   ctr->y = y/k;
00110 
00111   a = (x2-(double)x*x/k)/k;
00112   b = (xy-(double)x*y/k)/k;
00113   c = (y2-(double)y*y/k)/k;
00114   
00115   lambda2 = (a+c+sqrt((a-c)*(a-c)+4*b*b))/2; /* larger e.value */
00116 
00117   /* now find e.vector for lambda2 */
00118   a -= lambda2;
00119   c -= lambda2;
00120 
00121   if (fabs(a) >= fabs(c)) {
00122     l = sqrt(a*a+b*b);
00123     if (l!=0) {
00124       dir->x = -b/l;
00125       dir->y = a/l;
00126     }
00127   } else {
00128     l = sqrt(c*c+b*b);
00129     if (l!=0) {
00130       dir->x = -c/l;
00131       dir->y = b/l;
00132     }
00133   }
00134   if (l==0) {
00135     dir->x = dir->y = 0;   /* sometimes this can happen when k=4:
00136                               the two eigenvalues coincide */
00137   }
00138 }
00139 
00140 /* the type of (affine) quadratic forms, represented as symmetric 3x3
00141    matrices.  The value of the quadratic form at a vector (x,y) is v^t
00142    Q v, where v = (x,y,1)^t. */
00143 typedef double quadform_t[3][3];
00144 
00145 /* Apply quadratic form Q to vector w = (w.x,w.y) */
00146 static inline double quadform(quadform_t Q, dpoint_t w) {
00147   double v[3];
00148   int i, j;
00149   double sum;
00150 
00151   v[0] = w.x;
00152   v[1] = w.y;
00153   v[2] = 1;
00154   sum = 0.0;
00155 
00156   for (i=0; i<3; i++) {
00157     for (j=0; j<3; j++) {
00158       sum += v[i] * Q[i][j] * v[j];
00159     }
00160   }
00161   return sum;
00162 }
00163 
00164 /* calculate p1 x p2 */
00165 static inline int xprod(point_t p1, point_t p2) {
00166   return p1.x*p2.y - p1.y*p2.x;
00167 }
00168 
00169 /* calculate (p1-p0)x(p3-p2) */
00170 static inline double cprod(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3) {
00171   double x1, y1, x2, y2;
00172 
00173   x1 = p1.x - p0.x;
00174   y1 = p1.y - p0.y;
00175   x2 = p3.x - p2.x;
00176   y2 = p3.y - p2.y;
00177 
00178   return x1*y2 - x2*y1;
00179 }
00180 
00181 /* calculate (p1-p0)*(p2-p0) */
00182 static inline double iprod(dpoint_t p0, dpoint_t p1, dpoint_t p2) {
00183   double x1, y1, x2, y2;
00184 
00185   x1 = p1.x - p0.x;
00186   y1 = p1.y - p0.y;
00187   x2 = p2.x - p0.x;
00188   y2 = p2.y - p0.y;
00189 
00190   return x1*x2 + y1*y2;
00191 }
00192 
00193 /* calculate (p1-p0)*(p3-p2) */
00194 static inline double iprod1(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3) {
00195   double x1, y1, x2, y2;
00196 
00197   x1 = p1.x - p0.x;
00198   y1 = p1.y - p0.y;
00199   x2 = p3.x - p2.x;
00200   y2 = p3.y - p2.y;
00201 
00202   return x1*x2 + y1*y2;
00203 }
00204 
00205 /* calculate distance between two points */
00206 static inline double ddist(dpoint_t p, dpoint_t q) {
00207   return sqrt(sq(p.x-q.x)+sq(p.y-q.y));
00208 }
00209 
00210 /* calculate point of a bezier curve */
00211 static inline dpoint_t bezier(double t, dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3) {
00212   double s = 1-t;
00213   dpoint_t res;
00214 
00215   /* Note: a good optimizing compiler (such as gcc-3) reduces the
00216      following to 16 multiplications, using common subexpression
00217      elimination. */
00218 
00219   res.x = s*s*s*p0.x + 3*(s*s*t)*p1.x + 3*(t*t*s)*p2.x + t*t*t*p3.x;
00220   res.y = s*s*s*p0.y + 3*(s*s*t)*p1.y + 3*(t*t*s)*p2.y + t*t*t*p3.y;
00221 
00222   return res;
00223 }
00224 
00225 /* calculate the point t in [0..1] on the (convex) bezier curve
00226    (p0,p1,p2,p3) which is tangent to q1-q0. Return -1.0 if there is no
00227    solution in [0..1]. */
00228 static double tangent(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3, dpoint_t q0, dpoint_t q1) {
00229   double A, B, C;   /* (1-t)^2 A + 2(1-t)t B + t^2 C = 0 */
00230   double a, b, c;   /* a t^2 + b t + c = 0 */
00231   double d, s, r1, r2;
00232 
00233   A = cprod(p0, p1, q0, q1);
00234   B = cprod(p1, p2, q0, q1);
00235   C = cprod(p2, p3, q0, q1);
00236 
00237   a = A - 2*B + C;
00238   b = -2*A + 2*B;
00239   c = A;
00240   
00241   d = b*b - 4*a*c;
00242 
00243   if (a==0 || d<0) {
00244     return -1.0;
00245   }
00246 
00247   s = sqrt(d);
00248 
00249   r1 = (-b + s) / (2 * a);
00250   r2 = (-b - s) / (2 * a);
00251 
00252   if (r1 >= 0 && r1 <= 1) {
00253     return r1;
00254   } else if (r2 >= 0 && r2 <= 1) {
00255     return r2;
00256   } else {
00257     return -1.0;
00258   }
00259 }
00260 
00261 /* ---------------------------------------------------------------------- */
00262 /* Preparation: fill in the sum* fields of a path (used for later
00263    rapid summing). Return 0 on success, 1 with errno set on
00264    failure. */
00265 static int calc_sums(privpath_t *pp) {
00266   int i, x, y;
00267   int n = pp->len;
00268 
00269   SAFE_MALLOC(pp->sums, pp->len+1, sums_t);
00270 
00271   /* origin */
00272   pp->x0 = pp->pt[0].x;
00273   pp->y0 = pp->pt[0].y;
00274 
00275   /* preparatory computation for later fast summing */
00276   pp->sums[0].x2 = pp->sums[0].xy = pp->sums[0].y2 = pp->sums[0].x = pp->sums[0].y = 0;
00277   for (i=0; i<n; i++) {
00278     x = pp->pt[i].x - pp->x0;
00279     y = pp->pt[i].y - pp->y0;
00280     pp->sums[i+1].x = pp->sums[i].x + x;
00281     pp->sums[i+1].y = pp->sums[i].y + y;
00282     pp->sums[i+1].x2 = pp->sums[i].x2 + x*x;
00283     pp->sums[i+1].xy = pp->sums[i].xy + x*y;
00284     pp->sums[i+1].y2 = pp->sums[i].y2 + y*y;
00285   }
00286   return 0;  
00287 
00288  malloc_error:
00289   return 1;
00290 }
00291 
00292 /* ---------------------------------------------------------------------- */
00293 /* Stage 1: determine the straight subpaths (Sec. 2.2.1). Fill in the
00294    "lon" component of a path object (based on pt/len).  For each i,
00295    lon[i] is the furthest index such that a straight line can be drawn
00296    from i to lon[i]. Return 1 on error with errno set, else 0. */
00297 
00298 /* this algorithm depends on the fact that the existence of straight
00299    subpaths is a triplewise property. I.e., there exists a straight
00300    line through squares i0,...,in iff there exists a straight line
00301    through i,j,k, for all i0<=i<j<k<=in. (Proof?) */
00302 
00303 /* this implementation of calc_lon is O(n^2). It replaces an older
00304    O(n^3) version. A "constraint" means that future points must
00305    satisfy xprod(constraint[0], cur) >= 0 and xprod(constraint[1],
00306    cur) <= 0. */
00307 
00308 /* Remark for Potrace 1.1: the current implementation of calc_lon is
00309    more complex than the implementation found in Potrace 1.0, but it
00310    is considerably faster. The introduction of the "nc" data structure
00311    means that we only have to test the constraints for "corner"
00312    points. On a typical input file, this speeds up the calc_lon
00313    function by a factor of 31.2, thereby decreasing its time share
00314    within the overall Potrace algorithm from 72.6% to 7.82%, and
00315    speeding up the overall algorithm by a factor of 3.36. On another
00316    input file, calc_lon was sped up by a factor of 6.7, decreasing its
00317    time share from 51.4% to 13.61%, and speeding up the overall
00318    algorithm by a factor of 1.78. In any case, the savings are
00319    substantial. */
00320 
00321 /* returns 0 on success, 1 on error with errno set */
00322 static int calc_lon(privpath_t *pp) {
00323   point_t *pt = pp->pt;
00324   int n = pp->len;
00325   int i, j, k, k1;
00326   int ct[4], dir;
00327   point_t constraint[2];
00328   point_t cur;
00329   point_t off;
00330   int *pivk = NULL;  /* pivk[n] */
00331   int *nc = NULL;    /* nc[n]: next corner */
00332   point_t dk;  /* direction of k-k1 */
00333   int a, b, c, d;
00334 
00335   SAFE_MALLOC(pivk, n, int);
00336   SAFE_MALLOC(nc, n, int);
00337 
00338   /* initialize the nc data structure. Point from each point to the
00339      furthest future point to which it is connected by a vertical or
00340      horizontal segment. We take advantage of the fact that there is
00341      always a direction change at 0 (due to the path decomposition
00342      algorithm). But even if this were not so, there is no harm, as
00343      in practice, correctness does not depend on the word "furthest"
00344      above.  */
00345   k = 0;
00346   for (i=n-1; i>=0; i--) {
00347     if (pt[i].x != pt[k].x && pt[i].y != pt[k].y) {
00348       k = i+1;  /* necessarily i<n-1 in this case */
00349     }
00350     nc[i] = k;
00351   }
00352 
00353   SAFE_MALLOC(pp->lon, n, int);
00354 
00355   /* determine pivot points: for each i, let pivk[i] be the furthest k
00356      such that all j with i<j<k lie on a line connecting i,k. */
00357   
00358   for (i=n-1; i>=0; i--) {
00359     ct[0] = ct[1] = ct[2] = ct[3] = 0;
00360 
00361     /* keep track of "directions" that have occurred */
00362     dir = (3+3*(pt[mod(i+1,n)].x-pt[i].x)+(pt[mod(i+1,n)].y-pt[i].y))/2;
00363     ct[dir]++;
00364 
00365     constraint[0].x = 0;
00366     constraint[0].y = 0;
00367     constraint[1].x = 0;
00368     constraint[1].y = 0;
00369 
00370     /* find the next k such that no straight line from i to k */
00371     k = nc[i];
00372     k1 = i;
00373     while (1) {
00374       
00375       dir = (3+3*sign(pt[k].x-pt[k1].x)+sign(pt[k].y-pt[k1].y))/2;
00376       ct[dir]++;
00377 
00378       /* if all four "directions" have occurred, cut this path */
00379       if (ct[0] && ct[1] && ct[2] && ct[3]) {
00380         pivk[i] = k1;
00381         goto foundk;
00382       }
00383 
00384       cur.x = pt[k].x - pt[i].x;
00385       cur.y = pt[k].y - pt[i].y;
00386 
00387       /* see if current constraint is violated */
00388       if (xprod(constraint[0], cur) < 0 || xprod(constraint[1], cur) > 0) {
00389         goto constraint_viol;
00390       }
00391 
00392       /* else, update constraint */
00393       if (abs(cur.x) <= 1 && abs(cur.y) <= 1) {
00394         /* no constraint */
00395       } else {
00396         off.x = cur.x + ((cur.y>=0 && (cur.y>0 || cur.x<0)) ? 1 : -1);
00397         off.y = cur.y + ((cur.x<=0 && (cur.x<0 || cur.y<0)) ? 1 : -1);
00398         if (xprod(constraint[0], off) >= 0) {
00399           constraint[0] = off;
00400         }
00401         off.x = cur.x + ((cur.y<=0 && (cur.y<0 || cur.x<0)) ? 1 : -1);
00402         off.y = cur.y + ((cur.x>=0 && (cur.x>0 || cur.y<0)) ? 1 : -1);
00403         if (xprod(constraint[1], off) <= 0) {
00404           constraint[1] = off;
00405         }
00406       } 
00407       k1 = k;
00408       k = nc[k1];
00409       if (!cyclic(k,i,k1)) {
00410         break;
00411       }
00412     }
00413   constraint_viol:
00414     /* k1 was the last "corner" satisfying the current constraint, and
00415        k is the first one violating it. We now need to find the last
00416        point along k1..k which satisfied the constraint. */
00417     dk.x = sign(pt[k].x-pt[k1].x);
00418     dk.y = sign(pt[k].y-pt[k1].y);
00419     cur.x = pt[k1].x - pt[i].x;
00420     cur.y = pt[k1].y - pt[i].y;
00421     /* find largest integer j such that xprod(constraint[0], cur+j*dk)
00422        >= 0 and xprod(constraint[1], cur+j*dk) <= 0. Use bilinearity
00423        of xprod. */
00424     a = xprod(constraint[0], cur);
00425     b = xprod(constraint[0], dk);
00426     c = xprod(constraint[1], cur);
00427     d = xprod(constraint[1], dk);
00428     /* find largest integer j such that a+j*b>=0 and c+j*d<=0. This
00429        can be solved with integer arithmetic. */
00430     j = INFTY;
00431     if (b<0) {
00432       j = floordiv(a,-b);
00433     }
00434     if (d>0) {
00435       j = min(j, floordiv(-c,d));
00436     }
00437     pivk[i] = mod(k1+j,n);
00438   foundk:
00439     ;
00440   } /* for i */
00441 
00442   /* clean up: for each i, let lon[i] be the largest k such that for
00443      all i' with i<=i'<k, i'<k<=pivk[i']. */
00444 
00445   j=pivk[n-1];
00446   pp->lon[n-1]=j;
00447   for (i=n-2; i>=0; i--) {
00448     if (cyclic(i+1,pivk[i],j)) {
00449       j=pivk[i];
00450     }
00451     pp->lon[i]=j;
00452   }
00453 
00454   for (i=n-1; cyclic(mod(i+1,n),j,pp->lon[i]); i--) {
00455     pp->lon[i] = j;
00456   }
00457 
00458   free(pivk);
00459   free(nc);
00460   return 0;
00461 
00462  malloc_error:
00463   free(pivk);
00464   free(nc);
00465   return 1;
00466 }
00467 
00468 
00469 /* ---------------------------------------------------------------------- */
00470 /* Stage 2: calculate the optimal polygon (Sec. 2.2.2-2.2.4). */ 
00471 
00472 /* Auxiliary function: calculate the penalty of an edge from i to j in
00473    the given path. This needs the "lon" and "sum*" data. */
00474 
00475 static double penalty3(privpath_t *pp, int i, int j) {
00476   int n = pp->len;
00477   point_t *pt = pp->pt;
00478   sums_t *sums = pp->sums;
00479 
00480   /* assume 0<=i<j<=n  */
00481   double x, y, x2, xy, y2;
00482   double k;
00483   double a, b, c, s;
00484   double px, py, ex, ey;
00485 
00486   int r=0; /* rotations from i to j */
00487 
00488   if (j>=n) {
00489     j-=n;
00490     r+=1;
00491   }
00492   
00493   x = sums[j+1].x-sums[i].x+r*sums[n].x;
00494   y = sums[j+1].y-sums[i].y+r*sums[n].y;
00495   x2 = sums[j+1].x2-sums[i].x2+r*sums[n].x2;
00496   xy = sums[j+1].xy-sums[i].xy+r*sums[n].xy;
00497   y2 = sums[j+1].y2-sums[i].y2+r*sums[n].y2;
00498   k = j+1-i+r*n;
00499   
00500   px = (pt[i].x+pt[j].x)/2.0-pt[0].x;
00501   py = (pt[i].y+pt[j].y)/2.0-pt[0].y;
00502   ey = (pt[j].x-pt[i].x);
00503   ex = -(pt[j].y-pt[i].y);
00504 
00505   a = ((x2-2*x*px)/k+px*px);
00506   b = ((xy-x*py-y*px)/k+px*py);
00507   c = ((y2-2*y*py)/k+py*py);
00508   
00509   s = ex*ex*a + 2*ex*ey*b + ey*ey*c;
00510 
00511   return sqrt(s);
00512 }
00513 
00514 /* find the optimal polygon. Fill in the m and po components. Return 1
00515    on failure with errno set, else 0. Non-cyclic version: assumes i=0
00516    is in the polygon. Fixme: ### implement cyclic version. */
00517 static int bestpolygon(privpath_t *pp)
00518 {
00519   int i, j, m, k;     
00520   int n = pp->len;
00521   double *pen = NULL; /* pen[n+1]: penalty vector */
00522   int *prev = NULL;   /* prev[n+1]: best path pointer vector */
00523   int *clip0 = NULL;  /* clip0[n]: longest segment pointer, non-cyclic */
00524   int *clip1 = NULL;  /* clip1[n+1]: backwards segment pointer, non-cyclic */
00525   int *seg0 = NULL;    /* seg0[m+1]: forward segment bounds, m<=n */
00526   int *seg1 = NULL;   /* seg1[m+1]: backward segment bounds, m<=n */
00527   double thispen;
00528   double best;
00529   int c;
00530 
00531   SAFE_MALLOC(pen, n+1, double);
00532   SAFE_MALLOC(prev, n+1, int);
00533   SAFE_MALLOC(clip0, n, int);
00534   SAFE_MALLOC(clip1, n+1, int);
00535   SAFE_MALLOC(seg0, n+1, int);
00536   SAFE_MALLOC(seg1, n+1, int);
00537 
00538   /* calculate clipped paths */
00539   for (i=0; i<n; i++) {
00540     c = mod(pp->lon[mod(i-1,n)]-1,n);
00541     if (c == i) {
00542       c = mod(i+1,n);
00543     }
00544     if (c < i) {
00545       clip0[i] = n;
00546     } else {
00547       clip0[i] = c;
00548     }
00549   }
00550 
00551   /* calculate backwards path clipping, non-cyclic. j <= clip0[i] iff
00552      clip1[j] <= i, for i,j=0..n. */
00553   j = 1;
00554   for (i=0; i<n; i++) {
00555     while (j <= clip0[i]) {
00556       clip1[j] = i;
00557       j++;
00558     }
00559   }
00560 
00561   /* calculate seg0[j] = longest path from 0 with j segments */
00562   i = 0;
00563   for (j=0; i<n; j++) {
00564     seg0[j] = i;
00565     i = clip0[i];
00566   }
00567   seg0[j] = n;
00568   m = j;
00569 
00570   /* calculate seg1[j] = longest path to n with m-j segments */
00571   i = n;
00572   for (j=m; j>0; j--) {
00573     seg1[j] = i;
00574     i = clip1[i];
00575   }
00576   seg1[0] = 0;
00577 
00578   /* now find the shortest path with m segments, based on penalty3 */
00579   /* note: the outer 2 loops jointly have at most n interations, thus
00580      the worst-case behavior here is quadratic. In practice, it is
00581      close to linear since the inner loop tends to be short. */
00582   pen[0]=0;
00583   for (j=1; j<=m; j++) {
00584     for (i=seg1[j]; i<=seg0[j]; i++) {
00585       best = -1;
00586       for (k=seg0[j-1]; k>=clip1[i]; k--) {
00587         thispen = penalty3(pp, k, i) + pen[k];
00588         if (best < 0 || thispen < best) {
00589           prev[i] = k;
00590           best = thispen;
00591         }
00592       }
00593       pen[i] = best;
00594     }
00595   }
00596 
00597   pp->m = m;
00598   SAFE_MALLOC(pp->po, m, int);
00599 
00600   /* read off shortest path */
00601   for (i=n, j=m-1; i>0; j--) {
00602     i = prev[i];
00603     pp->po[j] = i;
00604   }
00605 
00606   free(pen);
00607   free(prev);
00608   free(clip0);
00609   free(clip1);
00610   free(seg0);
00611   free(seg1);
00612   return 0;
00613   
00614  malloc_error:
00615   free(pen);
00616   free(prev);
00617   free(clip0);
00618   free(clip1);
00619   free(seg0);
00620   free(seg1);
00621   return 1;
00622 }
00623 
00624 /* ---------------------------------------------------------------------- */
00625 /* Stage 3: vertex adjustment (Sec. 2.3.1). */
00626 
00627 /* Adjust vertices of optimal polygon: calculate the intersection of
00628    the two "optimal" line segments, then move it into the unit square
00629    if it lies outside. Return 1 with errno set on error; 0 on
00630    success. */
00631 
00632 static int adjust_vertices(privpath_t *pp) {
00633   int m = pp->m;
00634   int *po = pp->po;
00635   int n = pp->len;
00636   point_t *pt = pp->pt;
00637   int x0 = pp->x0;
00638   int y0 = pp->y0;
00639 
00640   dpoint_t *ctr = NULL;      /* ctr[m] */
00641   dpoint_t *dir = NULL;      /* dir[m] */
00642   quadform_t *q = NULL;      /* q[m] */
00643   double v[3];
00644   double d;
00645   int i, j, k, l;
00646   dpoint_t s;
00647   int r;
00648 
00649   SAFE_MALLOC(ctr, m, dpoint_t);
00650   SAFE_MALLOC(dir, m, dpoint_t);
00651   SAFE_MALLOC(q, m, quadform_t);
00652 
00653   r = privcurve_init(&pp->curve, m);
00654   if (r) {
00655     goto malloc_error;
00656   }
00657   
00658   /* calculate "optimal" point-slope representation for each line
00659      segment */
00660   for (i=0; i<m; i++) {
00661     j = po[mod(i+1,m)];
00662     j = mod(j-po[i],n)+po[i];
00663     pointslope(pp, po[i], j, &ctr[i], &dir[i]);
00664   }
00665 
00666   /* represent each line segment as a singular quadratic form; the
00667      distance of a point (x,y) from the line segment will be
00668      (x,y,1)Q(x,y,1)^t, where Q=q[i]. */
00669   for (i=0; i<m; i++) {
00670     d = sq(dir[i].x) + sq(dir[i].y);
00671     if (d == 0.0) {
00672       for (j=0; j<3; j++) {
00673         for (k=0; k<3; k++) {
00674           q[i][j][k] = 0;
00675         }
00676       }
00677     } else {
00678       v[0] = dir[i].y;
00679       v[1] = -dir[i].x;
00680       v[2] = - v[1] * ctr[i].y - v[0] * ctr[i].x;
00681       for (l=0; l<3; l++) {
00682         for (k=0; k<3; k++) {
00683           q[i][l][k] = v[l] * v[k] / d;
00684         }
00685       }
00686     }
00687   }
00688 
00689   /* now calculate the "intersections" of consecutive segments.
00690      Instead of using the actual intersection, we find the point
00691      within a given unit square which minimizes the square distance to
00692      the two lines. */
00693   for (i=0; i<m; i++) {
00694     quadform_t Q;
00695     dpoint_t w;
00696     double dx, dy;
00697     double det;
00698     double min, cand; /* minimum and candidate for minimum of quad. form */
00699     double xmin, ymin;  /* coordinates of minimum */
00700     int z;
00701 
00702     /* let s be the vertex, in coordinates relative to x0/y0 */
00703     s.x = pt[po[i]].x-x0;
00704     s.y = pt[po[i]].y-y0;
00705 
00706     /* intersect segments i-1 and i */
00707 
00708     j = mod(i-1,m);
00709 
00710     /* add quadratic forms */
00711     for (l=0; l<3; l++) {
00712       for (k=0; k<3; k++) {
00713         Q[l][k] = q[j][l][k] + q[i][l][k];
00714       }
00715     }
00716     
00717     while(1) {
00718       /* minimize the quadratic form Q on the unit square */
00719       /* find intersection */
00720 
00721 #ifdef HAVE_GCC_LOOP_BUG
00722       /* work around gcc bug #12243 */
00723       free(NULL);
00724 #endif
00725       
00726       det = Q[0][0]*Q[1][1] - Q[0][1]*Q[1][0];
00727       if (det != 0.0) {
00728         w.x = (-Q[0][2]*Q[1][1] + Q[1][2]*Q[0][1]) / det;
00729         w.y = ( Q[0][2]*Q[1][0] - Q[1][2]*Q[0][0]) / det;
00730         break;
00731       }
00732 
00733       /* matrix is singular - lines are parallel. Add another,
00734          orthogonal axis, through the center of the unit square */
00735       if (Q[0][0]>Q[1][1]) {
00736         v[0] = -Q[0][1];
00737         v[1] = Q[0][0];
00738       } else if (Q[1][1]) {
00739         v[0] = -Q[1][1];
00740         v[1] = Q[1][0];
00741       } else {
00742         v[0] = 1;
00743         v[1] = 0;
00744       }
00745       d = sq(v[0]) + sq(v[1]);
00746       v[2] = - v[1] * s.y - v[0] * s.x;
00747       for (l=0; l<3; l++) {
00748         for (k=0; k<3; k++) {
00749           Q[l][k] += v[l] * v[k] / d;
00750         }
00751       }
00752     }
00753     dx = fabs(w.x-s.x);
00754     dy = fabs(w.y-s.y);
00755     if (dx <= .5 && dy <= .5) {
00756       pp->curve.vertex[i].x = w.x+x0;
00757       pp->curve.vertex[i].y = w.y+y0;
00758       continue;
00759     }
00760 
00761     /* the minimum was not in the unit square; now minimize quadratic
00762        on boundary of square */
00763     min = quadform(Q, s);
00764     xmin = s.x;
00765     ymin = s.y;
00766 
00767     if (Q[0][0] == 0.0) {
00768       goto fixx;
00769     }
00770     for (z=0; z<2; z++) {   /* value of the y-coordinate */
00771       w.y = s.y-0.5+z;
00772       w.x = - (Q[0][1] * w.y + Q[0][2]) / Q[0][0];
00773       dx = fabs(w.x-s.x);
00774       cand = quadform(Q, w);
00775       if (dx <= .5 && cand < min) {
00776         min = cand;
00777         xmin = w.x;
00778         ymin = w.y;
00779       }
00780     }
00781   fixx:
00782     if (Q[1][1] == 0.0) {
00783       goto corners;
00784     }
00785     for (z=0; z<2; z++) {   /* value of the x-coordinate */
00786       w.x = s.x-0.5+z;
00787       w.y = - (Q[1][0] * w.x + Q[1][2]) / Q[1][1];
00788       dy = fabs(w.y-s.y);
00789       cand = quadform(Q, w);
00790       if (dy <= .5 && cand < min) {
00791         min = cand;
00792         xmin = w.x;
00793         ymin = w.y;
00794       }
00795     }
00796   corners:
00797     /* check four corners */
00798     for (l=0; l<2; l++) {
00799       for (k=0; k<2; k++) {
00800         w.x = s.x-0.5+l;
00801         w.y = s.y-0.5+k;
00802         cand = quadform(Q, w);
00803         if (cand < min) {
00804           min = cand;
00805           xmin = w.x;
00806           ymin = w.y;
00807         }
00808       }
00809     }
00810 
00811     pp->curve.vertex[i].x = xmin + x0;
00812     pp->curve.vertex[i].y = ymin + y0;
00813     continue;
00814   }
00815 
00816   free(ctr);
00817   free(dir);
00818   free(q);
00819   return 0;
00820 
00821  malloc_error:
00822   free(ctr);
00823   free(dir);
00824   free(q);
00825   return 1;
00826 }
00827 
00828 /* ---------------------------------------------------------------------- */
00829 /* Stage 4: smoothing and corner analysis (Sec. 2.3.3) */
00830 
00831 /* Always succeeds and returns 0 */
00832 static int smooth(privcurve_t *curve, int sign, double alphamax) {
00833   int m = curve->n;
00834 
00835   int i, j, k;
00836   double dd, denom, alpha;
00837   dpoint_t p2, p3, p4;
00838 
00839   if (sign == '-') {
00840     /* reverse orientation of negative paths */
00841     for (i=0, j=m-1; i<j; i++, j--) {
00842       dpoint_t tmp;
00843       tmp = curve->vertex[i];
00844       curve->vertex[i] = curve->vertex[j];
00845       curve->vertex[j] = tmp;
00846     }
00847   }
00848 
00849   /* examine each vertex and find its best fit */
00850   for (i=0; i<m; i++) {
00851     j = mod(i+1, m);
00852     k = mod(i+2, m);
00853     p4 = interval(1/2.0, curve->vertex[k], curve->vertex[j]);
00854 
00855     denom = ddenom(curve->vertex[i], curve->vertex[k]);
00856     if (denom != 0.0) {
00857       dd = dpara(curve->vertex[i], curve->vertex[j], curve->vertex[k]) / denom;
00858       dd = fabs(dd);
00859       alpha = dd>1 ? (1 - 1.0/dd) : 0;
00860       alpha = alpha / 0.75;
00861     } else {
00862       alpha = 4/3.0;
00863     }
00864     curve->alpha0[j] = alpha;    /* remember "original" value of alpha */
00865 
00866     if (alpha > alphamax) {  /* pointed corner */
00867       curve->tag[j] = POTRACE_CORNER;
00868       curve->c[j][1] = curve->vertex[j];
00869       curve->c[j][2] = p4;
00870     } else {
00871       if (alpha < 0.55) {
00872         alpha = 0.55;
00873       } else if (alpha > 1) {
00874         alpha = 1;
00875       }
00876       p2 = interval(.5+.5*alpha, curve->vertex[i], curve->vertex[j]);
00877       p3 = interval(.5+.5*alpha, curve->vertex[k], curve->vertex[j]);
00878       curve->tag[j] = POTRACE_CURVETO;
00879       curve->c[j][0] = p2;
00880       curve->c[j][1] = p3;
00881       curve->c[j][2] = p4;
00882     }
00883     curve->alpha[j] = alpha;    /* store the "cropped" value of alpha */
00884     curve->beta[j] = 0.5;
00885   }
00886   curve->alphacurve = 1;
00887 
00888   return 0;
00889 }
00890 
00891 /* ---------------------------------------------------------------------- */
00892 /* Stage 5: Curve optimization (Sec. 2.4) */
00893 
00894 /* a private type for the result of opti_penalty */
00895 struct opti_s {
00896   double pen;      /* penalty */
00897   dpoint_t c[2];   /* curve parameters */
00898   double t, s;     /* curve parameters */
00899   double alpha;    /* curve parameter */
00900 };
00901 typedef struct opti_s opti_t;
00902 
00903 /* calculate best fit from i+.5 to j+.5.  Assume i<j (cyclically).
00904    Return 0 and set badness and parameters (alpha, beta), if
00905    possible. Return 1 if impossible. */
00906 static int opti_penalty(privpath_t *pp, int i, int j, opti_t *res, double opttolerance, int *convc, double *areac) {
00907   int m = pp->curve.n;
00908   int k, k1, k2, conv, i1;
00909   double area, alpha, d, d1, d2;
00910   dpoint_t p0, p1, p2, p3, pt;
00911   double A, R, A1, A2, A3, A4;
00912   double s, t;
00913 
00914   /* check convexity, corner-freeness, and maximum bend < 179 degrees */
00915 
00916   if (i==j) {  /* sanity - a full loop can never be an opticurve */
00917     return 1;
00918   }
00919 
00920   k = i;
00921   i1 = mod(i+1, m);
00922   k1 = mod(k+1, m);
00923   conv = convc[k1];
00924   if (conv == 0) {
00925     return 1;
00926   }
00927   d = ddist(pp->curve.vertex[i], pp->curve.vertex[i1]);
00928   for (k=k1; k!=j; k=k1) {
00929     k1 = mod(k+1, m);
00930     k2 = mod(k+2, m);
00931     if (convc[k1] != conv) {
00932       return 1;
00933     }
00934     if (sign(cprod(pp->curve.vertex[i], pp->curve.vertex[i1], pp->curve.vertex[k1], pp->curve.vertex[k2])) != conv) {
00935       return 1;
00936     }
00937     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) {
00938       return 1;
00939     }
00940   }
00941 
00942   /* the curve we're working in: */
00943   p0 = pp->curve.c[mod(i,m)][2];
00944   p1 = pp->curve.vertex[mod(i+1,m)];
00945   p2 = pp->curve.vertex[mod(j,m)];
00946   p3 = pp->curve.c[mod(j,m)][2];
00947 
00948   /* determine its area */
00949   area = areac[j] - areac[i];
00950   area -= dpara(pp->curve.vertex[0], pp->curve.c[i][2], pp->curve.c[j][2])/2;
00951   if (i>=j) {
00952     area += areac[m];
00953   }
00954 
00955   /* find intersection o of p0p1 and p2p3. Let t,s such that o =
00956      interval(t,p0,p1) = interval(s,p3,p2). Let A be the area of the
00957      triangle (p0,o,p3). */
00958 
00959   A1 = dpara(p0, p1, p2);
00960   A2 = dpara(p0, p1, p3);
00961   A3 = dpara(p0, p2, p3);
00962   /* A4 = dpara(p1, p2, p3); */
00963   A4 = A1+A3-A2;    
00964   
00965   if (A2 == A1) {  /* this should never happen */
00966     return 1;
00967   }
00968 
00969   t = A3/(A3-A4);
00970   s = A2/(A2-A1);
00971   A = A2 * t / 2.0;
00972   
00973   if (A == 0.0) {  /* this should never happen */
00974     return 1;
00975   }
00976 
00977   R = area / A;  /* relative area */
00978   alpha = 2 - sqrt(4 - R / 0.3);  /* overall alpha for p0-o-p3 curve */
00979 
00980   res->c[0] = interval(t * alpha, p0, p1);
00981   res->c[1] = interval(s * alpha, p3, p2);
00982   res->alpha = alpha;
00983   res->t = t;
00984   res->s = s;
00985 
00986   p1 = res->c[0];
00987   p2 = res->c[1];  /* the proposed curve is now (p0,p1,p2,p3) */
00988 
00989   res->pen = 0;
00990 
00991   /* calculate penalty */
00992   /* check tangency with edges */
00993   for (k=mod(i+1,m); k!=j; k=k1) {
00994     k1 = mod(k+1,m);
00995     t = tangent(p0, p1, p2, p3, pp->curve.vertex[k], pp->curve.vertex[k1]);
00996     if (t<-.5) {
00997       return 1;
00998     }
00999     pt = bezier(t, p0, p1, p2, p3);
01000     d = ddist(pp->curve.vertex[k], pp->curve.vertex[k1]);
01001     if (d == 0.0) {  /* this should never happen */
01002       return 1;
01003     }
01004     d1 = dpara(pp->curve.vertex[k], pp->curve.vertex[k1], pt) / d;
01005     if (fabs(d1) > opttolerance) {
01006       return 1;
01007     }
01008     if (iprod(pp->curve.vertex[k], pp->curve.vertex[k1], pt) < 0 || iprod(pp->curve.vertex[k1], pp->curve.vertex[k], pt) < 0) {
01009       return 1;
01010     }
01011     res->pen += sq(d1);
01012   }
01013 
01014   /* check corners */
01015   for (k=i; k!=j; k=k1) {
01016     k1 = mod(k+1,m);
01017     t = tangent(p0, p1, p2, p3, pp->curve.c[k][2], pp->curve.c[k1][2]);
01018     if (t<-.5) {
01019       return 1;
01020     }
01021     pt = bezier(t, p0, p1, p2, p3);
01022     d = ddist(pp->curve.c[k][2], pp->curve.c[k1][2]);
01023     if (d == 0.0) {  /* this should never happen */
01024       return 1;
01025     }
01026     d1 = dpara(pp->curve.c[k][2], pp->curve.c[k1][2], pt) / d;
01027     d2 = dpara(pp->curve.c[k][2], pp->curve.c[k1][2], pp->curve.vertex[k1]) / d;
01028     d2 *= 0.75 * pp->curve.alpha[k1];
01029     if (d2 < 0) {
01030       d1 = -d1;
01031       d2 = -d2;
01032     }
01033     if (d1 < d2 - opttolerance) {
01034       return 1;
01035     }
01036     if (d1 < d2) {
01037       res->pen += sq(d1 - d2);
01038     }
01039   }
01040 
01041   return 0;
01042 }
01043 
01044 /* optimize the path p, replacing sequences of Bezier segments by a
01045    single segment when possible. Return 0 on success, 1 with errno set
01046    on failure. */
01047 static int opticurve(privpath_t *pp, double opttolerance) {
01048   int m = pp->curve.n;
01049   int *pt = NULL;     /* pt[m+1] */
01050   double *pen = NULL; /* pen[m+1] */
01051   int *len = NULL;    /* len[m+1] */
01052   opti_t *opt = NULL; /* opt[m+1] */
01053   int om;
01054   int i,j,r;
01055   opti_t o;
01056   dpoint_t p0;
01057   int i1;
01058   double area;
01059   double alpha;
01060   double *s = NULL;
01061   double *t = NULL;
01062 
01063   int *convc = NULL; /* conv[m]: pre-computed convexities */
01064   double *areac = NULL; /* cumarea[m+1]: cache for fast area computation */
01065 
01066   SAFE_MALLOC(pt, m+1, int);
01067   SAFE_MALLOC(pen, m+1, double);
01068   SAFE_MALLOC(len, m+1, int);
01069   SAFE_MALLOC(opt, m+1, opti_t);
01070   SAFE_MALLOC(convc, m, int);
01071   SAFE_MALLOC(areac, m+1, double);
01072 
01073   /* pre-calculate convexity: +1 = right turn, -1 = left turn, 0 = corner */
01074   for (i=0; i<m; i++) {
01075     if (pp->curve.tag[i] == POTRACE_CURVETO) {
01076       convc[i] = sign(dpara(pp->curve.vertex[mod(i-1,m)], pp->curve.vertex[i], pp->curve.vertex[mod(i+1,m)]));
01077     } else {
01078       convc[i] = 0;
01079     }
01080   }
01081 
01082   /* pre-calculate areas */
01083   area = 0.0;
01084   areac[0] = 0.0;
01085   p0 = pp->curve.vertex[0];
01086   for (i=0; i<m; i++) {
01087     i1 = mod(i+1, m);
01088     if (pp->curve.tag[i1] == POTRACE_CURVETO) {
01089       alpha = pp->curve.alpha[i1];
01090       area += 0.3*alpha*(4-alpha)*dpara(pp->curve.c[i][2], pp->curve.vertex[i1], pp->curve.c[i1][2])/2;
01091       area += dpara(p0, pp->curve.c[i][2], pp->curve.c[i1][2])/2;
01092     }
01093     areac[i+1] = area;
01094   }
01095 
01096   pt[0] = -1;
01097   pen[0] = 0;
01098   len[0] = 0;
01099 
01100   /* Fixme: we always start from a fixed point -- should find the best
01101      curve cyclically ### */
01102 
01103   for (j=1; j<=m; j++) {
01104     /* calculate best path from 0 to j */
01105     pt[j] = j-1;
01106     pen[j] = pen[j-1];
01107     len[j] = len[j-1]+1;
01108 
01109     for (i=j-2; i>=0; i--) {
01110       r = opti_penalty(pp, i, mod(j,m), &o, opttolerance, convc, areac);
01111       if (r) {
01112         break;
01113       }
01114       if (len[j] > len[i]+1 || (len[j] == len[i]+1 && pen[j] > pen[i] + o.pen)) {
01115         pt[j] = i;
01116         pen[j] = pen[i] + o.pen;
01117         len[j] = len[i] + 1;
01118         opt[j] = o;
01119       }
01120     }
01121   }
01122   om = len[m];
01123   r = privcurve_init(&pp->ocurve, om);
01124   if (r) {
01125     goto malloc_error;
01126   }
01127   SAFE_MALLOC(s, om, double);
01128   SAFE_MALLOC(t, om, double);
01129 
01130   j = m;
01131   for (i=om-1; i>=0; i--) {
01132     if (pt[j]==j-1) {
01133       pp->ocurve.tag[i]     = pp->curve.tag[mod(j,m)];
01134       pp->ocurve.c[i][0]    = pp->curve.c[mod(j,m)][0];
01135       pp->ocurve.c[i][1]    = pp->curve.c[mod(j,m)][1];
01136       pp->ocurve.c[i][2]    = pp->curve.c[mod(j,m)][2];
01137       pp->ocurve.vertex[i]  = pp->curve.vertex[mod(j,m)];
01138       pp->ocurve.alpha[i]   = pp->curve.alpha[mod(j,m)];
01139       pp->ocurve.alpha0[i]  = pp->curve.alpha0[mod(j,m)];
01140       pp->ocurve.beta[i]    = pp->curve.beta[mod(j,m)];
01141       s[i] = t[i] = 1.0;
01142     } else {
01143       pp->ocurve.tag[i] = POTRACE_CURVETO;
01144       pp->ocurve.c[i][0] = opt[j].c[0];
01145       pp->ocurve.c[i][1] = opt[j].c[1];
01146       pp->ocurve.c[i][2] = pp->curve.c[mod(j,m)][2];
01147       pp->ocurve.vertex[i] = interval(opt[j].s, pp->curve.c[mod(j,m)][2], pp->curve.vertex[mod(j,m)]);
01148       pp->ocurve.alpha[i] = opt[j].alpha;
01149       pp->ocurve.alpha0[i] = opt[j].alpha;
01150       s[i] = opt[j].s;
01151       t[i] = opt[j].t;
01152     }
01153     j = pt[j];
01154   }
01155 
01156   /* calculate beta parameters */
01157   for (i=0; i<om; i++) {
01158     i1 = mod(i+1,om);
01159     pp->ocurve.beta[i] = s[i] / (s[i] + t[i1]);
01160   }
01161   pp->ocurve.alphacurve = 1;
01162 
01163   free(pt);
01164   free(pen);
01165   free(len);
01166   free(opt);
01167   free(s);
01168   free(t);
01169   free(convc);
01170   free(areac);
01171   return 0;
01172 
01173  malloc_error:
01174   free(pt);
01175   free(pen);
01176   free(len);
01177   free(opt);
01178   free(s);
01179   free(t);
01180   free(convc);
01181   free(areac);
01182   return 1;
01183 }
01184 
01185 /* ---------------------------------------------------------------------- */
01186 
01187 #define TRY(x) if (x) goto try_error
01188 
01189 /* return 0 on success, 1 on error with errno set. */
01190 int process_path(path_t *plist, const potrace_param_t *param, progress_t *progress) {
01191   path_t *p;
01192   double nn = 0, cn = 0;
01193 
01194   if (progress->callback) {
01195     /* precompute task size for progress estimates */
01196     nn = 0;
01197     list_forall (p, plist) {
01198       nn += p->priv->len;
01199     }
01200     cn = 0;
01201   }
01202   
01203   /* call downstream function with each path */
01204   list_forall (p, plist) {
01205     TRY(calc_sums(p->priv));
01206     TRY(calc_lon(p->priv));
01207     TRY(bestpolygon(p->priv));
01208     TRY(adjust_vertices(p->priv));
01209     TRY(smooth(&p->priv->curve, p->sign, param->alphamax));
01210     if (param->opticurve) {
01211       TRY(opticurve(p->priv, param->opttolerance));
01212       p->priv->fcurve = &p->priv->ocurve;
01213     } else {
01214       p->priv->fcurve = &p->priv->curve;
01215     }
01216     privcurve_to_curve(p->priv->fcurve, &p->curve);
01217 
01218     if (progress->callback) {
01219       cn += p->priv->len;
01220       progress_update(cn/nn, progress);
01221     }
01222   }
01223 
01224   progress_update(1.0, progress);
01225 
01226   return 0;
01227 
01228  try_error:
01229   return 1;
01230 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines


portrait_painter
Author(s): Niklas Meinzer, Ina Baumgarten
autogenerated on Wed Dec 26 2012 16:00:43