← Back to branch summary

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

« back to all changes in this revision

Viewing changes to qwt-5.1.0/src/qwt_spline.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Fathi Boudra
  • Date: 2008-05-26 10:26:31 UTC
  • mfrom: (1.1.3 upstream) (2.1.1 lenny)
  • Revision ID: james.westby@ubuntu.com-20080526102631-bp95mfccnrb957nx
Tags: 5.1.1-1
New upstream release.

Show diffs side-by-side

added added

removed removed

Lines of Context:
qwt-5.1.0/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
 
QwtSpline::QwtSpline(const QwtSpline& other)
79
 
{
80
 
    d_data = new PrivateData(*other.d_data);
81
 
}
82
 
 
83
 
QwtSpline &QwtSpline::operator=( const QwtSpline &other)
84
 
{
85
 
    *d_data = *other.d_data;
86
 
    return *this;
87
 
}
88
 
 
89
 
//! Destructor
90
 
QwtSpline::~QwtSpline()
91
 
{
92
 
    delete d_data;
93
 
}
94
 
 
95
 
void QwtSpline::setSplineType(SplineType splineType)
96
 
{
97
 
    d_data->splineType = splineType;
98
 
}
99
 
 
100
 
QwtSpline::SplineType QwtSpline::splineType() const
101
 
{
102
 
    return d_data->splineType;
103
 
}
104
 
 
105
 
//! Determine the function table index corresponding to a value x
106
 
 
107
 
/*!
108
 
  \brief Calculate the spline coefficients
109
 
 
110
 
  Depending on the value of \a periodic, this function
111
 
  will determine the coefficients for a natural or a periodic
112
 
  spline and store them internally. 
113
 
  
114
 
  \param x
115
 
  \param y points
116
 
  \param size number of points
117
 
  \param periodic if true, calculate periodic spline
118
 
  \return true if successful
119
 
  \warning The sequence of x (but not y) values has to be strictly monotone
120
 
           increasing, which means <code>x[0] < x[1] < .... < x[n-1]</code>.
121
 
       If this is not the case, the function will return false
122
 
*/
123
 
#if QT_VERSION < 0x040000
124
 
bool QwtSpline::setPoints(const QwtArray<QwtDoublePoint>& points)
125
 
#else
126
 
bool QwtSpline::setPoints(const QPolygonF& points)
127
 
#endif
128
 
{
129
 
    const int size = points.size();
130
 
    if (size <= 2) 
131
 
    {
132
 
        reset();
133
 
        return false;
134
 
    }
135
 
 
136
 
#if QT_VERSION < 0x040000
137
 
    d_data->points = points.copy(); // Qt3: deep copy
138
 
#else
139
 
    d_data->points = points;
140
 
#endif
141
 
    
142
 
    d_data->a.resize(size-1);
143
 
    d_data->b.resize(size-1);
144
 
    d_data->c.resize(size-1);
145
 
 
146
 
    bool ok;
147
 
    if ( d_data->splineType == Periodic )
148
 
        ok = buildPeriodicSpline(points);
149
 
    else
150
 
        ok = buildNaturalSpline(points);
151
 
 
152
 
    if (!ok) 
153
 
        reset();
154
 
 
155
 
    return ok;
156
 
}
157
 
 
158
 
/*!
159
 
   Return points passed by setPoints
160
 
*/
161
 
#if QT_VERSION < 0x040000
162
 
QwtArray<QwtDoublePoint> QwtSpline::points() const
163
 
#else
164
 
QPolygonF QwtSpline::points() const
165
 
#endif
166
 
{
167
 
    return d_data->points;
168
 
}
169
 
 
170
 
 
171
 
//! Free allocated memory and set size to 0
172
 
void QwtSpline::reset()
173
 
{
174
 
    d_data->a.resize(0);
175
 
    d_data->b.resize(0);
176
 
    d_data->c.resize(0);
177
 
    d_data->points.resize(0);
178
 
}
179
 
 
180
 
//! True if valid
181
 
bool QwtSpline::isValid() const
182
 
