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

Last change on this file since 4271 was 4271, checked in by mervart, 7 years ago
File size: 9.0 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    {
251    }
252
253    double tolerance;
254};
255
256class QwtWeedingCurveFitter::Line
257{
258public:
259    Line( int i1 = 0, int i2 = 0 ):
260        from( i1 ),
261        to( i2 )
262    {
263    }
264
265    int from;
266    int to;
267};
268
269/*!
270   Constructor
271
272   \param tolerance Tolerance
273   \sa setTolerance(), tolerance()
274*/
275QwtWeedingCurveFitter::QwtWeedingCurveFitter( double tolerance )
276{
277    d_data = new PrivateData;
278    setTolerance( tolerance );
279}
280
281//! Destructor
282QwtWeedingCurveFitter::~QwtWeedingCurveFitter()
283{
284    delete d_data;
285}
286
287/*!
288 Assign the tolerance
289
290 The tolerance is the maximum distance, that is accaptable
291 between the original curve and the smoothed curve.
292
293 Increasing the tolerance will reduce the number of the
294 resulting points.
295
296 \param tolerance Tolerance
297
298 \sa tolerance()
299*/
300void QwtWeedingCurveFitter::setTolerance( double tolerance )
301{
302    d_data->tolerance = qMax( tolerance, 0.0 );
303}
304
305/*!
306  \return Tolerance
307  \sa setTolerance()
308*/
309double QwtWeedingCurveFitter::tolerance() const
310{
311    return d_data->tolerance;
312}
313
314/*!
315  \param points Series of data points
316  \return Curve points
317*/
318QPolygonF QwtWeedingCurveFitter::fitCurve( const QPolygonF &points ) const
319{
320    QStack<Line> stack;
321    stack.reserve( 500 );
322
323    const QPointF *p = points.data();
324    const int nPoints = points.size();
325
326    QVector<bool> usePoint( nPoints, false );
327
328    double distToSegment;
329
330    stack.push( Line( 0, nPoints - 1 ) );
331
332    while ( !stack.isEmpty() )
333    {
334        const Line r = stack.pop();
335
336        // initialize line segment
337        const double vecX = p[r.to].x() - p[r.from].x();
338        const double vecY = p[r.to].y() - p[r.from].y();
339
340        const double vecLength = qSqrt( vecX * vecX + vecY * vecY );
341
342        const double unitVecX = ( vecLength != 0.0 ) ? vecX / vecLength : 0.0;
343        const double unitVecY = ( vecLength != 0.0 ) ? vecY / vecLength : 0.0;
344
345        double maxDist = 0.0;
346        int nVertexIndexMaxDistance = r.from + 1;
347        for ( int i = r.from + 1; i < r.to; i++ )
348        {
349            //compare to anchor
350            const double fromVecX = p[i].x() - p[r.from].x();
351            const double fromVecY = p[i].y() - p[r.from].y();
352            const double fromVecLength =
353                qSqrt( fromVecX * fromVecX + fromVecY * fromVecY );
354
355            if ( fromVecX * unitVecX + fromVecY * unitVecY < 0.0 )
356            {
357                distToSegment = fromVecLength;
358            }
359            if ( fromVecX * unitVecX + fromVecY * unitVecY < 0.0 )
360            {
361                distToSegment = fromVecLength;
362            }
363            else
364            {
365                const double toVecX = p[i].x() - p[r.to].x();
366                const double toVecY = p[i].y() - p[r.to].y();
367                const double toVecLength = qSqrt( toVecX * toVecX + toVecY * toVecY );
368                const double s = toVecX * ( -unitVecX ) + toVecY * ( -unitVecY );
369                if ( s < 0.0 )
370                    distToSegment = toVecLength;
371                else
372                {
373                    distToSegment = qSqrt( qFabs( toVecLength * toVecLength - s * s ) );
374                }
375            }
376
377            if ( maxDist < distToSegment )
378            {
379                maxDist = distToSegment;
380                nVertexIndexMaxDistance = i;
381            }
382        }
383        if ( maxDist <= d_data->tolerance )
384        {
385            usePoint[r.from] = true;
386            usePoint[r.to] = true;
387        }
388        else
389        {
390            stack.push( Line( r.from, nVertexIndexMaxDistance ) );
391            stack.push( Line( nVertexIndexMaxDistance, r.to ) );
392        }
393    }
394
395    int cnt = 0;
396
397    QPolygonF stripped( nPoints );
398    for ( int i = 0; i < nPoints; i++ )
399    {
400        if ( usePoint[i] )
401            stripped[cnt++] = p[i];
402    }
403    stripped.resize( cnt );
404    return stripped;
405}
Note: See TracBrowser for help on using the repository browser.