Commit ec807c7
Add cubic Beziér curve fitting code from Graphics Gems
Changed files (7)
addedlibbirdfont/GraphicsGems/GGVecLib.caddedlibbirdfont/GraphicsGems/GraphicsGems.haddedlibbirdfont/GraphicsGems/fit_cubic.cmodifiedlibbirdfont/Path.valamodifiedlibbirdfont/PenTool.valamodifiedlibbirdfont/StrokeTool.valamodifiedscripts/bavala.py
Diff
new file mode 100644@@ -0,0 +1,458 @@+/* +2d and 3d Vector C Library +by Andrew Glassner+from "Graphics Gems", Academic Press, 1990+*/++/*+ * EULA: The Graphics Gems code is copyright-protected. In other words, + * you cannot claim the text of the code as your own and resell it. + * Using the code is permitted in any program, product, or library, + * non-commercial or commercial. Giving credit is not required, though + * is a nice gesture. The code comes as-is, and if there are any flaws + * or problems with any Gems code, nobody involved with Gems - authors,+ * editors, publishers, or webmasters - are to be held responsible. + * Basically, don't be a jerk, and remember that anything free comes + * with no guarantee. + */++#include <math.h>+#include "GraphicsGems.h"++/******************/+/* 2d Library */+/******************/++/* returns squared length of input vector */ +double V2SquaredLength(a) +Vector2 *a;+{ return((a->x * a->x)+(a->y * a->y));+ }+ +/* returns length of input vector */+double V2Length(a) +Vector2 *a;+{+ return(sqrt(V2SquaredLength(a)));+ }+ +/* negates the input vector and returns it */+Vector2 *V2Negate(v) +Vector2 *v;+{+ v->x = -v->x; v->y = -v->y;+ return(v);+ }++/* normalizes the input vector and returns it */+Vector2 *V2Normalize(v) +Vector2 *v;+{+double len = V2Length(v);+ if (len != 0.0) { v->x /= len; v->y /= len; }+ return(v);+ }+++/* scales the input vector to the new length and returns it */+Vector2 *V2Scale(v, newlen) +Vector2 *v;+double newlen;+{+double len = V2Length(v);+ if (len != 0.0) { v->x *= newlen/len; v->y *= newlen/len; }+ return(v);+ }++/* return vector sum c = a+b */+Vector2 *V2Add(a, b, c)+Vector2 *a, *b, *c;+{+ c->x = a->x+b->x; c->y = a->y+b->y;+ return(c);+ }+ +/* return vector difference c = a-b */+Vector2 *V2Sub(a, b, c)+Vector2 *a, *b, *c;+{+ c->x = a->x-b->x; c->y = a->y-b->y;+ return(c);+ }++/* return the dot product of vectors a and b */+double V2Dot(a, b) +Vector2 *a, *b; +{+ return((a->x*b->x)+(a->y*b->y));+ }++/* linearly interpolate between vectors by an amount alpha */+/* and return the resulting vector. */+/* When alpha=0, result=lo. When alpha=1, result=hi. */+Vector2 *V2Lerp(lo, hi, alpha, result) +Vector2 *lo, *hi, *result; +double alpha;+{+ result->x = LERP(alpha, lo->x, hi->x);+ result->y = LERP(alpha, lo->y, hi->y);+ return(result);+ }+++/* make a linear combination of two vectors and return the result. */+/* result = (a * ascl) + (b * bscl) */+Vector2 *V2Combine (a, b, result, ascl, bscl) +Vector2 *a, *b, *result;+double ascl, bscl;+{+ result->x = (ascl * a->x) + (bscl * b->x);+ result->y = (ascl * a->y) + (bscl * b->y);+ return(result);+ }++/* multiply two vectors together component-wise */+Vector2 *V2Mul (a, b, result) +Vector2 *a, *b, *result;+{+ result->x = a->x * b->x;+ result->y = a->y * b->y;+ return(result);+ }++/* return the distance between two points */+double V2DistanceBetween2Points(a, b)+Point2 *a, *b;+{+double dx = a->x - b->x;+double dy = a->y - b->y;+ return(sqrt((dx*dx)+(dy*dy)));+ }++/* return the vector perpendicular to the input vector a */+Vector2 *V2MakePerpendicular(a, ap)+Vector2 *a, *ap;+{+ ap->x = -a->y;+ ap->y = a->x;+ return(ap);+ }++/* create, initialize, and return a new vector */+Vector2 *V2New(x, y)+double x, y;+{+Vector2 *v = NEWTYPE(Vector2);+ v->x = x; v->y = y; + return(v);+ }+ ++/* create, initialize, and return a duplicate vector */+Vector2 *V2Duplicate(a)+Vector2 *a;+{+Vector2 *v = NEWTYPE(Vector2);+ v->x = a->x; v->y = a->y; + return(v);+ }+ +/* multiply a point by a matrix and return the transformed point */+Point2 *V2MulPointByMatrix(p, m)+Point2 *p;+Matrix3 *m;+{+double w;+Point2 ptmp;+ ptmp.x = (p->x * m->element[0][0]) + + (p->y * m->element[1][0]) + m->element[2][0];+ ptmp.y = (p->x * m->element[0][1]) + + (p->y * m->element[1][1]) + m->element[2][1];+ w = (p->x * m->element[0][2]) + + (p->y * m->element[1][2]) + m->element[2][2];+ if (w != 0.0) { ptmp.x /= w; ptmp.y /= w; }+ *p = ptmp;+ return(p);+ }++/* multiply together matrices c = ab */+/* note that c must not point to either of the input matrices */+Matrix3 *V2MatMul(a, b, c)+Matrix3 *a, *b, *c;+{+int i, j, k;+ for (i=0; i<3; i++) {+ for (j=0; j<3; j++) {+ c->element[i][j] = 0;+ for (k=0; k<3; k++) c->element[i][j] += + a->element[i][k] * b->element[k][j];+ }+ }+ return(c);+ }+++++/******************/+/* 3d Library */+/******************/+ +/* returns squared length of input vector */ +double V3SquaredLength(a) +Vector3 *a;+{+ return((a->x * a->x)+(a->y * a->y)+(a->z * a->z));+ }++/* returns length of input vector */+double V3Length(a) +Vector3 *a;+{+ return(sqrt(V3SquaredLength(a)));+ }++/* negates the input vector and returns it */+Vector3 *V3Negate(v) +Vector3 *v;+{+ v->x = -v->x; v->y = -v->y; v->z = -v->z;+ return(v);+ }++/* normalizes the input vector and returns it */+Vector3 *V3Normalize(v) +Vector3 *v;+{+double len = V3Length(v);+ if (len != 0.0) { v->x /= len; v->y /= len; v->z /= len; }+ return(v);+ }++/* scales the input vector to the new length and returns it */+Vector3 *V3Scale(v, newlen) +Vector3 *v;+double newlen;+{+double len = V3Length(v);+ if (len != 0.0) {+ v->x *= newlen/len; v->y *= newlen/len; v->z *= newlen/len;+ }+ return(v);+ }+++/* return vector sum c = a+b */+Vector3 *V3Add(a, b, c)+Vector3 *a, *b, *c;+{+ c->x = a->x+b->x; c->y = a->y+b->y; c->z = a->z+b->z;+ return(c);+ }+ +/* return vector difference c = a-b */+Vector3 *V3Sub(a, b, c)+Vector3 *a, *b, *c;+{+ c->x = a->x-b->x; c->y = a->y-b->y; c->z = a->z-b->z;+ return(c);+ }++/* return the dot product of vectors a and b */+double V3Dot(a, b) +Vector3 *a, *b; +{+ return((a->x*b->x)+(a->y*b->y)+(a->z*b->z));+ }++/* linearly interpolate between vectors by an amount alpha */+/* and return the resulting vector. */+/* When alpha=0, result=lo. When alpha=1, result=hi. */+Vector3 *V3Lerp(lo, hi, alpha, result) +Vector3 *lo, *hi, *result; +double alpha;+{+ result->x = LERP(alpha, lo->x, hi->x);+ result->y = LERP(alpha, lo->y, hi->y);+ result->z = LERP(alpha, lo->z, hi->z);+ return(result);+ }++/* make a linear combination of two vectors and return the result. */+/* result = (a * ascl) + (b * bscl) */+Vector3 *V3Combine (a, b, result, ascl, bscl) +Vector3 *a, *b, *result;+double ascl, bscl;+{+ result->x = (ascl * a->x) + (bscl * b->x);+ result->y = (ascl * a->y) + (bscl * b->y);+ result->z = (ascl * a->z) + (bscl * b->z);+ return(result);+ }+++/* multiply two vectors together component-wise and return the result */+Vector3 *V3Mul (a, b, result) +Vector3 *a, *b, *result;+{+ result->x = a->x * b->x;+ result->y = a->y * b->y;+ result->z = a->z * b->z;+ return(result);+ }++/* return the distance between two points */+double V3DistanceBetween2Points(a, b)+Point3 *a, *b;+{+double dx = a->x - b->x;+double dy = a->y - b->y;+double dz = a->z - b->z;+ return(sqrt((dx*dx)+(dy*dy)+(dz*dz)));+ }++/* return the cross product c = a cross b */+Vector3 *V3Cross(a, b, c)+Vector3 *a, *b, *c;+{+ c->x = (a->y*b->z) - (a->z*b->y);+ c->y = (a->z*b->x) - (a->x*b->z);+ c->z = (a->x*b->y) - (a->y*b->x);+ return(c);+ }++/* create, initialize, and return a new vector */+Vector3 *V3New(x, y, z)+double x, y, z;+{+Vector3 *v = NEWTYPE(Vector3);+ v->x = x; v->y = y; v->z = z;+ return(v);+ }++/* create, initialize, and return a duplicate vector */+Vector3 *V3Duplicate(a)+Vector3 *a;+{+Vector3 *v = NEWTYPE(Vector3);+ v->x = a->x; v->y = a->y; v->z = a->z;+ return(v);+ }++ +/* multiply a point by a matrix and return the transformed point */+Point3 *V3MulPointByMatrix(p, m)+Point3 *p;+Matrix4 *m;+{+double w;+Point3 ptmp;+ ptmp.x = (p->x * m->element[0][0]) + (p->y * m->element[1][0]) + + (p->z * m->element[2][0]) + m->element[3][0];+ ptmp.y = (p->x * m->element[0][1]) + (p->y * m->element[1][1]) + + (p->z * m->element[2][1]) + m->element[3][1];+ ptmp.z = (p->x * m->element[0][2]) + (p->y * m->element[1][2]) + + (p->z * m->element[2][2]) + m->element[3][2];+ w = (p->x * m->element[0][3]) + (p->y * m->element[1][3]) + + (p->z * m->element[2][3]) + m->element[3][3];+ if (w != 0.0) { ptmp.x /= w; ptmp.y /= w; ptmp.z /= w; }+ *p = ptmp;+ return(p);+ }++/* multiply together matrices c = ab */+/* note that c must not point to either of the input matrices */+Matrix4 *V3MatMul(a, b, c)+Matrix4 *a, *b, *c;+{+int i, j, k;+ for (i=0; i<4; i++) {+ for (j=0; j<4; j++) {+ c->element[i][j] = 0;+ for (k=0; k<4; k++) c->element[i][j] += + a->element[i][k] * b->element[k][j];+ }+ }+ return(c);+ }++/* binary greatest common divisor by Silver and Terzian. See Knuth */+/* both inputs must be >= 0 */+gcd(u, v)+int u, v;+{+int t, f;+ if ((u<0) || (v<0)) return(1); /* error if u<0 or v<0 */+ f = 1;+ while ((0 == (u%2)) && (0 == (v%2))) {+ u>>=1; v>>=1, f*=2;+ }+ if (u&01) { t = -v; goto B4; } else { t = u; }+ B3: if (t > 0) { t >>= 1; } else { t = -((-t) >> 1); }+ B4: if (0 == (t%2)) goto B3;++ if (t > 0) u = t; else v = -t;+ if (0 != (t = u - v)) goto B3;+ return(u*f);+ } ++/***********************/+/* Useful Routines */+/***********************/++/* return roots of ax^2+bx+c */+/* stable algebra derived from Numerical Recipes by Press et al.*/+int quadraticRoots(a, b, c, roots)+double a, b, c, *roots;+{+double d, q;+int count = 0;+ d = (b*b)-(4*a*c);+ if (d < 0.0) { *roots = *(roots+1) = 0.0; return(0); }+ q = -0.5 * (b + (SGN(b)*sqrt(d)));+ if (a != 0.0) { *roots++ = q/a; count++; }+ if (q != 0.0) { *roots++ = c/q; count++; }+ return(count);+ }+++/* generic 1d regula-falsi step. f is function to evaluate */+/* interval known to contain root is given in left, right */+/* returns new estimate */+double RegulaFalsi(f, left, right)+double (*f)(), left, right;+{+double d = (*f)(right) - (*f)(left);+ if (d != 0.0) return (right - (*f)(right)*(right-left)/d);+ return((left+right)/2.0);+ }++/* generic 1d Newton-Raphson step. f is function, df is derivative */+/* x is current best guess for root location. Returns new estimate */+double NewtonRaphson(f, df, x)+double (*f)(), (*df)(), x;+{+double d = (*df)(x);+ if (d != 0.0) return (x-((*f)(x)/d));+ return(x-1.0);+ }+++/* hybrid 1d Newton-Raphson/Regula Falsi root finder. */+/* input function f and its derivative df, an interval */+/* left, right known to contain the root, and an error tolerance */+/* Based on Blinn */+double findroot(left, right, tolerance, f, df)+double left, right, tolerance;+double (*f)(), (*df)();+{+double newx = left;+ while (ABS((*f)(newx)) > tolerance) {+ newx = NewtonRaphson(f, df, newx);+ if (newx < left || newx > right) + newx = RegulaFalsi(f, left, right);+ if ((*f)(newx) * (*f)(left) <= 0.0) right = newx; + else left = newx;+ }+ return(newx);+ } new file mode 100644@@ -0,0 +1,170 @@+/* + * GraphicsGems.h + * Version 1.0 - Andrew Glassner+ * from "Graphics Gems", Academic Press, 1990+ */++/*+ * EULA: The Graphics Gems code is copyright-protected. In other words, + * you cannot claim the text of the code as your own and resell it. + * Using the code is permitted in any program, product, or library, + * non-commercial or commercial. Giving credit is not required, though + * is a nice gesture. The code comes as-is, and if there are any flaws + * or problems with any Gems code, nobody involved with Gems - authors,+ * editors, publishers, or webmasters - are to be held responsible. + * Basically, don't be a jerk, and remember that anything free comes + * with no guarantee. + */+ +#ifndef GG_H++#define GG_H 1++/*********************/+/* 2d geometry types */+/*********************/++typedef struct Point2Struct { /* 2d point */+ double x, y;+ } Point2;+typedef Point2 Vector2;++typedef struct IntPoint2Struct { /* 2d integer point */+ int x, y;+ } IntPoint2;++typedef struct Matrix3Struct { /* 3-by-3 matrix */+ double element[3][3];+ } Matrix3;++typedef struct Box2dStruct { /* 2d box */+ Point2 min, max;+ } Box2;+ ++/*********************/+/* 3d geometry types */+/*********************/++typedef struct Point3Struct { /* 3d point */+ double x, y, z;+ } Point3;+typedef Point3 Vector3;++typedef struct IntPoint3Struct { /* 3d integer point */+ int x, y, z;+ } IntPoint3;+++typedef struct Matrix4Struct { /* 4-by-4 matrix */+ double element[4][4];+ } Matrix4;++typedef struct Box3dStruct { /* 3d box */+ Point3 min, max;+ } Box3;++++/***********************/+/* one-argument macros */+/***********************/++/* absolute value of a */+#define ABS(a) (((a)<0) ? -(a) : (a))++/* round a to nearest int */+#define ROUND(a) ((a)>0 ? (int)((a)+0.5) : -(int)(0.5-(a)))++/* take sign of a, either -1, 0, or 1 */+#define ZSGN(a) (((a)<0) ? -1 : (a)>0 ? 1 : 0) ++/* take binary sign of a, either -1, or 1 if >= 0 */+#define SGN(a) (((a)<0) ? -1 : 1)++/* shout if something that should be true isn't */+#define ASSERT(x) \+if (!(x)) fprintf(stderr," Assert failed: x\n");++/* square a */+#define SQR(a) ((a)*(a)) +++/***********************/+/* two-argument macros */+/***********************/++/* find minimum of a and b */+#define MIN(a,b) (((a)<(b))?(a):(b)) ++/* find maximum of a and b */+#define MAX(a,b) (((a)>(b))?(a):(b)) ++/* swap a and b (see Gem by Wyvill) */+#define SWAP(a,b) { a^=b; b^=a; a^=b; }++/* linear interpolation from l (when a=0) to h (when a=1)*/+/* (equal to (a*h)+((1-a)*l) */+#define LERP(a,l,h) ((l)+(((h)-(l))*(a)))++/* clamp the input to the specified range */+#define CLAMP(v,l,h) ((v)<(l) ? (l) : (v) > (h) ? (h) : v)+++/****************************/+/* memory allocation macros */+/****************************/++/* create a new instance of a structure (see Gem by Hultquist) */+#define NEWSTRUCT(x) (struct x *)(malloc((unsigned)sizeof(struct x)))++/* create a new instance of a type */+#define NEWTYPE(x) (x *)(malloc((unsigned)sizeof(x)))+++/********************/+/* useful constants */+/********************/++#define PI 3.141592 /* the venerable pi */+#define PITIMES2 6.283185 /* 2 * pi */+#define PIOVER2 1.570796 /* pi / 2 */+#define E 2.718282 /* the venerable e */+#define SQRT2 1.414214 /* sqrt(2) */+#define SQRT3 1.732051 /* sqrt(3) */+#define GOLDEN 1.618034 /* the golden ratio */+#define DTOR 0.017453 /* convert degrees to radians */+#define RTOD 57.29578 /* convert radians to degrees */+++/************/+/* booleans */+/************/++#define TRUE 1+#define FALSE 0+#define ON 1+#define OFF 0+typedef int boolean; /* boolean data type */+typedef boolean flag; /* flag data type */++extern double V2SquaredLength(), V2Length();+extern double V2Dot(), V2DistanceBetween2Points(); +extern Vector2 *V2Negate(), *V2Normalize(), *V2Scale(), *V2Add(), *V2Sub();+extern Vector2 *V2Lerp(), *V2Combine(), *V2Mul(), *V2MakePerpendicular();+extern Vector2 *V2New(), *V2Duplicate();+extern Point2 *V2MulPointByMatrix();+extern Matrix3 *V2MatMul();++extern double V3SquaredLength(), V3Length();+extern double V3Dot(), V3DistanceBetween2Points();+extern Vector3 *V3Normalize(), *V3Scale(), *V3Add(), *V3Sub();+extern Vector3 *V3Lerp(), *V3Combine(), *V3Mul(), *V3Cross();+extern Vector3 *V3New(), *V3Duplicate();+extern Point3 *V3MulPointByMatrix();+extern Matrix4 *V3MatMul();++extern double RegulaFalsi(), NewtonRaphson(), findroot();++#endif+new file mode 100644@@ -0,0 +1,570 @@+/*+An Algorithm for Automatically Fitting Digitized Curves+by Philip J. Schneider+from "Graphics Gems", Academic Press, 1990+Adapted to BirdFont by Johan Mattsson 2015+*/++/*+ * EULA: The Graphics Gems code is copyright-protected. In other words, + * you cannot claim the text of the code as your own and resell it. + * Using the code is permitted in any program, product, or library, + * non-commercial or commercial. Giving credit is not required, though + * is a nice gesture. The code comes as-is, and if there are any flaws + * or problems with any Gems code, nobody involved with Gems - authors,+ * editors, publishers, or webmasters - are to be held responsible. + * Basically, don't be a jerk, and remember that anything free comes + * with no guarantee. + */++/* fit_cubic.c */ +/* Piecewise cubic fitting code */++#include "GraphicsGems.h" +#include "birdfont.h"+#include <stdio.h>+#include <malloc.h>+#include <math.h>+#include <stdio.h>++typedef Point2 *BezierCurve;++/* Forward declarations */+void FitCurve();+static void FitCubic();+static double *Reparameterize();+static double NewtonRaphsonRootFind();+static Point2 BezierII();+static double B0(), B1(), B2(), B3();+static Vector2 ComputeLeftTangent();+static Vector2 ComputeRightTangent();+static Vector2 ComputeCenterTangent();+static double ComputeMaxError();+static double *ChordLengthParameterize();+static BezierCurve GenerateBezier();+static Vector2 V2AddII();+static Vector2 V2ScaleIII();+static Vector2 V2SubII();++#define MAXPOINTS 1000 /* The most points you can have */++BirdFontPath* simplified_path = NULL;++BirdFontPath* fit_bezier_curve_to_line (BirdFontPath* path, int start, int stop, double error) {+ int i, j, n;+ Point2* points;+ BirdFontPath* result; ++ if (simplified_path != NULL) {+ printf ("Simplifier is not thread safe.");+ return bird_font_path_new ();+ }++ n = stop - start + 1;+ + if (n > bird_font_path_size (path) || n < 0) {+ printf ("Index of of bounds in simplifier");+ return bird_font_path_new ();+ }+ + points = malloc (n * sizeof (Point2)); + + j = 0;+ for (i = start; i <= stop; i++) {+ points[j].x = bird_font_path_get_x_at (path, i);+ points[j].y = bird_font_path_get_y_at (path, i);+ j++;+ }+ + simplified_path = bird_font_path_new ();+ FitCurve(points, n, error);+ free (points);++ result = simplified_path;+ simplified_path = NULL;+ return result;+}++DrawBezierCurve(int n, Point2* curve) {+ if (n == 3) {+ bird_font_path_add_cubic_bezier_points (simplified_path,+ curve[0].x, curve[0].y, + curve[1].x, curve[1].y, + curve[2].x, curve[2].y,+ curve[3].x, curve[3].y);+ } else {+ printf ("Expecting three points\n");+ }+}++/*+ * FitCurve :+ * Fit a Bezier curve to a set of digitized points + */+void FitCurve(d, nPts, error)+ Point2 *d; /* Array of digitized points */+ int nPts; /* Number of digitized points */+ double error; /* User-defined error squared */+{+ Vector2 tHat1, tHat2; /* Unit tangent vectors at endpoints */++ tHat1 = ComputeLeftTangent(d, 0);+ tHat2 = ComputeRightTangent(d, nPts - 1);+ FitCubic(d, 0, nPts - 1, tHat1, tHat2, error);+}++/*+ * FitCubic :+ * Fit a Bezier curve to a (sub)set of digitized points+ */+static void FitCubic(d, first, last, tHat1, tHat2, error)+ Point2 *d; /* Array of digitized points */+ int first, last; /* Indices of first and last pts in region */+ Vector2 tHat1, tHat2; /* Unit tangent vectors at endpoints */+ double error; /* User-defined error squared */+{+ BezierCurve bezCurve; /*Control points of fitted Bezier curve*/+ double *u; /* Parameter values for point */+ double *uPrime; /* Improved parameter values */+ double maxError; /* Maximum fitting error */+ int splitPoint; /* Point to split point set at */+ int nPts; /* Number of points in subset */+ double iterationError; /*Error below which you try iterating */+ int maxIterations = 4; /* Max times to try iterating */+ Vector2 tHatCenter; /* Unit tangent vector at splitPoint */+ int i; ++ iterationError = error * error;+ nPts = last - first + 1;++ /* Use heuristic if region only has two points in it */+ if (nPts == 2) {+ double dist = V2DistanceBetween2Points(&d[last], &d[first]) / 3.0;++ bezCurve = (Point2 *)malloc(4 * sizeof(Point2));+ bezCurve[0] = d[first];+ bezCurve[3] = d[last];+ V2Add(&bezCurve[0], V2Scale(&tHat1, dist), &bezCurve[1]);+ V2Add(&bezCurve[3], V2Scale(&tHat2, dist), &bezCurve[2]);+ DrawBezierCurve(3, bezCurve);+ free((void *)bezCurve);+ return;+ }++ /* Parameterize points, and attempt to fit curve */+ u = ChordLengthParameterize(d, first, last);+ bezCurve = GenerateBezier(d, first, last, u, tHat1, tHat2);++ /* Find max deviation of points to fitted curve */+ maxError = ComputeMaxError(d, first, last, bezCurve, u, &splitPoint);+ if (maxError < error) {+ DrawBezierCurve(3, bezCurve);+ free((void *)u);+ free((void *)bezCurve);+ return;+ }+++ /* If error not too large, try some reparameterization */+ /* and iteration */+ if (maxError < iterationError) {+ for (i = 0; i < maxIterations; i++) {+ uPrime = Reparameterize(d, first, last, u, bezCurve);+ free((void *)bezCurve);+ bezCurve = GenerateBezier(d, first, last, uPrime, tHat1, tHat2);+ maxError = ComputeMaxError(d, first, last,+ bezCurve, uPrime, &splitPoint);+ if (maxError < error) {+ DrawBezierCurve(3, bezCurve);+ free((void *)u);+ free((void *)bezCurve);+ free((void *)uPrime);+ return;+ }+ free((void *)u);+ u = uPrime;+ }+ }++ /* Fitting failed -- split at max error point and fit recursively */+ free((void *)u);+ free((void *)bezCurve);+ tHatCenter = ComputeCenterTangent(d, splitPoint);+ FitCubic(d, first, splitPoint, tHat1, tHatCenter, error);+ V2Negate(&tHatCenter);+ FitCubic(d, splitPoint, last, tHatCenter, tHat2, error);+}+++/*+ * GenerateBezier :+ * Use least-squares method to find Bezier control points for region.+ *+ */+static BezierCurve GenerateBezier(d, first, last, uPrime, tHat1, tHat2)+ Point2 *d; /* Array of digitized points */+ int first, last; /* Indices defining region */+ double *uPrime; /* Parameter values for region */+ Vector2 tHat1, tHat2; /* Unit tangents at endpoints */+{+ int i;+ Vector2 A[MAXPOINTS][2]; /* Precomputed rhs for eqn */+ int nPts; /* Number of pts in sub-curve */+ double C[2][2]; /* Matrix C */+ double X[2]; /* Matrix X */+ double det_C0_C1, /* Determinants of matrices */+ det_C0_X,+ det_X_C1;+ double alpha_l, /* Alpha values, left and right */+ alpha_r;+ Vector2 tmp; /* Utility variable */+ BezierCurve bezCurve; /* RETURN bezier curve ctl pts */++ bezCurve = (Point2 *)malloc(4 * sizeof(Point2));+ nPts = last - first + 1;++ + /* Compute the A's */+ for (i = 0; i < nPts; i++) {+ Vector2 v1, v2;+ v1 = tHat1;+ v2 = tHat2;+ V2Scale(&v1, B1(uPrime[i]));+ V2Scale(&v2, B2(uPrime[i]));+ A[i][0] = v1;+ A[i][1] = v2;+ }++ /* Create the C and X matrices */+ C[0][0] = 0.0;+ C[0][1] = 0.0;+ C[1][0] = 0.0;+ C[1][1] = 0.0;+ X[0] = 0.0;+ X[1] = 0.0;++ for (i = 0; i < nPts; i++) {+ C[0][0] += V2Dot(&A[i][0], &A[i][0]);+ C[0][1] += V2Dot(&A[i][0], &A[i][1]);+/* C[1][0] += V2Dot(&A[i][0], &A[i][1]);*/ + C[1][0] = C[0][1];+ C[1][1] += V2Dot(&A[i][1], &A[i][1]);++ tmp = V2SubII(d[first + i],+ V2AddII(+ V2ScaleIII(d[first], B0(uPrime[i])),+ V2AddII(+ V2ScaleIII(d[first], B1(uPrime[i])),+ V2AddII(+ V2ScaleIII(d[last], B2(uPrime[i])),+ V2ScaleIII(d[last], B3(uPrime[i]))))));+ ++ X[0] += V2Dot(&A[i][0], &tmp);+ X[1] += V2Dot(&A[i][1], &tmp);+ }++ /* Compute the determinants of C and X */+ det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];+ det_C0_X = C[0][0] * X[1] - C[1][0] * X[0];+ det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1];++ /* Finally, derive alpha values */+ alpha_l = (det_C0_C1 == 0) ? 0.0 : det_X_C1 / det_C0_C1;+ alpha_r = (det_C0_C1 == 0) ? 0.0 : det_C0_X / det_C0_C1;++ /* If alpha negative, use the Wu/Barsky heuristic (see text) */+ /* (if alpha is 0, you get coincident control points that lead to+ * divide by zero in any subsequent NewtonRaphsonRootFind() call. */+ double segLength = V2DistanceBetween2Points(&d[last], &d[first]);+ double epsilon = 1.0e-6 * segLength;+ if (alpha_l < epsilon || alpha_r < epsilon)+ {+ /* fall back on standard (probably inaccurate) formula, and subdivide further if needed. */+ double dist = segLength / 3.0;+ bezCurve[0] = d[first];+ bezCurve[3] = d[last];+ V2Add(&bezCurve[0], V2Scale(&tHat1, dist), &bezCurve[1]);+ V2Add(&bezCurve[3], V2Scale(&tHat2, dist), &bezCurve[2]);+ return (bezCurve);+ }++ /* First and last control points of the Bezier curve are */+ /* positioned exactly at the first and last data points */+ /* Control points 1 and 2 are positioned an alpha distance out */+ /* on the tangent vectors, left and right, respectively */+ bezCurve[0] = d[first];+ bezCurve[3] = d[last];+ V2Add(&bezCurve[0], V2Scale(&tHat1, alpha_l), &bezCurve[1]);+ V2Add(&bezCurve[3], V2Scale(&tHat2, alpha_r), &bezCurve[2]);+ return (bezCurve);+}+++/*+ * Reparameterize:+ * Given set of points and their parameterization, try to find+ * a better parameterization.+ *+ */+static double *Reparameterize(d, first, last, u, bezCurve)+ Point2 *d; /* Array of digitized points */+ int first, last; /* Indices defining region */+ double *u; /* Current parameter values */+ BezierCurve bezCurve; /* Current fitted curve */+{+ int nPts = last-first+1; + int i;+ double *uPrime; /* New parameter values */++ uPrime = (double *)malloc(nPts * sizeof(double));+ for (i = first; i <= last; i++) {+ uPrime[i-first] = NewtonRaphsonRootFind(bezCurve, d[i], u[i-+ first]);+ }+ return (uPrime);+}++++/*+ * NewtonRaphsonRootFind :+ * Use Newton-Raphson iteration to find better root.+ */+static double NewtonRaphsonRootFind(Q, P, u)+ BezierCurve Q; /* Current fitted curve */+ Point2 P; /* Digitized point */+ double u; /* Parameter value for "P" */+{+ double numerator, denominator;+ Point2 Q1[3], Q2[2]; /* Q' and Q'' */+ Point2 Q_u, Q1_u, Q2_u; /*u evaluated at Q, Q', & Q'' */+ double uPrime; /* Improved u */+ int i;+ + /* Compute Q(u) */+ Q_u = BezierII(3, Q, u);+ + /* Generate control vertices for Q' */+ for (i = 0; i <= 2; i++) {+ Q1[i].x = (Q[i+1].x - Q[i].x) * 3.0;+ Q1[i].y = (Q[i+1].y - Q[i].y) * 3.0;+ }+ + /* Generate control vertices for Q'' */+ for (i = 0; i <= 1; i++) {+ Q2[i].x = (Q1[i+1].x - Q1[i].x) * 2.0;+ Q2[i].y = (Q1[i+1].y - Q1[i].y) * 2.0;+ }+ + /* Compute Q'(u) and Q''(u) */+ Q1_u = BezierII(2, Q1, u);+ Q2_u = BezierII(1, Q2, u);+ + /* Compute f(u)/f'(u) */+ numerator = (Q_u.x - P.x) * (Q1_u.x) + (Q_u.y - P.y) * (Q1_u.y);+ denominator = (Q1_u.x) * (Q1_u.x) + (Q1_u.y) * (Q1_u.y) ++ (Q_u.x - P.x) * (Q2_u.x) + (Q_u.y - P.y) * (Q2_u.y);+ if (denominator == 0.0f) return u;++ /* u = u - f(u)/f'(u) */+ uPrime = u - (numerator/denominator);+ return (uPrime);+}++ + +/*+ * Bezier :+ * Evaluate a Bezier curve at a particular parameter value+ * + */+static Point2 BezierII(degree, V, t)+ int degree; /* The degree of the bezier curve */+ Point2 *V; /* Array of control points */+ double t; /* Parametric value to find point for */+{+ int i, j; + Point2 Q; /* Point on curve at parameter t */+ Point2 *Vtemp; /* Local copy of control points */++ /* Copy array */+ Vtemp = (Point2 *)malloc((unsigned)((degree+1) + * sizeof (Point2)));+ for (i = 0; i <= degree; i++) {+ Vtemp[i] = V[i];+ }++ /* Triangle computation */+ for (i = 1; i <= degree; i++) { + for (j = 0; j <= degree-i; j++) {+ Vtemp[j].x = (1.0 - t) * Vtemp[j].x + t * Vtemp[j+1].x;+ Vtemp[j].y = (1.0 - t) * Vtemp[j].y + t * Vtemp[j+1].y;+ }+ }++ Q = Vtemp[0];+ free((void *)Vtemp);+ return Q;+}+++/*+ * B0, B1, B2, B3 :+ * Bezier multipliers+ */+static double B0(u)+ double u;+{+ double tmp = 1.0 - u;+ return (tmp * tmp * tmp);+}+++static double B1(u)+ double u;+{+ double tmp = 1.0 - u;+ return (3 * u * (tmp * tmp));+}++static double B2(u)+ double u;+{+ double tmp = 1.0 - u;+ return (3 * u * u * tmp);+}++static double B3(u)+ double u;+{+ return (u * u * u);+}++++/*+ * ComputeLeftTangent, ComputeRightTangent, ComputeCenterTangent :+ *Approximate unit tangents at endpoints and "center" of digitized curve+ */+static Vector2 ComputeLeftTangent(d, end)+ Point2 *d; /* Digitized points*/+ int end; /* Index to "left" end of region */+{+ Vector2 tHat1;+ tHat1 = V2SubII(d[end+1], d[end]);+ tHat1 = *V2Normalize(&tHat1);+ return tHat1;+}++static Vector2 ComputeRightTangent(d, end)+ Point2 *d; /* Digitized points */+ int end; /* Index to "right" end of region */+{+ Vector2 tHat2;+ tHat2 = V2SubII(d[end-1], d[end]);+ tHat2 = *V2Normalize(&tHat2);+ return tHat2;+}+++static Vector2 ComputeCenterTangent(d, center)+ Point2 *d; /* Digitized points */+ int center; /* Index to point inside region */+{+ Vector2 V1, V2, tHatCenter;++ V1 = V2SubII(d[center-1], d[center]);+ V2 = V2SubII(d[center], d[center+1]);+ tHatCenter.x = (V1.x + V2.x)/2.0;+ tHatCenter.y = (V1.y + V2.y)/2.0;+ tHatCenter = *V2Normalize(&tHatCenter);+ return tHatCenter;+}+++/*+ * ChordLengthParameterize :+ * Assign parameter values to digitized points + * using relative distances between points.+ */+static double *ChordLengthParameterize(d, first, last)+ Point2 *d; /* Array of digitized points */+ int first, last; /* Indices defining region */+{+ int i; + double *u; /* Parameterization */++ u = (double *)malloc((unsigned)(last-first+1) * sizeof(double));++ u[0] = 0.0;+ for (i = first+1; i <= last; i++) {+ u[i-first] = u[i-first-1] ++ V2DistanceBetween2Points(&d[i], &d[i-1]);+ }++ for (i = first + 1; i <= last; i++) {+ u[i-first] = u[i-first] / u[last-first];+ }++ return(u);+}+++++/*+ * ComputeMaxError :+ * Find the maximum squared distance of digitized points+ * to fitted curve.+*/+static double ComputeMaxError(d, first, last, bezCurve, u, splitPoint)+ Point2 *d; /* Array of digitized points */+ int first, last; /* Indices defining region */+ BezierCurve bezCurve; /* Fitted Bezier curve */+ double *u; /* Parameterization of points */+ int *splitPoint; /* Point of maximum error */+{+ int i;+ double maxDist; /* Maximum error */+ double dist; /* Current error */+ Point2 P; /* Point on curve */+ Vector2 v; /* Vector from point to curve */++ *splitPoint = (last - first + 1)/2;+ maxDist = 0.0;+ for (i = first + 1; i < last; i++) {+ P = BezierII(3, bezCurve, u[i-first]);+ v = V2SubII(P, d[i]);+ dist = V2SquaredLength(&v);+ if (dist >= maxDist) {+ maxDist = dist;+ *splitPoint = i;+ }+ }+ return (maxDist);+}+static Vector2 V2AddII(a, b)+ Vector2 a, b;+{+ Vector2 c;+ c.x = a.x + b.x; c.y = a.y + b.y;+ return (c);+}+static Vector2 V2ScaleIII(v, s)+ Vector2 v;+ double s;+{+ Vector2 result;+ result.x = v.x * s; result.y = v.y * s;+ return (result);+}++static Vector2 V2SubII(a, b)+ Vector2 a, b;+{+ Vector2 c;+ c.x = a.x - b.x; c.y = a.y - b.y;+ return (c);+}+@@ -17,6 +17,9 @@ using Math; namespace BirdFont { +[CCode (cname = "fit_bezier_curve_to_line")]+public extern static Path fit_bezier_curve_to_line (Path lines, int start, int stop, double error);+ public enum Direction { CLOCKWISE, COUNTER_CLOCKWISE@@ -2342,6 +2345,45 @@ public class Path { fast_stroke = StrokeTool.get_stroke_fast (this, stroke); return (!) fast_stroke; }+ + // Callback for path simplifier+ public void add_cubic_bezier_points (double x0, double y0, double x1, double y1,+ double x2, double y2, double x3, double y3) {+ + EditPoint start;+ EditPoint end;++ if (points.size > 0) {+ start = get_last_point ();+ } else {+ start = add (x0, y0);+ }+ + end = add (x3, y3);+ + start.set_point_type (PointType.CUBIC);+ end.set_point_type (PointType.CUBIC);+ + start.convert_to_curve ();+ end.convert_to_curve ();+ + start.get_right_handle ().move_to_coordinate (x1, y1);+ end.get_left_handle ().move_to_coordinate (x2, y2);+ }+ + public int size () {+ return points.size;+ }+ + public double get_x_at (int i) {+ return_val_if_fail (0 <= i < points.size, 0);+ return points.get (i).x;+ }++ public double get_y_at (int i) {+ return_val_if_fail (0 <= i < points.size, 0);+ return points.get (i).y;+ } } }@@ -1237,12 +1237,12 @@ public class PenTool : Tool { public static Path merge_open_paths (Path path1, Path path2) { Path union, merge; EditPoint first_point;- - return_if_fail (path1.points.size > 1);- return_if_fail (path2.points.size > 1); union = path1.copy (); merge = path2.copy ();+ + return_val_if_fail (path1.points.size >= 1, merge);+ return_val_if_fail (path2.points.size >= 1, union); merge.points.get (0).set_tie_handle (false); merge.points.get (0).set_reflective_handles (false);@@ -99,42 +99,198 @@ public class StrokeTool : Tool { } public static PathList get_stroke (Path path, double thickness) {- PathList o;+ PathList o, m; o = get_stroke_fast (path, thickness); o = get_all_parts (o); o = merge (o); + m = new PathList (); foreach (Path p in o.paths) {- simplify_stroke (p);+ m.add (simplify_stroke (p)); } - return o;+ return m; } - static void simplify_stroke (Path p) {- Gee.ArrayList<PointSelection> points = new Gee.ArrayList<PointSelection> ();- double error;+ static Path simplify_stroke (Path p) {+ Path simplified = new Path ();+ Path segment;+ EditPoint ep, ep_start, last, first;+ int start, stop;+ int j; - //p.remove_points_on_points (0.1); // FIXME:+ foreach (EditPoint e in p.points) {+ PenTool.convert_point_type (e, PointType.CUBIC);+ } - foreach (EditPoint ep in p.points) {+ foreach (EditPoint e in p.points) {+ if ((e.flags & EditPoint.CURVE) == 0) {+ p.set_new_start (e);+ break;+ }+ }+ + for (int i = 0; i < p.points.size; i++) {+ ep = p.points.get (i);+ if ((ep.flags & EditPoint.CURVE) > 0) {- points.add (new PointSelection (ep, p));+ start = i;+ ep_start = ep;+ for (j = start + 1; j < p.points.size; j++) {+ ep = p.points.get (j);+ if ((ep.flags & EditPoint.CURVE) == 0) {+ break;+ }+ }+ + stop = j;+ start -= 1;+ + if (start < 0) {+ warning ("start < 0");+ start = 0;+ }++ if (stop >= p.points.size) {+ warning ("stop >= p.points.size");+ stop = p.points.size - 1;+ }+ + if (ep_start.get_right_handle ().is_line () && ep.get_left_handle ().is_line ()) {+ simplified.add_point (ep_start.copy ());+ simplified.add_point (ep.copy ());+ continue;+ }+ + segment = fit_bezier_curve_to_line (p, start, stop, 0.0007);++ segment.get_first_point ().color = Color.yellow ();+ segment.get_last_point ().color = Color.pink ();+ + if (stroke_selected) { // FIXME: DELETE + ((!) BirdFont.get_current_font ().get_glyph ("l")).add_path (segment.copy ());+ }++ last = simplified.get_last_point ();+ first = segment.get_first_point ();+ + if (segment.get_last_point ().get_right_handle ().y > 140) {+ print (segment.get_last_point ().to_string ());+ print ("RIGHT\n");+ }+ + if (segment.get_last_point ().get_left_handle ().y > 140) {+ print (segment.get_last_point ().to_string ());+ print ("LEFT\n");+ }++ if (segment.get_first_point ().get_right_handle ().y > 140) {+ print (segment.get_first_point ().to_string ());+ print ("FRIGHT\n");+ }+ + if (segment.get_first_point ().get_left_handle ().y > 140) {+ print (segment.get_first_point ().to_string ());+ print ("FLEFT\n");+ }+ + if (simplified.points.size > 1) {+ simplified.delete_last_point ();+ }+ + if (ep_start.get_right_handle ().is_line ()) {+ segment.get_first_point ().get_right_handle ().convert_to_line ();+ segment.get_first_point ().recalculate_linear_handles ();+ }+ + first.get_left_handle ().convert_to_curve ();+ first.get_left_handle ().x = last.get_left_handle ().x;+ first.get_left_handle ().y = last.get_left_handle ().y;+ + if (first.get_left_handle ().y > 140) {+ print (first.to_string ());+ print ("LEFT2\n");+ }+ + simplified.append_path (segment.copy ());+ + last = simplified.get_last_point ();+ last.get_right_handle ().convert_to_line ();+ last.recalculate_linear_handles ();+ + i = stop;+ } else {+ if (ep.get_left_handle ().y > 140) {+ print (ep.to_string ());+ print ("LEFT3\n");+ }+ if (ep.get_right_handle ().y > 140) {+ print (ep.to_string ());+ print ("RIFGHIR 4\n");+ }+ + simplified.add_point (ep.copy ()); } } - if (stroke_selected) { // FIXME: DELETE - ((!) BirdFont.get_current_font ().get_glyph ("k")).add_path (p.copy ());+ simplified.close ();+ simplified.remove_points_on_points ();+ simplified.recalculate_linear_handles ();+ + if (p.get_first_point ().tie_handles) {+ first = simplified.get_first_point ();+ first.convert_to_curve ();+ first.tie_handles = true;+ first.process_tied_handle (); }- - foreach (PointSelection ps in points) {- error = PenTool.remove_point_simplify_path_fast (ps, 0.6, 0.8);++ if (p.get_last_point ().tie_handles) {+ last = simplified.get_last_point ();+ last.convert_to_curve ();+ last.tie_handles = true;+ last.process_tied_handle (); } - p.update_region_boundaries ();+ simplified.get_first_point ().color = Color.green ();+ simplified.get_last_point ().color = Color.blue ();++ foreach (EditPoint e in simplified.points) {+ if (e.get_right_handle ().y > 140) {+ print (e.to_string ());+ print ("R1 4\n");+ }++ if (e.get_left_handle ().y > 140) {+ print (e.to_string ());+ print ("L1 4\n");+ }+ }++ return simplified; } + static bool segment_is_line (Path p, int start, int stop) {+ EditPoint p0, p1, p2;+ + if (stop - start < 2) {+ warning ("Too small segment.");+ return true;+ }+ + for (int i = start; i <= stop - 2; i++) {+ p0 = p.points.get (i);+ p1 = p.points.get (i + 1);+ p2 = p.points.get (i + 2);+ + if (!is_line (p0.x, p0.y, p1.x, p1.y, p2.x, p2.y)) {+ return false;+ }+ }+ + return true;+ }+ /** Create one stroke from the outline and counter stroke and close the * open endings. * @@ -74,6 +74,7 @@ class Vala(object): self.vala = get_sources_path (src, '*.vala') self.c = get_sources_path (src, '*.c') # copy regular c sources+ self.c += get_sources_path (src, '*.h') # copy header files for the c sources self.cc = [join(build + '/' + src, f) for f in get_sources_name (src, '*.c') ] self.cc += [join(build + '/' + src, f.replace('.vala', '.c')) for f in get_sources_name (src, '*.vala')] self.obj = [self.build + '/' + self.src + '/' + f.replace('.c', '.o') for f in get_sources_name (src, '*.c')]@@ -111,7 +112,7 @@ class Vala(object): for f in self.c: yield {- 'name': 'copy_c',+ 'name': 'copy_c_' + f, 'actions': [ 'mkdir -p '+ self.build + '/' + self.src + '/', 'cp ' + f + ' ' + self.build + '/' + self.src + '/'