40742313779ee5e43be93a9191f1c86412cf183bKrzysztof Kosiński/* Bezier interpolation for inkscape drawing code.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Original code published in:
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * An Algorithm for Automatically Fitting Digitized Curves
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * by Philip J. Schneider
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * "Graphics Gems", Academic Press, 1990
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Philip J. Schneider
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Lauris Kaplinski <lauris@kaplinski.com>
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Peter Moulder <pmoulder@mail.csse.monash.edu.au>
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Copyright (C) 1990 Philip J. Schneider
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Copyright (C) 2001 Lauris Kaplinski
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Copyright (C) 2001 Ximian, Inc.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Copyright (C) 2003,2004 Monash University
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * This library is free software; you can redistribute it and/or
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * modify it either under the terms of the GNU Lesser General Public
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * License version 2.1 as published by the Free Software Foundation
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * (the "LGPL") or, at your option, under the terms of the Mozilla
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Public License Version 1.1 (the "MPL"). If you do not alter this
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * notice, a recipient may use your version of this file under either
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * the MPL or the LGPL.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * You should have received a copy of the LGPL along with this library
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * in the file COPYING-LGPL-2.1; if not, write to the Free Software
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * You should have received a copy of the MPL along with this library
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * in the file COPYING-MPL-1.1
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * The contents of this file are subject to the Mozilla Public License
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Version 1.1 (the "License"); you may not use this file except in
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * compliance with the License. You may obtain a copy of the License at
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * OF ANY KIND, either express or implied. See the LGPL or the MPL for
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * the specific language governing rights and limitations.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen/* Forward declarations */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstatic void generate_bezier(Point b[], Point const d[], double const u[], unsigned len,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Point const &tHat1, Point const &tHat2, double tolerance_sq);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Point const data[], double const u[], unsigned len,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstatic void estimate_bi(Point b[4], unsigned ei,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Point const data[], double const u[], unsigned len);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstatic void reparameterize(Point const d[], unsigned len, double u[], Point const bezCurve[]);
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstatic double NewtonRaphsonRootFind(Point const Q[], Point const &P, double u);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstatic Point darray_center_tangent(Point const d[], unsigned center, unsigned length);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstatic Point darray_right_tangent(Point const d[], unsigned const len);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstatic unsigned copy_without_nans_or_adjacent_duplicates(Point const src[], unsigned src_len, Point dest[]);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstatic void chord_length_parameterize(Point const d[], double u[], unsigned len);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstatic double compute_max_error_ratio(Point const d[], double const u[], unsigned len,
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelenstatic double compute_hook(Point const &a, Point const &b, double const u, Point const bezCurve[],
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double const tolerance);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * B0, B1, B2, B3 : Bezier multipliers
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#define B0(u) ( ( 1.0 - u ) * ( 1.0 - u ) * ( 1.0 - u ) )
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#define B1(u) ( 3 * u * ( 1.0 - u ) * ( 1.0 - u ) )
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen# define DOUBLE_ASSERT(x) assert( ( (x) > -SP_HUGE ) && ( (x) < SP_HUGE ) )
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen DOUBLE_ASSERT((b)[0][X]); DOUBLE_ASSERT((b)[0][Y]); \
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen DOUBLE_ASSERT((b)[1][X]); DOUBLE_ASSERT((b)[1][Y]); \
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen DOUBLE_ASSERT((b)[2][X]); DOUBLE_ASSERT((b)[2][Y]); \
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen DOUBLE_ASSERT((b)[3][X]); DOUBLE_ASSERT((b)[3][Y]); \
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Fit a single-segment Bezier curve to a set of digitized points.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \return Number of segments generated, or -1 on error.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenbezier_fit_cubic(Point *bezier, Point const *data, int len, double error)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen return bezier_fit_cubic_r(bezier, data, len, error, 1);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Fit a multi-segment Bezier curve to a set of digitized points, with
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * possible weedout of identical points and NaNs.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param max_beziers Maximum number of generated segments
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param Result array, must be large enough for n. segments * 4 elements.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \return Number of segments generated, or -1 on error.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenbezier_fit_cubic_r(Point bezier[], Point const data[], int const len, double const error, unsigned const max_beziers)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen unsigned uniqued_len = copy_without_nans_or_adjacent_duplicates(data, len, uniqued_data);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Call fit-cubic function with recursion. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen int const ret = bezier_fit_cubic_full(bezier, NULL, uniqued_data, uniqued_len,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Copy points from src to dest, filter out points containing NaN and
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * adjacent points with equal x and y.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \return length of dest
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstatic unsigned
981b809bc6ed10a21e89444d9447e5475801874fjohanengelencopy_without_nans_or_adjacent_duplicates(Point const src[], unsigned src_len, Point dest[])
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Fit a multi-segment Bezier curve to a set of digitized points, without
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * possible weedout of identical points and NaNs.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre data is uniqued, i.e. not exist i: data[i] == data[i + 1].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param max_beziers Maximum number of generated segments
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param Result array, must be large enough for n. segments * 4 elements.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenbezier_fit_cubic_full(Point bezier[], int split_points[],
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double const error, unsigned const max_beziers)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* We have 2 points, which can be fitted trivially. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double const dist = distance(bezier[0], bezier[3]) / 3.0;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Numerical problem, fall back to straight line segment. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Parameterize points, and attempt to fit curve */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen unsigned splitPoint; /* Point to split point set at. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double *u = new double[len];
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Zero-length path: every point in data[] is the same.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * (Clients aren't allowed to pass such data; handling the case is defensive
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * programming.)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen generate_bezier(bezier, data, u, len, tHat1, tHat2, error);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Find max deviation of points to fitted curve. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double maxErrorRatio = compute_max_error_ratio(data, u, len, bezier, tolerance, &splitPoint);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* If error not too large, then try some reparameterization and iteration. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if ( 0.0 <= maxErrorRatio && maxErrorRatio <= 3.0 ) {
c8589a6c7367d09fa756755cef0dd448c7328a71Johan B. C. Engelen int const maxIterations = 4; /* std::max times to try iterating */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen for (int i = 0; i < maxIterations; i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen generate_bezier(bezier, data, u, len, tHat1, tHat2, error);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen maxErrorRatio = compute_max_error_ratio(data, u, len, bezier, tolerance, &splitPoint);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Got spike even with unconstrained initial tangent. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen return bezier_fit_cubic_full(bezier, split_points, data, len, unconstrained_tangent, tHat2,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Got spike even with unconstrained final tangent. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen return bezier_fit_cubic_full(bezier, split_points, data, len, tHat1, unconstrained_tangent,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Fitting failed -- split at max error point and fit recursively
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen unsigned const rec_max_beziers1 = max_beziers - 1;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if(!(0 < splitPoint && splitPoint < unsigned(len - 1)))
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Unit tangent vector at splitPoint. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen recTHat2 = darray_center_tangent(data, splitPoint, len);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen int const nsegs1 = bezier_fit_cubic_full(bezier, split_points, data, splitPoint + 1,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen g_print("fit_cubic[1]: recursive call failed\n");
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen unsigned const rec_max_beziers2 = max_beziers - nsegs1;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen int const nsegs2 = bezier_fit_cubic_full(bezier + nsegs1*4,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen g_print("fit_cubic[2]: recursive call failed\n");
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen g_print("fit_cubic: success[nsegs: %d+%d=%d] on max_beziers:%u\n",
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Fill in \a bezier[] based on the given data and tangent requirements, using
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * a least-squares fit.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Each of tHat1 and tHat2 should be either a zero vector or a unit vector.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * If it is zero, then bezier[1 or 2] is estimated without constraint; otherwise,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * it bezier[1 or 2] is placed in the specified direction from bezier[0 or 3].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param tolerance_sq Used only for an initial guess as to tangent directions
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * when \a tHat1 or \a tHat2 is zero.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Point const data[], double const u[], unsigned const len,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen ? darray_right_tangent(data, len, tolerance_sq)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen estimate_lengths(bezier, data, u, len, est_tHat1, est_tHat2);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* We find that darray_right_tangent tends to produce better results
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen for our current freehand tool than full estimation. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen est_tHat1 = unit_vector(bezier[1] - bezier[0]);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen estimate_lengths(bezier, data, u, len, est_tHat1, est_tHat2);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Point const data[], double const uPrime[], unsigned const len,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Create the C and X matrices. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen C[0][0] = 0.0;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* First and last control points of the Bezier curve are positioned exactly at the first and
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen last data points. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen for (unsigned i = 0; i < len; i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Bezier control point coefficients. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* rhs for eqn */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Additional offset to the data point from the predicted point if we were to set bezier[1]
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen to bezier[0] and bezier[2] to bezier[3]. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* We've constructed a pair of equations in the form of a matrix product C * alpha = X.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Now solve for alpha. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Compute the determinants of C and X. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double const det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Apparently Kramer's rule. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double const det_C0_X = C[0][0] * X[1] - C[0][1] * X[0];
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double const det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1];
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* The matrix is under-determined. Try requiring alpha_l == alpha_r.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * One way of implementing the constraint alpha_l == alpha_r is to treat them as the same
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * variable in the equations. We can do this by adding the columns of C to form a single
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * column, to be multiplied by alpha to give the column vector X.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * We try each row in turn.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Let the below code handle this. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* If alpha negative, use the Wu/Barsky heuristic (see text). (If alpha is 0, you get
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen coincident control points that lead to divide by zero in any subsequent
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen NewtonRaphsonRootFind() call.) */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /// \todo Check whether this special-casing is necessary now that
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /// NewtonRaphsonRootFind handles non-positive denominator.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen alpha_l = alpha_r = distance(data[0], data[len-1]) / 3.0;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Control points 1 and 2 are positioned an alpha distance out on the tangent vectors, left and
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen right, respectively. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen return dot(p, p);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Point const data[], double const u[], unsigned const len)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double den = 0.;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen for (unsigned i = 0; i < len; ++i) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double const ui = u[i];
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double const b[4] = {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen for (unsigned d = 0; d < 2; ++d) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if (den != 0.) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen for (unsigned d = 0; d < 2; ++d) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen bezier[ei] = ( oi * bezier[0] + ei * bezier[3] ) / 3.;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Given set of points and their parameterization, try to find a better assignment of parameter
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * values for the points.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param d Array of digitized points.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param u Current parameter values.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param bezCurve Current fitted curve.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param len Number of values in both d and u arrays.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Also the size of the array that is allocated for return.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen unsigned const len,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Otherwise, consider including 0 and last in the below loop. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen u[i] = NewtonRaphsonRootFind(bezCurve, d[i], u[i]);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Use Newton-Raphson iteration to find better root.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param Q Current fitted curve
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param P Digitized point
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param u Parameter value for "P"
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \return Improved u
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstatic double
d37634d73670180f99a3e0ea583621373d90ec4fJohan EngelenNewtonRaphsonRootFind(Point const Q[], Point const &P, double const u)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Generate control vertices for Q'. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen for (unsigned i = 0; i < 3; i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Generate control vertices for Q''. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen for (unsigned i = 0; i < 2; i++) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Compute Q(u), Q'(u) and Q''(u). */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Compute f(u)/f'(u), where f is the derivative wrt u of distsq(u) = 0.5 * the square of the
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen distance from P to Q(u). Here we're using Newton-Raphson to find a stationary point in the
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen distsq(u), hopefully corresponding to a local minimum in distsq (and hence a local minimum
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen distance from P to Q(u)). */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double denominator = dot(Q1_u, Q1_u) + dot(diff, Q2_u);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* One iteration of Newton-Raphson:
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen improved_u = u - f(u)/f'(u) */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Using Newton-Raphson would move in the wrong direction (towards a local maximum rather
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen than local minimum), so we move an arbitrary amount in the right direction. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen } else if ( numerator < 0. ) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Deliberately asymmetrical, to reduce the chance of cycling. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Ensure that improved_u isn't actually worse. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen for (double proportion = .125; ; proportion += .125) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if ( lensq( bezier_pt(3, Q, improved_u) - P ) > diff_lensq ) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //g_warning("found proportion %g", proportion);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen improved_u = ( ( 1 - proportion ) * improved_u +
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Evaluate a Bezier curve at parameter value \a t.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param degree The degree of the Bezier curve: 3 for cubic, 2 for quadratic etc. Must be less
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param V The control points for the Bezier curve. Must have (\a degree+1)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \param t The "parameter" value, specifying whereabouts along the curve to
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * evaluate. Typically in the range [0.0, 1.0].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Let s = 1 - t.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * BezierII(1, V) gives (s, t) * V, i.e. t of the way
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * from V[0] to V[1].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * BezierII(2, V) gives (s**2, 2*s*t, t**2) * V.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * BezierII(3, V) gives (s**3, 3 s**2 t, 3s t**2, t**3) * V.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * The derivative of BezierII(i, V) with respect to t
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * is i * BezierII(i-1, V'), where for all j, V'[j] =
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * V[j + 1] - V[j].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenbezier_pt(unsigned const degree, Point const V[], double const t)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /** Pascal's triangle. */
76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński static int const pascal[4][4] = {{1, 0, 0, 0},
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double const s = 1.0 - t;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Calculate powers of t and s. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen ret += pascal[degree][i] * spow[degree - i] * tpow[i] * V[i];
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * ComputeLeftTangent, ComputeRightTangent, ComputeCenterTangent :
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Approximate unit tangents at endpoints and "center" of digitized curve
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Estimate the (forward) tangent at point d[first + 0.5].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Unlike the center and right versions, this calculates the tangent in
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * the way one might expect, i.e., wrt increasing index into d.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre (2 \<= len) and (d[0] != d[1]).
981b809bc6ed10a21e89444d9447e5475801874fjohanengelendarray_left_tangent(Point const d[], unsigned const len)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Estimates the (backward) tangent at d[last - 0.5].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \note The tangent is "backwards", i.e. it is with respect to
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * decreasing index rather than increasing index.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre 2 \<= len.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre d[len - 1] != d[len - 2].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre all[p in d] in_svg_plane(p).
981b809bc6ed10a21e89444d9447e5475801874fjohanengelendarray_right_tangent(Point const d[], unsigned const len)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Estimate the (forward) tangent at point d[0].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Unlike the center and right versions, this calculates the tangent in
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * the way one might expect, i.e., wrt increasing index into d.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre 2 \<= len.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre d[0] != d[1].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre all[p in d] in_svg_plane(p).
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \post is_unit_vector(ret).
981b809bc6ed10a21e89444d9447e5475801874fjohanengelendarray_left_tangent(Point const d[], unsigned const len, double const tolerance_sq)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen for (unsigned i = 1;;) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen return ( distsq == 0
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Estimates the (backward) tangent at d[last].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \note The tangent is "backwards", i.e. it is with respect to
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * decreasing index rather than increasing index.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre 2 \<= len.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre d[len - 1] != d[len - 2].
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre all[p in d] in_svg_plane(p).
981b809bc6ed10a21e89444d9447e5475801874fjohanengelendarray_right_tangent(Point const d[], unsigned const len, double const tolerance_sq)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if (i == 0) {
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen return ( distsq == 0
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Estimates the (backward) tangent at d[center], by averaging the two
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * segments connected to d[center] (and then normalizing the result).
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \note The tangent is "backwards", i.e. it is with respect to
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * decreasing index rather than increasing index.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre (0 \< center \< len - 1) and d is uniqued (at least in
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * the immediate vicinity of \a center).
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen unsigned const center,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen unsigned const len)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Rotate 90 degrees in an arbitrary direction. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Assign parameter values to digitized points using relative distances between points.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre Parameter array u must have space for \a len items.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenchord_length_parameterize(Point const d[], double u[], unsigned const len)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* First let u[i] equal the distance travelled along the path from d[0] to d[i]. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Then scale to [0.0 .. 1.0]. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if(!( tot_len != 0 ))
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* We could do better, but this probably never happens anyway. */
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * It's been reported that u[len - 1] can differ from 1.0 on some
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * systems (amd64), despite it having been calculated as x / x where x
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * is isFinite and non-zero.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //g_warning("u[len - 1] = %19g (= 1 + %19g), expecting exactly 1",
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen // u[len - 1], diff);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Find the maximum squared distance of digitized points to fitted curve, and (if this maximum
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * error is non-zero) set \a *splitPoint to the corresponding index.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre 2 \<= len.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre u[0] == 0.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \pre u[len - 1] == 1.0.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \post ((ret == 0.0)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * || ((*splitPoint \< len - 1)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * \&\& (*splitPoint != 0 || ret \< 0.0))).
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstatic double
981b809bc6ed10a21e89444d9447e5475801874fjohanengelencompute_max_error_ratio(Point const d[], double const u[], unsigned const len,
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelen Point const bezCurve[], double const tolerance,
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen unsigned *const splitPoint)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* I.e. assert that the error for the first & last points is zero.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Otherwise we should include those points in the below loop.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * The assertion is also necessary to ensure 0 < splitPoint < last.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen unsigned snap_end = 0;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Point const curr = bezier_pt(3, bezCurve, u[i]);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double const hook_ratio = compute_hook(prev, curr, .5 * (u[i - 1] + u[i]), bezCurve, tolerance);
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double const dist_ratio = sqrt(maxDistsq) / tolerance;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Whereas compute_max_error_ratio() checks for itself that each data point
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * is near some point on the curve, this function checks that each point on
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * the curve is near some data point (or near some point on the polyline
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * defined by the data points, or something like that: we allow for a
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * "reasonable curviness" from such a polyline). "Reasonable curviness"
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * means we draw a circle centred at the midpoint of a..b, of radius
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * proportional to the length |a - b|, and require that each point on the
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * segment of bezCurve between the parameters of a and b be within that circle.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * If any point P on the bezCurve segment is outside of that allowable
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * region (circle), then we return some metric that increases with the
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * distance from P to the circle.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * Given that this is a fairly arbitrary criterion for finding appropriate
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * places for sharp corners, we test only one point on bezCurve, namely
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * the point on bezCurve with parameter halfway between our estimated
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * parameters for a and b. (Alternatives are taking the farthest of a
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * few parameters between those of a and b, or even using a variant of
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * NewtonRaphsonFindRoot() for finding the maximum rather than minimum
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * distance.)
981b809bc6ed10a21e89444d9447e5475801874fjohanengelenstatic double
d37634d73670180f99a3e0ea583621373d90ec4fJohan Engelencompute_hook(Point const &a, Point const &b, double const u, Point const bezCurve[],
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double const allowed = distance(a, b) + tolerance;
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * effic: Hooks are very rare. We could start by comparing
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * distsq, only resorting to the more expensive L2 in cases of
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen * uncertainty.
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Local Variables:
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen c-file-style:"stroustrup"
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen indent-tabs-mode:nil
981b809bc6ed10a21e89444d9447e5475801874fjohanengelen fill-column:99
a4030d5ca449e7e384bc699cd249ee704faaeab0Chris Morgan// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :