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