[comp.graphics] Wanted: C Code for Graham Scan Convex Hull Algorithm...

rehm@cs.washington.edu (04/30/91)

Does anyone have C code for the Graham Scan Convex Hull Algorithm?

(I'm additionally interested in any implementation that avoids the use
of polar angles...)

howardl@sgi.com (Howard Look) (05/01/91)

In article <1991Apr30.053933.5561@beaver.cs.washington.edu> rehm@cs.washington.edu writes:
>Does anyone have C code for the Graham Scan Convex Hull Algorithm?
>
>(I'm additionally interested in any implementation that avoids the use
>of polar angles...)

I'm not really sure who Graham is, but here's code that works for me.
It was taken from Sedgewick's "Algorithms". The copy I was using has
been taken by its owner (the nerver) so I'm afraid I don't have any more
source information than that.

Hope this helps,
Howard.


/* 
Finds the convex hull of the given points. Returns the number
of points in the hull. The first and last point will NOT be duplicated.

Uses a wrapping algorithm as described in Sedgewick's Algorithms.
*/
int find_hull(Point3d *point, int num, Point3d *p)
{
	int min = 0;
	int i, hull_points = -1, done = FALSE;
	double current_angle, min_angle, angle;

	/* copy the points into the result and do the algorithm in place */
	for (i = 0; i < num; i++)
	{
		p[i][0] = point[i][0];
		p[i][1] = point[i][1];
	}
	
	/* find the index of the point lowest and leftmost */
	for (i = 0; i < num; i++)
	{
		if (p[i][1] < p[min][1])
			min = i;
		else if (p[i][1] == p[min][1])
		{
			if (p[i][0] < p[min][0])
				min = i;
		}
	}


	min_angle = 0.0;
	
	do
	{
		double t;

		hull_points++;

		/* swap points */
		if (min != hull_points)
		{
			t = p[hull_points][0];
			p[hull_points][0] = p[min][0];
			p[min][0] = t;
		
			t = p[hull_points][1];
			p[hull_points][1] = p[min][1];
			p[min][1] = t;
		}

		current_angle = min_angle;
		min_angle = 360.0;

		for (i = hull_points + 1; i <= num; i++)
		{
			/* first time through, dont stop */
			if ( (hull_points == 0) && (i == num))
				break;
				
			if (i != num)
				angle = theta(p[hull_points], p[i]);
			else
				angle = theta(p[hull_points], p[0]);

			if (( angle >= current_angle) && ( angle <= min_angle))
			{
				min = i;
				min_angle = angle;

				/* if the point of choice was the point
					we started with, we're done */
				done = (i == num);

			}

			/* also done if included all points are included */
			done = done || (hull_points == (num - 1));
		}

	} while (! done); /* until get back to start */

	return hull_points+1;
}	



/*
Doesn't return the angle between the horizontal, but does
return a number with the same order properties
*/
static double theta(Point3d p1, Point3d p2)
{
	double dx,dy,ax,ay,t;

	
	dx = p2[0] - p1[0];
	if (dx < 0.0)
		ax = -dx;
	else
		ax = dx;

	dy = p2[1] - p1[1];
	if (dy < 0.0)
		ay = -dy;
	else
		ay = dy;


	if ((dx == 0.0) && (dy == 0.0))
		t = 0.0;
	else
		t = dy/(ax + ay);

	if (dx < 0)
		t = 2.0 - t;
	else if (dy < 0)
		t = 4 + t;

	return t*90.0;
}



	

		
	



--
Howard Look   Silicon Graphics   howardl@sgi.com   (415) 335-1780
.__   One of these :) after being run over by one of these O-O