The Birdfont Source Code

--- a/dodo.py +++ b/dodo.py @@ -1,5 +1,5 @@ """ - Copyright (C) 2012, 2013, 2014 Eduardo Naufel Schettino and Johan Mattsson + Copyright (C) 2012, 2013, 2014 2015 Eduardo Naufel Schettino and Johan Mattsson This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -30,6 +30,7 @@ DOIT_CONFIG = { 'default_tasks': [ 'build', + 'libbirdgems', 'libbirdxml', 'libbirdfont', 'birdfont', @@ -62,6 +63,10 @@ 'posix', 'posixtypes' ] + + LIBBIRD_LIBS = [ + 'glib-2.0' + ] def task_build (): if not os.path.exists ("build/configured"): @@ -108,39 +113,44 @@ yield libbirdxml.gen_ln() - libbird = Vala(src='libbirdfont', build='build', library='birdfont', so_version=version.SO_VERSION, pkg_libs=LIBS, vala_deps=[libbirdxml]) + libbirdgems = Vala(src='libbirdgems', build='build', library='birdgems', so_version=version.LIBBIRDGEMS_SO_VERSION, pkg_libs=LIBBIRD_LIBS, vala_deps=[]) + def task_libbirdgems(): + yield libbirdgems.gen_c(valac_options) + yield libbirdgems.gen_o([]) + yield libbirdgems.gen_so('-L ./build -l m') + yield libbirdgems.gen_ln() + + + libbird = Vala(src='libbirdfont', build='build', library='birdfont', so_version=version.SO_VERSION, pkg_libs=LIBS, vala_deps=[libbirdgems, libbirdxml]) def task_libbirdfont(): yield libbird.gen_c(valac_options) yield libbird.gen_o(['-fPIC -I./build/', """-D 'GETTEXT_PACKAGE="birdfont"'"""]) - yield libbird.gen_so('-L ./build -l birdxml') + yield libbird.gen_so('-L ./build -l birdxml -L ./build -l birdgems') yield libbird.gen_ln() def task_birdfont(): - bird = Vala(src='birdfont', build='build', pkg_libs=LIBS, vala_deps=[libbird, libbirdxml]) + bird = Vala(src='birdfont', build='build', pkg_libs=LIBS, vala_deps=[libbird, libbirdxml, libbirdgems]) yield bird.gen_c(valac_options) yield bird.gen_bin(["""-D 'GETTEXT_PACKAGE="birdfont"' """]) def task_birdfont_autotrace(): - exp = Vala(src='birdfont-autotrace', build='build', pkg_libs=LIBS, vala_deps=[libbird, libbirdxml]) + exp = Vala(src='birdfont-autotrace', build='build', pkg_libs=LIBS, vala_deps=[libbird, libbirdxml, libbirdgems]) yield exp.gen_c(valac_options) yield exp.gen_bin(["""-D 'GETTEXT_PACKAGE="birdfont"' """]) def task_birdfont_export(): - exp = Vala(src='birdfont-export', build='build', pkg_libs=LIBS, vala_deps=[libbird, libbirdxml]) + exp = Vala(src='birdfont-export', build='build', pkg_libs=LIBS, vala_deps=[libbird, libbirdxml, libbirdgems]) yield exp.gen_c(valac_options) yield exp.gen_bin(["""-D 'GETTEXT_PACKAGE="birdfont"' """]) def task_birdfont_import(): - exp = Vala(src='birdfont-import', build='build', pkg_libs=LIBS, vala_deps=[libbird, libbirdxml]) + exp = Vala(src='birdfont-import', build='build', pkg_libs=LIBS, vala_deps=[libbird, libbirdxml, libbirdgems]) yield exp.gen_c(valac_options) yield exp.gen_bin(["""-D 'GETTEXT_PACKAGE="birdfont"' """]) - - - def task_compile_translations ():

--- a/install.py +++ b/install.py @@ -2,5 +2,5 @@ """ - Copyright (C) 2013 2014 Johan Mattsson + Copyright (C) 2013 2014 2015 Johan Mattsson This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -168,6 +168,24 @@ link (libdir, 'libbirdxml.' + version.LIBBIRDXML_SO_VERSION + '.dylib', ' libbirdxml.dylib') else: print ("Can't find libbirdxml.") + exit (1) + + + if os.path.isfile ('build/bin/libbirdgems.so.' + version.LIBBIRDGEMS_SO_VERSION): + install ('build/bin/libbirdgems.so.' + version.LIBBIRDGEMS_SO_VERSION, libdir, 644) + link (libdir, 'libbirdgems.so.' + version.LIBBIRDGEMS_SO_VERSION, ' libbirdgems.so.' + version.LIBBIRDXML_SO_VERSION_MAJOR) + link (libdir, 'libbirdgems.so.' + version.LIBBIRDGEMS_SO_VERSION, ' libbirdgems.so') + elif os.path.isfile ('build/libbirdgems.so.' + version.LIBBIRDGEMS_SO_VERSION): + install ('build/libbirdgems.so.' + version.LIBBIRDGEMS_SO_VERSION, libdir, 644) + link (libdir, 'libbirdgems.so.' + version.LIBBIRDGEMS_SO_VERSION, ' libbirdgems.so.' + version.LIBBIRDGEMS_SO_VERSION_MAJOR) + link (libdir, 'libbirdgems.so.' + version.LIBBIRDGEMS_SO_VERSION, ' libbirdgems.so') + elif os.path.isfile ('build/bin/libbirdgems.' + version.LIBBIRDGEMS_SO_VERSION + '.dylib'): + install ('build/bin/libbirdgems.' + version.LIBBIRDGEMS_SO_VERSION + '.dylib', libdir, 644) + link (libdir, 'libbirdgems.' + version.LIBBIRDGEMS_SO_VERSION + '.dylib', ' libbirdgems.dylib.' + version.LIBBIRDGEMS_SO_VERSION_MAJOR) + link (libdir, 'libbirdgems.' + version.LIBBIRDGEMS_SO_VERSION + '.dylib', ' libbirdgems.dylib') + else: + print ("Can't find libbirdgems.") + exit (1) #manpages

