00001
00002
00003
00004
00005
00006
00007
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
00025 QwtArray<double> a;
00026 QwtArray<double> b;
00027 QwtArray<double> c;
00028
00029
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
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();
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
00158
00160 void QwtSpline::reset()
00161 {
00162 d_data->a.resize(0);
00163 d_data->b.resize(0);
00164 d_data->c.resize(0);
00165 d_data->points.resize(0);
00166 }
00167
00168 bool QwtSpline::isValid() const
00169 {
00170 return d_data->a.size() > 0;
00171 }
00172
00177 double QwtSpline::value(double x) const
00178 {
00179 if (d_data->a.size() == 0)
00180 return 0.0;
00181
00182 const int i = lookup(x, d_data->points);
00183
00184 const double delta = x - d_data->points[i].x();
00185 return( ( ( ( d_data->a[i] * delta) + d_data->b[i] )
00186 * delta + d_data->c[i] ) * delta + d_data->points[i].y() );
00187 }
00188
00193 #if QT_VERSION < 0x040000
00194 bool QwtSpline::buildNaturalSpline(const QwtArray<QwtDoublePoint> &points)
00195 #else
00196 bool QwtSpline::buildNaturalSpline(const QPolygonF &points)
00197 #endif
00198 {
00199 int i;
00200
00201 #if QT_VERSION < 0x040000
00202 const QwtDoublePoint *p = points.data();
00203 #else
00204 const QPointF *p = points.data();
00205 #endif
00206 const int size = points.size();
00207
00208 double *a = d_data->a.data();
00209 double *b = d_data->b.data();
00210 double *c = d_data->c.data();
00211
00212
00213
00214 QwtArray<double> h(size-1);
00215 for (i = 0; i < size - 1; i++)
00216 {
00217 h[i] = p[i+1].x() - p[i].x();
00218 if (h[i] <= 0)
00219 return false;
00220 }
00221
00222 QwtArray<double> d(size-1);
00223 double dy1 = (p[1].y() - p[0].y()) / h[0];
00224 for (i = 1; i < size - 1; i++)
00225 {
00226 b[i] = c[i] = h[i];
00227 a[i] = 2.0 * (h[i-1] + h[i]);
00228
00229 const double dy2 = (p[i+1].y() - p[i].y()) / h[i];
00230 d[i] = 6.0 * ( dy1 - dy2);
00231 dy1 = dy2;
00232 }
00233
00234
00235
00236
00237
00238
00239 for(i = 1; i < size - 2;i++)
00240 {
00241 c[i] /= a[i];
00242 a[i+1] -= b[i] * c[i];
00243 }
00244
00245
00246 QwtArray<double> s(size);
00247 s[1] = d[1];
00248 for ( i = 2; i < size - 1; i++)
00249 s[i] = d[i] - c[i-1] * s[i-1];
00250
00251
00252 s[size - 2] = - s[size - 2] / a[size - 2];
00253 for (i = size -3; i > 0; i--)
00254 s[i] = - (s[i] + b[i] * s[i+1]) / a[i];
00255 s[size - 1] = s[0] = 0.0;
00256
00257
00258
00259
00260 for (i = 0; i < size - 1; i++)
00261 {
00262 a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00263 b[i] = 0.5 * s[i];
00264 c[i] = ( p[i+1].y() - p[i].y() ) / h[i]
00265 - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0;
00266 }
00267
00268 return true;
00269 }
00270
00275 #if QT_VERSION < 0x040000
00276 bool QwtSpline::buildPeriodicSpline(
00277 const QwtArray<QwtDoublePoint> &points)
00278 #else
00279 bool QwtSpline::buildPeriodicSpline(const QPolygonF &points)
00280 #endif
00281 {
00282 int i;
00283
00284 #if QT_VERSION < 0x040000
00285 const QwtDoublePoint *p = points.data();
00286 #else
00287 const QPointF *p = points.data();
00288 #endif
00289 const int size = points.size();
00290
00291 double *a = d_data->a.data();
00292 double *b = d_data->b.data();
00293 double *c = d_data->c.data();
00294
00295 QwtArray<double> d(size-1);
00296 QwtArray<double> h(size-1);
00297 QwtArray<double> s(size);
00298
00299
00300
00301
00302
00303 for (i = 0; i < size - 1; i++)
00304 {
00305 h[i] = p[i+1].x() - p[i].x();
00306 if (h[i] <= 0.0)
00307 return false;
00308 }
00309
00310 const int imax = size - 2;
00311 double htmp = h[imax];
00312 double dy1 = (p[0].y() - p[imax].y()) / htmp;
00313 for (i = 0; i <= imax; i++)
00314 {
00315 b[i] = c[i] = h[i];
00316 a[i] = 2.0 * (htmp + h[i]);
00317 const double dy2 = (p[i+1].y() - p[i].y()) / h[i];
00318 d[i] = 6.0 * ( dy1 - dy2);
00319 dy1 = dy2;
00320 htmp = h[i];
00321 }
00322
00323
00324
00325
00326
00327
00328 a[0] = sqrt(a[0]);
00329 c[0] = h[imax] / a[0];
00330 double sum = 0;
00331
00332 for( i = 0; i < imax - 1; i++)
00333 {
00334 b[i] /= a[i];
00335 if (i > 0)
00336 c[i] = - c[i-1] * b[i-1] / a[i];
00337 a[i+1] = sqrt( a[i+1] - qwtSqr(b[i]));
00338 sum += qwtSqr(c[i]);
00339 }
00340 b[imax-1] = (b[imax-1] - c[imax-2] * b[imax-2]) / a[imax-1];
00341 a[imax] = sqrt(a[imax] - qwtSqr(b[imax-1]) - sum);
00342
00343
00344
00345 s[0] = d[0] / a[0];
00346 sum = 0;
00347 for( i = 1; i < imax; i++)
00348 {
00349 s[i] = (d[i] - b[i-1] * s[i-1]) / a[i];
00350 sum += c[i-1] * s[i-1];
00351 }
00352 s[imax] = (d[imax] - b[imax-1] * s[imax-1] - sum) / a[imax];
00353
00354
00355
00356 s[imax] = - s[imax] / a[imax];
00357 s[imax-1] = -(s[imax-1] + b[imax-1] * s[imax]) / a[imax-1];
00358 for (i= imax - 2; i >= 0; i--)
00359 s[i] = - (s[i] + b[i] * s[i+1] + c[i] * s[imax]) / a[i];
00360
00361
00362
00363
00364 s[size-1] = s[0];
00365 for ( i=0; i < size-1; i++)
00366 {
00367 a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00368 b[i] = 0.5 * s[i];
00369 c[i] = ( p[i+1].y() - p[i].y() )
00370 / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0;
00371 }
00372
00373 return true;
00374 }