source: ntrip/trunk/BNC/qwt/qwt_curve_fitter.cpp @ 8127

Last change on this file since 8127 was 8127, checked in by stoecker, 2 years ago

update qwt and qwtpolar, many QT5 fixes (unfinished)

File size: 10.1 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_curve_fitter.h"
11#include "qwt_math.h"
12#include "qwt_spline.h"
13#include <qstack.h>
14#include <qvector.h>
15
16#if QT_VERSION < 0x040601
17#define qFabs(x) ::fabs(x)
18#endif
19
20//! Constructor
21QwtCurveFitter::QwtCurveFitter()
22{
23}
24
25//! Destructor
26QwtCurveFitter::~QwtCurveFitter()
27{
28}
29
30class QwtSplineCurveFitter::PrivateData
31{
32public:
33    PrivateData():
34        fitMode( QwtSplineCurveFitter::Auto ),
35        splineSize( 250 )
36    {
37    }
38
39    QwtSpline spline;
40    QwtSplineCurveFitter::FitMode fitMode;
41    int splineSize;
42};
43
44//! Constructor
45QwtSplineCurveFitter::QwtSplineCurveFitter()
46{
47    d_data = new PrivateData;
48}
49
50//! Destructor
51QwtSplineCurveFitter::~QwtSplineCurveFitter()
52{
53    delete d_data;
54}
55
56/*!
57  Select the algorithm used for building the spline
58
59  \param mode Mode representing a spline algorithm
60  \sa fitMode()
61*/
62void QwtSplineCurveFitter::setFitMode( FitMode mode )
63{
64    d_data->fitMode = mode;
65}
66
67/*!
68  \return Mode representing a spline algorithm
69  \sa setFitMode()
70*/
71QwtSplineCurveFitter::FitMode QwtSplineCurveFitter::fitMode() const
72{
73    return d_data->fitMode;
74}
75
76/*!
77  Assign a spline
78
79  \param spline Spline
80  \sa spline()
81*/
82void QwtSplineCurveFitter::setSpline( const QwtSpline &spline )
83{
84    d_data->spline = spline;
85    d_data->spline.reset();
86}
87
88/*!
89  \return Spline
90  \sa setSpline()
91*/
92const QwtSpline &QwtSplineCurveFitter::spline() const
93{
94    return d_data->spline;
95}
96
97/*!
98  \return Spline
99  \sa setSpline()
100*/
101QwtSpline &QwtSplineCurveFitter::spline()
102{
103    return d_data->spline;
104}
105
106/*!
107   Assign a spline size ( has to be at least 10 points )
108
109   \param splineSize Spline size
110   \sa splineSize()
111*/
112void QwtSplineCurveFitter::setSplineSize( int splineSize )
113{
114    d_data->splineSize = qMax( splineSize, 10 );
115}
116
117/*!
118  \return Spline size
119  \sa setSplineSize()
120*/
121int QwtSplineCurveFitter::splineSize() const
122{
123    return d_data->splineSize;
124}
125
126/*!
127  Find a curve which has the best fit to a series of data points
128
129  \param points Series of data points
130  \return Curve points
131*/
132QPolygonF QwtSplineCurveFitter::fitCurve( const QPolygonF &points ) const
133{
134    const int size = points.size();
135    if ( size <= 2 )
136        return points;
137
138    FitMode fitMode = d_data->fitMode;
139    if ( fitMode == Auto )
140    {
141        fitMode = Spline;
142
143        const QPointF *p = points.data();
144        for ( int i = 1; i < size; i++ )
145        {
146            if ( p[i].x() <= p[i-1].x() )
147            {
148                fitMode = ParametricSpline;
149                break;
150            }
151        };
152    }
153
154    if ( fitMode == ParametricSpline )
155        return fitParametric( points );
156    else
157        return fitSpline( points );
158}
159
160QPolygonF QwtSplineCurveFitter::fitSpline( const QPolygonF &points ) const
161{
162    d_data->spline.setPoints( points );
163    if ( !d_data->spline.isValid() )
164        return points;
165
166    QPolygonF fittedPoints( d_data->splineSize );
167
168    const double x1 = points[0].x();
169    const double x2 = points[int( points.size() - 1 )].x();
170    const double dx = x2 - x1;
171    const double delta = dx / ( d_data->splineSize - 1 );
172
173    for ( int i = 0; i < d_data->splineSize; i++ )
174    {
175        QPointF &p = fittedPoints[i];
176
177        const double v = x1 + i * delta;
178        const double sv = d_data->spline.value( v );
179
180        p.setX( v );
181        p.setY( sv );
182    }
183    d_data->spline.reset();
184
185    return fittedPoints;
186}
187
188QPolygonF QwtSplineCurveFitter::fitParametric( const QPolygonF &points ) const
189{
190    int i;
191    const int size = points.size();
192
193    QPolygonF fittedPoints( d_data->splineSize );
194    QPolygonF splinePointsX( size );
195    QPolygonF splinePointsY( size );
196
197    const QPointF *p = points.data();
198    QPointF *spX = splinePointsX.data();
199    QPointF *spY = splinePointsY.data();
200
201    double param = 0.0;
202    for ( i = 0; i < size; i++ )
203    {
204        const double x = p[i].x();
205        const double y = p[i].y();
206        if ( i > 0 )
207        {
208            const double delta = qSqrt( qwtSqr( x - spX[i-1].y() )
209                      + qwtSqr( y - spY[i-1].y() ) );
210            param += qMax( delta, 1.0 );
211        }
212        spX[i].setX( param );
213        spX[i].setY( x );
214        spY[i].setX( param );
215        spY[i].setY( y );
216    }
217
218    d_data->spline.setPoints( splinePointsX );
219    if ( !d_data->spline.isValid() )
220        return points;
221
222    const double deltaX =
223        splinePointsX[size - 1].x() / ( d_data->splineSize - 1 );
224    for ( i = 0; i < d_data->splineSize; i++ )
225    {
226        const double dtmp = i * deltaX;
227        fittedPoints[i].setX( d_data->spline.value( dtmp ) );
228    }
229
230    d_data->spline.setPoints( splinePointsY );
231    if ( !d_data->spline.isValid() )
232        return points;
233
234    const double deltaY =
235        splinePointsY[size - 1].x() / ( d_data->splineSize - 1 );
236    for ( i = 0; i < d_data->splineSize; i++ )
237    {
238        const double dtmp = i * deltaY;
239        fittedPoints[i].setY( d_data->spline.value( dtmp ) );
240    }
241
242    return fittedPoints;
243}
244
245class QwtWeedingCurveFitter::PrivateData
246{
247public:
248    PrivateData():
249        tolerance( 1.0 ),
250        chunkSize( 0 )
251    {
252    }
253
254    double tolerance;
255    uint chunkSize;
256};
257
258class QwtWeedingCurveFitter::Line
259{
260public:
261    Line( int i1 = 0, int i2 = 0 ):
262        from( i1 ),
263        to( i2 )
264    {
265    }
266
267    int from;
268    int to;
269};
270
271/*!
272   Constructor
273
274   \param tolerance Tolerance
275   \sa setTolerance(), tolerance()
276*/
277QwtWeedingCurveFitter::QwtWeedingCurveFitter( double tolerance )
278{
279    d_data = new PrivateData;
280    setTolerance( tolerance );
281}
282
283//! Destructor
284QwtWeedingCurveFitter::~QwtWeedingCurveFitter()
285{
286    delete d_data;
287}
288
289/*!
290 Assign the tolerance
291
292 The tolerance is the maximum distance, that is acceptable
293 between the original curve and the smoothed curve.
294
295 Increasing the tolerance will reduce the number of the
296 resulting points.
297
298 \param tolerance Tolerance
299
300 \sa tolerance()
301*/
302void QwtWeedingCurveFitter::setTolerance( double tolerance )
303{
304    d_data->tolerance = qMax( tolerance, 0.0 );
305}
306
307/*!
308  \return Tolerance
309  \sa setTolerance()
310*/
311double QwtWeedingCurveFitter::tolerance() const
312{
313    return d_data->tolerance;
314}
315
316/*!
317 Limit the number of points passed to a run of the algorithm
318
319 The runtime of the Douglas Peucker algorithm increases non linear
320 with the number of points. For a chunk size > 0 the polygon
321 is split into pieces passed to the algorithm one by one.
322
323 \param numPoints Maximum for the number of points passed to the algorithm
324
325 \sa chunkSize()
326*/
327void QwtWeedingCurveFitter::setChunkSize( uint numPoints )
328{
329    if ( numPoints > 0 )
330        numPoints = qMax( numPoints, 3U );
331
332    d_data->chunkSize = numPoints;
333}
334
335/*!
336 
337  \return Maximum for the number of points passed to a run
338          of the algorithm - or 0, when unlimited
339  \sa setChunkSize()
340*/
341uint QwtWeedingCurveFitter::chunkSize() const
342{
343    return d_data->chunkSize;
344}
345
346/*!
347  \param points Series of data points
348  \return Curve points
349*/
350QPolygonF QwtWeedingCurveFitter::fitCurve( const QPolygonF &points ) const
351{
352    QPolygonF fittedPoints;
353
354    if ( d_data->chunkSize == 0 )
355    {
356        fittedPoints = simplify( points );
357    }
358    else
359    {
360        for ( int i = 0; i < points.size(); i += d_data->chunkSize )
361        {
362            const QPolygonF p = points.mid( i, d_data->chunkSize );
363            fittedPoints += simplify( p );
364        }
365    }
366
367    return fittedPoints;
368}
369
370QPolygonF QwtWeedingCurveFitter::simplify( const QPolygonF &points ) const
371{
372    const double toleranceSqr = d_data->tolerance * d_data->tolerance;
373
374    QStack<Line> stack;
375    stack.reserve( 500 );
376
377    const QPointF *p = points.data();
378    const int nPoints = points.size();
379
380    QVector<bool> usePoint( nPoints, false );
381
382    stack.push( Line( 0, nPoints - 1 ) );
383
384    while ( !stack.isEmpty() )
385    {
386        const Line r = stack.pop();
387
388        // initialize line segment
389        const double vecX = p[r.to].x() - p[r.from].x();
390        const double vecY = p[r.to].y() - p[r.from].y();
391
392        const double vecLength = qSqrt( vecX * vecX + vecY * vecY );
393
394        const double unitVecX = ( vecLength != 0.0 ) ? vecX / vecLength : 0.0;
395        const double unitVecY = ( vecLength != 0.0 ) ? vecY / vecLength : 0.0;
396
397        double maxDistSqr = 0.0;
398        int nVertexIndexMaxDistance = r.from + 1;
399        for ( int i = r.from + 1; i < r.to; i++ )
400        {
401            //compare to anchor
402            const double fromVecX = p[i].x() - p[r.from].x();
403            const double fromVecY = p[i].y() - p[r.from].y();
404
405            double distToSegmentSqr;
406            if ( fromVecX * unitVecX + fromVecY * unitVecY < 0.0 )
407            {
408                distToSegmentSqr = fromVecX * fromVecX + fromVecY * fromVecY;
409            }
410            else
411            {
412                const double toVecX = p[i].x() - p[r.to].x();
413                const double toVecY = p[i].y() - p[r.to].y();
414                const double toVecLength = toVecX * toVecX + toVecY * toVecY;
415
416                const double s = toVecX * ( -unitVecX ) + toVecY * ( -unitVecY );
417                if ( s < 0.0 )
418                {
419                    distToSegmentSqr = toVecLength;
420                }
421                else
422                {
423                    distToSegmentSqr = qFabs( toVecLength - s * s );
424                }
425            }
426
427            if ( maxDistSqr < distToSegmentSqr )
428            {
429                maxDistSqr = distToSegmentSqr;
430                nVertexIndexMaxDistance = i;
431            }
432        }
433        if ( maxDistSqr <= toleranceSqr )
434        {
435            usePoint[r.from] = true;
436            usePoint[r.to] = true;
437        }
438        else
439        {
440            stack.push( Line( r.from, nVertexIndexMaxDistance ) );
441            stack.push( Line( nVertexIndexMaxDistance, r.to ) );
442        }
443    }
444
445    QPolygonF stripped;
446    for ( int i = 0; i < nPoints; i++ )
447    {
448        if ( usePoint[i] )
449            stripped += p[i];
450    }
451
452    return stripped;
453}
Note: See TracBrowser for help on using the repository browser.