Main Page | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Class Members | File Members | Related Pages

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 
00013 class QwtSpline::PrivateData
00014 {
00015 public:
00016     PrivateData():
00017         size(0),
00018         buffered(false)
00019     {
00020         a = b = c = NULL;
00021         xbuffer = ybuffer = x = y = NULL;
00022     }
00023 
00024     // coefficient vectors
00025     double *a;
00026     double *b;
00027     double *c;
00028     double *d;
00029 
00030     // values
00031     double *x;
00032     double *y;
00033     double *xbuffer;
00034     double *ybuffer;
00035     int size;
00036 
00037     //flags
00038     bool buffered;
00039 };
00040 
00042 QwtSpline::QwtSpline()
00043 {
00044     d_data = new PrivateData;
00045 }
00046 
00048 QwtSpline::~QwtSpline()
00049 {
00050     cleanup();
00051     delete d_data;
00052 }
00053 
00054 
00077 void QwtSpline::copyValues(bool tf)
00078 {
00079     cleanup();
00080     d_data->buffered = tf;
00081 }
00082 
00087 double QwtSpline::value(double x) const
00088 {
00089     if (!d_data->a)
00090         return 0.0;
00091 
00092     const int i = lookup(x);
00093 
00094     const double delta = x - d_data->x[i];
00095     return( ( ( ( d_data->a[i] * delta) + d_data->b[i] ) 
00096         * delta + d_data->c[i] ) * delta + d_data->y[i] );
00097 }
00098 
00100 int QwtSpline::lookup(double x) const
00101 {
00102     int i1, i2, i3;
00103     
00104     if (x <= d_data->x[0])
00105        i1 = 0;
00106     else if (x >= d_data->x[d_data->size - 2])
00107        i1 = d_data->size -2;
00108     else
00109     {
00110         i1 = 0;
00111         i2 = d_data->size -2;
00112         i3 = 0;
00113 
00114         while ( i2 - i1 > 1 )
00115         {
00116             i3 = i1 + ((i2 - i1) >> 1);
00117 
00118             if (d_data->x[i3] > x)
00119                i2 = i3;
00120             else
00121                i1 = i3;
00122 
00123         }
00124     }
00125     return i1;
00126 }
00127 
00128 
00148 bool QwtSpline::recalc(double *x, double *y, int n, bool periodic)
00149 {
00150     cleanup();
00151 
00152     if (n <= 2)
00153         return false;
00154 
00155     d_data->size = n;
00156 
00157     if (d_data->buffered)
00158     {
00159         d_data->xbuffer = new double[n];
00160         d_data->ybuffer = new double[n];
00161 
00162         for (int i = 0; i < n; i++)
00163         {
00164             d_data->xbuffer[i] = x[i];
00165             d_data->ybuffer[i] = y[i];
00166         }
00167         d_data->x = d_data->xbuffer;
00168         d_data->y = d_data->ybuffer;
00169     }
00170     else
00171     {
00172         d_data->x = x;
00173         d_data->y = y;
00174     }
00175     
00176     d_data->a = new double[n-1];
00177     d_data->b = new double[n-1];
00178     d_data->c = new double[n-1];
00179 
00180     bool ok;
00181     if(periodic)
00182        ok =  buildPerSpline();
00183     else
00184        ok =  buildNatSpline();
00185 
00186     if (!ok) 
00187         cleanup();
00188 
00189     return ok;
00190 }
00191 
00196 bool QwtSpline::buildNatSpline()
00197 {
00198     int i;
00199     
00200     double *d = new double[d_data->size-1];
00201     double *h = new double[d_data->size-1];
00202     double *s = new double[d_data->size];
00203 
00204     //
00205     //  set up tridiagonal equation system; use coefficient
00206     //  vectors as temporary buffers
00207     for (i = 0; i < d_data->size - 1; i++) 
00208     {
00209         h[i] = d_data->x[i+1] - d_data->x[i];
00210         if (h[i] <= 0)
00211         {
00212             delete[] h;
00213             delete[] s;
00214             delete[] d;
00215             return false;
00216         }
00217     }
00218     
00219     double dy1 = (d_data->y[1] - d_data->y[0]) / h[0];
00220     for (i = 1; i < d_data->size - 1; i++)
00221     {
00222         d_data->b[i] = d_data->c[i] = h[i];
00223         d_data->a[i] = 2.0 * (h[i-1] + h[i]);
00224 
00225         const double dy2 = (d_data->y[i+1] - d_data->y[i]) / h[i];
00226         d[i] = 6.0 * ( dy1 - dy2);
00227         dy1 = dy2;
00228     }
00229 
00230     //
00231     // solve it
00232     //
00233     
00234     // L-U Factorization
00235     for(i = 1; i < d_data->size - 2;i++)
00236     {
00237         d_data->c[i] /= d_data->a[i];
00238         d_data->a[i+1] -= d_data->b[i] * d_data->c[i]; 
00239     }
00240 
00241     // forward elimination
00242     s[1] = d[1];
00243     for(i=2;i<d_data->size - 1;i++)
00244        s[i] = d[i] - d_data->c[i-1] * s[i-1];
00245     
00246     // backward elimination
00247     s[d_data->size - 2] = - s[d_data->size - 2] / d_data->a[d_data->size - 2];
00248     for (i= d_data->size -3; i > 0; i--)
00249        s[i] = - (s[i] + d_data->b[i] * s[i+1]) / d_data->a[i];
00250 
00251     //
00252     // Finally, determine the spline coefficients
00253     //
00254     s[d_data->size - 1] = s[0] = 0.0;
00255     for (i = 0; i < d_data->size - 1; i++)
00256     {
00257         d_data->a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00258         d_data->b[i] = 0.5 * s[i];
00259         d_data->c[i] = ( d_data->y[i+1] - d_data->y[i] ) 
00260             / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0; 
00261     }
00262 
00263     delete[] d;
00264     delete[] s;
00265     delete[] h;
00266 
00267     return true;
00268 }
00269 
00274 bool QwtSpline::buildPerSpline()
00275 {
00276     int i;
00277     
00278     double *d = new double[d_data->size-1];
00279     double *h = new double[d_data->size-1];
00280     double *s = new double[d_data->size];
00281     
00282     //
00283     //  setup equation system; use coefficient
00284     //  vectors as temporary buffers
00285     //
00286     for (i=0; i<d_data->size - 1; i++)
00287     {
00288         h[i] = d_data->x[i+1] - d_data->x[i];
00289         if (h[i] <= 0.0)
00290         {
00291             delete[] h;
00292             delete[] s;
00293             delete[] d;
00294             return false;
00295         }
00296     }
00297     
00298     const int imax = d_data->size - 2;
00299     double htmp = h[imax];
00300     double dy1 = (d_data->y[0] - d_data->y[imax]) / htmp;
00301     for (i=0; i <= imax; i++)
00302     {
00303         d_data->b[i] = d_data->c[i] = h[i];
00304         d_data->a[i] = 2.0 * (htmp + h[i]);
00305         const double dy2 = (d_data->y[i+1] - d_data->y[i]) / h[i];
00306         d[i] = 6.0 * ( dy1 - dy2);
00307         dy1 = dy2;
00308         htmp = h[i];
00309     }
00310 
00311     //
00312     // solve it
00313     //
00314     
00315     // L-U Factorization
00316     d_data->a[0] = sqrt(d_data->a[0]);
00317     d_data->c[0] = h[imax] / d_data->a[0];
00318     double sum = 0;
00319 
00320     for(i=0;i<imax-1;i++)
00321     {
00322         d_data->b[i] /= d_data->a[i];
00323         if (i > 0)
00324            d_data->c[i] = - d_data->c[i-1] * d_data->b[i-1] / d_data->a[i];
00325         d_data->a[i+1] = sqrt( d_data->a[i+1] - qwtSqr(d_data->b[i]));
00326         sum += qwtSqr(d_data->c[i]);
00327     }
00328     d_data->b[imax-1] = (d_data->b[imax-1] - d_data->c[imax-2] * d_data->b[imax-2]) / d_data->a[imax-1];
00329     d_data->a[imax] = sqrt(d_data->a[imax] - qwtSqr(d_data->b[imax-1]) - sum);
00330     
00331 
00332     // forward elimination
00333     s[0] = d[0] / d_data->a[0];
00334     sum = 0;
00335     for(i=1;i<imax;i++)
00336     {
00337         s[i] = (d[i] - d_data->b[i-1] * s[i-1]) / d_data->a[i];
00338         sum += d_data->c[i-1] * s[i-1];
00339     }
00340     s[imax] = (d[imax] - d_data->b[imax-1]*s[imax-1] - sum) / d_data->a[imax];
00341     
00342     
00343     // backward elimination
00344     s[imax] = - s[imax] / d_data->a[imax];
00345     s[imax-1] = -(s[imax-1] + d_data->b[imax-1] * s[imax]) / d_data->a[imax-1];
00346     for (i= imax - 2; i >= 0; i--)
00347        s[i] = - (s[i] + d_data->b[i] * s[i+1] + d_data->c[i] * s[imax]) / d_data->a[i];
00348 
00349     //
00350     // Finally, determine the spline coefficients
00351     //
00352     s[d_data->size-1] = s[0];
00353     for (i=0;i<d_data->size-1;i++)
00354     {
00355         d_data->a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00356         d_data->b[i] = 0.5 * s[i];
00357         d_data->c[i] = ( d_data->y[i+1] - d_data->y[i] ) 
00358             / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0; 
00359     }
00360 
00361     delete[] d;
00362     delete[] s;
00363     delete[] h;
00364 
00365     return true;
00366 }
00367 
00368 
00370 void QwtSpline::cleanup()
00371 {
00372     delete[] d_data->a;
00373     delete[] d_data->b;
00374     delete[] d_data->c;
00375     delete[] d_data->xbuffer;
00376     delete[] d_data->ybuffer;
00377 
00378     d_data->a = d_data->b = d_data->c = NULL;
00379     d_data->xbuffer = d_data->ybuffer = d_data->x = d_data->y = NULL;
00380     d_data->size = 0;
00381 }

Generated on Mon Jan 30 22:16:26 2006 for Qwt User's Guide by  doxygen 1.4.4