KiCad PCB EDA Suite
trace.cpp File Reference
#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"

Go to the source code of this file.

Classes

struct  opti_s
 

Macros

#define INFTY
 
#define COS179   -0.999847695156 /* the cosine of 179 degrees */
 
#define SAFE_CALLOC(var, n, typ)
 
#define TRY(x)
 

Typedefs

typedef double quadform_t[3][3]
 
typedef struct opti_s opti_t
 

Functions

static point_t dorth_infty (dpoint_t p0, dpoint_t p2)
 
static double dpara (dpoint_t p0, dpoint_t p1, dpoint_t p2)
 
static double ddenom (dpoint_t p0, dpoint_t p2)
 
static int cyclic (int a, int b, int c)
 
static void pointslope (privpath_t *pp, int i, int j, dpoint_t *ctr, dpoint_t *dir)
 
static double quadform (quadform_t Q, dpoint_t w)
 
static int xprod (point_t p1, point_t p2)
 
static double cprod (dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3)
 
static double iprod (dpoint_t p0, dpoint_t p1, dpoint_t p2)
 
static double iprod1 (dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3)
 
static double ddist (dpoint_t p, dpoint_t q)
 
static dpoint_t bezier (double t, dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3)
 
static double tangent (dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3, dpoint_t q0, dpoint_t q1)
 
static int calc_sums (privpath_t *pp)
 
static int calc_lon (privpath_t *pp)
 
static double penalty3 (privpath_t *pp, int i, int j)
 
static int bestpolygon (privpath_t *pp)
 
static int adjust_vertices (privpath_t *pp)
 
static void reverse (privcurve_t *curve)
 
static void smooth (privcurve_t *curve, double alphamax)
 
static int opti_penalty (privpath_t *pp, int i, int j, opti_t *res, double opttolerance, int *convc, double *areac)
 
static int opticurve (privpath_t *pp, double opttolerance)
 
int process_path (path_t *plist, const potrace_param_t *param, progress_t *progress)
 

Macro Definition Documentation

#define COS179   -0.999847695156 /* the cosine of 179 degrees */

Definition at line 25 of file trace.cpp.

Referenced by opti_penalty().

#define INFTY
Value:
10000000 /* it suffices that this is longer than any
* path; it need not be really infinite */

Definition at line 23 of file trace.cpp.

Referenced by calc_lon().

#define SAFE_CALLOC (   var,
  n,
  typ 
)
Value:
if( ( var = (typ*) calloc( n, sizeof(typ) ) ) == NULL ) \
goto calloc_error

Definition at line 28 of file trace.cpp.

Referenced by adjust_vertices(), bestpolygon(), calc_lon(), calc_sums(), and opticurve().

#define TRY (   x)
Value:
if( x ) \
goto try_error

Definition at line 1485 of file trace.cpp.

Referenced by process_path().

Typedef Documentation

typedef struct opti_s opti_t

Definition at line 1111 of file trace.cpp.

typedef double quadform_t[3][3]

Definition at line 176 of file trace.cpp.

Function Documentation

static int adjust_vertices ( privpath_t pp)
static

Definition at line 769 of file trace.cpp.

References potrace_privpath_s::curve, potrace_privpath_s::len, potrace_privpath_s::m, min, mod, potrace_privpath_s::po, pointslope(), privcurve_init(), potrace_privpath_s::pt, quadform(), SAFE_CALLOC, sq, privcurve_s::vertex, point_s::x, potrace_dpoint_s::x, potrace_privpath_s::x0, point_s::y, potrace_dpoint_s::y, and potrace_privpath_s::y0.

Referenced by process_path().

770 {
771  int m = pp->m;
772  int* po = pp->po;
773  int n = pp->len;
774  point_t* pt = pp->pt;
775  int x0 = pp->x0;
776  int y0 = pp->y0;
777 
778  dpoint_t* ctr = NULL; /* ctr[m] */
779  dpoint_t* dir = NULL; /* dir[m] */
780  quadform_t* q = NULL; /* q[m] */
781  double v[3];
782  double d;
783  int i, j, k, l;
784  dpoint_t s;
785  int r;
786 
787  SAFE_CALLOC( ctr, m, dpoint_t );
788  SAFE_CALLOC( dir, m, dpoint_t );
789  SAFE_CALLOC( q, m, quadform_t );
790 
791  r = privcurve_init( &pp->curve, m );
792 
793  if( r )
794  {
795  goto calloc_error;
796  }
797 
798  /* calculate "optimal" point-slope representation for each line
799  * segment */
800  for( i = 0; i<m; i++ )
801  {
802  j = po[mod( i + 1, m )];
803  j = mod( j - po[i], n ) + po[i];
804  pointslope( pp, po[i], j, &ctr[i], &dir[i] );
805  }
806 
807  /* represent each line segment as a singular quadratic form; the
808  * distance of a point (x,y) from the line segment will be
809  * (x,y,1)Q(x,y,1)^t, where Q=q[i]. */
810  for( i = 0; i<m; i++ )
811  {
812  d = sq( dir[i].x ) + sq( dir[i].y );
813 
814  if( d == 0.0 )
815  {
816  for( j = 0; j<3; j++ )
817  {
818  for( k = 0; k<3; k++ )
819  {
820  q[i][j][k] = 0;
821  }
822  }
823  }
824  else
825  {
826  v[0] = dir[i].y;
827  v[1] = -dir[i].x;
828  v[2] = -v[1] * ctr[i].y - v[0] * ctr[i].x;
829 
830  for( l = 0; l<3; l++ )
831  {
832  for( k = 0; k<3; k++ )
833  {
834  q[i][l][k] = v[l] * v[k] / d;
835  }
836  }
837  }
838  }
839 
840  /* now calculate the "intersections" of consecutive segments.
841  * Instead of using the actual intersection, we find the point
842  * within a given unit square which minimizes the square distance to
843  * the two lines. */
844  for( i = 0; i<m; i++ )
845  {
846  quadform_t Q;
847  dpoint_t w;
848  double dx, dy;
849  double det;
850  double min, cand; /* minimum and candidate for minimum of quad. form */
851  double xmin, ymin; /* coordinates of minimum */
852  int z;
853 
854  /* let s be the vertex, in coordinates relative to x0/y0 */
855  s.x = pt[po[i]].x - x0;
856  s.y = pt[po[i]].y - y0;
857 
858  /* intersect segments i-1 and i */
859 
860  j = mod( i - 1, m );
861 
862  /* add quadratic forms */
863  for( l = 0; l<3; l++ )
864  {
865  for( k = 0; k<3; k++ )
866  {
867  Q[l][k] = q[j][l][k] + q[i][l][k];
868  }
869  }
870 
871  while( 1 )
872  {
873  /* minimize the quadratic form Q on the unit square */
874  /* find intersection */
875 
876 #ifdef HAVE_GCC_LOOP_BUG
877  /* work around gcc bug #12243 */
878  free( NULL );
879 #endif
880 
881  det = Q[0][0] * Q[1][1] - Q[0][1] * Q[1][0];
882 
883  if( det != 0.0 )
884  {
885  w.x = (-Q[0][2] * Q[1][1] + Q[1][2] * Q[0][1]) / det;
886  w.y = ( Q[0][2] * Q[1][0] - Q[1][2] * Q[0][0]) / det;
887  break;
888  }
889 
890  /* matrix is singular - lines are parallel. Add another,
891  * orthogonal axis, through the center of the unit square */
892  if( Q[0][0]>Q[1][1] )
893  {
894  v[0] = -Q[0][1];
895  v[1] = Q[0][0];
896  }
897  else if( Q[1][1] )
898  {
899  v[0] = -Q[1][1];
900  v[1] = Q[1][0];
901  }
902  else
903  {
904  v[0] = 1;
905  v[1] = 0;
906  }
907 
908  d = sq( v[0] ) + sq( v[1] );
909  v[2] = -v[1] * s.y - v[0] * s.x;
910 
911  for( l = 0; l<3; l++ )
912  {
913  for( k = 0; k<3; k++ )
914  {
915  Q[l][k] += v[l] * v[k] / d;
916  }
917  }
918  }
919 
920  dx = fabs( w.x - s.x );
921  dy = fabs( w.y - s.y );
922 
923  if( dx <= .5 && dy <= .5 )
924  {
925  pp->curve.vertex[i].x = w.x + x0;
926  pp->curve.vertex[i].y = w.y + y0;
927  continue;
928  }
929 
930  /* the minimum was not in the unit square; now minimize quadratic
931  * on boundary of square */
932  min = quadform( Q, s );
933  xmin = s.x;
934  ymin = s.y;
935 
936  if( Q[0][0] == 0.0 )
937  {
938  goto fixx;
939  }
940 
941  for( z = 0; z<2; z++ ) /* value of the y-coordinate */
942  {
943  w.y = s.y - 0.5 + z;
944  w.x = -(Q[0][1] * w.y + Q[0][2]) / Q[0][0];
945  dx = fabs( w.x - s.x );
946  cand = quadform( Q, w );
947 
948  if( dx <= .5 && cand < min )
949  {
950  min = cand;
951  xmin = w.x;
952  ymin = w.y;
953  }
954  }
955 
956 fixx:
957 
958  if( Q[1][1] == 0.0 )
959  {
960  goto corners;
961  }
962 
963  for( z = 0; z<2; z++ ) /* value of the x-coordinate */
964  {
965  w.x = s.x - 0.5 + z;
966  w.y = -(Q[1][0] * w.x + Q[1][2]) / Q[1][1];
967  dy = fabs( w.y - s.y );
968  cand = quadform( Q, w );
969 
970  if( dy <= .5 && cand < min )
971  {
972  min = cand;
973  xmin = w.x;
974  ymin = w.y;
975  }
976  }
977 
978 corners:
979 
980  /* check four corners */
981  for( l = 0; l<2; l++ )
982  {
983  for( k = 0; k<2; k++ )
984  {
985  w.x = s.x - 0.5 + l;
986  w.y = s.y - 0.5 + k;
987  cand = quadform( Q, w );
988 
989  if( cand < min )
990  {
991  min = cand;
992  xmin = w.x;
993  ymin = w.y;
994  }
995  }
996  }
997 
998  pp->curve.vertex[i].x = xmin + x0;
999  pp->curve.vertex[i].y = ymin + y0;
1000  continue;
1001  }
1002 
1003  free( ctr );
1004  free( dir );
1005  free( q );
1006  return 0;
1007 
1008 calloc_error:
1009  free( ctr );
1010  free( dir );
1011  free( q );
1012  return 1;
1013 }
long y
Definition: auxiliary.h:25
int privcurve_init(privcurve_t *curve, int n)
Definition: curve.cpp:98
static void pointslope(privpath_t *pp, int i, int j, dpoint_t *ctr, dpoint_t *dir)
Definition: trace.cpp:88
#define sq(a)
Definition: auxiliary.h:87
privcurve_t curve
Definition: curve.h:62
dpoint_t * vertex
Definition: curve.h:28
long x
Definition: auxiliary.h:24
point_t * pt
Definition: curve.h:53
#define SAFE_CALLOC(var, n, typ)
Definition: trace.cpp:28
static double quadform(quadform_t Q, dpoint_t w)
Definition: trace.cpp:179
#define mod(a, n)
Definition: greymap.cpp:24
double quadform_t[3][3]
Definition: trace.cpp:176
#define min(a, b)
Definition: auxiliary.h:85
static int bestpolygon ( privpath_t pp)
static

