Click here to Skip to main content
15,886,060 members
Search within:


Articles / Programming Languages / C++
 

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

Rate me:
Please Sign up or sign in to vote.
4.87/5 (39 votes)
3 Mar 20037 min read 421.6K   13.8K   126  
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.
 
#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(TPointContainer::const_iterator _i,TPointContainer::const_iterator _j, EBuildStep _s):i(_i),j(_j),s(_s){};

		TPointContainer::const_iterator i;
		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( TPointContainer::const_iterator i, 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;
	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_)

By viewing downloads associated with this article you agree to the Terms of Service and the article's licence.

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

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


Written By
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

Permalink
Advertise
Privacy
Cookies
Terms of Use
Layout: fixed | fluid

Posted 3 Jan 2002

Article Copyright 2002 by Jonathan de Halleux
Everything else Copyright © CodeProject, 1999-2024

Web03 2.8:2024-04-02:1