qwt_spline.cpp

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.h"
00011 #include "qwt_math.h"
00012 #include "qwt_array.h"
00013 
00014 class QwtSpline::PrivateData
00015 {
00016 public:
00017     PrivateData():
00018         splineType(QwtSpline::Natural)
00019     {
00020     }
00021 
00022     QwtSpline::SplineType splineType;
00023 
00024     // coefficient vectors
00025     QwtArray<double> a;
00026     QwtArray<double> b;
00027     QwtArray<double> c;
00028 
00029     // control points
00030 #if QT_VERSION < 0x040000
00031     QwtArray<QwtDoublePoint> points;
00032 #else
00033     QPolygonF points;
00034 #endif
00035 };
00036 
00037 #if QT_VERSION < 0x040000
00038 static int lookup(double x, const QwtArray<QwtDoublePoint> &values)
00039 #else
00040 static int lookup(double x, const QPolygonF &values)
00041 #endif
00042 {
00043 #if 0
00044 //qLowerBiund/qHigherBound ???
00045 #endif
00046     int i1;
00047     const int size = (int)values.size();
00048     
00049     if (x <= values[0].x())
00050        i1 = 0;
00051     else if (x >= values[size - 2].x())
00052        i1 = size - 2;
00053     else
00054     {
00055         i1 = 0;
00056         int i2 = size - 2;
00057         int i3 = 0;
00058 
00059         while ( i2 - i1 > 1 )
00060         {
00061             i3 = i1 + ((i2 - i1) >> 1);
00062 
00063             if (values[i3].x() > x)
00064                i2 = i3;
00065             else
00066                i1 = i3;
00067         }
00068     }
00069     return i1;
00070 }
00071 
00073 QwtSpline::QwtSpline()
00074 {
00075     d_data = new PrivateData;
00076 }
00077 
00078 QwtSpline::QwtSpline(const QwtSpline& other)
00079 {
00080     d_data = new PrivateData(*other.d_data);
00081 }
00082 
00083 QwtSpline &QwtSpline::operator=( const QwtSpline &other)
00084 {
00085     *d_data = *other.d_data;
00086     return *this;
00087 }
00088 
00090 QwtSpline::~QwtSpline()
00091 {
00092     delete d_data;
00093 }
00094 
00095 void QwtSpline::setSplineType(SplineType splineType)
00096 {
00097     d_data->splineType = splineType;
00098 }
00099 
00100 QwtSpline::SplineType QwtSpline::splineType() const
00101 {
00102     return d_data->splineType;
00103 }
00104 
00106 
00123 #if QT_VERSION < 0x040000
00124 bool QwtSpline::setPoints(const QwtArray<QwtDoublePoint>& points)
00125 #else
00126 bool QwtSpline::setPoints(const QPolygonF& points)
00127 #endif
00128 {
00129     const int size = points.size();
00130     if (size <= 2) 
00131     {
00132         reset();
00133         return false;
00134     }
00135 
00136 #if QT_VERSION < 0x040000
00137     d_data->points = points.copy(); // Qt3: deep copy
00138 #else
00139     d_data->points = points;
00140 #endif
00141     
00142     d_data->a.resize(size-1);
00143     d_data->b.resize(size-1);
00144     d_data->c.resize(size-1);
00145 
00146     bool ok;
00147     if ( d_data->splineType == Periodic )
00148         ok = buildPeriodicSpline(points);
00149     else
00150         ok = buildNaturalSpline(points);
00151 
00152     if (!ok) 
00153         reset();
00154 
00155     return ok;
00156 }
00157 
00161 #if QT_VERSION < 0x040000
00162 QwtArray<QwtDoublePoint> QwtSpline::points() const
00163 #else
00164 QPolygonF QwtSpline::points() const
00165 #endif
00166 {
00167     return d_data->points;
00168 }
00169 
00170 
00172 void QwtSpline::reset()
00173 {
00174     d_data->a.resize(0);
00175     d_data->b.resize(0);
00176     d_data->c.resize(0);
00177     d_data->points.resize(0);
00178 }
00179 
00181 bool QwtSpline::isValid() const
00182 {
00183     return d_data->a.size() > 0;
00184 }
00185 
00190 double QwtSpline::value(double x) const
00191 {
00192     if (d_data->a.size() == 0)
00193         return 0.0;
00194 
00195     const int i = lookup(x, d_data->points);
00196 
00197     const double delta = x - d_data->points[i].x();
00198     return( ( ( ( d_data->a[i] * delta) + d_data->b[i] ) 
00199         * delta + d_data->c[i] ) * delta + d_data->points[i].y() );
00200 }
00201 
00206 #if QT_VERSION < 0x040000
00207 bool QwtSpline::buildNaturalSpline(const QwtArray<QwtDoublePoint> &points)
00208 #else
00209 bool QwtSpline::buildNaturalSpline(const QPolygonF &points)
00210 #endif
00211 {
00212     int i;
00213     
00214 #if QT_VERSION < 0x040000
00215     const QwtDoublePoint *p = points.data();
00216 #else
00217     const QPointF *p = points.data();
00218 #endif
00219     const int size = points.size();
00220 
00221     double *a = d_data->a.data();
00222     double *b = d_data->b.data();
00223     double *c = d_data->c.data();
00224 
00225     //  set up tridiagonal equation system; use coefficient
00226     //  vectors as temporary buffers
00227     QwtArray<double> h(size-1);
00228     for (i = 0; i < size - 1; i++) 
00229     {
00230         h[i] = p[i+1].x() - p[i].x();
00231         if (h[i] <= 0)
00232             return false;
00233     }
00234     
00235     QwtArray<double> d(size-1);
00236     double dy1 = (p[1].y() - p[0].y()) / h[0];
00237     for (i = 1; i < size - 1; i++)
00238     {
00239         b[i] = c[i] = h[i];
00240         a[i] = 2.0 * (h[i-1] + h[i]);
00241 
00242         const double dy2 = (p[i+1].y() - p[i].y()) / h[i];
00243         d[i] = 6.0 * ( dy1 - dy2);
00244         dy1 = dy2;
00245     }
00246 
00247     //
00248     // solve it
00249     //
00250     
00251     // L-U Factorization
00252     for(i = 1; i < size - 2;i++)
00253     {
00254         c[i] /= a[i];
00255         a[i+1] -= b[i] * c[i]; 
00256     }
00257 
00258     // forward elimination
00259     QwtArray<double> s(size);
00260     s[1] = d[1];
00261     for ( i = 2; i < size - 1; i++)
00262        s[i] = d[i] - c[i-1] * s[i-1];
00263     
00264     // backward elimination
00265     s[size - 2] = - s[size - 2] / a[size - 2];
00266     for (i = size -3; i > 0; i--)
00267        s[i] = - (s[i] + b[i] * s[i+1]) / a[i];
00268     s[size - 1] = s[0] = 0.0;
00269 
00270     //
00271     // Finally, determine the spline coefficients
00272     //
00273     for (i = 0; i < size - 1; i++)
00274     {
00275         a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00276         b[i] = 0.5 * s[i];
00277         c[i] = ( p[i+1].y() - p[i].y() ) / h[i] 
00278             - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0; 
00279     }
00280 
00281     return true;
00282 }
00283 
00288 #if QT_VERSION < 0x040000
00289 bool QwtSpline::buildPeriodicSpline(
00290     const QwtArray<QwtDoublePoint> &points)
00291 #else
00292 bool QwtSpline::buildPeriodicSpline(const QPolygonF &points)
00293 #endif
00294 {
00295     int i;
00296     
00297 #if QT_VERSION < 0x040000
00298     const QwtDoublePoint *p = points.data();
00299 #else
00300     const QPointF *p = points.data();
00301 #endif
00302     const int size = points.size();
00303 
00304     double *a = d_data->a.data();
00305     double *b = d_data->b.data();
00306     double *c = d_data->c.data();
00307 
00308     QwtArray<double> d(size-1);
00309     QwtArray<double> h(size-1);
00310     QwtArray<double> s(size);
00311     
00312     //
00313     //  setup equation system; use coefficient
00314     //  vectors as temporary buffers
00315     //
00316     for (i = 0; i < size - 1; i++)
00317     {
00318         h[i] = p[i+1].x() - p[i].x();
00319         if (h[i] <= 0.0)
00320             return false;
00321     }
00322     
00323     const int imax = size - 2;
00324     double htmp = h[imax];
00325     double dy1 = (p[0].y() - p[imax].y()) / htmp;
00326     for (i = 0; i <= imax; i++)
00327     {
00328         b[i] = c[i] = h[i];
00329         a[i] = 2.0 * (htmp + h[i]);
00330         const double dy2 = (p[i+1].y() - p[i].y()) / h[i];
00331         d[i] = 6.0 * ( dy1 - dy2);
00332         dy1 = dy2;
00333         htmp = h[i];
00334     }
00335 
00336     //
00337     // solve it
00338     //
00339     
00340     // L-U Factorization
00341     a[0] = sqrt(a[0]);
00342     c[0] = h[imax] / a[0];
00343     double sum = 0;
00344 
00345     for( i = 0; i < imax - 1; i++)
00346     {
00347         b[i] /= a[i];
00348         if (i > 0)
00349            c[i] = - c[i-1] * b[i-1] / a[i];
00350         a[i+1] = sqrt( a[i+1] - qwtSqr(b[i]));
00351         sum += qwtSqr(c[i]);
00352     }
00353     b[imax-1] = (b[imax-1] - c[imax-2] * b[imax-2]) / a[imax-1];
00354     a[imax] = sqrt(a[imax] - qwtSqr(b[imax-1]) - sum);
00355     
00356 
00357     // forward elimination
00358     s[0] = d[0] / a[0];
00359     sum = 0;
00360     for( i = 1; i < imax; i++)
00361     {
00362         s[i] = (d[i] - b[i-1] * s[i-1]) / a[i];
00363         sum += c[i-1] * s[i-1];
00364     }
00365     s[imax] = (d[imax] - b[imax-1] * s[imax-1] - sum) / a[imax];
00366     
00367     
00368     // backward elimination
00369     s[imax] = - s[imax] / a[imax];
00370     s[imax-1] = -(s[imax-1] + b[imax-1] * s[imax]) / a[imax-1];
00371     for (i= imax - 2; i >= 0; i--)
00372        s[i] = - (s[i] + b[i] * s[i+1] + c[i] * s[imax]) / a[i];
00373 
00374     //
00375     // Finally, determine the spline coefficients
00376     //
00377     s[size-1] = s[0];
00378     for ( i=0; i < size-1; i++)
00379     {
00380         a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00381         b[i] = 0.5 * s[i];
00382         c[i] = ( p[i+1].y() - p[i].y() ) 
00383             / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0; 
00384     }
00385 
00386     return true;
00387 }

Generated on Mon Feb 26 21:22:38 2007 for Qwt User's Guide by  doxygen 1.4.6