--- a/libbirdfont/GraphicsGems/GGVecLib.c +++ /dev/null @@ -1,459 +1,1 @@ - /* - 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); - }

--- a/libbirdfont/GraphicsGems/GraphicsGems.h +++ /dev/null @@ -1,171 +1,1 @@ - /* - * 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 -

--- a/libbirdfont/GraphicsGems/fit_cubic.c +++ /dev/null @@ -1,571 +1,1 @@ - /* - 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); - } -

--- a/libbirdfont/Path.vala +++ b/libbirdfont/Path.vala @@ -16,9 +16,6 @@ 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, @@ -2369,22 +2366,8 @@ 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; } } }

--- a/libbirdfont/StrokeTool.vala +++ b/libbirdfont/StrokeTool.vala @@ -157,10 +157,10 @@ stop = p.points.size - 1; } - ep_start = p.points.get (start); + ep_start = p.points.get (start); ep = p.points.get (stop); - segment = fit_bezier_curve_to_line (p, start, stop, 0.0007); + segment = fit_bezier_path (p, start, stop, 0.0007); if (stroke_selected) { // FIXME: DELETE ((!) BirdFont.get_current_font ().get_glyph ("l")).add_path (segment.copy ()); @@ -227,6 +227,49 @@ rh.parent.recalculate_linear_handles (); } + return simplified; + } + + static Path fit_bezier_path (Path p, int start, int stop, double error) { + int index, size; + Path simplified; + double[] lines; + double[] result; + EditPoint ep; + + simplified = new Path (); + + return_val_if_fail (0 <= start < p.points.size, simplified); + return_val_if_fail (0 <= stop < p.points.size, simplified); + + size = stop - start + 1; + lines = new double[2 * size]; + + index = 0; + + for (int i = start; i <= stop; i++) { + ep = p.points.get (i); + lines[index] = ep.x; + index++; + + lines[index] = ep.y; + index++; + } + + return_val_if_fail (2 * size == index, new Path ()); + + Gems.fit_bezier_curve_to_line (lines, error, out result); + + return_val_if_fail (!is_null (result), simplified); + + for (int i = 0; i + 7 < result.length; i += 8) { + simplified.add_cubic_bezier_points ( + result[i], result[i + 1], + result[i + 2], result[i + 3], + result[i + 4], result[i + 5], + result[i + 6], result[i + 7]); + } + return simplified; }

--- /dev/null +++ b/libbirdgems/FitCubic.vala~ @@ -1,1 +1,4 @@ + class FitCubic { + + }

--- /dev/null +++ b/libbirdgems/GGVecLib.c @@ -1,1 +1,459 @@ + /* + 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); + }

--- /dev/null +++ b/libbirdgems/GraphicsGems.h @@ -1,1 +1,171 @@ + /* + * 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 +

--- /dev/null +++ b/libbirdgems/fit_cubic.c @@ -1,1 +1,618 @@ + /* + 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 <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 */ + + int simplified_path_buffer_size = 0; + int simplified_path_size = 0; + double* simplified_path = NULL; + + /** Generates an a bezier path and returns the length of the output array. */ + void fit_bezier_curve_to_line ( + double* lines, + int lines_size, + double error, + double** bezier_path, /** generated Bezier curves */ + int* bezier_path_size) + { + int i, j; + Point2* points; + int npoints; + + if (npoints % 2 != 0) { + fprintf (stderr, "Odd number of coordinates in fit_bezier_curve_to_line."); + return; + } + + if (lines == NULL || lines_size == 0) { + fprintf (stderr, "No lines in fit_bezier_curve_to_line."); + return; + } + + if (simplified_path != NULL) { + fprintf (stderr, "Path simplification is alreading running."); + return; + } + + if (bezier_path == NULL) { + fprintf (stderr, "No destination for output buffer in fit_bezier_curve_to_line"); + return; + } + + if (bezier_path_size == NULL) { + fprintf (stderr, "No destination for bezier_path_size in fit_bezier_curve_to_line"); + return; + } + + npoints = lines_size / 2; + + printf ("npoints: %d alloc: %d\n", npoints, npoints * sizeof (Point2)); + + points = malloc (npoints * sizeof (Point2)); + + j = 0; + for (i = 0; i < npoints; i++) { + points[i].x = lines[j]; + points[i].y = lines[j + 1]; + j += 2; + } + + printf ("simplified_path\n"); + simplified_path = malloc (8 * npoints * sizeof (double)); + simplified_path_buffer_size = 8 * npoints; + simplified_path_size = 0; + + FitCurve(points, npoints, error); + + *bezier_path = simplified_path; + *bezier_path_size = simplified_path_size; + + simplified_path = NULL; + + free (points); + } + + DrawBezierCurve(int n, Point2* curve) + { + int i; + + if (simplified_path_size + 8 > simplified_path_buffer_size) { + fprintf (stderr, "The bezier buffer is full (%d).\n", simplified_path_buffer_size); + return; + } + + if (n != 3) { + fprintf (stderr, "Expecting three points\n"); + return; + } + + i = simplified_path_size; + + simplified_path[i + 0] = curve[0].x; + simplified_path[i + 1] = curve[0].y; + simplified_path[i + 2] = curve[1].x; + simplified_path[i + 3] = curve[1].y; + simplified_path[i + 4] = curve[2].x; + simplified_path[i + 5] = curve[2].y; + simplified_path[i + 6] = curve[3].x; + simplified_path[i + 7] = curve[3].y; + + simplified_path_size = i + 8; + } + + /* + * 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); + } +

--- /dev/null +++ b/libbirdgems/fit_cubic.c~ @@ -1,1 +1,602 @@ + /* + 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 <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 */ + + int simplified_path_buffer_size = 0; + int simplified_path_size = 0; + double* simplified_path = NULL; + + /** Generates an a bezier path and returns the length of the output array. */ + int fit_bezier_curve_to_line (lines, lines_size, error, bezier_path) + double* lines, + int lines_size, + double error, + double** bezier_path /** generated Bezier curves */ + { + int i, j, n; + Point2* points + int npoints; + + npoints = lines_size / 2; + bezier_path = NULL; + + if (lines_size % 2 != 0) { + printf ("Odd number of coordinates."); + return 0; + } + + if (lines == NULL || lines_size == 0) { + printf ("No lines."); + return 0; + } + + if (simplified_path != NULL) { + printf ("Path simplification is alreading running."); + return 0; + } + + points = malloc (npoints * sizeof (Point2)); + + j = 0; + for (i = 0; i <= npoints; i++) { + points[i].x = lines[j]; + points[i].y = lines[j + 1]; + j += 2; + } + + simplified_path = malloc (4 * npoints * sizeof (double)); + simplified_path_buffer_size = 4 * npoints; + + FitCurve(points, n, error); + free (points); + + bezier_path = simplified_path; + simplified_path = NULL; + + return simplified_path_size; + } + + DrawBezierCurve(int n, Point2* curve) + { + int i; + + if (simplified_path_size + 4 > simplified_path_buffer_size) { + printf ("Out of memmory\n"); + return; + } + + if (n != 3) { + printf ("Expecting three points\n"); + } + + i = simplified_path_size; + + simplified_path[i++] = curve[0].x; + simplified_path[i++] = curve[0].y; + simplified_path[i++] = curve[1].x; + simplified_path[i++] = curve[1].y; + simplified_path[i++] = curve[2].x; + simplified_path[i++] = curve[2].y; + simplified_path[i++] = curve[3].x; + simplified_path[i++] = curve[3].y; + + simplified_path_size = i; + } + + /* + * 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); + } +

--- a/scripts/bavala.py +++ b/scripts/bavala.py @@ -82,13 +82,15 @@ if library: self.header = join(build, library) + '.h' - self.vapi = join(build, library) + '.vapi' + self.vapi = join(build, library) + '.vapi' # generated vapi file + self.other_vapi_files = get_sources_path (src, '*.vapi') # other vapi files self.so = join(build, src) + '.so.' + so_version self.so_link = join(build, src) + '.so' self.so_link_name = src + '.so' self.so_version = so_version self.so_name = 'lib' + library + '.so.' + so_version - + else: + self.other_vapi_files = [] def gen_c(self, opts): """translate code from vala to C and create .vapi""" @@ -107,25 +109,37 @@ dep_vapi = [d.vapi for d in self.vala_deps] action = cmd('valac', options, params, dep_vapi, self.vala) targets = self.cc[:] + if self.library: targets += [self.header, self.vapi] - + for f in self.c: yield { 'name': 'copy_c_' + f, 'actions': [ 'mkdir -p '+ self.build + '/' + self.src + '/', 'cp ' + f + ' ' + self.build + '/' + self.src + '/' + ], + } + + for f in self.other_vapi_files: + yield { + 'name': 'vapi_files_' + f, + 'actions': [ + 'mkdir -p '+ self.build + '/', + 'cp ' + f + ' ' + self.build + '/' ], } - print (action) - yield { - 'name': 'compile_c', - 'actions': [ action ], - 'file_dep': self.vala + dep_vapi, - 'targets': targets, - } + print (action) + + if not self.vala == []: + yield { + 'name': 'compile_c', + 'actions': [ action ], + 'file_dep': self.vala + dep_vapi, + 'targets': targets, + } def gen_o(self, opts):

--- a/scripts/build.py +++ b/scripts/build.py @@ -2,5 +2,5 @@ """ - Copyright (C) 2013, 2014 Johan Mattsson + Copyright (C) 2013, 2014 2015 Johan Mattsson This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as @@ -164,6 +164,61 @@ elif library.find ('.dylib') > -1: run ("""cd build/bin && ln -sf """ + library + " libbirdxml.dylib") + + def libbirdgems(prefix, cc, cflags, ldflags, valac, valaflags, library, nonNull = True): + #libbirdfont + run("mkdir -p build/libbirdgems") + run("mkdir -p build/bin") + + experimentalNonNull = "" + if nonNull: + experimentalNonNull = "--enable-experimental-non-null" + + run(valac + """\ + -C \ + """ + valaflags + """ \ + --pkg posix \ + --vapidir=./ \ + --basedir build/libbirdgems/ \ + """ + experimentalNonNull + """ \ + --enable-experimental \ + --library libbirdgems \ + -H build/libbirdxml/birdxml.h \ + libbirdxml/*.vala \ + """) + + if cc == "": + print ("Skipping compilation"); + else: + run(cc + " " + cflags + """ \ + -c build/libbirdgems/*.c \ + """) + + run("mv ./*.o build/libbirdgems/ ") + + if library.endswith (".dylib"): + sonameparam = "" # gcc on mac os does not have the soname parameter + else: + sonameparam = "-Wl,-soname," + library + + run(cc + " " + ldflags + """ \ + -shared \ + """ + sonameparam + """ \ + build/libbirdgems/*.o \ + $(pkg-config --libs glib-2.0) \ + $(pkg-config --libs gobject-2.0) \ + -o """ + library) + run("mv " + library + " build/bin/") + + if os.path.exists("build/bin/libbirdgems.so"): + run ("cd build/bin && unlink libbirdgems.so") + + # create link to the versioned library + if library.find ('.so') > -1: + run ("""cd build/bin && ln -sf """ + library + " libbirdgems.so") + elif library.find ('.dylib') > -1: + run ("""cd build/bin && ln -sf """ + library + " libbirdgems.dylib") + def birdfont_export(prefix, cc, cflags, ldflags, valac, valaflags, nonNull = True): # birdfont-export

--- a/scripts/linux_build.py +++ b/scripts/linux_build.py @@ -2,5 +2,5 @@ """ - Copyright (C) 2013 2014 Johan Mattsson + Copyright (C) 2013 2014 2015 Johan Mattsson This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as @@ -52,12 +52,16 @@ ldflags = options.ldflags + " " + "$(pkg-config --libs gdk-pixbuf-2.0)" library_cflags = options.cflags library_ldflags= options.ldflags + " -Wl,-soname," + "libbirdfont.so." + version.LIBBIRDXML_SO_VERSION + + birdgemslibrary_cflags = options.cflags + birdgemslibrary_ldflags= options.ldflags + " -Wl,-soname," + "birdgems.so." + version.LIBBIRDGEMS_SO_VERSION xmllibrary_cflags = options.cflags xmllibrary_ldflags= options.ldflags + " -Wl,-soname," + "libbirdxml.so." + version.SO_VERSION configfile.write_config (prefix) compile_translations() + build.libbirdgems(prefix, cc, birdgemslibrary_cflags, birdgemslibrary_ldflags, valac, valaflags, "birdgems.so." + version.LIBBIRDGEMS_SO_VERSION, False) build.libbirdxml(prefix, cc, xmllibrary_cflags, xmllibrary_ldflags, valac, valaflags, "libbirdxml.so." + version.LIBBIRDXML_SO_VERSION, False) build.libbirdfont(prefix, cc, library_cflags, library_ldflags, valac, valaflags, "libbirdfont.so." + version.SO_VERSION, False) build.birdfont_autotrace(prefix, cc, cflags, ldflags, valac, valaflags, False)

--- a/scripts/version.py +++ b/scripts/version.py @@ -21,4 +21,8 @@ LIBBIRDXML_SO_VERSION_MAJOR = '0' LIBBIRDXML_SO_VERSION_MINOR = '0' LIBBIRDXML_SO_VERSION = LIBBIRDXML_SO_VERSION_MAJOR + '.' + LIBBIRDXML_SO_VERSION_MINOR; + + LIBBIRDGEMS_SO_VERSION_MAJOR = '0' + LIBBIRDGEMS_SO_VERSION_MINOR = '0' + LIBBIRDGEMS_SO_VERSION = LIBBIRDGEMS_SO_VERSION_MAJOR + '.' + LIBBIRDGEMS_SO_VERSION_MINOR;

--- a/scripts/windows_build.py +++ b/scripts/windows_build.py @@ -30,6 +30,7 @@ from run import run compile_translations() + build.libbirdgems(prefix, cc, cflags, library_ldflags, valac, valaflags, "libbirdgems.dll") build.libbirdxml(prefix, cc, cflags, library_ldflags, valac, valaflags, "libbirdxml.dll") build.libbirdfont(prefix, cc, cflags, library_ldflags, valac, valaflags, "libbirdfont.dll")