Definition at line 628 of file trace.cpp.

References potrace_privpath_s::len, potrace_privpath_s::lon, potrace_privpath_s::m, mod, penalty3(), potrace_privpath_s::po, and SAFE_CALLOC.

Referenced by process_path().

629 {
630  int i, j, m, k;
631  int n = pp->len;
632  double* pen = NULL; /* pen[n+1]: penalty vector */
633  int* prev = NULL; /* prev[n+1]: best path pointer vector */
634  int* clip0 = NULL; /* clip0[n]: longest segment pointer, non-cyclic */
635  int* clip1 = NULL; /* clip1[n+1]: backwards segment pointer, non-cyclic */
636  int* seg0 = NULL; /* seg0[m+1]: forward segment bounds, m<=n */
637  int* seg1 = NULL; /* seg1[m+1]: backward segment bounds, m<=n */
638  double thispen;
639  double best;
640  int c;
641 
642  SAFE_CALLOC( pen, n + 1, double );
643  SAFE_CALLOC( prev, n + 1, int );
644  SAFE_CALLOC( clip0, n, int );
645  SAFE_CALLOC( clip1, n + 1, int );
646  SAFE_CALLOC( seg0, n + 1, int );
647  SAFE_CALLOC( seg1, n + 1, int );
648 
649  /* calculate clipped paths */
650  for( i = 0; i<n; i++ )
651  {
652  c = mod( pp->lon[mod( i - 1, n )] - 1, n );
653 
654  if( c == i )
655  {
656  c = mod( i + 1, n );
657  }
658 
659  if( c < i )
660  {
661  clip0[i] = n;
662  }
663  else
664  {
665  clip0[i] = c;
666  }
667  }
668 
669  /* calculate backwards path clipping, non-cyclic. j <= clip0[i] iff
670  * clip1[j] <= i, for i,j=0..n. */
671  j = 1;
672 
673  for( i = 0; i<n; i++ )
674  {
675  while( j <= clip0[i] )
676  {
677  clip1[j] = i;
678  j++;
679  }
680  }
681 
682  /* calculate seg0[j] = longest path from 0 with j segments */
683  i = 0;
684 
685  for( j = 0; i<n; j++ )
686  {
687  seg0[j] = i;
688  i = clip0[i];
689  }
690 
691  seg0[j] = n;
692  m = j;
693 
694  /* calculate seg1[j] = longest path to n with m-j segments */
695  i = n;
696 
697  for( j = m; j>0; j-- )
698  {
699  seg1[j] = i;
700  i = clip1[i];
701  }
702 
703  seg1[0] = 0;
704 
705  /* now find the shortest path with m segments, based on penalty3 */
706  /* note: the outer 2 loops jointly have at most n iterations, thus
707  * the worst-case behavior here is quadratic. In practice, it is
708  * close to linear since the inner loop tends to be short. */
709  pen[0] = 0;
710 
711  for( j = 1; j<=m; j++ )
712  {
713  for( i = seg1[j]; i<=seg0[j]; i++ )
714  {
715  best = -1;
716 
717  for( k = seg0[j - 1]; k>=clip1[i]; k-- )
718  {
719  thispen = penalty3( pp, k, i ) + pen[k];
720 
721  if( best < 0 || thispen < best )
722  {
723  prev[i] = k;
724  best = thispen;
725  }
726  }
727 
728  pen[i] = best;
729  }
730  }
731 
732  pp->m = m;
733  SAFE_CALLOC( pp->po, m, int );
734 
735  /* read off shortest path */
736  for( i = n, j = m - 1; i>0; j-- )
737  {
738  i = prev[i];
739  pp->po[j] = i;
740  }
741 
742  free( pen );
743  free( prev );
744  free( clip0 );
745  free( clip1 );
746  free( seg0 );
747  free( seg1 );
748  return 0;
749 
750 calloc_error:
751  free( pen );
752  free( prev );
753  free( clip0 );
754  free( clip1 );
755  free( seg0 );
756  free( seg1 );
757  return 1;
758 }
static double penalty3(privpath_t *pp, int i, int j)
Definition: trace.cpp:570
#define SAFE_CALLOC(var, n, typ)
Definition: trace.cpp:28
#define mod(a, n)
Definition: greymap.cpp:24
static dpoint_t bezier ( double  t,
dpoint_t  p0,
dpoint_t  p1,
dpoint_t  p2,
dpoint_t  p3 
)
inlinestatic

