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