← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/qwt/oneiric-proposed

« back to all changes in this revision

Viewing changes to qwt/src/qwt_spline.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Fathi Boudra
  • Date: 2009-07-01 16:08:21 UTC
  • mfrom: (1.1.5 upstream)
  • Revision ID: james.westby@ubuntu.com-20090701160821-gog44m4w7x2f4u6o
Tags: 5.2.0-1
* New upstream release.
* Add branch pull patch from svn r544.
* Bump Standards-Version to 3.8.2. No changes needed.
* Update installed manpages.
* Use qwt as base name for directory (drop version).
* Add missing warnings patch.
* Rewrite debian/rules: use dh.

Show diffs side-by-side

added added

removed removed

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