Definition at line 259 of file trace.cpp.

References potrace_dpoint_s::x, and potrace_dpoint_s::y.

Referenced by SCH_LEGACY_PLUGIN_CACHE::loadBezier(), and opti_penalty().

260 {
261  double s = 1 - t;
262  dpoint_t res;
263 
264  /* Note: a good optimizing compiler (such as gcc-3) reduces the
265  * following to 16 multiplications, using common subexpression
266  * elimination. */
267 
268  res.x = s * s * s * p0.x + 3 * (s * s * t) * p1.x + 3 * (t * t * s) * p2.x + t * t * t * p3.x;
269  res.y = s * s * s * p0.y + 3 * (s * s * t) * p1.y + 3 * (t * t * s) * p2.y + t * t * t * p3.y;
270 
271  return res;
272 }
static int calc_lon ( privpath_t pp)
static

Definition at line 386 of file trace.cpp.

References abs, cyclic(), floordiv(), INFTY, potrace_privpath_s::len, potrace_privpath_s::lon, min, mod, potrace_privpath_s::pt, SAFE_CALLOC, sign(), point_s::x, xprod(), and point_s::y.

Referenced by process_path().

387 {
388  point_t* pt = pp->pt;
389  int n = pp->len;
390  int i, j, k, k1;
391  int ct[4], dir;
392  point_t constraint[2];
393  point_t cur;
394  point_t off;
395  int* pivk = NULL; /* pivk[n] */
396  int* nc = NULL; /* nc[n]: next corner */
397  point_t dk; /* direction of k-k1 */
398  int a, b, c, d;
399 
400  SAFE_CALLOC( pivk, n, int );
401  SAFE_CALLOC( nc, n, int );
402 
403  /* initialize the nc data structure. Point from each point to the
404  * furthest future point to which it is connected by a vertical or
405  * horizontal segment. We take advantage of the fact that there is
406  * always a direction change at 0 (due to the path decomposition
407  * algorithm). But even if this were not so, there is no harm, as
408  * in practice, correctness does not depend on the word "furthest"
409  * above. */
410  k = 0;
411 
412  for( i = n - 1; i>=0; i-- )
413  {
414  if( pt[i].x != pt[k].x && pt[i].y != pt[k].y )
415  {
416  k = i + 1; /* necessarily i<n-1 in this case */
417  }
418 
419  nc[i] = k;
420  }
421 
422  SAFE_CALLOC( pp->lon, n, int );
423 
424  /* determine pivot points: for each i, let pivk[i] be the furthest k
425  * such that all j with i<j<k lie on a line connecting i,k. */
426 
427  for( i = n - 1; i>=0; i-- )
428  {
429  ct[0] = ct[1] = ct[2] = ct[3] = 0;
430 
431  /* keep track of "directions" that have occurred */
432  dir = ( 3 + 3 * (pt[mod( i + 1, n )].x - pt[i].x) + (pt[mod( i + 1, n )].y - pt[i].y) ) / 2;
433  ct[dir]++;
434 
435  constraint[0].x = 0;
436  constraint[0].y = 0;
437  constraint[1].x = 0;
438  constraint[1].y = 0;
439 
440  /* find the next k such that no straight line from i to k */
441  k = nc[i];
442  k1 = i;
443 
444  while( 1 )
445  {
446  dir = ( 3 + 3 * sign( pt[k].x - pt[k1].x ) + sign( pt[k].y - pt[k1].y ) ) / 2;
447  ct[dir]++;
448 
449  /* if all four "directions" have occurred, cut this path */
450  if( ct[0] && ct[1] && ct[2] && ct[3] )
451  {
452  pivk[i] = k1;
453  goto foundk;
454  }
455 
456  cur.x = pt[k].x - pt[i].x;
457  cur.y = pt[k].y - pt[i].y;
458 
459  /* see if current constraint is violated */
460  if( xprod( constraint[0], cur ) < 0 || xprod( constraint[1], cur ) > 0 )
461  {
462  goto constraint_viol;
463  }
464 
465  /* else, update constraint */
466  if( abs( cur.x ) <= 1 && abs( cur.y ) <= 1 )
467  {
468  /* no constraint */
469  }
470  else
471  {
472  off.x = cur.x + ( ( cur.y>=0 && (cur.y>0 || cur.x<0) ) ? 1 : -1 );
473  off.y = cur.y + ( ( cur.x<=0 && (cur.x<0 || cur.y<0) ) ? 1 : -1 );
474 
475  if( xprod( constraint[0], off ) >= 0 )
476  {
477  constraint[0] = off;
478  }
479 
480  off.x = cur.x + ( ( cur.y<=0 && (cur.y<0 || cur.x<0) ) ? 1 : -1 );
481  off.y = cur.y + ( ( cur.x>=0 && (cur.x>0 || cur.y<0) ) ? 1 : -1 );
482 
483  if( xprod( constraint[1], off ) <= 0 )
484  {
485  constraint[1] = off;
486  }
487  }
488 
489  k1 = k;
490  k = nc[k1];
491 
492  if( !cyclic( k, i, k1 ) )
493  {
494  break;
495  }
496  }
497 
498 constraint_viol:
499  /* k1 was the last "corner" satisfying the current constraint, and
500  * k is the first one violating it. We now need to find the last
501  * point along k1..k which satisfied the constraint. */
502  dk.x = sign( pt[k].x - pt[k1].x );
503  dk.y = sign( pt[k].y - pt[k1].y );
504  cur.x = pt[k1].x - pt[i].x;
505  cur.y = pt[k1].y - pt[i].y;
506  /* find largest integer j such that xprod(constraint[0], cur+j*dk)
507  * >= 0 and xprod(constraint[1], cur+j*dk) <= 0. Use bilinearity
508  * of xprod. */
509  a = xprod( constraint[0], cur );
510  b = xprod( constraint[0], dk );
511  c = xprod( constraint[1], cur );
512  d = xprod( constraint[1], dk );
513  /* find largest integer j such that a+j*b>=0 and c+j*d<=0. This
514  * can be solved with integer arithmetic. */
515  j = INFTY;
516 
517  if( b<0 )
518  {
519  j = floordiv( a, -b );
520  }
521 
522  if( d>0 )
523  {
524  j = min( j, floordiv( -c, d ) );
525  }
526 
527  pivk[i] = mod( k1 + j, n );
528 foundk:
529  ;
530  } /* for i */
531 
532  /* clean up: for each i, let lon[i] be the largest k such that for
533  * all i' with i<=i'<k, i'<k<=pivk[i']. */
534 
535  j = pivk[n - 1];
536  pp->lon[n - 1] = j;
537 
538  for( i = n - 2; i>=0; i-- )
539  {
540  if( cyclic( i + 1, pivk[i], j ) )
541  {
542  j = pivk[i];
543  }
544 
545  pp->lon[i] = j;
546  }
547 
548  for( i = n - 1; cyclic( mod( i + 1, n ), j, pp->lon[i] ); i-- )
549  {
550  pp->lon[i] = j;
551  }
552 
553  free( pivk );
554  free( nc );
555  return 0;
556 
557 calloc_error:
558  free( pivk );
559  free( nc );
560  return 1;
561 }
static int cyclic(int a, int b, int c)
Definition: trace.cpp:73
long y
Definition: auxiliary.h:25
#define INFTY
Definition: trace.cpp:23
static int xprod(point_t p1, point_t p2)
Definition: trace.cpp:203
#define abs(a)
Definition: auxiliary.h:84
long x
Definition: auxiliary.h:24
point_t * pt
Definition: curve.h:53
#define SAFE_CALLOC(var, n, typ)
Definition: trace.cpp:28
#define mod(a, n)
Definition: greymap.cpp:24
#define min(a, b)
Definition: auxiliary.h:85
int sign(T val)
Definition: math_util.h:44
static int floordiv(int a, int n)
Definition: auxiliary.h:70
static int calc_sums ( privpath_t pp)
static

