From 8b78a8791bc539bcffe7159f9d9714d577cb3d7d Mon Sep 17 00:00:00 2001 From: Michele Calgaro Date: Sun, 23 May 2021 20:48:35 +0900 Subject: Renaming of files in preparation for code style tools. Signed-off-by: Michele Calgaro --- karbon/tools/vcurvefit.cpp | 565 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 565 insertions(+) create mode 100644 karbon/tools/vcurvefit.cpp (limited to 'karbon/tools/vcurvefit.cpp') diff --git a/karbon/tools/vcurvefit.cpp b/karbon/tools/vcurvefit.cpp new file mode 100644 index 00000000..40d7b7e3 --- /dev/null +++ b/karbon/tools/vcurvefit.cpp @@ -0,0 +1,565 @@ +/* This file is part of the KDE project + Copyright (C) 2001, 2002, 2003 The Karbon Developers + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public + License as published by the Free Software Foundation; either + version 2 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public License + along with this library; see the file COPYING.LIB. If not, write to + the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, + * Boston, MA 02110-1301, USA. +*/ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/* + An Algorithm for Automatically Fitting Digitized Curves + by Philip J. Schneider + from "Graphics Gems", Academic Press, 1990 + + http://www.acm.org/pubs/tog/GraphicsGems/gems/FitCurves.c + http://www.acm.org/pubs/tog/GraphicsGems/gems/README +*/ + + +#include "vcurvefit.h" + +#define MAXPOINTS 1000 /* The most points you can have */ + + +class FitVector { + public: + FitVector(KoPoint &p){ + m_X=p.x(); + m_Y=p.y(); + } + + FitVector(){ + m_X=0; + m_Y=0; + } + + FitVector(KoPoint &a,KoPoint &b){ + m_X=a.x()-b.x(); + m_Y=a.y()-b.y(); + } + + void normalize(){ + double len=length(); + if(len==0.0f) + return; + m_X/=len; m_Y/=len; + } + + void negate(){ + m_X = -m_X; + m_Y = -m_Y; + } + + void scale(double s){ + double len = length(); + if(len==0.0f) + return; + m_X *= s/len; + m_Y *= s/len; + } + + double dot(FitVector &v){ + return ((m_X*v.m_X)+(m_Y*v.m_Y)); + } + + double length(){ + return (double) sqrt(m_X*m_X+m_Y*m_Y); + } + + KoPoint operator+(KoPoint &p){ + KoPoint b(p.x()+m_X,p.y()+m_Y); + return b; + } + + public: + double m_X,m_Y; +}; + +double distance(KoPoint *p1,KoPoint *p2){ + double dx = (p1->x()-p2->x()); + double dy = (p1->y()-p2->y()); + return sqrt( dx*dx + dy*dy ); +} + + +FitVector ComputeLeftTangent(TQPtrList &points,int end){ + FitVector tHat1(*points.at(end+1),*points.at(end)); + + tHat1.normalize(); + + return tHat1; +} + +FitVector ComputeRightTangent(TQPtrList &points,int end){ + FitVector tHat1(*points.at(end-1),*points.at(end)); + + tHat1.normalize(); + + return tHat1; +} + +/* + * ChordLengthParameterize : + * Assign parameter values to digitized points + * using relative distances between points. + */ +static double *ChordLengthParameterize(TQPtrList points,int first,int last) +{ + int i; + double *u; /* Parameterization */ + + u = new double[(last-first+1)]; + + u[0] = 0.0; + for (i = first+1; i <= last; i++) { + u[i-first] = u[i-first-1] + + distance(points.at(i), points.at(i-1)); + } + + for (i = first + 1; i <= last; i++) { + u[i-first] = u[i-first] / u[last-first]; + } + + return(u); +} + +static FitVector VectorAdd(FitVector a,FitVector b) +{ + FitVector c; + c.m_X = a.m_X + b.m_X; c.m_Y = a.m_Y + b.m_Y; + return (c); +} +static FitVector VectorScale(FitVector v,double s) +{ + FitVector result; + result.m_X = v.m_X * s; result.m_Y = v.m_Y * s; + return (result); +} + +static FitVector VectorSub(FitVector a,FitVector b) +{ + FitVector c; + c.m_X = a.m_X - b.m_X; c.m_Y = a.m_Y - b.m_Y; + return (c); +} + +static FitVector ComputeCenterTangent(TQPtrList points,int center) +{ + FitVector V1, V2, tHatCenter; + + FitVector cpointb = *points.at(center-1); + FitVector cpoint = *points.at(center); + FitVector cpointa = *points.at(center+1); + + V1 = VectorSub(cpointb,cpoint); + V2 = VectorSub(cpoint,cpointa); + tHatCenter.m_X= ((V1.m_X + V2.m_X)/2.0); + tHatCenter.m_Y= ((V1.m_Y + V2.m_Y)/2.0); + tHatCenter.normalize(); + return tHatCenter; +} + +/* + * B0, B1, B2, B3 : + * Bezier multipliers + */ +static double B0(double u) +{ + double tmp = 1.0 - u; + return (tmp * tmp * tmp); +} + + +static double B1(double u) +{ + double tmp = 1.0 - u; + return (3 * u * (tmp * tmp)); +} + +static double B2(double u) +{ + double tmp = 1.0 - u; + return (3 * u * u * tmp); +} + +static double B3(double u) +{ + return (u * u * u); +} + +/* + * GenerateBezier : + * Use least-squares method to find Bezier control points for region. + * + */ +KoPoint* GenerateBezier(TQPtrList &points, int first, int last, double *uPrime,FitVector tHat1,FitVector tHat2) +{ + int i; + FitVector A[MAXPOINTS][2]; /* Precomputed rhs for eqn */ + int nPts; /* Number of pts in sub-curve */ + double C[2][2]; /* Matrix C */ + double X[2]; /* Matrix X */ + double det_C0_C1, /* Determinants of matrices */ + det_C0_X, + det_X_C1; + double alpha_l, /* Alpha values, left and right */ + alpha_r; + FitVector tmp; /* Utility variable */ + KoPoint *curve; + + curve = new KoPoint[4]; + nPts = last - first + 1; + + + /* Compute the A's */ + for (i = 0; i < nPts; i++) { + FitVector v1, v2; + v1 = tHat1; + v2 = tHat2; + v1.scale(B1(uPrime[i])); + v2.scale(B2(uPrime[i])); + A[i][0] = v1; + A[i][1] = v2; + } + + /* Create the C and X matrices */ + C[0][0] = 0.0; + C[0][1] = 0.0; + C[1][0] = 0.0; + C[1][1] = 0.0; + X[0] = 0.0; + X[1] = 0.0; + + for (i = 0; i < nPts; i++) { + C[0][0] += (A[i][0]).dot(A[i][0]); + C[0][1] += A[i][0].dot(A[i][1]); + /* C[1][0] += V2Dot(&A[i][0], &A[i][1]);*/ + C[1][0] = C[0][1]; + C[1][1] += A[i][1].dot(A[i][1]); + + FitVector vfirstp1(*points.at(first+i)); + FitVector vfirst(*points.at(first)); + FitVector vlast(*points.at(last)); + + tmp = VectorSub(vfirstp1, + VectorAdd( + VectorScale(vfirst, B0(uPrime[i])), + VectorAdd( + VectorScale(vfirst, B1(uPrime[i])), + VectorAdd( + VectorScale(vlast, B2(uPrime[i])), + VectorScale(vlast, B3(uPrime[i])) )))); + + + X[0] += A[i][0].dot(tmp); + X[1] += A[i][1].dot(tmp); + } + + /* Compute the determinants of C and X */ + det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1]; + det_C0_X = C[0][0] * X[1] - C[0][1] * X[0]; + det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1]; + + /* Finally, derive alpha values */ + if (det_C0_C1 == 0.0) { + det_C0_C1 = (C[0][0] * C[1][1]) * 10e-12; + } + alpha_l = det_X_C1 / det_C0_C1; + alpha_r = det_C0_X / det_C0_C1; + + + /* If alpha negative, use the Wu/Barsky heuristic (see text) */ + /* (if alpha is 0, you get coincident control points that lead to + * divide by zero in any subsequent NewtonRaphsonRootFind() call. */ + if (alpha_l < 1.0e-6 || alpha_r < 1.0e-6) { + double dist = distance(points.at(last),points.at(first)) / + 3.0; + + curve[0] = *points.at(first); + curve[3] = *points.at(last); + + tHat1.scale(dist); + tHat2.scale(dist); + + curve[1] = tHat1 + curve[0]; + curve[2] = tHat2 + curve[3]; + return curve; + } + + /* First and last control points of the Bezier curve are */ + /* positioned exactly at the first and last data points */ + /* Control points 1 and 2 are positioned an alpha distance out */ + /* on the tangent vectors, left and right, respectively */ + curve[0] = *points.at(first); + curve[3] = *points.at(last); + + tHat1.scale(alpha_l); + tHat2.scale(alpha_r); + + curve[1] = tHat1 + curve[0]; + curve[2] = tHat2 + curve[3]; + + return (curve); +} + +/* + * Bezier : + * Evaluate a Bezier curve at a particular parameter value + * + */ +static KoPoint BezierII(int degree,KoPoint *V, double t) +{ + int i, j; + KoPoint Q; /* Point on curve at parameter t */ + KoPoint *Vtemp; /* Local copy of control points */ + + Vtemp = new KoPoint[degree+1]; + + for (i = 0; i <= degree; i++) { + Vtemp[i] = V[i]; + } + + /* Triangle computation */ + for (i = 1; i <= degree; i++) { + for (j = 0; j <= degree-i; j++) { + Vtemp[j].setX((1.0 - t) * Vtemp[j].x() + t * Vtemp[j+1].x()); + Vtemp[j].setY((1.0 - t) * Vtemp[j].y() + t * Vtemp[j+1].y()); + } + } + + Q = Vtemp[0]; + delete[] Vtemp; + return Q; +} + +/* + * ComputeMaxError : + * Find the maximum squared distance of digitized points + * to fitted curve. +*/ +static double ComputeMaxError(TQPtrList points,int first,int last,KoPoint *curve,double *u,int *splitPoint) +{ + int i; + double maxDist; /* Maximum error */ + double dist; /* Current error */ + KoPoint P; /* Point on curve */ + FitVector v; /* Vector from point to curve */ + + *splitPoint = (last - first + 1)/2; + maxDist = 0.0; + for (i = first + 1; i < last; i++) { + P = BezierII(3, curve, u[i-first]); + v = VectorSub(P, *points.at(i)); + dist = v.length(); + if (dist >= maxDist) { + maxDist = dist; + *splitPoint = i; + } + } + return (maxDist); +} + + +/* + * NewtonRaphsonRootFind : + * Use Newton-Raphson iteration to find better root. + */ +static double NewtonRaphsonRootFind(KoPoint *Q,KoPoint P,double u) +{ + double numerator, denominator; + KoPoint Q1[3], Q2[2]; /* Q' and Q'' */ + KoPoint TQ_u, Q1_u, Q2_u; /*u evaluated at Q, Q', & Q'' */ + double uPrime; /* Improved u */ + int i; + + /* Compute Q(u) */ + TQ_u = BezierII(3,Q, u); + + /* Generate control vertices for Q' */ + for (i = 0; i <= 2; i++) { + Q1[i].setX((Q[i+1].x() - Q[i].x()) * 3.0); + Q1[i].setY((Q[i+1].y() - Q[i].y()) * 3.0); + } + + /* Generate control vertices for Q'' */ + for (i = 0; i <= 1; i++) { + Q2[i].setX((Q1[i+1].x() - Q1[i].x()) * 2.0); + Q2[i].setY((Q1[i+1].y() - Q1[i].y()) * 2.0); + } + + /* Compute Q'(u) and Q''(u) */ + Q1_u = BezierII(2, Q1, u); + Q2_u = BezierII(1, Q2, u); + + /* Compute f(u)/f'(u) */ + numerator = (TQ_u.x() - P.x()) * (Q1_u.x()) + (TQ_u.y() - P.y()) * (Q1_u.y()); + denominator = (Q1_u.x()) * (Q1_u.x()) + (Q1_u.y()) * (Q1_u.y()) + + (TQ_u.x() - P.x()) * (Q2_u.x()) + (TQ_u.y() - P.y()) * (Q2_u.y()); + + /* u = u - f(u)/f'(u) */ + uPrime = u - (numerator/denominator); + return (uPrime); +} + +/* + * Reparameterize: + * Given set of points and their parameterization, try to find + * a better parameterization. + * + */ +static double *Reparameterize(TQPtrList points,int first,int last,double *u,KoPoint *curve) +{ + int nPts = last-first+1; + int i; + double *uPrime; /* New parameter values */ + + uPrime = new double[nPts]; + for (i = first; i <= last; i++) { + uPrime[i-first] = NewtonRaphsonRootFind(curve, *points.at(i), u[i- + first]); + } + return (uPrime); +} + +KoPoint *FitCubic(TQPtrList &points,int first,int last,FitVector tHat1,FitVector tHat2,float error,int &width){ + double *u; + double *uPrime; + double maxError; + int splitPoint; + int nPts; + double iterationError; + int maxIterations=4; + FitVector tHatCenter; + KoPoint *curve; + int i; + + width=0; + + + iterationError=error*error; + nPts = last-first+1; + + if(nPts == 2){ + double dist = distance(points.at(last), points.at(first)) / 3.0; + + curve = new KoPoint[4]; + + curve[0] = *points.at(first); + curve[3] = *points.at(last); + + tHat1.scale(dist); + tHat2.scale(dist); + + curve[1] = tHat1 + curve[0]; + curve[2] = tHat2 + curve[3]; + + width=4; + return curve; + } + + /* Parameterize points, and attempt to fit curve */ + u = ChordLengthParameterize(points, first, last); + curve = GenerateBezier(points, first, last, u, tHat1, tHat2); + + + /* Find max deviation of points to fitted curve */ + maxError = ComputeMaxError(points, first, last, curve, u, &splitPoint); + if (maxError < error) { + delete[] u; + width=4; + return curve; + } + + + /* If error not too large, try some reparameterization */ + /* and iteration */ + if (maxError < iterationError) { + for (i = 0; i < maxIterations; i++) { + uPrime = Reparameterize(points, first, last, u, curve); + curve = GenerateBezier(points, first, last, uPrime, tHat1, tHat2); + maxError = ComputeMaxError(points, first, last, + curve, uPrime, &splitPoint); + if (maxError < error) { + delete[] u; + width=4; + return curve; + } + delete[] u; + u = uPrime; + } + } + + /* Fitting failed -- split at max error point and fit recursively */ + delete[] u; + delete[] curve; + tHatCenter = ComputeCenterTangent(points, splitPoint); + + int w1,w2; + KoPoint *cu1=NULL, *cu2=NULL; + cu1 = FitCubic(points, first, splitPoint, tHat1, tHatCenter, error,w1); + + tHatCenter.negate(); + cu2 = FitCubic(points, splitPoint, last, tHatCenter, tHat2, error,w2); + + KoPoint *newcurve = new KoPoint[w1+w2]; + for(int i=0;i &points,float error){ + FitVector tHat1, tHat2; + + tHat1 = ComputeLeftTangent(points,0); + tHat2 = ComputeRightTangent(points,points.count()-1); + + int width=0; + KoPoint *curve; + curve = FitCubic(points,0,points.count()-1,tHat1,tHat2,error,width); + + VPath *path = new VPath(NULL); + + if(width>3){ + path->moveTo(curve[0]); + path->curveTo(curve[1],curve[2],curve[3]); + for(int i=4;icurveTo(curve[i+1],curve[i+2],curve[i+3]); + } + } + + + delete[] curve; + return path; +} + -- cgit v1.2.1