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