The Birdfont Source Code


All Repositories / birdfont.git / blob – RSS feed

fit_cubic.c in libbirdfont/GraphicsGems

This file is a part of the Birdfont project.

Contributing

Send patches or pull requests to johan.mattsson.m@gmail.com.
Clone this repository: git clone https://github.com/johanmattssonm/birdfont.git

Revisions

View the latest version of libbirdfont/GraphicsGems/fit_cubic.c.
Create a library of the Graphics Gems code
1 /* 2 An Algorithm for Automatically Fitting Digitized Curves 3 by Philip J. Schneider 4 from "Graphics Gems", Academic Press, 1990 5 Adapted to BirdFont by Johan Mattsson 2015 6 */ 7 8 /* 9 * EULA: The Graphics Gems code is copyright-protected. In other words, 10 * you cannot claim the text of the code as your own and resell it. 11 * Using the code is permitted in any program, product, or library, 12 * non-commercial or commercial. Giving credit is not required, though 13 * is a nice gesture. The code comes as-is, and if there are any flaws 14 * or problems with any Gems code, nobody involved with Gems - authors, 15 * editors, publishers, or webmasters - are to be held responsible. 16 * Basically, don't be a jerk, and remember that anything free comes 17 * with no guarantee. 18 */ 19 20 /* fit_cubic.c */ 21 /* Piecewise cubic fitting code */ 22 23 #include "GraphicsGems.h" 24 #include "birdfont.h" 25 #include <stdio.h> 26 #include <malloc.h> 27 #include <math.h> 28 #include <stdio.h> 29 30 typedef Point2 *BezierCurve; 31 32 /* Forward declarations */ 33 void FitCurve(); 34 static void FitCubic(); 35 static double *Reparameterize(); 36 static double NewtonRaphsonRootFind(); 37 static Point2 BezierII(); 38 static double B0(), B1(), B2(), B3(); 39 static Vector2 ComputeLeftTangent(); 40 static Vector2 ComputeRightTangent(); 41 static Vector2 ComputeCenterTangent(); 42 static double ComputeMaxError(); 43 static double *ChordLengthParameterize(); 44 static BezierCurve GenerateBezier(); 45 static Vector2 V2AddII(); 46 static Vector2 V2ScaleIII(); 47 static Vector2 V2SubII(); 48 49 #define MAXPOINTS 1000 /* The most points you can have */ 50 51 BirdFontPath* simplified_path = NULL; 52 53 BirdFontPath* fit_bezier_curve_to_line (BirdFontPath* path, int start, int stop, double error) { 54 int i, j, n; 55 Point2* points; 56 BirdFontPath* result; 57 58 if (simplified_path != NULL) { 59 printf ("Simplifier is not thread safe."); 60 return bird_font_path_new (); 61 } 62 63 n = stop - start + 1; 64 65 if (n > bird_font_path_size (path) || n < 0) { 66 printf ("Index of of bounds in simplifier"); 67 return bird_font_path_new (); 68 } 69 70 points = malloc (n * sizeof (Point2)); 71 72 j = 0; 73 for (i = start; i <= stop; i++) { 74 points[j].x = bird_font_path_get_x_at (path, i); 75 points[j].y = bird_font_path_get_y_at (path, i); 76 j++; 77 } 78 79 simplified_path = bird_font_path_new (); 80 FitCurve(points, n, error); 81 free (points); 82 83 result = simplified_path; 84 simplified_path = NULL; 85 return result; 86 } 87 88 DrawBezierCurve(int n, Point2* curve) { 89 if (n == 3) { 90 bird_font_path_add_cubic_bezier_points (simplified_path, 91 curve[0].x, curve[0].y, 92 curve[1].x, curve[1].y, 93 curve[2].x, curve[2].y, 94 curve[3].x, curve[3].y); 95 } else { 96 printf ("Expecting three points\n"); 97 } 98 } 99 100 /* 101 * FitCurve : 102 * Fit a Bezier curve to a set of digitized points 103 */ 104 void FitCurve(d, nPts, error) 105 Point2 *d; /* Array of digitized points */ 106 int nPts; /* Number of digitized points */ 107 double error; /* User-defined error squared */ 108 { 109 Vector2 tHat1, tHat2; /* Unit tangent vectors at endpoints */ 110 111 tHat1 = ComputeLeftTangent(d, 0); 112 tHat2 = ComputeRightTangent(d, nPts - 1); 113 FitCubic(d, 0, nPts - 1, tHat1, tHat2, error); 114 } 115 116 /* 117 * FitCubic : 118 * Fit a Bezier curve to a (sub)set of digitized points 119 */ 120 static void FitCubic(d, first, last, tHat1, tHat2, error) 121 Point2 *d; /* Array of digitized points */ 122 int first, last; /* Indices of first and last pts in region */ 123 Vector2 tHat1, tHat2; /* Unit tangent vectors at endpoints */ 124 double error; /* User-defined error squared */ 125 { 126 BezierCurve bezCurve; /*Control points of fitted Bezier curve*/ 127 double *u; /* Parameter values for point */ 128 double *uPrime; /* Improved parameter values */ 129 double maxError; /* Maximum fitting error */ 130 int splitPoint; /* Point to split point set at */ 131 int nPts; /* Number of points in subset */ 132 double iterationError; /*Error below which you try iterating */ 133 int maxIterations = 4; /* Max times to try iterating */ 134 Vector2 tHatCenter; /* Unit tangent vector at splitPoint */ 135 int i; 136 137 iterationError = error * error; 138 nPts = last - first + 1; 139 140 /* Use heuristic if region only has two points in it */ 141 if (nPts == 2) { 142 double dist = V2DistanceBetween2Points(&d[last], &d[first]) / 3.0; 143 144 bezCurve = (Point2 *)malloc(4 * sizeof(Point2)); 145 bezCurve[0] = d[first]; 146 bezCurve[3] = d[last]; 147 V2Add(&bezCurve[0], V2Scale(&tHat1, dist), &bezCurve[1]); 148 V2Add(&bezCurve[3], V2Scale(&tHat2, dist), &bezCurve[2]); 149 DrawBezierCurve(3, bezCurve); 150 free((void *)bezCurve); 151 return; 152 } 153 154 /* Parameterize points, and attempt to fit curve */ 155 u = ChordLengthParameterize(d, first, last); 156 bezCurve = GenerateBezier(d, first, last, u, tHat1, tHat2); 157 158 /* Find max deviation of points to fitted curve */ 159 maxError = ComputeMaxError(d, first, last, bezCurve, u, &splitPoint); 160 if (maxError < error) { 161 DrawBezierCurve(3, bezCurve); 162 free((void *)u); 163 free((void *)bezCurve); 164 return; 165 } 166 167 168 /* If error not too large, try some reparameterization */ 169 /* and iteration */ 170 if (maxError < iterationError) { 171 for (i = 0; i < maxIterations; i++) { 172 uPrime = Reparameterize(d, first, last, u, bezCurve); 173 free((void *)bezCurve); 174 bezCurve = GenerateBezier(d, first, last, uPrime, tHat1, tHat2); 175 maxError = ComputeMaxError(d, first, last, 176 bezCurve, uPrime, &splitPoint); 177 if (maxError < error) { 178 DrawBezierCurve(3, bezCurve); 179 free((void *)u); 180 free((void *)bezCurve); 181 free((void *)uPrime); 182 return; 183 } 184 free((void *)u); 185 u = uPrime; 186 } 187 } 188 189 /* Fitting failed -- split at max error point and fit recursively */ 190 free((void *)u); 191 free((void *)bezCurve); 192 tHatCenter = ComputeCenterTangent(d, splitPoint); 193 FitCubic(d, first, splitPoint, tHat1, tHatCenter, error); 194 V2Negate(&tHatCenter); 195 FitCubic(d, splitPoint, last, tHatCenter, tHat2, error); 196 } 197 198 199 /* 200 * GenerateBezier : 201 * Use least-squares method to find Bezier control points for region. 202 * 203 */ 204 static BezierCurve GenerateBezier(d, first, last, uPrime, tHat1, tHat2) 205 Point2 *d; /* Array of digitized points */ 206 int first, last; /* Indices defining region */ 207 double *uPrime; /* Parameter values for region */ 208 Vector2 tHat1, tHat2; /* Unit tangents at endpoints */ 209 { 210 int i; 211 Vector2 A[MAXPOINTS][2]; /* Precomputed rhs for eqn */ 212 int nPts; /* Number of pts in sub-curve */ 213 double C[2][2]; /* Matrix C */ 214 double X[2]; /* Matrix X */ 215 double det_C0_C1, /* Determinants of matrices */ 216 det_C0_X, 217 det_X_C1; 218 double alpha_l, /* Alpha values, left and right */ 219 alpha_r; 220 Vector2 tmp; /* Utility variable */ 221 BezierCurve bezCurve; /* RETURN bezier curve ctl pts */ 222 223 bezCurve = (Point2 *)malloc(4 * sizeof(Point2)); 224 nPts = last - first + 1; 225 226 227 /* Compute the A's */ 228 for (i = 0; i < nPts; i++) { 229 Vector2 v1, v2; 230 v1 = tHat1; 231 v2 = tHat2; 232 V2Scale(&v1, B1(uPrime[i])); 233 V2Scale(&v2, B2(uPrime[i])); 234 A[i][0] = v1; 235 A[i][1] = v2; 236 } 237 238 /* Create the C and X matrices */ 239 C[0][0] = 0.0; 240 C[0][1] = 0.0; 241 C[1][0] = 0.0; 242 C[1][1] = 0.0; 243 X[0] = 0.0; 244 X[1] = 0.0; 245 246 for (i = 0; i < nPts; i++) { 247 C[0][0] += V2Dot(&A[i][0], &A[i][0]); 248 C[0][1] += V2Dot(&A[i][0], &A[i][1]); 249 /* C[1][0] += V2Dot(&A[i][0], &A[i][1]);*/ 250 C[1][0] = C[0][1]; 251 C[1][1] += V2Dot(&A[i][1], &A[i][1]); 252 253 tmp = V2SubII(d[first + i], 254 V2AddII( 255 V2ScaleIII(d[first], B0(uPrime[i])), 256 V2AddII( 257 V2ScaleIII(d[first], B1(uPrime[i])), 258 V2AddII( 259 V2ScaleIII(d[last], B2(uPrime[i])), 260 V2ScaleIII(d[last], B3(uPrime[i])))))); 261 262 263 X[0] += V2Dot(&A[i][0], &tmp); 264 X[1] += V2Dot(&A[i][1], &tmp); 265 } 266 267 /* Compute the determinants of C and X */ 268 det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1]; 269 det_C0_X = C[0][0] * X[1] - C[1][0] * X[0]; 270 det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1]; 271 272 /* Finally, derive alpha values */ 273 alpha_l = (det_C0_C1 == 0) ? 0.0 : det_X_C1 / det_C0_C1; 274 alpha_r = (det_C0_C1 == 0) ? 0.0 : det_C0_X / det_C0_C1; 275 276 /* If alpha negative, use the Wu/Barsky heuristic (see text) */ 277 /* (if alpha is 0, you get coincident control points that lead to 278 * divide by zero in any subsequent NewtonRaphsonRootFind() call. */ 279 double segLength = V2DistanceBetween2Points(&d[last], &d[first]); 280 double epsilon = 1.0e-6 * segLength; 281 if (alpha_l < epsilon || alpha_r < epsilon) 282 { 283 /* fall back on standard (probably inaccurate) formula, and subdivide further if needed. */ 284 double dist = segLength / 3.0; 285 bezCurve[0] = d[first]; 286 bezCurve[3] = d[last]; 287 V2Add(&bezCurve[0], V2Scale(&tHat1, dist), &bezCurve[1]); 288 V2Add(&bezCurve[3], V2Scale(&tHat2, dist), &bezCurve[2]); 289 return (bezCurve); 290 } 291 292 /* First and last control points of the Bezier curve are */ 293 /* positioned exactly at the first and last data points */ 294 /* Control points 1 and 2 are positioned an alpha distance out */ 295 /* on the tangent vectors, left and right, respectively */ 296 bezCurve[0] = d[first]; 297 bezCurve[3] = d[last]; 298 V2Add(&bezCurve[0], V2Scale(&tHat1, alpha_l), &bezCurve[1]); 299 V2Add(&bezCurve[3], V2Scale(&tHat2, alpha_r), &bezCurve[2]); 300 return (bezCurve); 301 } 302 303 304 /* 305 * Reparameterize: 306 * Given set of points and their parameterization, try to find 307 * a better parameterization. 308 * 309 */ 310 static double *Reparameterize(d, first, last, u, bezCurve) 311 Point2 *d; /* Array of digitized points */ 312 int first, last; /* Indices defining region */ 313 double *u; /* Current parameter values */ 314 BezierCurve bezCurve; /* Current fitted curve */ 315 { 316 int nPts = last-first+1; 317 int i; 318 double *uPrime; /* New parameter values */ 319 320 uPrime = (double *)malloc(nPts * sizeof(double)); 321 for (i = first; i <= last; i++) { 322 uPrime[i-first] = NewtonRaphsonRootFind(bezCurve, d[i], u[i- 323 first]); 324 } 325 return (uPrime); 326 } 327 328 329 330 /* 331 * NewtonRaphsonRootFind : 332 * Use Newton-Raphson iteration to find better root. 333 */ 334 static double NewtonRaphsonRootFind(Q, P, u) 335 BezierCurve Q; /* Current fitted curve */ 336 Point2 P; /* Digitized point */ 337 double u; /* Parameter value for "P" */ 338 { 339 double numerator, denominator; 340 Point2 Q1[3], Q2[2]; /* Q' and Q'' */ 341 Point2 Q_u, Q1_u, Q2_u; /*u evaluated at Q, Q', & Q'' */ 342 double uPrime; /* Improved u */ 343 int i; 344 345 /* Compute Q(u) */ 346 Q_u = BezierII(3, Q, u); 347 348 /* Generate control vertices for Q' */ 349 for (i = 0; i <= 2; i++) { 350 Q1[i].x = (Q[i+1].x - Q[i].x) * 3.0; 351 Q1[i].y = (Q[i+1].y - Q[i].y) * 3.0; 352 } 353 354 /* Generate control vertices for Q'' */ 355 for (i = 0; i <= 1; i++) { 356 Q2[i].x = (Q1[i+1].x - Q1[i].x) * 2.0; 357 Q2[i].y = (Q1[i+1].y - Q1[i].y) * 2.0; 358 } 359 360 /* Compute Q'(u) and Q''(u) */ 361 Q1_u = BezierII(2, Q1, u); 362 Q2_u = BezierII(1, Q2, u); 363 364 /* Compute f(u)/f'(u) */ 365 numerator = (Q_u.x - P.x) * (Q1_u.x) + (Q_u.y - P.y) * (Q1_u.y); 366 denominator = (Q1_u.x) * (Q1_u.x) + (Q1_u.y) * (Q1_u.y) + 367 (Q_u.x - P.x) * (Q2_u.x) + (Q_u.y - P.y) * (Q2_u.y); 368 if (denominator == 0.0f) return u; 369 370 /* u = u - f(u)/f'(u) */ 371 uPrime = u - (numerator/denominator); 372 return (uPrime); 373 } 374 375 376 377 /* 378 * Bezier : 379 * Evaluate a Bezier curve at a particular parameter value 380 * 381 */ 382 static Point2 BezierII(degree, V, t) 383 int degree; /* The degree of the bezier curve */ 384 Point2 *V; /* Array of control points */ 385 double t; /* Parametric value to find point for */ 386 { 387 int i, j; 388 Point2 Q; /* Point on curve at parameter t */ 389 Point2 *Vtemp; /* Local copy of control points */ 390 391 /* Copy array */ 392 Vtemp = (Point2 *)malloc((unsigned)((degree+1) 393 * sizeof (Point2))); 394 for (i = 0; i <= degree; i++) { 395 Vtemp[i] = V[i]; 396 } 397 398 /* Triangle computation */ 399 for (i = 1; i <= degree; i++) { 400 for (j = 0; j <= degree-i; j++) { 401 Vtemp[j].x = (1.0 - t) * Vtemp[j].x + t * Vtemp[j+1].x; 402 Vtemp[j].y = (1.0 - t) * Vtemp[j].y + t * Vtemp[j+1].y; 403 } 404 } 405 406 Q = Vtemp[0]; 407 free((void *)Vtemp); 408 return Q; 409 } 410 411 412 /* 413 * B0, B1, B2, B3 : 414 * Bezier multipliers 415 */ 416 static double B0(u) 417 double u; 418 { 419 double tmp = 1.0 - u; 420 return (tmp * tmp * tmp); 421 } 422 423 424 static double B1(u) 425 double u; 426 { 427 double tmp = 1.0 - u; 428 return (3 * u * (tmp * tmp)); 429 } 430 431 static double B2(u) 432 double u; 433 { 434 double tmp = 1.0 - u; 435 return (3 * u * u * tmp); 436 } 437 438 static double B3(u) 439 double u; 440 { 441 return (u * u * u); 442 } 443 444 445 446 /* 447 * ComputeLeftTangent, ComputeRightTangent, ComputeCenterTangent : 448 *Approximate unit tangents at endpoints and "center" of digitized curve 449 */ 450 static Vector2 ComputeLeftTangent(d, end) 451 Point2 *d; /* Digitized points*/ 452 int end; /* Index to "left" end of region */ 453 { 454 Vector2 tHat1; 455 tHat1 = V2SubII(d[end+1], d[end]); 456 tHat1 = *V2Normalize(&tHat1); 457 return tHat1; 458 } 459 460 static Vector2 ComputeRightTangent(d, end) 461 Point2 *d; /* Digitized points */ 462 int end; /* Index to "right" end of region */ 463 { 464 Vector2 tHat2; 465 tHat2 = V2SubII(d[end-1], d[end]); 466 tHat2 = *V2Normalize(&tHat2); 467 return tHat2; 468 } 469 470 471 static Vector2 ComputeCenterTangent(d, center) 472 Point2 *d; /* Digitized points */ 473 int center; /* Index to point inside region */ 474 { 475 Vector2 V1, V2, tHatCenter; 476 477 V1 = V2SubII(d[center-1], d[center]); 478 V2 = V2SubII(d[center], d[center+1]); 479 tHatCenter.x = (V1.x + V2.x)/2.0; 480 tHatCenter.y = (V1.y + V2.y)/2.0; 481 tHatCenter = *V2Normalize(&tHatCenter); 482 return tHatCenter; 483 } 484 485 486 /* 487 * ChordLengthParameterize : 488 * Assign parameter values to digitized points 489 * using relative distances between points. 490 */ 491 static double *ChordLengthParameterize(d, first, last) 492 Point2 *d; /* Array of digitized points */ 493 int first, last; /* Indices defining region */ 494 { 495 int i; 496 double *u; /* Parameterization */ 497 498 u = (double *)malloc((unsigned)(last-first+1) * sizeof(double)); 499 500 u[0] = 0.0; 501 for (i = first+1; i <= last; i++) { 502 u[i-first] = u[i-first-1] + 503 V2DistanceBetween2Points(&d[i], &d[i-1]); 504 } 505 506 for (i = first + 1; i <= last; i++) { 507 u[i-first] = u[i-first] / u[last-first]; 508 } 509 510 return(u); 511 } 512 513 514 515 516 /* 517 * ComputeMaxError : 518 * Find the maximum squared distance of digitized points 519 * to fitted curve. 520 */ 521 static double ComputeMaxError(d, first, last, bezCurve, u, splitPoint) 522 Point2 *d; /* Array of digitized points */ 523 int first, last; /* Indices defining region */ 524 BezierCurve bezCurve; /* Fitted Bezier curve */ 525 double *u; /* Parameterization of points */ 526 int *splitPoint; /* Point of maximum error */ 527 { 528 int i; 529 double maxDist; /* Maximum error */ 530 double dist; /* Current error */ 531 Point2 P; /* Point on curve */ 532 Vector2 v; /* Vector from point to curve */ 533 534 *splitPoint = (last - first + 1)/2; 535 maxDist = 0.0; 536 for (i = first + 1; i < last; i++) { 537 P = BezierII(3, bezCurve, u[i-first]); 538 v = V2SubII(P, d[i]); 539 dist = V2SquaredLength(&v); 540 if (dist >= maxDist) { 541 maxDist = dist; 542 *splitPoint = i; 543 } 544 } 545 return (maxDist); 546 } 547 static Vector2 V2AddII(a, b) 548 Vector2 a, b; 549 { 550 Vector2 c; 551 c.x = a.x + b.x; c.y = a.y + b.y; 552 return (c); 553 } 554 static Vector2 V2ScaleIII(v, s) 555 Vector2 v; 556 double s; 557 { 558 Vector2 result; 559 result.x = v.x * s; result.y = v.y * s; 560 return (result); 561 } 562 563 static Vector2 V2SubII(a, b) 564 Vector2 a, b; 565 { 566 Vector2 c; 567 c.x = a.x - b.x; c.y = a.y - b.y; 568 return (c); 569 } 570 571