summaryrefslogtreecommitdiff
path: root/Build/source/texk/dvisvgm/dvisvgm-src/libs/potrace/trace.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/texk/dvisvgm/dvisvgm-src/libs/potrace/trace.c')
-rw-r--r--Build/source/texk/dvisvgm/dvisvgm-src/libs/potrace/trace.c1249
1 files changed, 0 insertions, 1249 deletions
diff --git a/Build/source/texk/dvisvgm/dvisvgm-src/libs/potrace/trace.c b/Build/source/texk/dvisvgm/dvisvgm-src/libs/potrace/trace.c
deleted file mode 100644
index b830764447b..00000000000
--- a/Build/source/texk/dvisvgm/dvisvgm-src/libs/potrace/trace.c
+++ /dev/null
@@ -1,1249 +0,0 @@
-/* Copyright (C) 2001-2019 Peter Selinger.
- This file is part of Potrace. It is free software and it is covered
- by the GNU General Public License. See the file COPYING for details. */
-
-/* transform jaggy paths into smooth curves */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdio.h>
-#include <math.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include "potracelib.h"
-#include "curve.h"
-#include "lists.h"
-#include "auxiliary.h"
-#include "trace.h"
-#include "progress.h"
-
-#define INFTY 10000000 /* it suffices that this is longer than any
- path; it need not be really infinite */
-#define COS179 -0.999847695156 /* the cosine of 179 degrees */
-
-/* ---------------------------------------------------------------------- */
-#define SAFE_CALLOC(var, n, typ) \
- if ((var = (typ *)calloc(n, sizeof(typ))) == NULL) goto calloc_error
-
-/* ---------------------------------------------------------------------- */
-/* auxiliary functions */
-
-/* return a direction that is 90 degrees counterclockwise from p2-p0,
- but then restricted to one of the major wind directions (n, nw, w, etc) */
-static inline point_t dorth_infty(dpoint_t p0, dpoint_t p2) {
- point_t r;
-
- r.y = sign(p2.x-p0.x);
- r.x = -sign(p2.y-p0.y);
-
- return r;
-}
-
-/* return (p1-p0)x(p2-p0), the area of the parallelogram */
-static inline double dpara(dpoint_t p0, dpoint_t p1, dpoint_t p2) {
- double x1, y1, x2, y2;
-
- x1 = p1.x-p0.x;
- y1 = p1.y-p0.y;
- x2 = p2.x-p0.x;
- y2 = p2.y-p0.y;
-
- return x1*y2 - x2*y1;
-}
-
-/* ddenom/dpara have the property that the square of radius 1 centered
- at p1 intersects the line p0p2 iff |dpara(p0,p1,p2)| <= ddenom(p0,p2) */
-static inline double ddenom(dpoint_t p0, dpoint_t p2) {
- point_t r = dorth_infty(p0, p2);
-
- return r.y*(p2.x-p0.x) - r.x*(p2.y-p0.y);
-}
-
-/* return 1 if a <= b < c < a, in a cyclic sense (mod n) */
-static inline int cyclic(int a, int b, int c) {
- if (a<=c) {
- return (a<=b && b<c);
- } else {
- return (a<=b || b<c);
- }
-}
-
-/* determine the center and slope of the line i..j. Assume i<j. Needs
- "sum" components of p to be set. */
-static void pointslope(privpath_t *pp, int i, int j, dpoint_t *ctr, dpoint_t *dir) {
- /* assume i<j */
-
- int n = pp->len;
- sums_t *sums = pp->sums;
-
- double x, y, x2, xy, y2;
- double k;
- double a, b, c, lambda2, l;
- int r=0; /* rotations from i to j */
-
- while (j>=n) {
- j-=n;
- r+=1;
- }
- while (i>=n) {
- i-=n;
- r-=1;
- }
- while (j<0) {
- j+=n;
- r-=1;
- }
- while (i<0) {
- i+=n;
- r+=1;
- }
-
- x = sums[j+1].x-sums[i].x+r*sums[n].x;
- y = sums[j+1].y-sums[i].y+r*sums[n].y;
- x2 = sums[j+1].x2-sums[i].x2+r*sums[n].x2;
- xy = sums[j+1].xy-sums[i].xy+r*sums[n].xy;
- y2 = sums[j+1].y2-sums[i].y2+r*sums[n].y2;
- k = j+1-i+r*n;
-
- ctr->x = x/k;
- ctr->y = y/k;
-
- a = (x2-(double)x*x/k)/k;
- b = (xy-(double)x*y/k)/k;
- c = (y2-(double)y*y/k)/k;
-
- lambda2 = (a+c+sqrt((a-c)*(a-c)+4*b*b))/2; /* larger e.value */
-
- /* now find e.vector for lambda2 */
- a -= lambda2;
- c -= lambda2;
-
- if (fabs(a) >= fabs(c)) {
- l = sqrt(a*a+b*b);
- if (l!=0) {
- dir->x = -b/l;
- dir->y = a/l;
- }
- } else {
- l = sqrt(c*c+b*b);
- if (l!=0) {
- dir->x = -c/l;
- dir->y = b/l;
- }
- }
- if (l==0) {
- dir->x = dir->y = 0; /* sometimes this can happen when k=4:
- the two eigenvalues coincide */
- }
-}
-
-/* the type of (affine) quadratic forms, represented as symmetric 3x3
- matrices. The value of the quadratic form at a vector (x,y) is v^t
- Q v, where v = (x,y,1)^t. */
-typedef double quadform_t[3][3];
-
-/* Apply quadratic form Q to vector w = (w.x,w.y) */
-static inline double quadform(quadform_t Q, dpoint_t w) {
- double v[3];
- int i, j;
- double sum;
-
- v[0] = w.x;
- v[1] = w.y;
- v[2] = 1;
- sum = 0.0;
-
- for (i=0; i<3; i++) {
- for (j=0; j<3; j++) {
- sum += v[i] * Q[i][j] * v[j];
- }
- }
- return sum;
-}
-
-/* calculate p1 x p2 */
-static inline int xprod(point_t p1, point_t p2) {
- return p1.x*p2.y - p1.y*p2.x;
-}
-
-/* calculate (p1-p0)x(p3-p2) */
-static inline double cprod(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3) {
- double x1, y1, x2, y2;
-
- x1 = p1.x - p0.x;
- y1 = p1.y - p0.y;
- x2 = p3.x - p2.x;
- y2 = p3.y - p2.y;
-
- return x1*y2 - x2*y1;
-}
-
-/* calculate (p1-p0)*(p2-p0) */
-static inline double iprod(dpoint_t p0, dpoint_t p1, dpoint_t p2) {
- double x1, y1, x2, y2;
-
- x1 = p1.x - p0.x;
- y1 = p1.y - p0.y;
- x2 = p2.x - p0.x;
- y2 = p2.y - p0.y;
-
- return x1*x2 + y1*y2;
-}
-
-/* calculate (p1-p0)*(p3-p2) */
-static inline double iprod1(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3) {
- double x1, y1, x2, y2;
-
- x1 = p1.x - p0.x;
- y1 = p1.y - p0.y;
- x2 = p3.x - p2.x;
- y2 = p3.y - p2.y;
-
- return x1*x2 + y1*y2;
-}
-
-/* calculate distance between two points */
-static inline double ddist(dpoint_t p, dpoint_t q) {
- return sqrt(sq(p.x-q.x)+sq(p.y-q.y));
-}
-
-/* calculate point of a bezier curve */
-static inline dpoint_t bezier(double t, dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3) {
- double s = 1-t;
- dpoint_t res;
-
- /* Note: a good optimizing compiler (such as gcc-3) reduces the
- following to 16 multiplications, using common subexpression
- elimination. */
-
- res.x = s*s*s*p0.x + 3*(s*s*t)*p1.x + 3*(t*t*s)*p2.x + t*t*t*p3.x;
- res.y = s*s*s*p0.y + 3*(s*s*t)*p1.y + 3*(t*t*s)*p2.y + t*t*t*p3.y;
-
- return res;
-}
-
-/* calculate the point t in [0..1] on the (convex) bezier curve
- (p0,p1,p2,p3) which is tangent to q1-q0. Return -1.0 if there is no
- solution in [0..1]. */
-static double tangent(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3, dpoint_t q0, dpoint_t q1) {
- double A, B, C; /* (1-t)^2 A + 2(1-t)t B + t^2 C = 0 */
- double a, b, c; /* a t^2 + b t + c = 0 */
- double d, s, r1, r2;
-
- A = cprod(p0, p1, q0, q1);
- B = cprod(p1, p2, q0, q1);
- C = cprod(p2, p3, q0, q1);
-
- a = A - 2*B + C;
- b = -2*A + 2*B;
- c = A;
-
- d = b*b - 4*a*c;
-
- if (a==0 || d<0) {
- return -1.0;
- }
-
- s = sqrt(d);
-
- r1 = (-b + s) / (2 * a);
- r2 = (-b - s) / (2 * a);
-
- if (r1 >= 0 && r1 <= 1) {
- return r1;
- } else if (r2 >= 0 && r2 <= 1) {
- return r2;
- } else {
- return -1.0;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-/* Preparation: fill in the sum* fields of a path (used for later
- rapid summing). Return 0 on success, 1 with errno set on
- failure. */
-static int calc_sums(privpath_t *pp) {
- int i, x, y;
- int n = pp->len;
-
- SAFE_CALLOC(pp->sums, pp->len+1, sums_t);
-
- /* origin */
- pp->x0 = pp->pt[0].x;
- pp->y0 = pp->pt[0].y;
-
- /* preparatory computation for later fast summing */
- pp->sums[0].x2 = pp->sums[0].xy = pp->sums[0].y2 = pp->sums[0].x = pp->sums[0].y = 0;
- for (i=0; i<n; i++) {
- x = pp->pt[i].x - pp->x0;
- y = pp->pt[i].y - pp->y0;
- pp->sums[i+1].x = pp->sums[i].x + x;
- pp->sums[i+1].y = pp->sums[i].y + y;
- pp->sums[i+1].x2 = pp->sums[i].x2 + (double)x*x;
- pp->sums[i+1].xy = pp->sums[i].xy + (double)x*y;
- pp->sums[i+1].y2 = pp->sums[i].y2 + (double)y*y;
- }
- return 0;
-
- calloc_error:
- return 1;
-}
-
-/* ---------------------------------------------------------------------- */
-/* Stage 1: determine the straight subpaths (Sec. 2.2.1). Fill in the
- "lon" component of a path object (based on pt/len). For each i,
- lon[i] is the furthest index such that a straight line can be drawn
- from i to lon[i]. Return 1 on error with errno set, else 0. */
-
-/* this algorithm depends on the fact that the existence of straight
- subpaths is a triplewise property. I.e., there exists a straight
- line through squares i0,...,in iff there exists a straight line
- through i,j,k, for all i0<=i<j<k<=in. (Proof?) */
-
-/* this implementation of calc_lon is O(n^2). It replaces an older
- O(n^3) version. A "constraint" means that future points must
- satisfy xprod(constraint[0], cur) >= 0 and xprod(constraint[1],
- cur) <= 0. */
-
-/* Remark for Potrace 1.1: the current implementation of calc_lon is
- more complex than the implementation found in Potrace 1.0, but it
- is considerably faster. The introduction of the "nc" data structure
- means that we only have to test the constraints for "corner"
- points. On a typical input file, this speeds up the calc_lon
- function by a factor of 31.2, thereby decreasing its time share
- within the overall Potrace algorithm from 72.6% to 7.82%, and
- speeding up the overall algorithm by a factor of 3.36. On another
- input file, calc_lon was sped up by a factor of 6.7, decreasing its
- time share from 51.4% to 13.61%, and speeding up the overall
- algorithm by a factor of 1.78. In any case, the savings are
- substantial. */
-
-/* returns 0 on success, 1 on error with errno set */
-static int calc_lon(privpath_t *pp) {
- point_t *pt = pp->pt;
- int n = pp->len;
- int i, j, k, k1;
- int ct[4], dir;
- point_t constraint[2];
- point_t cur;
- point_t off;
- int *pivk = NULL; /* pivk[n] */
- int *nc = NULL; /* nc[n]: next corner */
- point_t dk; /* direction of k-k1 */
- int a, b, c, d;
-
- SAFE_CALLOC(pivk, n, int);
- SAFE_CALLOC(nc, n, int);
-
- /* initialize the nc data structure. Point from each point to the
- furthest future point to which it is connected by a vertical or
- horizontal segment. We take advantage of the fact that there is
- always a direction change at 0 (due to the path decomposition
- algorithm). But even if this were not so, there is no harm, as
- in practice, correctness does not depend on the word "furthest"
- above. */
- k = 0;
- for (i=n-1; i>=0; i--) {
- if (pt[i].x != pt[k].x && pt[i].y != pt[k].y) {
- k = i+1; /* necessarily i<n-1 in this case */
- }
- nc[i] = k;
- }
-
- SAFE_CALLOC(pp->lon, n, int);
-
- /* determine pivot points: for each i, let pivk[i] be the furthest k
- such that all j with i<j<k lie on a line connecting i,k. */
-
- for (i=n-1; i>=0; i--) {
- ct[0] = ct[1] = ct[2] = ct[3] = 0;
-
- /* keep track of "directions" that have occurred */
- dir = (3+3*(pt[mod(i+1,n)].x-pt[i].x)+(pt[mod(i+1,n)].y-pt[i].y))/2;
- ct[dir]++;
-
- constraint[0].x = 0;
- constraint[0].y = 0;
- constraint[1].x = 0;
- constraint[1].y = 0;
-
- /* find the next k such that no straight line from i to k */
- k = nc[i];
- k1 = i;
- while (1) {
-
- dir = (3+3*sign(pt[k].x-pt[k1].x)+sign(pt[k].y-pt[k1].y))/2;
- ct[dir]++;
-
- /* if all four "directions" have occurred, cut this path */
- if (ct[0] && ct[1] && ct[2] && ct[3]) {
- pivk[i] = k1;
- goto foundk;
- }
-
- cur.x = pt[k].x - pt[i].x;
- cur.y = pt[k].y - pt[i].y;
-
- /* see if current constraint is violated */
- if (xprod(constraint[0], cur) < 0 || xprod(constraint[1], cur) > 0) {
- goto constraint_viol;
- }
-
- /* else, update constraint */
- if (abs(cur.x) <= 1 && abs(cur.y) <= 1) {
- /* no constraint */
- } else {
- off.x = cur.x + ((cur.y>=0 && (cur.y>0 || cur.x<0)) ? 1 : -1);
- off.y = cur.y + ((cur.x<=0 && (cur.x<0 || cur.y<0)) ? 1 : -1);
- if (xprod(constraint[0], off) >= 0) {
- constraint[0] = off;
- }
- off.x = cur.x + ((cur.y<=0 && (cur.y<0 || cur.x<0)) ? 1 : -1);
- off.y = cur.y + ((cur.x>=0 && (cur.x>0 || cur.y<0)) ? 1 : -1);
- if (xprod(constraint[1], off) <= 0) {
- constraint[1] = off;
- }
- }
- k1 = k;
- k = nc[k1];
- if (!cyclic(k,i,k1)) {
- break;
- }
- }
- constraint_viol:
- /* k1 was the last "corner" satisfying the current constraint, and
- k is the first one violating it. We now need to find the last
- point along k1..k which satisfied the constraint. */
- dk.x = sign(pt[k].x-pt[k1].x);
- dk.y = sign(pt[k].y-pt[k1].y);
- cur.x = pt[k1].x - pt[i].x;
- cur.y = pt[k1].y - pt[i].y;
- /* find largest integer j such that xprod(constraint[0], cur+j*dk)
- >= 0 and xprod(constraint[1], cur+j*dk) <= 0. Use bilinearity
- of xprod. */
- a = xprod(constraint[0], cur);
- b = xprod(constraint[0], dk);
- c = xprod(constraint[1], cur);
- d = xprod(constraint[1], dk);
- /* find largest integer j such that a+j*b>=0 and c+j*d<=0. This
- can be solved with integer arithmetic. */
- j = INFTY;
- if (b<0) {
- j = floordiv(a,-b);
- }
- if (d>0) {
- j = min(j, floordiv(-c,d));
- }
- pivk[i] = mod(k1+j,n);
- foundk:
- ;
- } /* for i */
-
- /* clean up: for each i, let lon[i] be the largest k such that for
- all i' with i<=i'<k, i'<k<=pivk[i']. */
-
- j=pivk[n-1];
- pp->lon[n-1]=j;
- for (i=n-2; i>=0; i--) {
- if (cyclic(i+1,pivk[i],j)) {
- j=pivk[i];
- }
- pp->lon[i]=j;
- }
-
- for (i=n-1; cyclic(mod(i+1,n),j,pp->lon[i]); i--) {
- pp->lon[i] = j;
- }
-
- free(pivk);
- free(nc);
- return 0;
-
- calloc_error:
- free(pivk);
- free(nc);
- return 1;
-}
-
-
-/* ---------------------------------------------------------------------- */
-/* Stage 2: calculate the optimal polygon (Sec. 2.2.2-2.2.4). */
-
-/* Auxiliary function: calculate the penalty of an edge from i to j in
- the given path. This needs the "lon" and "sum*" data. */
-
-static double penalty3(privpath_t *pp, int i, int j) {
- int n = pp->len;
- point_t *pt = pp->pt;
- sums_t *sums = pp->sums;
-
- /* assume 0<=i<j<=n */
- double x, y, x2, xy, y2;
- double k;
- double a, b, c, s;
- double px, py, ex, ey;
-
- int r = 0; /* rotations from i to j */
-
- if (j>=n) {
- j -= n;
- r = 1;
- }
-
- /* critical inner loop: the "if" gives a 4.6 percent speedup */
- if (r == 0) {
- x = sums[j+1].x - sums[i].x;
- y = sums[j+1].y - sums[i].y;
- x2 = sums[j+1].x2 - sums[i].x2;
- xy = sums[j+1].xy - sums[i].xy;
- y2 = sums[j+1].y2 - sums[i].y2;
- k = j+1 - i;
- } else {
- x = sums[j+1].x - sums[i].x + sums[n].x;
- y = sums[j+1].y - sums[i].y + sums[n].y;
- x2 = sums[j+1].x2 - sums[i].x2 + sums[n].x2;
- xy = sums[j+1].xy - sums[i].xy + sums[n].xy;
- y2 = sums[j+1].y2 - sums[i].y2 + sums[n].y2;
- k = j+1 - i + n;
- }
-
- px = (pt[i].x + pt[j].x) / 2.0 - pt[0].x;
- py = (pt[i].y + pt[j].y) / 2.0 - pt[0].y;
- ey = (pt[j].x - pt[i].x);
- ex = -(pt[j].y - pt[i].y);
-
- a = ((x2 - 2*x*px) / k + px*px);
- b = ((xy - x*py - y*px) / k + px*py);
- c = ((y2 - 2*y*py) / k + py*py);
-
- s = ex*ex*a + 2*ex*ey*b + ey*ey*c;
-
- return sqrt(s);
-}
-
-/* find the optimal polygon. Fill in the m and po components. Return 1
- on failure with errno set, else 0. Non-cyclic version: assumes i=0
- is in the polygon. Fixme: implement cyclic version. */
-static int bestpolygon(privpath_t *pp)
-{
- int i, j, m, k;
- int n = pp->len;
- double *pen = NULL; /* pen[n+1]: penalty vector */
- int *prev = NULL; /* prev[n+1]: best path pointer vector */
- int *clip0 = NULL; /* clip0[n]: longest segment pointer, non-cyclic */
- int *clip1 = NULL; /* clip1[n+1]: backwards segment pointer, non-cyclic */
- int *seg0 = NULL; /* seg0[m+1]: forward segment bounds, m<=n */
- int *seg1 = NULL; /* seg1[m+1]: backward segment bounds, m<=n */
- double thispen;
- double best;
- int c;
-
- SAFE_CALLOC(pen, n+1, double);
- SAFE_CALLOC(prev, n+1, int);
- SAFE_CALLOC(clip0, n, int);
- SAFE_CALLOC(clip1, n+1, int);
- SAFE_CALLOC(seg0, n+1, int);
- SAFE_CALLOC(seg1, n+1, int);
-
- /* calculate clipped paths */
- for (i=0; i<n; i++) {
- c = mod(pp->lon[mod(i-1,n)]-1,n);
- if (c == i) {
- c = mod(i+1,n);
- }
- if (c < i) {
- clip0[i] = n;
- } else {
- clip0[i] = c;
- }
- }
-
- /* calculate backwards path clipping, non-cyclic. j <= clip0[i] iff
- clip1[j] <= i, for i,j=0..n. */
- j = 1;
- for (i=0; i<n; i++) {
- while (j <= clip0[i]) {
- clip1[j] = i;
- j++;
- }
- }
-
- /* calculate seg0[j] = longest path from 0 with j segments */
- i = 0;
- for (j=0; i<n; j++) {
- seg0[j] = i;
- i = clip0[i];
- }
- seg0[j] = n;
- m = j;
-
- /* calculate seg1[j] = longest path to n with m-j segments */
- i = n;
- for (j=m; j>0; j--) {
- seg1[j] = i;
- i = clip1[i];
- }
- seg1[0] = 0;
-
- /* now find the shortest path with m segments, based on penalty3 */
- /* note: the outer 2 loops jointly have at most n iterations, thus
- the worst-case behavior here is quadratic. In practice, it is
- close to linear since the inner loop tends to be short. */
- pen[0]=0;
- for (j=1; j<=m; j++) {
- for (i=seg1[j]; i<=seg0[j]; i++) {
- best = -1;
- for (k=seg0[j-1]; k>=clip1[i]; k--) {
- thispen = penalty3(pp, k, i) + pen[k];
- if (best < 0 || thispen < best) {
- prev[i] = k;
- best = thispen;
- }
- }
- pen[i] = best;
- }
- }
-
- pp->m = m;
- SAFE_CALLOC(pp->po, m, int);
-
- /* read off shortest path */
- for (i=n, j=m-1; i>0; j--) {
- i = prev[i];
- pp->po[j] = i;
- }
-
- free(pen);
- free(prev);
- free(clip0);
- free(clip1);
- free(seg0);
- free(seg1);
- return 0;
-
- calloc_error:
- free(pen);
- free(prev);
- free(clip0);
- free(clip1);
- free(seg0);
- free(seg1);
- return 1;
-}
-
-/* ---------------------------------------------------------------------- */
-/* Stage 3: vertex adjustment (Sec. 2.3.1). */
-
-/* Adjust vertices of optimal polygon: calculate the intersection of
- the two "optimal" line segments, then move it into the unit square
- if it lies outside. Return 1 with errno set on error; 0 on
- success. */
-
-static int adjust_vertices(privpath_t *pp) {
- int m = pp->m;
- int *po = pp->po;
- int n = pp->len;
- point_t *pt = pp->pt;
- int x0 = pp->x0;
- int y0 = pp->y0;
-
- dpoint_t *ctr = NULL; /* ctr[m] */
- dpoint_t *dir = NULL; /* dir[m] */
- quadform_t *q = NULL; /* q[m] */
- double v[3];
- double d;
- int i, j, k, l;
- dpoint_t s;
- int r;
-
- SAFE_CALLOC(ctr, m, dpoint_t);
- SAFE_CALLOC(dir, m, dpoint_t);
- SAFE_CALLOC(q, m, quadform_t);
-
- r = privcurve_init(&pp->curve, m);
- if (r) {
- goto calloc_error;
- }
-
- /* calculate "optimal" point-slope representation for each line
- segment */
- for (i=0; i<m; i++) {
- j = po[mod(i+1,m)];
- j = mod(j-po[i],n)+po[i];
- pointslope(pp, po[i], j, &ctr[i], &dir[i]);
- }
-
- /* represent each line segment as a singular quadratic form; the
- distance of a point (x,y) from the line segment will be
- (x,y,1)Q(x,y,1)^t, where Q=q[i]. */
- for (i=0; i<m; i++) {
- d = sq(dir[i].x) + sq(dir[i].y);
- if (d == 0.0) {
- for (j=0; j<3; j++) {
- for (k=0; k<3; k++) {
- q[i][j][k] = 0;
- }
- }
- } else {
- v[0] = dir[i].y;
- v[1] = -dir[i].x;
- v[2] = - v[1] * ctr[i].y - v[0] * ctr[i].x;
- for (l=0; l<3; l++) {
- for (k=0; k<3; k++) {
- q[i][l][k] = v[l] * v[k] / d;
- }
- }
- }
- }
-
- /* now calculate the "intersections" of consecutive segments.
- Instead of using the actual intersection, we find the point
- within a given unit square which minimizes the square distance to
- the two lines. */
- for (i=0; i<m; i++) {
- quadform_t Q;
- dpoint_t w;
- double dx, dy;
- double det;
- double min, cand; /* minimum and candidate for minimum of quad. form */
- double xmin, ymin; /* coordinates of minimum */
- int z;
-
- /* let s be the vertex, in coordinates relative to x0/y0 */
- s.x = pt[po[i]].x-x0;
- s.y = pt[po[i]].y-y0;
-
- /* intersect segments i-1 and i */
-
- j = mod(i-1,m);
-
- /* add quadratic forms */
- for (l=0; l<3; l++) {
- for (k=0; k<3; k++) {
- Q[l][k] = q[j][l][k] + q[i][l][k];
- }
- }
-
- while(1) {
- /* minimize the quadratic form Q on the unit square */
- /* find intersection */
-
-#ifdef HAVE_GCC_LOOP_BUG
- /* work around gcc bug #12243 */
- free(NULL);
-#endif
-
- det = Q[0][0]*Q[1][1] - Q[0][1]*Q[1][0];
- if (det != 0.0) {
- w.x = (-Q[0][2]*Q[1][1] + Q[1][2]*Q[0][1]) / det;
- w.y = ( Q[0][2]*Q[1][0] - Q[1][2]*Q[0][0]) / det;
- break;
- }
-
- /* matrix is singular - lines are parallel. Add another,
- orthogonal axis, through the center of the unit square */
- if (Q[0][0]>Q[1][1]) {
- v[0] = -Q[0][1];
- v[1] = Q[0][0];
- } else if (Q[1][1]) {
- v[0] = -Q[1][1];
- v[1] = Q[1][0];
- } else {
- v[0] = 1;
- v[1] = 0;
- }
- d = sq(v[0]) + sq(v[1]);
- v[2] = - v[1] * s.y - v[0] * s.x;
- for (l=0; l<3; l++) {
- for (k=0; k<3; k++) {
- Q[l][k] += v[l] * v[k] / d;
- }
- }
- }
- dx = fabs(w.x-s.x);
- dy = fabs(w.y-s.y);
- if (dx <= .5 && dy <= .5) {
- pp->curve.vertex[i].x = w.x+x0;
- pp->curve.vertex[i].y = w.y+y0;
- continue;
- }
-
- /* the minimum was not in the unit square; now minimize quadratic
- on boundary of square */
- min = quadform(Q, s);
- xmin = s.x;
- ymin = s.y;
-
- if (Q[0][0] == 0.0) {
- goto fixx;
- }
- for (z=0; z<2; z++) { /* value of the y-coordinate */
- w.y = s.y-0.5+z;
- w.x = - (Q[0][1] * w.y + Q[0][2]) / Q[0][0];
- dx = fabs(w.x-s.x);
- cand = quadform(Q, w);
- if (dx <= .5 && cand < min) {
- min = cand;
- xmin = w.x;
- ymin = w.y;
- }
- }
- fixx:
- if (Q[1][1] == 0.0) {
- goto corners;
- }
- for (z=0; z<2; z++) { /* value of the x-coordinate */
- w.x = s.x-0.5+z;
- w.y = - (Q[1][0] * w.x + Q[1][2]) / Q[1][1];
- dy = fabs(w.y-s.y);
- cand = quadform(Q, w);
- if (dy <= .5 && cand < min) {
- min = cand;
- xmin = w.x;
- ymin = w.y;
- }
- }
- corners:
- /* check four corners */
- for (l=0; l<2; l++) {
- for (k=0; k<2; k++) {
- w.x = s.x-0.5+l;
- w.y = s.y-0.5+k;
- cand = quadform(Q, w);
- if (cand < min) {
- min = cand;
- xmin = w.x;
- ymin = w.y;
- }
- }
- }
-
- pp->curve.vertex[i].x = xmin + x0;
- pp->curve.vertex[i].y = ymin + y0;
- continue;
- }
-
- free(ctr);
- free(dir);
- free(q);
- return 0;
-
- calloc_error:
- free(ctr);
- free(dir);
- free(q);
- return 1;
-}
-
-/* ---------------------------------------------------------------------- */
-/* Stage 4: smoothing and corner analysis (Sec. 2.3.3) */
-
-/* reverse orientation of a path */
-static void reverse(privcurve_t *curve) {
- int m = curve->n;
- int i, j;
- dpoint_t tmp;
-
- for (i=0, j=m-1; i<j; i++, j--) {
- tmp = curve->vertex[i];
- curve->vertex[i] = curve->vertex[j];
- curve->vertex[j] = tmp;
- }
-}
-
-/* Always succeeds */
-static void smooth(privcurve_t *curve, double alphamax) {
- int m = curve->n;
-
- int i, j, k;
- double dd, denom, alpha;
- dpoint_t p2, p3, p4;
-
- /* examine each vertex and find its best fit */
- for (i=0; i<m; i++) {
- j = mod(i+1, m);
- k = mod(i+2, m);
- p4 = interval(1/2.0, curve->vertex[k], curve->vertex[j]);
-
- denom = ddenom(curve->vertex[i], curve->vertex[k]);
- if (denom != 0.0) {
- dd = dpara(curve->vertex[i], curve->vertex[j], curve->vertex[k]) / denom;
- dd = fabs(dd);
- alpha = dd>1 ? (1 - 1.0/dd) : 0;
- alpha = alpha / 0.75;
- } else {
- alpha = 4/3.0;
- }
- curve->alpha0[j] = alpha; /* remember "original" value of alpha */
-
- if (alpha >= alphamax) { /* pointed corner */
- curve->tag[j] = POTRACE_CORNER;
- curve->c[j][1] = curve->vertex[j];
- curve->c[j][2] = p4;
- } else {
- if (alpha < 0.55) {
- alpha = 0.55;
- } else if (alpha > 1) {
- alpha = 1;
- }
- p2 = interval(.5+.5*alpha, curve->vertex[i], curve->vertex[j]);
- p3 = interval(.5+.5*alpha, curve->vertex[k], curve->vertex[j]);
- curve->tag[j] = POTRACE_CURVETO;
- curve->c[j][0] = p2;
- curve->c[j][1] = p3;
- curve->c[j][2] = p4;
- }
- curve->alpha[j] = alpha; /* store the "cropped" value of alpha */
- curve->beta[j] = 0.5;
- }
- curve->alphacurve = 1;
-
- return;
-}
-
-/* ---------------------------------------------------------------------- */
-/* Stage 5: Curve optimization (Sec. 2.4) */
-
-/* a private type for the result of opti_penalty */
-struct opti_s {
- double pen; /* penalty */
- dpoint_t c[2]; /* curve parameters */
- double t, s; /* curve parameters */
- double alpha; /* curve parameter */
-};
-typedef struct opti_s opti_t;
-
-/* calculate best fit from i+.5 to j+.5. Assume i<j (cyclically).
- Return 0 and set badness and parameters (alpha, beta), if
- possible. Return 1 if impossible. */
-static int opti_penalty(privpath_t *pp, int i, int j, opti_t *res, double opttolerance, int *convc, double *areac) {
- int m = pp->curve.n;
- int k, k1, k2, conv, i1;
- double area, alpha, d, d1, d2;
- dpoint_t p0, p1, p2, p3, pt;
- double A, R, A1, A2, A3, A4;
- double s, t;
-
- /* check convexity, corner-freeness, and maximum bend < 179 degrees */
-
- if (i==j) { /* sanity - a full loop can never be an opticurve */
- return 1;
- }
-
- k = i;
- i1 = mod(i+1, m);
- k1 = mod(k+1, m);
- conv = convc[k1];
- if (conv == 0) {
- return 1;
- }
- d = ddist(pp->curve.vertex[i], pp->curve.vertex[i1]);
- for (k=k1; k!=j; k=k1) {
- k1 = mod(k+1, m);
- k2 = mod(k+2, m);
- if (convc[k1] != conv) {
- return 1;
- }
- if (sign(cprod(pp->curve.vertex[i], pp->curve.vertex[i1], pp->curve.vertex[k1], pp->curve.vertex[k2])) != conv) {
- return 1;
- }
- 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) {
- return 1;
- }
- }
-
- /* the curve we're working in: */
- p0 = pp->curve.c[mod(i,m)][2];
- p1 = pp->curve.vertex[mod(i+1,m)];
- p2 = pp->curve.vertex[mod(j,m)];
- p3 = pp->curve.c[mod(j,m)][2];
-
- /* determine its area */
- area = areac[j] - areac[i];
- area -= dpara(pp->curve.vertex[0], pp->curve.c[i][2], pp->curve.c[j][2])/2;
- if (i>=j) {
- area += areac[m];
- }
-
- /* find intersection o of p0p1 and p2p3. Let t,s such that o =
- interval(t,p0,p1) = interval(s,p3,p2). Let A be the area of the
- triangle (p0,o,p3). */
-
- A1 = dpara(p0, p1, p2);
- A2 = dpara(p0, p1, p3);
- A3 = dpara(p0, p2, p3);
- /* A4 = dpara(p1, p2, p3); */
- A4 = A1+A3-A2;
-
- if (A2 == A1) { /* this should never happen */
- return 1;
- }
-
- t = A3/(A3-A4);
- s = A2/(A2-A1);
- A = A2 * t / 2.0;
-
- if (A == 0.0) { /* this should never happen */
- return 1;
- }
-
- R = area / A; /* relative area */
- alpha = 2 - sqrt(4 - R / 0.3); /* overall alpha for p0-o-p3 curve */
-
- res->c[0] = interval(t * alpha, p0, p1);
- res->c[1] = interval(s * alpha, p3, p2);
- res->alpha = alpha;
- res->t = t;
- res->s = s;
-
- p1 = res->c[0];
- p2 = res->c[1]; /* the proposed curve is now (p0,p1,p2,p3) */
-
- res->pen = 0;
-
- /* calculate penalty */
- /* check tangency with edges */
- for (k=mod(i+1,m); k!=j; k=k1) {
- k1 = mod(k+1,m);
- t = tangent(p0, p1, p2, p3, pp->curve.vertex[k], pp->curve.vertex[k1]);
- if (t<-.5) {
- return 1;
- }
- pt = bezier(t, p0, p1, p2, p3);
- d = ddist(pp->curve.vertex[k], pp->curve.vertex[k1]);
- if (d == 0.0) { /* this should never happen */
- return 1;
- }
- d1 = dpara(pp->curve.vertex[k], pp->curve.vertex[k1], pt) / d;
- if (fabs(d1) > opttolerance) {
- return 1;
- }
- if (iprod(pp->curve.vertex[k], pp->curve.vertex[k1], pt) < 0 || iprod(pp->curve.vertex[k1], pp->curve.vertex[k], pt) < 0) {
- return 1;
- }
- res->pen += sq(d1);
- }
-
- /* check corners */
- for (k=i; k!=j; k=k1) {
- k1 = mod(k+1,m);
- t = tangent(p0, p1, p2, p3, pp->curve.c[k][2], pp->curve.c[k1][2]);
- if (t<-.5) {
- return 1;
- }
- pt = bezier(t, p0, p1, p2, p3);
- d = ddist(pp->curve.c[k][2], pp->curve.c[k1][2]);
- if (d == 0.0) { /* this should never happen */
- return 1;
- }
- d1 = dpara(pp->curve.c[k][2], pp->curve.c[k1][2], pt) / d;
- d2 = dpara(pp->curve.c[k][2], pp->curve.c[k1][2], pp->curve.vertex[k1]) / d;
- d2 *= 0.75 * pp->curve.alpha[k1];
- if (d2 < 0) {
- d1 = -d1;
- d2 = -d2;
- }
- if (d1 < d2 - opttolerance) {
- return 1;
- }
- if (d1 < d2) {
- res->pen += sq(d1 - d2);
- }
- }
-
- return 0;
-}
-
-/* optimize the path p, replacing sequences of Bezier segments by a
- single segment when possible. Return 0 on success, 1 with errno set
- on failure. */
-static int opticurve(privpath_t *pp, double opttolerance) {
- int m = pp->curve.n;
- int *pt = NULL; /* pt[m+1] */
- double *pen = NULL; /* pen[m+1] */
- int *len = NULL; /* len[m+1] */
- opti_t *opt = NULL; /* opt[m+1] */
- int om;
- int i,j,r;
- opti_t o;
- dpoint_t p0;
- int i1;
- double area;
- double alpha;
- double *s = NULL;
- double *t = NULL;
-
- int *convc = NULL; /* conv[m]: pre-computed convexities */
- double *areac = NULL; /* cumarea[m+1]: cache for fast area computation */
-
- SAFE_CALLOC(pt, m+1, int);
- SAFE_CALLOC(pen, m+1, double);
- SAFE_CALLOC(len, m+1, int);
- SAFE_CALLOC(opt, m+1, opti_t);
- SAFE_CALLOC(convc, m, int);
- SAFE_CALLOC(areac, m+1, double);
-
- /* pre-calculate convexity: +1 = right turn, -1 = left turn, 0 = corner */
- for (i=0; i<m; i++) {
- if (pp->curve.tag[i] == POTRACE_CURVETO) {
- convc[i] = sign(dpara(pp->curve.vertex[mod(i-1,m)], pp->curve.vertex[i], pp->curve.vertex[mod(i+1,m)]));
- } else {
- convc[i] = 0;
- }
- }
-
- /* pre-calculate areas */
- area = 0.0;
- areac[0] = 0.0;
- p0 = pp->curve.vertex[0];
- for (i=0; i<m; i++) {
- i1 = mod(i+1, m);
- if (pp->curve.tag[i1] == POTRACE_CURVETO) {
- alpha = pp->curve.alpha[i1];
- area += 0.3*alpha*(4-alpha)*dpara(pp->curve.c[i][2], pp->curve.vertex[i1], pp->curve.c[i1][2])/2;
- area += dpara(p0, pp->curve.c[i][2], pp->curve.c[i1][2])/2;
- }
- areac[i+1] = area;
- }
-
- pt[0] = -1;
- pen[0] = 0;
- len[0] = 0;
-
- /* Fixme: we always start from a fixed point -- should find the best
- curve cyclically */
-
- for (j=1; j<=m; j++) {
- /* calculate best path from 0 to j */
- pt[j] = j-1;
- pen[j] = pen[j-1];
- len[j] = len[j-1]+1;
-
- for (i=j-2; i>=0; i--) {
- r = opti_penalty(pp, i, mod(j,m), &o, opttolerance, convc, areac);
- if (r) {
- break;
- }
- if (len[j] > len[i]+1 || (len[j] == len[i]+1 && pen[j] > pen[i] + o.pen)) {
- pt[j] = i;
- pen[j] = pen[i] + o.pen;
- len[j] = len[i] + 1;
- opt[j] = o;
- }
- }
- }
- om = len[m];
- r = privcurve_init(&pp->ocurve, om);
- if (r) {
- goto calloc_error;
- }
- SAFE_CALLOC(s, om, double);
- SAFE_CALLOC(t, om, double);
-
- j = m;
- for (i=om-1; i>=0; i--) {
- if (pt[j]==j-1) {
- pp->ocurve.tag[i] = pp->curve.tag[mod(j,m)];
- pp->ocurve.c[i][0] = pp->curve.c[mod(j,m)][0];
- pp->ocurve.c[i][1] = pp->curve.c[mod(j,m)][1];
- pp->ocurve.c[i][2] = pp->curve.c[mod(j,m)][2];
- pp->ocurve.vertex[i] = pp->curve.vertex[mod(j,m)];
- pp->ocurve.alpha[i] = pp->curve.alpha[mod(j,m)];
- pp->ocurve.alpha0[i] = pp->curve.alpha0[mod(j,m)];
- pp->ocurve.beta[i] = pp->curve.beta[mod(j,m)];
- s[i] = t[i] = 1.0;
- } else {
- pp->ocurve.tag[i] = POTRACE_CURVETO;
- pp->ocurve.c[i][0] = opt[j].c[0];
- pp->ocurve.c[i][1] = opt[j].c[1];
- pp->ocurve.c[i][2] = pp->curve.c[mod(j,m)][2];
- pp->ocurve.vertex[i] = interval(opt[j].s, pp->curve.c[mod(j,m)][2], pp->curve.vertex[mod(j,m)]);
- pp->ocurve.alpha[i] = opt[j].alpha;
- pp->ocurve.alpha0[i] = opt[j].alpha;
- s[i] = opt[j].s;
- t[i] = opt[j].t;
- }
- j = pt[j];
- }
-
- /* calculate beta parameters */
- for (i=0; i<om; i++) {
- i1 = mod(i+1,om);
- pp->ocurve.beta[i] = s[i] / (s[i] + t[i1]);
- }
- pp->ocurve.alphacurve = 1;
-
- free(pt);
- free(pen);
- free(len);
- free(opt);
- free(s);
- free(t);
- free(convc);
- free(areac);
- return 0;
-
- calloc_error:
- free(pt);
- free(pen);
- free(len);
- free(opt);
- free(s);
- free(t);
- free(convc);
- free(areac);
- return 1;
-}
-
-/* ---------------------------------------------------------------------- */
-
-#define TRY(x) if (x) goto try_error
-
-/* return 0 on success, 1 on error with errno set. */
-int process_path(path_t *plist, const potrace_param_t *param, progress_t *progress) {
- path_t *p;
- double nn = 0, cn = 0;
-
- if (progress->callback) {
- /* precompute task size for progress estimates */
- nn = 0;
- list_forall (p, plist) {
- nn += p->priv->len;
- }
- cn = 0;
- }
-
- /* call downstream function with each path */
- list_forall (p, plist) {
- TRY(calc_sums(p->priv));
- TRY(calc_lon(p->priv));
- TRY(bestpolygon(p->priv));
- TRY(adjust_vertices(p->priv));
- if (p->sign == '-') { /* reverse orientation of negative paths */
- reverse(&p->priv->curve);
- }
- smooth(&p->priv->curve, param->alphamax);
- if (param->opticurve) {
- TRY(opticurve(p->priv, param->opttolerance));
- p->priv->fcurve = &p->priv->ocurve;
- } else {
- p->priv->fcurve = &p->priv->curve;
- }
- privcurve_to_curve(p->priv->fcurve, &p->curve);
-
- if (progress->callback) {
- cn += p->priv->len;
- progress_update(cn/nn, progress);
- }
- }
-
- progress_update(1.0, progress);
-
- return 0;
-
- try_error:
- return 1;
-}