2D Convex Hull Calculation

ConvexPolygon.cpp

// REPLACE WITH STACK ALLOCATOR !!!

#include "ConvexPolygon.h"
#include <sdla/SDLWrap.h>
#include <algorithm>


// This is a non-randomised algorithm that sorts the input points so that each new point added
// to the convex polygon is outside the one already defined. Normally this would be grossly
// inefficient (something randomised algorithms attempt to remedy), but the sorting based approach
// allows the MergePlanar phase to be simplified leading to an O(n log n) algorithm.
struct PointCompare
{
	bool operator () (const sdla::cVector& a, const sdla::cVector& b) const
	{
		if ((a.x < b.x) || (a.x == b.x && a.y < b.y))
			return (true);
		return (false);
	}
};


// Circular linked list node
// 16 + 8 = 24 bytes
// a 10k vert model would be ~240k
struct cConvexPolygon::Vertex
{
	// Construct from vector position
	Vertex(const sdla::cVector& v) :
		position(v)
	{
		// Make sure that the link forms its own loop
		prev = this;
		next = this;
	}

	// Vertex position
	sdla::cVector	position;

	// Links to previous and next vertices in the linked list
	Vertex*	prev;
	Vertex*	next;
};


cConvexPolygon::cConvexPolygon(const int max_nb_vertices) :

	m_Type(INVALID),
	m_FirstVertex(0),
	m_LastVertex(0)
	
{
	// Create the allocator with the worst-case entry count
	m_Allocator = new sdla::tSmallAllocator<Vertex>(max_nb_vertices);
}


cConvexPolygon::~cConvexPolygon(void)
{
	delete m_Allocator;
}


void cConvexPolygon::SortInput(std::vector<sdla::cVector>& vertices)
{
	// SGI STL uses the IntroSort algorithm which is, on average, at least as fast as QuickSort.
	// Worst case complexity for IntroSort is O(n log n), whereas QuickSort is quadratic.
	// Pretty nifty for a standard library function!
	std::sort(vertices.begin(), vertices.end(), PointCompare());
}


void cConvexPolygon::Merge(const sdla::cVector& p)
{
	switch (m_Type)
	{
		// Addition of the first point
		case (INVALID):
			m_FirstVertex = m_LastVertex = new (*m_Allocator) Vertex(p);
			m_Type = POINT;
			return;

		// Transition from point hull to linear hull is guaranteed when all
		// duplicate vertices have been removed.
		case (POINT):
			PushAfter(m_LastVertex, p);
			m_Type = LINEAR;
			return;

		// Linear hull case, possible change from linear to convex polygon
		case (LINEAR):
			MergeLinear(p);
			return;

		// Final state is a convex polygon
		case (PLANAR):
			MergePlanar(p);
			return;
	}
}


void cConvexPolygon::Draw(void) const
{
	// Enable state for drawing lines
	glDisable(GL_LIGHTING);
	glDisable(GL_DEPTH_TEST);

	// Draw the polygon!
	glBegin(GL_LINES);
	glLineWidth(3);
	Vertex* vptr = m_FirstVertex;
	while (vptr)
	{
		glVertex3fv(vptr->position);
		glVertex3fv(vptr->next->position);

		vptr = vptr->next;

		if (vptr == m_FirstVertex)
			break;
	}
	glEnd();

	// Restore state
	glEnable(GL_DEPTH_TEST);
	glEnable(GL_LIGHTING);
}


void cConvexPolygon::MergeLinear(const sdla::cVector& p)
{
	sdla::cVector q0 = m_FirstVertex->position;
	sdla::cVector q1 = m_LastVertex->position;

	switch (ColinearTest(p, q0, q1))
	{
		// Changing from linear hull to convex polygon
		case (POSITIVE):
			// Construct CCW polygon <p, q0, q1>
			PushBefore(m_FirstVertex, p);
			m_Type = PLANAR;
			return;

		// Changing from linear hull to convex polygon
		case (NEGATIVE):
			// Construct CCW polygon <p, q1, q0>
			Remove(m_FirstVertex);
			PushBefore(m_FirstVertex, p);
			PushAfter(m_LastVertex, q0);
			m_Type = PLANAR;
			return;

		// Point on the same line, in the order <p, q0, q1>
		case (COLINEAR_LEFT):
			// Construct the linear hull <p, q1>
			Remove(m_FirstVertex);
			PushBefore(m_FirstVertex, p);
			return;
		
		// Point on the same line, in the order <q0, q1, p>
		case (COLINEAR_RIGHT):
			// Construct the linear hull <q0, p>
			Remove(m_LastVertex);
			PushAfter(m_LastVertex, p);
			return;

		// Point on the same line, in the order <q0, p, q1>
		case (COLINEAR_CONTAIN):
			// Hull remains unchanged
			return;
	}
}


