gvsig-raster / libjni-potrace / trunk / libjni-potrace / resources / potrace-1.8 / src / trace.c @ 1780
History | View | Annotate | Download (31.2 KB)
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. */
|
4 |
|
5 |
/* $Id: trace.c 147 2007-04-09 00:44:09Z selinger $ */
|
6 |
/* transform jaggy paths into smooth curves */
|
7 |
|
8 |
#include <stdio.h> |
9 |
#include <math.h> |
10 |
#include <stdlib.h> |
11 |
#include <string.h> |
12 |
|
13 |
#include "potracelib.h" |
14 |
#include "curve.h" |
15 |
#include "lists.h" |
16 |
#include "auxiliary.h" |
17 |
#include "trace.h" |
18 |
#include "progress.h" |
19 |
|
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 */ |
23 |
|
24 |
/* ---------------------------------------------------------------------- */
|
25 |
#define SAFE_MALLOC(var, n, typ) \
|
26 |
if ((var = (typ *)malloc((n)*sizeof(typ))) == NULL) goto malloc_error |
27 |
|
28 |
/* ---------------------------------------------------------------------- */
|
29 |
/* auxiliary functions */
|
30 |
|
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; |
35 |
|
36 |
r.y = sign(p2.x-p0.x); |
37 |
r.x = -sign(p2.y-p0.y); |
38 |
|
39 |
return r;
|
40 |
} |
41 |
|
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;
|
45 |
|
46 |
x1 = p1.x-p0.x; |
47 |
y1 = p1.y-p0.y; |
48 |
x2 = p2.x-p0.x; |
49 |
y2 = p2.y-p0.y; |
50 |
|
51 |
return x1*y2 - x2*y1;
|
52 |
} |
53 |
|
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); |
58 |
|
59 |
return r.y*(p2.x-p0.x) - r.x*(p2.y-p0.y);
|
60 |
} |
61 |
|
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 |
} |
70 |
|
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 */
|
75 |
|
76 |
int n = pp->len;
|
77 |
sums_t *sums = pp->sums; |
78 |
|
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 */ |
83 |
|
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 |
} |
100 |
|
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;
|
107 |
|
108 |
ctr->x = x/k; |
109 |
ctr->y = y/k; |
110 |
|
111 |
a = (x2-(double)x*x/k)/k;
|
112 |
b = (xy-(double)x*y/k)/k;
|
113 |
c = (y2-(double)y*y/k)/k;
|
114 |
|
115 |
lambda2 = (a+c+sqrt((a-c)*(a-c)+4*b*b))/2; /* larger e.value */ |
116 |
|
117 |
/* now find e.vector for lambda2 */
|
118 |
a -= lambda2; |
119 |
c -= lambda2; |
120 |
|
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 |
} |
139 |
|
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]; |
144 |
|
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;
|
150 |
|
151 |
v[0] = w.x;
|
152 |
v[1] = w.y;
|
153 |
v[2] = 1; |
154 |
sum = 0.0; |
155 |
|
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 |
} |
163 |
|
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 |
} |
168 |
|
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;
|
172 |
|
173 |
x1 = p1.x - p0.x; |
174 |
y1 = p1.y - p0.y; |
175 |
x2 = p3.x - p2.x; |
176 |
y2 = p3.y - p2.y; |
177 |
|
178 |
return x1*y2 - x2*y1;
|
179 |
} |
180 |
|
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;
|
184 |
|
185 |
x1 = p1.x - p0.x; |
186 |
y1 = p1.y - p0.y; |
187 |
x2 = p2.x - p0.x; |
188 |
y2 = p2.y - p0.y; |
189 |
|
190 |
return x1*x2 + y1*y2;
|
191 |
} |
192 |
|
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;
|
196 |
|
197 |
x1 = p1.x - p0.x; |
198 |
y1 = p1.y - p0.y; |
199 |
x2 = p3.x - p2.x; |
200 |
y2 = p3.y - p2.y; |
201 |
|
202 |
return x1*x2 + y1*y2;
|
203 |
} |
204 |
|
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 |
} |
209 |
|
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; |
214 |
|
215 |
/* Note: a good optimizing compiler (such as gcc-3) reduces the
|
216 |
following to 16 multiplications, using common subexpression
|
217 |
elimination. */
|
218 |
|
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; |
221 |
|
222 |
return res;
|
223 |
} |
224 |
|
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;
|
232 |
|
233 |
A = cprod(p0, p1, q0, q1); |
234 |
B = cprod(p1, p2, q0, q1); |
235 |
C = cprod(p2, p3, q0, q1); |
236 |
|
237 |
a = A - 2*B + C;
|
238 |
b = -2*A + 2*B; |
239 |
c = A; |
240 |
|
241 |
d = b*b - 4*a*c;
|
242 |
|
243 |
if (a==0 || d<0) { |
244 |
return -1.0; |
245 |
} |
246 |
|
247 |
s = sqrt(d); |
248 |
|
249 |
r1 = (-b + s) / (2 * a);
|
250 |
r2 = (-b - s) / (2 * a);
|
251 |
|
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 |
} |
260 |
|
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;
|
268 |
|
269 |
SAFE_MALLOC(pp->sums, pp->len+1, sums_t);
|
270 |
|
271 |
/* origin */
|
272 |
pp->x0 = pp->pt[0].x;
|
273 |
pp->y0 = pp->pt[0].y;
|
274 |
|
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; |
287 |
|
288 |
malloc_error:
|
289 |
return 1; |
290 |
} |
291 |
|
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. */
|
297 |
|
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?) */
|
302 |
|
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. */
|
307 |
|
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. */
|
320 |
|
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;
|
334 |
|
335 |
SAFE_MALLOC(pivk, n, int);
|
336 |
SAFE_MALLOC(nc, n, int);
|
337 |
|
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 |
} |
352 |
|
353 |
SAFE_MALLOC(pp->lon, n, int);
|
354 |
|
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. */
|
357 |
|
358 |
for (i=n-1; i>=0; i--) { |
359 |
ct[0] = ct[1] = ct[2] = ct[3] = 0; |
360 |
|
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]++; |
364 |
|
365 |
constraint[0].x = 0; |
366 |
constraint[0].y = 0; |
367 |
constraint[1].x = 0; |
368 |
constraint[1].y = 0; |
369 |
|
370 |
/* find the next k such that no straight line from i to k */
|
371 |
k = nc[i]; |
372 |
k1 = i; |
373 |
while (1) { |
374 |
|
375 |
dir = (3+3*sign(pt[k].x-pt[k1].x)+sign(pt[k].y-pt[k1].y))/2; |
376 |
ct[dir]++; |
377 |
|
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 |
} |
383 |
|
384 |
cur.x = pt[k].x - pt[i].x; |
385 |
cur.y = pt[k].y - pt[i].y; |
386 |
|
387 |
/* see if current constraint is violated */
|
388 |
if (xprod(constraint[0], cur) < 0 || xprod(constraint[1], cur) > 0) { |
389 |
goto constraint_viol;
|
390 |
} |
391 |
|
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 */
|
441 |
|
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']. */
|
444 |
|
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 |
} |
453 |
|
454 |
for (i=n-1; cyclic(mod(i+1,n),j,pp->lon[i]); i--) { |
455 |
pp->lon[i] = j; |
456 |
} |
457 |
|
458 |
free(pivk); |
459 |
free(nc); |
460 |
return 0; |
461 |
|
462 |
malloc_error:
|
463 |
free(pivk); |
464 |
free(nc); |
465 |
return 1; |
466 |
} |
467 |
|
468 |
|
469 |
/* ---------------------------------------------------------------------- */
|
470 |
/* Stage 2: calculate the optimal polygon (Sec. 2.2.2-2.2.4). */
|
471 |
|
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. */
|
474 |
|
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; |
479 |
|
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;
|
485 |
|
486 |
int r=0; /* rotations from i to j */ |
487 |
|
488 |
if (j>=n) {
|
489 |
j-=n; |
490 |
r+=1;
|
491 |
} |
492 |
|
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;
|
499 |
|
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); |
504 |
|
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);
|
508 |
|
509 |
s = ex*ex*a + 2*ex*ey*b + ey*ey*c;
|
510 |
|
511 |
return sqrt(s);
|
512 |
} |
513 |
|
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;
|
530 |
|
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); |
537 |
|
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 |
} |
550 |
|
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 |
} |
560 |
|
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; |
569 |
|
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; |
577 |
|
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 |
} |
596 |
|
597 |
pp->m = m; |
598 |
SAFE_MALLOC(pp->po, m, int);
|
599 |
|
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 |
} |
605 |
|
606 |
free(pen); |
607 |
free(prev); |
608 |
free(clip0); |
609 |
free(clip1); |
610 |
free(seg0); |
611 |
free(seg1); |
612 |
return 0; |
613 |
|
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 |
} |
623 |
|
624 |
/* ---------------------------------------------------------------------- */
|
625 |
/* Stage 3: vertex adjustment (Sec. 2.3.1). */
|
626 |
|
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. */
|
631 |
|
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;
|
639 |
|
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;
|
648 |
|
649 |
SAFE_MALLOC(ctr, m, dpoint_t); |
650 |
SAFE_MALLOC(dir, m, dpoint_t); |
651 |
SAFE_MALLOC(q, m, quadform_t); |
652 |
|
653 |
r = privcurve_init(&pp->curve, m); |
654 |
if (r) {
|
655 |
goto malloc_error;
|
656 |
} |
657 |
|
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 |
} |
665 |
|
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 |
} |
688 |
|
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;
|
701 |
|
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; |
705 |
|
706 |
/* intersect segments i-1 and i */
|
707 |
|
708 |
j = mod(i-1,m);
|
709 |
|
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 |
} |
716 |
|
717 |
while(1) { |
718 |
/* minimize the quadratic form Q on the unit square */
|
719 |
/* find intersection */
|
720 |
|
721 |
#ifdef HAVE_GCC_LOOP_BUG
|
722 |
/* work around gcc bug #12243 */
|
723 |
free(NULL);
|
724 |
#endif
|
725 |
|
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 |
} |
732 |
|
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 |
} |
760 |
|
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; |
766 |
|
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 |
} |
810 |
|
811 |
pp->curve.vertex[i].x = xmin + x0; |
812 |
pp->curve.vertex[i].y = ymin + y0; |
813 |
continue;
|
814 |
} |
815 |
|
816 |
free(ctr); |
817 |
free(dir); |
818 |
free(q); |
819 |
return 0; |
820 |
|
821 |
malloc_error:
|
822 |
free(ctr); |
823 |
free(dir); |
824 |
free(q); |
825 |
return 1; |
826 |
} |
827 |
|
828 |
/* ---------------------------------------------------------------------- */
|
829 |
/* Stage 4: smoothing and corner analysis (Sec. 2.3.3) */
|
830 |
|
831 |
/* Always succeeds and returns 0 */
|
832 |
static int smooth(privcurve_t *curve, int sign, double alphamax) { |
833 |
int m = curve->n;
|
834 |
|
835 |
int i, j, k;
|
836 |
double dd, denom, alpha;
|
837 |
dpoint_t p2, p3, p4; |
838 |
|
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 |
} |
848 |
|
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]); |
854 |
|
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 */
|
865 |
|
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;
|
887 |
|
888 |
return 0; |
889 |
} |
890 |
|
891 |
/* ---------------------------------------------------------------------- */
|
892 |
/* Stage 5: Curve optimization (Sec. 2.4) */
|
893 |
|
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; |
902 |
|
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;
|
913 |
|
914 |
/* check convexity, corner-freeness, and maximum bend < 179 degrees */
|
915 |
|
916 |
if (i==j) { /* sanity - a full loop can never be an opticurve */ |
917 |
return 1; |
918 |
} |
919 |
|
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 |
} |
941 |
|
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];
|
947 |
|
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 |
} |
954 |
|
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). */
|
958 |
|
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; |
964 |
|
965 |
if (A2 == A1) { /* this should never happen */ |
966 |
return 1; |
967 |
} |
968 |
|
969 |
t = A3/(A3-A4); |
970 |
s = A2/(A2-A1); |
971 |
A = A2 * t / 2.0; |
972 |
|
973 |
if (A == 0.0) { /* this should never happen */ |
974 |
return 1; |
975 |
} |
976 |
|
977 |
R = area / A; /* relative area */
|
978 |
alpha = 2 - sqrt(4 - R / 0.3); /* overall alpha for p0-o-p3 curve */ |
979 |
|
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; |
985 |
|
986 |
p1 = res->c[0];
|
987 |
p2 = res->c[1]; /* the proposed curve is now (p0,p1,p2,p3) */ |
988 |
|
989 |
res->pen = 0;
|
990 |
|
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 |
} |
1013 |
|
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 |
} |
1040 |
|
1041 |
return 0; |
1042 |
} |
1043 |
|
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; |
1062 |
|
1063 |
int *convc = NULL; /* conv[m]: pre-computed convexities */ |
1064 |
double *areac = NULL; /* cumarea[m+1]: cache for fast area computation */ |
1065 |
|
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); |
1072 |
|
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 |
} |
1081 |
|
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 |
} |
1095 |
|
1096 |
pt[0] = -1; |
1097 |
pen[0] = 0; |
1098 |
len[0] = 0; |
1099 |
|
1100 |
/* Fixme: we always start from a fixed point -- should find the best
|
1101 |
curve cyclically ### */
|
1102 |
|
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; |
1108 |
|
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);
|
1129 |
|
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 |
} |
1155 |
|
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;
|
1162 |
|
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; |
1172 |
|
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 |
} |
1184 |
|
1185 |
/* ---------------------------------------------------------------------- */
|
1186 |
|
1187 |
#define TRY(x) if (x) goto try_error |
1188 |
|
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; |
1193 |
|
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 |
} |
1202 |
|
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); |
1217 |
|
1218 |
if (progress->callback) {
|
1219 |
cn += p->priv->len; |
1220 |
progress_update(cn/nn, progress); |
1221 |
} |
1222 |
} |
1223 |
|
1224 |
progress_update(1.0, progress); |
1225 |
|
1226 |
return 0; |
1227 |
|
1228 |
try_error:
|
1229 |
return 1; |
1230 |
} |