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