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