Definition at line 324 of file trace.cpp.

References potrace_privpath_s::len, potrace_privpath_s::pt, SAFE_CALLOC, potrace_privpath_s::sums, point_s::x, sums_s::x, potrace_privpath_s::x0, sums_s::x2, sums_s::xy, point_s::y, sums_s::y, potrace_privpath_s::y0, and sums_s::y2.

Referenced by process_path().

325 {
326  int i, x, y;
327  int n = pp->len;
328 
329  SAFE_CALLOC( pp->sums, pp->len + 1, sums_t );
330 
331  /* origin */
332  pp->x0 = pp->pt[0].x;
333  pp->y0 = pp->pt[0].y;
334 
335  /* preparatory computation for later fast summing */
336  pp->sums[0].x2 = pp->sums[0].xy = pp->sums[0].y2 = pp->sums[0].x = pp->sums[0].y = 0;
337 
338  for( i = 0; i<n; i++ )
339  {
340  x = pp->pt[i].x - pp->x0;
341  y = pp->pt[i].y - pp->y0;
342  pp->sums[i + 1].x = pp->sums[i].x + x;
343  pp->sums[i + 1].y = pp->sums[i].y + y;
344  pp->sums[i + 1].x2 = pp->sums[i].x2 + x * x;
345  pp->sums[i + 1].xy = pp->sums[i].xy + x * y;
346  pp->sums[i + 1].y2 = pp->sums[i].y2 + y * y;
347  }
348 
349  return 0;
350 
351 calloc_error:
352  return 1;
353 }
long y
Definition: auxiliary.h:25
long x
Definition: auxiliary.h:24
double y2
Definition: curve.h:41
double x
Definition: curve.h:37
double xy
Definition: curve.h:40
double y
Definition: curve.h:38
point_t * pt
Definition: curve.h:53
#define SAFE_CALLOC(var, n, typ)
Definition: trace.cpp:28
sums_t * sums
Definition: curve.h:57
Definition: curve.h:35
double x2
Definition: curve.h:39
static double cprod ( dpoint_t  p0,
dpoint_t  p1,
dpoint_t  p2,
dpoint_t  p3 
)
inlinestatic

Definition at line 210 of file trace.cpp.

References potrace_dpoint_s::x, and potrace_dpoint_s::y.

Referenced by opti_penalty(), and tangent().

211 {
212  double x1, y1, x2, y2;
213 
214  x1 = p1.x - p0.x;
215  y1 = p1.y - p0.y;
216  x2 = p3.x - p2.x;
217  y2 = p3.y - p2.y;
218 
219  return x1 * y2 - x2 * y1;
220 }
static int cyclic ( int  a,
int  b,
int  c 
)
inlinestatic

Definition at line 73 of file trace.cpp.

Referenced by calc_lon().

74 {
75  if( a<=c )
76  {
77  return a<=b && b<c;
78  }
79  else
80  {
81  return a<=b || b<c;
82  }
83 }
static double ddenom ( dpoint_t  p0,
dpoint_t  p2 
)
inlinestatic

Definition at line 64 of file trace.cpp.

References dorth_infty(), point_s::x, potrace_dpoint_s::x, point_s::y, and potrace_dpoint_s::y.

Referenced by smooth().

65 {
66  point_t r = dorth_infty( p0, p2 );
67 
68  return r.y * (p2.x - p0.x) - r.x * (p2.y - p0.y);
69 }
long y
Definition: auxiliary.h:25
long x
Definition: auxiliary.h:24
static point_t dorth_infty(dpoint_t p0, dpoint_t p2)
Definition: trace.cpp:37
static double ddist ( dpoint_t  p,
dpoint_t  q 
)
inlinestatic

Definition at line 252 of file trace.cpp.

References sq, potrace_dpoint_s::x, and potrace_dpoint_s::y.

Referenced by opti_penalty().

253 {
254  return sqrt( sq( p.x - q.x ) + sq( p.y - q.y ) );
255 }
#define sq(a)
Definition: auxiliary.h:87
static point_t dorth_infty ( dpoint_t  p0,
dpoint_t  p2 
)
inlinestatic

Definition at line 37 of file trace.cpp.

References sign(), point_s::x, potrace_dpoint_s::x, point_s::y, and potrace_dpoint_s::y.

Referenced by ddenom().

38 {
39  point_t r;
40 
41  r.y = sign( p2.x - p0.x );
42  r.x = -sign( p2.y - p0.y );
43 
44  return r;
45 }
long y
Definition: auxiliary.h:25
long x
Definition: auxiliary.h:24
int sign(T val)
Definition: math_util.h:44
static double dpara ( dpoint_t  p0,
dpoint_t  p1,
dpoint_t  p2 
)
inlinestatic

Definition at line 49 of file trace.cpp.

References potrace_dpoint_s::x, and potrace_dpoint_s::y.

Referenced by opti_penalty(), opticurve(), and smooth().

50 {
51  double x1, y1, x2, y2;
52 
53  x1 = p1.x - p0.x;
54  y1 = p1.y - p0.y;
55  x2 = p2.x - p0.x;
56  y2 = p2.y - p0.y;
57 
58  return x1 * y2 - x2 * y1;
59 }
static double iprod ( dpoint_t  p0,
dpoint_t  p1,
dpoint_t  p2 
)
inlinestatic

Definition at line 224 of file trace.cpp.

References potrace_dpoint_s::x, and potrace_dpoint_s::y.

Referenced by opti_penalty().

225 {
226  double x1, y1, x2, y2;
227 
228  x1 = p1.x - p0.x;
229  y1 = p1.y - p0.y;
230  x2 = p2.x - p0.x;
231  y2 = p2.y - p0.y;
232 
233  return x1 * x2 + y1 * y2;
234 }
static double iprod1 ( dpoint_t  p0,
dpoint_t  p1,
dpoint_t  p2,
dpoint_t  p3 
)
inlinestatic

Definition at line 238 of file trace.cpp.

References potrace_dpoint_s::x, and potrace_dpoint_s::y.

Referenced by opti_penalty().

239 {
240  double x1, y1, x2, y2;
241 
242  x1 = p1.x - p0.x;
243  y1 = p1.y - p0.y;
244  x2 = p3.x - p2.x;
245  y2 = p3.y - p2.y;
246 
247  return x1 * x2 + y1 * y2;
248 }
static int opti_penalty ( privpath_t pp,
int  i,
int  j,
opti_t res,
double  opttolerance,
int *  convc,
double *  areac 
)
static

Definition at line 1116 of file trace.cpp.

References privcurve_s::alpha, opti_s::alpha, bezier(), privcurve_s::c, opti_s::c, COS179, cprod(), potrace_privpath_s::curve, ddist(), dpara(), interval(), iprod(), iprod1(), mod, privcurve_s::n, opti_s::pen, opti_s::s, sign(), sq, opti_s::t, tangent(), and privcurve_s::vertex.

Referenced by opticurve().

