The Birdfont Source Code


All Repositories / birdfont.git / blob – RSS feed

fit_cubic.c in libbirdgems

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