void cConvexPolygon::MergePlanar(const sdla::cVector& p)
{
	Vertex*	U;
	Vertex*	L;

	// Begin both searches from the last inserted point
	// Due to the pre-process sort the point to be added is guaranteed to be
	// outside the convex polygon so there is no need to check for this

	// Find upper tangent point by scanning CCW along each edge for the first non visible edge
	for (U = m_FirstVertex; ; U = U->next)
	{
		// Check to see if current edge is visible to this point
		PointPlacement placement = ColinearTest(p, U->position, U->next->position);

		// Edge is visible, skip to the next one
		if (placement == NEGATIVE)
			continue;

		// Edge not visible, the first point of the segment is now the upper tangent point
		if (placement == POSITIVE || placement == COLINEAR_LEFT)
			break;

		// Theoretically not possible with a distinct and sorted pointset.
		// Might happen due to FPU inaccuracies.
		// Assume P is on the hull polygon so return.
		return;
	}

	// Find lower tangent point by scanning CW along each edge for the first non visible edge
	for (L = m_FirstVertex; ; L = L->prev)
	{
		// Check to see if current edge is visible to this point
		PointPlacement placement = ColinearTest(p, L->prev->position, L->position);

		// Edge is visible, skip to the next one
		if (placement == NEGATIVE)
			continue;

		// Edge not visible, the second point of the segment is now the lower tangent point
		if (placement == POSITIVE || placement == COLINEAR_RIGHT)
			break;

		// Theoretically not possible with a distinct and sorted pointset.
		// Might happen due to FPU inaccuracies.
		// Assume P is on the hull polygon so return.
		return;
	}

	// Allocate the vertex to be inserted
	Vertex* P = new (*m_Allocator) Vertex(p);

	// Completely remove the visible edges by linking the lower tangent vertex to P...
	L->next = P;
	P->prev = L;

	// ...and P to the upper tangent vertex
	P->next = U;
	U->prev = P;

	// Ensure the last inserted vertex is always the first in the list
	m_FirstVertex = P;
	m_LastVertex = U;
}


cConvexPolygon::PointPlacement cConvexPolygon::ColinearTest(const sdla::cVector& p, const sdla::cVector& q0, const sdla::cVector& q1)
{
	// Line direction
	sdla::cVector d = q1 - q0;

	// The line normal is the negative perpendicular to the line direction
	sdla::cVector n = -sdla::cVector(-d.y, d.x, 0);

	// Vector from line segment start point to p
	sdla::cVector a = p - q0;

	// p is on the line if n and a are co-linear
	float nda = n * a;
	if (nda > 0)
		return (POSITIVE);
	if (nda < 0)
		return (NEGATIVE);

	float dda = d * a;

	// d and a are directed in opposite directions so p is on the left of q0
	if (dda < 0)
		return (COLINEAR_LEFT);

	// p is on the right if d.a is larger than the square magnitude of the direction
	// (only true if d and a are directed equally, which is true at this point)
	if (dda > d * d)
		return (COLINEAR_RIGHT);

	return (COLINEAR_CONTAIN);
}


void cConvexPolygon::PushBefore(Vertex* before, const sdla::cVector& v)
{
	// Allocate the vertex
	Vertex* vertex = new (*m_Allocator) Vertex(v);

	// Link in the existing chain
	before->prev->next = vertex;
	vertex->next = before;
	vertex->prev = before->prev;
	before->prev = vertex;

	// Assign as first?
	if (before == m_FirstVertex)
		m_FirstVertex = vertex;
}


void cConvexPolygon::PushAfter(Vertex* after, const sdla::cVector& v)
{
	// Allocate the vertex
	Vertex* vertex = new (*m_Allocator) Vertex(v);

	// Link in the existing chain
	after->next->prev = vertex;
	vertex->prev = after;
	vertex->next = after->next;
	after->next = vertex;

	// Assign as last?
	if (after == m_LastVertex)
		m_LastVertex = vertex;
}


void cConvexPolygon::Remove(Vertex* vertex)
{
	// Remove from the linked list
	vertex->prev->next = vertex->next;
	vertex->next->prev = vertex->prev;

	// Update the first and last vertices
	if (m_FirstVertex == vertex) m_FirstVertex = vertex->next;
	else if (m_LastVertex == vertex) m_LastVertex = vertex->prev;

	// Release the memory
	sdla::cVector v = vertex->position;
}

ConvexPolygon.h

#pragma once


#include <vector>
#include <sdla/Vector.h>
#include <sdla/SmallAllocator.h>


class cConvexPolygon
{
public:
	enum Type
	{
		// Not yet constructed
		INVALID,

		// The hull is one point
		POINT,

		// The hull is two points denoting a line
		LINEAR,

		// The hull is a convex polygon
		PLANAR
	};

	// Viewing the line segment <q0, q1>, a given point p can be placed anywhere in R2
	// and fall under the following placement classifications.
	enum PointPlacement
	{
		// Given the left-facing normal of <q0, q1>, p is on the positive side of the segment
		POSITIVE,

		// Given the left-facing normal of <q0, q1>, p is on the negative side of the segment
		NEGATIVE,

		// p is on the line through which the segment <q0, q1> exists, and to the left of q0
		COLINEAR_LEFT,

		// p is on the line through which the segment <q0, q1> exists, and to the right of q1
		COLINEAR_RIGHT,

		// p is on the line through which the segment <q0, q1> exists, and within the segment bounds
		COLINEAR_CONTAIN
	};

	// Default constructor
	explicit cConvexPolygon(const int max_nb_vertices);
	~cConvexPolygon(void);

	// Sort the input with the required less-than operator for correct polygon construction
	static void SortInput(std::vector<sdla::cVector>& vertices);

	// Merge a new point with the convex polygon
	void	Merge(const sdla::cVector& p);

	// Draw the polygon as a set of lines
	void	Draw(void) const;

private:
	struct Vertex;

	void			MergeLinear(const sdla::cVector& p);
	void			MergePlanar(const sdla::cVector& p);
	PointPlacement	ColinearTest(const sdla::cVector& p, const sdla::cVector& q0, const sdla::cVector& q1);
	void			PushBefore(Vertex* before, const sdla::cVector& v);
	void			PushAfter(Vertex* after, const sdla::cVector& v);
	void			Remove(Vertex* vertex);

	// Type of convex polygon
	Type	m_Type;

	// Allocator for type vertex type
	sdla::tSmallAllocator<Vertex>*	m_Allocator;

	// Polygon first and last point
	Vertex*	m_FirstVertex;
	Vertex*	m_LastVertex;
};