karbon

vcurvefit.cc

00001 /* This file is part of the KDE project
00002    Copyright (C) 2001, 2002, 2003 The Karbon Developers
00003 
00004    This library is free software; you can redistribute it and/or
00005    modify it under the terms of the GNU Library General Public
00006    License as published by the Free Software Foundation; either
00007    version 2 of the License, or (at your option) any later version.
00008 
00009    This library is distributed in the hope that it will be useful,
00010    but WITHOUT ANY WARRANTY; without even the implied warranty of
00011    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012    Library General Public License for more details.
00013 
00014    You should have received a copy of the GNU Library General Public License
00015    along with this library; see the file COPYING.LIB.  If not, write to
00016    the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
00017  * Boston, MA 02110-1301, USA.
00018 */
00019 
00020 #include <karbon_part.h>
00021 #include <karbon_view.h>
00022 #include <core/vcolor.h>
00023 #include <core/vcomposite.h>
00024 #include <core/vfill.h>
00025 #include <core/vstroke.h>
00026 #include <core/vglobal.h>
00027 #include <render/vpainter.h>
00028 #include <render/vpainterfactory.h>
00029 #include <commands/vshapecmd.h>
00030 
00031 /*
00032     An Algorithm for Automatically Fitting Digitized Curves
00033     by Philip J. Schneider
00034     from "Graphics Gems", Academic Press, 1990
00035 
00036     http://www.acm.org/pubs/tog/GraphicsGems/gems/FitCurves.c
00037     http://www.acm.org/pubs/tog/GraphicsGems/gems/README
00038 */
00039 
00040 
00041 #include "vcurvefit.h"
00042 
00043 #define MAXPOINTS   1000        /* The most points you can have */
00044 
00045 
00046 class FitVector {
00047     public:
00048     FitVector(KoPoint &p){
00049         m_X=p.x();
00050         m_Y=p.y();
00051     }
00052     
00053     FitVector(){
00054         m_X=0;
00055         m_Y=0;
00056     }
00057 
00058     FitVector(KoPoint &a,KoPoint &b){
00059         m_X=a.x()-b.x();
00060         m_Y=a.y()-b.y();
00061     }
00062 
00063     void normalize(){
00064         double len=length();
00065         if(len==0.0f)
00066             return;
00067         m_X/=len; m_Y/=len;
00068     }
00069 
00070     void negate(){
00071         m_X = -m_X;
00072         m_Y = -m_Y;
00073     }
00074 
00075     void scale(double s){
00076         double len = length();
00077         if(len==0.0f)
00078             return;
00079         m_X *= s/len;
00080         m_Y *= s/len;
00081     }
00082 
00083     double dot(FitVector &v){
00084         return ((m_X*v.m_X)+(m_Y*v.m_Y));
00085     }
00086 
00087     double length(){
00088         return (double) sqrt(m_X*m_X+m_Y*m_Y); 
00089     }
00090 
00091     KoPoint operator+(KoPoint &p){
00092         KoPoint b(p.x()+m_X,p.y()+m_Y);
00093         return b;
00094     }
00095 
00096     public:
00097         double m_X,m_Y;
00098 };
00099 
00100 double distance(KoPoint *p1,KoPoint *p2){
00101     double dx = (p1->x()-p2->x());
00102     double dy = (p1->y()-p2->y());
00103     return sqrt( dx*dx + dy*dy );
00104 }
00105 
00106 
00107 FitVector ComputeLeftTangent(QPtrList<KoPoint> &points,int end){
00108     FitVector tHat1(*points.at(end+1),*points.at(end));
00109 
00110     tHat1.normalize();
00111 
00112     return tHat1;
00113 }
00114 
00115 FitVector ComputeRightTangent(QPtrList<KoPoint> &points,int end){
00116     FitVector tHat1(*points.at(end-1),*points.at(end));
00117 
00118     tHat1.normalize();
00119 
00120     return tHat1;
00121 }
00122 
00123 /*
00124  *  ChordLengthParameterize :
00125  *  Assign parameter values to digitized points 
00126  *  using relative distances between points.
00127  */
00128 static double *ChordLengthParameterize(QPtrList<KoPoint> points,int first,int last)
00129 {
00130     int     i;  
00131     double  *u;         /*  Parameterization        */
00132 
00133     u = new double[(last-first+1)];
00134 
00135     u[0] = 0.0;
00136     for (i = first+1; i <= last; i++) {
00137         u[i-first] = u[i-first-1] +
00138                 distance(points.at(i), points.at(i-1));
00139     }
00140 
00141     for (i = first + 1; i <= last; i++) {
00142         u[i-first] = u[i-first] / u[last-first];
00143     }
00144 
00145     return(u);
00146 }
00147 
00148 static FitVector VectorAdd(FitVector a,FitVector b)
00149 {
00150     FitVector   c;
00151     c.m_X = a.m_X + b.m_X;  c.m_Y = a.m_Y + b.m_Y;
00152     return (c);
00153 }
00154 static FitVector VectorScale(FitVector v,double s)
00155 {
00156     FitVector result;
00157     result.m_X = v.m_X * s; result.m_Y = v.m_Y * s;
00158     return (result);
00159 }
00160 
00161 static FitVector VectorSub(FitVector a,FitVector b)
00162 {
00163     FitVector   c;
00164     c.m_X = a.m_X - b.m_X; c.m_Y = a.m_Y - b.m_Y;
00165     return (c);
00166 }
00167 
00168 static FitVector ComputeCenterTangent(QPtrList<KoPoint> points,int center)
00169 {
00170     FitVector V1, V2, tHatCenter;
00171     
00172     FitVector cpointb = *points.at(center-1);
00173     FitVector cpoint = *points.at(center);
00174     FitVector cpointa = *points.at(center+1);
00175 
00176     V1 = VectorSub(cpointb,cpoint);
00177     V2 = VectorSub(cpoint,cpointa);
00178     tHatCenter.m_X= ((V1.m_X + V2.m_X)/2.0);
00179     tHatCenter.m_Y= ((V1.m_Y + V2.m_Y)/2.0);
00180     tHatCenter.normalize();
00181     return tHatCenter;
00182 }
00183 
00184 /*
00185  *  B0, B1, B2, B3 :
00186  *  Bezier multipliers
00187  */
00188 static double B0(double u)
00189 {
00190     double tmp = 1.0 - u;
00191     return (tmp * tmp * tmp);
00192 }
00193 
00194 
00195 static double B1(double u)
00196 {
00197     double tmp = 1.0 - u;
00198     return (3 * u * (tmp * tmp));
00199 }
00200 
00201 static double B2(double u)
00202 {
00203     double tmp = 1.0 - u;
00204     return (3 * u * u * tmp);
00205 }
00206 
00207 static double B3(double u)
00208 {
00209     return (u * u * u);
00210 }
00211 
00212 /*
00213  *  GenerateBezier :
00214  *  Use least-squares method to find Bezier control points for region.
00215  *
00216  */
00217 KoPoint* GenerateBezier(QPtrList<KoPoint> &points, int first, int last, double *uPrime,FitVector tHat1,FitVector tHat2)
00218 {
00219     int     i;
00220     FitVector   A[MAXPOINTS][2];    /* Precomputed rhs for eqn  */
00221     int     nPts;           /* Number of pts in sub-curve */
00222     double  C[2][2];            /* Matrix C     */
00223     double  X[2];           /* Matrix X         */
00224     double  det_C0_C1,      /* Determinants of matrices */
00225             det_C0_X,
00226             det_X_C1;
00227     double  alpha_l,        /* Alpha values, left and right */
00228             alpha_r;
00229     FitVector   tmp;            /* Utility variable     */
00230     KoPoint *curve;
00231     
00232     curve = new KoPoint[4];
00233     nPts = last - first + 1;
00234 
00235  
00236     /* Compute the A's  */
00237     for (i = 0; i < nPts; i++) {
00238         FitVector v1, v2;
00239         v1 = tHat1;
00240         v2 = tHat2;
00241         v1.scale(B1(uPrime[i]));
00242         v2.scale(B2(uPrime[i]));
00243         A[i][0] = v1;
00244         A[i][1] = v2;
00245     }
00246 
00247     /* Create the C and X matrices  */
00248     C[0][0] = 0.0;
00249     C[0][1] = 0.0;
00250     C[1][0] = 0.0;
00251     C[1][1] = 0.0;
00252     X[0]    = 0.0;
00253     X[1]    = 0.0;
00254 
00255     for (i = 0; i < nPts; i++) {
00256         C[0][0] += (A[i][0]).dot(A[i][0]);
00257         C[0][1] += A[i][0].dot(A[i][1]);
00258         /* C[1][0] += V2Dot(&A[i][0], &A[i][1]);*/  
00259         C[1][0] = C[0][1];
00260         C[1][1] += A[i][1].dot(A[i][1]);
00261 
00262         FitVector vfirstp1(*points.at(first+i));
00263         FitVector vfirst(*points.at(first));
00264         FitVector vlast(*points.at(last));
00265 
00266         tmp = VectorSub(vfirstp1,
00267             VectorAdd(
00268               VectorScale(vfirst, B0(uPrime[i])),
00269                 VectorAdd(
00270                     VectorScale(vfirst, B1(uPrime[i])),
00271                             VectorAdd(
00272                             VectorScale(vlast, B2(uPrime[i])),
00273                                 VectorScale(vlast, B3(uPrime[i])) ))));
00274     
00275 
00276     X[0] += A[i][0].dot(tmp);
00277     X[1] += A[i][1].dot(tmp);
00278     }
00279 
00280     /* Compute the determinants of C and X  */
00281     det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
00282     det_C0_X  = C[0][0] * X[1]    - C[0][1] * X[0];
00283     det_X_C1  = X[0]    * C[1][1] - X[1]    * C[0][1];
00284 
00285     /* Finally, derive alpha values */
00286     if (det_C0_C1 == 0.0) {
00287         det_C0_C1 = (C[0][0] * C[1][1]) * 10e-12;
00288     }
00289     alpha_l = det_X_C1 / det_C0_C1;
00290     alpha_r = det_C0_X / det_C0_C1;
00291 
00292 
00293     /*  If alpha negative, use the Wu/Barsky heuristic (see text) */
00294     /* (if alpha is 0, you get coincident control points that lead to
00295      * divide by zero in any subsequent NewtonRaphsonRootFind() call. */
00296     if (alpha_l < 1.0e-6 || alpha_r < 1.0e-6) {
00297         double  dist = distance(points.at(last),points.at(first)) /
00298                     3.0;
00299 
00300         curve[0] = *points.at(first);
00301         curve[3] = *points.at(last);
00302 
00303         tHat1.scale(dist);
00304         tHat2.scale(dist);
00305 
00306         curve[1] = tHat1 + curve[0];
00307         curve[2] = tHat2 + curve[3];
00308             return curve;
00309     }
00310 
00311     /*  First and last control points of the Bezier curve are */
00312     /*  positioned exactly at the first and last data points */
00313     /*  Control points 1 and 2 are positioned an alpha distance out */
00314     /*  on the tangent vectors, left and right, respectively */
00315     curve[0] = *points.at(first);
00316     curve[3] = *points.at(last);
00317 
00318     tHat1.scale(alpha_l);
00319     tHat2.scale(alpha_r);
00320 
00321     curve[1] = tHat1 + curve[0];
00322     curve[2] = tHat2 + curve[3];
00323     
00324     return (curve);
00325 }
00326 
00327 /*
00328  *  Bezier :
00329  *      Evaluate a Bezier curve at a particular parameter value
00330  * 
00331  */
00332 static KoPoint BezierII(int degree,KoPoint *V, double t)
00333 {
00334     int     i, j;       
00335     KoPoint     Q;          /* Point on curve at parameter t    */
00336     KoPoint     *Vtemp;     /* Local copy of control points     */
00337 
00338     Vtemp = new KoPoint[degree+1];
00339     
00340     for (i = 0; i <= degree; i++) {
00341         Vtemp[i] = V[i];
00342     }
00343 
00344     /* Triangle computation */
00345     for (i = 1; i <= degree; i++) { 
00346         for (j = 0; j <= degree-i; j++) {
00347             Vtemp[j].setX((1.0 - t) * Vtemp[j].x() + t * Vtemp[j+1].x());
00348             Vtemp[j].setY((1.0 - t) * Vtemp[j].y() + t * Vtemp[j+1].y());
00349         }
00350     }
00351 
00352     Q = Vtemp[0];
00353     delete[] Vtemp;
00354     return Q;
00355 }
00356 
00357 /*
00358  *  ComputeMaxError :
00359  *  Find the maximum squared distance of digitized points
00360  *  to fitted curve.
00361 */
00362 static double ComputeMaxError(QPtrList<KoPoint> points,int first,int last,KoPoint *curve,double *u,int *splitPoint)
00363 {
00364     int     i;
00365     double  maxDist;        /*  Maximum error       */
00366     double  dist;       /*  Current error       */
00367     KoPoint P;          /*  Point on curve      */
00368     FitVector   v;          /*  Vector from point to curve  */
00369 
00370     *splitPoint = (last - first + 1)/2;
00371     maxDist = 0.0;
00372     for (i = first + 1; i < last; i++) {
00373         P = BezierII(3, curve, u[i-first]);
00374         v = VectorSub(P, *points.at(i));
00375         dist = v.length();
00376         if (dist >= maxDist) {
00377             maxDist = dist;
00378             *splitPoint = i;
00379         }
00380     }
00381     return (maxDist);
00382 }
00383 
00384 
00385 /*
00386  *  NewtonRaphsonRootFind :
00387  *  Use Newton-Raphson iteration to find better root.
00388  */
00389 static double NewtonRaphsonRootFind(KoPoint *Q,KoPoint P,double u)
00390 {
00391     double      numerator, denominator;
00392     KoPoint         Q1[3], Q2[2];   /*  Q' and Q''          */
00393     KoPoint     Q_u, Q1_u, Q2_u; /*u evaluated at Q, Q', & Q''  */
00394     double      uPrime;     /*  Improved u          */
00395     int         i;
00396     
00397     /* Compute Q(u) */
00398     Q_u = BezierII(3,Q, u);
00399     
00400     /* Generate control vertices for Q' */
00401     for (i = 0; i <= 2; i++) {
00402         Q1[i].setX((Q[i+1].x() - Q[i].x()) * 3.0);
00403         Q1[i].setY((Q[i+1].y() - Q[i].y()) * 3.0);
00404     }
00405     
00406     /* Generate control vertices for Q'' */
00407     for (i = 0; i <= 1; i++) {
00408         Q2[i].setX((Q1[i+1].x() - Q1[i].x()) * 2.0);
00409         Q2[i].setY((Q1[i+1].y() - Q1[i].y()) * 2.0);
00410     }
00411     
00412     /* Compute Q'(u) and Q''(u) */
00413     Q1_u = BezierII(2, Q1, u);
00414     Q2_u = BezierII(1, Q2, u);
00415     
00416     /* Compute f(u)/f'(u) */
00417     numerator = (Q_u.x() - P.x()) * (Q1_u.x()) + (Q_u.y() - P.y()) * (Q1_u.y());
00418     denominator = (Q1_u.x()) * (Q1_u.x()) + (Q1_u.y()) * (Q1_u.y()) +
00419                   (Q_u.x() - P.x()) * (Q2_u.x()) + (Q_u.y() - P.y()) * (Q2_u.y());
00420     
00421     /* u = u - f(u)/f'(u) */
00422     uPrime = u - (numerator/denominator);
00423     return (uPrime);
00424 }
00425 
00426 /*
00427  *  Reparameterize:
00428  *  Given set of points and their parameterization, try to find
00429  *   a better parameterization.
00430  *
00431  */
00432 static double *Reparameterize(QPtrList<KoPoint> points,int first,int last,double *u,KoPoint *curve)
00433 {
00434     int     nPts = last-first+1;    
00435     int     i;
00436     double  *uPrime;        /*  New parameter values    */
00437 
00438     uPrime = new double[nPts];
00439     for (i = first; i <= last; i++) {
00440         uPrime[i-first] = NewtonRaphsonRootFind(curve, *points.at(i), u[i-
00441                     first]);
00442     }
00443     return (uPrime);
00444 }
00445 
00446 KoPoint *FitCubic(QPtrList<KoPoint> &points,int first,int last,FitVector tHat1,FitVector tHat2,float error,int &width){
00447     double *u;
00448     double *uPrime;
00449     double maxError;
00450     int splitPoint;
00451     int nPts;
00452     double iterationError;
00453     int maxIterations=4;
00454     FitVector tHatCenter;
00455     KoPoint *curve;
00456     int i;
00457 
00458     width=0;
00459 
00460    
00461     iterationError=error*error;
00462     nPts = last-first+1;
00463 
00464     if(nPts == 2){
00465             double dist = distance(points.at(last), points.at(first)) / 3.0;
00466 
00467         curve = new KoPoint[4];
00468         
00469         curve[0] = *points.at(first);
00470         curve[3] = *points.at(last);
00471 
00472         tHat1.scale(dist);
00473         tHat2.scale(dist);
00474     
00475         curve[1] = tHat1 + curve[0];
00476         curve[2] = tHat2 + curve[3];
00477 
00478         width=4;    
00479         return curve;
00480     }
00481     
00482     /*  Parameterize points, and attempt to fit curve */
00483     u = ChordLengthParameterize(points, first, last);
00484     curve = GenerateBezier(points, first, last, u, tHat1, tHat2);
00485 
00486 
00487     /*  Find max deviation of points to fitted curve */
00488     maxError = ComputeMaxError(points, first, last, curve, u, &splitPoint);
00489     if (maxError < error) {
00490         delete[] u;
00491         width=4;    
00492         return curve;
00493     }
00494 
00495 
00496     /*  If error not too large, try some reparameterization  */
00497     /*  and iteration */
00498     if (maxError < iterationError) {
00499         for (i = 0; i < maxIterations; i++) {
00500             uPrime = Reparameterize(points, first, last, u, curve);
00501             curve = GenerateBezier(points, first, last, uPrime, tHat1, tHat2);
00502             maxError = ComputeMaxError(points, first, last,
00503                     curve, uPrime, &splitPoint);
00504             if (maxError < error) {
00505                 delete[] u;
00506                 width=4;    
00507                 return curve;
00508             }
00509             delete[] u;
00510             u = uPrime;
00511         }
00512     }
00513 
00514     /* Fitting failed -- split at max error point and fit recursively */
00515     delete[] u;
00516     delete[] curve;
00517     tHatCenter = ComputeCenterTangent(points, splitPoint);
00518 
00519     int w1,w2;
00520     KoPoint *cu1=NULL, *cu2=NULL;
00521     cu1 = FitCubic(points, first, splitPoint, tHat1, tHatCenter, error,w1);
00522 
00523     tHatCenter.negate();
00524     cu2 = FitCubic(points, splitPoint, last, tHatCenter, tHat2, error,w2);
00525 
00526     KoPoint *newcurve = new KoPoint[w1+w2];
00527     for(int i=0;i<w1;i++){
00528         newcurve[i]=cu1[i];
00529     }
00530     for(int i=0;i<w2;i++){
00531         newcurve[i+w1]=cu2[i];
00532     }
00533     
00534     delete[] cu1;
00535     delete[] cu2;
00536     width=w1+w2;
00537     return newcurve;
00538 }
00539 
00540 
00541 VPath *bezierFit(QPtrList<KoPoint> &points,float error){
00542     FitVector tHat1, tHat2;
00543 
00544     tHat1 = ComputeLeftTangent(points,0);
00545     tHat2 = ComputeRightTangent(points,points.count()-1);
00546     
00547     int width=0;
00548     KoPoint *curve;
00549     curve = FitCubic(points,0,points.count()-1,tHat1,tHat2,error,width);
00550     
00551     VPath *path = new VPath(NULL);
00552 
00553     if(width>3){
00554         path->moveTo(curve[0]);
00555         path->curveTo(curve[1],curve[2],curve[3]);
00556         for(int i=4;i<width;i+=4){
00557             path->curveTo(curve[i+1],curve[i+2],curve[i+3]);    
00558         }
00559     }
00560 
00561 
00562     delete[] curve;
00563     return path;
00564 }
00565 
KDE Home | KDE Accessibility Home | Description of Access Keys