1123 {
1124  int m = pp->curve.n;
1125  int k, k1, k2, conv, i1;
1126  double area, alpha, d, d1, d2;
1127  dpoint_t p0, p1, p2, p3, pt;
1128  double A, R, A1, A2, A3, A4;
1129  double s, t;
1130 
1131  /* check convexity, corner-freeness, and maximum bend < 179 degrees */
1132 
1133  if( i==j ) /* sanity - a full loop can never be an opticurve */
1134  {
1135  return 1;
1136  }
1137 
1138  k = i;
1139  i1 = mod( i + 1, m );
1140  k1 = mod( k + 1, m );
1141  conv = convc[k1];
1142 
1143  if( conv == 0 )
1144  {
1145  return 1;
1146  }
1147 
1148  d = ddist( pp->curve.vertex[i], pp->curve.vertex[i1] );
1149 
1150  for( k = k1; k!=j; k = k1 )
1151  {
1152  k1 = mod( k + 1, m );
1153  k2 = mod( k + 2, m );
1154 
1155  if( convc[k1] != conv )
1156  {
1157  return 1;
1158  }
1159 
1160  if( sign( cprod( pp->curve.vertex[i], pp->curve.vertex[i1], pp->curve.vertex[k1],
1161  pp->curve.vertex[k2] ) ) != conv )
1162  {
1163  return 1;
1164  }
1165 
1166  if( iprod1( pp->curve.vertex[i], pp->curve.vertex[i1], pp->curve.vertex[k1],
1167  pp->curve.vertex[k2] ) <
1168  d * ddist( pp->curve.vertex[k1], pp->curve.vertex[k2] ) * COS179 )
1169  {
1170  return 1;
1171  }
1172  }
1173 
1174  /* the curve we're working in: */
1175  p0 = pp->curve.c[mod( i, m )][2];
1176  p1 = pp->curve.vertex[mod( i + 1, m )];
1177  p2 = pp->curve.vertex[mod( j, m )];
1178  p3 = pp->curve.c[mod( j, m )][2];
1179 
1180  /* determine its area */
1181  area = areac[j] - areac[i];
1182  area -= dpara( pp->curve.vertex[0], pp->curve.c[i][2], pp->curve.c[j][2] ) / 2;
1183 
1184  if( i>=j )
1185  {
1186  area += areac[m];
1187  }
1188 
1189  /* find intersection o of p0p1 and p2p3. Let t,s such that o =
1190  * interval(t,p0,p1) = interval(s,p3,p2). Let A be the area of the
1191  * triangle (p0,o,p3). */
1192 
1193  A1 = dpara( p0, p1, p2 );
1194  A2 = dpara( p0, p1, p3 );
1195  A3 = dpara( p0, p2, p3 );
1196  /* A4 = dpara(p1, p2, p3); */
1197  A4 = A1 + A3 - A2;
1198 
1199  if( A2 == A1 ) /* this should never happen */
1200  {
1201  return 1;
1202  }
1203 
1204  t = A3 / (A3 - A4);
1205  s = A2 / (A2 - A1);
1206  A = A2 * t / 2.0;
1207 
1208  if( A == 0.0 ) /* this should never happen */
1209  {
1210  return 1;
1211  }
1212 
1213  R = area / A; /* relative area */
1214  alpha = 2 - sqrt( 4 - R / 0.3 ); /* overall alpha for p0-o-p3 curve */
1215 
1216  res->c[0] = interval( t * alpha, p0, p1 );
1217  res->c[1] = interval( s * alpha, p3, p2 );
1218  res->alpha = alpha;
1219  res->t = t;
1220  res->s = s;
1221 
1222  p1 = res->c[0];
1223  p2 = res->c[1]; /* the proposed curve is now (p0,p1,p2,p3) */
1224 
1225  res->pen = 0;
1226 
1227  /* calculate penalty */
1228  /* check tangency with edges */
1229  for( k = mod( i + 1, m ); k!=j; k = k1 )
1230  {
1231  k1 = mod( k + 1, m );
1232  t = tangent( p0, p1, p2, p3, pp->curve.vertex[k], pp->curve.vertex[k1] );
1233 
1234  if( t<-.5 )
1235  {
1236  return 1;
1237  }
1238 
1239  pt = bezier( t, p0, p1, p2, p3 );
1240  d = ddist( pp->curve.vertex[k], pp->curve.vertex[k1] );
1241 
1242  if( d == 0.0 ) /* this should never happen */
1243  {
1244  return 1;
1245  }
1246 
1247  d1 = dpara( pp->curve.vertex[k], pp->curve.vertex[k1], pt ) / d;
1248 
1249  if( fabs( d1 ) > opttolerance )
1250  {
1251  return 1;
1252  }
1253 
1254  if( iprod( pp->curve.vertex[k], pp->curve.vertex[k1],
1255  pt ) < 0 || iprod( pp->curve.vertex[k1], pp->curve.vertex[k], pt ) < 0 )
1256  {
1257  return 1;
1258  }
1259 
1260  res->pen += sq( d1 );
1261  }
1262 
1263  /* check corners */
1264  for( k = i; k!=j; k = k1 )
1265  {
1266  k1 = mod( k + 1, m );
1267  t = tangent( p0, p1, p2, p3, pp->curve.c[k][2], pp->curve.c[k1][2] );
1268 
1269  if( t<-.5 )
1270  {
1271  return 1;
1272  }
1273 
1274  pt = bezier( t, p0, p1, p2, p3 );
1275  d = ddist( pp->curve.c[k][2], pp->curve.c[k1][2] );
1276 
1277  if( d == 0.0 ) /* this should never happen */
1278  {
1279  return 1;
1280  }
1281 
1282  d1 = dpara( pp->curve.c[k][2], pp->curve.c[k1][2], pt ) / d;
1283  d2 = dpara( pp->curve.c[k][2], pp->curve.c[k1][2], pp->curve.vertex[k1] ) / d;
1284  d2 *= 0.75 * pp->curve.alpha[k1];
1285 
1286  if( d2 < 0 )
1287  {
1288  d1 = -d1;
1289  d2 = -d2;
1290  }
1291 
1292  if( d1 < d2 - opttolerance )
1293  {
1294  return 1;
1295  }
1296 
1297  if( d1 < d2 )
1298  {
1299  res->pen += sq( d1 - d2 );
1300  }
1301  }
1302 
1303  return 0;
1304 }
static double dpara(dpoint_t p0, dpoint_t p1, dpoint_t p2)
Definition: trace.cpp:49
static double iprod1(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3)
Definition: trace.cpp:238
double * alpha
Definition: curve.h:29
#define sq(a)
Definition: auxiliary.h:87
static double ddist(dpoint_t p, dpoint_t q)
Definition: trace.cpp:252
static dpoint_t bezier(double t, dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3)
Definition: trace.cpp:259
privcurve_t curve
Definition: curve.h:62
dpoint_t * vertex
Definition: curve.h:28
static double iprod(dpoint_t p0, dpoint_t p1, dpoint_t p2)
Definition: trace.cpp:224
static dpoint_t interval(double lambda, dpoint_t a, dpoint_t b)
Definition: auxiliary.h:43
double alpha
Definition: trace.cpp:1109
static double tangent(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3, dpoint_t q0, dpoint_t q1)
Definition: trace.cpp:278
dpoint_t c[2]
Definition: trace.cpp:1107
double s
Definition: trace.cpp:1108
static double cprod(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3)
Definition: trace.cpp:210
double pen
Definition: trace.cpp:1106
#define COS179
Definition: trace.cpp:25
double t
Definition: trace.cpp:1108
int n
Definition: curve.h:20
#define mod(a, n)
Definition: greymap.cpp:24
int sign(T val)
Definition: math_util.h:44
dpoint_t(* c)[3]
Definition: curve.h:22
static int opticurve ( privpath_t pp,
double  opttolerance 
)
static

Definition at line 1310 of file trace.cpp.

References privcurve_s::alpha, opti_s::alpha, privcurve_s::alpha0, privcurve_s::alphacurve, privcurve_s::beta, privcurve_s::c, opti_s::c, potrace_privpath_s::curve, dpara(), interval(), mod, privcurve_s::n, potrace_privpath_s::ocurve, opti_penalty(), opti_s::pen, POTRACE_CURVETO, privcurve_init(), opti_s::s, SAFE_CALLOC, sign(), opti_s::t, privcurve_s::tag, and privcurve_s::vertex.

Referenced by process_path().

