Click here to Skip to main content
Click here to Skip to main content

A C++ implementation of Douglas-Peucker Line Approximation Algorithm

, 3 Mar 2003
Rate this:
Please Sign up or sign in to vote.
DP Line approximation algorithm is a well-known method to approximate 2D lines. It is quite fast, O(nlog_2(n)) for a n-points line and can drastically compress a data curve. Here, a fully OOP implementation is given.

Sample Image - dphull.jpg

Introduction

When working with mathematical simulations or engineering problems, it is not unusual to handle curves that contains thousands of points. Usually, displaying all the points is not useful, a number of them will be rendered on the same pixel since the screen precision is finite. Hence, you use a lot of resource for nothing!

This article presents a fast 2D-line approximation algorithm based on the Douglas-Peucker algorithm (see [1]), well-known in the cartography community. It computes a hull, scaled by a tolerance factor, around the curve by choosing a minimum of key points. This algorithm has several advantages:

  • It is fast!: For a n point curve, computation time of the approximation is proportional to nlog_2(n),
  • You don't need a priori knowledge on the curve,
  • The compression ratio can be easily tuned using the tolerance parameter.

The class has been integrated to a plotting library: Plot Graphic Library.

Douglas-Peucker (DP) Line Simplification Algorithm

The DP line simplification algorithm was originally proposed in [1]. John Hershberger and Jack Snoeyink have implemented it in C in [2] as a package named DPSimp:

DPsimp is a Douglas-Peucker line simplification algorithm implementation by John Hershberger and Jack Snoeyink. They analyze the line simplification algorithm reported by Douglas and Peucker and show that its worst case is quadratic in n, the number of input points. The algorithm is based on path hulls, that uses the geometric structure of the problem to attain a worst-case running time proportional to nlog_2(n), which is the best case of the Douglas algorithm.

The algorithm is based on a recursive construction of path hull, as depicted in the picture below. They did all the hard work (writing the base code in C), I ported it to C++ templates.

Modifications to DPSimp

DPSimp was using a recursive call in the DP method. This could lead to a stack overflow when the algorithm would go deep in recursion. To avoid this problem an internal stack has been added to the class to mimic the recursive function call without stack overflow.

Concepts and class design

Let the points denote all the points of the original curve and the keys the points from the original curve that are kept for the approximation.

The idea is that the user provides a container for the points, denoted PointContainer, and for the keys, denoted KeyContainer, and the link between those containers will be the line simplification class, denoted LineApproximator.

How do we build a class hierarchy without restricting ourselves to particular container? In fact, one user might store it's points in vector< pair<T,T>> and another one in 2 separate vector<T>. Of course, the same argument applies to the KeyContainer.

