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
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
00226
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
00249
00250
00251
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
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
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
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
00314
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
00338
00339
00340
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
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
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
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 }