1311 {
1312  int m = pp->curve.n;
1313  int* pt = NULL; /* pt[m+1] */
1314  double* pen = NULL; /* pen[m+1] */
1315  int* len = NULL; /* len[m+1] */
1316  opti_t* opt = NULL; /* opt[m+1] */
1317  int om;
1318  int i, j, r;
1319  opti_t o;
1320  dpoint_t p0;
1321  int i1;
1322  double area;
1323  double alpha;
1324  double* s = NULL;
1325  double* t = NULL;
1326 
1327  int* convc = NULL; /* conv[m]: pre-computed convexities */
1328  double* areac = NULL; /* cumarea[m+1]: cache for fast area computation */
1329 
1330  SAFE_CALLOC( pt, m + 1, int );
1331  SAFE_CALLOC( pen, m + 1, double );
1332  SAFE_CALLOC( len, m + 1, int );
1333  SAFE_CALLOC( opt, m + 1, opti_t );
1334  SAFE_CALLOC( convc, m, int );
1335  SAFE_CALLOC( areac, m + 1, double );
1336 
1337  /* pre-calculate convexity: +1 = right turn, -1 = left turn, 0 = corner */
1338  for( i = 0; i<m; i++ )
1339  {
1340  if( pp->curve.tag[i] == POTRACE_CURVETO )
1341  {
1342  convc[i] =
1343  sign( dpara( pp->curve.vertex[mod( i - 1, m )], pp->curve.vertex[i],
1344  pp->curve.vertex[mod( i + 1, m )] ) );
1345  }
1346  else
1347  {
1348  convc[i] = 0;
1349  }
1350  }
1351 
1352  /* pre-calculate areas */
1353  area = 0.0;
1354  areac[0] = 0.0;
1355  p0 = pp->curve.vertex[0];
1356 
1357  for( i = 0; i<m; i++ )
1358  {
1359  i1 = mod( i + 1, m );
1360 
1361  if( pp->curve.tag[i1] == POTRACE_CURVETO )
1362  {
1363  alpha = pp->curve.alpha[i1];
1364  area += 0.3 * alpha * (4 - alpha) * dpara( pp->curve.c[i][2],
1365  pp->curve.vertex[i1],
1366  pp->curve.c[i1][2] ) / 2;
1367  area += dpara( p0, pp->curve.c[i][2], pp->curve.c[i1][2] ) / 2;
1368  }
1369 
1370  areac[i + 1] = area;
1371  }
1372 
1373  pt[0] = -1;
1374  pen[0] = 0;
1375  len[0] = 0;
1376 
1377  /* Fixme: we always start from a fixed point -- should find the best
1378  * curve cyclically */
1379 
1380  for( j = 1; j<=m; j++ )
1381  {
1382  /* calculate best path from 0 to j */
1383  pt[j] = j - 1;
1384  pen[j] = pen[j - 1];
1385  len[j] = len[j - 1] + 1;
1386 
1387  for( i = j - 2; i>=0; i-- )
1388  {
1389  r = opti_penalty( pp, i, mod( j, m ), &o, opttolerance, convc, areac );
1390 
1391  if( r )
1392  {
1393  break;
1394  }
1395 
1396  if( len[j] > len[i] + 1 || (len[j] == len[i] + 1 && pen[j] > pen[i] + o.pen) )
1397  {
1398  pt[j] = i;
1399  pen[j] = pen[i] + o.pen;
1400  len[j] = len[i] + 1;
1401  opt[j] = o;
1402  }
1403  }
1404  }
1405 
1406  om = len[m];
1407  r = privcurve_init( &pp->ocurve, om );
1408 
1409  if( r )
1410  {
1411  goto calloc_error;
1412  }
1413 
1414  SAFE_CALLOC( s, om, double );
1415  SAFE_CALLOC( t, om, double );
1416 
1417  j = m;
1418 
1419  for( i = om - 1; i>=0; i-- )
1420  {
1421  if( pt[j]==j - 1 )
1422  {
1423  pp->ocurve.tag[i] = pp->curve.tag[mod( j, m )];
1424  pp->ocurve.c[i][0] = pp->curve.c[mod( j, m )][0];
1425  pp->ocurve.c[i][1] = pp->curve.c[mod( j, m )][1];
1426  pp->ocurve.c[i][2] = pp->curve.c[mod( j, m )][2];
1427  pp->ocurve.vertex[i] = pp->curve.vertex[mod( j, m )];
1428  pp->ocurve.alpha[i] = pp->curve.alpha[mod( j, m )];
1429  pp->ocurve.alpha0[i] = pp->curve.alpha0[mod( j, m )];
1430  pp->ocurve.beta[i] = pp->curve.beta[mod( j, m )];
1431  s[i] = t[i] = 1.0;
1432  }
1433  else
1434  {
1435  pp->ocurve.tag[i] = POTRACE_CURVETO;
1436  pp->ocurve.c[i][0] = opt[j].c[0];
1437  pp->ocurve.c[i][1] = opt[j].c[1];
1438  pp->ocurve.c[i][2] = pp->curve.c[mod( j, m )][2];
1439  pp->ocurve.vertex[i] = interval( opt[j].s, pp->curve.c[mod( j,
1440  m )][2],
1441  pp->curve.vertex[mod( j, m )] );
1442  pp->ocurve.alpha[i] = opt[j].alpha;
1443  pp->ocurve.alpha0[i] = opt[j].alpha;
1444  s[i] = opt[j].s;
1445  t[i] = opt[j].t;
1446  }
1447 
1448  j = pt[j];
1449  }
1450 
1451  /* calculate beta parameters */
1452  for( i = 0; i<om; i++ )
1453  {
1454  i1 = mod( i + 1, om );
1455  pp->ocurve.beta[i] = s[i] / (s[i] + t[i1]);
1456  }
1457 
1458  pp->ocurve.alphacurve = 1;
1459 
1460  free( pt );
1461  free( pen );
1462  free( len );
1463  free( opt );
1464  free( s );
1465  free( t );
1466  free( convc );
1467  free( areac );
1468  return 0;
1469 
1470 calloc_error:
1471  free( pt );
1472  free( pen );
1473  free( len );
1474  free( opt );
1475  free( s );
1476  free( t );
1477  free( convc );
1478  free( areac );
1479  return 1;
1480 }
static double dpara(dpoint_t p0, dpoint_t p1, dpoint_t p2)
Definition: trace.cpp:49
int privcurve_init(privcurve_t *curve, int n)
Definition: curve.cpp:98
double * alpha
Definition: curve.h:29
privcurve_t curve
Definition: curve.h:62
int alphacurve
Definition: curve.h:27
dpoint_t * vertex
Definition: curve.h:28
static dpoint_t interval(double lambda, dpoint_t a, dpoint_t b)
Definition: auxiliary.h:43
privcurve_t ocurve
Definition: curve.h:63
double alpha
Definition: trace.cpp:1109
dpoint_t c[2]
Definition: trace.cpp:1107
double s
Definition: trace.cpp:1108
double * beta
Definition: curve.h:31
#define POTRACE_CURVETO
Definition: potracelib.h:78
static int opti_penalty(privpath_t *pp, int i, int j, opti_t *res, double opttolerance, int *convc, double *areac)
Definition: trace.cpp:1116
#define SAFE_CALLOC(var, n, typ)
Definition: trace.cpp:28
double pen
Definition: trace.cpp:1106
double * alpha0
Definition: curve.h:30
int * tag
Definition: curve.h:21
double t
Definition: trace.cpp:1108
int n
Definition: curve.h:20
#define mod(a, n)
Definition: greymap.cpp:24
int sign(T val)
Definition: math_util.h:44
dpoint_t(* c)[3]
Definition: curve.h:22
static double penalty3 ( privpath_t pp,
int  i,
int  j 
)
static

Definition at line 570 of file trace.cpp.

References potrace_privpath_s::len, potrace_privpath_s::pt, potrace_privpath_s::sums, point_s::x, sums_s::x, sums_s::x2, sums_s::xy, point_s::y, sums_s::y, and sums_s::y2.

