ground/gcs/src/libs/qwt/src/qwt_spline.cpp

   1 /* -*- mode: C++ ; c-file-style: "stroustrup" -*- *****************************
   2  * Qwt Widget Library
   3  * Copyright (C) 1997   Josef Wilgen
   4  * Copyright (C) 2002   Uwe Rathmann
   5  *
   6  * This library is free software; you can redistribute it and/or
   7  * modify it under the terms of the Qwt License, Version 1.0
   8  *****************************************************************************/
   9
  10 #include "qwt_spline.h"
  11 #include "qwt_math.h"
  12
  13 class QwtSpline::PrivateData
  14 {
  15 public:
  16     PrivateData():
  17         splineType( QwtSpline::Natural )
  18     {
  19     }
  20
  21     QwtSpline::SplineType splineType;
  22
  23     // coefficient vectors
  24     QVector<double> a;
  25     QVector<double> b;
  26     QVector<double> c;
  27
  28     // control points
  29     QPolygonF points;
  30 };
  31
  32 static int lookup( double x, const QPolygonF &values )
  33 {
  34 #if 0
  35 //qLowerBound/qHigherBound ???
  36 #endif
  37     int i1;
  38     const int size = values.size();
  39
  40     if ( x <= values[0].x() )
  41         i1 = 0;
  42     else if ( x >= values[size - 2].x() )
  43         i1 = size - 2;
  44     else
  45     {
  46         i1 = 0;
  47         int i2 = size - 2;
  48         int i3 = 0;
  49
  50         while ( i2 - i1 > 1 )
  51         {
  52             i3 = i1 + ( ( i2 - i1 ) >> 1 );
  53
  54             if ( values[i3].x() > x )
  55                 i2 = i3;
  56             else
  57                 i1 = i3;
  58         }
  59     }
  60     return i1;
  61 }
  62
  63 //! Constructor
  64 QwtSpline::QwtSpline()
  65 {
  66     d_data = new PrivateData;
  67 }
  68
  69 /*!
  70    Copy constructor
  71    \param other Spline used for initialization
  72 */
  73 QwtSpline::QwtSpline( const QwtSpline& other )
  74 {
  75     d_data = new PrivateData( *other.d_data );
  76 }
  77
  78 /*!
  79    Assignment operator
  80    \param other Spline used for initialization
  81    \return *this
  82 */
  83 QwtSpline &QwtSpline::operator=( const QwtSpline & other )
  84 {
  85     *d_data = *other.d_data;
  86     return *this;
  87 }
  88
  89 //! Destructor
  90 QwtSpline::~QwtSpline()
  91 {
  92     delete d_data;
  93 }
  94
  95 /*!
  96    Select the algorithm used for calculating the spline
  97
  98    \param splineType Spline type
  99    \sa splineType()
 100 */
 101 void QwtSpline::setSplineType( SplineType splineType )
 102 {
 103     d_data->splineType = splineType;
 104 }
 105
 106 /*!
 107    \return the spline type
 108    \sa setSplineType()
 109 */
 110 QwtSpline::SplineType QwtSpline::splineType() const
 111 {
 112     return d_data->splineType;
 113 }
 114
 115 /*!
 116   \brief Calculate the spline coefficients
 117
 118   Depending on the value of \a periodic, this function
 119   will determine the coefficients for a natural or a periodic
 120   spline and store them internally.
 121
 122   \param points Points
 123   \return true if successful
 124   \warning The sequence of x (but not y) values has to be strictly monotone
 125            increasing, which means <code>points[i].x() < points[i+1].x()</code>.
 126        If this is not the case, the function will return false
 127 */
 128 bool QwtSpline::setPoints( const QPolygonF& points )
 129 {
 130     const int size = points.size();
 131     if ( size <= 2 )
 132     {
 133         reset();
 134         return false;
 135     }
 136
 137     d_data->points = points;
 138
 139     d_data->a.resize( size - 1 );
 140     d_data->b.resize( size - 1 );
 141     d_data->c.resize( size - 1 );
 142
 143     bool ok;
 144     if ( d_data->splineType == Periodic )
 145         ok = buildPeriodicSpline( points );
 146     else
 147         ok = buildNaturalSpline( points );
 148
 149     if ( !ok )
 150         reset();
 151
 152     return ok;
 153 }
 154
 155 /*!
 156    \return Points, that have been by setPoints()
 157 */
 158 QPolygonF QwtSpline::points() const
 159 {
 160     return d_data->points;
 161 }
 162
 163 //! \return A coefficients
 164 const QVector<double> &QwtSpline::coefficientsA() const
 165 {
 166     return d_data->a;
 167 }
 168
 169 //! \return B coefficients
 170 const QVector<double> &QwtSpline::coefficientsB() const
 171 {
 172     return d_data->b;
 173 }
 174
 175 //! \return C coefficients
 176 const QVector<double> &QwtSpline::coefficientsC() const
 177 {
 178     return d_data->c;
 179 }
 180
 181
 182 //! Free allocated memory and set size to 0
 183 void QwtSpline::reset()
 184 {
 185     d_data->a.resize( 0 );
 186     d_data->b.resize( 0 );
 187     d_data->c.resize( 0 );
 188     d_data->points.resize( 0 );
 189 }
 190
 191 //! True if valid
 192 bool QwtSpline::isValid() const
 193 {
 194     return d_data->a.size() > 0;
 195 }
 196
 197 /*!
 198   Calculate the interpolated function value corresponding
 199   to a given argument x.
 200
 201   \param x Coordinate
 202   \return Interpolated coordinate
 203 */
 204 double QwtSpline::value( double x ) const
 205 {
 206     if ( d_data->a.size() == 0 )
 207         return 0.0;
 208
 209     const int i = lookup( x, d_data->points );
 210
 211     const double delta = x - d_data->points[i].x();
 212     return( ( ( ( d_data->a[i] * delta ) + d_data->b[i] )
 213         * delta + d_data->c[i] ) * delta + d_data->points[i].y() );
 214 }
 215
 216 /*!
 217   \brief Determines the coefficients for a natural spline
 218   \return true if successful
 219 */
 220 bool QwtSpline::buildNaturalSpline( const QPolygonF &points )
 221 {
 222     int i;
 223
 224     const QPointF *p = points.data();
 225     const int size = points.size();
 226
 227     double *a = d_data->a.data();
 228     double *b = d_data->b.data();
 229     double *c = d_data->c.data();
 230
 231     //  set up tridiagonal equation system; use coefficient
 232     //  vectors as temporary buffers
 233     QVector<double> h( size - 1 );
 234     for ( i = 0; i < size - 1; i++ )
 235     {
 236         h[i] = p[i+1].x() - p[i].x();
 237         if ( h[i] <= 0 )
 238             return false;
 239     }
 240
 241     QVector<double> d( size - 1 );
 242     double dy1 = ( p[1].y() - p[0].y() ) / h[0];
 243     for ( i = 1; i < size - 1; i++ )
 244     {
 245         b[i] = c[i] = h[i];
 246         a[i] = 2.0 * ( h[i-1] + h[i] );
 247
 248         const double dy2 = ( p[i+1].y() - p[i].y() ) / h[i];
 249         d[i] = 6.0 * ( dy1 - dy2 );
 250         dy1 = dy2;
 251     }
 252
 253     //
 254     // solve it
 255     //
 256
 257     // L-U Factorization
 258     for ( i = 1; i < size - 2; i++ )
 259     {
 260         c[i] /= a[i];
 261         a[i+1] -= b[i] * c[i];
 262     }
 263
 264     // forward elimination
 265     QVector<double> s( size );
 266     s[1] = d[1];
 267     for ( i = 2; i < size - 1; i++ )
 268         s[i] = d[i] - c[i-1] * s[i-1];
 269
 270     // backward elimination
 271     s[size - 2] = - s[size - 2] / a[size - 2];
 272     for ( i = size - 3; i > 0; i-- )
 273         s[i] = - ( s[i] + b[i] * s[i+1] ) / a[i];
 274     s[size - 1] = s[0] = 0.0;
 275
 276     //
 277     // Finally, determine the spline coefficients
 278     //
 279     for ( i = 0; i < size - 1; i++ )
 280     {
 281         a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i] );
 282         b[i] = 0.5 * s[i];
 283         c[i] = ( p[i+1].y() - p[i].y() ) / h[i]
 284             - ( s[i+1] + 2.0 * s[i] ) * h[i] / 6.0;
 285     }
 286
 287     return true;
 288 }
 289
 290 /*!
 291   \brief Determines the coefficients for a periodic spline
 292   \return true if successful
 293 */
 294 bool QwtSpline::buildPeriodicSpline( const QPolygonF &points )
 295 {
 296     int i;
 297
 298     const QPointF *p = points.data();
 299     const int size = points.size();
 300
 301     double *a = d_data->a.data();
 302     double *b = d_data->b.data();
 303     double *c = d_data->c.data();
 304
 305     QVector<double> d( size - 1 );
 306     QVector<double> h( size - 1 );
 307     QVector<double> s( size );
 308
 309     //
 310     //  setup equation system; use coefficient
 311     //  vectors as temporary buffers
 312     //
 313     for ( i = 0; i < size - 1; i++ )
 314     {
 315         h[i] = p[i+1].x() - p[i].x();
 316         if ( h[i] <= 0.0 )
 317             return false;
 318     }
 319
 320     const int imax = size - 2;
 321     double htmp = h[imax];
 322     double dy1 = ( p[0].y() - p[imax].y() ) / htmp;
 323     for ( i = 0; i <= imax; i++ )
 324     {
 325         b[i] = c[i] = h[i];
 326         a[i] = 2.0 * ( htmp + h[i] );
 327         const double dy2 = ( p[i+1].y() - p[i].y() ) / h[i];
 328         d[i] = 6.0 * ( dy1 - dy2 );
 329         dy1 = dy2;
 330         htmp = h[i];
 331     }
 332
 333     //
 334     // solve it
 335     //
 336
 337     // L-U Factorization
 338     a[0] = qSqrt( a[0] );
 339     c[0] = h[imax] / a[0];
 340     double sum = 0;
 341
 342     for ( i = 0; i < imax - 1; i++ )
 343     {
 344         b[i] /= a[i];
 345         if ( i > 0 )
 346             c[i] = - c[i-1] * b[i-1] / a[i];
 347         a[i+1] = qSqrt( a[i+1] - qwtSqr( b[i] ) );
 348         sum += qwtSqr( c[i] );
 349     }
 350     b[imax-1] = ( b[imax-1] - c[imax-2] * b[imax-2] ) / a[imax-1];
 351     a[imax] = qSqrt( a[imax] - qwtSqr( b[imax-1] ) - sum );
 352
 353
 354     // forward elimination
 355     s[0] = d[0] / a[0];
 356     sum = 0;
 357     for ( i = 1; i < imax; i++ )
 358     {
 359         s[i] = ( d[i] - b[i-1] * s[i-1] ) / a[i];
 360         sum += c[i-1] * s[i-1];
 361     }
 362     s[imax] = ( d[imax] - b[imax-1] * s[imax-1] - sum ) / a[imax];
 363
 364
 365     // backward elimination
 366     s[imax] = - s[imax] / a[imax];
 367     s[imax-1] = -( s[imax-1] + b[imax-1] * s[imax] ) / a[imax-1];
 368     for ( i = imax - 2; i >= 0; i-- )
 369         s[i] = - ( s[i] + b[i] * s[i+1] + c[i] * s[imax] ) / a[i];
 370
 371     //
 372     // Finally, determine the spline coefficients
 373     //
 374     s[size-1] = s[0];
 375     for ( i = 0; i < size - 1; i++ )
 376     {
 377         a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i] );
 378         b[i] = 0.5 * s[i];
 379         c[i] = ( p[i+1].y() - p[i].y() )
 380             / h[i] - ( s[i+1] + 2.0 * s[i] ) * h[i] / 6.0;
 381     }
 382
 383     return true;
 384 }