diff options
Diffstat (limited to '')
-rw-r--r-- | src/path/splinefit/bezier-fit.cpp | 90 | ||||
-rw-r--r-- | src/path/splinefit/bezier-fit.h | 20 | ||||
-rw-r--r-- | src/path/splinefit/splinefit.c | 1427 | ||||
-rw-r--r-- | src/path/splinefit/splinefit.h | 78 | ||||
-rw-r--r-- | src/path/splinefit/splinefont.c | 1174 | ||||
-rw-r--r-- | src/path/splinefit/splinefont.h | 191 | ||||
-rw-r--r-- | src/path/splinefit/splinerefigure.c | 117 | ||||
-rw-r--r-- | src/path/splinefit/splinerefigure.h | 9 |
8 files changed, 3106 insertions, 0 deletions
diff --git a/src/path/splinefit/bezier-fit.cpp b/src/path/splinefit/bezier-fit.cpp new file mode 100644 index 0000000..111f538 --- /dev/null +++ b/src/path/splinefit/bezier-fit.cpp @@ -0,0 +1,90 @@ +// SPDX-License-Identifier: GPL-2.0-or-later + +#include <iostream> +#include <vector> +#include "bezier-fit.h" +#include <2geom/bezier-utils.h> +#include <2geom/point.h> + +extern "C" { + #include "splinefit.h" + #include "splinefont.h" +} + +int bezier_fit(Geom::Point bezier[4], const std::vector<InputPoint>& data) { + + if (data.size() <= 2) return 0; + + int order2 = false; // not 2nd order, so cubic + // "Fitting cubic Bézier curves" + // https://raphlinus.github.io/curves/2021/03/11/bezier-fitting.html + mergetype mt = mt_levien; + auto len = data.size(); + + std::vector<FitPoint> fit; + for (int i = 0; i < len; ++i) { + fit.push_back({}); + auto& fp = fit.back(); + fp.p.x = data[i].x(); + fp.p.y = data[i].y(); + fp.t = data[i].t; + fp.ut.x = fp.ut.y = 0; + } + + // transform data into spline set format + + auto input = (SplineSet*)chunkalloc(sizeof(SplineSet)); + + for (int i = 0; i < len; ++i) { + auto& d = data[i]; + auto sp = SplinePointCreate(d.x(), d.y()); + if (d.have_slope) { + sp->nextcp.x = d.front.x(); + sp->nextcp.y = d.front.y(); + sp->nonextcp = false; + sp->prevcp.x = d.back.x(); + sp->prevcp.y = d.back.y(); + sp->noprevcp = false; + } + + if (i == 0) { + input->first = input->last = sp; + } + else { + SplineMake(input->last, sp, order2); + input->last = sp; + } + } + + Spline* spline = ApproximateSplineFromPointsSlopes(input->first, input->last, fit.data(), fit.size(), order2, mt); + bool ok = spline != nullptr; + + if (!spline) { + std::vector<Geom::Point> inp; + inp.reserve(data.size()); + for (auto& pt : data) { + inp.push_back(pt); + } + ok = bezier_fit_cubic(bezier, inp.data(), inp.size(), 0.5) > 0; + } + + if (spline) { + bezier[0].x() = spline->from->me.x; + bezier[0].y() = spline->from->me.y; + + bezier[1].x() = spline->from->nextcp.x; + bezier[1].y() = spline->from->nextcp.y; + + bezier[2].x() = spline->to->prevcp.x; + bezier[2].y() = spline->to->prevcp.y; + + bezier[3].x() = spline->to->me.x; + bezier[3].y() = spline->to->me.y; + } + + SplinePointListFree(input); + //TODO: verify that all C structs are freed up + // SplineFree(spline); + + return ok; +} diff --git a/src/path/splinefit/bezier-fit.h b/src/path/splinefit/bezier-fit.h new file mode 100644 index 0000000..7eb0440 --- /dev/null +++ b/src/path/splinefit/bezier-fit.h @@ -0,0 +1,20 @@ +// SPDX-License-Identifier: GPL-2.0-or-later + +#include <vector> +#include "2geom/point.h" + +struct InputPoint : Geom::Point { + InputPoint() {} + InputPoint(const Geom::Point& pt) : Point(pt) {} + InputPoint(const Geom::Point& pt, double t) : Point(pt), t(t) {} + InputPoint(const Geom::Point& pt, const Geom::Point& front, const Geom::Point& back, double t) + : Point(pt), front(front), back(back), t(t), have_slope(true) {} + + Geom::Point front; + Geom::Point back; + double t = 0; + bool have_slope = false; +}; + +// Fit cubic Bezier to input points; use slope of the first and last points to find a fit +int bezier_fit(Geom::Point bezier[4], const std::vector<InputPoint>& data); diff --git a/src/path/splinefit/splinefit.c b/src/path/splinefit/splinefit.c new file mode 100644 index 0000000..9e3039a --- /dev/null +++ b/src/path/splinefit/splinefit.c @@ -0,0 +1,1427 @@ +// SPDX-License-Identifier: GPL-2.0-or-later +/* -*- coding: utf-8 -*- */ +/* Copyright (C) 2000-2012 by George Williams, 2021 by Linus Romer */ +/* + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + + * The name of the author may not be used to endorse or promote products + * derived from this software without specific prior written permission. + + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED + * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO + * EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; + * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, + * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR + * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF + * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/* Incorporated into Inkscape sources; splinefit.c as of 2023-02-15. + + * Note: The only changes to this file in Inkscape codebase are modification of includes + * and adding this note. Formatting is intact to facilitate future merges. + */ + +#include <memory.h> +#include <stdbool.h> +#include <stdint.h> +#include "splinefit.h" + +#include <math.h> + +static Spline *IsLinearApprox(SplinePoint *from, SplinePoint *to, + FitPoint *mid, int cnt, int order2) { + bigreal vx, vy, slope; + int i; + + vx = to->me.x-from->me.x; vy = to->me.y-from->me.y; + if ( vx==0 && vy==0 ) { + for ( i=0; i<cnt; ++i ) + if ( mid[i].p.x != from->me.x || mid[i].p.y != from->me.y ) +return( NULL ); + } else if ( fabs(vx)>fabs(vy) ) { + slope = vy/vx; + for ( i=0; i<cnt; ++i ) + if ( !RealWithin(mid[i].p.y,from->me.y+slope*(mid[i].p.x-from->me.x),.7) ) +return( NULL ); + } else { + slope = vx/vy; + for ( i=0; i<cnt; ++i ) + if ( !RealWithin(mid[i].p.x,from->me.x+slope*(mid[i].p.y-from->me.y),.7) ) +return( NULL ); + } + from->nextcp = from->me; + to->prevcp = to->me; +return( SplineMake(from,to,order2) ); +} + +/* This routine should almost never be called now. It uses a flawed algorithm */ +/* which won't produce the best results. It gets called only when the better */ +/* approach doesn't work (singular matrices, etc.) */ +/* Old comment, back when I was confused... */ +/* Least squares tells us that: + | S(xi*ti^3) | | S(ti^6) S(ti^5) S(ti^4) S(ti^3) | | a | + | S(xi*ti^2) | = | S(ti^5) S(ti^4) S(ti^3) S(ti^2) | * | b | + | S(xi*ti) | | S(ti^4) S(ti^3) S(ti^2) S(ti) | | c | + | S(xi) | | S(ti^3) S(ti^2) S(ti) n | | d | + and the definition of a spline tells us: + | x1 | = | 1 1 1 1 | * (a b c d) + | x0 | = | 0 0 0 1 | * (a b c d) +So we're a bit over specified. Let's use the last two lines of least squares +and the 2 from the spline defn. So d==x0. Now we've got three unknowns +and only three equations... + +For order2 splines we've got + | S(xi*ti^2) | | S(ti^4) S(ti^3) S(ti^2) | | b | + | S(xi*ti) | = | S(ti^3) S(ti^2) S(ti) | * | c | + | S(xi) | | S(ti^2) S(ti) n | | d | + and the definition of a spline tells us: + | x1 | = | 1 1 1 | * (b c d) + | x0 | = | 0 0 1 | * (b c d) +=> + d = x0 + b+c = x1-x0 + S(ti^2)*b + S(ti)*c = S(xi)-n*x0 + S(ti^2)*b + S(ti)*(x1-x0-b) = S(xi)-n*x0 + [ S(ti^2)-S(ti) ]*b = S(xi)-S(ti)*(x1-x0) - n*x0 +*/ +static int _ApproximateSplineFromPoints(SplinePoint *from, SplinePoint *to, + FitPoint *mid, int cnt, BasePoint *nextcp, BasePoint *prevcp, + int order2) { + bigreal tt, ttn; + int i, j, ret; + bigreal vx[3], vy[3], m[3][3]; + bigreal ts[7], xts[4], yts[4]; + BasePoint nres, pres; + int nrescnt=0, prescnt=0; + bigreal nmin, nmax, pmin, pmax, test, ptest; + bigreal bx, by, cx, cy; + + memset(&nres,0,sizeof(nres)); memset(&pres,0,sizeof(pres)); + + /* Add the initial and end points */ + ts[0] = 2; for ( i=1; i<7; ++i ) ts[i] = 1; + xts[0] = from->me.x+to->me.x; yts[0] = from->me.y+to->me.y; + xts[3] = xts[2] = xts[1] = to->me.x; yts[3] = yts[2] = yts[1] = to->me.y; + nmin = pmin = 0; nmax = pmax = (to->me.x-from->me.x)*(to->me.x-from->me.x)+(to->me.y-from->me.y)*(to->me.y-from->me.y); + for ( i=0; i<cnt; ++i ) { + xts[0] += mid[i].p.x; + yts[0] += mid[i].p.y; + ++ts[0]; + tt = mid[i].t; + xts[1] += tt*mid[i].p.x; + yts[1] += tt*mid[i].p.y; + ts[1] += tt; + ts[2] += (ttn=tt*tt); + xts[2] += ttn*mid[i].p.x; + yts[2] += ttn*mid[i].p.y; + ts[3] += (ttn*=tt); + xts[3] += ttn*mid[i].p.x; + yts[3] += ttn*mid[i].p.y; + ts[4] += (ttn*=tt); + ts[5] += (ttn*=tt); + ts[6] += (ttn*=tt); + + test = (mid[i].p.x-from->me.x)*(to->me.x-from->me.x) + (mid[i].p.y-from->me.y)*(to->me.y-from->me.y); + if ( test<nmin ) nmin=test; + if ( test>nmax ) nmax=test; + test = (mid[i].p.x-to->me.x)*(from->me.x-to->me.x) + (mid[i].p.y-to->me.y)*(from->me.y-to->me.y); + if ( test<pmin ) pmin=test; + if ( test>pmax ) pmax=test; + } + pmin *= 1.2; pmax *= 1.2; nmin *= 1.2; nmax *= 1.2; + + for ( j=0; j<3; ++j ) { + if ( order2 ) { + if ( RealNear(ts[j+2],ts[j+1]) ) + continue; + /* This produces really bad results!!!! But I don't see what I can do to improve it */ + bx = (xts[j]-ts[j+1]*(to->me.x-from->me.x) - ts[j]*from->me.x) / (ts[j+2]-ts[j+1]); + by = (yts[j]-ts[j+1]*(to->me.y-from->me.y) - ts[j]*from->me.y) / (ts[j+2]-ts[j+1]); + cx = to->me.x-from->me.x-bx; + cy = to->me.y-from->me.y-by; + + nextcp->x = from->me.x + cx/2; + nextcp->y = from->me.y + cy/2; + *prevcp = *nextcp; + } else { + vx[0] = xts[j+1]-ts[j+1]*from->me.x; + vx[1] = xts[j]-ts[j]*from->me.x; + vx[2] = to->me.x-from->me.x; /* always use the defn of spline */ + + vy[0] = yts[j+1]-ts[j+1]*from->me.y; + vy[1] = yts[j]-ts[j]*from->me.y; + vy[2] = to->me.y-from->me.y; + + m[0][0] = ts[j+4]; m[0][1] = ts[j+3]; m[0][2] = ts[j+2]; + m[1][0] = ts[j+3]; m[1][1] = ts[j+2]; m[1][2] = ts[j+1]; + m[2][0] = 1; m[2][1] = 1; m[2][2] = 1; + + /* Remove a terms from rows 0 and 1 */ + vx[0] -= ts[j+4]*vx[2]; + vy[0] -= ts[j+4]*vy[2]; + m[0][0] = 0; m[0][1] -= ts[j+4]; m[0][2] -= ts[j+4]; + vx[1] -= ts[j+3]*vx[2]; + vy[1] -= ts[j+3]*vy[2]; + m[1][0] = 0; m[1][1] -= ts[j+3]; m[1][2] -= ts[j+3]; + + if ( fabs(m[1][1])<fabs(m[0][1]) ) { + bigreal temp; + temp = vx[1]; vx[1] = vx[0]; vx[0] = temp; + temp = vy[1]; vy[1] = vy[0]; vy[0] = temp; + temp = m[1][1]; m[1][1] = m[0][1]; m[0][1] = temp; + temp = m[1][2]; m[1][2] = m[0][2]; m[0][2] = temp; + } + /* remove b terms from rows 0 and 2 (first normalize row 1 so m[1][1] is 1*/ + vx[1] /= m[1][1]; + vy[1] /= m[1][1]; + m[1][2] /= m[1][1]; + m[1][1] = 1; + vx[0] -= m[0][1]*vx[1]; + vy[0] -= m[0][1]*vy[1]; + m[0][2] -= m[0][1]*m[1][2]; m[0][1] = 0; + vx[2] -= m[2][1]*vx[1]; + vy[2] -= m[2][1]*vy[1]; + m[2][2] -= m[2][1]*m[1][2]; m[2][1] = 0; + + vx[0] /= m[0][2]; /* This is cx */ + vy[0] /= m[0][2]; /* This is cy */ + /*m[0][2] = 1;*/ + + vx[1] -= m[1][2]*vx[0]; /* This is bx */ + vy[1] -= m[1][2]*vy[0]; /* This is by */ + /* m[1][2] = 0; */ + vx[2] -= m[2][2]*vx[0]; /* This is ax */ + vy[2] -= m[2][2]*vy[0]; /* This is ay */ + /* m[2][2] = 0; */ + + nextcp->x = from->me.x + vx[0]/3; + nextcp->y = from->me.y + vy[0]/3; + prevcp->x = nextcp->x + (vx[0]+vx[1])/3; + prevcp->y = nextcp->y + (vy[0]+vy[1])/3; + } + + test = (nextcp->x-from->me.x)*(to->me.x-from->me.x) + + (nextcp->y-from->me.y)*(to->me.y-from->me.y); + ptest = (prevcp->x-to->me.x)*(from->me.x-to->me.x) + + (prevcp->y-to->me.y)*(from->me.y-to->me.y); + if ( order2 && + (test<nmin || test>nmax || ptest<pmin || ptest>pmax)) + continue; + if ( test>=nmin && test<=nmax ) { + nres.x += nextcp->x; nres.y += nextcp->y; + ++nrescnt; + } + if ( test>=pmin && test<=pmax ) { + pres.x += prevcp->x; pres.y += prevcp->y; + ++prescnt; + } + if ( nrescnt==1 && prescnt==1 ) + break; + } + + ret = 0; + if ( nrescnt>0 ) { + ret |= 1; + nextcp->x = nres.x/nrescnt; + nextcp->y = nres.y/nrescnt; + } else + *nextcp = from->nextcp; + if ( prescnt>0 ) { + ret |= 2; + prevcp->x = pres.x/prescnt; + prevcp->y = pres.y/prescnt; + } else + *prevcp = to->prevcp; + if ( order2 && ret!=3 ) { + nextcp->x = (nextcp->x + prevcp->x)/2; + nextcp->y = (nextcp->y + prevcp->y)/2; + } + if ( order2 ) + *prevcp = *nextcp; +return( ret ); +} + +static void TestForLinear(SplinePoint *from,SplinePoint *to) { + BasePoint off, cpoff, cpoff2; + bigreal len, co, co2; + + /* Did we make a line? */ + off.x = to->me.x-from->me.x; off.y = to->me.y-from->me.y; + len = sqrt(off.x*off.x + off.y*off.y); + if ( len!=0 ) { + off.x /= len; off.y /= len; + cpoff.x = from->nextcp.x-from->me.x; cpoff.y = from->nextcp.y-from->me.y; + len = sqrt(cpoff.x*cpoff.x + cpoff.y*cpoff.y); + if ( len!=0 ) { + cpoff.x /= len; cpoff.y /= len; + } + cpoff2.x = to->prevcp.x-from->me.x; cpoff2.y = to->prevcp.y-from->me.y; + len = sqrt(cpoff2.x*cpoff2.x + cpoff2.y*cpoff2.y); + if ( len!=0 ) { + cpoff2.x /= len; cpoff2.y /= len; + } + co = cpoff.x*off.y - cpoff.y*off.x; co2 = cpoff2.x*off.y - cpoff2.y*off.x; + if ( co<.05 && co>-.05 && co2<.05 && co2>-.05 ) { + from->nextcp = from->me; + to->prevcp = to->me; + } else { + Spline temp; + memset(&temp,0,sizeof(temp)); + temp.from = from; temp.to = to; + SplineRefigure(&temp); + if ( SplineIsLinear(&temp)) { + from->nextcp = from->me; + to->prevcp = to->me; + } + } + } +} + +/* Find a spline which best approximates the list of intermediate points we */ +/* are given. No attempt is made to use fixed slope angles */ +/* given a set of points (x,y,t) */ +/* find the bezier spline which best fits those points */ + +/* OK, we know the end points, so all we really need are the control points */ + /* For cubics.... */ +/* Pf = point from */ +/* CPf = control point, from nextcp */ +/* CPt = control point, to prevcp */ +/* Pt = point to */ +/* S(t) = Pf + 3*(CPf-Pf)*t + 3*(CPt-2*CPf+Pf)*t^2 + (Pt-3*CPt+3*CPf-Pf)*t^3 */ +/* S(t) = Pf - 3*Pf*t + 3*Pf*t^2 - Pf*t^3 + Pt*t^3 + */ +/* 3*(t-2*t^2+t^3)*CPf + */ +/* 3*(t^2-t^3)*CPt */ +/* We want to minimize Σ [S(ti)-Pi]^2 */ +/* There are four variables CPf.x, CPf.y, CPt.x, CPt.y */ +/* When we take the derivative of the error term above with each of these */ +/* variables, we find that the two coordinates are separate. So I shall only */ +/* work through the equations once, leaving off the coordinate */ +/* d error/dCPf = Σ 2*3*(t-2*t^2+t^3) * [S(ti)-Pi] = 0 */ +/* d error/dCPt = Σ 2*3*(t^2-t^3) * [S(ti)-Pi] = 0 */ + /* For quadratics.... */ +/* CP = control point, there's only one */ +/* S(t) = Pf + 2*(CP-Pf)*t + (Pt-2*CP+Pf)*t^2 */ +/* S(t) = Pf - 2*Pf*t + Pf*t^2 + Pt*t^2 + */ +/* 2*(t-2*t^2)*CP */ +/* We want to minimize Σ [S(ti)-Pi]^2 */ +/* There are two variables CP.x, CP.y */ +/* d error/dCP = Σ 2*2*(t-2*t^2) * [S(ti)-Pi] = 0 */ +/* Σ (t-2*t^2) * [Pf - 2*Pf*t + Pf*t^2 + Pt*t^2 - Pi + */ +/* 2*(t-2*t^2)*CP] = 0 */ +/* CP * (Σ 2*(t-2*t^2)*(t-2*t^2)) = Σ (t-2*t^2) * [Pf - 2*Pf*t + Pf*t^2 + Pt*t^2 - Pi] */ + +/* Σ (t-2*t^2) * [Pf - 2*Pf*t + Pf*t^2 + Pt*t^2 - Pi] */ +/* CP = ----------------------------------------------------- */ +/* Σ 2*(t-2*t^2)*(t-2*t^2) */ +Spline *ApproximateSplineFromPoints(SplinePoint *from, SplinePoint *to, + FitPoint *mid, int cnt, int order2) { + int ret; + Spline *spline; + BasePoint nextcp, prevcp; + int i; + + if ( order2 ) { + bigreal xconst, yconst, term /* Same for x and y */; + xconst = yconst = term = 0; + for ( i=0; i<cnt; ++i ) { + bigreal t = mid[i].t, t2 = t*t; + bigreal tfactor = (t-2*t2); + term += 2*tfactor*tfactor; + xconst += tfactor*(from->me.x*(1-2*t+t2) + to->me.x*t2 - mid[i].p.x); + yconst += tfactor*(from->me.y*(1-2*t+t2) + to->me.y*t2 - mid[i].p.y); + } + if ( term!=0 ) { + BasePoint cp; + cp.x = xconst/term; cp.y = yconst/term; + from->nextcp = to->prevcp = cp; +return( SplineMake2(from,to)); + } + } else { + bigreal xconst[2], yconst[2], f_term[2], t_term[2] /* Same for x and y */; + bigreal tfactor[2], determinant; + xconst[0] = xconst[1] = yconst[0] = yconst[1] = + f_term[0] = f_term[1] = t_term[0] = t_term[1] = 0; + for ( i=0; i<cnt; ++i ) { + bigreal t = mid[i].t, t2 = t*t, t3=t*t2; + bigreal xc = (from->me.x*(1-3*t+3*t2-t3) + to->me.x*t3 - mid[i].p.x); + bigreal yc = (from->me.y*(1-3*t+3*t2-t3) + to->me.y*t3 - mid[i].p.y); + tfactor[0] = (t-2*t2+t3); tfactor[1]=(t2-t3); + xconst[0] += tfactor[0]*xc; + xconst[1] += tfactor[1]*xc; + yconst[0] += tfactor[0]*yc; + yconst[1] += tfactor[1]*yc; + f_term[0] += 3*tfactor[0]*tfactor[0]; + f_term[1] += 3*tfactor[0]*tfactor[1]; + /*t_term[0] += 3*tfactor[1]*tfactor[0];*/ + t_term[1] += 3*tfactor[1]*tfactor[1]; + } + t_term[0] = f_term[1]; + determinant = f_term[1]*t_term[0] - f_term[0]*t_term[1]; + if ( determinant!=0 ) { + to->prevcp.x = -(xconst[0]*f_term[1]-xconst[1]*f_term[0])/determinant; + to->prevcp.y = -(yconst[0]*f_term[1]-yconst[1]*f_term[0])/determinant; + if ( f_term[0]!=0 ) { + from->nextcp.x = (-xconst[0]-t_term[0]*to->prevcp.x)/f_term[0]; + from->nextcp.y = (-yconst[0]-t_term[0]*to->prevcp.y)/f_term[0]; + } else { + from->nextcp.x = (-xconst[1]-t_term[1]*to->prevcp.x)/f_term[1]; + from->nextcp.y = (-yconst[1]-t_term[1]*to->prevcp.y)/f_term[1]; + } +return( SplineMake3(from,to)); + } + } + + if ( (spline = IsLinearApprox(from,to,mid,cnt,order2))!=NULL ) +return( spline ); + + ret = _ApproximateSplineFromPoints(from,to,mid,cnt,&nextcp,&prevcp,order2); + + if ( ret&1 ) { + from->nextcp = nextcp; + } else { + from->nextcp = from->me; + } + if ( ret&2 ) { + to->prevcp = prevcp; + } else { + to->prevcp = to->me; + } + TestForLinear(from,to); + spline = SplineMake(from,to,order2); +return( spline ); +} + +static bigreal ClosestSplineSolve(Spline1D *sp,bigreal sought,bigreal close_to_t) { + /* We want to find t so that spline(t) = sought */ + /* find the value which is closest to close_to_t */ + /* on error return closetot */ + extended ts[3]; + int i; + bigreal t, best, test; + + _CubicSolve(sp,sought,ts); + best = 9e20; t= close_to_t; + for ( i=0; i<3; ++i ) if ( ts[i]>-.0001 && ts[i]<1.0001 ) { + if ( (test=ts[i]-close_to_t)<0 ) test = -test; + if ( test<best ) { + best = test; + t = ts[i]; + } + } + +return( t ); +} + +struct dotbounds { + BasePoint unit; + BasePoint base; + bigreal len; + bigreal min,max; /* If min<0 || max>len the spline extends beyond its endpoints */ +}; + +static bigreal SigmaDeltas(Spline *spline, FitPoint *mid, int cnt, DBounds *b, struct dotbounds *db) { + int i; + bigreal xdiff, ydiff, sum, temp, t; + SplinePoint *to = spline->to, *from = spline->from; + extended ts[2], x,y; + struct dotbounds db2; + bigreal dot; + int near_vert, near_horiz; + + if ( (xdiff = to->me.x-from->me.x)<0 ) xdiff = -xdiff; + if ( (ydiff = to->me.y-from->me.y)<0 ) ydiff = -ydiff; + near_vert = ydiff>2*xdiff; + near_horiz = xdiff>2*ydiff; + + sum = 0; + for ( i=0; i<cnt; ++i ) { + if ( near_vert ) { + t = ClosestSplineSolve(&spline->splines[1],mid[i].p.y,mid[i].t); + } else if ( near_horiz ) { + t = ClosestSplineSolve(&spline->splines[0],mid[i].p.x,mid[i].t); + } else { + t = (ClosestSplineSolve(&spline->splines[1],mid[i].p.y,mid[i].t) + + ClosestSplineSolve(&spline->splines[0],mid[i].p.x,mid[i].t))/2; + } + temp = mid[i].p.x - ( ((spline->splines[0].a*t+spline->splines[0].b)*t+spline->splines[0].c)*t + spline->splines[0].d ); + sum += temp*temp; + temp = mid[i].p.y - ( ((spline->splines[1].a*t+spline->splines[1].b)*t+spline->splines[1].c)*t + spline->splines[1].d ); + sum += temp*temp; + } + + /* Ok, we've got distances from a set of points on the old spline to the */ + /* new one. Let's do the reverse: find the extrema of the new and see how*/ + /* close they get to the bounding box of the old */ + /* And get really unhappy if the spline extends beyond the end-points */ + db2.min = 0; db2.max = db->len; + SplineFindExtrema(&spline->splines[0],&ts[0],&ts[1]); + for ( i=0; i<2; ++i ) { + if ( ts[i]!=-1 ) { + x = ((spline->splines[0].a*ts[i]+spline->splines[0].b)*ts[i]+spline->splines[0].c)*ts[i] + spline->splines[0].d; + y = ((spline->splines[1].a*ts[i]+spline->splines[1].b)*ts[i]+spline->splines[1].c)*ts[i] + spline->splines[1].d; + if ( x<b->minx ) + sum += (x-b->minx)*(x-b->minx); + else if ( x>b->maxx ) + sum += (x-b->maxx)*(x-b->maxx); + dot = (x-db->base.x)*db->unit.x + (y-db->base.y)*db->unit.y; + if ( dot<db2.min ) db2.min = dot; + if ( dot>db2.max ) db2.max = dot; + } + } + SplineFindExtrema(&spline->splines[1],&ts[0],&ts[1]); + for ( i=0; i<2; ++i ) { + if ( ts[i]!=-1 ) { + x = ((spline->splines[0].a*ts[i]+spline->splines[0].b)*ts[i]+spline->splines[0].c)*ts[i] + spline->splines[0].d; + y = ((spline->splines[1].a*ts[i]+spline->splines[1].b)*ts[i]+spline->splines[1].c)*ts[i] + spline->splines[1].d; + if ( y<b->miny ) + sum += (y-b->miny)*(y-b->miny); + else if ( y>b->maxy ) + sum += (y-b->maxy)*(y-b->maxy); + dot = (x-db->base.x)*db->unit.x + (y-db->base.y)*db->unit.y; + if ( dot<db2.min ) db2.min = dot; + if ( dot>db2.max ) db2.max = dot; + } + } + + /* Big penalty for going beyond the range of the desired spline */ + if ( db->min==0 && db2.min<0 ) + sum += 10000 + db2.min*db2.min; + else if ( db2.min<db->min ) + sum += 100 + (db2.min-db->min)*(db2.min-db->min); + if ( db->max==db->len && db2.max>db->len ) + sum += 10000 + (db2.max-db->max)*(db2.max-db->max); + else if ( db2.max>db->max ) + sum += 100 + (db2.max-db->max)*(db2.max-db->max); + +return( sum ); +} + +static void ApproxBounds(DBounds *b, FitPoint *mid, int cnt, struct dotbounds *db) { + int i; + bigreal dot; + + b->minx = b->maxx = mid[0].p.x; + b->miny = b->maxy = mid[0].p.y; + db->min = 0; db->max = db->len; + for ( i=1; i<cnt; ++i ) { + if ( mid[i].p.x>b->maxx ) b->maxx = mid[i].p.x; + if ( mid[i].p.x<b->minx ) b->minx = mid[i].p.x; + if ( mid[i].p.y>b->maxy ) b->maxy = mid[i].p.y; + if ( mid[i].p.y<b->miny ) b->miny = mid[i].p.y; + dot = (mid[i].p.x-db->base.x)*db->unit.x + (mid[i].p.y-db->base.y)*db->unit.y; + if ( dot<db->min ) db->min = dot; + if ( dot>db->max ) db->max = dot; + } +} + +/* pf == point from (start point) */ +/* Δf == slope from (cp(from) - from) */ +/* pt == point to (end point, t==1) */ +/* Δt == slope to (cp(to) - to) */ + +/* A spline from pf to pt with slope vectors rf*Δf, rt*Δt is: */ +/* p(t) = pf + [ 3*rf*Δf ]*t + 3*[pt-pf+rt*Δt-2*rf*Δf] *t^2 + */ +/* [2*pf-2*pt+3*rf*Δf-3*rt*Δt]*t^3 */ + +/* So I want */ +/* d Σ (p(t(i))-p(i))^2/ d rf == 0 */ +/* d Σ (p(t(i))-p(i))^2/ d rt == 0 */ +/* now... */ +/* d Σ (p(t(i))-p(i))^2/ d rf == 0 */ +/* => Σ 3*t*Δf*(1-2*t+t^2)* + * [pf-pi+ 3*(pt-pf)*t^2 + 2*(pf-pt)*t^3] + + * 3*[t - 2*t^2 + t^3]*Δf*rf + + * 3*[t^2-t^3]*Δt*rt */ +/* and... */ +/* d Σ (p(t(i))-p(i))^2/ d rt == 0 */ +/* => Σ 3*t^2*Δt*(1-t)* + * [pf-pi+ 3*(pt-pf)*t^2 + 2*(pf-pt)*t^3] + + * 3*[t - 2*t^2 + t^3]*Δf*rf + + * 3*[t^2-t^3]*Δt*rt */ + +/* Now for a long time I looked at that and saw four equations and two unknowns*/ +/* That was I was trying to solve for x and y separately, and that doesn't work. */ +/* There are really just two equations and each sums over both x and y components */ + +/* Old comment: */ +/* I used to do a least squares approach adding two more to the above set of equations */ +/* which held the slopes constant. But that didn't work very well. So instead*/ +/* Then I tried doing the approximation, and then forcing the control points */ +/* to be in line (with the original slopes), getting a better approximation */ +/* to "t" for each data point and then calculating an error array, approximating*/ +/* it, and using that to fix up the final result */ +/* Then I tried checking various possible cp lengths in the desired directions*/ +/* finding the best one or two, and doing a 2D binary search using that as a */ +/* starting point. */ +/* And sometimes a least squares approach will give us the right answer, so */ +/* try that too. */ +/* This still isn't as good as I'd like it... But I haven't been able to */ +/* improve it further yet */ +/* The mergetype mt is either of: */ +/* mt_matrix; original, fast, all-purpose (relies on matrix calculations) */ +/* mt_levien; by Raph Levien (implemented by Linus Romer), fast, accurate, use only if mid is on spline */ +/* mt_bruteforce; slow, all-purpose, normally more accurate than mt_matrix.*/ +/* The mt_levien algorithm is explained here: */ +/* raphlinus.github.io/curves/2021/03/11/bezier-fitting.html */ +/* The notation used here is a bit different: Instead of theta1, theta2, */ +/* delta1, delta2, momentx, area we use alpha,beta,a,b,m,f: */ +/* Here is to complete math that we are using: */ +/* Signed area of the cubic bezier spline a .. controls b and c .. d to the x-axis */ +/* f = ((xb-xa)*(10*ya+6*yb+3*yc+yd)+(xc-xb)*(4*ya+6*yb+6*yc+4*yd)+(xd-xc)*(ya+3*yb+6*yc+10*yd))/20; */ +/* simplified for the normed case */ +/* f = 3/20*(2*a*sin(alpha)+2*b*sin(beta)-a*b*sin(alpha+beta)); */ +/* solved for b */ +/* b = (20*f-6*a*sin(alpha))/(6*sin(beta)-3*a*sin(alpha+beta)). */ +/* Signed area of the cubic bezier spline a .. controls b and c .. d to the x-axis */ +/* from point a up to the bezier point at time t */ +/* f(t) = ((((1-t)*xa+xb*t)-xa)*(10*ya+6*((1-t)*ya+yb*t)+3*((1-t)^2*ya+2*(1-t)*t*yb+t^2*yc) */ +/* +((1-t)^3*ya+3*(1-t)^2*t*yb+3*(1-t)*t^2*yc+t^3*yd))+(((1-t)^2*xa+2*(1-t)*t*xb+t^2*xc) */ +/* -((1-t)*xa+xb*t))*(4*ya+6*((1-t)*ya+yb*t)+6*((1-t)^2*ya+2*(1-t)*t*yb+t^2*yc) */ +/* +4*((1-t)^3*ya+3*(1-t)^2*t*yb+3*(1-t)*t^2*yc+t^3*yd))+(((1-t)^3*xa+3*(1-t)^2*t*xb */ +/* +3*(1-t)*t^2*xc+t^3*xd)-((1-t)^2*xa+2*(1-t)*t*xb+t^2*xc))*(ya+3*((1-t)*ya+yb*t) */ +/* +6*((1-t)^2*ya+2*(1-t)*t*yb+t^2*yc)+10*((1-t)^3*ya+3*(1-t)^2*t*yb+3*(1-t)*t^2*yc+t^3*yd)))/20; */ +/* simplified for the normed case: */ +/* f(t) = -(3*(30*a*b*sin(beta-alpha)*t^6+15*b^2*sin(2*beta)*t^6-20*b*sin(beta)*t^6 */ +/* -15*a^2*sin(2*alpha)*t^6+20*a*sin(alpha)*t^6+6*a*b*sin(beta+alpha)*t^5 */ +/* -90*a*b*sin(beta-alpha)*t^5-30*b^2*sin(2*beta)*t^5+48*b*sin(beta)*t^5 */ +/* +60*a^2*sin(2*alpha)*t^5-72*a*sin(alpha)*t^5-15*a*b*sin(beta+alpha)*t^4 */ +/* +90*a*b*sin(beta-alpha)*t^4+15*b^2*sin(2*beta)*t^4-30*b*sin(beta)*t^4 */ +/* -90*a^2*sin(2*alpha)*t^4+90*a*sin(alpha)*t^4+10*a*b*sin(beta+alpha)*t^3 */ +/* -30*a*b*sin(beta-alpha)*t^3+60*a^2*sin(2*alpha)*t^3-40*a*sin(alpha)*t^3 */ +/* -15*a^2*sin(2*alpha)*t^2))/20. */ +/* First moment about y-axis = \int x dA = \int x dA/dt dt for a cubic bezier */ +/* path a .. controls b and c .. d */ +/* m = (280*xd^2*yd-105*xc*xd*yd-30*xb*xd*yd-5*xa*xd*yd-45*xc^2*yd-45*xb*xc*yd */ +/* -12*xa*xc*yd-18*xb^2*yd-15*xa*xb*yd-5*xa^2*yd+105*xd^2*yc+45*xc*xd*yc */ +/* -3*xa*xd*yc-27*xb*xc*yc-18*xa*xc*yc-27*xb^2*yc-45*xa*xb*yc-30*xa^2*yc */ +/* +30*xd^2*yb+45*xc*xd*yb+18*xb*xd*yb+3*xa*xd*yb+27*xc^2*yb+27*xb*xc*yb */ +/* -45*xa*xb*yb-105*xa^2*yb+5*xd^2*ya+15*xc*xd*ya+12*xb*xd*ya+5*xa*xd*ya */ +/* +18*xc^2*ya+45*xb*xc*ya+30*xa*xc*ya+45*xb^2*ya+105*xa*xb*ya-280*xa^2*ya)/840; */ +/* simplified for the normed case */ +/* m = (9*a*cos(alpha)*b^2*cos(beta)*sin(beta)-15*b^2*cos(beta)*sin(beta) */ +/* -9*a^2*cos(alpha)^2*b*sin(beta)-9*a*cos(alpha)*b*sin(beta)+50*b*sin(beta) */ +/* +9*a*sin(alpha)*b^2*cos(beta)^2-9*a^2*cos(alpha)*sin(alpha)*b*cos(beta) */ +/* -33*a*sin(alpha)*b*cos(beta)+15*a^2*cos(alpha)*sin(alpha)+34*a*sin(alpha))/280; */ +/* normed case combined with the formula for b depending on the area (see above): */ +/* m = (34*a*sin(alpha)+50*(20*f-6*a*sin(alpha))/(6*sin(beta)-3*a*sin(beta+alpha))*sin(beta) */ +/* +15*a^2*sin(alpha)*cos(alpha)-15*(20*f-6*a*sin(alpha))/(6*sin(beta) */ +/* -3*a*sin(beta+alpha))^2*sin(beta)*cos(beta)-a*(20*f-6*a*sin(alpha))/(6*sin(beta) */ +/* -3*a*sin(beta+alpha))*(33*sin(alpha)*cos(beta)+9*cos(alpha)*sin(beta)) */ +/* -9*a^2*(20*f-6*a*sin(alpha))/(6*sin(beta)-3*a*sin(beta+alpha))*sin(alpha+beta)*cos(alpha) */ +/* +9*a*(20*f-6*a*sin(alpha))/(6*sin(beta)-3*a*sin(beta+alpha))^2*sin(alpha+beta)*cos(beta))/280; */ +/* and reduced to a quartic equation with sa = sin(alpha), sb = sin(beta), ca = cos(alpha), cb = cos(beta) */ +/* 0 = -9*ca*(((2*sb*cb*ca+sa*(2*cb*cb-1))*ca-2*sb*cb)*ca-cb*cb*sa) * a^4 */ +/* + 12*((((cb*(30*f*cb-sb)-15*f)*ca+2*sa-cb*sa*(cb+30*f*sb))*ca+cb*(sb-15*f*cb))*ca-sa*cb*cb) * a^3 */ +/* + 12*((((70*m+15*f)*sb^2+cb*(9*sb-70*cb*m-5*cb*f))*ca-5*sa*sb*(3*sb-4*cb*(7*m+f)))*ca-cb*(9*sb-70*cb*m-5*cb*f)) * a^2 */ +/* + 16*(((12*sa-5*ca*(42*m-17*f))*sb-70*cb*(3*m-f)*sa-75*ca*cb*f*f)*sb-75*cb^2*f^2*sa) * a */ +/* + 80*sb*(42*sb*m-25*f*(sb-cb*f)); */ +/* this quartic equation reduces to a quadratic for the special case beta = pi - alpha or beta = -alpha */ +/* 0 = -9*ca*sa^2 * a^3 */ +/* + 6*sa*(4*sa+5*ca*f) * a^2 */ +/* + 10*((42*m-25*f)*sa-25*ca*f^2). */ +/* The derivative of the first moment (not the quartic) = 0 results in a quartic as well: */ +/* 0 = -9*ca*sa*sab^3 * a^4 */ +/* -3*sab^2*(9*ca*sa*sb-(17*sa+30*ca*f)*sab+15*cb*sa^2) * a^3 */ +/* +18*sab*sb*(21*ca*sa*sb-(17*sa+30*ca*f)*sab+15*cb*sa^2) * a^2 */ +/* -4*(144*ca*sa*sb^3+((-78*sa-135*ca*f)*sab+108*cb*sa^2)*sb^2+(-125*f*sab^2-45*cb*f*sa*sab)*sb+150*cb*f^2*sab^2) * a */ +/* +8*sb*((24*sa+45*ca*f)*sb^2+(15*cb*f*sa-125*f*sab)*sb+100*cb*f^2*sab) */ +/* this quartic equation reduces to a linear for the special case beta = pi - alpha or beta = -alpha */ +/* 0 = -3*ca*sa * a */ +/* +4*sa+5*ca*f */ +#define TRY_CNT 2 +#define DECIMATION 5 +Spline *ApproximateSplineFromPointsSlopes(SplinePoint *from, SplinePoint *to, + FitPoint *mid, int cnt, int order2, enum mergetype mt) { + BasePoint tounit, fromunit, ftunit; + bigreal flen,tlen,ftlen,dot; + Spline *spline, temp; + BasePoint nextcp; + int bettern, betterp; + bigreal offn, offp, incrn, incrp, trylen; + int nocnt = 0, totcnt; + bigreal curdiff, bestdiff[TRY_CNT]; + int i,j,besti[TRY_CNT],bestj[TRY_CNT],k,l; + bigreal fdiff, tdiff, fmax, tmax, fdotft, tdotft; + DBounds b; + struct dotbounds db; + bigreal offn_, offp_, finaldiff; + bigreal pt_pf_x, pt_pf_y, determinant; + bigreal consts[2], rt_terms[2], rf_terms[2]; + + /* If all the selected points are at the same spot, and one of the */ + /* end-points is also at that spot, then just copy the control point */ + /* But our caller seems to have done that for us */ + + /* If the two end-points are corner points then allow the slope to vary */ + /* Or if one end-point is a tangent but the point defining the tangent's */ + /* slope is being removed then allow the slope to vary */ + /* Except if the slope is horizontal or vertical then keep it fixed */ + if ( ( !from->nonextcp && ( from->nextcp.x==from->me.x || from->nextcp.y==from->me.y)) || + (!to->noprevcp && ( to->prevcp.x==to->me.x || to->prevcp.y==to->me.y)) ) + /* Preserve the slope */; + else if ( ((from->pointtype == pt_corner && from->nonextcp) || + (from->pointtype == pt_tangent && + ((from->nonextcp && from->noprevcp) || !from->noprevcp))) && + ((to->pointtype == pt_corner && to->noprevcp) || + (to->pointtype == pt_tangent && + ((to->nonextcp && to->noprevcp) || !to->nonextcp))) ) { + from->pointtype = to->pointtype = pt_corner; +return( ApproximateSplineFromPoints(from,to,mid,cnt,order2) ); + } + + /* If we are going to honour the slopes of a quadratic spline, there is */ + /* only one possibility */ + if ( order2 ) { + if ( from->nonextcp ) + from->nextcp = from->next->to->me; + if ( to->noprevcp ) + to->prevcp = to->prev->from->me; + fromunit.x = from->nextcp.x-from->me.x; fromunit.y = from->nextcp.y-from->me.y; + tounit.x = to->prevcp.x-to->me.x; tounit.y = to->prevcp.y-to->me.y; + + if ( !IntersectLines(&nextcp,&from->nextcp,&from->me,&to->prevcp,&to->me) || + (nextcp.x-from->me.x)*fromunit.x + (nextcp.y-from->me.y)*fromunit.y < 0 || + (nextcp.x-to->me.x)*tounit.x + (nextcp.y-to->me.y)*tounit.y < 0 ) { + /* If the slopes don't intersect then use a line */ + /* (or if the intersection is patently absurd) */ + from->nextcp = from->me; + to->prevcp = to->me; + TestForLinear(from,to); + } else { + from->nextcp = to->prevcp = nextcp; + } +return( SplineMake2(from,to)); + } + /* From here down we are only working with cubic splines */ + + if ( cnt==0 ) { + /* Just use what we've got, no data to improve it */ + /* But we do sometimes get some cps which are too big */ + bigreal len = sqrt((to->me.x-from->me.x)*(to->me.x-from->me.x) + (to->me.y-from->me.y)*(to->me.y-from->me.y)); + if ( len==0 ) { + from->nextcp = from->me; + to->prevcp = to->me; + } else { + BasePoint noff, poff; + bigreal nlen, plen; + noff.x = from->nextcp.x-from->me.x; noff.y = from->nextcp.y-from->me.y; + poff.x = to->me.x-to->prevcp.x; poff.y = to->me.y-to->prevcp.y; + nlen = sqrt(noff.x*noff.x + noff.y+noff.y); + plen = sqrt(poff.x*poff.x + poff.y+poff.y); + if ( nlen>len/3 ) { + noff.x *= len/3/nlen; noff.y *= len/3/nlen; + from->nextcp.x = from->me.x + noff.x; + from->nextcp.y = from->me.y + noff.y; + } + if ( plen>len/3 ) { + poff.x *= len/3/plen; poff.y *= len/3/plen; + to->prevcp.x = to->me.x + poff.x; + to->prevcp.y = to->me.y + poff.y; + } + } +return( SplineMake3(from,to)); + } + + if ( to->prev!=NULL && (( to->noprevcp && to->nonextcp ) || to->prev->knownlinear )) { + tounit.x = to->prev->from->me.x-to->me.x; tounit.y = to->prev->from->me.y-to->me.y; +// g_message("1 tu: %f %f, %d %d %d", tounit.x, tounit.y, to->noprevcp, to->nonextcp, to->prev->knownlinear); + } else if ( !to->noprevcp || to->pointtype == pt_corner ) { + tounit.x = to->prevcp.x-to->me.x; tounit.y = to->prevcp.y-to->me.y; +// g_message("2 tu: %f %f", tounit.x, tounit.y); + } else { + tounit.x = to->me.x-to->nextcp.x; tounit.y = to->me.y-to->nextcp.y; +// g_message("3 tu: %f %f", tounit.x, tounit.y); + } + tlen = sqrt(tounit.x*tounit.x + tounit.y*tounit.y); + if ( from->next!=NULL && (( from->noprevcp && from->nonextcp ) || from->next->knownlinear) ) { + fromunit.x = from->next->to->me.x-from->me.x; fromunit.y = from->next->to->me.y-from->me.y; +// g_message("1 fu"); + } else if ( !from->nonextcp || from->pointtype == pt_corner ) { + fromunit.x = from->nextcp.x-from->me.x; fromunit.y = from->nextcp.y-from->me.y; +// g_message("2 fu"); + } else { + fromunit.x = from->me.x-from->prevcp.x; fromunit.y = from->me.y-from->prevcp.y; +// g_message("3 fu"); + } + flen = sqrt(fromunit.x*fromunit.x + fromunit.y*fromunit.y); +// g_message("tu: %f %f, fu: %f, %f", tounit.x, tounit.y, fromunit.x, fromunit.y); + if ( (tlen==0 || flen==0) && (from->next==NULL || to->prev==NULL) ) { + memset(&temp,0,sizeof(temp)); + temp.from = from; temp.to = to; + SplineRefigure(&temp); + from->next = to->prev = NULL; + } + if ( tlen==0 ) { + if ( to->prev!=NULL ) { + temp = *to->prev; + } + if ( (to->pointtype==pt_curve || to->pointtype==pt_hvcurve) && + to->next && !to->nonextcp ) { + tounit.x = to->me.x-to->nextcp.x; tounit.y = to->me.y-to->nextcp.y; + } else { + tounit.x = -( (3*temp.splines[0].a*.9999+2*temp.splines[0].b)*.9999+temp.splines[0].c ); + tounit.y = -( (3*temp.splines[1].a*.9999+2*temp.splines[1].b)*.9999+temp.splines[1].c ); + } + tlen = sqrt(tounit.x*tounit.x + tounit.y*tounit.y); + } + tounit.x /= tlen; tounit.y /= tlen; + + if ( flen==0 ) { + if ( from->next!=NULL ) { + temp = *from->next; + } + if ( (from->pointtype==pt_curve || from->pointtype==pt_hvcurve) && + from->prev && !from->noprevcp ) { + fromunit.x = from->me.x-from->prevcp.x; fromunit.y = from->me.y-from->prevcp.y; + } else { + fromunit.x = ( (3*temp.splines[0].a*.0001+2*temp.splines[0].b)*.0001+temp.splines[0].c ); + fromunit.y = ( (3*temp.splines[1].a*.0001+2*temp.splines[1].b)*.0001+temp.splines[1].c ); + } + flen = sqrt(fromunit.x*fromunit.x + fromunit.y*fromunit.y); + } + fromunit.x /= flen; fromunit.y /= flen; + + ftunit.x = (to->me.x-from->me.x); ftunit.y = (to->me.y-from->me.y); + ftlen = sqrt(ftunit.x*ftunit.x + ftunit.y*ftunit.y); + if ( ftlen!=0 ) { + ftunit.x /= ftlen; ftunit.y /= ftlen; + } + + if ( (dot=fromunit.x*tounit.y - fromunit.y*tounit.x)<.0001 && dot>-.0001 && + (dot=ftunit.x*tounit.y - ftunit.y*tounit.x)<.0001 && dot>-.0001 ) { + /* It's a line. Slopes are parallel, and parallel to vector between (from,to) */ + from->nextcp = from->me; to->prevcp = to->me; +return( SplineMake3(from,to)); + } + /* This is the generic case, where a generic part is approximated by a cubic */ + /* bezier spline. */ + if ( ( ftlen == 0 ) && ( mt != mt_matrix ) ) + mt = mt_matrix; + if ( mt == mt_levien ) { + bigreal f,m,xa,ya,xb,yb,xc,yc,xd,yd,sasa,sab; + int numberOfSolutions; + SplinePoint *frompoint,*topoint; + f = 0; /* area */ + m = 0; /* first area moment about y (along x) */ + frompoint = from; + if ( from->next==NULL ) + topoint=to; + else + topoint=from->next->to; + for ( ; ; frompoint = topoint->next->from, topoint = topoint->next->to ) { + /* normalizing transformation (chord to x axis and length 1) */ + xa = ((frompoint->me.x-from->me.x)*ftunit.x+(frompoint->me.y-from->me.y)*ftunit.y)/ftlen; + ya = (-(frompoint->me.x-from->me.x)*ftunit.y+(frompoint->me.y-from->me.y)*ftunit.x)/ftlen; + xb = ((frompoint->nextcp.x-from->me.x)*ftunit.x+(frompoint->nextcp.y-from->me.y)*ftunit.y)/ftlen; + yb = (-(frompoint->nextcp.x-from->me.x)*ftunit.y+(frompoint->nextcp.y-from->me.y)*ftunit.x)/ftlen; + xc = ((topoint->prevcp.x-from->me.x)*ftunit.x+(topoint->prevcp.y-from->me.y)*ftunit.y)/ftlen; + yc = (-(topoint->prevcp.x-from->me.x)*ftunit.y+(topoint->prevcp.y-from->me.y)*ftunit.x)/ftlen; + xd = ((topoint->me.x-from->me.x)*ftunit.x+(topoint->me.y-from->me.y)*ftunit.y)/ftlen; + yd = (-(topoint->me.x-from->me.x)*ftunit.y+(topoint->me.y-from->me.y)*ftunit.x)/ftlen; + f += ((xb-xa)*(10*ya+6*yb+3*yc+yd)+(xc-xb)*(4*ya+6*yb+6*yc+4*yd)+(xd-xc)*(ya+3*yb+6*yc+10*yd))/20; + m += (280*xd*xd*yd-105*xc*xd*yd-30*xb*xd*yd-5*xa*xd*yd-45*xc*xc*yd-45*xb*xc*yd-12*xa*xc*yd-18*xb*xb*yd + -15*xa*xb*yd-5*xa*xa*yd+105*xd*xd*yc+45*xc*xd*yc-3*xa*xd*yc-27*xb*xc*yc-18*xa*xc*yc-27*xb*xb*yc + -45*xa*xb*yc-30*xa*xa*yc+30*xd*xd*yb+45*xc*xd*yb+18*xb*xd*yb+3*xa*xd*yb+27*xc*xc*yb+27*xb*xc*yb + -45*xa*xb*yb-105*xa*xa*yb+5*xd*xd*ya+15*xc*xd*ya+12*xb*xd*ya+5*xa*xd*ya+18*xc*xc*ya+45*xb*xc*ya + +30*xa*xc*ya+45*xb*xb*ya+105*xa*xb*ya-280*xa*xa*ya)/840; + if ( topoint==to ) + break; + } + BasePoint aunit = (BasePoint) { BPDot(ftunit, fromunit), BPCross(ftunit, fromunit) }; /* normed direction at "from" */ + BasePoint bunit = (BasePoint) { BPDot(BPRev(ftunit), tounit),BPCross(ftunit, tounit) }; /* normed direction at "to" */ + if ( aunit.y < 0 ) { /* normalize aunit.y to >= 0: */ + aunit.y = -aunit.y; + bunit.y = -bunit.y; + m = -m; + f = -f; + } + /* calculate the Tunni point (where the tangents at "from" and "to" intersect) */ + bigreal aMax = 100; /* maximal value that the handle a can reach up to the Tunni point, 100 is really long */ + bigreal bMax = 100; /* maximal value that the handle b can reach up to the Tunni point, 100 is really long */ + sab = aunit.y*bunit.x+aunit.x*bunit.y; + if (sab != 0) { /* if handles not parallel */ + aMax = bunit.y/sab; + bMax = aunit.y/sab; + if ( aMax < 0 ) { + aMax = 100; + } + if ( bMax < 0 ) { + bMax = 100; + } + } + /* start approximation by solving the quartic equation */ + sasa = aunit.y*aunit.y; /* reducing the multiplications */ + Quartic quad; + if ( (aunit.x == -bunit.x && aunit.y == bunit.y) || (aunit.x == bunit.x && aunit.y == -bunit.y) ) { /* handles are parallel */ + quad.a = 0; + quad.b = 0; + quad.c = -9*aunit.x*sasa; + quad.d = 6*aunit.y*(4*aunit.y+5*aunit.x*f); + quad.e = 10*((42*m-25*f)*aunit.y-25*aunit.x*f*f); + } else { /* generic situation */ + quad.a = -9*aunit.x*(((2*bunit.y*bunit.x*aunit.x+aunit.y + *(2*bunit.x*bunit.x-1))*aunit.x-2*bunit.y*bunit.x) + *aunit.x-bunit.x*bunit.x*aunit.y); + quad.b = 12*((((bunit.x*(30*f*bunit.x-bunit.y)-15*f) + *aunit.x+2*aunit.y-bunit.x*aunit.y*(bunit.x+30*f*bunit.y)) + *aunit.x+bunit.x*(bunit.y-15*f*bunit.x)) + *aunit.x-aunit.y*bunit.x*bunit.x); + quad.c = 12*((((70*m+15*f)*bunit.y*bunit.y+bunit.x + *(9*bunit.y-70*bunit.x*m-5*bunit.x*f)) + *aunit.x-5*aunit.y*bunit.y*(3*bunit.y-4*bunit.x + *(7*m+f)))*aunit.x-bunit.x*(9*bunit.y-70*bunit.x*m-5*bunit.x*f)); + quad.d = 16*(((12*aunit.y-5*aunit.x*(42*m-17*f))*bunit.y + -70*bunit.x*(3*m-f)*aunit.y-75*aunit.x*bunit.x*f*f) + *bunit.y-75*bunit.x*bunit.x*f*f*aunit.y); + quad.e = 80*bunit.y*(42*bunit.y*m-25*f*(bunit.y-bunit.x*f)); + } + extended solutions[4] = {-999999,-999999,-999999,-999999}; + _QuarticSolve(&quad,solutions); + extended abSolutions[10][2]; /* there are at most 4+4+1+1=10 solutions of pairs of a and b (quartic=0,derivative=0,b=0.01,a=0.01) */ + numberOfSolutions = 0; + extended a,b; + for( int i = 0; i < 4; i++ ){ + a = solutions[i]; + if ( a >= 0 && a < aMax ) { + b = (20*f-6*a*aunit.y)/(3*(2*bunit.y-a*sab)); + if ( b >= 0 && b < bMax ) { + abSolutions[numberOfSolutions][0] = a; + abSolutions[numberOfSolutions++][1] = b; + } + } + } + /* and now again for the derivative (of m not of the upper quartic): */ + if ( (aunit.x == -bunit.x && aunit.y == bunit.y) || (aunit.x == bunit.x && aunit.y == -bunit.y) ) { /* handles are parallel */ + quad.a = 0; + quad.b = 0; + quad.c = 0; + quad.d = -3*aunit.x*aunit.y; + quad.e = 4*aunit.y+5*aunit.x*f; + } else { /* generic situation */ + bigreal sbsb = bunit.y*bunit.y; + bigreal sabsab = sab*sab; + quad.a = -9*aunit.x*aunit.y*sabsab*sab; + quad.b = -3*sabsab*(9*aunit.x*aunit.y*bunit.y-(17*aunit.y + +30*aunit.x*f)*sab+15*bunit.x*sasa); + quad.c = 18*sab*bunit.y*(21*aunit.x*aunit.y*bunit.y-(17*aunit.y + +30*aunit.x*f)*sab+15*bunit.x*sasa); + quad.d = -4*(144*aunit.x*aunit.y*sbsb*bunit.y+((-78*aunit.y- + 135*aunit.x*f)*sab+108*bunit.x*sasa)*sbsb+(-125*f*sabsab + -45*bunit.x*f*aunit.y*sab)*bunit.y+150*bunit.x*f*f*sabsab); + quad.e = 8*bunit.y*((24*aunit.y+45*aunit.x*f)*sbsb + +(15*bunit.x*f*aunit.y-125*f*sab)*bunit.y+100*bunit.x*f*f*sab); + } + for( int i = 0; i < 4; i++ ) /* overwriting (reusing) */ + solutions[i] = -999999; + _QuarticSolve(&quad,solutions); + for( int i = 0; i < 4; i++ ){ + a = solutions[i]; + if ( a >= 0 && a < aMax ) { + b = (20*f-6*a*aunit.y)/(3*(2*bunit.y-a*sab)); + if ( b >= 0 && b < bMax ) { + abSolutions[numberOfSolutions][0] = a; + abSolutions[numberOfSolutions++][1] = b; + } + } + } + /* Add the solution of b = 0.01 (approximately 0 but above because of direction). */ +#if 0 /* those solutions lead to subpar approximations */ + /* This solution is not part of the original algorithm by Raph Levien. */ + a = (2000*f-6*bunit.y)/(600*aunit.y-3*sab); + if ( a >= 0 && a < aMax ) { + abSolutions[numberOfSolutions][0] = a; + abSolutions[numberOfSolutions++][1] = 0.01; + } + /* Add the solution of a = 0.01 (approximately 0 but above because of direction). */ + /* This solution is not part of the original algorithm by Raph Levien. */ + b = (2000*f-6*aunit.y)/(600*bunit.y-3*sab); + if ( b >= 0 && b < bMax ) { + abSolutions[numberOfSolutions][0] = 0.01; + abSolutions[numberOfSolutions++][1] = b; + } + if ( numberOfSolutions == 0) { /* add solutions that extend up to the Tunni point */ + /* try solution with a = aMax and b area-equal*/ + b = (20*f-6*aMax*aunit.y)/(3*(2*bunit.y-aMax*sab)); + if ( b >= 0 && b < bMax ) { + abSolutions[numberOfSolutions][0] = aMax; + abSolutions[numberOfSolutions++][1] = b; + } + /* try solution with b = bMax and a area-equal*/ + a = (20*f-6*bMax*bunit.y)/(3*(2*aunit.y-bMax*sab)); + if ( a >= 0 && a < aMax ) { + abSolutions[numberOfSolutions][0] = a; + abSolutions[numberOfSolutions++][1] = bMax; + } + } +#endif + if ( numberOfSolutions == 0) { + /* no solutions found, quit */ +return NULL; + + /* solution with a = aMax and b = bMax*/ + abSolutions[numberOfSolutions][0] = aMax; + abSolutions[numberOfSolutions++][1] = bMax; + } + if ( numberOfSolutions == 1) { + from->nextcp.x = from->me.x+ftlen*fromunit.x*abSolutions[0][0]; + from->nextcp.y = from->me.y+ftlen*fromunit.y*abSolutions[0][0]; + to->prevcp.x = to->me.x+ftlen*tounit.x*abSolutions[0][1]; + to->prevcp.y = to->me.y+ftlen*tounit.y*abSolutions[0][1]; + } else { /* compare L2 errors to choose the best solution */ + bigreal bestError = 1e30; + bigreal t,error,errorsum,dist; + BasePoint prevcp,nextcp,coeff1,coeff2,coeff3; + int last_best_j; + for (int k=0; k<numberOfSolutions; k++) { + nextcp.x = from->me.x+ftlen*fromunit.x*abSolutions[k][0]; + nextcp.y = from->me.y+ftlen*fromunit.y*abSolutions[k][0]; + prevcp.x = to->me.x+ftlen*tounit.x*abSolutions[k][1]; + prevcp.y = to->me.y+ftlen*tounit.y*abSolutions[k][1]; + /* Calculate the error of the cubic bezier path from,nextcp,prevcp,to: */ + /* In order to do that we calculate 99 points on the bezier path. */ + coeff3.x = -from->me.x+3*nextcp.x-3*prevcp.x+to->me.x; + coeff3.y = -from->me.y+3*nextcp.y-3*prevcp.y+to->me.y; + coeff2.x = 3*from->me.x-6*nextcp.x+3*prevcp.x; + coeff2.y = 3*from->me.y-6*nextcp.y+3*prevcp.y; + coeff1.x = -3*from->me.x+3*nextcp.x; + coeff1.y = -3*from->me.y+3*nextcp.y; + BasePoint approx[99]; + for (int i=0; i<99; i++) { + t = (i+1)/100.0; + approx[i].x = from->me.x+t*(coeff1.x+t*(coeff2.x+t*coeff3.x)); + approx[i].y = from->me.y+t*(coeff1.y+t*(coeff2.y+t*coeff3.y)); + } + /* Now we calculate the error by determining the minimal quadratic distance to the mid points. */ + errorsum = 0.0; + last_best_j = 0; + for (int i=0; i<cnt; i++) { /* Going through the mid points */ + error = 1e30; + /* For this mid point, find the distance to the closest one of the */ + /* 99 points on the approximate cubic bezier. */ + /* To not favour approximations which trace the original multiple times */ + /* by going back and forth, only consider monotonic mappings. */ + /* I.e., start from the point that was closest to the previous mid point: */ + for (int j=last_best_j; j<99; j++) { + dist = (mid[i].p.x-approx[j].x)*(mid[i].p.x-approx[j].x) + +(mid[i].p.y-approx[j].y)*(mid[i].p.y-approx[j].y); + if (dist < error) { + error = dist; + last_best_j = j; + } + } + errorsum += error; + if (errorsum > bestError) + break; + } + if (errorsum < bestError) { + bestError = errorsum; + from->nextcp = nextcp; + to->prevcp = prevcp; + } + } + } + return( SplineMake3(from,to)); + } else if ( mt == mt_bruteforce ) { + bigreal best_error = 1e30; + bigreal t,error,errorsum,dist; + BasePoint prevcp,coeff1,coeff2,coeff3; + bigreal best_fromhandle = 0.0; + bigreal best_tohandle = 0.0; + BasePoint approx[99]; /* The 99 points on the approximate cubic bezier */ + /* We make 2 runs: The first run to narrow the variation range, the second run to finetune */ + /* The optimal length of the two handles are determined by brute force. */ + for (int run=0; run<2; ++run) { + for (int fromhandle=((run==0)?1:-29); fromhandle<=((run==0)?60:29); ++fromhandle) { + for (int tohandle=((run==0)?1:-29); tohandle<=((run==0)?60:29); ++tohandle) { + nextcp.x = from->me.x+ftlen*fromunit.x*( (run==0)?fromhandle:best_fromhandle+fromhandle/30.0 )/60.0; + nextcp.y = from->me.y+ftlen*fromunit.y*( (run==0)?fromhandle:best_fromhandle+fromhandle/30.0 )/60.0; + prevcp.x = to->me.x+ftlen*tounit.x*( (run==0)?tohandle:best_tohandle+tohandle/30.0 )/60.0; + prevcp.y = to->me.y+ftlen*tounit.y*( (run==0)?tohandle:best_tohandle+tohandle/30.0 )/60.0; + /* Calculate the error of the cubic bezier path from,nextcp,prevcp,to: */ + /* In order to do that we calculate 99 points on the bezier path. */ + coeff3.x = -from->me.x+3*nextcp.x-3*prevcp.x+to->me.x; + coeff3.y = -from->me.y+3*nextcp.y-3*prevcp.y+to->me.y; + coeff2.x = 3*from->me.x-6*nextcp.x+3*prevcp.x; + coeff2.y = 3*from->me.y-6*nextcp.y+3*prevcp.y; + coeff1.x = -3*from->me.x+3*nextcp.x; + coeff1.y = -3*from->me.y+3*nextcp.y; + for (int i=0; i<99; ++i) { + t = (i+1)/100.0; + approx[i].x = from->me.x+t*(coeff1.x+t*(coeff2.x+t*coeff3.x)); + approx[i].y = from->me.y+t*(coeff1.y+t*(coeff2.y+t*coeff3.y)); + } + /* Now we calculate the error by determining the minimal quadratic distance to the mid points. */ + errorsum = 0.0; + for (int i=0; i<cnt; ++i) { /* Going through the mid points */ + error = (mid[i].p.x-approx[0].x)*(mid[i].p.x-approx[0].x) + +(mid[i].p.y-approx[0].y)*(mid[i].p.y-approx[0].y); + /* Above we have just initialized the error and */ + /* now we are going through the remaining 98 of */ + /* 99 points on the approximate cubic bezier: */ + for (int j=1; j<99; ++j) { + dist = (mid[i].p.x-approx[j].x)*(mid[i].p.x-approx[j].x) + +(mid[i].p.y-approx[j].y)*(mid[i].p.y-approx[j].y); + if (dist < error) + error = dist; + } + errorsum += error; + if (errorsum > best_error) + break; + } + if (errorsum < best_error) { + best_error = errorsum; + if (run == 0) { + best_fromhandle = fromhandle; + best_tohandle = tohandle; + } + from->nextcp = nextcp; + to->prevcp = prevcp; + } + } + } + } + return( SplineMake3(from,to)); + } + else { /* mergetype mt_matrix (original algorithm) */ + pt_pf_x = to->me.x - from->me.x; + pt_pf_y = to->me.y - from->me.y; + consts[0] = consts[1] = rt_terms[0] = rt_terms[1] = rf_terms[0] = rf_terms[1] = 0; + for ( i=0; i<cnt; ++i ) { + bigreal t = mid[i].t, t2 = t*t, t3=t2*t; + bigreal factor_from = t-2*t2+t3; + bigreal factor_to = t2-t3; + bigreal const_x = from->me.x-mid[i].p.x + 3*pt_pf_x*t2 - 2*pt_pf_x*t3; + bigreal const_y = from->me.y-mid[i].p.y + 3*pt_pf_y*t2 - 2*pt_pf_y*t3; + bigreal temp1 = 3*(t-2*t2+t3); + bigreal rf_term_x = temp1*fromunit.x; + bigreal rf_term_y = temp1*fromunit.y; + bigreal temp2 = 3*(t2-t3); + bigreal rt_term_x = -temp2*tounit.x; + bigreal rt_term_y = -temp2*tounit.y; + + consts[0] += factor_from*( fromunit.x*const_x + fromunit.y*const_y ); + consts[1] += factor_to *( -tounit.x*const_x + -tounit.y*const_y); + rf_terms[0] += factor_from*( fromunit.x*rf_term_x + fromunit.y*rf_term_y); + rf_terms[1] += factor_to*( -tounit.x*rf_term_x + -tounit.y*rf_term_y); + rt_terms[0] += factor_from*( fromunit.x*rt_term_x + fromunit.y*rt_term_y); + rt_terms[1] += factor_to*( -tounit.x*rt_term_x + -tounit.y*rt_term_y); + } + + /* I've only seen singular matrices (determinant==0) when cnt==1 */ + /* but even with cnt==1 the determinant is usually non-0 (16 times out of 17)*/ + determinant = (rt_terms[0]*rf_terms[1]-rt_terms[1]*rf_terms[0]); + if ( determinant!=0 ) { + bigreal rt, rf; + rt = (consts[1]*rf_terms[0]-consts[0]*rf_terms[1])/determinant; + if ( rf_terms[0]!=0 ) + rf = -(consts[0]+rt*rt_terms[0])/rf_terms[0]; + else /* if ( rf_terms[1]!=0 ) This can't happen, otherwise the determinant would be 0 */ + rf = -(consts[1]+rt*rt_terms[1])/rf_terms[1]; + /* If we get bad values (ones that point diametrically opposed to what*/ + /* we need), then fix that factor at 0, and see what we get for the */ + /* other */ + if ( rf>=0 && rt>0 && rf_terms[0]!=0 && + (rf = -consts[0]/rf_terms[0])>0 ) { + rt = 0; + } else if ( rf<0 && rt<=0 && rt_terms[1]!=0 && + (rt = -consts[1]/rt_terms[1])<0 ) { + rf = 0; + } + if ( rt<=0 && rf>=0 ) { + from->nextcp.x = from->me.x + rf*fromunit.x; + from->nextcp.y = from->me.y + rf*fromunit.y; + to->prevcp.x = to->me.x - rt*tounit.x; + to->prevcp.y = to->me.y - rt*tounit.y; +return( SplineMake3(from,to)); + } + } + + trylen = (to->me.x-from->me.x)*fromunit.x + (to->me.y-from->me.y)*fromunit.y; + if ( trylen>flen ) flen = trylen; + + trylen = (from->me.x-to->me.x)*tounit.x + (from->me.y-to->me.y)*tounit.y; + if ( trylen>tlen ) tlen = trylen; + + for ( i=0; i<cnt; ++i ) { + trylen = (mid[i].p.x-from->me.x)*fromunit.x + (mid[i].p.y-from->me.y)*fromunit.y; + if ( trylen>flen ) flen = trylen; + trylen = (mid[i].p.x-to->me.x)*tounit.x + (mid[i].p.y-to->me.y)*tounit.y; + if ( trylen>tlen ) tlen = trylen; + } + + fdotft = fromunit.x*ftunit.x + fromunit.y*ftunit.y; + fmax = fdotft>0 ? ftlen/fdotft : 1e10; + tdotft = -tounit.x*ftunit.x - tounit.y*ftunit.y; + tmax = tdotft>0 ? ftlen/tdotft : 1e10; + /* At fmax, tmax the control points will stretch beyond the other endpoint*/ + /* when projected along the line between the two endpoints */ + + db.base = from->me; + db.unit = ftunit; + db.len = ftlen; + ApproxBounds(&b,mid,cnt,&db); + + for ( k=0; k<TRY_CNT; ++k ) { + bestdiff[k] = 1e20; + besti[k] = -1; bestj[k] = -1; + } + fdiff = flen/DECIMATION; + tdiff = tlen/DECIMATION; + from->nextcp = from->me; + memset(&temp,0,sizeof(Spline)); + temp.from = from; temp.to = to; + for ( i=1; i<DECIMATION; ++i ) { + from->nextcp.x += fdiff*fromunit.x; from->nextcp.y += fdiff*fromunit.y; + to->prevcp = to->me; + for ( j=1; j<DECIMATION; ++j ) { + to->prevcp.x += tdiff*tounit.x; to->prevcp.y += tdiff*tounit.y; + SplineRefigure(&temp); + curdiff = SigmaDeltas(&temp,mid,cnt,&b,&db); + for ( k=0; k<TRY_CNT; ++k ) { + if ( curdiff<bestdiff[k] ) { + for ( l=TRY_CNT-1; l>k; --l ) { + bestdiff[l] = bestdiff[l-1]; + besti[l] = besti[l-1]; + bestj[l] = bestj[l-1]; + } + bestdiff[k] = curdiff; + besti[k] = i; bestj[k]=j; + break; + } + } + } + } + + finaldiff = 1e20; + offn_ = offp_ = -1; + spline = SplineMake(from,to,false); + for ( k=-1; k<TRY_CNT; ++k ) { + if ( k<0 ) { + BasePoint nextcp, prevcp; + bigreal temp1, temp2; + int ret = _ApproximateSplineFromPoints(from,to,mid,cnt,&nextcp,&prevcp,false); + /* sometimes least squares gives us the right answer */ + if ( !(ret&1) || !(ret&2)) + continue; + temp1 = (prevcp.x-to->me.x)*tounit.x + (prevcp.y-to->me.y)*tounit.y; + temp2 = (nextcp.x-from->me.x)*fromunit.x + (nextcp.y-from->me.y)*fromunit.y; + if ( temp1<=0 || temp2<=0 ) /* A nice solution... but the control points are diametrically opposed to what they should be */ + continue; + tlen = temp1; flen = temp2; + } else { + if ( bestj[k]<0 || besti[k]<0 ) + continue; + tlen = bestj[k]*tdiff; flen = besti[k]*fdiff; + } + to->prevcp.x = to->me.x + tlen*tounit.x; to->prevcp.y = to->me.y + tlen*tounit.y; + from->nextcp.x = from->me.x + flen*fromunit.x; from->nextcp.y = from->me.y + flen*fromunit.y; + SplineRefigure(spline); + + bettern = betterp = false; + incrn = tdiff/2.0; incrp = fdiff/2.0; + offn = flen; offp = tlen; + nocnt = 0; + curdiff = SigmaDeltas(spline,mid,cnt,&b,&db); + totcnt = 0; + for (;;) { + bigreal fadiff, fsdiff; + bigreal tadiff, tsdiff; + + from->nextcp.x = from->me.x + (offn+incrn)*fromunit.x; from->nextcp.y = from->me.y + (offn+incrn)*fromunit.y; + to->prevcp.x = to->me.x + offp*tounit.x; to->prevcp.y = to->me.y + offp*tounit.y; + SplineRefigure(spline); + fadiff = SigmaDeltas(spline,mid,cnt,&b,&db); + from->nextcp.x = from->me.x + (offn-incrn)*fromunit.x; from->nextcp.y = from->me.y + (offn-incrn)*fromunit.y; + SplineRefigure(spline); + fsdiff = SigmaDeltas(spline,mid,cnt,&b,&db); + from->nextcp.x = from->me.x + offn*fromunit.x; from->nextcp.y = from->me.y + offn*fromunit.y; + if ( offn-incrn<=0 ) + fsdiff += 1e10; + + to->prevcp.x = to->me.x + (offp+incrp)*tounit.x; to->prevcp.y = to->me.y + (offp+incrp)*tounit.y; + SplineRefigure(spline); + tadiff = SigmaDeltas(spline,mid,cnt,&b,&db); + to->prevcp.x = to->me.x + (offp-incrp)*tounit.x; to->prevcp.y = to->me.y + (offp-incrp)*tounit.y; + SplineRefigure(spline); + tsdiff = SigmaDeltas(spline,mid,cnt,&b,&db); + to->prevcp.x = to->me.x + offp*tounit.x; to->prevcp.y = to->me.y + offp*tounit.y; + if ( offp-incrp<=0 ) + tsdiff += 1e10; + + if ( offn>=incrn && fsdiff<curdiff && + (fsdiff<fadiff && fsdiff<tsdiff && fsdiff<tadiff)) { + offn -= incrn; + if ( bettern>0 ) + incrn /= 2; + bettern = -1; + nocnt = 0; + curdiff = fsdiff; + } else if ( offn+incrn<fmax && fadiff<curdiff && + (fadiff<=fsdiff && fadiff<tsdiff && fadiff<tadiff)) { + offn += incrn; + if ( bettern<0 ) + incrn /= 2; + bettern = 1; + nocnt = 0; + curdiff = fadiff; + } else if ( offp>=incrp && tsdiff<curdiff && + (tsdiff<=fsdiff && tsdiff<=fadiff && tsdiff<tadiff)) { + offp -= incrp; + if ( betterp>0 ) + incrp /= 2; + betterp = -1; + nocnt = 0; + curdiff = tsdiff; + } else if ( offp+incrp<tmax && tadiff<curdiff && + (tadiff<=fsdiff && tadiff<=fadiff && tadiff<=tsdiff)) { + offp += incrp; + if ( betterp<0 ) + incrp /= 2; + betterp = 1; + nocnt = 0; + curdiff = tadiff; + } else { + if ( ++nocnt > 6 ) + break; + incrn /= 2; + incrp /= 2; + } + if ( curdiff<1 ) + break; + if ( incrp<tdiff/1024 || incrn<fdiff/1024 ) + break; + if ( ++totcnt>200 ) + break; + if ( offn<0 || offp<0 ) { + IError("Approximation got inverse control points"); + break; + } + } + if ( curdiff<finaldiff ) { + finaldiff = curdiff; + offn_ = offn; + offp_ = offp; + } + } + + to->noprevcp = offp_==0; + from->nonextcp = offn_==0; + to->prevcp.x = to->me.x + offp_*tounit.x; to->prevcp.y = to->me.y + offp_*tounit.y; + from->nextcp.x = from->me.x + offn_*fromunit.x; from->nextcp.y = from->me.y + offn_*fromunit.y; + /* I used to check for a spline very close to linear (and if so, make it */ + /* linear). But in when stroking a path with an elliptical pen we transform*/ + /* the coordinate system and our normal definitions of "close to linear" */ + /* don't apply */ + /*TestForLinear(from,to);*/ + SplineRefigure(spline); + +return( spline ); + } +} +#undef TRY_CNT +#undef DECIMATION + +SplinePoint *_ApproximateSplineSetFromGen(SplinePoint *from, SplinePoint *to, + bigreal start_t, bigreal end_t, + bigreal toler, int toler_is_sumsq, + GenPointsP genp, void *tok, + int order2, int depth) { + int cnt, i, maxerri=0, created = false; + bigreal errsum=0, maxerr=0, d, mid_t; + FitPoint *fp; + SplinePoint *mid, *r; + + cnt = (*genp)(tok, start_t, end_t, &fp); + if ( cnt < 2 ) + return NULL; + + // Rescale zero to one + for ( i=1; i<(cnt-1); ++i ) + fp[i].t = (fp[i].t-fp[0].t)/(fp[cnt-1].t-fp[0].t); + fp[0].t = 0.0; + fp[cnt-1].t = 1.0; + + from->nextcp.x = from->me.x + fp[0].ut.x; + from->nextcp.y = from->me.y + fp[0].ut.y; + from->nonextcp = false; + if ( to!=NULL ) + to->me = fp[cnt-1].p; + else { + to = SplinePointCreate(fp[cnt-1].p.x, fp[cnt-1].p.y); + created = true; + } + to->prevcp.x = to->me.x - fp[cnt-1].ut.x; + to->prevcp.y = to->me.y - fp[cnt-1].ut.y; + to->noprevcp = false; + ApproximateSplineFromPointsSlopes(from,to,fp+1,cnt-2,order2,mt_matrix); + + for ( i=0; i<cnt; ++i ) { + d = SplineMinDistanceToPoint(from->next, &fp[i].p); + errsum += d*d; + if ( d>maxerr ) { + maxerr = d; + maxerri = i; + } + } + // printf(" Error sum %lf, max error %lf at depth %d\n", errsum, maxerr, depth); + + if ( (toler_is_sumsq ? errsum : maxerr) > toler && depth < 6 ) { + mid_t = fp[maxerri].t * (end_t-start_t) + start_t; + free(fp); + SplineFree(from->next); + from->next = NULL; + to->prev = NULL; + mid = _ApproximateSplineSetFromGen(from, NULL, start_t, mid_t, toler, + toler_is_sumsq, genp, tok, order2, + depth+1); + if ( mid ) { + r = _ApproximateSplineSetFromGen(mid, to, mid_t, end_t, toler, + toler_is_sumsq, genp, tok, + order2, depth+1); + if ( r ) + return r; + else { + if ( created ) + SplinePointFree(to); + else + to->prev = NULL; + SplinePointFree(mid); + SplineFree(from->next); + from->next = NULL; + return NULL; + } + } else { + if ( created ) + SplinePointFree(to); + return NULL; + } + } else if ( (toler_is_sumsq ? errsum : maxerr) > toler ) { + TRACE("%s %lf exceeds %lf at maximum depth %d\n", + toler_is_sumsq ? "Sum of squared errors" : "Maximum error length", + toler_is_sumsq ? errsum : maxerr, toler, depth); + } + free(fp); + return to; +} + +SplinePoint *ApproximateSplineSetFromGen(SplinePoint *from, SplinePoint *to, + bigreal start_t, bigreal end_t, + bigreal toler, int toler_is_sumsq, + GenPointsP genp, void *tok, + int order2) { + return _ApproximateSplineSetFromGen(from, to, start_t, end_t, toler, + toler_is_sumsq, genp, tok, order2, 0); +} diff --git a/src/path/splinefit/splinefit.h b/src/path/splinefit/splinefit.h new file mode 100644 index 0000000..22e21fc --- /dev/null +++ b/src/path/splinefit/splinefit.h @@ -0,0 +1,78 @@ +// SPDX-License-Identifier: GPL-2.0-or-later +/* Copyright (C) 2000-2012 by George Williams, 2019 by Skef Iterum, 2021 by Linus Romer */ +/* + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + + * The name of the author may not be used to endorse or promote products + * derived from this software without specific prior written permission. + + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED + * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO + * EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; + * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, + * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR + * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF + * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef FONTFORGE_SPLINEFIT_H +#define FONTFORGE_SPLINEFIT_H + +// #include <fontforge-config.h> +#include "splinefont.h" + +typedef struct fitpoint { + BasePoint p; + BasePoint ut; + bigreal t; +} FitPoint; + +#define FITPOINT_EMPTY { {0.0, 0.0}, {0.0, 0.0}, 0.0 } + +enum mergetype { mt_matrix, mt_levien, mt_bruteforce }; + +extern Spline *ApproximateSplineFromPoints(SplinePoint *from, SplinePoint *to, + FitPoint *mid, int cnt, int order2); +extern Spline *ApproximateSplineFromPointsSlopes(SplinePoint *from, SplinePoint *to, + FitPoint *mid, int cnt, int order2, enum mergetype mt); + +/* ApproximateSplineSetFromGen() fits a one or more splines to data + * generated by calls to genp, within the tolerance toler. The data + * are treated as noiseless although some noise may be harmless. When + * an interval cannot be fit with a single spline it is divided at + * the point of largest error, where the position and slope are fixed + * according to the FitPoint entries. + * + * A GenPointsP function should allocate a free-able array of FitPoints, + * set *fpp to point to it, and return the number of points allocated. + * The points should correspond to values chosen between the interval + * t_start and t_end, but just what this interval corresponds to depends + * on vinfo and is therefore opaque to the tracing system. It is up to the + * implementation how many points to return on the interval and how they are + * "t-spaced". + * + * The function can also return 0 points, in which case *fpp should be set + * to NULL. This will cancel the tracing. This semantic is useful for + * communicating information found during tracing back to the caller + * of ApproximateSplineSetFromGen(), presumably via the vinfo structure, + * so that it can retrace with a different choice of interval. + */ +typedef int (*GenPointsP)(void *vinfo, bigreal t_start, bigreal t_end, FitPoint **fpp); + +extern SplinePoint *ApproximateSplineSetFromGen(SplinePoint *from, SplinePoint *to, + bigreal start_t, bigreal end_t, + bigreal toler, int toler_is_sumsq, + GenPointsP genp, void *tok, int order2); + +#endif /* FONTFORGE_SPLINEFIT_H */ diff --git a/src/path/splinefit/splinefont.c b/src/path/splinefit/splinefont.c new file mode 100644 index 0000000..4681628 --- /dev/null +++ b/src/path/splinefit/splinefont.c @@ -0,0 +1,1174 @@ +// SPDX-License-Identifier: GPL-2.0-or-later + +#include <math.h> +#include <stdbool.h> +#include <stdint.h> +#include "splinefont.h" +#include "splinefit.h" + +#define FONTFORGE_CONFIG_USE_DOUBLE 1 + +bigreal BPDot(BasePoint v1, BasePoint v2) { + return v1.x * v2.x + v1.y * v2.y; +} + +bigreal BPCross(BasePoint v1, BasePoint v2) { + return v1.x * v2.y - v1.y * v2.x; +} + +BasePoint BPRev(BasePoint v) { + return (BasePoint) { -v.x, -v.y }; +} + +int RealWithin(real a,real b,real fudge) { + return( b>=a-fudge && b<=a+fudge ); +} + +BOOL RealNear(real a,real b) { + real d = a-b; +#ifdef FONTFORGE_CONFIG_USE_DOUBLE + // These tighter equals-zero tests are retained for code tuned when + // passing zero as a constant + if ( a==0 ) + return b>-1e-8 && b<1e-8; + if ( b==0 ) + return a>-1e-8 && a<1e-8; + + return d>-1e-6 && d<1e-6; +#else /* For floats */ + return d>-1e-5 && d<1e-5 +#endif +} + +int RealApprox(real a,real b) { + + if ( a==0 ) { + if ( b<.0001 && b>-.0001 ) +return( true ); + } else if ( b==0 ) { + if ( a<.0001 && a>-.0001 ) +return( true ); + } else { + a /= b; + if ( a>=.95 && a<=1.05 ) +return( true ); + } +return( false ); +} + +void LineListFree(LineList *ll) { + LineList *next; + + while ( ll!=NULL ) { + next = ll->next; + chunkfree(ll,sizeof(LineList)); + ll = next; + } +} + +void LinearApproxFree(LinearApprox *la) { + LinearApprox *next; + + while ( la!=NULL ) { + next = la->next; + LineListFree(la->lines); + chunkfree(la,sizeof(LinearApprox)); + la = next; + } +} + +void SplineFree(Spline *spline) { + LinearApproxFree(spline->approx); + chunkfree(spline,sizeof(Spline)); +} + +SplinePoint *SplinePointCreate(real x, real y) { + SplinePoint *sp; + if ( (sp=chunkalloc(sizeof(SplinePoint)))!=NULL ) { + sp->me.x = x; sp->me.y = y; + sp->nextcp = sp->prevcp = sp->me; + sp->nonextcp = sp->noprevcp = true; + sp->nextcpdef = sp->prevcpdef = false; + sp->ttfindex = sp->nextcpindex = 0xfffe; + sp->name = NULL; + } + return( sp ); +} + +void SplinePointsFree(SplinePointList *spl) { + Spline *first, *spline, *next; + int nonext; + + if ( spl==NULL ) + return; + if ( spl->first!=NULL ) { + nonext = spl->first->next==NULL; // If there is no spline, we set a flag. + first = NULL; + // We start on the first spline if it exists. + for ( spline = spl->first->next; spline!=NULL && spline!=first; spline = next ) { + next = spline->to->next; // Cache the location of the next spline. + SplinePointFree(spline->to); // Free the destination point. + SplineFree(spline); // Free the spline. + if ( first==NULL ) first = spline; // We want to avoid repeating the circuit. + } + // If the path is open or has no splines, free the starting point. + if ( spl->last!=spl->first || nonext ) + SplinePointFree(spl->first); + } +} + +void SplinePointListFree(SplinePointList *spl) { + + if ( spl==NULL ) return; + SplinePointsFree(spl); + // free(spl->spiros); + free(spl->contour_name); + chunkfree(spl,sizeof(SplinePointList)); +} + +void SplineRefigure2(Spline *spline) { + SplinePoint *from = spline->from, *to = spline->to; + Spline1D *xsp = &spline->splines[0], *ysp = &spline->splines[1]; + Spline old; + +#ifdef DEBUG + if ( RealNear(from->me.x,to->me.x) && RealNear(from->me.y,to->me.y)) + IError("Zero length spline created"); +#endif + if ( spline->acceptableextrema ) + old = *spline; + + if ( ( from->nextcp.x==from->me.x && from->nextcp.y==from->me.y && from->nextcpindex>=0xfffe ) + || ( to->prevcp.x==to->me.x && to->prevcp.y==to->me.y && from->nextcpindex>=0xfffe ) ) { + from->nonextcp = to->noprevcp = true; + from->nextcp = from->me; + to->prevcp = to->me; + } else { + from->nonextcp = to->noprevcp = false; + if ( from->nextcp.x==from->me.x && from->nextcp.y==from->me.y ) + to->prevcp = from->me; + else if ( to->prevcp.x==to->me.x && to->prevcp.y==to->me.y ) + from->nextcp = to->me; + } + + if ( from->nonextcp && to->noprevcp ) + /* Ok */; + else if ( from->nextcp.x!=to->prevcp.x || from->nextcp.y!=to->prevcp.y ) { + if ( RealNear(from->nextcp.x,to->prevcp.x) && + RealNear(from->nextcp.y,to->prevcp.y)) { + from->nextcp.x = to->prevcp.x = (from->nextcp.x+to->prevcp.x)/2; + from->nextcp.y = to->prevcp.y = (from->nextcp.y+to->prevcp.y)/2; + } else { + IError("Invalid 2nd order spline in SplineRefigure2" ); +#ifndef GWW_TEST + /* I don't want these to go away when I'm debugging. I want to */ + /* know how I got them */ + from->nextcp.x = to->prevcp.x = (from->nextcp.x+to->prevcp.x)/2; + from->nextcp.y = to->prevcp.y = (from->nextcp.y+to->prevcp.y)/2; +#endif + } + } + + xsp->d = from->me.x; ysp->d = from->me.y; + if ( from->nonextcp && to->noprevcp ) { + spline->islinear = true; + xsp->c = to->me.x-from->me.x; + ysp->c = to->me.y-from->me.y; + xsp->a = xsp->b = 0; + ysp->a = ysp->b = 0; + } else { + /* from p. 393 (Operator Details, curveto) PostScript Lang. Ref. Man. (Red book) */ + xsp->c = 2*(from->nextcp.x-from->me.x); + ysp->c = 2*(from->nextcp.y-from->me.y); + xsp->b = to->me.x-from->me.x-xsp->c; + ysp->b = to->me.y-from->me.y-ysp->c; + xsp->a = 0; + ysp->a = 0; + if ( RealNear(xsp->c,0)) xsp->c=0; + if ( RealNear(ysp->c,0)) ysp->c=0; + if ( RealNear(xsp->b,0)) xsp->b=0; + if ( RealNear(ysp->b,0)) ysp->b=0; + spline->islinear = false; + if ( ysp->b==0 && xsp->b==0 ) + spline->islinear = true; /* This seems extremely unlikely... */ + if ( from->nextcpselected || to->prevcpselected ) { + // The convention for tracking selection of quadratic control + // points is to use nextcpselected except at the tail of the + // list, where it's prevcpselected on the first point. + from->nextcpselected = true; + to->prevcpselected = false; + } + } + if ( isnan(ysp->b) || isnan(xsp->b) ) + IError("NaN value in spline creation"); + LinearApproxFree(spline->approx); + spline->approx = NULL; + spline->knowncurved = false; + spline->knownlinear = spline->islinear; + SplineIsLinear(spline); + spline->isquadratic = !spline->knownlinear; + spline->order2 = true; + + if ( spline->acceptableextrema ) { + /* I don't check "d", because changes to that reflect simple */ + /* translations which will not affect the shape of the spline */ + /* (I don't check "a" because it is always 0 in a quadratic spline) */ + if ( !RealNear(old.splines[0].b,spline->splines[0].b) || + !RealNear(old.splines[0].c,spline->splines[0].c) || + !RealNear(old.splines[1].b,spline->splines[1].b) || + !RealNear(old.splines[1].c,spline->splines[1].c) ) + spline->acceptableextrema = false; + } +} + +Spline *SplineMake(SplinePoint *from, SplinePoint *to, int order2) { + if (order2 > 0) +return( SplineMake2(from,to)); + else +return( SplineMake3(from,to)); +} + +Spline *SplineMake2(SplinePoint *from, SplinePoint *to) { + Spline *spline = chunkalloc(sizeof(Spline)); + + spline->from = from; spline->to = to; + from->next = to->prev = spline; + spline->order2 = true; + SplineRefigure2(spline); +return( spline ); +} + +Spline *SplineMake3(SplinePoint *from, SplinePoint *to) { + Spline *spline = chunkalloc(sizeof(Spline)); + + spline->from = from; spline->to = to; + from->next = to->prev = spline; + SplineRefigure3(spline); +return( spline ); +} + +void SplinePointFree(SplinePoint *sp) { + // chunkfree(sp->hintmask,sizeof(HintMask)); + free(sp->name); + chunkfree(sp,sizeof(SplinePoint)); +} + +void SplineRefigure(Spline *spline) { + if ( spline==NULL ) +return; + if ( spline->order2 ) + SplineRefigure2(spline); + else + SplineRefigure3(spline); +} + +# define RE_NearZero .00000001 +# define RE_Factor (1024.0*1024.0*1024.0*1024.0*1024.0*2.0) /* 52 bits => divide by 2^51 */ + +int Within16RoundingErrors(bigreal v1, bigreal v2) { + bigreal temp=v1*v2; + bigreal re; + + if ( temp<0 ) /* Ok, if the two values are on different sides of 0 there */ +return( false ); /* is no way they can be within a rounding error of each other */ + else if ( temp==0 ) { + if ( v1==0 ) +return( v2<RE_NearZero && v2>-RE_NearZero ); + else +return( v1<RE_NearZero && v1>-RE_NearZero ); + } else if ( v1>0 ) { + if ( v1>v2 ) { /* Rounding error from the biggest absolute value */ + re = v1/ (RE_Factor/16); +return( v1-v2 < re ); + } else { + re = v2/ (RE_Factor/16); +return( v2-v1 < re ); + } + } else { + if ( v1<v2 ) { + re = v1/ (RE_Factor/16); /* This will be a negative number */ +return( v1-v2 > re ); + } else { + re = v2/ (RE_Factor/16); +return( v2-v1 > re ); + } + } +} + +/* An IEEE double has 52 bits of precision. So one unit of rounding error will be */ +/* the number divided by 2^51 */ +# define D_RE_Factor (1024.0*1024.0*1024.0*1024.0*1024.0*2.0) +/* But that's not going to work near 0, so, since the t values we care about */ +/* are [0,1], let's use 1.0/D_RE_Factor */ + +double CheckExtremaForSingleBitErrors(const Spline1D *sp, double t, double othert) { + double u1, um1; + double slope, slope1, slopem1; + int err; + double diff, factor; + + if ( t<0 || t>1 ) +return( t ); + + factor = t*0x40000/D_RE_Factor; + if ( (diff = t-othert)<0 ) diff= -diff; + if ( factor>diff/4 && diff!=0 ) /* This little check is to insure we don't skip beyond the well of this extremum into the next */ + factor = diff/4; + + slope = (3*(double) sp->a*t+2*sp->b)*t+sp->c; + if ( slope<0 ) slope = -slope; + + for ( err = 0x40000; err!=0; err>>=1 ) { + u1 = t+factor; + slope1 = (3*(double) sp->a*u1+2*sp->b)*u1+sp->c; + if ( slope1<0 ) slope1 = -slope1; + + um1 = t-factor; + slopem1 = (3*(double) sp->a*um1+2*sp->b)*um1+sp->c; + if ( slopem1<0 ) slopem1 = -slopem1; + + if ( slope1<slope && slope1<=slopem1 && u1<=1.0 ) { + t = u1; + } else if ( slopem1<slope && slopem1<=slope1 && um1>=0.0 ) { + t = um1; + } + factor /= 2.0; + } + /* that seems as good as it gets */ + +return( t ); +} + +void SplineFindExtrema(const Spline1D *sp, extended *_t1, extended *_t2 ) { + extended t1= -1, t2= -1; + extended b2_fourac; + + /* Find the extreme points on the curve */ + /* Set to -1 if there are none or if they are outside the range [0,1] */ + /* Order them so that t1<t2 */ + /* If only one valid extremum it will be t1 */ + /* (Does not check the end points unless they have derivative==0) */ + /* (Does not check to see if d/dt==0 points are inflection points (rather than extrema) */ + if ( sp->a!=0 ) { + /* cubic, possibly 2 extrema (possibly none) */ + b2_fourac = 4*(extended) sp->b*sp->b - 12*(extended) sp->a*sp->c; + if ( b2_fourac>=0 ) { + b2_fourac = sqrt(b2_fourac); + t1 = (-2*sp->b - b2_fourac) / (6*sp->a); + t2 = (-2*sp->b + b2_fourac) / (6*sp->a); + t1 = CheckExtremaForSingleBitErrors(sp,t1,t2); + t2 = CheckExtremaForSingleBitErrors(sp,t2,t1); + if ( t1>t2 ) { extended temp = t1; t1 = t2; t2 = temp; } + else if ( t1==t2 ) t2 = -1; + if ( RealNear(t1,0)) t1=0; else if ( RealNear(t1,1)) t1=1; + if ( RealNear(t2,0)) t2=0; else if ( RealNear(t2,1)) t2=1; + if ( t2<=0 || t2>=1 ) t2 = -1; + if ( t1<=0 || t1>=1 ) { t1 = t2; t2 = -1; } + } + } else if ( sp->b!=0 ) { + /* Quadratic, at most one extremum */ + t1 = -sp->c/(2.0*(extended) sp->b); + if ( t1<=0 || t1>=1 ) t1 = -1; + } else /*if ( sp->c!=0 )*/ { + /* linear, no extrema */ + } + *_t1 = t1; *_t2 = t2; +} + +int IntersectLines(BasePoint *inter, + BasePoint *line1_1, BasePoint *line1_2, + BasePoint *line2_1, BasePoint *line2_2) { + // A lot of functions call this with the same address as an input and the output. + // In order to avoid unexpected behavior, we delay writing to the output until the end. + bigreal s1, s2; + BasePoint _output; + BasePoint * output = &_output; + if ( line1_1->x == line1_2->x ) { + // Line 1 is vertical. + output->x = line1_1->x; + if ( line2_1->x == line2_2->x ) { + // Line 2 is vertical. + if ( line2_1->x!=line1_1->x ) + return( false ); /* Parallel vertical lines */ + output->y = (line1_1->y+line2_1->y)/2; + } else { + output->y = line2_1->y + (output->x-line2_1->x) * (line2_2->y - line2_1->y)/(line2_2->x - line2_1->x); + } + *inter = *output; + return( true ); + } else if ( line2_1->x == line2_2->x ) { + // Line 2 is vertical, but we know that line 1 is not. + output->x = line2_1->x; + output->y = line1_1->y + (output->x-line1_1->x) * (line1_2->y - line1_1->y)/(line1_2->x - line1_1->x); + *inter = *output; + return( true ); + } else { + // Both lines are oblique. + s1 = (line1_2->y - line1_1->y)/(line1_2->x - line1_1->x); + s2 = (line2_2->y - line2_1->y)/(line2_2->x - line2_1->x); + if ( RealNear(s1,s2)) { + if ( !RealNear(line1_1->y + (line2_1->x-line1_1->x) * s1,line2_1->y)) + return( false ); + output->x = (line1_2->x+line2_2->x)/2; + output->y = (line1_2->y+line2_2->y)/2; + } else { + output->x = (s1*line1_1->x - s2*line2_1->x - line1_1->y + line2_1->y)/(s1-s2); + output->y = line1_1->y + (output->x-line1_1->x) * s1; + } + *inter = *output; + return( true ); + } +} + +static int MinMaxWithin(Spline *spline) { + extended dx, dy; + int which; + extended t1, t2; + extended w; + /* We know that this "spline" is basically one dimensional. As long as its*/ + /* extrema are between the start and end points on that line then we can */ + /* treat it as a line. If the extrema are way outside the line segment */ + /* then it's a line that backtracks on itself */ + + if ( (dx = spline->to->me.x - spline->from->me.x)<0 ) dx = -dx; + if ( (dy = spline->to->me.y - spline->from->me.y)<0 ) dy = -dy; + which = dx<dy; + SplineFindExtrema(&spline->splines[which],&t1,&t2); + if ( t1==-1 ) +return( true ); + w = ((spline->splines[which].a*t1 + spline->splines[which].b)*t1 + + spline->splines[which].c)*t1 + spline->splines[which].d; + if ( RealNear(w, (&spline->to->me.x)[which]) || RealNear(w, (&spline->from->me.x)[which]) ) + /* Close enough */; + else if ( (w<(&spline->to->me.x)[which] && w<(&spline->from->me.x)[which]) || + (w>(&spline->to->me.x)[which] && w>(&spline->from->me.x)[which]) ) +return( false ); /* Outside */ + + w = ((spline->splines[which].a*t2 + spline->splines[which].b)*t2 + + spline->splines[which].c)*t2 + spline->splines[which].d; + if ( RealNear(w, (&spline->to->me.x)[which]) || RealNear(w, (&spline->from->me.x)[which]) ) + /* Close enough */; + else if ( (w<(&spline->to->me.x)[which] && w<(&spline->from->me.x)[which]) || + (w>(&spline->to->me.x)[which] && w>(&spline->from->me.x)[which]) ) +return( false ); /* Outside */ + +return( true ); +} + +int SplineIsLinear(Spline *spline) { + bigreal t1,t2, t3,t4; + int ret; + + if ( spline->knownlinear ) +return( true ); + if ( spline->knowncurved ) +return( false ); + + if ( spline->splines[0].a==0 && spline->splines[0].b==0 && + spline->splines[1].a==0 && spline->splines[1].b==0 ) +return( true ); + + /* Something is linear if the control points lie on the line between the */ + /* two base points */ + + /* Vertical lines */ + if ( RealNear(spline->from->me.x,spline->to->me.x) ) { + ret = RealNear(spline->from->me.x,spline->from->nextcp.x) && + RealNear(spline->from->me.x,spline->to->prevcp.x); + if ( ret && ! ((spline->from->nextcp.y >= spline->from->me.y && + spline->from->nextcp.y <= spline->to->me.y && + spline->to->prevcp.y >= spline->from->me.y && + spline->to->prevcp.y <= spline->to->me.y ) || + (spline->from->nextcp.y <= spline->from->me.y && + spline->from->nextcp.y >= spline->to->me.y && + spline->to->prevcp.y <= spline->from->me.y && + spline->to->prevcp.y >= spline->to->me.y )) ) + ret = MinMaxWithin(spline); + /* Horizontal lines */ + } else if ( RealNear(spline->from->me.y,spline->to->me.y) ) { + ret = RealNear(spline->from->me.y,spline->from->nextcp.y) && + RealNear(spline->from->me.y,spline->to->prevcp.y); + if ( ret && ! ((spline->from->nextcp.x >= spline->from->me.x && + spline->from->nextcp.x <= spline->to->me.x && + spline->to->prevcp.x >= spline->from->me.x && + spline->to->prevcp.x <= spline->to->me.x) || + (spline->from->nextcp.x <= spline->from->me.x && + spline->from->nextcp.x >= spline->to->me.x && + spline->to->prevcp.x <= spline->from->me.x && + spline->to->prevcp.x >= spline->to->me.x)) ) + ret = MinMaxWithin(spline); + } else { + ret = true; + t1 = (spline->from->nextcp.y-spline->from->me.y)/(spline->to->me.y-spline->from->me.y); + t2 = (spline->from->nextcp.x-spline->from->me.x)/(spline->to->me.x-spline->from->me.x); + t3 = (spline->to->me.y-spline->to->prevcp.y)/(spline->to->me.y-spline->from->me.y); + t4 = (spline->to->me.x-spline->to->prevcp.x)/(spline->to->me.x-spline->from->me.x); + ret = (Within16RoundingErrors(t1,t2) || (RealApprox(t1,0) && RealApprox(t2,0))) && + (Within16RoundingErrors(t3,t4) || (RealApprox(t3,0) && RealApprox(t4,0))); + if ( ret ) { + if ( t1<0 || t2<0 || t3<0 || t4<0 || + t1>1 || t2>1 || t3>1 || t4>1 ) + ret = MinMaxWithin(spline); + } + } + spline->knowncurved = !ret; + spline->knownlinear = ret; + if ( ret ) { + /* A few places that if the spline is knownlinear then its splines[?] */ + /* are linear. So give the linear version and not that suggested by */ + /* the control points */ + spline->splines[0].a = spline->splines[0].b = 0; + spline->splines[0].d = spline->from->me.x; + spline->splines[0].c = spline->to->me.x-spline->from->me.x; + spline->splines[1].a = spline->splines[1].b = 0; + spline->splines[1].d = spline->from->me.y; + spline->splines[1].c = spline->to->me.y-spline->from->me.y; + } +return( ret ); +} + +static bigreal FindZero5(bigreal w[7],bigreal tlow, bigreal thigh) { + /* Somewhere between tlow and thigh there is a value of t where w(t)==0 */ + /* It is conceiveable that there might be 3 such ts if there are some high frequency effects */ + /* but I ignore that for now */ + bigreal t, test; + int bot_negative; + + t = tlow; + test = ((((w[5]*t+w[4])*t+w[3])*t+w[2])*t+w[1])*t + w[0]; + bot_negative = test<0; + + for (;;) { + t = (thigh+tlow)/2; + if ( thigh==t || tlow==t ) +return( t ); /* close as we can get */ + test = ((((w[5]*t+w[4])*t+w[3])*t+w[2])*t+w[1])*t + w[0]; + if ( test==0 ) +return( t ); + if ( bot_negative ) { + if ( test<0 ) + tlow = t; + else + thigh = t; + } else { + if ( test<0 ) + thigh = t; + else + tlow = t; + } + } +} + +static bigreal FindZero3(bigreal w[7],bigreal tlow, bigreal thigh) { + /* Somewhere between tlow and thigh there is a value of t where w(t)==0 */ + /* It is conceiveable that there might be 3 such ts if there are some high frequency effects */ + /* but I ignore that for now */ + bigreal t, test; + int bot_negative; + + t = tlow; + test = ((w[3]*t+w[2])*t+w[1])*t + w[0]; + bot_negative = test<0; + + for (;;) { + t = (thigh+tlow)/2; + if ( thigh==t || tlow==t ) +return( t ); /* close as we can get */ + test = ((w[3]*t+w[2])*t+w[1])*t + w[0]; + if ( test==0 ) +return( t ); + if ( bot_negative ) { + if ( test<0 ) + tlow = t; + else + thigh = t; + } else { + if ( test<0 ) + thigh = t; + else + tlow = t; + } + } +} + +bigreal SplineMinDistanceToPoint(Spline *s, BasePoint *p) { + /* So to find the minimum distance we want the sqrt( (sx(t)-px)^2 + (sy(t)-py)^2 ) */ + /* Same minima as (sx(t)-px)^2 + (sy(t)-py)^2, which is easier to deal with */ + bigreal w[7]; + Spline1D *x = &s->splines[0], *y = &s->splines[1]; + bigreal off[2], best; + + off[0] = (x->d-p->x); off[1] = (y->d-p->y); + + w[6] = (x->a*x->a) + (y->a*y->a); + w[5] = 2*(x->a*x->b + y->a*y->b); + w[4] = (x->b*x->b) + 2*(x->a*x->c) + (y->b*y->b) + 2*(y->a*y->c); + w[3] = 2* (x->b*x->c + x->a*off[0] + y->b*y->c + y->a*off[1]); + w[2] = (x->c*x->c) + 2*(x->b*off[0]) + (y->c*y->c) + 2*y->b*off[1]; + w[1] = 2*(x->c*off[0] + y->c*off[1]); + w[0] = off[0]*off[0] + off[1]*off[1]; + + /* Take derivative */ + w[0] = w[1]; + w[1] = 2*w[2]; + w[2] = 3*w[3]; + w[3] = 4*w[4]; + w[4] = 5*w[5]; + w[5] = 6*w[6]; + w[6] = 0; + + if ( w[5]!=0 ) { + bigreal tzeros[8], t, incr, test, lasttest, zerot; + int i, zcnt=0; + /* Well, we've got a 5th degree poly and no way to play cute tricks. */ + /* brute force it */ + incr = 1.0/1024; + lasttest = w[0]; + for ( t = incr; t<=1.0; t += incr ) { + test = ((((w[5]*t+w[4])*t+w[3])*t+w[2])*t+w[1])*t + w[0]; + if ( test==0 ) + tzeros[zcnt++] = t; + else { + if ( lasttest!=0 && (test>0) != (lasttest>0) ) { + zerot = FindZero5(w,t-incr,t); + if ( zerot>0 ) + tzeros[zcnt++] = zerot; + } + } + lasttest = test; + } + best = off[0]*off[0] + off[1]*off[1]; /* t==0 */ + test = (x->a+x->b+x->c+off[0])*(x->a+x->b+x->c+off[0]) + + (y->a+y->b+y->c+off[1])*(y->a+y->b+y->c+off[1]); /* t==1 */ + if ( best>test ) best = test; + for ( i=0; i<zcnt; ++i ) { + bigreal tx, ty; + tx = ((x->a*tzeros[i]+x->b)*tzeros[i]+x->c)*tzeros[i] + off[0]; + ty = ((y->a*tzeros[i]+y->b)*tzeros[i]+y->c)*tzeros[i] + off[1]; + test = tx*tx + ty*ty; + if ( best>test ) best = test; + } +return( sqrt(best)); + } else if ( w[4]==0 && w[3]!=0 ) { + /* Started with a quadratic -- now, find 0s of a cubic */ + /* We could find the extrema, so we have a bunch of monotonics */ + /* Or we could brute force it as above */ + bigreal tzeros[8], test, zerot; + bigreal quad[3], disc, e[5], t1, t2; + int i, zcnt=0, ecnt; + + quad[2] = 3*w[3]; quad[1] = 2*w[2]; quad[0] = w[1]; + disc = (-quad[1]*quad[1] - 4*quad[2]*quad[0]); + e[0] = 0; + if ( disc<0 ) { + e[1] = 1.0; + ecnt = 2; + } else + disc = sqrt(disc); + t1 = (-w[1] - disc) / (2*w[2]); + t2 = (-w[1] + disc) / (2*w[2]); + if ( t1>t2 ) { + bigreal temp = t1; + t1 = t2; + t2 = temp; + } + ecnt=1; + if ( t1>0 && t1<1 ) + e[ecnt++] = t1; + if ( t2>0 && t2<1 && t1!=t2 ) + e[ecnt++] = t2; + e[ecnt++] = 1.0; + for ( i=1; i<ecnt; ++i ) { + zerot = FindZero3(w,e[i-1],e[i]); + if ( zerot>0 ) + tzeros[zcnt++] = zerot; + } + best = off[0]*off[0] + off[1]*off[1]; /* t==0 */ + test = (x->b+x->c+off[0])*(x->b+x->c+off[0]) + + (y->b+y->c+off[1])*(y->b+y->c+off[1]); /* t==1 */ + if ( best>test ) best = test; + for ( i=0; i<zcnt; ++i ) { + bigreal tx, ty; + tx = (x->b*tzeros[i]+x->c)*tzeros[i] + off[0]; + ty = (y->b*tzeros[i]+y->c)*tzeros[i] + off[1]; + test = tx*tx + ty*ty; + if ( best>test ) best = test; + } +return( sqrt(best)); + } else if ( w[2]==0 && w[1]!=0 ) { + /* Started with a line */ + bigreal t = -w[0]/w[1], test, best; + best = off[0]*off[0] + off[1]*off[1]; /* t==0 */ + test = (x->c+off[0])*(x->c+off[0]) + (y->c+off[1])*(y->c+off[1]); /* t==1 */ + if ( best>test ) best = test; + if ( t>0 && t<1 ) { + test = (x->c*t+off[0])*(x->c*t+off[0]) + (y->c*t+off[1])*(y->c*t+off[1]); + if ( best>test ) best = test; + } +return(sqrt(best)); + } else if ( w[4]!=0 && w[3]!=0 && w[2]!=0 && w[1]!=0 ) { + IError( "Impossible condition in SplineMinDistanceToPoint"); + } else { + /* It's a point, minimum distance is the only distance */ +return( sqrt(off[0]*off[0] + off[1]*off[1]) ); + } +return( -1 ); +} + +/* This returns all real solutions, even those out of bounds */ +/* I use -999999 as an error flag, since we're really only interested in */ +/* solns near 0 and 1 that should be ok. -1 is perhaps a little too close */ +/* Sigh. When solutions are near 0, the rounding errors are appalling. */ +int _CubicSolve(const Spline1D *sp,bigreal sought, extended ts[3]) { + extended d, xN, yN, delta2, temp, delta, h, t2, t3, theta; + extended sa=sp->a, sb=sp->b, sc=sp->c, sd=sp->d-sought; + int i=0; + + ts[0] = ts[1] = ts[2] = -999999; + if ( sd==0 && sa!=0 ) { + /* one of the roots is 0, the other two are the soln of a quadratic */ + ts[0] = 0; + if ( sc==0 ) { + ts[1] = -sb/(extended) sa; /* two zero roots */ + } else { + temp = sb*(extended) sb-4*(extended) sa*sc; + if ( RealNear(temp,0)) + ts[1] = -sb/(2*(extended) sa); + else if ( temp>=0 ) { + temp = sqrt(temp); + ts[1] = (-sb+temp)/(2*(extended) sa); + ts[2] = (-sb-temp)/(2*(extended) sa); + } + } + } else if ( sa!=0 ) { + /* http://www.m-a.org.uk/eb/mg/mg077ch.pdf */ + /* this nifty solution to the cubic neatly avoids complex arithmetic */ + xN = -sb/(3*(extended) sa); + yN = ((sa*xN + sb)*xN+sc)*xN + sd; + + delta2 = (sb*(extended) sb-3*(extended) sa*sc)/(9*(extended) sa*sa); + /*if ( RealWithin(delta2,0,.00000001) ) delta2 = 0;*/ + + /* the descriminant is yN^2-h^2, but delta might be <0 so avoid using h */ + d = yN*yN - 4*sa*sa*delta2*delta2*delta2; + if ( ((yN>.01 || yN<-.01) && RealNear(d/yN,0)) || ((yN<=.01 && yN>=-.01) && RealNear(d,0)) ) + d = 0; + if ( d>0 ) { + temp = sqrt(d); + t2 = (-yN-temp)/(2*sa); + t2 = (t2==0) ? 0 : (t2<0) ? -pow(-t2,1./3.) : pow(t2,1./3.); + t3 = (-yN+temp)/(2*sa); + t3 = t3==0 ? 0 : (t3<0) ? -pow(-t3,1./3.) : pow(t3,1./3.); + ts[0] = xN + t2 + t3; + } else if ( d<0 ) { + if ( delta2>=0 ) { + delta = sqrt(delta2); + h = 2*sa*delta2*delta; + temp = -yN/h; + if ( temp>=-1.0001 && temp<=1.0001 ) { + if ( temp<-1 ) temp = -1; else if ( temp>1 ) temp = 1; + theta = acos(temp)/3; + ts[i++] = xN+2*delta*cos(theta); + ts[i++] = xN+2*delta*cos(2.0943951+theta); /* 2*pi/3 */ + ts[i++] = xN+2*delta*cos(4.1887902+theta); /* 4*pi/3 */ + } + } + } else if ( /* d==0 && */ delta2!=0 ) { + delta = yN/(2*sa); + delta = delta==0 ? 0 : delta>0 ? pow(delta,1./3.) : -pow(-delta,1./3.); + ts[i++] = xN + delta; /* this root twice, but that's irrelevant to me */ + ts[i++] = xN - 2*delta; + } else if ( /* d==0 && */ delta2==0 ) { + if ( xN>=-0.0001 && xN<=1.0001 ) ts[0] = xN; + } + } else if ( sb!=0 ) { + extended d = sc*(extended) sc-4*(extended) sb*sd; + if ( d<0 && RealNear(d,0)) d=0; + if ( d<0 ) +return(false); /* All roots imaginary */ + d = sqrt(d); + ts[0] = (-sc-d)/(2*(extended) sb); + ts[1] = (-sc+d)/(2*(extended) sb); + } else if ( sc!=0 ) { + ts[0] = -sd/(extended) sc; + } else { + /* If it's a point then either everything is a solution, or nothing */ + } +return( ts[0]!=-999999 ); +} + +int _QuarticSolve(Quartic *q,extended ts[4]) { + extended extrema[5]; + Spline1D sp; + int ecnt = 0, i, zcnt; + + /* Two special cases */ + if ( q->a==0 ) { /* It's really a cubic */ + sp.a = q->b; + sp.b = q->c; + sp.c = q->d; + sp.d = q->e; + ts[3] = -999999; +return( _CubicSolve(&sp,0,ts)); + } else if ( q->e==0 ) { /* we can factor out a zero root */ + sp.a = q->a; + sp.b = q->b; + sp.c = q->c; + sp.d = q->d; + ts[0] = 0; +return( _CubicSolve(&sp,0,ts+1)+1); + } + + sp.a = 4*q->a; + sp.b = 3*q->b; + sp.c = 2*q->c; + sp.d = q->d; + if ( _CubicSolve(&sp,0,extrema)) { + ecnt = 1; + if ( extrema[1]!=-999999 ) { + ecnt = 2; + if ( extrema[1]<extrema[0] ) { + extended temp = extrema[1]; extrema[1] = extrema[0]; extrema[0]=temp; + } + if ( extrema[2]!=-999999 ) { + ecnt = 3; + if ( extrema[2]<extrema[0] ) { + extended temp = extrema[2]; extrema[2] = extrema[0]; extrema[0]=temp; + } + if ( extrema[2]<extrema[1] ) { + extended temp = extrema[2]; extrema[2] = extrema[1]; extrema[1]=temp; + } + } + } + } + for ( i=ecnt-1; i>=0 ; --i ) + extrema[i+1] = extrema[i]; + /* Upper and lower bounds within which we'll search */ + extrema[0] = -999; + extrema[ecnt+1] = 999; + ecnt += 2; + /* divide into monotonic sections & use binary search to find zeroes */ + for ( i=zcnt=0; i<ecnt-1; ++i ) { + extended top, bottom, val; + extended topt, bottomt, t; + topt = extrema[i+1]; + bottomt = extrema[i]; + top = (((q->a*topt+q->b)*topt+q->c)*topt+q->d)*topt+q->e; + bottom = (((q->a*bottomt+q->b)*bottomt+q->c)*bottomt+q->d)*bottomt+q->e; + if ( top<bottom ) { + extended temp = top; top = bottom; bottom = temp; + temp = topt; topt = bottomt; bottomt = temp; + } + if ( bottom>.001 ) /* this monotonic is all above 0 */ + continue; + if ( top<-.001 ) /* this monotonic is all below 0 */ + continue; + if ( bottom>0 ) { + ts[zcnt++] = bottomt; + continue; + } + if ( top<0 ) { + ts[zcnt++] = topt; + continue; + } + for (;;) { + t = (topt+bottomt)/2; + if ( isnan(t) ) { + break; + } else if ( t==topt || t==bottomt ) { + ts[zcnt++] = t; + break; + } + + val = (((q->a*t+q->b)*t+q->c)*t+q->d)*t+q->e; + if ( val>-.0001 && val<.0001 ) { + ts[zcnt++] = t; + break; + } else if ( val>0 ) { + top = val; + topt = t; + } else { + bottom = val; + bottomt = t; + } + } + } + for ( i=zcnt; i<4; ++i ) + ts[i] = -999999; +return( zcnt ); +} + +/* calculating the actual length of a spline is hard, this gives a very */ +/* rough (but quick) approximation */ +static bigreal SplineLenApprox(Spline *spline) { + bigreal len, slen, temp; + + if ( (temp = spline->to->me.x-spline->from->me.x)<0 ) temp = -temp; + len = temp; + if ( (temp = spline->to->me.y-spline->from->me.y)<0 ) temp = -temp; + len += temp; + if ( !spline->to->noprevcp || !spline->from->nonextcp ) { + if ( (temp = spline->from->nextcp.x-spline->from->me.x)<0 ) temp = -temp; + slen = temp; + if ( (temp = spline->from->nextcp.y-spline->from->me.y)<0 ) temp = -temp; + slen += temp; + if ( (temp = spline->to->prevcp.x-spline->from->nextcp.x)<0 ) temp = -temp; + slen += temp; + if ( (temp = spline->to->prevcp.y-spline->from->nextcp.y)<0 ) temp = -temp; + slen += temp; + if ( (temp = spline->to->me.x-spline->to->prevcp.x)<0 ) temp = -temp; + slen += temp; + if ( (temp = spline->to->me.y-spline->to->prevcp.y)<0 ) temp = -temp; + slen += temp; + len = (len + slen)/2; + } +return( len ); +} + +FitPoint *SplinesFigureFPsBetween(SplinePoint *from, SplinePoint *to, + int *tot) { + int cnt, i, j, pcnt; + bigreal len, slen, lbase; + SplinePoint *np; + FitPoint *fp; + bigreal _lens[10], *lens = _lens; + int _cnts[10], *cnts = _cnts; + /* I used just to give every spline 10 points. But that gave much more */ + /* weight to small splines than to big ones */ + + cnt = 0; + for ( np = from->next->to; ; np = np->next->to ) { + ++cnt; + if ( np==to ) + break; + } + if ( cnt>10 ) { + lens = malloc(cnt*sizeof(bigreal)); + cnts = malloc(cnt*sizeof(int)); + } + cnt = 0; len = 0; + for ( np = from->next->to; ; np = np->next->to ) { + lens[cnt] = SplineLenApprox(np->prev); + len += lens[cnt]; + ++cnt; + if ( np==to ) + break; + } + if ( len!=0 ) { + pcnt = 0; + for ( i=0; i<cnt; ++i ) { + int pnts = rint( (10*cnt*lens[i])/len ); + if ( pnts<2 ) pnts = 2; + cnts[i] = pnts; + pcnt += pnts; + } + } else + pcnt = 2*cnt; + + fp = malloc((pcnt+1)*sizeof(FitPoint)); i = 0; + if ( len==0 ) { + for ( i=0; i<=pcnt; ++i ) { + fp[i].t = i/(pcnt); + fp[i].p.x = from->me.x; + fp[i].p.y = from->me.y; + } + } else { + lbase = 0; + for ( i=cnt=0, np = from->next->to; ; np = np->next->to, ++cnt ) { + slen = SplineLenApprox(np->prev); + for ( j=0; j<cnts[cnt]; ++j ) { + bigreal t = j/(bigreal) cnts[cnt]; + fp[i].t = (lbase+ t*slen)/len; + fp[i].p.x = ((np->prev->splines[0].a*t+np->prev->splines[0].b)*t+np->prev->splines[0].c)*t + np->prev->splines[0].d; + fp[i++].p.y = ((np->prev->splines[1].a*t+np->prev->splines[1].b)*t+np->prev->splines[1].c)*t + np->prev->splines[1].d; + } + lbase += slen; + if ( np==to ) + break; + } + } + if ( cnts!=_cnts ) free(cnts); + if ( lens!=_lens ) free(lens); + + *tot = i; + +return( fp ); +} + +static int SplinePointCategory(SplinePoint *sp) { + enum pointtype pt; + + pt = pt_corner; + if ( sp->next==NULL && sp->prev==NULL ) + ; + else if ( (sp->next!=NULL && sp->next->to->me.x==sp->me.x && sp->next->to->me.y==sp->me.y) || + (sp->prev!=NULL && sp->prev->from->me.x==sp->me.x && sp->prev->from->me.y==sp->me.y )) + ; + else if ( sp->next==NULL ) { + pt = sp->noprevcp ? pt_corner : pt_curve; + } else if ( sp->prev==NULL ) { + pt = sp->nonextcp ? pt_corner : pt_curve; + } else if ( sp->nonextcp && sp->noprevcp ) { + ; + } else { + BasePoint ndir, ncdir, ncunit, pdir, pcdir, pcunit; + bigreal nlen, nclen, plen, pclen; + bigreal cross, bounds; + + ncdir.x = sp->nextcp.x - sp->me.x; ncdir.y = sp->nextcp.y - sp->me.y; + pcdir.x = sp->prevcp.x - sp->me.x; pcdir.y = sp->prevcp.y - sp->me.y; + ndir.x = ndir.y = pdir.x = pdir.y = 0; + if ( sp->next!=NULL ) { + ndir.x = sp->next->to->me.x - sp->me.x; ndir.y = sp->next->to->me.y - sp->me.y; + } + if ( sp->prev!=NULL ) { + pdir.x = sp->prev->from->me.x - sp->me.x; pdir.y = sp->prev->from->me.y - sp->me.y; + } + nclen = sqrt(ncdir.x*ncdir.x + ncdir.y*ncdir.y); + pclen = sqrt(pcdir.x*pcdir.x + pcdir.y*pcdir.y); + nlen = sqrt(ndir.x*ndir.x + ndir.y*ndir.y); + plen = sqrt(pdir.x*pdir.x + pdir.y*pdir.y); + ncunit = ncdir; pcunit = pcdir; + if ( nclen!=0 ) { ncunit.x /= nclen; ncunit.y /= nclen; } + if ( pclen!=0 ) { pcunit.x /= pclen; pcunit.y /= pclen; } + if ( nlen!=0 ) { ndir.x /= nlen; ndir.y /= nlen; } + if ( plen!=0 ) { pdir.x /= plen; pdir.y /= plen; } + + /* find out which side has the shorter control vector. Cross that vector */ + /* with the normal of the unit vector on the other side. If the */ + /* result is less than 1 em-unit then we've got colinear control points */ + /* (within the resolution of the integer grid) */ + /* Not quite... they could point in the same direction */ + if ( sp->pointtype==pt_curve ) + bounds = 4.0; + else + bounds = 1.0; + if ( nclen!=0 && pclen!=0 && + ((nclen>=pclen && (cross = pcdir.x*ncunit.y - pcdir.y*ncunit.x)<bounds && cross>-bounds ) || + (pclen>nclen && (cross = ncdir.x*pcunit.y - ncdir.y*pcunit.x)<bounds && cross>-bounds )) && + ncdir.x*pcdir.x + ncdir.y*pcdir.y < 0 ) + pt = pt_curve; + /* Cross product of control point with unit vector normal to line in */ + /* opposite direction should be less than an em-unit for a tangent */ + else if ( ( nclen==0 && pclen!=0 + && (cross = pcdir.x*ndir.y-pcdir.y*ndir.x)<bounds + && cross>-bounds && (pcdir.x*ndir.x+pcdir.y*ndir.y)<0 ) + || + ( pclen==0 && nclen!=0 + && (cross = ncdir.x*pdir.y-ncdir.y*pdir.x)<bounds + && cross>-bounds && (ncdir.x*pdir.x+ncdir.y*pdir.y)<0 ) ) + pt = pt_tangent; + + if (pt == pt_curve && + ((sp->nextcp.x==sp->me.x && sp->prevcp.x==sp->me.x && sp->nextcp.y!=sp->me.y) || + (sp->nextcp.y==sp->me.y && sp->prevcp.y==sp->me.y && sp->nextcp.x!=sp->me.x))) + pt = pt_hvcurve; + } + return pt; +} + +static enum pointtype SplinePointDowngrade(int current, int geom) { + enum pointtype np = current; + + if ( current==pt_curve && geom!=pt_curve ) { + if ( geom==pt_hvcurve ) + np = pt_curve; + else + np = pt_corner; + } else if ( current==pt_hvcurve && geom!=pt_hvcurve ) { + if ( geom==pt_curve ) + np = pt_curve; + else + np = pt_corner; + } else if ( current==pt_tangent && geom!=pt_tangent ) { + np = pt_corner; + } + + return np; +} + +// Assumes flag combinations are already verified. Only returns false +// when called with check_compat +int _SplinePointCategorize(SplinePoint *sp, int flags) { + enum pointtype geom, dg, cur; + + if ( flags & pconvert_flag_none ) + // No points selected for conversion -- keep type as is + return true; + if ( flags & pconvert_flag_smooth && sp->pointtype == pt_corner ) + // Convert only "smooth" points, not corners + return true; + + geom = SplinePointCategory(sp); + dg = SplinePointDowngrade(sp->pointtype, geom); + + if ( flags & pconvert_flag_incompat && sp->pointtype == dg ) + // Only convert points incompatible with current type + return true; + + if ( flags & pconvert_flag_by_geom ) { + if ( ! ( flags & pconvert_flag_hvcurve ) && geom == pt_hvcurve ) + sp->pointtype = pt_curve; + else + sp->pointtype = geom; + } else if ( flags & pconvert_flag_downgrade ) { + sp->pointtype = dg; + } else if ( flags & pconvert_flag_force_type ) { + if ( sp->pointtype != dg ) { + cur = sp->pointtype; + sp->pointtype = dg; + /* SPChangePointType(sp,cur); */ + } + } else if ( flags & pconvert_flag_check_compat ) { + if ( sp->pointtype != dg ) + return false; + } + return true; +} + +void SplinePointCategorize(SplinePoint *sp) { + _SplinePointCategorize(sp, pconvert_flag_all|pconvert_flag_by_geom); +} + +static void SplinePointReCategorize(SplinePoint *sp,int oldpt) { + SplinePointCategorize(sp); + if ( sp->pointtype!=oldpt ) { + if ( sp->pointtype==pt_curve && oldpt==pt_hvcurve && + ((sp->nextcp.x == sp->me.x && sp->nextcp.y != sp->me.y ) || + (sp->nextcp.y == sp->me.y && sp->nextcp.x != sp->me.x ))) + sp->pointtype = pt_hvcurve; + } +} + +void SplinesRemoveBetween(SplinePoint *from, SplinePoint *to, int type) { + int tot; + FitPoint *fp; + SplinePoint *np, oldfrom; + int oldfpt = from->pointtype, oldtpt = to->pointtype; + Spline *sp; + int order2 = from->next->order2; + + oldfrom = *from; + fp = SplinesFigureFPsBetween(from,to,&tot); + + if ( type==1 ) + ApproximateSplineFromPointsSlopes(from,to,fp,tot-1,order2,mt_levien); + else + ApproximateSplineFromPoints(from,to,fp,tot-1,order2); + + /* Have to do the frees after the approximation because the approx */ + /* uses the splines to determine slopes */ + for ( sp = oldfrom.next; ; ) { + np = sp->to; + SplineFree(sp); + if ( np==to ) + break; + sp = np->next; + // SplinePointMDFree(sc,np); + } + + free(fp); + + SplinePointReCategorize(from,oldfpt); + SplinePointReCategorize(to,oldtpt); +} diff --git a/src/path/splinefit/splinefont.h b/src/path/splinefit/splinefont.h new file mode 100644 index 0000000..83db934 --- /dev/null +++ b/src/path/splinefit/splinefont.h @@ -0,0 +1,191 @@ +// SPDX-License-Identifier: GPL-2.0-or-later + +#ifndef _SEEN_SPLINEFONT_H_ +#define _SEEN_SPLINEFONT_H_ + +#include <glib.h> + +typedef double real; +typedef double bigreal; +typedef double extended; +typedef int BOOL; + +#define chunkalloc(size) calloc(1,size) +#define chunkfree(item,size) free(item) + +typedef struct basepoint { + real x; + real y; +} BasePoint; + +typedef struct ipoint { + int x; + int y; +} IPoint; + +enum pointtype { pt_curve, pt_corner, pt_tangent, pt_hvcurve }; +typedef struct splinepoint { + BasePoint me; + BasePoint nextcp; /* control point */ + BasePoint prevcp; /* control point */ + unsigned int nonextcp:1; + unsigned int noprevcp:1; + unsigned int nextcpdef:1; + unsigned int prevcpdef:1; + unsigned int selected:1; /* for UI */ + unsigned int nextcpselected: 2; /* Is the next BCP selected */ + unsigned int prevcpselected: 2; /* Is the prev BCP selected */ + unsigned int pointtype:2; + unsigned int isintersection: 1; + unsigned int flexy: 1; /* When "freetype_markup" is on in charview.c:DrawPoint */ + unsigned int flexx: 1; /* flexy means select nextcp, and flexx means draw circle around nextcp */ + unsigned int roundx: 1; /* For true type hinting */ + unsigned int roundy: 1; /* For true type hinting */ + unsigned int dontinterpolate: 1; /* in ttf, don't imply point by interpolating between cps */ + unsigned int ticked: 1; + unsigned int watched: 1; + /* 1 bits left... */ + uint16_t ptindex; /* Temporary value used by metafont routine */ + uint16_t ttfindex; /* Truetype point index */ + /* Special values 0xffff => point implied by averaging control points */ + /* 0xfffe => point created with no real number yet */ + /* (or perhaps point in context where no number is possible as in a glyph with points & refs) */ + uint16_t nextcpindex; /* Truetype point index */ + struct spline *next; + struct spline *prev; + /* Inkscape: not used; HintMask *hintmask; */ + char* name; +} SplinePoint; + + +typedef struct spline1d { + real a, b, c, d; +} Spline1D; + +typedef struct spline { + unsigned int islinear: 1; /* No control points */ + unsigned int isquadratic: 1; /* probably read in from ttf */ + unsigned int isticked: 1; + unsigned int isneeded: 1; /* Used in remove overlap */ + unsigned int isunneeded: 1; /* Used in remove overlap */ + unsigned int exclude: 1; /* Used in remove overlap variant: exclude */ + unsigned int ishorvert: 1; + unsigned int knowncurved: 1; /* We know that it curves */ + unsigned int knownlinear: 1; /* it might have control points, but still traces out a line */ + /* If neither knownlinear nor curved then we haven't checked */ + unsigned int order2: 1; /* It's a bezier curve with only one cp */ + unsigned int touched: 1; + unsigned int leftedge: 1; + unsigned int rightedge: 1; + unsigned int acceptableextrema: 1; /* This spline has extrema, but we don't care */ + SplinePoint *from; + SplinePoint *to; + Spline1D splines[2]; /* splines[0] is the x spline, splines[1] is y */ + struct linearapprox *approx; + /* Possible optimizations: + Precalculate bounding box + Precalculate min/max/ points of inflection + */ +} Spline; + +typedef struct splinepointlist { + SplinePoint *first, *last; + struct splinepointlist *next; + /* Not used: spiro_cp *spiros; */ + uint16_t spiro_cnt, spiro_max; + /* These could be bit fields, but bytes are easier to access and we */ + /* don't need the space (yet) */ + uint8_t ticked; + uint8_t beziers_need_optimizer; /* If the spiros have changed in spiro mode, then reverting to bezier mode might, someday, run a simplifier */ + uint8_t is_clip_path; /* In type3/svg fonts */ + int start_offset; // This indicates which point is the canonical first for purposes of outputting to U. F. O.. + char *contour_name; +} SplinePointList, SplineSet; + +typedef struct dbounds { + real minx, maxx; + real miny, maxy; +} DBounds; + +typedef struct quartic { + bigreal a,b,c,d,e; +} Quartic; + + +int RealWithin(real a,real b,real fudge); +BOOL RealNear(real a, real b); + +Spline *SplineMake(SplinePoint *from, SplinePoint *to, int order2); +Spline *SplineMake2(SplinePoint *from, SplinePoint *to); +Spline *SplineMake3(SplinePoint *from, SplinePoint *to); +SplinePoint *SplinePointCreate(real x, real y); + +void SplineRefigure3(Spline *spline); +void SplineRefigure(Spline *spline); +int SplineIsLinear(Spline *spline); +void SplineFindExtrema(const Spline1D *sp, extended *_t1, extended *_t2 ); +bigreal SplineMinDistanceToPoint(Spline *s, BasePoint *p); + +void SplinePointFree(SplinePoint *sp); +void SplineFree(Spline *spline); +void SplinePointListFree(SplinePointList *spl); + +bigreal BPDot(BasePoint v1, BasePoint v2); +bigreal BPCross(BasePoint v1, BasePoint v2); +BasePoint BPRev(BasePoint v); + +int _CubicSolve(const Spline1D *sp,bigreal sought, extended ts[3]); +int _QuarticSolve(Quartic *q,extended ts[4]); +int IntersectLines(BasePoint *inter, + BasePoint *line1_1, BasePoint *line1_2, + BasePoint *line2_1, BasePoint *line2_2); + +#define IError(msg) g_warning(msg) +#define TRACE g_message + +enum linelist_flags { cvli_onscreen=0x1, cvli_clipped=0x2 }; + +typedef struct linelist { + IPoint here; + struct linelist *next; + /* The first two fields are constant for the linelist, the next ones */ + /* refer to a particular screen. If some portion of the line from */ + /* this point to the next one is on the screen then set cvli_onscreen */ + /* if this point needs to be clipped then set cvli_clipped */ + /* asend and asstart are the actual screen locations where this point */ + /* intersects the clip edge. */ + enum linelist_flags flags; + IPoint asend, asstart; +} LineList; + +typedef struct linearapprox { + real scale; + unsigned int oneline: 1; + unsigned int onepoint: 1; + unsigned int any: 1; /* refers to a particular screen */ + struct linelist *lines; + struct linearapprox *next; +} LinearApprox; + +void LinearApproxFree(LinearApprox *la); + +int Within16RoundingErrors(bigreal v1, bigreal v2); + +enum pconvert_flags { + // Point selection (mutually exclusive) + pconvert_flag_none = 0x01, + pconvert_flag_all = 0x02, + pconvert_flag_smooth = 0x04, + pconvert_flag_incompat = 0x08, + // Conversion modes (mutually exclusive) + pconvert_flag_by_geom = 0x100, + pconvert_flag_force_type = 0x200, + pconvert_flag_downgrade = 0x400, + pconvert_flag_check_compat = 0x0800, + // Additional + pconvert_flag_hvcurve = 0x4000 +}; + +void SplinesRemoveBetween(SplinePoint *from, SplinePoint *to, int type); + +#endif // _SEEN_SPLINEFONT_H_ diff --git a/src/path/splinefit/splinerefigure.c b/src/path/splinefit/splinerefigure.c new file mode 100644 index 0000000..2f6d09b --- /dev/null +++ b/src/path/splinefit/splinerefigure.c @@ -0,0 +1,117 @@ +// SPDX-License-Identifier: GPL-2.0-or-later +/* Copyright (C) 2000-2012 by George Williams */ +/* + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + + * The name of the author may not be used to endorse or promote products + * derived from this software without specific prior written permission. + + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR IMPLIED + * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO + * EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; + * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, + * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR + * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF + * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#include <stdbool.h> +#include <stdint.h> + +#include "splinerefigure.h" + +#include "splinefont.h" + +#include <math.h> +#include <stdio.h> +#ifdef HAVE_IEEEFP_H +# include <ieeefp.h> /* Solaris defines isnan in ieeefp rather than math.h */ +#endif + +#include <stdbool.h> + +/* The slight errors introduced by the optimizer turn out to have nasty */ +/* side effects. An error on the order of 7e-8 in splines[1].b caused */ +/* the rasterizer to have kaniptions */ +void SplineRefigure3(Spline *spline) { + SplinePoint *from = spline->from, *to = spline->to; + Spline1D *xsp = &spline->splines[0], *ysp = &spline->splines[1]; + Spline old; + + spline->isquadratic = false; + if ( spline->acceptableextrema ) + old = *spline; + xsp->d = from->me.x; ysp->d = from->me.y; + // Set noprevcp and nonextcp based on point values but then make sure both + // have the same value + from->nonextcp = from->nextcp.x==from->me.x && from->nextcp.y == from->me.y; + to->noprevcp = to->prevcp.x==to->me.x && to->prevcp.y == to->me.y; + if ( !from->nonextcp || !to->noprevcp ) + from->nonextcp = to->noprevcp = false; + if ( from->nonextcp && to->noprevcp ) { + spline->islinear = true; + xsp->c = to->me.x-from->me.x; + ysp->c = to->me.y-from->me.y; + xsp->a = xsp->b = 0; + ysp->a = ysp->b = 0; + } else { + /* from p. 393 (Operator Details, curveto) PostScript Lang. Ref. Man. (Red book) */ + xsp->c = 3*(from->nextcp.x-from->me.x); + ysp->c = 3*(from->nextcp.y-from->me.y); + xsp->b = 3*(to->prevcp.x-from->nextcp.x)-xsp->c; + ysp->b = 3*(to->prevcp.y-from->nextcp.y)-ysp->c; + xsp->a = to->me.x-from->me.x-xsp->c-xsp->b; + ysp->a = to->me.y-from->me.y-ysp->c-ysp->b; + if ( RealNear(xsp->c,0)) xsp->c=0; + if ( RealNear(ysp->c,0)) ysp->c=0; + if ( RealNear(xsp->b,0)) xsp->b=0; + if ( RealNear(ysp->b,0)) ysp->b=0; + if ( RealNear(xsp->a,0)) xsp->a=0; + if ( RealNear(ysp->a,0)) ysp->a=0; + if ( xsp->a!=0 && ( Within16RoundingErrors(xsp->a+from->me.x,from->me.x) || + Within16RoundingErrors(xsp->a+to->me.x,to->me.x))) + xsp->a = 0; + if ( ysp->a!=0 && ( Within16RoundingErrors(ysp->a+from->me.y,from->me.y) || + Within16RoundingErrors(ysp->a+to->me.y,to->me.y))) + ysp->a = 0; + SplineIsLinear(spline); + spline->islinear = false; + if ( ysp->a==0 && xsp->a==0 ) { + if ( ysp->b==0 && xsp->b==0 ) + spline->islinear = true; /* This seems extremely unlikely... */ + else + spline->isquadratic = true; /* Only likely if we read in a TTF */ + } + } + if ( !isfinite(ysp->a) || !isfinite(xsp->a) || !isfinite(ysp->c) || !isfinite(xsp->c) || !isfinite(ysp->d) || !isfinite(xsp->d)) + IError("NaN value in spline creation"); + LinearApproxFree(spline->approx); + spline->approx = NULL; + spline->knowncurved = false; + spline->knownlinear = spline->islinear; + SplineIsLinear(spline); + spline->order2 = false; + + if ( spline->acceptableextrema ) { + /* I don't check "d", because changes to that reflect simple */ + /* translations which will not affect the shape of the spline */ + if ( !RealNear(old.splines[0].a,spline->splines[0].a) || + !RealNear(old.splines[0].b,spline->splines[0].b) || + !RealNear(old.splines[0].c,spline->splines[0].c) || + !RealNear(old.splines[1].a,spline->splines[1].a) || + !RealNear(old.splines[1].b,spline->splines[1].b) || + !RealNear(old.splines[1].c,spline->splines[1].c) ) + spline->acceptableextrema = false; + } +} diff --git a/src/path/splinefit/splinerefigure.h b/src/path/splinefit/splinerefigure.h new file mode 100644 index 0000000..110e2f9 --- /dev/null +++ b/src/path/splinefit/splinerefigure.h @@ -0,0 +1,9 @@ +// SPDX-License-Identifier: GPL-2.0-or-later +#ifndef FONTFORGE_SPLINEREFIGURE_H +#define FONTFORGE_SPLINEREFIGURE_H + +#include "splinefont.h" + +extern void SplineRefigure3(Spline *spline); + +#endif /* FONTFORGE_SPLINEREFIGURE_H */ |