Referenced by bestpolygon().

571 {
572  int n = pp->len;
573  point_t* pt = pp->pt;
574  sums_t* sums = pp->sums;
575 
576  /* assume 0<=i<j<=n */
577  double x, y, x2, xy, y2;
578  double k;
579  double a, b, c, s;
580  double px, py, ex, ey;
581 
582  int r = 0; /* rotations from i to j */
583 
584  if( j>=n )
585  {
586  j -= n;
587  r = 1;
588  }
589 
590  /* critical inner loop: the "if" gives a 4.6 percent speedup */
591  if( r == 0 )
592  {
593  x = sums[j + 1].x - sums[i].x;
594  y = sums[j + 1].y - sums[i].y;
595  x2 = sums[j + 1].x2 - sums[i].x2;
596  xy = sums[j + 1].xy - sums[i].xy;
597  y2 = sums[j + 1].y2 - sums[i].y2;
598  k = j + 1 - i;
599  }
600  else
601  {
602  x = sums[j + 1].x - sums[i].x + sums[n].x;
603  y = sums[j + 1].y - sums[i].y + sums[n].y;
604  x2 = sums[j + 1].x2 - sums[i].x2 + sums[n].x2;
605  xy = sums[j + 1].xy - sums[i].xy + sums[n].xy;
606  y2 = sums[j + 1].y2 - sums[i].y2 + sums[n].y2;
607  k = j + 1 - i + n;
608  }
609 
610  px = (pt[i].x + pt[j].x) / 2.0 - pt[0].x;
611  py = (pt[i].y + pt[j].y) / 2.0 - pt[0].y;
612  ey = (pt[j].x - pt[i].x);
613  ex = -(pt[j].y - pt[i].y);
614 
615  a = ( (x2 - 2 * x * px) / k + px * px );
616  b = ( (xy - x * py - y * px) / k + px * py );
617  c = ( (y2 - 2 * y * py) / k + py * py );
618 
619  s = ex * ex * a + 2 * ex * ey * b + ey * ey * c;
620 
621  return sqrt( s );
622 }
long y
Definition: auxiliary.h:25
long x
Definition: auxiliary.h:24
double y2
Definition: curve.h:41
double x
Definition: curve.h:37
double xy
Definition: curve.h:40
double y
Definition: curve.h:38
point_t * pt
Definition: curve.h:53
sums_t * sums
Definition: curve.h:57
Definition: curve.h:35
double x2
Definition: curve.h:39
static void pointslope ( privpath_t pp,
int  i,
int  j,
dpoint_t ctr,
dpoint_t dir 
)
static

Definition at line 88 of file trace.cpp.

References potrace_privpath_s::len, potrace_privpath_s::sums, sums_s::x, potrace_dpoint_s::x, sums_s::x2, sums_s::xy, sums_s::y, potrace_dpoint_s::y, and sums_s::y2.

Referenced by adjust_vertices().

89 {
90  /* assume i<j */
91 
92  int n = pp->len;
93  sums_t* sums = pp->sums;
94 
95  double x, y, x2, xy, y2;
96  double k;
97  double a, b, c, lambda2, l;
98  int r = 0; /* rotations from i to j */
99 
100  while( j>=n )
101  {
102  j -= n;
103  r += 1;
104  }
105 
106  while( i>=n )
107  {
108  i -= n;
109  r -= 1;
110  }
111 
112  while( j<0 )
113  {
114  j += n;
115  r -= 1;
116  }
117 
118  while( i<0 )
119  {
120  i += n;
121  r += 1;
122  }
123 
124  x = sums[j + 1].x - sums[i].x + r * sums[n].x;
125  y = sums[j + 1].y - sums[i].y + r * sums[n].y;
126  x2 = sums[j + 1].x2 - sums[i].x2 + r * sums[n].x2;
127  xy = sums[j + 1].xy - sums[i].xy + r * sums[n].xy;
128  y2 = sums[j + 1].y2 - sums[i].y2 + r * sums[n].y2;
129  k = j + 1 - i + r * n;
130 
131  ctr->x = x / k;
132  ctr->y = y / k;
133 
134  a = (x2 - (double) x * x / k) / k;
135  b = (xy - (double) x * y / k) / k;
136  c = (y2 - (double) y * y / k) / k;
137 
138  lambda2 = ( a + c + sqrt( (a - c) * (a - c) + 4 * b * b ) ) / 2; /* larger e.value */
139 
140  /* now find e.vector for lambda2 */
141  a -= lambda2;
142  c -= lambda2;
143 
144  if( fabs( a ) >= fabs( c ) )
145  {
146  l = sqrt( a * a + b * b );
147 
148  if( l!=0 )
149  {
150  dir->x = -b / l;
151  dir->y = a / l;
152  }
153  }
154  else
155  {
156  l = sqrt( c * c + b * b );
157 
158  if( l!=0 )
159  {
160  dir->x = -c / l;
161  dir->y = b / l;
162  }
163  }
164 
165  if( l==0 )
166  {
167  dir->x = dir->y = 0; /* sometimes this can happen when k=4:
168  * the two eigenvalues coincide */
169  }
170 }
double y2
Definition: curve.h:41
double x
Definition: curve.h:37
double xy
Definition: curve.h:40
double y
Definition: curve.h:38
sums_t * sums
Definition: curve.h:57
Definition: curve.h:35
double x2
Definition: curve.h:39
int process_path ( path_t plist,
const potrace_param_t param,
progress_t progress 
)

Definition at line 1489 of file trace.cpp.

References adjust_vertices(), potrace_param_s::alphamax, bestpolygon(), calc_lon(), calc_sums(), progress_s::callback, potrace_privpath_s::curve, potrace_path_s::curve, potrace_privpath_s::fcurve, potrace_privpath_s::len, list_forall, potrace_privpath_s::ocurve, potrace_param_s::opticurve, opticurve(), potrace_param_s::opttolerance, potrace_path_s::priv, privcurve_to_curve(), progress_update(), reverse(), potrace_path_s::sign, smooth(), and TRY.

Referenced by potrace_trace().

1490 {
1491  path_t* p;
1492  double nn = 0, cn = 0;
1493 
1494  if( progress->callback )
1495  {
1496  /* precompute task size for progress estimates */
1497  nn = 0;
1498  list_forall( p, plist ) {
1499  nn += p->priv->len;
1500  }
1501  cn = 0;
1502  }
1503 
1504  /* call downstream function with each path */
1505  list_forall( p, plist ) {
1506  TRY( calc_sums( p->priv ) );
1507  TRY( calc_lon( p->priv ) );
1508  TRY( bestpolygon( p->priv ) );
1509  TRY( adjust_vertices( p->priv ) );
1510 
1511  if( p->sign == '-' ) /* reverse orientation of negative paths */
1512  {
1513  reverse( &p->priv->curve );
1514  }
1515 
1516  smooth( &p->priv->curve, param->alphamax );
1517 
1518  if( param->opticurve )
1519  {
1520  TRY( opticurve( p->priv, param->opttolerance ) );
1521  p->priv->fcurve = &p->priv->ocurve;
1522  }
1523  else
1524  {
1525  p->priv->fcurve = &p->priv->curve;
1526  }
1527 
1528  privcurve_to_curve( p->priv->fcurve, &p->curve );
1529 
1530  if( progress->callback )
1531  {
1532  cn += p->priv->len;
1533  progress_update( cn / nn, progress );
1534  }
1535  }
1536 
1537  progress_update( 1.0, progress );
1538 
1539  return 0;
1540 
1541 try_error:
1542  return 1;
1543 }
static int calc_lon(privpath_t *pp)
Definition: trace.cpp:386
#define TRY(x)
Definition: trace.cpp:1485
static int bestpolygon(privpath_t *pp)
Definition: trace.cpp:628
double alphamax
Definition: potracelib.h:42
privcurve_t curve
Definition: curve.h:62
static int calc_sums(privpath_t *pp)
Definition: trace.cpp:324
privcurve_t ocurve
Definition: curve.h:63
static void smooth(privcurve_t *curve, double alphamax)
Definition: trace.cpp:1036
struct potrace_privpath_s * priv
Definition: potracelib.h:103
privcurve_t * fcurve
Definition: curve.h:64
static int opticurve(privpath_t *pp, double opttolerance)
Definition: trace.cpp:1310
void privcurve_to_curve(privcurve_t *pc, potrace_curve_t *c)
Definition: curve.cpp:122
#define list_forall(elt, list)
Definition: lists.h:42
double opttolerance
Definition: potracelib.h:44
static void progress_update(double d, progress_t *prog)
Definition: progress.h:29
static void reverse(privcurve_t *curve)
Definition: trace.cpp:1020
void(* callback)(double progress, void *privdata)
Definition: progress.h:17
static int adjust_vertices(privpath_t *pp)
Definition: trace.cpp:769
potrace_curve_t curve
Definition: potracelib.h:96
static double quadform ( quadform_t  Q,
dpoint_t  w 
)
inlinestatic

