Application: gvSIG desktop » gvSIG raster

/* Copyright (C) 2001-2007 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. */

/* $Id: trace.c 147 2007-04-09 00:44:09Z selinger $ */

/* transform jaggy paths into smooth curves */

#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_MALLOC(var, n, typ) \

  if ((var = (typ *)malloc((n)*sizeof(typ))) == NULL) goto malloc_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_MALLOC(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 + x*x;

    pp->sums[i+1].xy = pp->sums[i].xy + x*y;

    pp->sums[i+1].y2 = pp->sums[i].y2 + y*y;

  }

  return 0;

 malloc_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_MALLOC(pivk, n, int);

  SAFE_MALLOC(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_MALLOC(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;

 malloc_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;

  }

  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;

  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_MALLOC(pen, n+1, double);

  SAFE_MALLOC(prev, n+1, int);

  SAFE_MALLOC(clip0, n, int);

  SAFE_MALLOC(clip1, n+1, int);

  SAFE_MALLOC(seg0, n+1, int);

  SAFE_MALLOC(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 interations, 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_MALLOC(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;

 malloc_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_MALLOC(ctr, m, dpoint_t);

  SAFE_MALLOC(dir, m, dpoint_t);

  SAFE_MALLOC(q, m, quadform_t);

  r = privcurve_init(&pp->curve, m);

  if (r) {

    goto malloc_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;

 malloc_error:

  free(ctr);

  free(dir);

  free(q);

  return 1;

}

/* ---------------------------------------------------------------------- */

/* Stage 4: smoothing and corner analysis (Sec. 2.3.3) */

/* Always succeeds and returns 0 */

static int smooth(privcurve_t *curve, int sign, double alphamax) {

  int m = curve->n;

  int i, j, k;

  double dd, denom, alpha;

  dpoint_t p2, p3, p4;

  if (sign == '-') {

    /* reverse orientation of negative paths */

    for (i=0, j=m-1; i<j; i++, j--) {

      dpoint_t tmp;

      tmp = curve->vertex[i];

      curve->vertex[i] = curve->vertex[j];

      curve->vertex[j] = tmp;

    }

  }

  /* 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 0;

}

/* ---------------------------------------------------------------------- */

/* 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]);