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