SourceBrowser

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