{
183
 
    return d_data->a.size() > 0;
184
 
}
185
 
 
186
 
/*!
187
 
  Calculate the interpolated function value corresponding 
188
 
  to a given argument x.
189
 
*/
190
 
double QwtSpline::value(double x) const
191
 
{
192
 
    if (d_data->a.size() == 0)
193
 
        return 0.0;
194
 
 
195
 
    const int i = lookup(x, d_data->points);
196
 
 
197
 
    const double delta = x - d_data->points[i].x();
198
 
    return( ( ( ( d_data->a[i] * delta) + d_data->b[i] ) 
199
 
        * delta + d_data->c[i] ) * delta + d_data->points[i].y() );
200
 
}
201
 
 
202
 
/*!
203
 
  \brief Determines the coefficients for a natural spline
204
 
  \return true if successful
205
 
*/
206
 
#if QT_VERSION < 0x040000
207
 
bool QwtSpline::buildNaturalSpline(const QwtArray<QwtDoublePoint> &points)
208
 
#else
209
 
bool QwtSpline::buildNaturalSpline(const QPolygonF &points)
210
 
#endif
211
 
{
212
 
    int i;
213
 
    
214
 
#if QT_VERSION < 0x040000
215
 
    const QwtDoublePoint *p = points.data();
216
 
#else
217
 
    const QPointF *p = points.data();
218
 
#endif
219
 
    const int size = points.size();
220
 
 
221
 
    double *a = d_data->a.data();
222
 
    double *b = d_data->b.data();
223
 
    double *c = d_data->c.data();
224
 
 
225
 
    //  set up tridiagonal equation system; use coefficient
226
 
    //  vectors as temporary buffers
227
 
    QwtArray<double> h(size-1);
228
 
    for (i = 0; i < size - 1; i++) 
229
 
    {
230
 
        h[i] = p[i+1].x() - p[i].x();
231
 
        if (h[i] <= 0)
232
 
            return false;
233
 
    }
234
 
    
235
 
    QwtArray<double> d(size-1);
236
 
    double dy1 = (p[1].y() - p[0].y()) / h[0];
237
 
    for (i = 1; i < size - 1; i++)
238
 
    {
239
 
        b[i] = c[i] = h[i];
240
 
        a[i] = 2.0 * (h[i-1] + h[i]);
241
 
 
242
 
        const double dy2 = (p[i+1].y() - p[i].y()) / h[i];
243
 
        d[i] = 6.0 * ( dy1 - dy2);
244
 
        dy1 = dy2;
245
 
    }
246
 
 
247
 
    //
248
 
    // solve it
249
 
    //
250
 
    
251
 
    // L-U Factorization
252
 
    for(i = 1; i < size - 2;i++)
253
 
    {
254
 
        c[i] /= a[i];
255
 
        a[i+1] -= b[i] * c[i]; 
256
 
    }
257
 
 
258
 
    // forward elimination
259
 
    QwtArray<double> s(size);
260
 
    s[1] = d[1];
261
 
    for ( i = 2; i < size - 1; i++)
262
 
       s[i] = d[i] - c[i-1] * s[i-1];
263
 
    
264
 
    // backward elimination
265
 
    s[size - 2] = - s[size - 2] / a[size - 2];
266
 
    for (i = size -3; i > 0; i--)
267
 
       s[i] = - (s[i] + b[i] * s[i+1]) / a[i];
268
 
    s[size - 1] = s[0] = 0.0;
269
 
 
270
 
    //
271
 
    // Finally, determine the spline coefficients
272
 
    //
273
 
    for (i = 0; i < size - 1; i++)
274
 
    {
275
 
        a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
276
 
        b[i] = 0.5 * s[i];
277
 
        c[i] = ( p[i+1].y() - p[i].y() ) / h[i] 
278
 
            - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0; 
279
 
    }
280
 
 
281
 
    return true;
282
 
}
283
 
 
284
 
/*!
285
 
  \brief Determines the coefficients for a periodic spline
286
 
  \return true if successful
287
 
*/
288
 
#if QT_VERSION < 0x040000
289
 
