/* Copyright (C) 2001-2015 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 + 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;

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


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