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