Warning, file /office/calligra/libs/flake/KoCurveFit.cpp was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).
0001 /* This file is part of the KDE project 0002 Copyright (C) 2001-2003 Rob Buis <buis@kde.org> 0003 Copyright (C) 2007, 2009 Jan Hambrecht <jaham@gmx.net> 0004 0005 This library is free software; you can redistribute it and/or 0006 modify it under the terms of the GNU Library General Public 0007 License as published by the Free Software Foundation; either 0008 version 2 of the License, or (at your option) any later version. 0009 0010 This library is distributed in the hope that it will be useful, 0011 but WITHOUT ANY WARRANTY; without even the implied warranty of 0012 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0013 Library General Public License for more details. 0014 0015 You should have received a copy of the GNU Library General Public License 0016 along with this library; see the file COPYING.LIB. If not, write to 0017 the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 0018 * Boston, MA 02110-1301, USA. 0019 */ 0020 0021 #include "KoCurveFit.h" 0022 #include <KoPathShape.h> 0023 #include <QVector> 0024 #include <math.h> 0025 0026 /// our equivalent to zero 0027 const qreal Zero = 10e-12; 0028 0029 /* 0030 An Algorithm for Automatically Fitting Digitized Curves 0031 by Philip J. Schneider 0032 from "Graphics Gems", Academic Press, 1990 0033 0034 http://tog.acm.org/resources/GraphicsGems/gems/FitCurves.c 0035 http://tog.acm.org/resources/GraphicsGems/gems/README 0036 */ 0037 0038 class FitVector { 0039 public: 0040 FitVector(const QPointF &p) { 0041 m_X = p.x(); 0042 m_Y = p.y(); 0043 } 0044 0045 FitVector() { 0046 m_X = 0; 0047 m_Y = 0; 0048 } 0049 0050 FitVector(const QPointF &a, const QPointF &b) { 0051 m_X = a.x() - b.x(); 0052 m_Y = a.y() - b.y(); 0053 } 0054 0055 void normalize() { 0056 qreal len = length(); 0057 if (qFuzzyCompare(len, qreal(0.0))) { 0058 return; 0059 } 0060 m_X /= len; m_Y /= len; 0061 } 0062 0063 void negate() { 0064 m_X = -m_X; 0065 m_Y = -m_Y; 0066 } 0067 0068 void scale(qreal s) { 0069 qreal len = length(); 0070 if (qFuzzyCompare(len, qreal(0.0))) { 0071 return; 0072 } 0073 m_X *= s / len; 0074 m_Y *= s / len; 0075 } 0076 0077 qreal dot(const FitVector &v) const { 0078 return ((m_X*v.m_X) + (m_Y*v.m_Y)); 0079 } 0080 0081 qreal length() const { 0082 return (qreal) sqrt(m_X*m_X + m_Y*m_Y); 0083 } 0084 0085 QPointF operator+(const QPointF &p) { 0086 QPointF b(p.x() + m_X, p.y() + m_Y); 0087 return b; 0088 } 0089 0090 public: 0091 qreal m_X, m_Y; 0092 }; 0093 0094 qreal distance(const QPointF &p1, const QPointF &p2) 0095 { 0096 qreal dx = (p1.x() - p2.x()); 0097 qreal dy = (p1.y() - p2.y()); 0098 return sqrt(dx*dx + dy*dy); 0099 } 0100 0101 0102 FitVector ComputeLeftTangent(const QVector<QPointF> &points, int end) 0103 { 0104 FitVector tHat1(points.at(end + 1), points.at(end)); 0105 0106 tHat1.normalize(); 0107 0108 return tHat1; 0109 } 0110 0111 FitVector ComputeRightTangent(const QVector<QPointF> &points, int end) 0112 { 0113 FitVector tHat1(points.at(end - 1), points.at(end)); 0114 0115 tHat1.normalize(); 0116 0117 return tHat1; 0118 } 0119 0120 /* 0121 * ChordLengthParameterize : 0122 * Assign parameter values to digitized points 0123 * using relative distances between points. 0124 */ 0125 static qreal *ChordLengthParameterize(const QVector<QPointF> &points, int first, int last) 0126 { 0127 int i; 0128 qreal *u; /* Parameterization */ 0129 0130 u = new qreal[(last-first+1)]; 0131 0132 u[0] = 0.0; 0133 for (i = first + 1; i <= last; ++i) { 0134 u[i-first] = u[i-first-1] + 0135 distance(points.at(i), points.at(i - 1)); 0136 } 0137 0138 qreal denominator = u[last-first]; 0139 if (qFuzzyCompare(denominator, qreal(0.0))) { 0140 denominator = Zero; 0141 } 0142 0143 for (i = first + 1; i <= last; ++i) { 0144 u[i-first] = u[i-first] / denominator; 0145 } 0146 0147 return(u); 0148 } 0149 0150 static FitVector VectorAdd(FitVector a, FitVector b) 0151 { 0152 FitVector c; 0153 c.m_X = a.m_X + b.m_X; c.m_Y = a.m_Y + b.m_Y; 0154 return (c); 0155 } 0156 static FitVector VectorScale(FitVector v, qreal s) 0157 { 0158 FitVector result; 0159 result.m_X = v.m_X * s; result.m_Y = v.m_Y * s; 0160 return (result); 0161 } 0162 0163 static FitVector VectorSub(FitVector a, FitVector b) 0164 { 0165 FitVector c; 0166 c.m_X = a.m_X - b.m_X; c.m_Y = a.m_Y - b.m_Y; 0167 return (c); 0168 } 0169 0170 static FitVector ComputeCenterTangent(const QVector<QPointF> &points, int center) 0171 { 0172 FitVector V1, V2, tHatCenter; 0173 0174 FitVector cpointb(points.at(center - 1)); 0175 FitVector cpoint(points.at(center)); 0176 FitVector cpointa(points.at(center + 1)); 0177 0178 V1 = VectorSub(cpointb, cpoint); 0179 V2 = VectorSub(cpoint, cpointa); 0180 tHatCenter.m_X = ((V1.m_X + V2.m_X) / 2.0); 0181 tHatCenter.m_Y = ((V1.m_Y + V2.m_Y) / 2.0); 0182 tHatCenter.normalize(); 0183 return tHatCenter; 0184 } 0185 0186 /* 0187 * B0, B1, B2, B3 : 0188 * Bezier multipliers 0189 */ 0190 static qreal B0(qreal u) 0191 { 0192 qreal tmp = 1.0 - u; 0193 return (tmp * tmp * tmp); 0194 } 0195 0196 0197 static qreal B1(qreal u) 0198 { 0199 qreal tmp = 1.0 - u; 0200 return (3 * u *(tmp * tmp)); 0201 } 0202 0203 static qreal B2(qreal u) 0204 { 0205 qreal tmp = 1.0 - u; 0206 return (3 * u * u * tmp); 0207 } 0208 0209 static qreal B3(qreal u) 0210 { 0211 return (u * u * u); 0212 } 0213 0214 /* 0215 * GenerateBezier : 0216 * Use least-squares method to find Bezier control points for region. 0217 * 0218 */ 0219 QPointF* GenerateBezier(const QVector<QPointF> &points, int first, int last, qreal *uPrime, FitVector tHat1, FitVector tHat2) 0220 { 0221 int i; 0222 int nPts; /* Number of pts in sub-curve */ 0223 qreal C[2][2]; /* Matrix C */ 0224 qreal X[2]; /* Matrix X */ 0225 qreal det_C0_C1, /* Determinants of matrices */ 0226 det_C0_X, 0227 det_X_C1; 0228 qreal alpha_l, /* Alpha values, left and right */ 0229 alpha_r; 0230 FitVector tmp; /* Utility variable */ 0231 QPointF *curve; 0232 0233 curve = new QPointF[4]; 0234 nPts = last - first + 1; 0235 0236 /* Precomputed rhs for eqn */ 0237 // FitVector A[nPts][2] 0238 QVector< QVector<FitVector> > A(nPts, QVector<FitVector>(2)); 0239 0240 /* Compute the A's */ 0241 for (i = 0; i < nPts; ++i) { 0242 FitVector v1, v2; 0243 v1 = tHat1; 0244 v2 = tHat2; 0245 v1.scale(B1(uPrime[i])); 0246 v2.scale(B2(uPrime[i])); 0247 A[i][0] = v1; 0248 A[i][1] = v2; 0249 } 0250 0251 /* Create the C and X matrices */ 0252 C[0][0] = 0.0; 0253 C[0][1] = 0.0; 0254 C[1][0] = 0.0; 0255 C[1][1] = 0.0; 0256 X[0] = 0.0; 0257 X[1] = 0.0; 0258 0259 for (i = 0; i < nPts; ++i) { 0260 C[0][0] += (A[i][0]).dot(A[i][0]); 0261 C[0][1] += A[i][0].dot(A[i][1]); 0262 /* C[1][0] += V2Dot(&A[i][0], &A[i][1]);*/ 0263 C[1][0] = C[0][1]; 0264 C[1][1] += A[i][1].dot(A[i][1]); 0265 0266 FitVector vfirstp1(points.at(first + i)); 0267 FitVector vfirst(points.at(first)); 0268 FitVector vlast(points.at(last)); 0269 0270 tmp = VectorSub(vfirstp1, 0271 VectorAdd( 0272 VectorScale(vfirst, B0(uPrime[i])), 0273 VectorAdd( 0274 VectorScale(vfirst, B1(uPrime[i])), 0275 VectorAdd( 0276 VectorScale(vlast, B2(uPrime[i])), 0277 VectorScale(vlast, B3(uPrime[i])))))); 0278 0279 0280 X[0] += A[i][0].dot(tmp); 0281 X[1] += A[i][1].dot(tmp); 0282 } 0283 0284 /* Compute the determinants of C and X */ 0285 det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1]; 0286 det_C0_X = C[0][0] * X[1] - C[0][1] * X[0]; 0287 det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1]; 0288 0289 /* Finally, derive alpha values */ 0290 if (qFuzzyCompare(det_C0_C1, qreal(0.0))) { 0291 det_C0_C1 = (C[0][0] * C[1][1]) * 10e-12; 0292 if (qFuzzyCompare(det_C0_C1, qreal(0.0))) { 0293 det_C0_C1 = Zero; 0294 } 0295 } 0296 alpha_l = det_X_C1 / det_C0_C1; 0297 alpha_r = det_C0_X / det_C0_C1; 0298 0299 0300 /* If alpha negative, use the Wu/Barsky heuristic (see text) */ 0301 /* (if alpha is 0, you get coincident control points that lead to 0302 * divide by zero in any subsequent NewtonRaphsonRootFind() call. */ 0303 if (alpha_l < 1.0e-6 || alpha_r < 1.0e-6) { 0304 qreal dist = distance(points.at(last), points.at(first)) / 3.0; 0305 0306 curve[0] = points.at(first); 0307 curve[3] = points.at(last); 0308 0309 tHat1.scale(dist); 0310 tHat2.scale(dist); 0311 0312 curve[1] = tHat1 + curve[0]; 0313 curve[2] = tHat2 + curve[3]; 0314 return curve; 0315 } 0316 0317 /* First and last control points of the Bezier curve are */ 0318 /* positioned exactly at the first and last data points */ 0319 /* Control points 1 and 2 are positioned an alpha distance out */ 0320 /* on the tangent vectors, left and right, respectively */ 0321 curve[0] = points.at(first); 0322 curve[3] = points.at(last); 0323 0324 tHat1.scale(alpha_l); 0325 tHat2.scale(alpha_r); 0326 0327 curve[1] = tHat1 + curve[0]; 0328 curve[2] = tHat2 + curve[3]; 0329 0330 return (curve); 0331 } 0332 0333 /* 0334 * Bezier : 0335 * Evaluate a Bezier curve at a particular parameter value 0336 * 0337 */ 0338 static QPointF BezierII(int degree, QPointF *V, qreal t) 0339 { 0340 int i, j; 0341 QPointF Q; /* Point on curve at parameter t */ 0342 QPointF *Vtemp; /* Local copy of control points */ 0343 0344 Vtemp = new QPointF[degree+1]; 0345 0346 for (i = 0; i <= degree; ++i) { 0347 Vtemp[i] = V[i]; 0348 } 0349 0350 /* Triangle computation */ 0351 for (i = 1; i <= degree; ++i) { 0352 for (j = 0; j <= degree - i; ++j) { 0353 Vtemp[j].setX((1.0 - t) * Vtemp[j].x() + t * Vtemp[j+1].x()); 0354 Vtemp[j].setY((1.0 - t) * Vtemp[j].y() + t * Vtemp[j+1].y()); 0355 } 0356 } 0357 0358 Q = Vtemp[0]; 0359 delete[] Vtemp; 0360 return Q; 0361 } 0362 0363 /* 0364 * ComputeMaxError : 0365 * Find the maximum squared distance of digitized points 0366 * to fitted curve. 0367 */ 0368 static qreal ComputeMaxError(const QVector<QPointF> &points, int first, int last, QPointF *curve, qreal *u, int *splitPoint) 0369 { 0370 int i; 0371 qreal maxDist; /* Maximum error */ 0372 qreal dist; /* Current error */ 0373 QPointF P; /* Point on curve */ 0374 FitVector v; /* Vector from point to curve */ 0375 0376 *splitPoint = (last - first + 1) / 2; 0377 maxDist = 0.0; 0378 for (i = first + 1; i < last; ++i) { 0379 P = BezierII(3, curve, u[i-first]); 0380 v = VectorSub(P, points.at(i)); 0381 dist = v.length(); 0382 if (dist >= maxDist) { 0383 maxDist = dist; 0384 *splitPoint = i; 0385 } 0386 } 0387 return (maxDist); 0388 } 0389 0390 0391 /* 0392 * NewtonRaphsonRootFind : 0393 * Use Newton-Raphson iteration to find better root. 0394 */ 0395 static qreal NewtonRaphsonRootFind(QPointF *Q, QPointF P, qreal u) 0396 { 0397 qreal numerator, denominator; 0398 QPointF Q1[3], Q2[2]; /* Q' and Q'' */ 0399 QPointF Q_u, Q1_u, Q2_u; /*u evaluated at Q, Q', & Q'' */ 0400 qreal uPrime; /* Improved u */ 0401 int i; 0402 0403 /* Compute Q(u) */ 0404 Q_u = BezierII(3, Q, u); 0405 0406 /* Generate control vertices for Q' */ 0407 for (i = 0; i <= 2; ++i) { 0408 Q1[i].setX((Q[i+1].x() - Q[i].x()) * 3.0); 0409 Q1[i].setY((Q[i+1].y() - Q[i].y()) * 3.0); 0410 } 0411 0412 /* Generate control vertices for Q'' */ 0413 for (i = 0; i <= 1; ++i) { 0414 Q2[i].setX((Q1[i+1].x() - Q1[i].x()) * 2.0); 0415 Q2[i].setY((Q1[i+1].y() - Q1[i].y()) * 2.0); 0416 } 0417 0418 /* Compute Q'(u) and Q''(u) */ 0419 Q1_u = BezierII(2, Q1, u); 0420 Q2_u = BezierII(1, Q2, u); 0421 0422 /* Compute f(u)/f'(u) */ 0423 numerator = (Q_u.x() - P.x()) * (Q1_u.x()) + (Q_u.y() - P.y()) * (Q1_u.y()); 0424 denominator = (Q1_u.x()) * (Q1_u.x()) + (Q1_u.y()) * (Q1_u.y()) + 0425 (Q_u.x() - P.x()) * (Q2_u.x()) + (Q_u.y() - P.y()) * (Q2_u.y()); 0426 0427 if (qFuzzyCompare(denominator, qreal(0.0))) { 0428 denominator = Zero; 0429 } 0430 0431 /* u = u - f(u)/f'(u) */ 0432 uPrime = u - (numerator / denominator); 0433 return (uPrime); 0434 } 0435 0436 /* 0437 * Reparameterize: 0438 * Given set of points and their parameterization, try to find 0439 * a better parameterization. 0440 * 0441 */ 0442 static qreal *Reparameterize(const QVector<QPointF> &points, int first, int last, qreal *u, QPointF *curve) 0443 { 0444 int nPts = last - first + 1; 0445 int i; 0446 qreal *uPrime; /* New parameter values */ 0447 0448 uPrime = new qreal[nPts]; 0449 for (i = first; i <= last; ++i) { 0450 uPrime[i-first] = NewtonRaphsonRootFind(curve, points.at(i), u[i-first]); 0451 } 0452 return (uPrime); 0453 } 0454 0455 QPointF *FitCubic(const QVector<QPointF> &points, int first, int last, FitVector tHat1, FitVector tHat2, float error, int &width) 0456 { 0457 qreal *u; 0458 qreal *uPrime; 0459 qreal maxError; 0460 int splitPoint; 0461 int nPts; 0462 qreal iterationError; 0463 int maxIterations = 4; 0464 FitVector tHatCenter; 0465 QPointF *curve; 0466 int i; 0467 0468 width = 0; 0469 0470 0471 iterationError = error * error; 0472 nPts = last - first + 1; 0473 0474 if (nPts == 2) { 0475 qreal dist = distance(points.at(last), points.at(first)) / 3.0; 0476 0477 curve = new QPointF[4]; 0478 0479 curve[0] = points.at(first); 0480 curve[3] = points.at(last); 0481 0482 tHat1.scale(dist); 0483 tHat2.scale(dist); 0484 0485 curve[1] = tHat1 + curve[0]; 0486 curve[2] = tHat2 + curve[3]; 0487 0488 width = 4; 0489 return curve; 0490 } 0491 0492 /* Parameterize points, and attempt to fit curve */ 0493 u = ChordLengthParameterize(points, first, last); 0494 curve = GenerateBezier(points, first, last, u, tHat1, tHat2); 0495 0496 0497 /* Find max deviation of points to fitted curve */ 0498 maxError = ComputeMaxError(points, first, last, curve, u, &splitPoint); 0499 if (maxError < error) { 0500 delete[] u; 0501 width = 4; 0502 return curve; 0503 } 0504 0505 0506 /* If error not too large, try some reparameterization */ 0507 /* and iteration */ 0508 if (maxError < iterationError) { 0509 for (i = 0; i < maxIterations; ++i) { 0510 uPrime = Reparameterize(points, first, last, u, curve); 0511 delete[] curve; 0512 curve = GenerateBezier(points, first, last, uPrime, tHat1, tHat2); 0513 maxError = ComputeMaxError(points, first, last, 0514 curve, uPrime, &splitPoint); 0515 if (maxError < error) { 0516 delete[] u; 0517 delete[] uPrime; 0518 width = 4; 0519 return curve; 0520 } 0521 delete[] u; 0522 u = uPrime; 0523 } 0524 } 0525 0526 /* Fitting failed -- split at max error point and fit recursively */ 0527 delete[] u; 0528 delete[] curve; 0529 tHatCenter = ComputeCenterTangent(points, splitPoint); 0530 0531 int w1, w2; 0532 QPointF *cu1 = 0, *cu2 = 0; 0533 cu1 = FitCubic(points, first, splitPoint, tHat1, tHatCenter, error, w1); 0534 0535 tHatCenter.negate(); 0536 cu2 = FitCubic(points, splitPoint, last, tHatCenter, tHat2, error, w2); 0537 0538 QPointF *newcurve = new QPointF[w1+w2]; 0539 for (int i = 0; i < w1; ++i) { 0540 newcurve[i] = cu1[i]; 0541 } 0542 for (int i = 0; i < w2; ++i) { 0543 newcurve[i+w1] = cu2[i]; 0544 } 0545 0546 delete[] cu1; 0547 delete[] cu2; 0548 width = w1 + w2; 0549 return newcurve; 0550 } 0551 0552 0553 KoPathShape * bezierFit(const QVector<QPointF> &points, float error) 0554 { 0555 FitVector tHat1, tHat2; 0556 0557 tHat1 = ComputeLeftTangent(points, 0); 0558 tHat2 = ComputeRightTangent(points, points.count() - 1); 0559 0560 int width = 0; 0561 QPointF *curve; 0562 curve = FitCubic(points, 0, points.count() - 1, tHat1, tHat2, error, width); 0563 0564 KoPathShape * path = new KoPathShape(); 0565 0566 if (width > 3) { 0567 path->moveTo(curve[0]); 0568 path->curveTo(curve[1], curve[2], curve[3]); 0569 for (int i = 4; i < width; i += 4) { 0570 path->curveTo(curve[i+1], curve[i+2], curve[i+3]); 0571 } 0572 } 0573 0574 delete[] curve; 0575 return path; 0576 } 0577