bool QwtSpline::buildPeriodicSpline(
290
 
    const QwtArray<QwtDoublePoint> &points)
291
 
#else
292
 
bool QwtSpline::buildPeriodicSpline(const QPolygonF &points)
293
 
#endif
294
 
{
295
 
    int i;
296
 
    
297
 
#if QT_VERSION < 0x040000
298
 
    const QwtDoublePoint *p = points.data();
299
 
#else
300
 
    const QPointF *p = points.data();
301
 
#endif
302
 
    const int size = points.size();
303
 
 
304
 
    double *a = d_data->a.data();
305
 
    double *b = d_data->b.data();
306
 
    double *c = d_data->c.data();
307
 
 
308
 
    QwtArray<double> d(size-1);
309
 
    QwtArray<double> h(size-1);
310
 
    QwtArray<double> s(size);
311
 
    
312
 
    //
313
 
    //  setup equation system; use coefficient
314
 
    //  vectors as temporary buffers
315
 
    //
316
 
    for (i = 0; i < size - 1; i++)
317
 
    {
318
 
        h[i] = p[i+1].x() - p[i].x();
319
 
        if (h[i] <= 0.0)
320
 
            return false;
321
 
    }
322
 
    
323
 
    const int imax = size - 2;
324
 
    double htmp = h[imax];
325
 
    double dy1 = (p[0].y() - p[imax].y()) / htmp;
326
 
    for (i = 0; i <= imax; i++)
327
 
    {
328
 
        b[i] = c[i] = h[i];
329
 
        a[i] = 2.0 * (htmp + h[i]);
330
 
        const double dy2 = (p[i+1].y() - p[i].y()) / h[i];
331
 
        d[i] = 6.0 * ( dy1 - dy2);
332
 
        dy1 = dy2;
333
 
        htmp = h[i];
334
 
    }
335
 
 
336
 
    //
337
 
    // solve it
338
 
    //
339
 
    
340
 
    // L-U Factorization
341
 
    a[0] = sqrt(a[0]);
342
 
    c[0] = h[imax] / a[0];
343
 
    double sum = 0;
344
 
 
345
 
    for( i = 0; i < imax - 1; i++)
346
 
    {
347
 
        b[i] /= a[i];
348
 
        if (i > 0)
349
 
           c[i] = - c[i-1] * b[i-1] / a[i];
350
 
        a[i+1] = sqrt( a[i+1] - qwtSqr(b[i]));
351
 
        sum += qwtSqr(c[i]);
352
 
    }
353
 
    b[imax-1] = (b[imax-1] - c[imax-2] * b[imax-2]) / a[imax-1];
354
 
    a[imax] = sqrt(a[imax] - qwtSqr(b[imax-1]) - sum);
355
 
    
356
 
 
357
 
    // forward elimination
358
 
    s[0] = d[0] / a[0];
359
 
    sum = 0;
360
 
    for( i = 1; i < imax; i++)
361
 
    {
362
 
        s[i] = (d[i] - b[i-1] * s[i-1]) / a[i];
363
 
        sum += c[i-1] * s[i-1];
364
 
    }
365
 
    s[imax] = (d[imax] - b[imax-1] * s[imax-1] - sum) / a[imax];
366
 
    
367
 
    
368
 
    // backward elimination
369
 
    s[imax] = - s[imax] / a[imax];
370
 
    s[imax-1] = -(s[imax-1] + b[imax-1] * s[imax]) / a[imax-1];
371
 
    for (i= imax - 2; i >= 0; i--)
372
 
       s[i] = - (s[i] + b[i] * s[i+1] + c[i] * s[imax]) / a[i];
373
 
 
374
 
    //
375
 
    // Finally, determine the spline coefficients
376
 
    //
377
 
    s[size-1] = s[0];
378
 
    for ( i=0; i < size-1; i++)
379
 
    {
380
 
        a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
381
 
        b[i] = 0.5 * s[i];
382
 
        c[i] = ( p[i+1].y() - p[i].y() ) 
383
 
            / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0; 
384
 
    }
385
 
 
386
 
    return true;
387
 
}