source: ntrip/branches/BNC_2.12/qwt/qwt_spline.cpp@ 8142

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