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: trace.c,v 1.10 2005/02/27 01:54:08 selinger Exp $ */
6
/* transform jaggy paths into smooth curves */
13
#include "potracelib.h"
16
#include "auxiliary.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
/* some useful macros. Note: the "mod" macro works correctly for
26
negative a. Also note that the test for a>=n, while redundant,
27
speeds up the mod function by 70% in the average case (significant
28
since the program spends about 16% of its time here - or 40%
29
without the test). The "floordiv" macro returns the largest integer
30
<= a/n, and again this works correctly for negative a, as long as
31
a,n are integers and n>0. */
33
#define SAFE_MALLOC(var, n, typ) \
34
if ((var = (typ *)malloc((n)*sizeof(typ))) == NULL) goto malloc_error
36
/* ---------------------------------------------------------------------- */
37
/* auxiliary functions */
39
/* return a direction that is 90 degrees counterclockwise from p2-p0,
40
but then restricted to one of the major wind directions (n, nw, w, etc) */
41
static inline point_t dorth_infty(dpoint_t p0, dpoint_t p2) {
44
r.y = sign(p2.x-p0.x);
45
r.x = -sign(p2.y-p0.y);
50
/* return (p1-p0)x(p2-p0), the area of the parallelogram */
51
static inline double dpara(dpoint_t p0, dpoint_t p1, dpoint_t p2) {
52
double x1, y1, x2, y2;
62
/* ddenom/dpara have the property that the square of radius 1 centered
63
at p1 intersects the line p0p2 iff |dpara(p0,p1,p2)| <= ddenom(p0,p2) */
64
static inline double ddenom(dpoint_t p0, dpoint_t p2) {
65
point_t r = dorth_infty(p0, p2);
67
return r.y*(p2.x-p0.x) - r.x*(p2.y-p0.y);
70
/* return 1 if a <= b < c < a, in a cyclic sense (mod n) */
71
static inline int cyclic(int a, int b, int c) {
79
/* determine the center and slope of the line i..j. Assume i<j. Needs
80
"sum" components of p to be set. */
81
static void pointslope(privpath_t *pp, int i, int j, dpoint_t *ctr, dpoint_t *dir) {
85
sums_t *sums = pp->sums;
87
double x, y, x2, xy, y2;
89
double a, b, c, lambda2, l;
90
int r=0; /* rotations from i to j */
109
x = sums[j+1].x-sums[i].x+r*sums[n].x;
110
y = sums[j+1].y-sums[i].y+r*sums[n].y;
111
x2 = sums[j+1].x2-sums[i].x2+r*sums[n].x2;
112
xy = sums[j+1].xy-sums[i].xy+r*sums[n].xy;
113
y2 = sums[j+1].y2-sums[i].y2+r*sums[n].y2;
119
a = (x2-(double)x*x/k)/k;
120
b = (xy-(double)x*y/k)/k;
121
c = (y2-(double)y*y/k)/k;
123
lambda2 = (a+c+sqrt((a-c)*(a-c)+4*b*b))/2; /* larger e.value */
125
/* now find e.vector for lambda2 */
129
if (fabs(a) >= fabs(c)) {
143
dir->x = dir->y = 0; /* sometimes this can happen when k=4:
144
the two eigenvalues coincide */
148
/* the type of (affine) quadratic forms, represented as symmetric 3x3
149
matrices. The value of the quadratic form at a vector (x,y) is v^t
150
Q v, where v = (x,y,1)^t. */
151
typedef double quadform_t[3][3];
153
/* Apply quadratic form Q to vector w = (w.x,w.y) */
154
static inline double quadform(quadform_t Q, dpoint_t w) {
164
for (i=0; i<3; i++) {
165
for (j=0; j<3; j++) {
166
sum += v[i] * Q[i][j] * v[j];
172
/* calculate p1 x p2 */
173
static inline int xprod(point_t p1, point_t p2) {
174
return p1.x*p2.y - p1.y*p2.x;
177
/* calculate (p1-p0)x(p3-p2) */
178
static inline double cprod(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3) {
179
double x1, y1, x2, y2;
186
return x1*y2 - x2*y1;
189
/* calculate (p1-p0)*(p2-p0) */
190
static inline double iprod(dpoint_t p0, dpoint_t p1, dpoint_t p2) {
191
double x1, y1, x2, y2;
198
return x1*x2 + y1*y2;
201
/* calculate (p1-p0)*(p3-p2) */
202
static inline double iprod1(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3) {
203
double x1, y1, x2, y2;
210
return x1*x2 + y1*y2;
213
/* calculate distance between two points */
214
static inline double ddist(dpoint_t p, dpoint_t q) {
215
return sqrt(sq(p.x-q.x)+sq(p.y-q.y));
218
/* calculate point of a bezier curve */
219
static inline dpoint_t bezier(double t, dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3) {
223
/* Note: a good optimizing compiler (such as gcc-3) reduces the
224
following to 16 multiplications, using common subexpression
227
res.x = s*s*s*p0.x + 3*(s*s*t)*p1.x + 3*(t*t*s)*p2.x + t*t*t*p3.x;
228
res.y = s*s*s*p0.y + 3*(s*s*t)*p1.y + 3*(t*t*s)*p2.y + t*t*t*p3.y;
233
/* calculate the point t in [0..1] on the (convex) bezier curve
234
(p0,p1,p2,p3) which is tangent to q1-q0. Return -1.0 if there is no
235
solution in [0..1]. */
236
static double tangent(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3, dpoint_t q0, dpoint_t q1) {
237
double A, B, C; /* (1-t)^2 A + 2(1-t)t B + t^2 C = 0 */
238
double a, b, c; /* a t^2 + b t + c = 0 */
241
A = cprod(p0, p1, q0, q1);
242
B = cprod(p1, p2, q0, q1);
243
C = cprod(p2, p3, q0, q1);
257
r1 = (-b + s) / (2 * a);
258
r2 = (-b - s) / (2 * a);
260
if (r1 >= 0 && r1 <= 1) {
262
} else if (r2 >= 0 && r2 <= 1) {
269
/* ---------------------------------------------------------------------- */
270
/* Preparation: fill in the sum* fields of a path (used for later
271
rapid summing). Return 0 on success, 1 with errno set on
273
static int calc_sums(privpath_t *pp) {
277
SAFE_MALLOC(pp->sums, pp->len+1, sums_t);
280
pp->x0 = pp->pt[0].x;
281
pp->y0 = pp->pt[0].y;
283
/* preparatory computation for later fast summing */
284
pp->sums[0].x2 = pp->sums[0].xy = pp->sums[0].y2 = pp->sums[0].x = pp->sums[0].y = 0;
285
for (i=0; i<n; i++) {
286
x = pp->pt[i].x - pp->x0;
287
y = pp->pt[i].y - pp->y0;
288
pp->sums[i+1].x = pp->sums[i].x + x;
289
pp->sums[i+1].y = pp->sums[i].y + y;
290
pp->sums[i+1].x2 = pp->sums[i].x2 + x*x;
291
pp->sums[i+1].xy = pp->sums[i].xy + x*y;
292
pp->sums[i+1].y2 = pp->sums[i].y2 + y*y;
300
/* ---------------------------------------------------------------------- */
301
/* Stage 1: determine the straight subpaths (Sec. 2.2.1). Fill in the
302
"lon" component of a path object (based on pt/len). For each i,
303
lon[i] is the furthest index such that a straight line can be drawn
304
from i to lon[i]. Return 1 on error with errno set, else 0. */
306
/* this algorithm depends on the fact that the existence of straight
307
subpaths is a triplewise property. I.e., there exists a straight
308
line through squares i0,...,in iff there exists a straight line
309
through i,j,k, for all i0<=i<j<k<=in. (Proof?) */
311
/* this implementation of calc_lon is O(n^2). It replaces an older
312
O(n^3) version. A "constraint" means that future points must
313
satisfy xprod(constraint[0], cur) >= 0 and xprod(constraint[1],
316
/* Remark for potrace 1.1: the current implementation of calc_lon is
317
more complex than the implementation found in potrace 1.0, but it
318
is considerably faster. The introduction of the "nc" data structure
319
means that we only have to test the constraints for "corner"
320
points. On a typical input file, this speeds up the calc_lon
321
function by a factor of 31.2, thereby decreasing its time share
322
within the overall potrace algorithm from 72.6% to 7.82%, and
323
speeding up the overall algorithm by a factor of 3.36. On another
324
input file, calc_lon was sped up by a factor of 6.7, decreasing its
325
time share from 51.4% to 13.61%, and speeding up the overall
326
algorithm by a factor of 1.78. In any case, the savings are
329
/* returns 0 on success, 1 on error with errno set */
330
static int calc_lon(privpath_t *pp) {
331
point_t *pt = pp->pt;
335
point_t constraint[2];
338
int *pivk = NULL; /* pivk[n] */
339
int *nc = NULL; /* nc[n]: next corner */
340
point_t dk; /* direction of k-k1 */
343
SAFE_MALLOC(pivk, n, int);
344
SAFE_MALLOC(nc, n, int);
346
/* initialize the nc data structure. Point from each point to the
347
furthest future point to which it is connected by a vertical or
348
horizontal segment. We take advantage of the fact that there is
349
always a direction change at 0 (due to the path decomposition
350
algorithm). But even if this were not so, there is no harm, as
351
in practice, correctness does not depend on the word "furthest"
354
for (i=n-1; i>=0; i--) {
355
if (pt[i].x != pt[k].x && pt[i].y != pt[k].y) {
356
k = i+1; /* necessarily i<n-1 in this case */
361
SAFE_MALLOC(pp->lon, n, int);
363
/* determine pivot points: for each i, let pivk[i] be the furthest k
364
such that all j with i<j<k lie on a line connecting i,k. */
366
for (i=n-1; i>=0; i--) {
367
ct[0] = ct[1] = ct[2] = ct[3] = 0;
369
/* keep track of "directions" that have occurred */
370
dir = (3+3*(pt[mod(i+1,n)].x-pt[i].x)+(pt[mod(i+1,n)].y-pt[i].y))/2;
378
/* find the next k such that no straight line from i to k */
383
dir = (3+3*sign(pt[k].x-pt[k1].x)+sign(pt[k].y-pt[k1].y))/2;
386
/* if all four "directions" have occurred, cut this path */
387
if (ct[0] && ct[1] && ct[2] && ct[3]) {
392
cur.x = pt[k].x - pt[i].x;
393
cur.y = pt[k].y - pt[i].y;
395
/* see if current constraint is violated */
396
if (xprod(constraint[0], cur) < 0 || xprod(constraint[1], cur) > 0) {
397
goto constraint_viol;
400
/* else, update constraint */
401
if (abs(cur.x) <= 1 && abs(cur.y) <= 1) {
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[0], off) >= 0) {
409
off.x = cur.x + ((cur.y<=0 && (cur.y<0 || cur.x<0)) ? 1 : -1);
410
off.y = cur.y + ((cur.x>=0 && (cur.x>0 || cur.y<0)) ? 1 : -1);
411
if (xprod(constraint[1], off) <= 0) {
417
if (!cyclic(k,i,k1)) {
422
/* k1 was the last "corner" satisfying the current constraint, and
423
k is the first one violating it. We now need to find the last
424
point along k1..k which satisfied the constraint. */
425
dk.x = sign(pt[k].x-pt[k1].x);
426
dk.y = sign(pt[k].y-pt[k1].y);
427
cur.x = pt[k1].x - pt[i].x;
428
cur.y = pt[k1].y - pt[i].y;
429
/* find largest integer j such that xprod(constraint[0], cur+j*dk)
430
>= 0 and xprod(constraint[1], cur+j*dk) <= 0. Use bilinearity
432
a = xprod(constraint[0], cur);
433
b = xprod(constraint[0], dk);
434
c = xprod(constraint[1], cur);
435
d = xprod(constraint[1], dk);
436
/* find largest integer j such that a+j*b>=0 and c+j*d<=0. This
437
can be solved with integer arithmetic. */
443
j = min(j, floordiv(-c,d));
445
pivk[i] = mod(k1+j,n);
450
/* clean up: for each i, let lon[i] be the largest k such that for
451
all i' with i<=i'<k, i'<k<=pivk[i']. */
455
for (i=n-2; i>=0; i--) {
456
if (cyclic(i+1,pivk[i],j)) {
462
for (i=n-1; cyclic(mod(i+1,n),j,pp->lon[i]); i--) {
477
/* ---------------------------------------------------------------------- */
478
/* Stage 2: calculate the optimal polygon (Sec. 2.2.2-2.2.4). */
480
/* Auxiliary function: calculate the penalty of an edge from i to j in
481
the given path. This needs the "lon" and "sum*" data. */
483
static double penalty3(privpath_t *pp, int i, int j) {
485
point_t *pt = pp->pt;
486
sums_t *sums = pp->sums;
488
/* assume 0<=i<j<=n */
489
double x, y, x2, xy, y2;
492
double px, py, ex, ey;
494
int r=0; /* rotations from i to j */
501
x = sums[j+1].x-sums[i].x+r*sums[n].x;
502
y = sums[j+1].y-sums[i].y+r*sums[n].y;
503
x2 = sums[j+1].x2-sums[i].x2+r*sums[n].x2;
504
xy = sums[j+1].xy-sums[i].xy+r*sums[n].xy;
505
y2 = sums[j+1].y2-sums[i].y2+r*sums[n].y2;
508
px = (pt[i].x+pt[j].x)/2.0-pt[0].x;
509
py = (pt[i].y+pt[j].y)/2.0-pt[0].y;
510
ey = (pt[j].x-pt[i].x);
511
ex = -(pt[j].y-pt[i].y);
513
a = ((x2-2*x*px)/k+px*px);
514
b = ((xy-x*py-y*px)/k+px*py);
515
c = ((y2-2*y*py)/k+py*py);
517
s = ex*ex*a + 2*ex*ey*b + ey*ey*c;
522
/* find the optimal polygon. Fill in the m and po components. Return 1
523
on failure with errno set, else 0. Non-cyclic version: assumes i=0
524
is in the polygon. Fixme: ### implement cyclic version. */
525
static int bestpolygon(privpath_t *pp)
529
double *pen = NULL; /* pen[n+1]: penalty vector */
530
int *prev = NULL; /* prev[n+1]: best path pointer vector */
531
int *clip0 = NULL; /* clip0[n]: longest segment pointer, non-cyclic */
532
int *clip1 = NULL; /* clip1[n+1]: backwards segment pointer, non-cyclic */
533
int *seg0 = NULL; /* seg0[m+1]: forward segment bounds, m<=n */
534
int *seg1 = NULL; /* seg1[m+1]: backward segment bounds, m<=n */
539
SAFE_MALLOC(pen, n+1, double);
540
SAFE_MALLOC(prev, n+1, int);
541
SAFE_MALLOC(clip0, n, int);
542
SAFE_MALLOC(clip1, n+1, int);
543
SAFE_MALLOC(seg0, n+1, int);
544
SAFE_MALLOC(seg1, n+1, int);
546
/* calculate clipped paths */
547
for (i=0; i<n; i++) {
548
c = mod(pp->lon[mod(i-1,n)]-1,n);
559
/* calculate backwards path clipping, non-cyclic. j <= clip0[i] iff
560
clip1[j] <= i, for i,j=0..n. */
562
for (i=0; i<n; i++) {
563
while (j <= clip0[i]) {
569
/* calculate seg0[j] = longest path from 0 with j segments */
571
for (j=0; i<n; j++) {
578
/* calculate seg1[j] = longest path to n with m-j segments */
580
for (j=m; j>0; j--) {
586
/* now find the shortest path with m segments, based on penalty3 */
587
/* note: the outer 2 loops jointly have at most n interations, thus
588
the worst-case behavior here is quadratic. In practice, it is
589
close to linear since the inner loop tends to be short. */
591
for (j=1; j<=m; j++) {
592
for (i=seg1[j]; i<=seg0[j]; i++) {
594
for (k=seg0[j-1]; k>=clip1[i]; k--) {
595
thispen = penalty3(pp, k, i) + pen[k];
596
if (best < 0 || thispen < best) {
606
SAFE_MALLOC(pp->po, m, int);
608
/* read off shortest path */
609
for (i=n, j=m-1; i>0; j--) {
632
/* ---------------------------------------------------------------------- */
633
/* Stage 3: vertex adjustment (Sec. 2.3.1). */
635
/* Adjust vertices of optimal polygon: calculate the intersection of
636
the two "optimal" line segments, then move it into the unit square
637
if it lies outside. Return 1 with errno set on error; 0 on
640
static int adjust_vertices(privpath_t *pp) {
644
point_t *pt = pp->pt;
648
dpoint_t *ctr = NULL; /* ctr[m] */
649
dpoint_t *dir = NULL; /* dir[m] */
650
quadform_t *q = NULL; /* q[m] */
657
SAFE_MALLOC(ctr, m, dpoint_t);
658
SAFE_MALLOC(dir, m, dpoint_t);
659
SAFE_MALLOC(q, m, quadform_t);
661
r = privcurve_init(&pp->curve, m);
666
/* calculate "optimal" point-slope representation for each line
668
for (i=0; i<m; i++) {
670
j = mod(j-po[i],n)+po[i];
671
pointslope(pp, po[i], j, &ctr[i], &dir[i]);
674
/* represent each line segment as a singular quadratic form; the
675
distance of a point (x,y) from the line segment will be
676
(x,y,1)Q(x,y,1)^t, where Q=q[i]. */
677
for (i=0; i<m; i++) {
678
d = sq(dir[i].x) + sq(dir[i].y);
680
for (j=0; j<3; j++) {
681
for (k=0; k<3; k++) {
688
v[2] = - v[1] * ctr[i].y - v[0] * ctr[i].x;
689
for (l=0; l<3; l++) {
690
for (k=0; k<3; k++) {
691
q[i][l][k] = v[l] * v[k] / d;
697
/* now calculate the "intersections" of consecutive segments.
698
Instead of using the actual intersection, we find the point
699
within a given unit square which minimizes the square distance to
701
for (i=0; i<m; i++) {
706
double min, cand; /* minimum and candidate for minimum of quad. form */
707
double xmin, ymin; /* coordinates of minimum */
710
/* let s be the vertex, in coordinates relative to x0/y0 */
711
s.x = pt[po[i]].x-x0;
712
s.y = pt[po[i]].y-y0;
714
/* intersect segments i-1 and i */
718
/* add quadratic forms */
719
for (l=0; l<3; l++) {
720
for (k=0; k<3; k++) {
721
Q[l][k] = q[j][l][k] + q[i][l][k];
726
/* minimize the quadratic form Q on the unit square */
727
/* find intersection */
729
#ifdef HAVE_GCC_LOOP_BUG
730
/* work around gcc bug #12243 */
734
det = Q[0][0]*Q[1][1] - Q[0][1]*Q[1][0];
736
w.x = (-Q[0][2]*Q[1][1] + Q[1][2]*Q[0][1]) / det;
737
w.y = ( Q[0][2]*Q[1][0] - Q[1][2]*Q[0][0]) / det;
741
/* matrix is singular - lines are parallel. Add another,
742
orthogonal axis, through the center of the unit square */
743
if (Q[0][0]>Q[1][1]) {
746
} else if (Q[1][1]) {
753
d = sq(v[0]) + sq(v[1]);
754
v[2] = - v[1] * s.y - v[0] * s.x;
755
for (l=0; l<3; l++) {
756
for (k=0; k<3; k++) {
757
Q[l][k] += v[l] * v[k] / d;
763
if (dx <= .5 && dy <= .5) {
764
pp->curve.vertex[i].x = w.x+x0;
765
pp->curve.vertex[i].y = w.y+y0;
769
/* the minimum was not in the unit square; now minimize quadratic
770
on boundary of square */
771
min = quadform(Q, s);
775
if (Q[0][0] == 0.0) {
778
for (z=0; z<2; z++) { /* value of the y-coordinate */
780
w.x = - (Q[0][1] * w.y + Q[0][2]) / Q[0][0];
782
cand = quadform(Q, w);
783
if (dx <= .5 && cand < min) {
790
if (Q[1][1] == 0.0) {
793
for (z=0; z<2; z++) { /* value of the x-coordinate */
795
w.y = - (Q[1][0] * w.x + Q[1][2]) / Q[1][1];
797
cand = quadform(Q, w);
798
if (dy <= .5 && cand < min) {
805
/* check four corners */
806
for (l=0; l<2; l++) {
807
for (k=0; k<2; k++) {
810
cand = quadform(Q, w);
819
pp->curve.vertex[i].x = xmin + x0;
820
pp->curve.vertex[i].y = ymin + y0;
836
/* ---------------------------------------------------------------------- */
837
/* Stage 4: smoothing and corner analysis (Sec. 2.3.3) */
839
/* Always succeeds and returns 0 */
840
static int smooth(privcurve_t *curve, int sign, double alphamax) {
844
double dd, denom, alpha;
848
/* reverse orientation of negative paths */
849
for (i=0, j=m-1; i<j; i++, j--) {
851
tmp = curve->vertex[i];
852
curve->vertex[i] = curve->vertex[j];
853
curve->vertex[j] = tmp;
857
/* examine each vertex and find its best fit */
858
for (i=0; i<m; i++) {
861
p4 = interval(1/2.0, curve->vertex[k], curve->vertex[j]);
863
denom = ddenom(curve->vertex[i], curve->vertex[k]);
865
dd = dpara(curve->vertex[i], curve->vertex[j], curve->vertex[k]) / denom;
867
alpha = dd>1 ? (1 - 1.0/dd) : 0;
868
alpha = alpha / 0.75;
872
curve->alpha0[j] = alpha; /* remember "original" value of alpha */
874
if (alpha > alphamax) { /* pointed corner */
875
curve->tag[j] = POTRACE_CORNER;
876
curve->c[j][1] = curve->vertex[j];
881
} else if (alpha > 1) {
884
p2 = interval(.5+.5*alpha, curve->vertex[i], curve->vertex[j]);
885
p3 = interval(.5+.5*alpha, curve->vertex[k], curve->vertex[j]);
886
curve->tag[j] = POTRACE_CURVETO;
891
curve->alpha[j] = alpha; /* store the "cropped" value of alpha */
892
curve->beta[j] = 0.5;
898
/* ---------------------------------------------------------------------- */
899
/* Stage 5: Curve optimization (Sec. 2.4) */
901
/* a private type for the result of opti_penalty */
903
double pen; /* penalty */
904
dpoint_t c[2]; /* curve parameters */
905
double t, s; /* curve parameters */
906
double alpha; /* curve parameter */
908
typedef struct opti_s opti_t;
910
/* calculate best fit from i+.5 to j+.5. Assume i<j (cyclically).
911
Return 0 and set badness and parameters (alpha, beta), if
912
possible. Return 1 if impossible. */
913
static int opti_penalty(privpath_t *pp, int i, int j, opti_t *res, double opttolerance, int *convc, double *areac) {
915
int k, k1, k2, conv, i1;
916
double area, alpha, d, d1, d2;
917
dpoint_t p0, p1, p2, p3, pt;
918
double A, R, A1, A2, A3, A4;
921
/* check convexity, corner-freeness, and maximum bend < 179 degrees */
923
if (i==j) { /* sanity - a full loop can never be an opticurve */
934
d = ddist(pp->curve.vertex[i], pp->curve.vertex[i1]);
935
for (k=k1; k!=j; k=k1) {
938
if (convc[k1] != conv) {
941
if (sign(cprod(pp->curve.vertex[i], pp->curve.vertex[i1], pp->curve.vertex[k1], pp->curve.vertex[k2])) != conv) {
944
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) {
949
/* the curve we're working in: */
950
p0 = pp->curve.c[mod(i,m)][2];
951
p1 = pp->curve.vertex[mod(i+1,m)];
952
p2 = pp->curve.vertex[mod(j,m)];
953
p3 = pp->curve.c[mod(j,m)][2];
955
/* determine its area */
956
area = areac[j] - areac[i];
957
area -= dpara(pp->curve.vertex[0], pp->curve.c[i][2], pp->curve.c[j][2])/2;
962
/* find intersection o of p0p1 and p2p3. Let t,s such that o =
963
interval(t,p0,p1) = interval(s,p3,p2). Let A be the area of the
964
triangle (p0,o,p3). */
966
A1 = dpara(p0, p1, p2);
967
A2 = dpara(p0, p1, p3);
968
A3 = dpara(p0, p2, p3);
969
/* A4 = dpara(p1, p2, p3); */
972
if (A2 == A1) { /* this should never happen */
980
if (A == 0.0) { /* this should never happen */
984
R = area / A; /* relative area */
985
alpha = 2 - sqrt(4 - R / 0.3); /* overall alpha for p0-o-p3 curve */
987
res->c[0] = interval(t * alpha, p0, p1);
988
res->c[1] = interval(s * alpha, p3, p2);
994
p2 = res->c[1]; /* the proposed curve is now (p0,p1,p2,p3) */
998
/* calculate penalty */
999
/* check tangency with edges */
1000
for (k=mod(i+1,m); k!=j; k=k1) {
1002
t = tangent(p0, p1, p2, p3, pp->curve.vertex[k], pp->curve.vertex[k1]);
1006
pt = bezier(t, p0, p1, p2, p3);
1007
d = ddist(pp->curve.vertex[k], pp->curve.vertex[k1]);
1008
if (d == 0.0) { /* this should never happen */
1011
d1 = dpara(pp->curve.vertex[k], pp->curve.vertex[k1], pt) / d;
1012
if (fabs(d1) > opttolerance) {
1015
if (iprod(pp->curve.vertex[k], pp->curve.vertex[k1], pt) < 0 || iprod(pp->curve.vertex[k1], pp->curve.vertex[k], pt) < 0) {
1022
for (k=i; k!=j; k=k1) {
1024
t = tangent(p0, p1, p2, p3, pp->curve.c[k][2], pp->curve.c[k1][2]);
1028
pt = bezier(t, p0, p1, p2, p3);
1029
d = ddist(pp->curve.c[k][2], pp->curve.c[k1][2]);
1030
if (d == 0.0) { /* this should never happen */
1033
d1 = dpara(pp->curve.c[k][2], pp->curve.c[k1][2], pt) / d;
1034
d2 = dpara(pp->curve.c[k][2], pp->curve.c[k1][2], pp->curve.vertex[k1]) / d;
1035
d2 *= 0.75 * pp->curve.alpha[k1];
1040
if (d1 < d2 - opttolerance) {
1044
res->pen += sq(d1 - d2);
1051
/* optimize the path p, replacing sequences of Bezier segments by a
1052
single segment when possible. Return 0 on success, 1 with errno set
1054
static int opticurve(privpath_t *pp, double opttolerance) {
1055
int m = pp->curve.n;
1056
int *pt = NULL; /* pt[m+1] */
1057
double *pen = NULL; /* pen[m+1] */
1058
int *len = NULL; /* len[m+1] */
1059
opti_t *opt = NULL; /* opt[m+1] */
1070
int *convc = NULL; /* conv[m]: pre-computed convexities */
1071
double *areac = NULL; /* cumarea[m+1]: cache for fast area computation */
1073
SAFE_MALLOC(pt, m+1, int);
1074
SAFE_MALLOC(pen, m+1, double);
1075
SAFE_MALLOC(len, m+1, int);
1076
SAFE_MALLOC(opt, m+1, opti_t);
1077
SAFE_MALLOC(convc, m, int);
1078
SAFE_MALLOC(areac, m+1, double);
1080
/* pre-calculate convexity: +1 = right turn, -1 = left turn, 0 = corner */
1081
for (i=0; i<m; i++) {
1082
if (pp->curve.tag[i] == POTRACE_CURVETO) {
1083
convc[i] = sign(dpara(pp->curve.vertex[mod(i-1,m)], pp->curve.vertex[i], pp->curve.vertex[mod(i+1,m)]));
1089
/* pre-calculate areas */
1092
p0 = pp->curve.vertex[0];
1093
for (i=0; i<m; i++) {
1095
if (pp->curve.tag[i1] == POTRACE_CURVETO) {
1096
alpha = pp->curve.alpha[i1];
1097
area += 0.3*alpha*(4-alpha)*dpara(pp->curve.c[i][2], pp->curve.vertex[i1], pp->curve.c[i1][2])/2;
1098
area += dpara(p0, pp->curve.c[i][2], pp->curve.c[i1][2])/2;
1107
/* Fixme: we always start from a fixed point -- should find the best
1108
curve cyclically ### */
1110
for (j=1; j<=m; j++) {
1111
/* calculate best path from 0 to j */
1114
len[j] = len[j-1]+1;
1116
for (i=j-2; i>=0; i--) {
1117
r = opti_penalty(pp, i, mod(j,m), &o, opttolerance, convc, areac);
1121
if (len[j] > len[i]+1 || (len[j] == len[i]+1 && pen[j] > pen[i] + o.pen)) {
1123
pen[j] = pen[i] + o.pen;
1124
len[j] = len[i] + 1;
1130
r = privcurve_init(&pp->ocurve, om);
1134
SAFE_MALLOC(s, om, double);
1135
SAFE_MALLOC(t, om, double);
1138
for (i=om-1; i>=0; i--) {
1140
pp->ocurve.tag[i] = pp->curve.tag[mod(j,m)];
1141
pp->ocurve.c[i][0] = pp->curve.c[mod(j,m)][0];
1142
pp->ocurve.c[i][1] = pp->curve.c[mod(j,m)][1];
1143
pp->ocurve.c[i][2] = pp->curve.c[mod(j,m)][2];
1144
pp->ocurve.vertex[i] = pp->curve.vertex[mod(j,m)];
1145
pp->ocurve.alpha[i] = pp->curve.alpha[mod(j,m)];
1146
pp->ocurve.alpha0[i] = pp->curve.alpha0[mod(j,m)];
1147
pp->ocurve.beta[i] = pp->curve.beta[mod(j,m)];
1150
pp->ocurve.tag[i] = POTRACE_CURVETO;
1151
pp->ocurve.c[i][0] = opt[j].c[0];
1152
pp->ocurve.c[i][1] = opt[j].c[1];
1153
pp->ocurve.c[i][2] = pp->curve.c[mod(j,m)][2];
1154
pp->ocurve.vertex[i] = interval(opt[j].s, pp->curve.c[mod(j,m)][2], pp->curve.vertex[mod(j,m)]);
1155
pp->ocurve.alpha[i] = opt[j].alpha;
1156
pp->ocurve.alpha0[i] = opt[j].alpha;
1163
/* calculate beta parameters */
1164
for (i=0; i<om; i++) {
1166
pp->ocurve.beta[i] = s[i] / (s[i] + t[i1]);
1191
/* ---------------------------------------------------------------------- */
1193
#define TRY(x) if (x) goto try_error
1195
/* return 0 on success, 1 on error with errno set. */
1196
int process_path(path_t *plist, const potrace_param_t *param, progress_t *progress) {
1198
double nn = 0, cn = 0;
1200
if (progress->callback) {
1201
/* precompute task size for progress estimates */
1203
list_forall (p, plist) {
1209
/* call downstream function with each path */
1210
list_forall (p, plist) {
1211
TRY(calc_sums(p->priv));
1212
TRY(calc_lon(p->priv));
1213
TRY(bestpolygon(p->priv));
1214
TRY(adjust_vertices(p->priv));
1215
TRY(smooth(&p->priv->curve, p->sign, param->alphamax));
1216
if (param->opticurve) {
1217
TRY(opticurve(p->priv, param->opttolerance));
1218
p->priv->fcurve = &p->priv->ocurve;
1220
p->priv->fcurve = &p->priv->curve;
1222
privcurve_to_curve(p->priv->fcurve, &p->curve);
1224
if (progress->callback) {
1226
progress_update(cn/nn, progress);
1230
progress_update(1.0, progress);
added
removed
3rdparty/potrace/trace.c
