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