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