1
#define __SP_BEZIER_UTILS_C__
4
* Bezier interpolation for inkscape drawing code.
7
* Original code published in:
8
* An Algorithm for Automatically Fitting Digitized Curves
9
* by Philip J. Schneider
10
* "Graphics Gems", Academic Press, 1990
14
* Lauris Kaplinski <lauris@kaplinski.com>
15
* Peter Moulder <pmoulder@mail.csse.monash.edu.au>
17
* Copyright (C) 1990 Philip J. Schneider
18
* Copyright (C) 2001 Lauris Kaplinski
19
* Copyright (C) 2001 Ximian, Inc.
20
* Copyright (C) 2003,2004 Monash University
22
* This library is free software; you can redistribute it and/or
23
* modify it either under the terms of the GNU Lesser General Public
24
* License version 2.1 as published by the Free Software Foundation
25
* (the "LGPL") or, at your option, under the terms of the Mozilla
26
* Public License Version 1.1 (the "MPL"). If you do not alter this
27
* notice, a recipient may use your version of this file under either
28
* the MPL or the LGPL.
30
* You should have received a copy of the LGPL along with this library
31
* in the file COPYING-LGPL-2.1; if not, write to the Free Software
32
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
33
* You should have received a copy of the MPL along with this library
34
* in the file COPYING-MPL-1.1
36
* The contents of this file are subject to the Mozilla Public License
37
* Version 1.1 (the "License"); you may not use this file except in
38
* compliance with the License. You may obtain a copy of the License at
39
* http://www.mozilla.org/MPL/
41
* This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
42
* OF ANY KIND, either express or implied. See the LGPL or the MPL for
43
* the specific language governing rights and limitations.
48
#define noBEZIER_DEBUG
54
#include <2geom/bezier-utils.h>
56
#include <2geom/isnan.h>
61
typedef Point BezierCurve[];
63
/* Forward declarations */
64
static void generate_bezier(Point b[], Point const d[], double const u[], unsigned len,
65
Point const &tHat1, Point const &tHat2, double tolerance_sq);
66
static void estimate_lengths(Point bezier[],
67
Point const data[], double const u[], unsigned len,
68
Point const &tHat1, Point const &tHat2);
69
static void estimate_bi(Point b[4], unsigned ei,
70
Point const data[], double const u[], unsigned len);
71
static void reparameterize(Point const d[], unsigned len, double u[], BezierCurve const bezCurve);
72
static double NewtonRaphsonRootFind(BezierCurve const Q, Point const &P, double u);
73
static Point darray_center_tangent(Point const d[], unsigned center, unsigned length);
74
static Point darray_right_tangent(Point const d[], unsigned const len);
75
static unsigned copy_without_nans_or_adjacent_duplicates(Point const src[], unsigned src_len, Point dest[]);
76
static void chord_length_parameterize(Point const d[], double u[], unsigned len);
77
static double compute_max_error_ratio(Point const d[], double const u[], unsigned len,
78
BezierCurve const bezCurve, double tolerance,
79
unsigned *splitPoint);
80
static double compute_hook(Point const &a, Point const &b, double const u, BezierCurve const bezCurve,
81
double const tolerance);
84
static Point const unconstrained_tangent(0, 0);
88
* B0, B1, B2, B3 : Bezier multipliers
91
#define B0(u) ( ( 1.0 - u ) * ( 1.0 - u ) * ( 1.0 - u ) )
92
#define B1(u) ( 3 * u * ( 1.0 - u ) * ( 1.0 - u ) )
93
#define B2(u) ( 3 * u * u * ( 1.0 - u ) )
94
#define B3(u) ( u * u * u )
97
# define DOUBLE_ASSERT(x) assert( ( (x) > -SP_HUGE ) && ( (x) < SP_HUGE ) )
98
# define BEZIER_ASSERT(b) do { \
99
DOUBLE_ASSERT((b)[0][X]); DOUBLE_ASSERT((b)[0][Y]); \
100
DOUBLE_ASSERT((b)[1][X]); DOUBLE_ASSERT((b)[1][Y]); \
101
DOUBLE_ASSERT((b)[2][X]); DOUBLE_ASSERT((b)[2][Y]); \
102
DOUBLE_ASSERT((b)[3][X]); DOUBLE_ASSERT((b)[3][Y]); \
105
# define DOUBLE_ASSERT(x) do { } while(0)
106
# define BEZIER_ASSERT(b) do { } while(0)
111
* Fit a single-segment Bezier curve to a set of digitized points.
113
* \return Number of segments generated, or -1 on error.
116
bezier_fit_cubic(Point *bezier, Point const *data, int len, double error)
118
return bezier_fit_cubic_r(bezier, data, len, error, 1);
122
* Fit a multi-segment Bezier curve to a set of digitized points, with
123
* possible weedout of identical points and NaNs.
125
* \param max_beziers Maximum number of generated segments
126
* \param Result array, must be large enough for n. segments * 4 elements.
128
* \return Number of segments generated, or -1 on error.
131
bezier_fit_cubic_r(Point bezier[], Point const data[], int const len, double const error, unsigned const max_beziers)
136
max_beziers >= (1ul << (31 - 2 - 1 - 3)))
139
Point *uniqued_data = new Point[len];
140
unsigned uniqued_len = copy_without_nans_or_adjacent_duplicates(data, len, uniqued_data);
142
if ( uniqued_len < 2 ) {
143
delete[] uniqued_data;
147
/* Call fit-cubic function with recursion. */
148
int const ret = bezier_fit_cubic_full(bezier, NULL, uniqued_data, uniqued_len,
149
unconstrained_tangent, unconstrained_tangent,
151
delete[] uniqued_data;
156
* Copy points from src to dest, filter out points containing NaN and
157
* adjacent points with equal x and y.
158
* \return length of dest
161
copy_without_nans_or_adjacent_duplicates(Point const src[], unsigned src_len, Point dest[])
165
if ( si == src_len ) {
168
if (!IS_NAN(src[si][X]) &&
169
!IS_NAN(src[si][Y])) {
170
dest[0] = Point(src[si]);
177
for (; si < src_len; ++si) {
178
Point const src_pt = Point(src[si]);
179
if ( src_pt != dest[di]
180
&& !IS_NAN(src_pt[X])
181
&& !IS_NAN(src_pt[Y])) {
185
unsigned dest_len = di + 1;
186
assert( dest_len <= src_len );
191
* Fit a multi-segment Bezier curve to a set of digitized points, without
192
* possible weedout of identical points and NaNs.
194
* \pre data is uniqued, i.e. not exist i: data[i] == data[i + 1].
195
* \param max_beziers Maximum number of generated segments
196
* \param Result array, must be large enough for n. segments * 4 elements.
199
bezier_fit_cubic_full(Point bezier[], int split_points[],
200
Point const data[], int const len,
201
Point const &tHat1, Point const &tHat2,
202
double const error, unsigned const max_beziers)
204
int const maxIterations = 4; /* std::max times to try iterating */
206
if(!(bezier != NULL) ||
209
!(max_beziers >= 1) ||
213
if ( len < 2 ) return 0;
216
/* We have 2 points, which can be fitted trivially. */
218
bezier[3] = data[len - 1];
219
double const dist = distance(bezier[0], bezier[3]) / 3.0;
221
/* Numerical problem, fall back to straight line segment. */
222
bezier[1] = bezier[0];
223
bezier[2] = bezier[3];
225
bezier[1] = ( is_zero(tHat1)
226
? ( 2 * bezier[0] + bezier[3] ) / 3.
227
: bezier[0] + dist * tHat1 );
228
bezier[2] = ( is_zero(tHat2)
229
? ( bezier[0] + 2 * bezier[3] ) / 3.
230
: bezier[3] + dist * tHat2 );
232
BEZIER_ASSERT(bezier);
236
/* Parameterize points, and attempt to fit curve */
237
unsigned splitPoint; /* Point to split point set at. */
240
double *u = new double[len];
241
chord_length_parameterize(data, u, len);
242
if ( u[len - 1] == 0.0 ) {
243
/* Zero-length path: every point in data[] is the same.
245
* (Clients aren't allowed to pass such data; handling the case is defensive
252
generate_bezier(bezier, data, u, len, tHat1, tHat2, error);
253
reparameterize(data, len, u, bezier);
255
/* Find max deviation of points to fitted curve. */
256
double const tolerance = sqrt(error + 1e-9);
257
double maxErrorRatio = compute_max_error_ratio(data, u, len, bezier, tolerance, &splitPoint);
259
if ( fabs(maxErrorRatio) <= 1.0 ) {
260
BEZIER_ASSERT(bezier);
265
/* If error not too large, then try some reparameterization and iteration. */
266
if ( 0.0 <= maxErrorRatio && maxErrorRatio <= 3.0 ) {
267
for (int i = 0; i < maxIterations; i++) {
268
generate_bezier(bezier, data, u, len, tHat1, tHat2, error);
269
reparameterize(data, len, u, bezier);
270
maxErrorRatio = compute_max_error_ratio(data, u, len, bezier, tolerance, &splitPoint);
271
if ( fabs(maxErrorRatio) <= 1.0 ) {
272
BEZIER_ASSERT(bezier);
279
is_corner = (maxErrorRatio < 0);
283
assert(splitPoint < unsigned(len));
284
if (splitPoint == 0) {
285
if (is_zero(tHat1)) {
286
/* Got spike even with unconstrained initial tangent. */
289
return bezier_fit_cubic_full(bezier, split_points, data, len, unconstrained_tangent, tHat2,
292
} else if (splitPoint == unsigned(len - 1)) {
293
if (is_zero(tHat2)) {
294
/* Got spike even with unconstrained final tangent. */
297
return bezier_fit_cubic_full(bezier, split_points, data, len, tHat1, unconstrained_tangent,
303
if ( 1 < max_beziers ) {
305
* Fitting failed -- split at max error point and fit recursively
307
unsigned const rec_max_beziers1 = max_beziers - 1;
309
Point recTHat2, recTHat1;
311
if(!(0 < splitPoint && splitPoint < unsigned(len - 1)))
313
recTHat1 = recTHat2 = unconstrained_tangent;
315
/* Unit tangent vector at splitPoint. */
316
recTHat2 = darray_center_tangent(data, splitPoint, len);
317
recTHat1 = -recTHat2;
319
int const nsegs1 = bezier_fit_cubic_full(bezier, split_points, data, splitPoint + 1,
320
tHat1, recTHat2, error, rec_max_beziers1);
323
g_print("fit_cubic[1]: recursive call failed\n");
327
assert( nsegs1 != 0 );
328
if (split_points != NULL) {
329
split_points[nsegs1 - 1] = splitPoint;
331
unsigned const rec_max_beziers2 = max_beziers - nsegs1;
332
int const nsegs2 = bezier_fit_cubic_full(bezier + nsegs1*4,
333
( split_points == NULL
335
: split_points + nsegs1 ),
336
data + splitPoint, len - splitPoint,
337
recTHat1, tHat2, error, rec_max_beziers2);
340
g_print("fit_cubic[2]: recursive call failed\n");
346
g_print("fit_cubic: success[nsegs: %d+%d=%d] on max_beziers:%u\n",
347
nsegs1, nsegs2, nsegs1 + nsegs2, max_beziers);
349
return nsegs1 + nsegs2;
357
* Fill in \a bezier[] based on the given data and tangent requirements, using
358
* a least-squares fit.
360
* Each of tHat1 and tHat2 should be either a zero vector or a unit vector.
361
* If it is zero, then bezier[1 or 2] is estimated without constraint; otherwise,
362
* it bezier[1 or 2] is placed in the specified direction from bezier[0 or 3].
364
* \param tolerance_sq Used only for an initial guess as to tangent directions
365
* when \a tHat1 or \a tHat2 is zero.
368
generate_bezier(Point bezier[],
369
Point const data[], double const u[], unsigned const len,
370
Point const &tHat1, Point const &tHat2,
371
double const tolerance_sq)
373
bool const est1 = is_zero(tHat1);
374
bool const est2 = is_zero(tHat2);
375
Point est_tHat1( est1
376
? darray_left_tangent(data, len, tolerance_sq)
378
Point est_tHat2( est2
379
? darray_right_tangent(data, len, tolerance_sq)
381
estimate_lengths(bezier, data, u, len, est_tHat1, est_tHat2);
382
/* We find that darray_right_tangent tends to produce better results
383
for our current freehand tool than full estimation. */
385
estimate_bi(bezier, 1, data, u, len);
386
if (bezier[1] != bezier[0]) {
387
est_tHat1 = unit_vector(bezier[1] - bezier[0]);
389
estimate_lengths(bezier, data, u, len, est_tHat1, est_tHat2);
395
estimate_lengths(Point bezier[],
396
Point const data[], double const uPrime[], unsigned const len,
397
Point const &tHat1, Point const &tHat2)
399
double C[2][2]; /* Matrix C. */
400
double X[2]; /* Matrix X. */
402
/* Create the C and X matrices. */
410
/* First and last control points of the Bezier curve are positioned exactly at the first and
413
bezier[3] = data[len - 1];
415
for (unsigned i = 0; i < len; i++) {
416
/* Bezier control point coefficients. */
417
double const b0 = B0(uPrime[i]);
418
double const b1 = B1(uPrime[i]);
419
double const b2 = B2(uPrime[i]);
420
double const b3 = B3(uPrime[i]);
423
Point const a1 = b1 * tHat1;
424
Point const a2 = b2 * tHat2;
426
C[0][0] += dot(a1, a1);
427
C[0][1] += dot(a1, a2);
429
C[1][1] += dot(a2, a2);
431
/* Additional offset to the data point from the predicted point if we were to set bezier[1]
432
to bezier[0] and bezier[2] to bezier[3]. */
433
Point const shortfall
435
- ( ( b0 + b1 ) * bezier[0] )
436
- ( ( b2 + b3 ) * bezier[3] ) );
437
X[0] += dot(a1, shortfall);
438
X[1] += dot(a2, shortfall);
441
/* We've constructed a pair of equations in the form of a matrix product C * alpha = X.
442
Now solve for alpha. */
443
double alpha_l, alpha_r;
445
/* Compute the determinants of C and X. */
446
double const det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
447
if ( det_C0_C1 != 0 ) {
448
/* Apparently Kramer's rule. */
449
double const det_C0_X = C[0][0] * X[1] - C[0][1] * X[0];
450
double const det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1];
451
alpha_l = det_X_C1 / det_C0_C1;
452
alpha_r = det_C0_X / det_C0_C1;
454
/* The matrix is under-determined. Try requiring alpha_l == alpha_r.
456
* One way of implementing the constraint alpha_l == alpha_r is to treat them as the same
457
* variable in the equations. We can do this by adding the columns of C to form a single
458
* column, to be multiplied by alpha to give the column vector X.
460
* We try each row in turn.
462
double const c0 = C[0][0] + C[0][1];
464
alpha_l = alpha_r = X[0] / c0;
466
double const c1 = C[1][0] + C[1][1];
468
alpha_l = alpha_r = X[1] / c1;
470
/* Let the below code handle this. */
471
alpha_l = alpha_r = 0.;
476
/* If alpha negative, use the Wu/Barsky heuristic (see text). (If alpha is 0, you get
477
coincident control points that lead to divide by zero in any subsequent
478
NewtonRaphsonRootFind() call.) */
479
/// \todo Check whether this special-casing is necessary now that
480
/// NewtonRaphsonRootFind handles non-positive denominator.
481
if ( alpha_l < 1.0e-6 ||
484
alpha_l = alpha_r = distance(data[0], data[len-1]) / 3.0;
487
/* Control points 1 and 2 are positioned an alpha distance out on the tangent vectors, left and
488
right, respectively. */
489
bezier[1] = alpha_l * tHat1 + bezier[0];
490
bezier[2] = alpha_r * tHat2 + bezier[3];
495
static double lensq(Point const p) {
500
estimate_bi(Point bezier[4], unsigned const ei,
501
Point const data[], double const u[], unsigned const len)
503
if(!(1 <= ei && ei <= 2))
505
unsigned const oi = 3 - ei;
506
double num[2] = {0., 0.};
508
for (unsigned i = 0; i < len; ++i) {
509
double const ui = u[i];
510
double const b[4] = {
517
for (unsigned d = 0; d < 2; ++d) {
518
num[d] += b[ei] * (b[0] * bezier[0][d] +
519
b[oi] * bezier[oi][d] +
520
b[3] * bezier[3][d] +
523
den -= b[ei] * b[ei];
527
for (unsigned d = 0; d < 2; ++d) {
528
bezier[ei][d] = num[d] / den;
531
bezier[ei] = ( oi * bezier[0] + ei * bezier[3] ) / 3.;
536
* Given set of points and their parameterization, try to find a better assignment of parameter
537
* values for the points.
539
* \param d Array of digitized points.
540
* \param u Current parameter values.
541
* \param bezCurve Current fitted curve.
542
* \param len Number of values in both d and u arrays.
543
* Also the size of the array that is allocated for return.
546
reparameterize(Point const d[],
549
BezierCurve const bezCurve)
553
unsigned const last = len - 1;
554
assert( bezCurve[0] == d[0] );
555
assert( bezCurve[3] == d[last] );
556
assert( u[0] == 0.0 );
557
assert( u[last] == 1.0 );
558
/* Otherwise, consider including 0 and last in the below loop. */
560
for (unsigned i = 1; i < last; i++) {
561
u[i] = NewtonRaphsonRootFind(bezCurve, d[i], u[i]);
566
* Use Newton-Raphson iteration to find better root.
568
* \param Q Current fitted curve
569
* \param P Digitized point
570
* \param u Parameter value for "P"
575
NewtonRaphsonRootFind(BezierCurve const Q, Point const &P, double const u)
580
/* Generate control vertices for Q'. */
582
for (unsigned i = 0; i < 3; i++) {
583
Q1[i] = 3.0 * ( Q[i+1] - Q[i] );
586
/* Generate control vertices for Q''. */
588
for (unsigned i = 0; i < 2; i++) {
589
Q2[i] = 2.0 * ( Q1[i+1] - Q1[i] );
592
/* Compute Q(u), Q'(u) and Q''(u). */
593
Point const Q_u = bezier_pt(3, Q, u);
594
Point const Q1_u = bezier_pt(2, Q1, u);
595
Point const Q2_u = bezier_pt(1, Q2, u);
597
/* Compute f(u)/f'(u), where f is the derivative wrt u of distsq(u) = 0.5 * the square of the
598
distance from P to Q(u). Here we're using Newton-Raphson to find a stationary point in the
599
distsq(u), hopefully corresponding to a local minimum in distsq (and hence a local minimum
600
distance from P to Q(u)). */
601
Point const diff = Q_u - P;
602
double numerator = dot(diff, Q1_u);
603
double denominator = dot(Q1_u, Q1_u) + dot(diff, Q2_u);
606
if ( denominator > 0. ) {
607
/* One iteration of Newton-Raphson:
608
improved_u = u - f(u)/f'(u) */
609
improved_u = u - ( numerator / denominator );
611
/* Using Newton-Raphson would move in the wrong direction (towards a local maximum rather
612
than local minimum), so we move an arbitrary amount in the right direction. */
613
if ( numerator > 0. ) {
614
improved_u = u * .98 - .01;
615
} else if ( numerator < 0. ) {
616
/* Deliberately asymmetrical, to reduce the chance of cycling. */
617
improved_u = .031 + u * .98;
623
if (!IS_FINITE(improved_u)) {
625
} else if ( improved_u < 0.0 ) {
627
} else if ( improved_u > 1.0 ) {
631
/* Ensure that improved_u isn't actually worse. */
633
double const diff_lensq = lensq(diff);
634
for (double proportion = .125; ; proportion += .125) {
635
if ( lensq( bezier_pt(3, Q, improved_u) - P ) > diff_lensq ) {
636
if ( proportion > 1.0 ) {
637
//g_warning("found proportion %g", proportion);
641
improved_u = ( ( 1 - proportion ) * improved_u +
649
DOUBLE_ASSERT(improved_u);
654
* Evaluate a Bezier curve at parameter value \a t.
656
* \param degree The degree of the Bezier curve: 3 for cubic, 2 for quadratic etc. Must be less
658
* \param V The control points for the Bezier curve. Must have (\a degree+1)
660
* \param t The "parameter" value, specifying whereabouts along the curve to
661
* evaluate. Typically in the range [0.0, 1.0].
664
* BezierII(1, V) gives (s, t) * V, i.e. t of the way
666
* BezierII(2, V) gives (s**2, 2*s*t, t**2) * V.
667
* BezierII(3, V) gives (s**3, 3 s**2 t, 3s t**2, t**3) * V.
669
* The derivative of BezierII(i, V) with respect to t
670
* is i * BezierII(i-1, V'), where for all j, V'[j] =
674
bezier_pt(unsigned const degree, Point const V[], double const t)
676
/** Pascal's triangle. */
677
static int const pascal[4][4] = {{1},
682
double const s = 1.0 - t;
684
/* Calculate powers of t and s. */
687
spow[0] = 1.0; spow[1] = s;
688
tpow[0] = 1.0; tpow[1] = t;
689
for (unsigned i = 1; i < degree; ++i) {
690
spow[i + 1] = spow[i] * s;
691
tpow[i + 1] = tpow[i] * t;
694
Point ret = spow[degree] * V[0];
695
for (unsigned i = 1; i <= degree; ++i) {
696
ret += pascal[degree][i] * spow[degree - i] * tpow[i] * V[i];
702
* ComputeLeftTangent, ComputeRightTangent, ComputeCenterTangent :
703
* Approximate unit tangents at endpoints and "center" of digitized curve
707
* Estimate the (forward) tangent at point d[first + 0.5].
709
* Unlike the center and right versions, this calculates the tangent in
710
* the way one might expect, i.e., wrt increasing index into d.
711
* \pre (2 \<= len) and (d[0] != d[1]).
714
darray_left_tangent(Point const d[], unsigned const len)
717
assert( d[0] != d[1] );
718
return unit_vector( d[1] - d[0] );
722
* Estimates the (backward) tangent at d[last - 0.5].
724
* \note The tangent is "backwards", i.e. it is with respect to
725
* decreasing index rather than increasing index.
728
* \pre d[len - 1] != d[len - 2].
729
* \pre all[p in d] in_svg_plane(p).
732
darray_right_tangent(Point const d[], unsigned const len)
735
unsigned const last = len - 1;
736
unsigned const prev = last - 1;
737
assert( d[last] != d[prev] );
738
return unit_vector( d[prev] - d[last] );
742
* Estimate the (forward) tangent at point d[0].
744
* Unlike the center and right versions, this calculates the tangent in
745
* the way one might expect, i.e., wrt increasing index into d.
749
* \pre all[p in d] in_svg_plane(p).
750
* \post is_unit_vector(ret).
753
darray_left_tangent(Point const d[], unsigned const len, double const tolerance_sq)
756
assert( 0 <= tolerance_sq );
757
for (unsigned i = 1;;) {
758
Point const pi(d[i]);
759
Point const t(pi - d[0]);
760
double const distsq = dot(t, t);
761
if ( tolerance_sq < distsq ) {
762
return unit_vector(t);
767
? darray_left_tangent(d, len)
774
* Estimates the (backward) tangent at d[last].
776
* \note The tangent is "backwards", i.e. it is with respect to
777
* decreasing index rather than increasing index.
780
* \pre d[len - 1] != d[len - 2].
781
* \pre all[p in d] in_svg_plane(p).
784
darray_right_tangent(Point const d[], unsigned const len, double const tolerance_sq)
787
assert( 0 <= tolerance_sq );
788
unsigned const last = len - 1;
789
for (unsigned i = last - 1;; i--) {
790
Point const pi(d[i]);
791
Point const t(pi - d[last]);
792
double const distsq = dot(t, t);
793
if ( tolerance_sq < distsq ) {
794
return unit_vector(t);
798
? darray_right_tangent(d, len)
805
* Estimates the (backward) tangent at d[center], by averaging the two
806
* segments connected to d[center] (and then normalizing the result).
808
* \note The tangent is "backwards", i.e. it is with respect to
809
* decreasing index rather than increasing index.
811
* \pre (0 \< center \< len - 1) and d is uniqued (at least in
812
* the immediate vicinity of \a center).
815
darray_center_tangent(Point const d[],
816
unsigned const center,
819
assert( center != 0 );
820
assert( center < len - 1 );
823
if ( d[center + 1] == d[center - 1] ) {
824
/* Rotate 90 degrees in an arbitrary direction. */
825
Point const diff = d[center] - d[center - 1];
828
ret = d[center - 1] - d[center + 1];
836
* Assign parameter values to digitized points using relative distances between points.
838
* \pre Parameter array u must have space for \a len items.
841
chord_length_parameterize(Point const d[], double u[], unsigned const len)
846
/* First let u[i] equal the distance travelled along the path from d[0] to d[i]. */
848
for (unsigned i = 1; i < len; i++) {
849
double const dist = distance(d[i], d[i-1]);
850
u[i] = u[i-1] + dist;
853
/* Then scale to [0.0 .. 1.0]. */
854
double tot_len = u[len - 1];
855
if(!( tot_len != 0 ))
857
if (IS_FINITE(tot_len)) {
858
for (unsigned i = 1; i < len; ++i) {
862
/* We could do better, but this probably never happens anyway. */
863
for (unsigned i = 1; i < len; ++i) {
864
u[i] = i / (double) ( len - 1 );
869
* It's been reported that u[len - 1] can differ from 1.0 on some
870
* systems (amd64), despite it having been calculated as x / x where x
871
* is isFinite and non-zero.
873
if (u[len - 1] != 1) {
874
double const diff = u[len - 1] - 1;
875
if (fabs(diff) > 1e-13) {
876
assert(0); // No warnings in 2geom
877
//g_warning("u[len - 1] = %19g (= 1 + %19g), expecting exactly 1",
878
// u[len - 1], diff);
884
assert( u[0] == 0.0 && u[len - 1] == 1.0 );
885
for (unsigned i = 1; i < len; i++) {
886
assert( u[i] >= u[i-1] );
895
* Find the maximum squared distance of digitized points to fitted curve, and (if this maximum
896
* error is non-zero) set \a *splitPoint to the corresponding index.
900
* \pre u[len - 1] == 1.0.
901
* \post ((ret == 0.0)
902
* || ((*splitPoint \< len - 1)
903
* \&\& (*splitPoint != 0 || ret \< 0.0))).
906
compute_max_error_ratio(Point const d[], double const u[], unsigned const len,
907
BezierCurve const bezCurve, double const tolerance,
908
unsigned *const splitPoint)
911
unsigned const last = len - 1;
912
assert( bezCurve[0] == d[0] );
913
assert( bezCurve[3] == d[last] );
914
assert( u[0] == 0.0 );
915
assert( u[last] == 1.0 );
916
/* I.e. assert that the error for the first & last points is zero.
917
* Otherwise we should include those points in the below loop.
918
* The assertion is also necessary to ensure 0 < splitPoint < last.
921
double maxDistsq = 0.0; /* Maximum error */
922
double max_hook_ratio = 0.0;
923
unsigned snap_end = 0;
924
Point prev = bezCurve[0];
925
for (unsigned i = 1; i <= last; i++) {
926
Point const curr = bezier_pt(3, bezCurve, u[i]);
927
double const distsq = lensq( curr - d[i] );
928
if ( distsq > maxDistsq ) {
932
double const hook_ratio = compute_hook(prev, curr, .5 * (u[i - 1] + u[i]), bezCurve, tolerance);
933
if (max_hook_ratio < hook_ratio) {
934
max_hook_ratio = hook_ratio;
940
double const dist_ratio = sqrt(maxDistsq) / tolerance;
942
if (max_hook_ratio <= dist_ratio) {
945
assert(0 < snap_end);
946
ret = -max_hook_ratio;
947
*splitPoint = snap_end - 1;
950
|| ( ( *splitPoint < last )
951
&& ( *splitPoint != 0 || ret < 0. ) ) );
956
* Whereas compute_max_error_ratio() checks for itself that each data point
957
* is near some point on the curve, this function checks that each point on
958
* the curve is near some data point (or near some point on the polyline
959
* defined by the data points, or something like that: we allow for a
960
* "reasonable curviness" from such a polyline). "Reasonable curviness"
961
* means we draw a circle centred at the midpoint of a..b, of radius
962
* proportional to the length |a - b|, and require that each point on the
963
* segment of bezCurve between the parameters of a and b be within that circle.
964
* If any point P on the bezCurve segment is outside of that allowable
965
* region (circle), then we return some metric that increases with the
966
* distance from P to the circle.
968
* Given that this is a fairly arbitrary criterion for finding appropriate
969
* places for sharp corners, we test only one point on bezCurve, namely
970
* the point on bezCurve with parameter halfway between our estimated
971
* parameters for a and b. (Alternatives are taking the farthest of a
972
* few parameters between those of a and b, or even using a variant of
973
* NewtonRaphsonFindRoot() for finding the maximum rather than minimum
977
compute_hook(Point const &a, Point const &b, double const u, BezierCurve const bezCurve,
978
double const tolerance)
980
Point const P = bezier_pt(3, bezCurve, u);
981
double const dist = distance((a+b)*.5, P);
982
if (dist < tolerance) {
985
double const allowed = distance(a, b) + tolerance;
986
return dist / allowed;
988
* effic: Hooks are very rare. We could start by comparing
989
* distsq, only resorting to the more expensive L2 in cases of
999
c-file-style:"stroustrup"
1000
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
1001
indent-tabs-mode:nil
1005
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
added
removed
inkscape-0.47pre1/src/2geom/bezier-utils.cpp