Definition at line 179 of file trace.cpp.

References potrace_dpoint_s::x, and potrace_dpoint_s::y.

Referenced by adjust_vertices().

180 {
181  double v[3];
182  int i, j;
183  double sum;
184 
185  v[0] = w.x;
186  v[1] = w.y;
187  v[2] = 1;
188  sum = 0.0;
189 
190  for( i = 0; i<3; i++ )
191  {
192  for( j = 0; j<3; j++ )
193  {
194  sum += v[i] * Q[i][j] * v[j];
195  }
196  }
197 
198  return sum;
199 }
static void reverse ( privcurve_t curve)
static

Definition at line 1020 of file trace.cpp.

References privcurve_s::n, and privcurve_s::vertex.

Referenced by AUTOPLACER::choose_side_for_fields(), process_path(), PNS::LINE::Reverse(), SHAPE_LINE_CHAIN::Reverse(), and PNS::SHOVE::runOptimizer().

1021 {
1022  int m = curve->n;
1023  int i, j;
1024  dpoint_t tmp;
1025 
1026  for( i = 0, j = m - 1; i<j; i++, j-- )
1027  {
1028  tmp = curve->vertex[i];
1029  curve->vertex[i] = curve->vertex[j];
1030  curve->vertex[j] = tmp;
1031  }
1032 }
dpoint_t * vertex
Definition: curve.h:28
int n
Definition: curve.h:20
static void smooth ( privcurve_t curve,
double  alphamax 
)
static

Definition at line 1036 of file trace.cpp.

References privcurve_s::alpha, privcurve_s::alpha0, privcurve_s::alphacurve, privcurve_s::beta, privcurve_s::c, ddenom(), dpara(), interval(), mod, privcurve_s::n, POTRACE_CORNER, POTRACE_CURVETO, privcurve_s::tag, and privcurve_s::vertex.

Referenced by process_path().

1037 {
1038  int m = curve->n;
1039 
1040  int i, j, k;
1041  double dd, denom, alpha;
1042  dpoint_t p2, p3, p4;
1043 
1044  /* examine each vertex and find its best fit */
1045  for( i = 0; i<m; i++ )
1046  {
1047  j = mod( i + 1, m );
1048  k = mod( i + 2, m );
1049  p4 = interval( 1 / 2.0, curve->vertex[k], curve->vertex[j] );
1050 
1051  denom = ddenom( curve->vertex[i], curve->vertex[k] );
1052 
1053  if( denom != 0.0 )
1054  {
1055  dd = dpara( curve->vertex[i], curve->vertex[j], curve->vertex[k] ) / denom;
1056  dd = fabs( dd );
1057  alpha = dd>1 ? (1 - 1.0 / dd) : 0;
1058  alpha = alpha / 0.75;
1059  }
1060  else
1061  {
1062  alpha = 4 / 3.0;
1063  }
1064 
1065  curve->alpha0[j] = alpha; /* remember "original" value of alpha */
1066 
1067  if( alpha >= alphamax ) /* pointed corner */
1068  {
1069  curve->tag[j] = POTRACE_CORNER;
1070  curve->c[j][1] = curve->vertex[j];
1071  curve->c[j][2] = p4;
1072  }
1073  else
1074  {
1075  if( alpha < 0.55 )
1076  {
1077  alpha = 0.55;
1078  }
1079  else if( alpha > 1 )
1080  {
1081  alpha = 1;
1082  }
1083 
1084  p2 = interval( .5 + .5 * alpha, curve->vertex[i], curve->vertex[j] );
1085  p3 = interval( .5 + .5 * alpha, curve->vertex[k], curve->vertex[j] );
1086  curve->tag[j] = POTRACE_CURVETO;
1087  curve->c[j][0] = p2;
1088  curve->c[j][1] = p3;
1089  curve->c[j][2] = p4;
1090  }
1091 
1092  curve->alpha[j] = alpha; /* store the "cropped" value of alpha */
1093  curve->beta[j] = 0.5;
1094  }
1095 
1096  curve->alphacurve = 1;
1097 }
static double dpara(dpoint_t p0, dpoint_t p1, dpoint_t p2)
Definition: trace.cpp:49
double * alpha
Definition: curve.h:29
int alphacurve
Definition: curve.h:27
#define POTRACE_CORNER
Definition: potracelib.h:79
dpoint_t * vertex
Definition: curve.h:28
static dpoint_t interval(double lambda, dpoint_t a, dpoint_t b)
Definition: auxiliary.h:43
double * beta
Definition: curve.h:31
#define POTRACE_CURVETO
Definition: potracelib.h:78
static double ddenom(dpoint_t p0, dpoint_t p2)
Definition: trace.cpp:64
double * alpha0
Definition: curve.h:30
int * tag
Definition: curve.h:21
int n
Definition: curve.h:20
#define mod(a, n)
Definition: greymap.cpp:24
dpoint_t(* c)[3]
Definition: curve.h:22
static double tangent ( dpoint_t  p0,
dpoint_t  p1,
dpoint_t  p2,
dpoint_t  p3,
dpoint_t  q0,
dpoint_t  q1 
)
static

Definition at line 278 of file trace.cpp.

References cprod().

Referenced by opti_penalty().

280 {
281  double A, B, C; /* (1-t)^2 A + 2(1-t)t B + t^2 C = 0 */
282  double a, b, c; /* a t^2 + b t + c = 0 */
283  double d, s, r1, r2;
284 
285  A = cprod( p0, p1, q0, q1 );
286  B = cprod( p1, p2, q0, q1 );
287  C = cprod( p2, p3, q0, q1 );
288 
289  a = A - 2 * B + C;
290  b = -2 * A + 2 * B;
291  c = A;
292 
293  d = b * b - 4 * a * c;
294 
295  if( a==0 || d<0 )
296  {
297  return -1.0;
298  }
299 
300  s = sqrt( d );
301 
302  r1 = (-b + s) / (2 * a);
303  r2 = (-b - s) / (2 * a);
304 
305  if( r1 >= 0 && r1 <= 1 )
306  {
307  return r1;
308  }
309  else if( r2 >= 0 && r2 <= 1 )
310  {
311  return r2;
312  }
313  else
314  {
315  return -1.0;
316  }
317 }
static double cprod(dpoint_t p0, dpoint_t p1, dpoint_t p2, dpoint_t p3)
Definition: trace.cpp:210
static int xprod ( point_t  p1,
point_t  p2 
)
inlinestatic

Definition at line 203 of file trace.cpp.

References point_s::x, and point_s::y.

Referenced by calc_lon().

204 {
205  return p1.x * p2.y - p1.y * p2.x;
206 }
long y
Definition: auxiliary.h:25
long x
Definition: auxiliary.h:24