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
|
---|
21 | QwtCurveFitter::QwtCurveFitter()
|
---|
22 | {
|
---|
23 | }
|
---|
24 |
|
---|
25 | //! Destructor
|
---|
26 | QwtCurveFitter::~QwtCurveFitter()
|
---|
27 | {
|
---|
28 | }
|
---|
29 |
|
---|
30 | class QwtSplineCurveFitter::PrivateData
|
---|
31 | {
|
---|
32 | public:
|
---|
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
|
---|
45 | QwtSplineCurveFitter::QwtSplineCurveFitter()
|
---|
46 | {
|
---|
47 | d_data = new PrivateData;
|
---|
48 | }
|
---|
49 |
|
---|
50 | //! Destructor
|
---|
51 | QwtSplineCurveFitter::~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 | */
|
---|
62 | void QwtSplineCurveFitter::setFitMode( FitMode mode )
|
---|
63 | {
|
---|
64 | d_data->fitMode = mode;
|
---|
65 | }
|
---|
66 |
|
---|
67 | /*!
|
---|
68 | \return Mode representing a spline algorithm
|
---|
69 | \sa setFitMode()
|
---|
70 | */
|
---|
71 | QwtSplineCurveFitter::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 | */
|
---|
82 | void 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 | */
|
---|
92 | const QwtSpline &QwtSplineCurveFitter::spline() const
|
---|
93 | {
|
---|
94 | return d_data->spline;
|
---|
95 | }
|
---|
96 |
|
---|
97 | /*!
|
---|
98 | \return Spline
|
---|
99 | \sa setSpline()
|
---|
100 | */
|
---|
101 | QwtSpline &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 | */
|
---|
112 | void QwtSplineCurveFitter::setSplineSize( int splineSize )
|
---|
113 | {
|
---|
114 | d_data->splineSize = qMax( splineSize, 10 );
|
---|
115 | }
|
---|
116 |
|
---|
117 | /*!
|
---|
118 | \return Spline size
|
---|
119 | \sa setSplineSize()
|
---|
120 | */
|
---|
121 | int 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 | */
|
---|
132 | QPolygonF 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 |
|
---|
160 | QPolygonF 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 |
|
---|
188 | QPolygonF 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 |
|
---|
245 | class QwtWeedingCurveFitter::PrivateData
|
---|
246 | {
|
---|
247 | public:
|
---|
248 | PrivateData():
|
---|
249 | tolerance( 1.0 ),
|
---|
250 | chunkSize( 0 )
|
---|
251 | {
|
---|
252 | }
|
---|
253 |
|
---|
254 | double tolerance;
|
---|
255 | uint chunkSize;
|
---|
256 | };
|
---|
257 |
|
---|
258 | class QwtWeedingCurveFitter::Line
|
---|
259 | {
|
---|
260 | public:
|
---|
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 | */
|
---|
277 | QwtWeedingCurveFitter::QwtWeedingCurveFitter( double tolerance )
|
---|
278 | {
|
---|
279 | d_data = new PrivateData;
|
---|
280 | setTolerance( tolerance );
|
---|
281 | }
|
---|
282 |
|
---|
283 | //! Destructor
|
---|
284 | QwtWeedingCurveFitter::~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 | */
|
---|
302 | void QwtWeedingCurveFitter::setTolerance( double tolerance )
|
---|
303 | {
|
---|
304 | d_data->tolerance = qMax( tolerance, 0.0 );
|
---|
305 | }
|
---|
306 |
|
---|
307 | /*!
|
---|
308 | \return Tolerance
|
---|
309 | \sa setTolerance()
|
---|
310 | */
|
---|
311 | double 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 | */
|
---|
327 | void QwtWeedingCurveFitter::setChunkSize( uint numPoints )
|
---|
328 | {
|
---|
329 | if ( numPoints > 0 )
|
---|
330 | numPoints = qMax( numPoints, 3U );
|
---|
331 |
|
---|
332 | d_data->chunkSize = numPoints;
|
---|
333 | }
|
---|
334 |
|
---|
335 | /*!
|
---|
336 | \return Maximum for the number of points passed to a run
|
---|
337 | of the algorithm - or 0, when unlimited
|
---|
338 | \sa setChunkSize()
|
---|
339 | */
|
---|
340 | uint QwtWeedingCurveFitter::chunkSize() const
|
---|
341 | {
|
---|
342 | return d_data->chunkSize;
|
---|
343 | }
|
---|
344 |
|
---|
345 | /*!
|
---|
346 | \param points Series of data points
|
---|
347 | \return Curve points
|
---|
348 | */
|
---|
349 | QPolygonF QwtWeedingCurveFitter::fitCurve( const QPolygonF &points ) const
|
---|
350 | {
|
---|
351 | if ( points.isEmpty() )
|
---|
352 | return points;
|
---|
353 |
|
---|
354 | QPolygonF fittedPoints;
|
---|
355 | if ( d_data->chunkSize == 0 )
|
---|
356 | {
|
---|
357 | fittedPoints = simplify( points );
|
---|
358 | }
|
---|
359 | else
|
---|
360 | {
|
---|
361 | for ( int i = 0; i < points.size(); i += d_data->chunkSize )
|
---|
362 | {
|
---|
363 | const QPolygonF p = points.mid( i, d_data->chunkSize );
|
---|
364 | fittedPoints += simplify( p );
|
---|
365 | }
|
---|
366 | }
|
---|
367 |
|
---|
368 | return fittedPoints;
|
---|
369 | }
|
---|
370 |
|
---|
371 | QPolygonF QwtWeedingCurveFitter::simplify( const QPolygonF &points ) const
|
---|
372 | {
|
---|
373 | const double toleranceSqr = d_data->tolerance * d_data->tolerance;
|
---|
374 |
|
---|
375 | QStack<Line> stack;
|
---|
376 | stack.reserve( 500 );
|
---|
377 |
|
---|
378 | const QPointF *p = points.data();
|
---|
379 | const int nPoints = points.size();
|
---|
380 |
|
---|
381 | QVector<bool> usePoint( nPoints, false );
|
---|
382 |
|
---|
383 | stack.push( Line( 0, nPoints - 1 ) );
|
---|
384 |
|
---|
385 | while ( !stack.isEmpty() )
|
---|
386 | {
|
---|
387 | const Line r = stack.pop();
|
---|
388 |
|
---|
389 | // initialize line segment
|
---|
390 | const double vecX = p[r.to].x() - p[r.from].x();
|
---|
391 | const double vecY = p[r.to].y() - p[r.from].y();
|
---|
392 |
|
---|
393 | const double vecLength = qSqrt( vecX * vecX + vecY * vecY );
|
---|
394 |
|
---|
395 | const double unitVecX = ( vecLength != 0.0 ) ? vecX / vecLength : 0.0;
|
---|
396 | const double unitVecY = ( vecLength != 0.0 ) ? vecY / vecLength : 0.0;
|
---|
397 |
|
---|
398 | double maxDistSqr = 0.0;
|
---|
399 | int nVertexIndexMaxDistance = r.from + 1;
|
---|
400 | for ( int i = r.from + 1; i < r.to; i++ )
|
---|
401 | {
|
---|
402 | //compare to anchor
|
---|
403 | const double fromVecX = p[i].x() - p[r.from].x();
|
---|
404 | const double fromVecY = p[i].y() - p[r.from].y();
|
---|
405 |
|
---|
406 | double distToSegmentSqr;
|
---|
407 | if ( fromVecX * unitVecX + fromVecY * unitVecY < 0.0 )
|
---|
408 | {
|
---|
409 | distToSegmentSqr = fromVecX * fromVecX + fromVecY * fromVecY;
|
---|
410 | }
|
---|
411 | else
|
---|
412 | {
|
---|
413 | const double toVecX = p[i].x() - p[r.to].x();
|
---|
414 | const double toVecY = p[i].y() - p[r.to].y();
|
---|
415 | const double toVecLength = toVecX * toVecX + toVecY * toVecY;
|
---|
416 |
|
---|
417 | const double s = toVecX * ( -unitVecX ) + toVecY * ( -unitVecY );
|
---|
418 | if ( s < 0.0 )
|
---|
419 | {
|
---|
420 | distToSegmentSqr = toVecLength;
|
---|
421 | }
|
---|
422 | else
|
---|
423 | {
|
---|
424 | distToSegmentSqr = qFabs( toVecLength - s * s );
|
---|
425 | }
|
---|
426 | }
|
---|
427 |
|
---|
428 | if ( maxDistSqr < distToSegmentSqr )
|
---|
429 | {
|
---|
430 | maxDistSqr = distToSegmentSqr;
|
---|
431 | nVertexIndexMaxDistance = i;
|
---|
432 | }
|
---|
433 | }
|
---|
434 | if ( maxDistSqr <= toleranceSqr )
|
---|
435 | {
|
---|
436 | usePoint[r.from] = true;
|
---|
437 | usePoint[r.to] = true;
|
---|
438 | }
|
---|
439 | else
|
---|
440 | {
|
---|
441 | stack.push( Line( r.from, nVertexIndexMaxDistance ) );
|
---|
442 | stack.push( Line( nVertexIndexMaxDistance, r.to ) );
|
---|
443 | }
|
---|
444 | }
|
---|
445 |
|
---|
446 | QPolygonF stripped;
|
---|
447 | for ( int i = 0; i < nPoints; i++ )
|
---|
448 | {
|
---|
449 | if ( usePoint[i] )
|
---|
450 | stripped += p[i];
|
---|
451 | }
|
---|
452 |
|
---|
453 | return stripped;
|
---|
454 | }
|
---|