Remaining files transfered

1 files changed, 1464 insertions, 0 deletions

diff --git a/potrace/trace.cpp b/potrace/trace.cpp
new file mode 100644
index 0000000..fdf4afa
--- /dev/null
+++ b/potrace/trace.cpp
@@ -0,0 +1,1464 @@
+/* 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 <cmath>
+#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] );
+        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      seg_count = pp->curve.n;   // segment count in pp->curve
+    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      j, r;
+    opti_t   curve_prms;
+    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_MALLOC( pt, seg_count + 1, int );
+    SAFE_MALLOC( pen, seg_count + 1, double );
+    SAFE_MALLOC( len, seg_count + 1, int );
+    SAFE_MALLOC( opt, seg_count + 1, opti_t );
+    SAFE_MALLOC( convc, seg_count, int );
+    SAFE_MALLOC( areac, seg_count + 1, double );
+
+    /* pre-calculate convexity: +1 = right turn, -1 = left turn, 0 = corner */
+    for( int i = 0; i<seg_count; i++ )
+    {
+        if( pp->curve.tag[i] == POTRACE_CURVETO )
+        {
+            convc[i] =
+                sign( dpara( pp->curve.vertex[mod( i - 1, seg_count )],
+                             pp->curve.vertex[i],
+                             pp->curve.vertex[mod( i + 1, seg_count )] ) );
+        }
+        else
+        {
+            convc[i] = 0;
+        }
+    }
+
+    /* pre-calculate areas */
+    area     = 0.0;
+    areac[0] = 0.0;
+    p0 = pp->curve.vertex[0];
+    for( int i = 0; i<seg_count; i++ )
+    {
+        i1 = mod( i + 1, seg_count );
+        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;
+
+    // Avoid not initialized value for opt[j] (should not occur, but...)
+    for( j = 0; j<=seg_count; j++ )
+    {
+        opt[j].pen = 0.0;   // penalty
+        opt[j].c[0].x = opt[j].c[1].x = opt[j].c[0].y = opt[j].c[1].y = 0;
+        opt[j].t = opt[j].s = opt[j].alpha = 0.0; // curve parameters
+    }
+
+
+    /* Fixme: we always start from a fixed point -- should find the best
+     *  curve cyclically ### */
+
+    for( j = 1; j<=seg_count; j++ )
+    {
+        /* calculate best path from 0 to j */
+        pt[j]  = j - 1;
+        pen[j] = pen[j - 1];
+        len[j] = len[j - 1] + 1;
+
+        for( int i = j - 2; i>=0; i-- )
+        {
+            r = opti_penalty( pp, i, mod( j, seg_count ), &curve_prms, opttolerance, convc, areac );
+            if( r )
+            {
+                break;
+            }
+            if( len[j] > len[i] + 1 || (len[j] == len[i] + 1 && pen[j] > pen[i] + curve_prms.pen) )
+            {
+                pt[j]  = i;
+                pen[j] = pen[i] + curve_prms.pen;
+                len[j] = len[i] + 1;
+                opt[j] = curve_prms;
+            }
+        }
+    }
+
+    om = len[seg_count];
+    r  = privcurve_init( &pp->ocurve, om );
+    if( r )
+    {
+        goto malloc_error;
+    }
+    SAFE_MALLOC( s, om, double );
+    SAFE_MALLOC( t, om, double );
+
+    j = seg_count;
+    for( int i = om - 1; i>=0; i-- )
+    {
+        if( pt[j]==j - 1 )
+        {
+            pp->ocurve.tag[i]    = pp->curve.tag[mod( j, seg_count )];
+            pp->ocurve.c[i][0]   = pp->curve.c[mod( j, seg_count )][0];
+            pp->ocurve.c[i][1]   = pp->curve.c[mod( j, seg_count )][1];
+            pp->ocurve.c[i][2]   = pp->curve.c[mod( j, seg_count )][2];
+            pp->ocurve.vertex[i] = pp->curve.vertex[mod( j, seg_count )];
+            pp->ocurve.alpha[i]  = pp->curve.alpha[mod( j, seg_count )];
+            pp->ocurve.alpha0[i] = pp->curve.alpha0[mod( j, seg_count )];
+            pp->ocurve.beta[i]   = pp->curve.beta[mod( j, seg_count )];
+            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, seg_count )][2];
+            pp->ocurve.vertex[i] = interval( opt[j].s, pp->curve.c[mod( j, seg_count )][2],
+                                             pp->curve.vertex[mod( j, seg_count )] );
+            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( int 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;
+
+malloc_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 ) );
+        TRY( smooth( &p->priv->curve, p->sign, 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;
+}