Click here to Skip to main content
Click here to Skip to main content
Add your own
alternative version

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

, 3 Mar 2003
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.
// 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<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;};
	const TPointContainer::const_iterator& GetpElt(int i)	const	{	ASSERT(i<m_ppElt.size()); return m_ppElt[i];};
	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(const TPointContainer::const_iterator& p)		{	ASSERT(m_iTop>=0); ASSERT(m_iTop<m_ppElt.size()); m_ppElt[m_iTop]=p;};
	void SetBotElt(const TPointContainer::const_iterator& p)		{	ASSERT(m_iBot>=0); ASSERT(m_iBot<m_ppElt.size()); m_ppElt[m_iBot]=p;};

	void Split(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(const TLine<T,TPointContainer,TKeyContainer>::SHomog& line, 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(const TPointContainer::const_iterator& e1, 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(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(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(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(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, 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_)

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

Share

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).

| Advertise | Privacy | Terms of Use | Mobile
Web03 | 2.8.150327.1 | Last Updated 4 Mar 2003
Article Copyright 2002 by Jonathan de Halleux
Everything else Copyright © CodeProject, 1999-2015
Layout: fixed | fluid