qwt_spline_cubic.cpp
Go to the documentation of this file.
00001 /* -*- mode: C++ ; c-file-style: "stroustrup" -*- *****************************
00002  * Qwt Widget Library
00003  * Copyright (C) 1997   Josef Wilgen
00004  * Copyright (C) 2002   Uwe Rathmann
00005  *
00006  * This library is free software; you can redistribute it and/or
00007  * modify it under the terms of the Qwt License, Version 1.0
00008  *****************************************************************************/
00009 
00010 #include "qwt_spline_cubic.h"
00011 #include <qdebug.h>
00012 
00013 #define SLOPES_INCREMENTAL 0
00014 #define KAHAN 0
00015 
00016 namespace QwtSplineCubicP
00017 {
00018     class KahanSum 
00019     {
00020     public:
00021         inline KahanSum( double value = 0.0 ):
00022             d_sum( value ),
00023             d_carry( 0.0 )
00024         {
00025         }
00026 
00027         inline void reset()
00028         {
00029             d_sum = d_carry = 0.0;
00030         }
00031 
00032         inline double value() const 
00033         {
00034             return d_sum;
00035         }
00036 
00037         inline void add( double value )
00038         {
00039             const double y = value - d_carry;
00040             const double t = d_sum + y;
00041 
00042             d_carry = ( t - d_sum ) - y;
00043             d_sum = t;
00044         }
00045 
00046         static inline double sum3( double d1, double d2, double d3 )
00047         {
00048             KahanSum sum( d1 );
00049             sum.add( d2 );
00050             sum.add( d3 );
00051 
00052             return sum.value();
00053         }
00054 
00055         static inline double sum4( double d1, double d2, double d3, double d4 )
00056         {
00057             KahanSum sum( d1 );
00058             sum.add( d2 );
00059             sum.add( d3 );
00060             sum.add( d4 );
00061 
00062             return sum.value();
00063         }   
00064 
00065 
00066     private:
00067         double d_sum; 
00068         double d_carry; // The carry from the previous operation
00069     };
00070 
00071     class CurvatureStore
00072     {
00073     public:
00074         inline void setup( int size )
00075         {
00076             d_curvatures.resize( size );
00077             d_cv = d_curvatures.data();
00078         }
00079 
00080         inline void storeFirst( double, 
00081             const QPointF &, const QPointF &, double b1, double )
00082         {
00083             d_cv[0] = 2.0 * b1;
00084         }
00085 
00086         inline void storeNext( int index, double,
00087             const QPointF &, const QPointF &, double, double b2 )
00088         {
00089             d_cv[index] = 2.0 * b2;
00090         }
00091 
00092         inline void storeLast( double,
00093             const QPointF &, const QPointF &, double, double b2 )
00094         {
00095             d_cv[d_curvatures.size() - 1] = 2.0 * b2;
00096         }
00097 
00098         inline void storePrevious( int index, double,
00099             const QPointF &, const QPointF &, double b1, double )
00100         {
00101             d_cv[index] = 2.0 * b1;
00102         }
00103 
00104         inline void closeR()
00105         {
00106             d_cv[0] = d_cv[d_curvatures.size() - 1];
00107         }
00108 
00109         QVector<double> curvatures() const { return d_curvatures; }
00110 
00111     private:
00112         QVector<double> d_curvatures;
00113         double *d_cv;
00114     };
00115 
00116     class SlopeStore
00117     {
00118     public:
00119         inline void setup( int size )
00120         {
00121             d_slopes.resize( size );
00122             d_m = d_slopes.data();
00123         }
00124         
00125         inline void storeFirst( double h, 
00126             const QPointF &p1, const QPointF &p2, double b1, double b2 )
00127         {  
00128             const double s = ( p2.y() - p1.y() ) / h;
00129             d_m[0] = s - h * ( 2.0 * b1 + b2 ) / 3.0;
00130 #if KAHAN
00131             d_sum.add( d_m[0] );
00132 #endif
00133         }
00134         
00135         inline void storeNext( int index, double h,
00136             const QPointF &p1, const QPointF &p2, double b1, double b2 )
00137         {
00138 #if SLOPES_INCREMENTAL
00139             Q_UNUSED( p1 )
00140             Q_UNUSED( p2 )
00141 #if KAHAN
00142             d_sum.add( ( b1 + b2 ) * h );
00143             d_m[index] = d_sum.value();
00144 #else
00145             d_m[index] = d_m[index-1] + ( b1 + b2 ) * h;
00146 #endif
00147 #else
00148             const double s = ( p2.y() - p1.y() ) / h;
00149             d_m[index] = s + h * ( b1 + 2.0 * b2 ) / 3.0;
00150 #endif
00151         }
00152 
00153         inline void storeLast( double h,
00154             const QPointF &p1, const QPointF &p2, double b1, double b2 )
00155         {   
00156             const double s = ( p2.y() - p1.y() ) / h;
00157             d_m[d_slopes.size() - 1] = s + h * ( b1 + 2.0 * b2 ) / 3.0;
00158 #if KAHAN
00159             d_sum.add( d_m[d_slopes.size() - 1] );
00160 #endif
00161         }
00162 
00163         inline void storePrevious( int index, double h,
00164             const QPointF &p1, const QPointF &p2, double b1, double b2 )
00165         {
00166 #if SLOPES_INCREMENTAL
00167             Q_UNUSED( p1 )
00168             Q_UNUSED( p2 )
00169 #if KAHAN
00170             d_sum.add( -( b1 + b2 ) * h );
00171             d_m[index] = d_sum.value();
00172 #else
00173             d_m[index] = d_m[index+1] - ( b1 + b2 ) * h;
00174 #endif
00175 
00176 #else
00177             const double s = ( p2.y() - p1.y() ) / h;
00178             d_m[index] = s - h * ( 2.0 * b1 + b2 ) / 3.0;
00179 #endif
00180         }
00181 
00182         inline void closeR()
00183         {
00184             d_m[0] = d_m[d_slopes.size() - 1];
00185         }
00186 
00187         QVector<double> slopes() const { return d_slopes; }
00188     
00189     private:
00190         QVector<double> d_slopes;
00191         double *d_m;
00192 #if SLOPES_INCREMENTAL
00193         KahanSum d_sum;
00194 #endif
00195     };
00196 };
00197 
00198 namespace QwtSplineCubicP
00199 {
00200     class Equation2
00201     {
00202     public:
00203         inline Equation2()
00204         {
00205         }
00206 
00207         inline Equation2( double p0, double q0, double r0 ):
00208             p( p0 ),
00209             q( q0 ),
00210             r( r0 )
00211         {
00212         }
00213 
00214         inline void setup( double p0 , double q0, double r0 )
00215         {
00216             p = p0;
00217             q = q0;
00218             r = r0;
00219         }
00220             
00221         inline Equation2 normalized() const
00222         {
00223             Equation2 c;
00224             c.p = 1.0;
00225             c.q = q / p;
00226             c.r = r / p;
00227 
00228             return c;
00229         }
00230 
00231         inline double resolved1( double x2 ) const
00232         {
00233             return ( r - q * x2 ) / p;
00234         }
00235 
00236         inline double resolved2( double x1 ) const
00237         {
00238             return ( r - p * x1 ) / q;
00239         }
00240 
00241         inline double resolved1( const Equation2 &eq ) const
00242         {
00243             // find x1
00244             double k = q / eq.q;
00245             return ( r - k * eq.r ) / ( p - k * eq.p );
00246         }
00247 
00248         inline double resolved2( const Equation2 &eq ) const
00249         {
00250             // find x2
00251             const double k = p / eq.p;
00252             return ( r - k * eq.r ) / ( q - k * eq.q );
00253         }
00254 
00255         // p * x1 + q * x2 = r
00256         double p, q, r;
00257     };
00258 
00259     class Equation3
00260     {
00261     public:
00262         inline Equation3()
00263         {
00264         }
00265 
00266         inline Equation3( const QPointF &p1, const QPointF &p2, const QPointF &p3 )
00267         {
00268             const double h1 = p2.x() - p1.x();
00269             const double s1  = ( p2.y() - p1.y() ) / h1;
00270 
00271             const double h2 = p3.x() - p2.x();
00272             const double s2  = ( p3.y() - p2.y() ) / h2;
00273 
00274             p = h1;
00275             q = 2 * ( h1 + h2 );
00276             u = h2;
00277             r = 3 * ( s2 - s1 );
00278         }
00279 
00280         inline Equation3( double cp, double cq, double du, double dr ):
00281             p( cp ),
00282             q( cq ),
00283             u( du ),
00284             r( dr )
00285         {
00286         }
00287     
00288         inline bool operator==( const Equation3 &c ) const
00289         {
00290             return ( p == c.p ) && ( q == c.q ) && 
00291                 ( u == c.u ) && ( r == c.r );
00292         }
00293 
00294         inline void setup( double cp, double cq, double du, double dr )
00295         {
00296             p = cp;
00297             q = cq;
00298             u = du;
00299             r = dr;
00300         }
00301 
00302         inline Equation3 normalized() const
00303         {
00304             Equation3 c;
00305             c.p = 1.0;
00306             c.q = q / p;
00307             c.u = u / p;
00308             c.r = r / p;
00309 
00310             return c;
00311         }
00312 
00313         inline Equation2 substituted1( const Equation3 &eq ) const
00314         {
00315             // eliminate x1
00316             const double k = p / eq.p;
00317             return Equation2( q - k * eq.q, u - k * eq.u, r - k * eq.r );
00318         }
00319 
00320         inline Equation2 substituted2( const Equation3 &eq ) const
00321         {
00322             // eliminate x2
00323 
00324             const double k = q / eq.q;
00325             return Equation2( p - k * eq.p, u - k * eq.u, r - k * eq.r );
00326         }
00327 
00328         inline Equation2 substituted3( const Equation3 &eq ) const
00329         {
00330             // eliminate x3
00331 
00332             const double k = u / eq.u;
00333             return Equation2( p - k * eq.p, q - k * eq.q, r - k * eq.r );
00334         }
00335 
00336         inline Equation2 substituted1( const Equation2 &eq ) const
00337         {
00338             // eliminate x1
00339             const double k = p / eq.p;
00340             return Equation2( q - k * eq.q, u, r - k * eq.r );
00341         }
00342 
00343         inline Equation2 substituted3( const Equation2 &eq ) const
00344         {
00345             // eliminate x3
00346 
00347             const double k = u / eq.q;
00348             return Equation2( p, q - k * eq.p, r - k * eq.r );
00349         }
00350 
00351 
00352         inline double resolved1( double x2, double x3 ) const
00353         {
00354             return ( r - q * x2 - u * x3 ) / p;
00355         }
00356 
00357         inline double resolved2( double x1, double x3 ) const
00358         {
00359             return ( r - u * x3 - p * x1 ) / q;
00360         }
00361 
00362         inline double resolved3( double x1, double x2 ) const
00363         {
00364             return ( r - p * x1 - q * x2 ) / u;
00365         }
00366 
00367         // p * x1 + q * x2 + u * x3 = r
00368         double p, q, u, r;
00369     };
00370 };
00371          
00372 QDebug operator<<( QDebug debug, const QwtSplineCubicP::Equation2 &eq )
00373 {
00374     debug.nospace() << "EQ2(" << eq.p << ", " << eq.q << ", " << eq.r << ")";
00375     return debug.space();
00376 }
00377 
00378 QDebug operator<<( QDebug debug, const QwtSplineCubicP::Equation3 &eq )
00379 {
00380     debug.nospace() << "EQ3(" << eq.p << ", " 
00381         << eq.q << ", " << eq.u << ", " << eq.r << ")";
00382     return debug.space();
00383 }
00384 
00385 namespace QwtSplineCubicP
00386 {
00387     template <class T>
00388     class EquationSystem
00389     {
00390     public:
00391         void setStartCondition( double p, double q, double u, double r )
00392         {
00393             d_conditionsEQ[0].setup( p, q, u, r );
00394         }
00395 
00396         void setEndCondition( double p, double q, double u, double r )
00397         {
00398             d_conditionsEQ[1].setup( p, q, u, r );
00399         }
00400 
00401         const T &store() const 
00402         { 
00403             return d_store;
00404         }
00405 
00406         void resolve( const QPolygonF &p ) 
00407         {
00408             const int n = p.size();
00409             if ( n < 3 )
00410                 return;
00411 
00412             if ( d_conditionsEQ[0].p == 0.0 ||
00413                 ( d_conditionsEQ[0].q == 0.0 && d_conditionsEQ[0].u != 0.0 ) )
00414             {
00415                 return;
00416             }
00417 
00418             if ( d_conditionsEQ[1].u == 0.0 ||
00419                 ( d_conditionsEQ[1].q == 0.0 && d_conditionsEQ[1].p != 0.0 ) )
00420             {
00421                 return;
00422             }
00423 
00424             const double h0 = p[1].x() - p[0].x();
00425             const double h1 = p[2].x() - p[1].x();
00426             const double hn = p[n-1].x() - p[n-2].x();
00427 
00428             d_store.setup( n );
00429             
00430 
00431             if ( n == 3 )
00432             {
00433                 // For certain conditions the first/last point does not 
00434                 // necessarily meet the spline equation and we would
00435                 // have many solutions. In this case we resolve using
00436                 // the spline equation - as for all other conditions.
00437 
00438                 const Equation3 eqSpline0( p[0], p[1], p[2] ); // ???
00439                 const Equation2 eq0 = d_conditionsEQ[0].substituted1( eqSpline0 );
00440 
00441                 // The equation system can be solved without substitution
00442                 // from the start/end conditions and eqSpline0 ( = eqSplineN ).
00443 
00444                 double b1;
00445 
00446                 if ( d_conditionsEQ[0].normalized() == d_conditionsEQ[1].normalized() )
00447                 {
00448                     // When we have 3 points only and start/end conditions
00449                     // for 3 points mean the same condition the system
00450                     // is under-determined and has many solutions.
00451                     // We chose b1 = 0.0
00452 
00453                     b1 = 0.0;
00454                 }
00455                 else 
00456                 {
00457                     const Equation2 eq = d_conditionsEQ[1].substituted1( eqSpline0 );
00458                     b1 = eq0.resolved1( eq );
00459                 }
00460 
00461                 const double b2 = eq0.resolved2( b1 );
00462                 const double b0 = eqSpline0.resolved1( b1, b2 );
00463 
00464                 d_store.storeFirst( h0, p[0], p[1], b0, b1 );
00465                 d_store.storeNext( 1, h0, p[0], p[1], b0, b1 );
00466                 d_store.storeNext( 2, h1, p[1], p[2], b1, b2 );
00467 
00468                 return; 
00469             }
00470 
00471             const Equation3 eqSplineN( p[n-3], p[n-2], p[n-1] );
00472             const Equation2 eqN = d_conditionsEQ[1].substituted3( eqSplineN );
00473 
00474             Equation2 eq = eqN;
00475             if ( n > 4 )
00476             {
00477                 const Equation3 eqSplineR( p[n-4], p[n-3], p[n-2] );
00478                 eq = eqSplineR.substituted3( eq );
00479                 eq = substituteSpline( p, eq );
00480             }
00481 
00482             const Equation3 eqSpline0( p[0], p[1], p[2] ); 
00483 
00484             double b0, b1;
00485             if ( d_conditionsEQ[0].u == 0.0 )
00486             {
00487                 eq = eqSpline0.substituted3( eq );
00488 
00489                 const Equation3 &eq0 = d_conditionsEQ[0];
00490                 b0 = Equation2( eq0.p, eq0.q, eq0.r ).resolved1( eq );
00491                 b1 = eq.resolved2( b0 );
00492             }
00493             else
00494             {
00495                 const Equation2 eqX = d_conditionsEQ[0].substituted3( eq );
00496                 const Equation2 eqY = eqSpline0.substituted3( eq );
00497 
00498                 b0 = eqY.resolved1( eqX );
00499                 b1 = eqY.resolved2( b0 );
00500             }
00501 
00502             d_store.storeFirst( h0, p[0], p[1], b0, b1 );
00503             d_store.storeNext( 1, h0, p[0], p[1], b0, b1 );
00504 
00505             const double bn2 = resolveSpline( p, b1 );
00506 
00507             const double bn1 = eqN.resolved2( bn2 );
00508             const double bn0 = d_conditionsEQ[1].resolved3( bn2, bn1 );
00509 
00510             const double hx = p[n-2].x() - p[n-3].x();
00511             d_store.storeNext( n - 2, hx, p[n-3], p[n-2], bn2, bn1 );
00512             d_store.storeNext( n - 1, hn, p[n-2], p[n-1], bn1, bn0 );
00513         }
00514 
00515     private:
00516         Equation2 substituteSpline( const QPolygonF &points, const Equation2 &eq )
00517         {
00518             const int n = points.size();
00519 
00520             d_eq.resize( n - 2 );
00521             d_eq[n-3] = eq;
00522 
00523             // eq[i].resolved2( b[i-1] ) => b[i]
00524 
00525             double slope2 = ( points[n-3].y() - points[n-4].y() ) / eq.p;
00526 
00527             for ( int i = n - 4; i > 1; i-- )
00528             {
00529                 const Equation2 &eq2 = d_eq[i+1];
00530                 Equation2 &eq1 = d_eq[i];
00531 
00532                 eq1.p = points[i].x() - points[i-1].x();
00533                 const double slope1 = ( points[i].y() - points[i-1].y() ) / eq1.p;
00534 
00535                 const double v = eq2.p / eq2.q;
00536 
00537                 eq1.q = 2.0 * ( eq1.p + eq2.p ) - v * eq2.p;
00538                 eq1.r = 3.0 * ( slope2 - slope1 ) - v * eq2.r;
00539 
00540                 slope2 = slope1;
00541             }
00542 
00543             return d_eq[2];
00544         }
00545 
00546         double resolveSpline( const QPolygonF &points, double b1 )
00547         {
00548             const int n = points.size();
00549             const QPointF *p = points.constData();
00550 
00551             for ( int i = 2; i < n - 2; i++ )
00552             {
00553                 // eq[i].resolved2( b[i-1] ) => b[i]
00554                 const double b2 = d_eq[i].resolved2( b1 );
00555                 d_store.storeNext( i, d_eq[i].p, p[i-1], p[i], b1, b2 );
00556 
00557                 b1 = b2;
00558             }
00559 
00560             return b1;
00561         }
00562 
00563     private:
00564         Equation3 d_conditionsEQ[2];
00565         QVector<Equation2> d_eq;
00566         T d_store;
00567     };
00568 
00569     template <class T>
00570     class EquationSystem2
00571     {
00572     public:
00573         const T &store() const
00574         {
00575             return d_store;
00576         }
00577 
00578         void resolve( const QPolygonF &p )
00579         {
00580             const int n = p.size();
00581 
00582             if ( p[n-1].y() != p[0].y() )
00583             {
00584                 // TODO ???
00585             }
00586 
00587             const double h0 = p[1].x() - p[0].x();
00588             const double s0 = ( p[1].y() - p[0].y() ) / h0;
00589 
00590             if ( n == 3 )
00591             {
00592                 const double h1 = p[2].x() - p[1].x();
00593                 const double s1 = ( p[2].y() - p[1].y() ) / h1;
00594 
00595                 const double b = 3.0 * ( s0 - s1 ) / ( h0 + h1 );
00596 
00597                 d_store.setup( 3 );
00598                 d_store.storeLast( h1, p[1], p[2], -b, b );
00599                 d_store.storePrevious( 1, h1, p[1], p[2], -b, b );
00600                 d_store.closeR();
00601 
00602                 return;
00603             }
00604 
00605             const double hn = p[n-1].x() - p[n-2].x();
00606 
00607             Equation2 eqn, eqX;
00608             substitute( p, eqn, eqX );
00609 
00610             const double b0 = eqn.resolved2( eqX );
00611             const double bn = eqn.resolved1( b0 );
00612 
00613             d_store.setup( n );
00614             d_store.storeLast( hn, p[n-2], p[n-1], bn, b0 );
00615             d_store.storePrevious( n - 2, hn, p[n-2], p[n-1], bn, b0 );
00616 
00617             resolveSpline( p, b0, bn );
00618 
00619             d_store.closeR();
00620         }
00621 
00622     private:
00623 
00624         void substitute( const QPolygonF &points, Equation2 &eqn, Equation2 &eqX )
00625         {
00626             const int n = points.size();
00627 
00628             const double hn = points[n-1].x() - points[n-2].x();
00629 
00630             const Equation3 eqSpline0( points[0], points[1], points[2] );
00631             const Equation3 eqSplineN( 
00632                 QPointF( points[0].x() - hn, points[n-2].y() ), points[0], points[1] );
00633 
00634             d_eq.resize( n - 1 );
00635             
00636             double dq = 0;
00637             double dr = 0;
00638 
00639             d_eq[1] = eqSpline0;
00640 
00641             double slope1 = ( points[2].y() - points[1].y() ) / d_eq[1].u;
00642 
00643             // a) p1 * b[0] + q1 * b[1] + u1 * b[2] = r1
00644             // b) p2 * b[n-2] + q2 * b[0] + u2 * b[1] = r2
00645             // c) pi * b[i-1] + qi * b[i] + ui * b[i+1] = ri
00646             //
00647             // Using c) we can substitute b[i] ( starting from 2 ) by b[i+1] 
00648             // until we reach n-1. As we know, that b[0] == b[n-1] we found 
00649             // an equation where only 2 coefficients ( for b[n-2], b[0] ) are left unknown.
00650             // Each step we have an equation that depends on b[0], b[i] and b[i+1] 
00651             // that can also be used to substitute b[i] in b). Ding so we end up with another
00652             // equation depending on b[n-2], b[0] only.
00653             // Finally 2 equations with 2 coefficients can be solved.
00654 
00655             for ( int i = 2; i < n - 1; i++ )
00656             {
00657                 const Equation3 &eq1 = d_eq[i-1];
00658                 Equation3 &eq2 = d_eq[i];
00659 
00660                 dq += eq1.p * eq1.p / eq1.q;
00661                 dr += eq1.p * eq1.r / eq1.q;
00662 
00663                 eq2.u = points[i+1].x() - points[i].x();
00664                 const double slope2 = ( points[i+1].y() - points[i].y() ) / eq2.u;
00665 
00666                 const double k = eq1.u / eq1.q;
00667 
00668                 eq2.p = -eq1.p * k;
00669                 eq2.q = 2.0 * ( eq1.u + eq2.u ) - eq1.u * k;
00670                 eq2.r = 3.0 * ( slope2 - slope1 ) - eq1.r * k;
00671 
00672                 slope1 = slope2;
00673             }
00674 
00675 
00676             // b[0] * d_p[n-2] + b[n-2] * d_q[n-2] + b[n-1] * pN = d_r[n-2]
00677             eqn.setup( d_eq[n-2].q, d_eq[n-2].p + eqSplineN.p , d_eq[n-2].r );
00678 
00679             // b[n-2] * pN + b[0] * ( qN - dq ) + b[n-2] * d_p[n-2] = rN - dr
00680             eqX.setup( d_eq[n-2].p + eqSplineN.p, eqSplineN.q - dq, eqSplineN.r - dr );
00681         }
00682 
00683         void resolveSpline( const QPolygonF &points, double b0, double bi )
00684         {
00685             const int n = points.size();
00686 
00687             for ( int i = n - 3; i >= 1; i-- )
00688             {
00689                 const Equation3 &eq = d_eq[i];
00690 
00691                 const double b = eq.resolved2( b0, bi );
00692                 d_store.storePrevious( i, eq.u, points[i], points[i+1], b, bi );
00693 
00694                 bi = b;
00695             }
00696         }
00697 
00698         void resolveSpline2( const QPolygonF &points,
00699             double b0, double bi, QVector<double> &m )
00700         {
00701             const int n = points.size();
00702 
00703             bi = d_eq[0].resolved3( b0, bi );
00704 
00705             for ( int i = 1; i < n - 2; i++ )
00706             {
00707                 const Equation3 &eq = d_eq[i];
00708 
00709                 const double b = eq.resolved3( b0, bi );
00710                 m[i+1] = m[i] + ( b + bi ) * d_eq[i].u;
00711 
00712                 bi = b;
00713             }
00714         }
00715 
00716         void resolveSpline3( const QPolygonF &points,
00717             double b0, double b1, QVector<double> &m )
00718         {
00719             const int n = points.size();
00720 
00721             double h0 = ( points[1].x() - points[0].x() );
00722             double s0 = ( points[1].y() - points[0].y() ) / h0;
00723 
00724             m[1] = m[0] + ( b0 + b1 ) * h0;
00725 
00726             for ( int i = 1; i < n - 1; i++ )
00727             {
00728                 const double h1 = ( points[i+1].x() - points[i].x() );
00729                 const double s1 = ( points[i+1].y() - points[i].y() ) / h1;
00730 
00731                 const double r = 3.0 * ( s1 - s0 );
00732 
00733                 const double b2 = ( r - h0 * b0 - 2.0 * ( h0 + h1 ) * b1 ) / h1;
00734                 m[i+1] = m[i] + ( b1 + b2 ) * h1;
00735 
00736                 h0 = h1;
00737                 s0 = s1;
00738                 b0 = b1;
00739                 b1 = b2;
00740             }
00741         }
00742 
00743         void resolveSpline4( const QPolygonF &points,
00744             double b2, double b1, QVector<double> &m )
00745         {
00746             const int n = points.size();
00747 
00748             double h2 = ( points[n-1].x() - points[n-2].x() );
00749             double s2 = ( points[n-1].y() - points[n-2].y() ) / h2;
00750 
00751             for ( int i = n - 2; i > 1; i-- )
00752             {
00753                 const double h1 = ( points[i].x() - points[i-1].x() );
00754                 const double s1 = ( points[i].y() - points[i-1].y() ) / h1;
00755 
00756                 const double r = 3.0 * ( s2 - s1 );
00757                 const double k = 2.0 * ( h1 + h2 );
00758 
00759                 const double b0 = ( r - h2 * b2 - k * b1 ) / h1;
00760 
00761                 m[i-1] = m[i] - ( b0 + b1 ) * h1;
00762 
00763                 h2 = h1;
00764                 s2 = s1;
00765                 b2 = b1;
00766                 b1 = b0;
00767             }
00768         }
00769 
00770     public:
00771         QVector<Equation3> d_eq;
00772         T d_store;
00773     };
00774 }
00775 
00776 static void qwtSetupEndEquations( 
00777     int conditionBegin, double valueBegin, int conditionEnd, double valueEnd, 
00778     const QPolygonF &points, QwtSplineCubicP::Equation3 eq[2] )
00779 {
00780     const int n = points.size();
00781 
00782     const double h0 = points[1].x() - points[0].x();
00783     const double s0 = ( points[1].y() - points[0].y() ) / h0;
00784 
00785     const double hn = ( points[n-1].x() - points[n-2].x() );
00786     const double sn = ( points[n-1].y() - points[n-2].y() ) / hn;
00787 
00788     switch( conditionBegin )
00789     {
00790         case QwtSpline::Clamped1:
00791         {
00792             // first derivative at end points given
00793 
00794             // 3 * a1 * h + b1 = b2
00795             // a1 * h * h + b1 * h + c1 = s
00796 
00797             // c1 = slopeBegin
00798             // => b1 * ( 2 * h / 3.0 ) + b2 * ( h / 3.0 ) = s - slopeBegin
00799 
00800             // c2 = slopeEnd
00801             // => b1 * ( 1.0 / 3.0 ) + b2 * ( 2.0 / 3.0 ) = ( slopeEnd - s ) / h;
00802 
00803             eq[0].setup( 2 * h0 / 3.0, h0 / 3.0, 0.0, s0 - valueBegin );
00804             break;
00805         }
00806         case QwtSpline::Clamped2:
00807         {
00808             // second derivative at end points given
00809 
00810             // b0 = 0.5 * cvStart
00811             // => b0 * 1.0 + b1 * 0.0 = 0.5 * cvStart
00812 
00813             // b1 = 0.5 * cvEnd
00814             // => b0 * 0.0 + b1 * 1.0 = 0.5 * cvEnd
00815 
00816             eq[0].setup( 1.0, 0.0, 0.0, 0.5 * valueBegin ); 
00817             break;
00818         }
00819         case QwtSpline::Clamped3:
00820         {
00821             // third derivative at end point given
00822 
00823             // 3 * a * h0 + b[0] = b[1]
00824 
00825             // a = marg_0 / 6.0
00826             // => b[0] * 1.0 + b[1] * ( -1.0 ) = -0.5 * v0 * h0
00827 
00828             // a = marg_n / 6.0
00829             // => b[n-2] * 1.0 + b[n-1] * ( -1.0 ) = -0.5 * v1 * h5
00830 
00831             eq[0].setup( 1.0, -1.0, 0.0, -0.5 * valueBegin * h0 ); 
00832 
00833             break;
00834         }
00835         case QwtSpline::LinearRunout:
00836         {
00837             const double r0 = qBound( 0.0, valueBegin, 1.0 );
00838             if ( r0 == 0.0 )
00839             {
00840                 // clamping s0
00841                 eq[0].setup( 2 * h0 / 3.0, h0 / 3.0, 0.0, 0.0 );
00842             }
00843             else
00844             {
00845                 eq[0].setup( 1.0 + 2.0 / r0, 2.0 + 1.0 / r0, 0.0, 0.0 );
00846             }
00847             break;
00848         }
00849         case QwtSplineC2::NotAKnot:
00850         case QwtSplineC2::CubicRunout:
00851         {
00852             // building one cubic curve from 3 points
00853 
00854             double v0;
00855 
00856             if ( conditionBegin == QwtSplineC2::CubicRunout )
00857             {
00858                 // first/last point are the endpoints of the curve
00859 
00860                 // b0 = 2 * b1 - b2
00861                 // => 1.0 * b0 - 2 * b1 + 1.0 * b2 = 0.0
00862 
00863                 v0 = 1.0;
00864             }
00865             else
00866             {
00867                 // first/last points are on the curve, 
00868                 // the imaginary endpoints have the same distance as h0/hn
00869 
00870                 v0 = h0 / ( points[2].x() - points[1].x() );
00871             }
00872 
00873             eq[0].setup( 1.0, -( 1.0 + v0 ), v0, 0.0 ); 
00874             break;
00875         }
00876         default:
00877         {
00878             // a natural spline, where the
00879             // second derivative at end points set to 0.0
00880             eq[0].setup( 1.0, 0.0, 0.0, 0.0 ); 
00881             break;
00882         }
00883     }
00884 
00885     switch( conditionEnd )
00886     {
00887         case QwtSpline::Clamped1:
00888         {
00889             // first derivative at end points given
00890             eq[1].setup( 0.0, 1.0 / 3.0 * hn, 2.0 / 3.0 * hn, valueEnd - sn );
00891             break;
00892         }
00893         case QwtSpline::Clamped2:
00894         {
00895             // second derivative at end points given
00896             eq[1].setup( 0.0, 0.0, 1.0, 0.5 * valueEnd ); 
00897             break;
00898         }
00899         case QwtSpline::Clamped3:
00900         {
00901             // third derivative at end point given
00902             eq[1].setup( 0.0, 1.0, -1.0, -0.5 * valueEnd * hn ); 
00903             break;
00904         }
00905         case QwtSpline::LinearRunout:
00906         {
00907             const double rn = qBound( 0.0, valueEnd, 1.0 );
00908             if ( rn == 0.0 )
00909             {
00910                 // clamping sn
00911                 eq[1].setup( 0.0, 1.0 / 3.0 * hn, 2.0 / 3.0 * hn, 0.0 );
00912             }
00913             else
00914             {
00915                 eq[1].setup( 0.0, 2.0 + 1.0 / rn, 1.0 + 2.0 / rn, 0.0 );
00916             }
00917 
00918             break;
00919         }
00920         case QwtSplineC2::NotAKnot:
00921         case QwtSplineC2::CubicRunout:
00922         {
00923             // building one cubic curve from 3 points
00924 
00925             double vn;
00926 
00927             if ( conditionEnd == QwtSplineC2::CubicRunout )
00928             {
00929                 // last point is the endpoints of the curve
00930                 vn = 1.0;
00931             }
00932             else
00933             {
00934                 // last points on the curve, 
00935                 // the imaginary endpoints have the same distance as hn
00936 
00937                 vn = hn / ( points[n-2].x() - points[n-3].x() );
00938             }
00939 
00940             eq[1].setup( vn, -( 1.0 + vn ), 1.0, 0.0 ); 
00941             break;
00942         }
00943         default:
00944         {
00945             // a natural spline, where the
00946             // second derivative at end points set to 0.0
00947             eq[1].setup( 0.0, 0.0, 1.0, 0.0 ); 
00948             break;
00949         }
00950     }
00951 }
00952 
00953 class QwtSplineCubic::PrivateData
00954 {
00955 public:
00956     PrivateData()
00957     {
00958     }
00959 };
00960 
00965 QwtSplineCubic::QwtSplineCubic()
00966 {
00967     d_data = new PrivateData;
00968 
00969     // a natural spline
00970 
00971     setBoundaryCondition( QwtSpline::AtBeginning, QwtSpline::Clamped2 );
00972     setBoundaryValue( QwtSpline::AtBeginning, 0.0 );
00973 
00974     setBoundaryCondition( QwtSpline::AtEnd, QwtSpline::Clamped2 );
00975     setBoundaryValue( QwtSpline::AtEnd, 0.0 );
00976 }
00977 
00979 QwtSplineCubic::~QwtSplineCubic()
00980 {
00981     delete d_data;
00982 }
00983 
00990 uint QwtSplineCubic::locality() const
00991 {
00992     return 0;
00993 }
00994 
01007 QVector<double> QwtSplineCubic::slopes( const QPolygonF &points ) const
01008 {
01009     using namespace QwtSplineCubicP;
01010 
01011     if ( points.size() <= 2 )
01012         return QVector<double>();
01013 
01014     if ( ( boundaryType() == QwtSpline::PeriodicPolygon )
01015         || ( boundaryType() == QwtSpline::ClosedPolygon ) )
01016     {
01017         EquationSystem2<SlopeStore> eqs;
01018         eqs.resolve( points );
01019 
01020         return eqs.store().slopes();
01021     }
01022 
01023     if ( points.size() == 3 )
01024     {
01025         if ( boundaryCondition( QwtSpline::AtBeginning ) == QwtSplineCubic::NotAKnot
01026             || boundaryCondition( QwtSpline::AtEnd ) == QwtSplineCubic::NotAKnot )
01027         {
01028 #if 0
01029             const double h0 = points[1].x() - points[0].x();
01030             const double h1 = points[2].x() - points[1].x();
01031 
01032             const double s0 = ( points[1].y() - points[0].y() ) / h0;
01033             const double s1 = ( points[2].y() - points[1].y() ) / h1;
01034 
01035             /*
01036               the system is under-determined and we only 
01037               compute a quadratic spline.                   
01038              */
01039 
01040             const double b = ( s1 - s0 ) / ( h0 + h1 );
01041 
01042             QVector<double> m( 3 );
01043             m[0] = s0 - h0 * b;
01044             m[1] = s1 - h1 * b;
01045             m[2] = s1 + h1 * b;
01046 
01047             return m;
01048 #else
01049             return QVector<double>();
01050 #endif
01051         }
01052     }
01053 
01054     Equation3 eq[2];
01055     qwtSetupEndEquations( 
01056         boundaryCondition( QwtSpline::AtBeginning ), 
01057         boundaryValue( QwtSpline::AtBeginning ), 
01058         boundaryCondition( QwtSpline::AtEnd ), 
01059         boundaryValue( QwtSpline::AtEnd ), 
01060         points, eq );
01061 
01062     EquationSystem<SlopeStore> eqs;
01063     eqs.setStartCondition( eq[0].p, eq[0].q, eq[0].u, eq[0].r );
01064     eqs.setEndCondition( eq[1].p, eq[1].q, eq[1].u, eq[1].r );
01065     eqs.resolve( points );
01066 
01067     return eqs.store().slopes();
01068 }
01069 
01079 QVector<double> QwtSplineCubic::curvatures( const QPolygonF &points ) const
01080 {
01081     using namespace QwtSplineCubicP;
01082 
01083     if ( points.size() <= 2 )
01084         return QVector<double>();
01085 
01086     if ( ( boundaryType() == QwtSpline::PeriodicPolygon )
01087         || ( boundaryType() == QwtSpline::ClosedPolygon ) )
01088     {
01089         EquationSystem2<CurvatureStore> eqs;
01090         eqs.resolve( points );
01091 
01092         return eqs.store().curvatures();
01093     }
01094 
01095     if ( points.size() == 3 )
01096     {
01097         if ( boundaryCondition( QwtSpline::AtBeginning ) == QwtSplineC2::NotAKnot 
01098             || boundaryCondition( QwtSpline::AtEnd ) == QwtSplineC2::NotAKnot )
01099         {
01100             return QVector<double>();
01101         }
01102     }
01103 
01104     Equation3 eq[2];
01105     qwtSetupEndEquations(
01106         boundaryCondition( QwtSpline::AtBeginning ),
01107         boundaryValue( QwtSpline::AtBeginning ),
01108         boundaryCondition( QwtSpline::AtEnd ),
01109         boundaryValue( QwtSpline::AtEnd ),
01110         points, eq );
01111 
01112     EquationSystem<CurvatureStore> eqs;
01113     eqs.setStartCondition( eq[0].p, eq[0].q, eq[0].u, eq[0].r );
01114     eqs.setEndCondition( eq[1].p, eq[1].q, eq[1].u, eq[1].r );
01115     eqs.resolve( points );
01116 
01117     return eqs.store().curvatures();
01118 }
01119 
01131 QPainterPath QwtSplineCubic::painterPath( const QPolygonF &points ) const
01132 {   
01133     // as QwtSplineCubic can calcuate slopes directly we can
01134     // use the implementation of QwtSplineC1 without any performance loss.
01135 
01136     return QwtSplineC1::painterPath( points );
01137 }
01138     
01150 QVector<QLineF> QwtSplineCubic::bezierControlLines( const QPolygonF &points ) const
01151 {   
01152     // as QwtSplineCubic can calcuate slopes directly we can
01153     // use the implementation of QwtSplineC1 without any performance loss.
01154 
01155     return QwtSplineC1::bezierControlLines( points ); 
01156 }   
01157 
01168 QVector<QwtSplinePolynomial> QwtSplineCubic::polynomials( const QPolygonF &points ) const
01169 {
01170     return QwtSplineC2::polynomials( points );
01171 }
01172 


plotjuggler
Author(s): Davide Faconti
autogenerated on Fri Sep 1 2017 02:41:56