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 MAC
26 #include <malloc/malloc.h>
27 #else
28 #include <malloc.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