1 /**********************************************************
2 * Version $Id$
3 *********************************************************/
4 /* Plot Graphic Library
5
6 Copyright (C) 2002 Pelikhan, Inc. All rights reserved
7 Go to http://eauminerale.syldavie.csam.ucl.ac.be/peli/pgl/pgl.html for the latest PGL binaries
8 This software is NOT freeware.
9
10 This software is provided "as is", with no warranty.
11 Read the license agreement provided with the files
12
13 This software, and all accompanying files, data and materials, are distributed "AS IS" and with no warranties of any kind,
14 whether express or implied. his disclaimer of warranty constitutes an essential part of the agreement.
15 In no event shall Pelikhan, or its principals, shareholders, officers, employees, affiliates, contractors, subsidiaries,
16 or parent organizations, be liable for any incidental, consequential, or punitive damages whatsoever relating to the use of PGL,
17 or your relationship with Pelikhan.
18
19 Author: Jonathan de Halleux, dehalleux@auto.ucl.ac.be
20 */
21 #if !defined(AFX_LINEAPPROXIMATOR_H__F5E6E8DC_1185_4AC0_A061_7B3309700E9D__INCLUDED_)
22 #define AFX_LINEAPPROXIMATOR_H__F5E6E8DC_1185_4AC0_A061_7B3309700E9D__INCLUDED_
23
24 #if _MSC_VER > 1000
25 #pragma once
26 #endif // _MSC_VER > 1000
27
28 #include <crtdbg.h>
29 #ifndef ASSERT
30 #ifdef _DEBUG
31 #define ASSERT(a) _ASSERT(a);
32 #else
33 #define ASSERT(a)
34 #endif
35 #endif
36 #ifdef min
37 #undef min
38 #endif
39 #ifdef max
40 #undef max
41 #endif
42
43 #include <iostream>
44 #include <vector>
45 #include <list>
46 #include <limits>
47 #include <algorithm>
48
49 /*!
50 \defgroup LAGroup Line approximation algorithms
51 */
52 namespace hull
53 {
54
55 /*! \brief A point (x,y)
56
57 \param T float or double
58
59 A pair of (x,y) values.
60 \ingroup LAGroup
61 */
62 template<typename T>
63 class TPoint
64 {
65 public:
66 //! Default constructor
x(_x)67 TPoint( T _x = 0, T _y=0):x(_x), y(_y){};
68 //! Assignement constructor
69 TPoint& operator = ( const TPoint<T>& t) { if (&t != this){ x=t.x; y=t.y;} return *this;};
70
71 //! x value
72 T x;
73 //! y value
74 T y;
75 };
76
77 /*! \brief A point (x,y) with references
78
79 \param T float or double
80
81 \ingroup LAGroup
82 */
83 template<typename T>
84 class TPointRef
85 {
86 public:
87 //! Default contructor
TPointRef()88 TPointRef():x(xDummy),y(yDummy){};
89 /*! \brief Constructor with 2 values
90
91 \param _x value to get reference to.
92 \param _y value to get reference to.
93 */
TPointRef(T & _x,T & _y)94 TPointRef(T& _x, T& _y):x(_x),y(_y){};
95 //! Assignement constructor
96 TPointRef<T>& operator = (const TPointRef<T>& t) { if (this != &t){ x=t.x;y=t.y;} return *this;};
97 //! Assignement constructor with point
98 TPointRef<T>& operator = (TPoint<T> t) { x=t.x;y=t.y; return *this;};
99 //! x, reference to a value
100 T& x;
101 //! y, reference to a value
102 T& y;
103 private:
104 static T xDummy;
105 static T yDummy;
106 };
107
108 template<typename T> T TPointRef<T>::xDummy=0;
109 template<typename T> T TPointRef<T>::yDummy=0;
110
111 /*! \brief Virtual base class for Line approximator classes
112
113 \par Template arguments
114
115 \param T float or double
116 \param TPointContainer a container of points (structure with x,y) with random access iterators
117 \param TKeyContainer a container of TPointContainer::const_iterator (structure with x,y) with single direction iterators
118
119
120 \par Notations:
121
122 - points : data to interpolate
123 - keys : selected points from available set that interpolate within the desired tolerance
124
125 \par Containers:
126
127 - Points are stored in a #PointContainer. #PointContainer is a container such that
128 #- has random access iterator,
129 #- value_type is a structure/class with x,y members
130 - Keys are stored in a #KeyContainer. #KeyContainer is a container such that:
131 #- has single directional iterator,
132 #- value_type is PointContainer::const_iterator
133
134 \par Data normalization:
135
136 To avoid numerical instability when computing cross product and so, data is normalized ( see #NormalizePoints ).
137 To enable/disable normalization, use #SetNormalization.
138
139 \ingroup LAGroup
140 */
141 template<typename T, typename TPointContainer, typename TKeyContainer>
142 class TLine
143 {
144 public:
145
146 //! \name Structures and typedefs
147 //@{
148 //! Point container
149 typedef TPointContainer PointContainer;
150 //! Key container
151 typedef TKeyContainer KeyContainer;
152 //! 2D homogenous point
153 struct SHomog
154 {
155 public:
156 //! Default constructor
157 SHomog(T _x=0, T _y=0, T _w=1) { x=_x; y=_y; w=_w;};
158
159 T x;
160 T y;
161 T w;
162 };
163
164 /*! \brief returns square of euclidian distance
165
166 */
SqrDist(const TPointContainer::const_iterator & p,const TPointContainer::const_iterator & q)167 static T SqrDist( const TPointContainer::const_iterator& p, const TPointContainer::const_iterator& q)
168 {
169 T dx=p->x-q->x;
170 T dy=p->y-q->y;
171 return dx*dx+dy*dy;
172 }
173
174 /*! \brief Cross product
175
176 \param p an iterator with x,y members
177 \param q an iterator with x,y members
178 \result r cross product of p,q
179
180 The euclidian cross product of p,q:
181 \f[ \vec r = \vec p \times \vec q = \left( \begin{array}{c} p_x q_y - p_y q_x \\ -q_y + p_y \\ q_x - p_x \end{array}\right)\f]
182 */
CrossProduct(const TPointContainer::const_iterator & p,const TPointContainer::const_iterator & q,SHomog & r)183 static void CrossProduct( const TPointContainer::const_iterator& p, const TPointContainer::const_iterator& q, SHomog& r)
184 {
185 r.w = p->x * q->y - p->y * q->x;
186 r.x = - q->y + p->y;
187 r.y = q->x - p->x;
188 };
189
190 /*! \brief Dot product
191
192 \param p an iterator with x,y members
193 \param q an 2D homogenous point
194 \result dot product of p,q
195
196 The euclidian dot product of p,q:
197 \f[ <\vec p, \vec q> = q_w + p_x q_x + p_y q_y \f]
198 */
DotProduct(const TPointContainer::const_iterator & p,const SHomog & q)199 static T DotProduct( const TPointContainer::const_iterator& p, const SHomog& q)
200 {
201 return q.w + p->x*q.x + p->y*q.y;
202 };
203
204 /*! \brief Dot product
205
206 \param p an iterator with x,y members
207 \param q an iterator with x,y members
208 \result dot product of p,q
209
210 The euclidian dot product of p,q:
211 \f[ < \vec p,\vec q> = p_x q_x + p_y q_y \f]
212 */
DotProduct(const TPointContainer::const_iterator & p,const TPointContainer::const_iterator & q)213 static T DotProduct( const TPointContainer::const_iterator& p, const TPointContainer::const_iterator& q)
214 {
215 return p->x*q->x + p->y*q->y;
216 };
217
218 /*! \brief Linear combination of points
219
220 \param a p multiplier
221 \param p a point
222 \param b q multiplier
223 \param q a point
224 \param r linear combination of p,q
225
226 The linear combination is:
227 \f[ \vec r = a \vec p + b \vec q \f]
228 */
LinComb(T a,const TPointContainer::const_iterator & p,T b,const TPointContainer::const_iterator & q,const TPointContainer::const_iterator & r)229 static void LinComb( T a, const TPointContainer::const_iterator& p, T b, const TPointContainer::const_iterator& q, const TPointContainer::const_iterator& r)
230 {
231 r->x = a * p->x + b * q->x;
232 r->y = a * p->y + b * q->y;
233 };
234
235 //! Internal limit structure
236 struct SLimits
237 {
238 T dMinX;
239 T dMaxX;
240 T dMinY;
241 T dMaxY;
GetCenterXSLimits242 T GetCenterX() { return static_cast<T>((dMaxX+dMinX)/2.0);};
GetCenterYSLimits243 T GetCenterY() { return static_cast<T>((dMaxY+dMinY)/2.0);};
GetWidthSLimits244 T GetWidth() { return static_cast<T>(dMaxX-dMinX);};
GetHeightSLimits245 T GetHeight() { return static_cast<T>(dMaxY-dMinY);};
246 };
247
248 //! T container
249 typedef std::vector<T> TVector;
250 //@}
251
252 //! Default constructor
TLine()253 TLine():m_bNormalization(true){};
254
255 //! \name Point handling
256 //@{
257 //! returns number of points
GetPointSize()258 size_t GetPointSize() const { return m_cPoints.size();};
259 //! return vector of points
GetPoints()260 PointContainer& GetPoints() { return m_cPoints;};
261 //! return vector of points, const
GetPoints()262 const PointContainer& GetPoints() const { return m_cPoints;};
263 //@}
264
265 //! \name Key handling
266 //@{
267 //! returns number of keys
GetKeySize()268 size_t GetKeySize() const { return m_cKeys.size();};
269 //! return keys
GetKeys()270 KeyContainer& GetKeys() { return m_cKeys;};
271 //! return keys, const
GetKeys()272 const KeyContainer& GetKeys() const { return m_cKeys;};
273 //@}
274
275 //! \name Helpers
276 //@{
277 //! Setting points from vectors
278 void SetPoints( const std::vector<T>& vX, const std::vector<T>& vY);
279 //! Returning keys to vectors
280 void GetKeys( std::vector<T>& vX, std::vector<T>& vY) const;
281 //@}
282
283 //! \name Bounding box
284 //@{
285 //! compute the bounding box
286 void ComputeBoundingBox();
287 //! return the point bounding box
GetBoundingBox()288 const SLimits& GetBoundingBox() const { return m_limits;};
289 //@}
290
291 //! \name Normalization
292 //@{
293 /*! \brief Point normalization
294
295 Let \f$(x_i,y_i)\f$, the original points and \f$(\hat x_x, \hat y_i)\f$ the normalized points:
296 \f[
297 \hat x_i = \frac{x_i - \bar x]}{\max_i (x_i-x_j)}
298 \f]
299 where \f$\bar x\f$, \f$\bar y\f$ denote respectively the mean value of the \f$x_i\f$ and \f$y_i\f$.
300
301 \sa DeNormalizePoints
302 */
303 void NormalizePoints();
304
305 /*! \brief Roll back normalization
306
307 \sa NormalizePoints
308 */
309 void DeNormalizePoints();
310 //! enabled, disable normalization
311 void SetNormalization( bool bEnabled = true) { m_bNormalization = true;};
312 //! returns true if normalizing
IsNormalization()313 bool IsNormalization() const { return m_bNormalization;};
314 //@}
315
316 //! \name Double points checking and loop...
317 //@{
318 /*! \brief Discard double points
319
320 For each pair of points \f$p_i, p_{i+1}\f$, discards \f$p_{i+1}\f$ if
321 \f[ \| p_i - p_{i+1} \| < \varepsilon \f]
322 */
323 void DiscardDoublePoints();
324 //! Test for loops
325 void FindLoop(size_t uStartPoint, size_t& uEndPoint);
326 //@}
327
328 protected:
329 //! \name Attributes
330 //@{
331 TPointContainer m_cPoints;
332 TKeyContainer m_cKeys;
333 SLimits m_limits;
334 bool m_bNormalization;
335 //@}
336 };
337
338 namespace priv
339 {
340 template<class Container, class Pred>
341 struct PredX : public std::binary_function< Container::iterator, Container::iterator, bool>
342 {
operatorPredX343 bool operator()( const Container::value_type& p1, const Container::value_type& p2)
344 { return m_pred(p1.x, p2.x); };
345 protected:
346 Pred m_pred;
347 };
348
349 template<class Container, class Pred>
350 struct PredY : public std::binary_function< Container::iterator, Container::iterator, bool>
351 {
operatorPredY352 bool operator()( const Container::value_type& p1, const Container::value_type& p2)
353 { return m_pred(p1.y, p2.y); };
354 protected:
355 Pred m_pred;
356 };
357 };
358
359 template <typename T, typename TPointContainer, typename TKeyContainer>
ComputeBoundingBox()360 void TLine<T, TPointContainer,TKeyContainer>::ComputeBoundingBox()
361 {
362 if (m_cPoints.size() < 2)
363 return;
364
365 PointContainer::const_iterator it = (*((const TPointContainer*)&m_cPoints)).begin();
366
367
368 //finding minimum and maximum...
369 m_limits.dMinX=std::min_element( m_cPoints.begin(), m_cPoints.end(), priv::PredX<TPointContainer, std::less<T> >() )->x ;
370 m_limits.dMaxX=std::max_element( m_cPoints.begin(), m_cPoints.end(), priv::PredX<TPointContainer, std::less<T> >() )->x ;
371 m_limits.dMinY=std::min_element( m_cPoints.begin(), m_cPoints.end(), priv::PredY<TPointContainer, std::less<T> >() )->y ;
372 m_limits.dMaxY=std::max_element( m_cPoints.begin(), m_cPoints.end(), priv::PredY<TPointContainer, std::less<T> >() )->y ;
373
374 if ( fabs( m_limits.GetWidth() ) < std::numeric_limits<T>::epsilon() )
375 {
376 m_limits.dMaxX = m_limits.dMinX+1;
377 }
378 if ( fabs( m_limits.GetHeight() ) < std::numeric_limits<T>::epsilon() )
379 {
380 m_limits.dMaxY = m_limits.dMinY+1;
381 }
382 }
383
384 namespace priv
385 {
386 template<typename T>
387 struct Rect
388 {
RectRect389 Rect( T xm, T ym, T dx, T dy)
390 :m_xm(xm),m_ym(ym),m_dx(dx),m_dy(dy){};
391 protected:
392 T m_xm;
393 T m_ym;
394 T m_dx;
395 T m_dy;
396 };
397
398 template<typename T>
399 struct NormalizePoint : public std::unary_function< TPoint<T>& , int>, public Rect<T>
400 {
NormalizePointNormalizePoint401 NormalizePoint( T xm, T ym, T dx, T dy)
402 : Rect<T>(xm,ym,dx,dy){};
operatorNormalizePoint403 int operator() ( TPoint<T>& point)
404 {
405 point.x=(point.x-m_xm)/m_dx;
406 point.y=(point.y-m_ym)/m_dy;
407 return 0;
408 };
409 };
410
411 template<typename T>
412 struct DeNormalizePoint : public std::unary_function< TPoint<T>& , int>, public Rect<T>
413 {
DeNormalizePointDeNormalizePoint414 DeNormalizePoint( T xm, T ym, T dx, T dy)
415 : Rect<T>(xm,ym,dx,dy){};
operatorDeNormalizePoint416 int operator() ( TPoint<T>& point)
417 {
418 point.x=m_xm+point.x*m_dx;
419 point.y=m_ym+point.y*m_dy;
420 return 0;
421 };
422 };
423 };
424
425 template <typename T, typename TPointContainer, typename TKeyContainer>
NormalizePoints()426 void TLine<T, TPointContainer, TKeyContainer>::NormalizePoints()
427 {
428 T xm,ym,dx,dy;
429 // normalizing...
430 xm=m_limits.GetCenterX();
431 ym=m_limits.GetCenterY();
432 dx=m_limits.GetWidth();
433 dy=m_limits.GetHeight();
434
435 std::for_each(
436 m_cPoints.begin(),
437 m_cPoints.end(),
438 priv::NormalizePoint<T>( xm, ym, dx, dy)
439 );
440 }
441
442 template <typename T, typename TPointContainer, typename TKeyContainer>
DeNormalizePoints()443 void TLine<T, TPointContainer, TKeyContainer>::DeNormalizePoints()
444 {
445 T xm,ym,dx,dy;
446 // normalizing...
447 xm=m_limits.GetCenterX();
448 ym=m_limits.GetCenterY();
449 dx=m_limits.GetWidth();
450 dy=m_limits.GetHeight();
451
452 std::for_each(
453 m_cPoints.begin(),
454 m_cPoints.end(),
455 priv::DeNormalizePoint<T>( xm, ym, dx, dy)
456 );
457 }
458
459 template <typename T, typename TPointContainer, typename TKeyContainer>
DiscardDoublePoints()460 void TLine<T, TPointContainer, TKeyContainer>::DiscardDoublePoints()
461 {
462 // creating a list...
463 TPointContainer& pc=GetPoints();
464 TPointContainer::iterator it, it1;
465 T epsilon2=std::numeric_limits<T>::epsilon();
466
467 #ifdef _DEBUG
468 size_t count=0;
469 #endif
470
471 it1=it=pc.begin();
472 ++it1;
473 while (it !=pc.end() && it1!=pc.end())
474 {
475 if ( SqrDist(it, it1) < epsilon2 )
476 {
477 it1=pc.erase(it1);
478 }
479 else
480 {
481 ++it; ++it1;
482 }
483 }
484
485 #ifdef _DEBUG
486 TRACE( _T("Numer of (double) points erased: %d\n"), count);
487 #endif
488 };
489
490 template <typename T, typename TPointContainer, typename TKeyContainer>
SetPoints(const std::vector<T> & vX,const std::vector<T> & vY)491 void TLine<T, TPointContainer, TKeyContainer>::SetPoints( const std::vector<T>& vX, const std::vector<T>& vY)
492 {
493 TPointContainer& pc=GetPoints();
494 const size_t n = __min( vX.size(), vY.size());
495
496 pc.resize(n);
497 for (size_t i=0;i<n;i++)
498 {
499 pc[i]= TPoint<T>( vX[i], vY[i]);
500 }
501 };
502
503 template <typename T, typename TPointContainer, typename TKeyContainer>
GetKeys(std::vector<T> & vX,std::vector<T> & vY)504 void TLine<T, TPointContainer, TKeyContainer>::GetKeys( std::vector<T>& vX, std::vector<T>& vY) const
505 {
506 const TKeyContainer& kc=GetKeys();
507 TKeyContainer::const_iterator it;
508 size_t i;
509
510 vX.resize(kc.size());
511 vY.resize(kc.size());
512
513 for (it=kc.begin(), i=0;it!=kc.end();it++, i++)
514 {
515 vX[i]=(*it)->x;
516 vY[i]=(*it)->y;
517 }
518 };
519
520 template <typename T, typename TPointContainer, typename TKeyContainer>
FindLoop(size_t uStartPoint,size_t & uEndPoint)521 void TLine<T, TPointContainer, TKeyContainer>::FindLoop(size_t uStartPoint, size_t& uEndPoint)
522 {
523 };
524
525
526 /*! \brief Base class for line approximators
527
528
529
530 \ingroup LAGroup
531 */
532 template <typename T, typename TPointContainer, typename TKeyContainer>
533 class TLineApproximator : virtual public TLine<T,TPointContainer, TKeyContainer>
534 {
535 public:
536
537 //! \name Constructor
538 //@{
TLineApproximator()539 TLineApproximator(): m_dTol(0)
540 {m_limits.dMinX=m_limits.dMinY=0;m_limits.dMaxX=m_limits.dMaxY=1;};
~TLineApproximator()541 ~TLineApproximator(){};
542 //@}
543
544
545 //! \name Tolerance
546 //@{
547 //! sets the tolerance
SetTol(double dTol)548 void SetTol( double dTol) { m_dTol = __max( dTol, 0);};
549 //! return current tolerance
GetTol()550 double GetTol() const { return m_dTol;};
551 //@}
552
553 //! \name Simplification functions
554 //@{
555 //! Initialize simplification
ClearKeys()556 void ClearKeys() { m_cKeys.clear();};
557 //! Compute the keys
558 void Simplify();
559 /*! Shrink to compression level
560
561 \param dScale scaling to apply [0...1]
562 \param dScaleTol [optional] tolerance with respect to dScale, default is 0.05
563 \param eTolRight [optional] first estimate on right tolerance
564 \param nMaxIter [optional] maximum number of iterations, default is 250
565 \return number of estimations
566 */
567 size_t ShrinkNorm( T dScale, T dScaleTol = 0.05, T eTolRight = -1, size_t nMaxIter = 250);
568
569 /*! Shrink to a specified number of points
570
571 \param n desired number of points in the approximate curve
572 \param nTol [optional] tolerance with respect to n, default is 10
573 \param eTolRight [optional] first estimate on right tolerance
574 \param nMaxIter [optional] maximum number of iterations, default is 250
575 \return number of estimations
576 */
577 size_t Shrink( size_t nDesiredPoints, size_t nTol = 10, T eTolRight = -1, size_t nMaxIter = 250);
578 //@}
579
580 protected:
581 //! \name Virtual functions
582 //@{
583 /*! \brief Virtual approximation function
584
585 This function must be overriden in inherited classes. To implement your own algorithm,
586 override this function.
587 */
ComputeKeys()588 virtual void ComputeKeys() { ClearKeys();};
589 //@}
590
591 private:
592 T m_dTol;
593 };
594
595
596 template <typename T, typename TPointContainer, typename TKeyContainer>
ShrinkNorm(T dScale,T dScaleTol,T eTolRight,size_t nMaxIter)597 size_t TLineApproximator<T,TPointContainer,TKeyContainer>::ShrinkNorm( T dScale, T dScaleTol, T eTolRight ,size_t nMaxIter)
598 {
599 // number of points wanted...
600 size_t uWantedPoints= __min(m_cPoints.size(), __max(2, static_cast<size_t>(floor(m_cPoints.size()*dScale))));
601 size_t uTol = __min(m_cPoints.size(), __max(0, static_cast<size_t>(floor(m_cPoints.size()*dScaleTol)) ));
602
603 return Shrink( uWantedPoints, uTol, eTolRight, nMaxIter);
604 }
605
606 namespace priv
607 {
608 // (C) Copyright Gennadiy Rozental 2001-2002.
609 // Permission to copy, use, modify, sell and distribute this software
610 // is granted provided this copyright notice appears in all copies.
611 // This software is provided "as is" without express or implied warranty,
612 // and with no claim as to its suitability for any purpose.
613
614 // See http://www.boost.org for most recent version including documentation.
615 //
616 // File : $RCSfile: LineApproximator.h,v $
617 //
618 // Version : $Id$
619 //
620 // Description : defines algoirthms for comparing 2 floating point values
621 // ***************************************************************************
622 template<typename FPT>
623 inline FPT
fpt_abs(FPT arg)624 fpt_abs( FPT arg )
625 {
626 return arg < 0 ? -arg : arg;
627 }
628
629 // both f1 and f2 are unsigned here
630 template<typename FPT>
631 inline FPT
safe_fpt_division(FPT uf1,FPT uf2)632 safe_fpt_division( FPT uf1, FPT uf2 )
633 {
634 return ( uf1 < 1 && uf1 > uf2 * std::numeric_limits<FPT>::max())
635 ? std::numeric_limits<FPT>::max() :
636 ((uf2 > 1 && uf1 < uf2 * std::numeric_limits<FPT>::min() ||
637 uf1 == 0) ? 0 :
638 uf1/uf2 );
639 }
640
641 template<typename FPT>
642 class close_at_tolerance
643 {
644 public:
645 explicit close_at_tolerance( FPT tolerance, bool strong_or_weak = true )
p_tolerance(tolerance)646 : p_tolerance( tolerance ),m_strong_or_weak( strong_or_weak ) { };
647
648 explicit close_at_tolerance( int number_of_rounding_errors, bool strong_or_weak = true )
649 : p_tolerance( std::numeric_limits<FPT>::epsilon() * number_of_rounding_errors/2 ),
650 m_strong_or_weak( strong_or_weak ) {}
651
operator()652 bool operator()( FPT left, FPT right ) const
653 {
654 FPT diff = fpt_abs( left - right );
655 FPT d1 = safe_fpt_division( diff, fpt_abs( right ) );
656 FPT d2 = safe_fpt_division( diff, fpt_abs( left ) );
657
658 return m_strong_or_weak ? (d1 <= p_tolerance.get() && d2 <= p_tolerance.get())
659 : (d1 <= p_tolerance.get() || d2 <= p_tolerance.get());
660 }
661
662 // Data members
663 class p_tolerance_class
664 {
665 private:
666 FPT f;
667 public:
f(_f)668 p_tolerance_class(FPT _f=0):f(_f){};
get()669 FPT get() const{ return f;};
670 };
671 p_tolerance_class p_tolerance;
672 private:
673 bool m_strong_or_weak;
674 };
675
676 template <typename T>
IsEqual(T x,T y)677 inline bool IsEqual(T x, T y)
678 {
679 static close_at_tolerance<T> comp( std::numeric_limits<T>::epsilon()/2*10);
680 return comp(fpt_abs(x),fpt_abs(y));
681 };
682
683 template <typename T>
IsEmptyInterval(T x,T y)684 inline bool IsEmptyInterval(T x, T y)
685 {
686 return ( x>=y || IsEqual(x,y) );
687 }
688
689 };
690
691 template<typename T, typename TPointContainer, typename TKeyContainer>
Shrink(size_t nDesiredPoints,size_t nTol,T eTolRight,size_t nMaxIter)692 size_t TLineApproximator<T,TPointContainer,TKeyContainer>::Shrink( size_t nDesiredPoints, size_t nTol, T eTolRight, size_t nMaxIter)
693 {
694 if (m_cPoints.size()<2)
695 return 0;
696
697 // number of points wanted...
698 T dWantedPoints= __min(m_cPoints.size(), __max(2, nDesiredPoints));
699 T uMinWantedPoints = __min(m_cPoints.size(), __max(2, nDesiredPoints-nTol ));
700 T uMaxWantedPoints = __min(m_cPoints.size(), __max(2, nDesiredPoints+nTol ));
701
702 T eLeft, eRight, eMiddle;
703 T dResultLeft, dResultRight;
704 size_t iter=0;
705
706 // compute limits
707 ComputeBoundingBox();
708
709 // normalize if needed
710 if (m_bNormalization)
711 NormalizePoints();
712
713 // first estimation
714 eLeft = 0;
715 SetTol(eLeft);
716
717 ComputeKeys();
718
719 dResultLeft = m_cKeys.size();
720 iter++;
721 // test if success
722 if ( (m_cKeys.size()<=uMaxWantedPoints) && (m_cKeys.size() >= uMinWantedPoints) )
723 goto PostProcess;
724
725 // second estimation
726 if (eTolRight<=0)
727 eRight=__max( m_limits.GetWidth(), m_limits.GetHeight());
728 else
729 eRight=eTolRight;
730 SetTol(eRight);
731
732 ComputeKeys();
733
734 dResultRight = m_cKeys.size();
735
736 // test if optimization possible
737 // if (dResultLeft<uMinWantedPoints || dResultRight>uMaxWantedPoints)
738 // throw _T("TLineApproximator<T>::Shrink failed: Desired compression ratio not possible in the tolerance domain.");
739
740 iter++;
741 // test if success
742 if ( ((m_cKeys.size()<=uMaxWantedPoints) && (m_cKeys.size() >= uMinWantedPoints)) || (dResultLeft == dResultRight) )
743 goto PostProcess;
744
745 // main loop, dichotomy
746 do
747 {
748 // test middle
749 eMiddle=(eLeft +eRight)/2;
750 SetTol(eMiddle);
751
752 // computing new DP...
753 ComputeKeys();
754
755 // updating...
756 if ( (m_cKeys.size()-dWantedPoints)*( dResultLeft-dResultRight) < 0 )
757 {
758 eRight=eMiddle;
759 dResultRight=m_cKeys.size();
760 }
761 else
762 {
763 eLeft=eMiddle;
764 dResultLeft=m_cKeys.size();
765 }
766
767 iter++;
768 } while ( ((m_cKeys.size()>uMaxWantedPoints) || (m_cKeys.size() < uMinWantedPoints)) /* checking that we are in the acceptable compression */
769 && !priv::IsEmptyInterval(eLeft,eRight) /* interval is non empty */
770 && (dResultRight != dResultLeft)
771 && iter<nMaxIter /* checking for maximum number of iterations */);
772
773 PostProcess:
774 if (m_bNormalization)
775 DeNormalizePoints();
776
777 return iter;
778 }
779
780
781 template <typename T, typename TPointContainer, typename TKeyContainer>
Simplify()782 void TLineApproximator<T,TPointContainer, TKeyContainer>::Simplify()
783 {
784 if (m_cPoints.size()<2)
785 return;
786
787 // compute limits
788 ComputeBoundingBox();
789
790 // preprocess...
791 if (m_bNormalization)
792 NormalizePoints();
793
794 ComputeKeys();
795
796 if (m_bNormalization)
797 DeNormalizePoints();
798 }
799
800 }; // namespace hull
801
802 #endif // !defined(AFX_LINEAPPROXIMATOR_H__F5E6E8DC_1185_4AC0_A061_7B3309700E9D__INCLUDED_)
803