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