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