A possible answer is templating. Passing the PointContainer and KeyContainer as template arguments for LineApproximator allows to build the approximation class without specifying the containers, since the class is built at compilation time (We could write interface for those containers but in fact, I'm too lazy for that Smile | :) ).

With this idea in mind, here are the specifications of the container:

PointContainer

Let

  • Point a structure, template or class that has x,y of type T as member,
  • PointRef a structure, template or class that has x,y of type T& as member,

PointContainer behaves like a vector of Point:

  • has clear(), size() and resize() methods,
  • has random access iterator,
  • const_iterator points to a structure similar to Point,
  • iterator points to a structure similar to PointRef
  • operator[] const returns a Point,
  • operator[] returns a PointRef

A simple example of valid PointContainer is

vector<Point>

However, a hybrid container has been developed to handle the case where x and y are in separate containers (See below).

KeyContainer

KeyContainer behaves like a list of PointContainer::const_iterator:

  • has size(), clear() methods,
  • support push_back method

A simple example of valid KeyContainer is

vector<PointContainer::const_iterator>

Again, a hydrid container to handle the case where the keys must be outputted in separate containers is provided.

Templates

All the classes are templated to support float and double version of the algorithm.

The template TDPHull is the user interface to the DP algorithm. However, it relies on a series of subclasses detailed below:

Name Description Use
TLine<T, TPointContainer, TKeyContainer> 2D Line template Points and keys
TLineApproximator<T, TPointContainer, TKeyContainer> 2D Line approximator base class Default interface to approximation algorithms
TDPHull<T, TPointContainer, TKeyContainer> Implementing Douglas-Peukler algorithm User front end
TPathHull<T, TPointContainer, TKeyContainer> Path hull Internal in TDPHull<T>
TPoint<T> A pair of T: x, y Template for 2D point TLineApproximator<T>
TPointRef<T> A pair of T&: x, y Template for 2D point TLineApproximator<T>
SHomog 2D Homogenous point, T x,y,w Internal structure to TLineApproximator<T>

How to use TDPHull?

In the following examples, we adopt the following notations

using namespace hull; // all classes are in the hull namespace,
using namespace std; // using some STL containers

// defining the point container
typedef vector<TPoint<float>> VectorPointContainer; 
// defining the key container
typedef vector<MyPointContainer::const_iterator> ListKeyContainer; 
// defining the line approximation class
typedef TDPHull<float, VectorPointContainer, ListKeyContainer> CDPHullF; 
CDPHullF dp; // a DPHull object, note that you can also work with doubles

The approximation methods throw exception, so you should always enclose them in a try-catch statement.

Normalization

The data points are, by default, normalized before approximation. This is in order to reduce numerical errors in the gradients computations. This happens when the data is badly scaled: using big numbers close together will lead to disastrous loss of significant numbers.

However, if you feel confident about your data, you can disable it by using SetNormalization.

Handling points and keys

Get a reference to the point container and modify it:

// getting reference to container,
TDPHull<float>::PointContainer& pc = dp.GetPoints(); 
for (UINT i=0;i<pc.size();i++)
{
    pc[i].x=...;
    pc[i].y=...;
}

If you are using normalization (default behavior), do not forget to re-compute the data bounding box after your changes:

dp.ComputeBoundingBox();

Approximation tuning

You can control the compression ratio by different ways:

  • Setting the tolerance
    // Setting tolerance of approximation
    dp.SetTol(0.5);
  • Setting a compression ratio, an acceptable compression threshold:
    // dp will find the tolerance corresponding to 10 % of
    // the original points, with 5% of possible error.
    try
    {
        dp.ShrinkNorm(0.1,0.05);
    }
    catch(TCHAR* str)
    {
        // catch and process errors throw by dp
        // example: output to cerr
        cerr<<str<<endl;
    }
    

    The method uses dichotomy to accelerate convergence.

  • Setting the desired number of points, an acceptable number of points threshold:
    // dp will find the tolerance corresponding to 100
    // in the approximated curve, with 5 points of possible error.
    try
    {
        dp.Shrink(100,5);
    }
    catch(TCHAR* str)
    {
        // catch and process errors throw by dp
        ...
    }

Simplifaction

The easiest part of the job:

try
{
    dp.Simplify();
}
catch(TCHAR* str)
{
    // catch and process errors throw by dp
    ...
}

or by using ShrinkNorm, Shrink methods.

Accessing the approximated curve

The keys are stored as PointContainer::const_iterator. You can access the key container by using GetKeys:

// getting conts reference to keys
const TDPHull<float>::KeyContainer& kc = dp.GetKeys(); 
TDPHull<float>::KeyContainer::const_iterator it; // iterator for keys
for (it = kc.begin(); it != kc.end(); it++)
{
    // it is an const_iterator pointing to a PointContainer::const_iterator
    xi=(*it)->x; 
    yi=(*it)->y;
}

Implementing your own algorithm

All you have to do is inherit a class from TLineApproximator and override the function ComputeKeys.

Hydrid containers

You can implement your own containers for points and keys as long as they meet the requirements.

Separate containers for x,y

It is not unusual to have x,y stored in separate containers and moreover these containers are not of the same type. To tackle this problem, two wrapper templates have been written: TPointDoubleContainer and TKeyDoubleContainer which serve as an interface between the approximation algorithms and the containers:

CVectorX vX;    // the vector of x coordinates
CVectorY vY;    // the vector of y coordinates

// defining the hydrid container
typedef TPointDoubleContainer<float, CVectorX, 
                   CVectorY> HybridPointContainer; 

// a list of key x coordinates
CListX lKeyX;    
// a vector of key y coordinates, any container
// that has push_back will do the job :)
CVectorY vKeyY;    
// defining the hydrid container
typedef TKeyDoubleContainer<float, CListX, CVectorY> HybridKeyContainer; 

// creating approximator
TDPHull< float, HybridPointContainer, HydridKeyContainer> dp; 
// attaching point containers
dp.GetPoints().SetContainers( &vX, &vY);    
// attaching key containers
dp.GetKeys().SetContainers( &lKeyX, &vKeyY);    
// dp is now ready to work

Using the demo

The demo shows a real time approximation of a curve by different algorithms.

  • You can stop/start the animation using the toolbar buttons,
  • You can modify the shrinking ration with the scroll bar,
  • You can load your own data with the menu, File->Load Data Set. The file must be formatted with a pair x,y per line.

Using it in your project

Insert the following files in your project and you're done.

LineApproximator.h,DPHull.h, PathHull.h

Known issues

  • The original curve must not self-intersect. This means also that the start and end points must be different, no closed curve !
  • Sick dataset and stack overflow: solved. The problem was due to recursion leading to stack overflow. It is solved now.

Update history

  • 04-03-2003
    • Got rid of DP recursion by adding an internal function call stack. Hence, the stack overflow problem is solved!!!
    • Applied stuff I learned in Effective STL to the classes: using algorithms, functors, etc...
    • Changed class T to typename T
    • Better floating point comparison using boost::close_at_tolerance
  • 15-11-2002
    • More and more templating,
    • Detecting when curve is closed
    • Hybrid containers
    • Fixed bug in compute limits
    • Added LOTS of ASSERT, so Debug version is significantly slower than release build
  • 7-11-2002
    • Fixed a bug in the SLimits structure (GetCenterY)
    • TPathHull working with iterators rather that SPoint*
    • Added exception to handle errors
    • Fixed a bug in TDPHull::ComputeKeys. Was using pc.end() rather that pc.end()--
  • 6-11-2002
    • Added base class TLineApproximator
    • Added proposed algorithm by S.Rog: see TKeyFramer, TGlobalKeyFramer, TDispKeyFramer, TVectKeyFramer
    • Updated demo
  • 3-11-2002
    • Added data normalization for better numerical behavior. Avoids the algorithm to crash when using badly conditioned data. Problem submitted by Corey W.
    • Templated version
    • Got rid of macro and rewrote in more OOP style

References

  1. D. H. Douglas and T. K. Peucker. Algorithms for the reduction of the number of points required to represent a line or its caricature. The Canadian Cartographer, 10(2):112--122, 1973.
  2. J. Hershberger and J. Snoeyink. Speeding up the Douglas-Peucker line simplification algorithm. In Proc. 5th Intl. Symp. Spatial Data Handling. IGU Commission on GIS, pages 134--143, 1992. (home page).

License

This article has no explicit license attached to it but may contain usage terms in the article text or the download files themselves. If in doubt please contact the author via the discussion board below.

A list of licenses authors might use can be found here

About the Author

Jonathan de Halleux
Engineer
United States United States
Jonathan de Halleux is Civil Engineer in Applied Mathematics. He finished his PhD in 2004 in the rainy country of Belgium. After 2 years in the Common Language Runtime (i.e. .net), he is now working at Microsoft Research on Pex (http://research.microsoft.com/pex).

Comments and Discussions

 
GeneralMy vote of 5 Pinmembermanoj kumar choubey15-Feb-12 23:16 
QuestionMakefile? Pinmemberdaviddoria25-Jul-11 3:31 
Generalhow to convert a recursive DP algorithm to non-recursive? Pinmemberunravel-the-sky4-May-11 3:19 
Generalthis template version has bugs Pinmemberyhu2213001-Sep-10 3:33 
GeneralRe: this template version has bugs PinmemberElmar de Koning8-Oct-10 1:17 
GeneralGood job Pinmemberstcornflakes19-Dec-07 0:30 
Questionwhat is the license PinmemberAndreas Fabri13-Nov-07 21:58 
Generalcollinear points [modified] PinmemberRomain Michel12-Mar-07 4:33 
GeneralRe: collinear points PinmemberRalf Wirtz3-Apr-07 20:34 
GeneralRe: collinear points PinmemberRomain Michel11-Apr-07 5:01 
GeneralRe: collinear points Pinmemberdonkey30124-Aug-09 16:41 
QuestionBUG?? PinmemberManiez5-Mar-07 1:55 
AnswerRe: BUG?? Pinmemberdonkey30124-Aug-09 16:44 
Generalanother simple implementation in ANSI C PinmemberAlbertino14-Feb-07 12:04 
Questionany plans for 3d? Pinmembergok25-Jan-06 11:42 
GeneralOriginal C code PinmemberReunion18-Jul-05 19:15 
GeneralRe: Original C code PinmemberJohn M. Drescher18-Jul-05 21:09 
GeneralRe: Original C code PinmemberReunion19-Jul-05 4:04 
GeneralLoop when calling Shrink function PinmemberChristian Staffe8-Aug-04 21:59 
GeneralRe: Loop when calling Shrink function PinmemberCorey W25-Nov-04 11:52 
GeneralRe: Loop when calling Shrink function PinmemberChristian Staffe26-Nov-04 9:20 
GeneralRe: Loop when calling Shrink function PinmemberCorey W26-Nov-04 17:44 
GeneralRe: Loop when calling Shrink function PinmemberCorey W26-Nov-04 17:58 
Generaldouglas line simplification Pinmemberjehovah10-Jun-04 10:24 
QuestionHow do you pronounce Peucker? PinmemberJeff Bogan19-Apr-04 7:55 
AnswerRe: How do you pronounce Peucker? PinsussHans Peucker22-Jun-04 3:46 
Generalusing douglas Pinmembermih221-Jan-04 23:54 
Generalproblem with dphull Pinmembermih221-Jan-04 23:44 
GeneralCompiling with VS.Net 2003 PinmemberCorey W10-Oct-03 11:55 
GeneralRe: Compiling with VS.Net 2003 Pinmembercentomo25-Nov-04 8:54 
GeneralRe: Compiling with VS.Net 2003 PinmemberCorey W25-Nov-04 11:49 
GeneralRe: Compiling with VS.Net 2003 PinmemberChristian Staffe18-May-05 11:44 
I have ported it to Visual Studio 2003 but without much testing. The updated version of the four include files (DPHull.h, LineApproximator.h, PathHull.h and Containers.h) can be found below. Most of the job was limited to the insertion of "typename" keywords at required places.
 
LineApproximator.h
------------------
 
<code>
/* Plot Graphic Library
 
     Copyright (C) 2002 Pelikhan, Inc. All rights reserved
     Go to http://eauminerale.syldavie.csam.ucl.ac.be/peli/pgl/pgl.html for the latest PGL binaries
     This software is NOT freeware.
 
     This software is provided "as is", with no warranty.
     Read the license agreement provided with the files
 
     This software, and all accompanying files, data and materials, are distributed "AS IS" and with no warranties of any kind,
     whether express or implied.   his disclaimer of warranty constitutes an essential part of the agreement.  
     In no event shall Pelikhan, or its principals, shareholders, officers, employees, affiliates, contractors, subsidiaries,
     or parent organizations, be liable for any incidental, consequential, or punitive damages whatsoever relating to the use of PGL,
     or your relationship with Pelikhan.
 
     Author: Jonathan de Halleux, dehalleux@auto.ucl.ac.be
*/
#if !defined(AFX_LINEAPPROXIMATOR_H__F5E6E8DC_1185_4AC0_A061_7B3309700E9D__INCLUDED_)
#define AFX_LINEAPPROXIMATOR_H__F5E6E8DC_1185_4AC0_A061_7B3309700E9D__INCLUDED_
 
#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
 
#include <crtdbg.h>
#ifndef ASSERT
     #ifdef _DEBUG
          #define ASSERT(a)     _ASSERT(a);
     #else
          #define ASSERT(a)
     #endif
#endif
#ifdef min
     #undef min
#endif
#ifdef max
     #undef max
#endif
    
#include <iostream>
#include <vector>
#include <list>
#include <limits>
#include <algorithm>
 
/*!
   \defgroup LAGroup Line approximation algorithms
*/
namespace hull
{
 
/*! \brief A point (x,y)
 
\param T float or double
 
   A pair of (x,y) values.
\ingroup LAGroup
*/
template<typename T>
class TPoint
{
public:
     //! Default constructor
     TPoint( T _x = 0, T _y=0):x(_x), y(_y){};
     //! Assignement constructor
     TPoint& operator = ( const TPoint<T>& t)     {     if (&t != this){ x=t.x; y=t.y;} return *this;};
 
     //! x value
     T x;
     //! y value
     T y;
};
 
/*! \brief A point (x,y) with references
 
\param T float or double
 
   \ingroup LAGroup
*/
template<typename T>
class TPointRef
{
public:
     //! Default contructor
     TPointRef():x(xDummy),y(yDummy){};
     /*! \brief Constructor with 2 values
 
     \param _x value to get reference to.
     \param _y value to get reference to.
     */
     TPointRef(T& _x, T& _y):x(_x),y(_y){};
     //! Assignement constructor
     TPointRef<T>& operator = (const TPointRef<T>& t)     {     if (this != &t){ x=t.x;y=t.y;} return *this;};
     //! Assignement constructor with point
     TPointRef<T>& operator = (TPoint<T> t)          {     x=t.x;y=t.y; return *this;};
     //! x, reference to a value
     T& x;
     //! y, reference to a value
     T& y;
private:
     static T xDummy;
     static T yDummy;
};
 
template<typename T> T TPointRef<T>::xDummy=0;
template<typename T> T TPointRef<T>::yDummy=0;
 
/*! \brief Virtual base class for Line approximator classes
 
\par Template arguments
 
\param T float or double
\param TPointContainer a container of points (structure with x,y) with random access iterators
\param TKeyContainer a container of TPointContainer::const_iterator (structure with x,y) with single direction iterators
 

\par Notations:
 
   - points : data to interpolate
   - keys : selected points from available set that interpolate within the desired tolerance
 
\par Containers:
 
   - Points are stored in a #PointContainer. #PointContainer is a container such that
     #- has random access iterator,
     #- value_type is a structure/class with x,y members
   - Keys are stored in a #KeyContainer. #KeyContainer is a container such that:
     #- has single directional iterator,
     #- value_type is PointContainer::const_iterator
 
\par Data normalization:
 
   To avoid numerical instability when computing cross product and so, data is normalized ( see #NormalizePoints ).
   To enable/disable normalization, use #SetNormalization.
 
   \ingroup LAGroup
*/
template<typename T, typename TPointContainer, typename TKeyContainer>
class TLine
{
public:
 
     //! \name Structures and typedefs
     //@{
     //! Point container
     typedef TPointContainer PointContainer;
     //! Key container
     typedef TKeyContainer KeyContainer;
     //! 2D homogenous point
     struct SHomog
     {
     public:
          //! Default constructor
          SHomog(T _x=0, T _y=0, T _w=1)     { x=_x; y=_y; w=_w;};    
 
          T x;
          T y;
          T w;
     };
 
     /*! \brief returns square of euclidian distance
    
     */
     static T SqrDist( typename const TPointContainer::const_iterator& p, typename const TPointContainer::const_iterator& q)
     {
          T dx=p->x-q->x;
          T dy=p->y-q->y;
          return dx*dx+dy*dy;
     }
 
     /*! \brief Cross product
 
     \param p an iterator with x,y members
     \param q an iterator with x,y members
     \result r cross product of p,q
 
        The euclidian cross product of p,q:
     \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]
     */
     static void CrossProduct( typename const TPointContainer::const_iterator& p, typename const TPointContainer::const_iterator& q, SHomog& r)
     {
               r.w = p->x * q->y - p->y * q->x;
               r.x = - q->y + p->y;
               r.y = q->x - p->x;
     };
         
     /*! \brief Dot product
 
     \param p an iterator with x,y members
     \param q an 2D homogenous point
     \result dot product of p,q
 
     The euclidian dot product of p,q:
     \f[ <\vec p, \vec q> = q_w + p_x q_x + p_y q_y \f]
     */
     static T DotProduct( typename const TPointContainer::const_iterator& p, const SHomog& q)
     {
          return q.w + p->x*q.x + p->y*q.y;
     };
    
     /*! \brief Dot product
 
     \param p an iterator with x,y members
     \param q an iterator with x,y members
     \result dot product of p,q
 
     The euclidian dot product of p,q:
     \f[ < \vec p,\vec q> = p_x q_x + p_y q_y \f]
     */
     static T DotProduct( typename const TPointContainer::const_iterator& p, typename const TPointContainer::const_iterator& q)
     {
          return p->x*q->x + p->y*q->y;
     };
 
     /*! \brief Linear combination of points
 
     \param a p multiplier
     \param p a point
     \param b q multiplier
     \param q a point
     \param r linear combination of p,q
 
     The linear combination is:
     \f[ \vec r = a \vec p + b \vec q \f]
     */
     static void LinComb( T a, typename const TPointContainer::const_iterator& p, T b, typename const TPointContainer::const_iterator& q, typename const TPointContainer::const_iterator& r)
     {
          r->x = a * p->x + b * q->x;
          r->y = a * p->y + b * q->y;
     };
 
     //! Internal limit structure
     struct SLimits
     {
          T dMinX;
          T dMaxX;
          T dMinY;
          T dMaxY;
          T GetCenterX()     {     return static_cast<T>((dMaxX+dMinX)/2.0);};
          T GetCenterY()     {     return static_cast<T>((dMaxY+dMinY)/2.0);};
          T GetWidth()     {     return static_cast<T>(dMaxX-dMinX);};
          T GetHeight()     {     return static_cast<T>(dMaxY-dMinY);};
     };
 
     //! T container
     typedef std::vector<T> TVector;
     //@}
 
     //! Default constructor
     TLine():m_bNormalization(true){};
 
     //! \name Point handling
     //@{
     //! returns number of points
     size_t GetPointSize() const                    {     return m_cPoints.size();};
     //! return vector of points
     PointContainer& GetPoints()                    {     return m_cPoints;};
     //! return vector of points, const
     const PointContainer& GetPoints() const     {     return m_cPoints;};
     //@}
 
     //! \name Key handling
     //@{
     //! returns number of keys
     size_t GetKeySize() const                         {     return m_cKeys.size();};
     //! return keys
     KeyContainer& GetKeys()                         {     return m_cKeys;};
     //! return keys, const
     const KeyContainer& GetKeys() const          {     return m_cKeys;};
     //@}
 
     //! \name Helpers
     //@{
     //! Setting points from vectors
     void SetPoints( const std::vector<T>& vX, const std::vector<T>& vY);
     //! Returning keys to vectors
     void GetKeys( std::vector<T>& vX, std::vector<T>& vY) const;
     //@}
 
     //! \name Bounding box
     //@{
     //! compute the bounding box
     void ComputeBoundingBox();
     //! return the point bounding box
     const SLimits& GetBoundingBox() const          {     return m_limits;};
     //@}
 
     //! \name Normalization
     //@{
     /*! \brief Point normalization
    
     Let \f$(x_i,y_i)\f$, the original points and \f$(\hat x_x, \hat y_i)\f$ the normalized points:
     \f[
     \hat x_i = \frac{x_i - \bar x]}{\max_i (x_i-x_j)}
     \f]
     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$.
 
     \sa DeNormalizePoints
     */
     void NormalizePoints();
 
     /*! \brief Roll back normalization
 
     \sa NormalizePoints
     */
     void DeNormalizePoints();
     //! enabled, disable normalization
     void SetNormalization( bool bEnabled = true)     {     m_bNormalization = true;};
     //! returns true if normalizing
     bool IsNormalization() const                         {     return m_bNormalization;};
     //@}
 
     //! \name Double points checking and loop...
     //@{
     /*! \brief Discard double points
 
     For each pair of points \f$p_i, p_{i+1}\f$, discards \f$p_{i+1}\f$ if
     \f[ \| p_i - p_{i+1} \| < \varepsilon \f]
     */
     void DiscardDoublePoints();
     //! Test for loops
     void FindLoop(size_t uStartPoint, size_t& uEndPoint);
     //@}
 
protected:
     //! \name Attributes
     //@{
     TPointContainer m_cPoints;
     TKeyContainer m_cKeys;
     SLimits m_limits;
     bool m_bNormalization;
     //@}
};
 
namespace priv
{
     template<class Container, class Pred>
     struct PredX : public std::binary_function< typename Container::iterator, typename Container::iterator, bool>
     {
          bool operator()( typename const Container::value_type& p1, typename const Container::value_type& p2)    
          {     return m_pred(p1.x, p2.x); };
     protected:
          Pred m_pred;
     };
 
     template<class Container, class Pred>
     struct PredY : public std::binary_function< typename Container::iterator, typename Container::iterator, bool>
     {
          bool operator()( typename const Container::value_type& p1, typename const Container::value_type& p2)    
          {     return m_pred(p1.y, p2.y); };
     protected:
          Pred m_pred;
     };
};
 
template <typename T, typename TPointContainer, typename TKeyContainer>
void TLine<T, TPointContainer,TKeyContainer>::ComputeBoundingBox()
{
     if (m_cPoints.size() < 2)
          return;
 
     PointContainer::const_iterator it = (*((const TPointContainer*)&m_cPoints)).begin();
 

     //finding minimum and maximum...
     m_limits.dMinX=std::min_element( m_cPoints.begin(), m_cPoints.end(), priv::PredX<TPointContainer, std::less<T> >() )->x ;
     m_limits.dMaxX=std::max_element( m_cPoints.begin(), m_cPoints.end(), priv::PredX<TPointContainer, std::less<T> >() )->x ;
     m_limits.dMinY=std::min_element( m_cPoints.begin(), m_cPoints.end(), priv::PredY<TPointContainer, std::less<T> >() )->y ;
     m_limits.dMaxY=std::max_element( m_cPoints.begin(), m_cPoints.end(), priv::PredY<TPointContainer, std::less<T> >() )->y ;
 
     if ( fabs( m_limits.GetWidth() ) < std::numeric_limits<T>::epsilon() )
     {
          m_limits.dMaxX = m_limits.dMinX+1;
     }
     if ( fabs( m_limits.GetHeight() ) < std::numeric_limits<T>::epsilon() )
     {
          m_limits.dMaxY = m_limits.dMinY+1;
     }
}
 
namespace priv
{
     template<typename T>
     struct Rect
     {
          Rect( T xm, T ym, T dx, T dy)
               :m_xm(xm),m_ym(ym),m_dx(dx),m_dy(dy){};
     protected:
          T m_xm;
          T m_ym;
          T m_dx;
          T m_dy;
     };
 
     template<typename T>
     struct NormalizePoint : public std::unary_function< TPoint<T>& , int>, public Rect<T>
     {
          NormalizePoint( T xm, T ym, T dx, T dy)
               : Rect<T>(xm,ym,dx,dy){};
          int operator() ( TPoint<T>& point)
          {    
               point.x=(point.x-m_xm)/m_dx;
               point.y=(point.y-m_ym)/m_dy;
               return 0;
          };
     };
 
     template<typename T>
     struct DeNormalizePoint : public std::unary_function< TPoint<T>& , int>, public Rect<T>
     {
          DeNormalizePoint( T xm, T ym, T dx, T dy)
               : Rect<T>(xm,ym,dx,dy){};
          int operator() ( TPoint<T>& point)
          {    
               point.x=m_xm+point.x*m_dx;
               point.y=m_ym+point.y*m_dy;
               return 0;
          };
     };
};
 
template <typename T, typename TPointContainer, typename TKeyContainer>
void TLine<T, TPointContainer, TKeyContainer>::NormalizePoints()
{
     T xm,ym,dx,dy;
     // normalizing...
     xm=m_limits.GetCenterX();
     ym=m_limits.GetCenterY();
     dx=m_limits.GetWidth();
     dy=m_limits.GetHeight();
 
     std::for_each(
          m_cPoints.begin(),
          m_cPoints.end(),
          priv::NormalizePoint<T>( xm, ym, dx, dy)
          );
}
 
template <typename T, typename TPointContainer, typename TKeyContainer>
void TLine<T, TPointContainer, TKeyContainer>::DeNormalizePoints()
{
     T xm,ym,dx,dy;
     // normalizing...
     xm=m_limits.GetCenterX();
     ym=m_limits.GetCenterY();
     dx=m_limits.GetWidth();
     dy=m_limits.GetHeight();
 
     std::for_each(
          m_cPoints.begin(),
          m_cPoints.end(),
          priv::DeNormalizePoint<T>( xm, ym, dx, dy)
          );
}
 
template <typename T, typename TPointContainer, typename TKeyContainer>
void TLine<T, TPointContainer, TKeyContainer>::DiscardDoublePoints()
{
     // creating a list...
     TPointContainer& pc=GetPoints();
     TPointContainer::iterator it, it1;
     T epsilon2=std::numeric_limits<T>::epsilon();
 
#ifdef _DEBUG
     size_t count=0;
#endif
    
     it1=it=pc.begin();
     ++it1;
     while (it !=pc.end() && it1!=pc.end())
     {
          if ( SqrDist(it, it1) < epsilon2 )
          {
               it1=pc.erase(it1);
          }
          else
          {
               ++it; ++it1;
          }
     }
 
#ifdef _DEBUG
     TRACE( _T("Numer of (double) points erased: %d\n"), count);
#endif
};
 
template <typename T, typename TPointContainer, typename TKeyContainer>
void TLine<T, TPointContainer, TKeyContainer>::SetPoints( const std::vector<T>& vX, const std::vector<T>& vY)
{
     TPointContainer& pc=GetPoints();
     const size_t n = __min( vX.size(), vY.size());
 
     pc.resize(n);
     for (size_t i=0;i<n;i++)
     {
          pc[i]= TPoint<T>( vX[i], vY[i]);
     }
};
 
template <typename T, typename TPointContainer, typename TKeyContainer>
void TLine<T, TPointContainer, TKeyContainer>::GetKeys( std::vector<T>& vX, std::vector<T>& vY) const
{
     const TKeyContainer& kc=GetKeys();
     TKeyContainer::const_iterator it;
     size_t i;
 
     vX.resize(kc.size());
     vY.resize(kc.size());
 
     for (it=kc.begin(), i=0;it!=kc.end();it++, i++)
     {
          vX[i]=(*it)->x;
          vY[i]=(*it)->y;
     }
};
 
template <typename T, typename TPointContainer, typename TKeyContainer>
void TLine<T, TPointContainer, TKeyContainer>::FindLoop(size_t uStartPoint, size_t& uEndPoint)
{
};
 

/*! \brief Base class for line approximators
 
    
 
   \ingroup LAGroup
*/
template <typename T, typename TPointContainer, typename TKeyContainer>
class TLineApproximator : virtual public TLine<T,TPointContainer, TKeyContainer>
{
public:
 
     //! \name Constructor
     //@{
     TLineApproximator(): m_dTol(0)
     {m_limits.dMinX=m_limits.dMinY=0;m_limits.dMaxX=m_limits.dMaxY=1;};
     ~TLineApproximator(){};
     //@}
 

     //! \name Tolerance
     //@{
     //! sets the tolerance
     void SetTol( double dTol)                    {     m_dTol = __max( dTol, 0);};
     //! return current tolerance
     double GetTol() const                         {     return m_dTol;};
     //@}
 
     //! \name Simplification functions
     //@{
     //! Initialize simplification
     void ClearKeys()                                        {     m_cKeys.clear();};
     //! Compute the keys
     void Simplify();
     /*! Shrink to compression level
    
     \param dScale scaling to apply [0...1]
     \param dScaleTol [optional] tolerance with respect to dScale, default is 0.05
     \param eTolRight [optional] first estimate on right tolerance
     \param nMaxIter [optional] maximum number of iterations, default is 250
     \return number of estimations
     */
     size_t ShrinkNorm( T dScale, T dScaleTol = 0.05, T eTolRight = -1, size_t nMaxIter = 250);
 
     /*! Shrink to a specified number of points
    
     \param n desired number of points in the approximate curve
     \param nTol [optional] tolerance with respect to n, default is 10
     \param eTolRight [optional] first estimate on right tolerance
     \param nMaxIter [optional] maximum number of iterations, default is 250
     \return number of estimations
     */
     size_t Shrink( size_t nDesiredPoints, size_t nTol = 10, T eTolRight = -1, size_t nMaxIter = 250);
     //@}
 
protected:
     //! \name Virtual functions
     //@{
     /*! \brief Virtual approximation function
 
     This function must be overriden in inherited classes. To implement your own algorithm,
     override this function.
     */
     virtual void ComputeKeys()          {     ClearKeys();};
     //@}
 
private:
     T m_dTol;
};
 

template <typename T, typename TPointContainer, typename TKeyContainer>
size_t TLineApproximator<T,TPointContainer,TKeyContainer>::ShrinkNorm( T dScale, T dScaleTol, T eTolRight ,size_t nMaxIter)
{
     // number of points wanted...
     size_t uWantedPoints= __min(m_cPoints.size(), __max(2, static_cast<size_t>(floor(m_cPoints.size()*dScale))));
     size_t uTol = __min(m_cPoints.size(), __max(0, static_cast<size_t>(floor(m_cPoints.size()*dScaleTol)) ));
 
     return Shrink( uWantedPoints, uTol, eTolRight, nMaxIter);
}
 
namespace priv
{
//   (C) Copyright Gennadiy Rozental 2001-2002.
//   Permission to copy, use, modify, sell and distribute this software
//   is granted provided this copyright notice appears in all copies.
//   This software is provided "as is" without express or implied warranty,
//   and with no claim as to its suitability for any purpose.
 
//   See http://www.boost.org for most recent version including documentation.
//
//   File            : $RCSfile: floating_point_comparison.hpp,v $
//
//   Version      : $Id: floating_point_comparison.hpp,v 1.6 2002/09/16 08:47:29 rogeeff Exp $
//
//   Description : defines algoirthms for comparing 2 floating point values
// ***************************************************************************
     template<typename FPT>
     inline FPT
     fpt_abs( FPT arg )
     {
           return arg < 0 ? -arg : arg;
     }
 
     // both f1 and f2 are unsigned here
     template<typename FPT>
     inline FPT
     safe_fpt_division( FPT uf1, FPT uf2 )
     {
           return   ( uf1 < 1 && uf1 > uf2 * std::numeric_limits<FPT>::max())  
               ? std::numeric_limits<FPT>::max() :
                    ((uf2 > 1 && uf1 < uf2 * std::numeric_limits<FPT>::min() ||
                       uf1 == 0)                                                                     ? 0                                             :
                                                                                                              uf1/uf2 );
     }
 
     template<typename FPT>
     class close_at_tolerance
     {
     public:
           explicit close_at_tolerance( FPT tolerance, bool strong_or_weak = true )
               : p_tolerance( tolerance ),m_strong_or_weak( strong_or_weak ) { };
    
     explicit      close_at_tolerance( int number_of_rounding_errors, bool strong_or_weak = true )
           : p_tolerance( std::numeric_limits<FPT>::epsilon() * number_of_rounding_errors/2 ),
        m_strong_or_weak( strong_or_weak ) {}
    
           bool            operator()( FPT left, FPT right ) const
           {
                 FPT diff = fpt_abs( left - right );
           FPT d1   = safe_fpt_division( diff, fpt_abs( right ) );
           FPT d2   = safe_fpt_division( diff, fpt_abs( left ) );
                
                 return m_strong_or_weak ? (d1 <= p_tolerance.get() && d2 <= p_tolerance.get())
                                               : (d1 <= p_tolerance.get() || d2 <= p_tolerance.get());
           }
    
           // Data members
          class p_tolerance_class
          {
          private:
               FPT f;
          public:
               p_tolerance_class(FPT _f=0):f(_f){};
               FPT   get() const{     return f;};
          };
          p_tolerance_class p_tolerance;    
     private:
     bool            m_strong_or_weak;
     };
 
     template <typename T>
     inline bool IsEqual(T x, T y)         
     {
        static close_at_tolerance<T> comp( std::numeric_limits<T>::epsilon()/2*10);
        return comp(fpt_abs(x),fpt_abs(y));
     };
  
     template <typename T>
     inline bool IsEmptyInterval(T x, T y)
     {
          return ( x>=y || IsEqual(x,y) );
     }
 
};
 
template<typename T, typename TPointContainer, typename TKeyContainer>
size_t TLineApproximator<T,TPointContainer,TKeyContainer>::Shrink( size_t nDesiredPoints, size_t nTol, T eTolRight, size_t nMaxIter)
{
     if (m_cPoints.size()<2)
          return 0;
    
     // number of points wanted...
     T dWantedPoints= __min(m_cPoints.size(), __max(2, nDesiredPoints));
     T uMinWantedPoints = __min(m_cPoints.size(), __max(2, nDesiredPoints-nTol ));
     T uMaxWantedPoints = __min(m_cPoints.size(), __max(2, nDesiredPoints+nTol ));
 
     T eLeft, eRight, eMiddle;
     T dResultLeft, dResultRight;
     size_t iter=0;
 
     // compute limits
     ComputeBoundingBox();
 
     // normalize if needed
     if (m_bNormalization)
          NormalizePoints();
 
     // first estimation
     eLeft = 0;
     SetTol(eLeft);
 
     ComputeKeys();
 
     dResultLeft =   m_cKeys.size();
     iter++;
     // test if success
     if ( (m_cKeys.size()<=uMaxWantedPoints) && (m_cKeys.size() >= uMinWantedPoints)   )
          goto PostProcess;
 
     // second estimation
     if (eTolRight<=0)
          eRight=__max( m_limits.GetWidth(), m_limits.GetHeight());
     else
          eRight=eTolRight;
     SetTol(eRight);
 
     ComputeKeys();
 
     dResultRight =   m_cKeys.size();
 
     // test if optimization possible
//     if (dResultLeft<uMinWantedPoints ||   dResultRight>uMaxWantedPoints)
//          throw _T("TLineApproximator<T>::Shrink failed: Desired compression ratio not possible in the tolerance domain.");
 
     iter++;
     // test if success
     if ( ((m_cKeys.size()<=uMaxWantedPoints) && (m_cKeys.size() >= uMinWantedPoints)) || (dResultLeft == dResultRight) )
          goto PostProcess;
 
     // main loop, dichotomy
     do
     {
          // test middle
          eMiddle=(eLeft +eRight)/2;
          SetTol(eMiddle);
 
          // computing new DP...
          ComputeKeys();
         
          // updating...
          if ( (m_cKeys.size()-dWantedPoints)*( dResultLeft-dResultRight) < 0 )
          {
               eRight=eMiddle;
               dResultRight=m_cKeys.size();
          }
          else
          {
               eLeft=eMiddle;
               dResultLeft=m_cKeys.size();
          }
 
          iter++;
     } while ( ((m_cKeys.size()>uMaxWantedPoints) || (m_cKeys.size() < uMinWantedPoints)) /* checking that we are in the acceptable compression */
          && !priv::IsEmptyInterval(eLeft,eRight) /* interval is non empty */
          && (dResultRight != dResultLeft)
          && iter<nMaxIter /* checking for maximum number of iterations */);
 
PostProcess:
     if (m_bNormalization)
          DeNormalizePoints();
 
     return iter;
}
 

template <typename T, typename TPointContainer, typename TKeyContainer>
void TLineApproximator<T,TPointContainer, TKeyContainer>::Simplify()
{
     if (m_cPoints.size()<2)
          return;
 
     // compute limits
     ComputeBoundingBox();
    
     // preprocess...
     if (m_bNormalization)
          NormalizePoints();
 
     ComputeKeys();
 
     if (m_bNormalization)
          DeNormalizePoints();
}
 
}; // namespace hull
 
#endif // !defined(AFX_LINEAPPROXIMATOR_H__F5E6E8DC_1185_4AC0_A061_7B3309700E9D__INCLUDED_)
</code>
 
DPHull.h
--------
 
<code>
 
#if !defined(AFX_DPHULL1_H__6CE88E63_3AC7_4E18_87FB_CACF5BE62BE4__INCLUDED_)
#define AFX_DPHULL1_H__6CE88E63_3AC7_4E18_87FB_CACF5BE62BE4__INCLUDED_
 
#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
 
#ifndef UINT
     #define UINT unsigned int
#endif
 
#include <tchar.h>
#include <vector>
#include <stack>
#include "PathHull.h"
#include "KeyExporter.h"
#include "Containers.h"
 
namespace hull
{
 
/*! \brief Douglas-Peukler Appromixation algorithm
 
\ingroup LAGroup
*/
template <typename T, typename TPointContainer, typename TKeyContainer>
class TDPHull : public TLineApproximator<T,TPointContainer,TKeyContainer>
{
public:
     //! Build step types
     enum EBuildStep
     {
          //! Output to vertex
          BuildStepOutputVertex,
          //! Call DP
          BuildStepDP,
          //! Call Build
          BuildStepBuild,
          //! Is a return key
          BuildStepReturnKey
     };
     //! A build step structure
     struct SBuildStep
     {
          SBuildStep(typename TPointContainer::const_iterator _i,typename TPointContainer::const_iterator _j, EBuildStep _s):i(_i),j(_j),s(_s){};
 
          typename TPointContainer::const_iterator i;
          typename TPointContainer::const_iterator j;
          EBuildStep s;
     };
 
public:
     //! \name Constructors
     //@{
     //! Default constructor
     explicit TDPHull(){};
     //! Destructor
     virtual ~TDPHull(){};
     //@}
 
protected:
     //! A path hull
     typedef TPathHull<T,TPointContainer,TKeyContainer> PathHull;
     typedef std::stack< SBuildStep > BuildStack;
 
     //! Computes the keys
     virtual void ComputeKeys();
 
     //! \name Hull methods
     //@{
     /*! \brief Adds a function call to the stack
 
     \param i left point iterator
     \param j right point iterator
     \param buildStep step description
     */
     void AddBuildStep( typename TPointContainer::const_iterator i, typename TPointContainer::const_iterator j, EBuildStep buildStep)
     {     m_stack.push( BuildStack::value_type(i,j,buildStep) );     }
     /*! \brief Apply Douglas-Peucker algo
 
     \pre m_stack not empty
     \pre m_stack.top().s == BuildStepDP
     */
     void DP();
 
     /*! \brief Builds the path hull
 
     \pre m_stack not empty
     \pre m_stack.top().s == BuildStepBuild
     */
     void Build();
 
     /*! \brief Stores key
 
     \pre m_stack size = 2
     \pre m_stack.top().s == BuildStepReturnKey
     \pre m_stack.top()(twice).s == BuildStepOutputVertex
     */
     void OutputVertex()
     {
          ASSERT(!m_stack.empty());
          ASSERT(m_stack.top().s==BuildStepReturnKey);
 
          GetKeys().push_back(m_stack.top().i);
          m_stack.pop();
          ASSERT(!m_stack.empty());
          ASSERT(m_stack.top().s==BuildStepOutputVertex);
          m_stack.pop();
     };
     //@}
 
protected:
     //! \name Attributes
     //@{
     TPathHull<T,TPointContainer,TKeyContainer> m_phRight;
     TPathHull<T,TPointContainer,TKeyContainer> m_phLeft;
     typename TPointContainer::const_iterator m_pPHtag;
 
     BuildStack m_stack;
     //@}
};
 
template<typename T, typename TPointContainer, typename TKeyContainer>
void TDPHull<T,TPointContainer,TKeyContainer>::DP()
     {
          T ld, rd, len_sq;
          SHomog l;
          register TPointContainer::const_iterator le;
          register TPointContainer::const_iterator re;
         
          ASSERT( !m_stack.empty() );
          ASSERT( m_stack.top().s == BuildStepDP );
          TPointContainer::const_iterator i(m_stack.top().i);
          TPointContainer::const_iterator j(m_stack.top().j);
          m_stack.pop();
 
          CrossProduct(i, j, l);
          len_sq = l.x * l.x + l.y * l.y;
         
          if (j - i < 8)
          {          /* chain small */
               rd   = 0.0;
               for (le = i + 1; le < j; ++le)
               {
                    ld = DotProduct(le, l);
                    if (ld < 0) ld = - ld;
                    if (ld > rd)
                    {
                         rd = ld;
                         re = le;
                    }
               }
               if (rd * rd > GetTol() * len_sq)
               {
                    AddBuildStep( re, j, BuildStepDP );
                    AddBuildStep( i, re, BuildStepOutputVertex );
                    AddBuildStep( i, re, BuildStepDP );
                    return;
//                    OutputVertex(DP(i, re));
//                    return(DP(re, j));
               }
               else
               {
                    AddBuildStep(j,j,BuildStepReturnKey);
                    return;
//                    return j;
               }
          }
          else
          {    
               /* chain large */
               int sbase, sbrk, mid, lo, m1, brk, m2, hi;
               T d1, d2;
               if ((m_phLeft.GetTop() - m_phLeft.GetBot()) > 8)
               {
                    /* left hull large */
                    lo = m_phLeft.GetBot();
                    hi = m_phLeft.GetTop() - 1;
                    sbase = m_phLeft.SlopeSign(hi, lo, l);
                    do
                    {
                         brk = (lo + hi) / 2;
                         if (sbase == (sbrk = m_phLeft.SlopeSign(brk, brk+1, l)))
                              if (sbase == (m_phLeft.SlopeSign(lo, brk+1, l)))
                                   lo = brk + 1;
                              else
                                   hi = brk;
                    }
                    while (sbase == sbrk && lo < hi);
                   
                    m1 = brk;
                    while (lo < m1)
                    {
                         mid = (lo + m1) / 2;
                         if (sbase == (m_phLeft.SlopeSign(mid, mid+1, l)))
                              lo = mid + 1;
                         else
                              m1 = mid;
                    }
                   
                    m2 = brk;
                    while (m2 < hi)
                    {
                         mid = (m2 + hi) / 2;
                         if (sbase == (m_phLeft.SlopeSign(mid, mid+1, l)))
                              hi = mid;
                         else
                              m2 = mid + 1;
                    };
                   
                   
                    if ((d1 = DotProduct(m_phLeft.GetpElt(lo), l)) < 0) d1 = - d1;
                    if ((d2 = DotProduct(m_phLeft.GetpElt(m2), l)) < 0) d2 = - d2;
                    ld = (d1 > d2 ? (le = m_phLeft.GetpElt(lo), d1) : (le = m_phLeft.GetpElt(m2), d2));
               }
               else
               {               /* Few SPoints in left hull */
                    ld = 0.0;
                    for (mid = m_phLeft.GetBot(); mid < m_phLeft.GetTop(); mid++)
                    {
                         if ((d1 = DotProduct(m_phLeft.GetpElt(mid), l)) < 0) d1 = - d1;
                         if (d1 > ld)
                         {
                              ld = d1;
                              le = m_phLeft.GetpElt(mid);
                         }
                    }
               }
              
               if ((m_phRight.GetTop() - m_phRight.GetBot()) > 8)
               {               /* right hull large */
                    lo = m_phRight.GetBot(); hi = m_phRight.GetTop() - 1;
                    sbase = m_phRight.SlopeSign(hi, lo, l);
                    do
                    {
                         brk = (lo + hi) / 2;
                         if (sbase == (sbrk = m_phRight.SlopeSign(brk, brk+1, l)))
                              if (sbase == (m_phRight.SlopeSign(lo, brk+1, l)))
                                   lo = brk + 1;
                              else
                                   hi = brk;
                    }
                    while (sbase == sbrk && lo < hi);
                   
                    m1 = brk;
                    while (lo < m1)
                    {
                         mid = (lo + m1) / 2;
                         if (sbase == (m_phRight.SlopeSign(mid, mid+1, l)))
                              lo = mid + 1;
                         else
                              m1 = mid;
                    }
                   
                    m2 = brk;
                    while (m2 < hi)
                    {
                         mid = (m2 + hi) / 2;
                         if (sbase == (m_phRight.SlopeSign(mid, mid+1, l)))
                              hi = mid;
                         else
                              m2 = mid + 1;
                    };
                   
                    if ((d1 = DotProduct(m_phRight.GetpElt(lo), l)) < 0) d1 = - d1;
                    if ((d2 = DotProduct(m_phRight.GetpElt(m2), l)) < 0) d2 = - d2;
                    rd = (d1 > d2 ? (re = m_phRight.GetpElt(lo), d1) : (re = m_phRight.GetpElt(m2), d2));
               }
               else
               {               /* Few SPoints in righthull */
                    rd = 0.0;
                    for (mid = m_phRight.GetBot(); mid < m_phRight.GetTop(); mid++)
                    {
                         if ((d1 = DotProduct(m_phRight.GetpElt(mid), l)) < 0) d1 = - d1;
                         if (d1 > rd)
                         {
                              rd = d1;
                              re = m_phRight.GetpElt(mid);
                         }
                    }
               }
      }
    
    
     if (ld > rd)
          if (ld * ld > GetTol() * len_sq)
          {
               /* split left */
               register int tmpo;
              
               while ((m_phLeft.GetHp() >= 0)
                    && ( (tmpo = m_phLeft.GetOps()[m_phLeft.GetHp()] ),
                    ((re = m_phLeft.GetpHelt(m_phLeft.GetHp())) != le) || (tmpo != PathHull::StackPushOp)))
               {
                    m_phLeft.DownHp();
                    switch (tmpo)
                    {
                    case PathHull::StackPushOp:
                         m_phLeft.DownTop();
                         m_phLeft.UpBot();
                         break;
                    case PathHull::StackTopOp:
                         m_phLeft.UpTop();
                         m_phLeft.SetTopElt(re);
                         break;
                    case PathHull::StackBotOp:
                         m_phLeft.DownBot();
                         m_phLeft.SetBotElt(re);
                         break;
                    }
               }
              
               AddBuildStep( le, j, BuildStepDP );
               AddBuildStep( le, j, BuildStepBuild);
               AddBuildStep( i, le, BuildStepOutputVertex);
               AddBuildStep( i, le, BuildStepDP);
               AddBuildStep( i, le, BuildStepBuild);
               return;
//               Build(i, le);
//               OutputVertex(DP(i, le));
//               Build(le, j);
//               return DP(le, j);
          }
          else
          {
               AddBuildStep(j,j,BuildStepReturnKey);
               return;
//               return(j);
          }
          else                    /* extreme on right */
               if (rd * rd > GetTol() * len_sq)
               {                    /* split right or both */
//                    if (m_pPHtag == re)
//                    {
//                         AddBuildStep( i, re, BuildStepBuild);
//                         Build(i, re);
//                    }
//                    else
                    if (m_pPHtag != re)
                    {    
                         /* split right */
                         register int tmpo;
                        
                         while ((m_phRight.GetHp() >= 0)
                              && ((tmpo = m_phRight.GetOps()[m_phRight.GetHp()]),
                              ((le = m_phRight.GetpHelt(m_phRight.GetHp())) != re) || (tmpo != PathHull::StackPushOp)))
                         {
                              m_phRight.DownHp();
                              switch (tmpo)
                              {
                              case PathHull::StackPushOp:
                                   m_phRight.DownTop();
                                   m_phRight.UpBot();
                                   break;
                              case PathHull::StackTopOp:
                                   m_phRight.UpTop();
                                   m_phRight.SetTopElt(le);
                                   break;
                              case PathHull::StackBotOp:
                                   m_phRight.DownBot();
                                   m_phRight.SetBotElt(le);
                                   break;
                              }
                         }
                    }
 
                    AddBuildStep( re, j ,BuildStepDP );
                    AddBuildStep( re, j, BuildStepBuild );
                    AddBuildStep( i, re, BuildStepOutputVertex );
                    AddBuildStep( i, re, BuildStepDP );
                    if (m_pPHtag == re)
                         AddBuildStep( i, re, BuildStepBuild);
 
                    return;
 
//                    OutputVertex(DP(i, re));
//                    Build(re, j);
//                    return(DP(re, j));
               }
               else
                    AddBuildStep( j,j, BuildStepReturnKey);
//                    return(j);    
     }
 
template<typename T, typename TPointContainer, typename TKeyContainer>
void TDPHull<T,TPointContainer,TKeyContainer>::Build()
{
          TPointContainer::const_iterator k;
          int topflag, botflag;
         
          ASSERT( !m_stack.empty() );
          ASSERT( m_stack.top().s == BuildStepBuild );
          TPointContainer::const_iterator i(m_stack.top().i);
          TPointContainer::const_iterator j(m_stack.top().j);
          m_stack.pop();
 
          m_pPHtag = i + (j - i) / 2;     /* Assign tag vertex */
         
         
          m_phLeft.Init(m_pPHtag, m_pPHtag - 1); /* \va{left} hull */
          for (k = m_pPHtag - 2; k >= i; --k)
          {
               topflag = m_phLeft.LeftOfTop(k);
               botflag = m_phLeft.LeftOfBot(k);
               if ((topflag || botflag) && !(topflag && botflag))
               {
                    while (topflag)
                    {
                         m_phLeft.PopTop();
                         topflag = m_phLeft.LeftOfTop(k);
                    }
                    while (botflag)
                    {
                         m_phLeft.PopBot();
                         botflag = m_phLeft.LeftOfBot(k);
                    }
                    m_phLeft.Push(k);
               }
          }
         
          m_phRight.Init(m_pPHtag, m_pPHtag + 1); /* \va{right} hull */
          for (k = m_pPHtag + 2; k <= j; ++k)
          {
               topflag = m_phRight.LeftOfTop(k);
               botflag = m_phRight.LeftOfBot(k);
               if ((topflag || botflag) && !(topflag && botflag))
               {
                    while (topflag)
                    {
                         m_phRight.PopTop();
                         topflag = m_phRight.LeftOfTop(k);
                    }
                    while (botflag)
                    {
                         m_phRight.PopBot();
                         botflag = m_phRight.LeftOfBot(k);
                    }
                    m_phRight.Push(k);
               }
          }
     };
 

template<typename T, typename TPointContainer, typename TKeyContainer>
void TDPHull<T,TPointContainer,TKeyContainer>::ComputeKeys()
{
     using namespace std;
     static const T epsilon2 = numeric_limits<T>::epsilon()*numeric_limits<T>::epsilon();
     TLineApproximator<T,TPointContainer,TKeyContainer>::ComputeKeys();
     const TPointContainer& pc=GetPoints();
     TPointContainer::const_iterator pcend=pc.end();
     --pcend;
     T len_sq;
     SHomog l;
 
     /////////////////////////////////////////////////////////////////////////////
     CrossProduct(pc.begin(), pcend, l);
 
     len_sq = l.x * l.x + l.y * l.y;
     if (len_sq < epsilon2)
          throw _T("TDPHull<T,TKeyExporter>::DP failed: Start and end points are equal or at a distance < epsilon\n");
 
     ////////////////////////////////////////////////////////////////////////
     // prepraring path if needed...
     m_phLeft.SetMaxSize(pc.size()+1);
     m_phRight.SetMaxSize(pc.size()+1);
 
     /////////////////////////////////////////////////////////////////////////
     // building hull
//     Build(pc.begin(), pcend);     /* Build the initial path hull */    
//     OutputVertex( pc.begin() );
//     OutputVertex( DP(pc.begin(), pcend) ); /* Simplify */
 
     AddBuildStep( pc.begin(), pcend, BuildStepBuild );
     Build();
     AddBuildStep( pc.begin(), pc.begin(), BuildStepOutputVertex );
     AddBuildStep( pc.begin(), pc.begin(), BuildStepReturnKey );
     OutputVertex();
 
     AddBuildStep( pc.begin(), pc.begin(), BuildStepOutputVertex );
     AddBuildStep( pc.begin(), pcend, BuildStepDP );
     while (!m_stack.empty())
     {
//          std::cerr<<"stack size: "<<m_stack.size()<<std::endl;
          switch( m_stack.top().s)
          {
          case BuildStepDP:
               DP();
               break;
          case BuildStepBuild:
               Build();
               break;
          case BuildStepReturnKey:
               OutputVertex();    
               break;
          case BuildStepOutputVertex:
               ASSERT(false);
               break;
          default:
               ASSERT(false);
          }
     }
 
     //////////////////////////////////////////////////////////////////////////
     // cleaning path...
     m_phLeft.SetMaxSize(0);
     m_phRight.SetMaxSize(0);
};
 
/*! \brief A single precision DPHull
 
The classical DPHull object:
     - PointContainer: vector<TPoint<float>>
     - KeyContainer: list<PointContainer::const_iterator>
 
   \ingroup LAGroup
*/
typedef TDPHull<float, std::vector< TPoint<float> >, std::vector< std::vector< TPoint< float > >::const_iterator > > CDPHullF;
 
/*! \brief A double precision DPHull
 
The classical DPHull object:
     - PointContainer: vector<TPoint<double>>
     - KeyContainer: list<PointContainer::const_iterator>
 
   \ingroup LAGroup
*/
typedef TDPHull<double, std::vector< TPoint< double > >, std::vector< std::vector< TPoint< double   > >::const_iterator > > CDPHullD;
 
};
 
#endif // !defined(AFX_DPHULL1_H__6CE88E63_3AC7_4E18_87FB_CACF5BE62BE4__INCLUDED_)
</code>
 
PathHull.h
----------
 
<code>
// PathHull.h: interface for the CPathHull class.
//
//////////////////////////////////////////////////////////////////////
 
#if !defined(AFX_PATHHULL_H__50C639BA_585B_4272_9AF4_4632128D8938__INCLUDED_)
#define AFX_PATHHULL_H__50C639BA_585B_4272_9AF4_4632128D8938__INCLUDED_
 
#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
 
#include "LineApproximator.h"
 
namespace hull
{
 

#define DP_SGN(a) (a >= 0)
 
/*! \brief A path
     \ingroup
*/
template<typename T, typename TPointContainer, typename TKeyContainer>
class TPathHull  
{
public:
     enum EStackOp
     {
          StackPushOp=0,
          StackTopOp=1,
          StackBotOp=2
     };
 
     typedef std::vector<signed char> OpContainer;
     typedef std::vector<typename TPointContainer::const_iterator> PCItContainer;
 
     TPathHull(): m_iHullMax(0)
     {};
     virtual ~TPathHull()
     {};
 
     void SetMaxSize(int iHullMax);
 
     int GetHp() const                                                            {     return m_iHp;};
     int GetBot() const                                                            {     return m_iBot;};
     int GetTop() const                                                            {     return m_iTop;};
     typename const TPointContainer::const_iterator& GetpElt(int i)     const     {     ASSERT(i<m_ppElt.size()); return m_ppElt[i];};
     typename const TPointContainer::const_iterator& GetpHelt(int i)     const     {     ASSERT(i<m_ppHelt.size()); return m_ppHelt[i];};
     OpContainer& GetOps()                                                       {     return m_pOp;};
 
     void SetHp(int hp)                                                            {     m_iHp=hp;};
 
     void UpHp()          {     m_iHp++;};
     void UpTop()     {     m_iTop++;};
     void UpBot()     {     m_iBot++;};
     void DownHp()     {     m_iHp--;};
     void DownTop()     {     m_iTop--;};
     void DownBot()     {     m_iBot--;};
 
     void SetTopElt(typename const TPointContainer::const_iterator& p)          {     ASSERT(m_iTop>=0); ASSERT(m_iTop<m_ppElt.size()); m_ppElt[m_iTop]=p;};
     void SetBotElt(typename const TPointContainer::const_iterator& p)          {     ASSERT(m_iBot>=0); ASSERT(m_iBot<m_ppElt.size()); m_ppElt[m_iBot]=p;};
 
     void Split(typename const TPointContainer::const_iterator& e)
     {
          TPointContainer::const_iterator tmpe;
          int tmpo;
         
          ASSERT(m_iHp<m_ppHelt.size());
          ASSERT(m_iHp<m_pOp.size());
          while ((m_iHp >= 0)
               && ((tmpo = m_pOp[m_iHp]),
               ((tmpe = m_ppHelt[m_iHp]) != e) || (tmpo != StackPushOp)))
          {
               m_iHp--;
               switch (tmpo)
               {
               case StackPushOp:
                    m_iTop--;
                    m_iBot++;
                    break;
               case StackTopOp:
                    ASSERT(m_iTop-1>=0);
                    ASSERT(m_iTop+1<m_ppElt.size());
                    m_ppElt[++m_iTop] = tmpe;
                    break;
               case StackBotOp:
                    ASSERT(m_iBot-1>=0);
                    ASSERT(m_iBot-1<m_ppElt.size());
                    m_ppElt[--m_iBot] = tmpe;
                    break;
               }
          }
     }
 
     void FindExtreme(typename const TLine<T,TPointContainer,TKeyContainer>::SHomog& line, typename TPointContainer::const_iterator* e, T& dist)
     {
     int sbase, sbrk, mid,lo, m1, brk, m2, hi;
     T d1, d2;
    
     if ((m_iTop - m_iBot) > 8)
      {
          lo = m_iBot; hi = m_iTop - 1;
          sbase = SlopeSign(hi, lo, line);
          do
          {
               brk = (lo + hi) / 2;
               if (sbase == (sbrk = SlopeSign(brk, brk+1, line)))
                    if (sbase == (SlopeSign(lo, brk+1, line)))
                         lo = brk + 1;
                    else
                         hi = brk;
          }
          while (sbase == sbrk);
         
          m1 = brk;
          while (lo < m1)
          {
               mid = (lo + m1) / 2;
               if (sbase == (SlopeSign(mid, mid+1, line)))
                    lo = mid + 1;
               else
                    m1 = mid;
          }
         
          m2 = brk;
          while (m2 < hi)
          {
               mid = (m2 + hi) / 2;
               if (sbase == (SlopeSign(mid, mid+1, line)))
                    hi = mid;
               else
                    m2 = mid + 1;
          }
         
          ASSERT(lo>=0);
          ASSERT(lo <m_ppElt.size());
          if ((d1 = TLine<T,TPointContainer,TKeyContainer>::DotProduct(*m_ppElt[lo], line)) < 0)
               d1 = - d1;
 
          ASSERT(m2>=0);
          ASSERT(m2 <m_ppElt.size());
          if ((d2 = TLine<T,TPointContainer,TKeyContainer>::DotProduct(*m_ppElt[m2], line)) < 0)
               d2 = - d2;
         
          dist = (d1 > d2 ? (*e = m_ppElt[lo], d1) : (*e = m_ppElt[m2], d2));
      }
     else                    /* Few DP_POINTs in hull */
      {
          dist = 0.0;
          for (mid = m_iBot; mid < m_iTop; mid++)
          {
              
               ASSERT(mid>=0);
               ASSERT(mid<m_ppElt->size());
               if ((d1 = TLine<T,TPointContainer,TKeyContainer>::::DotProduct(*m_ppElt[mid], line)) < 0)
                    d1 = - d1;
               if (d1 > *dist)
               {
                    dist = d1;
                    *e = m_ppElt[mid];
               }
          }
      }
}
 
     void Init(typename const TPointContainer::const_iterator& e1, typename const TPointContainer::const_iterator& e2)
     {
          /* Initialize path hull and history   */
          ASSERT(m_iHullMax>=0);
          ASSERT(m_iHullMax+1<m_ppElt.size());
          m_ppElt[m_iHullMax] = e1;
          m_ppElt[m_iTop = m_iHullMax + 1] =
               m_ppElt[m_iBot = m_iHullMax - 1] =
               m_ppHelt[m_iHp = 0] = e2;
          m_pOp[0] = StackPushOp;
     }
 
     void Push(typename const TPointContainer::const_iterator& e)
     {
          ASSERT(m_iTop+1 >= 0);
          ASSERT(m_iTop+1 < m_ppElt.size());
          ASSERT(m_iBot-1 >= 0);
          ASSERT(m_iBot-1 < m_ppElt.size());
          ASSERT(m_iHp+1 >= 0);
          ASSERT(m_iHp+1 < m_ppHelt.size());
          ASSERT(m_iHp+1 < m_pOp.size());
 
          /* Push element $e$ onto path hull $h$ */
          m_ppElt[++m_iTop] = m_ppElt[--m_iBot] = m_ppHelt[++m_iHp] = e;
          m_pOp[m_iHp] = StackPushOp;
     }
 
     void PopTop()
     {    
          ASSERT(m_iTop >= 0);
          ASSERT(m_iTop < m_ppElt.size());
          ASSERT(m_iHp+1 >= 0);
          ASSERT(m_iHp+1 < m_ppHelt.size());
          ASSERT(m_iHp+1 < m_pOp.size());
 
          m_ppHelt[++m_iHp] = m_ppElt[m_iTop--];
          m_pOp[m_iHp] = StackTopOp;
     }
 
     void PopBot()
     {
          ASSERT(m_iBot >= 0);
          ASSERT(m_iBot < m_ppElt.size());
          ASSERT(m_iHp+1 >= 0);
          ASSERT(m_iHp+1 < m_ppHelt.size());
          ASSERT(m_iHp+1 < m_pOp.size());
 
          /* Pop from bottom */
          m_ppHelt[++m_iHp] = m_ppElt[m_iBot++];
          m_pOp[m_iHp] = StackBotOp;
     }
    
     void Add(typename const TPointContainer::const_iterator& p)
     {
          int topflag, botflag;
         
          topflag = LeftOfTop(p);
          botflag = LeftOfBot(p);
         
          if (topflag || botflag)
          {
               while (topflag)
               {
                    PopTop();
                    topflag = LeftOfTop(p);
               }
               while (botflag)
               {
                    PopBot();
                    botflag = LeftOfBot(p);
               }
               Push(p);
          }
     }
 
     int LeftOfTop(typename const TPointContainer::const_iterator& c)
     {
          ASSERT(m_iTop >= 1);
          ASSERT(m_iTop < m_ppElt.size());
 
          /* Determine if point c is left of line a to b */
          return (((*m_ppElt[m_iTop]).x - (*c).x)*((*m_ppElt[m_iTop-1]).y - (*c).y)
               >= ((*m_ppElt[m_iTop-1]).x - (*c).x)*((*m_ppElt[m_iTop]).y - (*c).y));
     }
 
     int LeftOfBot(typename const TPointContainer::const_iterator& c)
     {
          ASSERT(m_iBot >= 0);
          ASSERT(m_iBot+1 < m_ppElt.size());
 
          /* Determine if point c is left of line a to b */
          return (((*m_ppElt[m_iBot+1]).x - (*c).x)*((*m_ppElt[m_iBot]).y - (*c).y)
               >= ((*m_ppElt[m_iBot]).x - (*c).x)*((*m_ppElt[m_iBot+1]).y - (*c).y));
     }
 

     int SlopeSign(int p, int q, typename const TLine<T,TPointContainer,TKeyContainer>::SHomog& l)
     {
          ASSERT(p >= 0);
          ASSERT(p < m_ppElt.size());
          ASSERT(q >= 0);
          ASSERT(q < m_ppElt.size());
 
          /* Return the sign of the projection
                       of $h[q] - h[p]$ onto the normal
                       to line $l$ */
          return (int) (DP_SGN(
               (l.x)*((*m_ppElt[q]).x - (*m_ppElt[p]).x)
               + (l.y)*((*m_ppElt[q]).y - (*m_ppElt[p]).y) ) ) ;
     };
 
protected:
     /// Maxium number of elements in hull
     int m_iHullMax;
    
     /// internal values
     int m_iTop;
     int m_iBot;
     int m_iHp;
     OpContainer m_pOp;    
     PCItContainer m_ppElt;
     PCItContainer m_ppHelt;
};
 
template <typename T,typename TPointContainer,typename TKeyContainer>
void TPathHull<T,TPointContainer,TKeyContainer>::SetMaxSize(int iHullMax)
{
     if (m_iHullMax == iHullMax)
          return;
 
     m_iHullMax=iHullMax;
    
     m_pOp.resize(3*m_iHullMax);
     m_ppElt.resize(2*m_iHullMax);
     m_ppHelt.resize(3*m_iHullMax);
}
 

};
 
#endif // !defined(AFX_PATHHULL_H__50C639BA_585B_4272_9AF4_4632128D8938__INCLUDED_)
</code>
 
Containers.h
------------
 
<code>
// Containers.h: interface for the CContainers class.
//
//////////////////////////////////////////////////////////////////////
 
#if !defined(AFX_CONTAINERS_H__33EAD324_F802_4BC8_BEC8_39E69C1C21A7__INCLUDED_)
#define AFX_CONTAINERS_H__33EAD324_F802_4BC8_BEC8_39E69C1C21A7__INCLUDED_
 
#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
 
#include <vector>
#include <list>
 
namespace hull
{
 
/*!
     \defgroup DCGroup Double container groups
     \ingroup LAGroup
*/
 
/*! \brief Virtual base class for Point double container
 
\ingroup DCGroup    
*/
template<typename T, typename TVectorX, typename TVectorY>
class TPointDoubleContainerBase
{
public:
     //! \name Typedefs
     //@{
     typedef TPoint<T> value_type;
     typedef TPointRef<T> reference;
     typedef TVectorX container_x;
     typedef TVectorY container_y;
     //@}
 
     /*! \brief Constructor
 
     \param pX pointer to x coordinate container,
     \param pY pointer to y coordinate container
     */
     TPointDoubleContainerBase(container_x* pX = NULL, container_y* pY =NULL):m_pX(pX),m_pY(pY){};
     virtual ~TPointDoubleContainerBase(){};
    
     //! return the size of the point container
     size_t size() const               {     ASSERT(m_pX); ASSERT(m_pY); return __min( m_pX->size(), m_pY->size());};
     /*! \brief resize the container
    
     \param size the new size of the vector.
     */
     void resize(size_t size)     {     ASSERT(m_pX); ASSERT(m_pY); m_pX->resize(size); m_pY->resize(size);};
 
     //! return the value at position i
     value_type operator[](UINT i) const     {     return value_type((*m_pX)[i],(*m_pY)[i]);};
     //! return a reference to position i
     reference operator[](UINT i)          {     return reference((*m_pX)[i],(*m_pY)[i]);};
 
     //! Sets the container pointers
     void SetContainers( container_x* pX, container_y* pY)     {     m_pX=pX; m_pY=pY;};
     //! return the x coordinate pointer, const
     const container_x* GetXContainer() const          {     return m_pX;};
     //! return the y coordinate pointer, const
     const container_y* GetYContainer() const          {     return m_pY;};
     //! return the x coordinate pointer
     container_x* GetXContainer()                         {     return m_pX;};
     //! return the y coordinate pointer
     container_y* GetYContainer()                         {     return m_pY;};
 
protected:
     container_x* m_pX;
     container_y* m_pY;
};
 
/*! \brief Base class for point iterator
 
\param T float or double
\param TVectorX x coordinate container with random access iterator
\param TVectorY y coordinate container with random access iterator
 
     A random access iterator base class.
 
   \ingroup DCGroup
*/
template<typename T,typename TVectorX, typename TVectorY>
class TPointIt
{
public:
     TPointIt(const TPointDoubleContainerBase<T,TVectorX,TVectorY>* pV = NULL, UINT index = 0):m_pV(pV),m_index(index){};
     virtual ~TPointIt(){};
 
     TPointIt<T,TVectorX,TVectorY>& operator = (const TPointIt<T,TVectorX,TVectorY>& t)    
     {    
          if (&t!=this)
          {
               m_pV=t.m_pV;
               m_index=t.m_index;
          };
          return *this;
     };
 
     friend UINT operator - (const TPointIt<T,TVectorX,TVectorY>& t1, const TPointIt<T,TVectorX,TVectorY>& t2){     return t1.m_index-t2.m_index; };
 
     friend bool operator == ( const TPointIt<T,TVectorX,TVectorY>& t1, const TPointIt<T,TVectorX,TVectorY>& t2)     {     return t1.m_index==t2.m_index && t1.m_pV==t2.m_pV;};
     friend bool operator != ( const TPointIt<T,TVectorX,TVectorY>& t1, const TPointIt<T,TVectorX,TVectorY>& t2)     {     return   t1.m_index!=t2.m_index || t1.m_pV!=t2.m_pV;};
     friend bool operator < ( const TPointIt<T,TVectorX,TVectorY>& t1, const TPointIt<T,TVectorX,TVectorY>& t2)     {     return   t1.m_index<t2.m_index;};
     friend bool operator > ( const TPointIt<T,TVectorX,TVectorY>& t1, const TPointIt<T,TVectorX,TVectorY>& t2)     {     return   t1.m_index>t2.m_index;};
     friend bool operator <= ( const TPointIt<T,TVectorX,TVectorY>& t1, const TPointIt<T,TVectorX,TVectorY>& t2)     {     return   t1.m_index<=t2.m_index;};
     friend bool operator >= ( const TPointIt<T,TVectorX,TVectorY>& t1, const TPointIt<T,TVectorX,TVectorY>& t2)     {     return   t1.m_index>=t2.m_index;};
 
protected:
     long m_index;
     const TPointDoubleContainerBase<T,TVectorX,TVectorY>* m_pV;
};
 
/*! \brief PointDoubleContainer iterator
 
\param T float or double
\param TVectorX x coordinate container with random access iterator
\param TVectorY y coordinate container with random access iterator
 
   A random access iterator.
 
   \ingroup DCGroup
*/
template<typename T, typename TVectorX, typename TVectorY>
class TPointRandIt : public TPointIt<T,TVectorX,TVectorY>
{
public:
     TPointRandIt(const TPointDoubleContainerBase<T,TVectorX,TVectorY>* pV = NULL, UINT index = 0):TPointIt<T,TVectorX,TVectorY>(pV,index)
     {};
     TPointRandIt<T,TVectorX,TVectorY>& operator = (const TPointRandIt<T,TVectorX,TVectorY>& t)    
     {    
          if (&t!=this)
          {
               TPointIt<T,TVectorX,TVectorY>::operator=(t);
          };
          return *this;
     };
 

     TPoint<T> operator*() const              
     {    
          ASSERT(m_pV);
          ASSERT(m_index<m_pV->size());
          return TPoint<T>((*m_pV->GetXContainer())[m_index], (*m_pV->GetYContainer())[m_index]);
     };
 
     TPointRef<T> operator*()                             
     {    
          ASSERT(m_pV);
          ASSERT(m_index<m_pV->size());
          return TPointRef<T>( (*m_pV->GetXContainer())[m_index], (*m_pV->GetYContainer())[m_index]);
     };
 
     TPointRandIt<T,TVectorX,TVectorY>& operator++(int)                    {     return TPointRandIt<T,TVectorX,TVectorY>( m_pV, m_index+1); };
     TPointRandIt<T,TVectorX,TVectorY>& operator--(int)                    {     return TPointRandIt<T,TVectorX,TVectorY>( m_pV, m_index-1); };
     TPointRandIt<T,TVectorX,TVectorY>& operator++()                         {     m_index++; return *this; };
     TPointRandIt<T,TVectorX,TVectorY>& operator--()                         {     m_index--; return *this; };
     TPointRandIt<T,TVectorX,TVectorY>& operator+=(UINT i)               {     m_index+=i; return *this; };
     TPointRandIt<T,TVectorX,TVectorY>& operator-=(UINT i)               {     m_index-=i;   return *this; };
     friend TPointRandIt<T,TVectorX,TVectorY> operator + (const TPointRandIt<T,TVectorX,TVectorY>& t, UINT i)     {     return TPointRandIt<T,TVectorX,TVectorY>( t.m_pV, t.m_index+i);};
     friend TPointRandIt<T,TVectorX,TVectorY> operator - (const TPointRandIt<T,TVectorX,TVectorY>& t, UINT i)     {     return TPointRandIt<T,TVectorX,TVectorY>( t.m_pV, t.m_index-i);};
 
};
 
/*! \brief PointDoubleContainer const_iterator
 
\param T float or double
\param TVectorX x coordinate container with random access iterator
\param TVectorY y coordinate container with random access iterator
 
     A const random access iterator.
 
   \ingroup DCGroup
*/
template<typename T, typename TVectorX, typename TVectorY>
class TPointConstRandIt : public TPointIt<T,TVectorX,TVectorY>
{
public:
     TPointConstRandIt(const TPointDoubleContainerBase<T,TVectorX,TVectorY>* pV = NULL, UINT index = 0):TPointIt<T,TVectorX,TVectorY>(pV,index){};
 
     TPointConstRandIt& operator = (const TPointConstRandIt<T,TVectorX,TVectorY>& t)    
     {    
          if (&t!=this)
          {
               TPointIt<T,TVectorX,TVectorY>::operator=(t);
          };
          return *this;
     };
 
     TPointConstRandIt& operator = (const TPointRandIt<T,TVectorX,TVectorY>& t)    
     {    
          TPointIt<T,TVectorX,TVectorY>::operator=(t);
          return *this;
     };
 
     TPoint<T> operator*() const                                  
     {              
          ASSERT(m_pV);
          ASSERT(m_index<m_pV->size());    
          return TPoint<T>((*m_pV->GetXContainer())[m_index], (*m_pV->GetYContainer())[m_index]);
     };
 
     TPointConstRandIt<T,TVectorX,TVectorY>& operator++(int)                    {     return TPointConstRandIt<T,TVectorX,TVectorY>( t.m_pV, m_index+1); };
     TPointConstRandIt<T,TVectorX,TVectorY>& operator--(int)                    {     return TPointConstRandIt<T,TVectorX,TVectorY>( t.m_pV, m_index-1); };
     TPointConstRandIt<T,TVectorX,TVectorY>& operator++()                         {     m_index++; return *this; };
     TPointConstRandIt<T,TVectorX,TVectorY>& operator--()                         {     m_index--; return *this; };
     TPointConstRandIt<T,TVectorX,TVectorY>& operator+=(UINT i)               {     m_index+=i; return *this; };
     TPointConstRandIt<T,TVectorX,TVectorY>& operator-=(UINT i)               {     m_index-=i;   return *this; };
     friend TPointConstRandIt<T,TVectorX,TVectorY> operator + (const TPointConstRandIt<T,TVectorX,TVectorY>& t, UINT i)     {     return TPointConstRandIt<T,TVectorX,TVectorY>( t.m_pV, t.m_index+i);};
     friend TPointConstRandIt<T,TVectorX,TVectorY> operator - (const TPointConstRandIt<T,TVectorX,TVectorY>& t, UINT i)     {     return TPointConstRandIt<T,TVectorX,TVectorY>( t.m_pV, t.m_index-i);};
};
 
/*! \brief A container linking two separate containers into a point container
 
\param T float or double
\param TVectorX x coordinate container with random access iterator
\param TVectorY y coordinate container with random access iterator
 
   \ingroup DCGroup
*/
template<typename T, typename TVectorX, typename TVectorY>
class TPointDoubleContainer : virtual public TPointDoubleContainerBase<T,TVectorX,TVectorY>
{
public:
     TPointDoubleContainer(TVectorX* pX = NULL, TVectorY* pY =NULL):TPointDoubleContainerBase<T,TVectorX,TVectorY>(pX,pY){};
     virtual~TPointDoubleContainer(){};
 
     typedef TPointConstRandIt<T,TVectorX,TVectorY> const_iterator;
     typedef TPointRandIt<T,TVectorX,TVectorY> iterator;
 
     iterator begin()     {     return iterator(this,0); };
     iterator end()          {     return iterator(this,size()); };
 
     const_iterator begin() const     {     return const_iterator(this,0); };
     const_iterator end() const          {     return const_iterator(this,size()); };
};
 
/*! \brief A container to export point to two containers
 
\param T float or double
\param TListX x coordinate container with single direction iterator
\param TListY y coordinate container with single direction iterator
\param TPointContainer The point container
 
     \ingroup DCGroup
*/
template<typename T, typename TListX, typename TListY, typename TPointContainer>
class TKeyDoubleContainer
{
public:
     typedef typename TPointContainer::const_iterator value_type;
     typedef TListX container_x;
     typedef TListY container_y;
 
     TKeyDoubleContainer( container_x* pListX = NULL, container_y* pListY = NULL):m_pListX(pListX), m_pListY(pListY){};
 
     void SetContainers( container_x* pListX, container_y* pListY)          {     m_pListX=pListX; m_pListY=pListY;};
     container_x* GetXContainer()                                                  {     return m_pListX;};
     container_y* GetYContainer()                                                  {     return m_pListY;};
     const container_x* GetXContainer() const                                   {     return m_pListX;};
     const container_y* GetYContainer() const                                   {     return m_pListY;};
 
     void clear()          {     ASSERT(m_pListX); ASSERT(m_pListY); m_pListX->clear(); m_pListY->clear();};
     size_t size() const     {     ASSERT(m_pListX); ASSERT(m_pListY); return __min( m_pListX->size(), m_pListY->size());};
 
     void push_back( const value_type& p)     { ASSERT(m_pListX); ASSERT(m_pListY); m_pListX->push_back((*p).x); m_pListY->push_back((*p).y);};
     void push_front( const value_type& p)     { ASSERT(m_pListX); ASSERT(m_pListY); m_pListX->push_front((*p).x); m_pListY->push_front((*p).y);};
 
protected:
     container_x* m_pListX;
     container_y* m_pListY;
};
 

 
};
 
#endif // !defined(AFX_CONTAINERS_H__33EAD324_F802_4BC8_BEC8_39E69C1C21A7__INCLUDED_)
</code>

GeneralConfused Pinmemberashnyr25-Sep-03 22:48 
Generalanother implementation Pinsussfaraz05-Aug-03 2:00 
GeneralRe: another implementation PinmemberJonathan de Halleux5-Aug-03 5:24 
GeneralWell done PinmemberBitterman4-Mar-03 13:19 
GeneralSQL PinmemberJonathan de Halleux5-Mar-03 8:14 
GeneralRe: SQL PinmemberCarsten Breum21-Mar-03 1:44 
GeneralDiving on hold PinmemberJonathan de Halleux27-Mar-03 6:43 
GeneralRe: Diving on hold PinmemberCarsten Breum27-Mar-03 20:30 
GeneralRe: Diving on hold PinmemberJonathan de Halleux27-Mar-03 22:18 
GeneralRe: Diving on hold PinmemberCeox14-Oct-03 23:37 
GeneralRe: SQL PinmemberJonathan de Halleux28-Jul-03 13:26 
GeneralExamples PinmemberChopper4-Mar-03 2:53 
GeneralRe: Examples PinmemberCorey W4-Mar-03 4:27 
GeneralRe: Examples PinmemberChopper5-Mar-03 4:46 
GeneralSick datatypes PinmemberJonathan de Halleux12-Mar-03 7:27 
GeneralRe: Sick datatypes PinmemberCorey W12-Mar-03 7:34 
GeneralCool PinmemberJonathan de Halleux13-Mar-03 0:46 
GeneralRe: Cool PinmemberCorey W13-Mar-03 5:23 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Mobile
Web02 | 2.8.140721.1 | Last Updated 4 Mar 2003
Article Copyright 2002 by Jonathan de Halleux
Everything else Copyright © CodeProject, 1999-2014
Terms of Service
Layout: fixed | fluid