/***************************************************************************** * Routine to implement convex hull creation out of given set of 2D points. * * Algorithm is based on two articels: * * 1. An Efficient Algorithm For Determining The convex Hull of a Finite Set. * * By R.L. Graham, Information processing letters (1972) 132-133. * * 2. A Reevolution of an Efficient Algorithm For Determining The Convex Hull * * of a Finite Planar Set. By K.R. Anderson, Information Processing * * Letters, January 1978, Vol. 7, Num. 1, 53-55. * * * * Written by: Gershon Elber Ver 0.1, Dec 1987 * *****************************************************************************/ #include "irit_sm.h" #include "imalloc.h" #include "extra_fn.h" #define MAX_STACK 1000 /* Max. number of points */ #define XAXIS 0 #define YAXIS 1 static RealType Ymin, Xmin;/* Global minimum Y point - used by CompareAngle. */ static int GlblCHError; static int IsConvex(RealType p1x, RealType p1y, RealType p2x, RealType p2y, RealType p3x, RealType p3y); static void ResetLocalStack(void); static void PushIndex(int i); static int PopIndex(void); #if defined(ultrix) && defined(mips) static int CompareAngle(VoidPtr Pt1, VoidPtr Pt2); #else static int CompareAngle(const VoidPtr Pt1, const VoidPtr Pt2); #endif /* ultrix && mips (no const support) */ /***************************************************************************** * DESCRIPTION: M * Convex Hull computation of a set of points. The Convex Hull is returned M * in place, updating NumOfPoints. M * Algorithm is based on two articles: M * 1. An Efficient Algorithm For Determining The convex Hull of a Finite Set. V * By R.L. Graham, Information processing letters (1972) 132-133. V * 2. A Reevolution of an Efficient Algorithm For Determining The Convex Hull V * of a Finite Planar Set. By K.R. Anderson, Information Processing V * Letters, January 1978, Vol. 7, Num. 1, 53-55. V * * * PARAMETERS: M * DTPts: The set of point to compute their convex hull. M * NumOfPoints: Number of points in set DTPts. M * * * RETURN VALUE: M * int: TRUE if successful, FALSE otherwise. M * * * KEYWORDS: M * ConvexHull M *****************************************************************************/ int ConvexHull(ConvexHullPoint *DTPts, int *NumOfPoints) { int i, j, YminIndex, Index1, Index2, Index3, CHIndex; RealType p1x, p1y, p2x, p2y, p3x, p3y; ConvexHullPoint *CHPts; Xmin = DTPts[0].Pt[0]; Ymin = DTPts[0].Pt[1]; /* Find Y minimum: */ YminIndex = 0; for (i = 1; i < *NumOfPoints; i++) if (Ymin >= DTPts[i].Pt[1]) { /* If Y levels are the same pick the one leftmost on X axis: */ if (Ymin == DTPts[i].Pt[1] && Xmin < DTPts[i].Pt[0]) continue; Xmin = DTPts[i].Pt[0]; Ymin = DTPts[i].Pt[1]; YminIndex = i; } /* Make sure first point in array is at Y minimum. */ SWAP(RealType, DTPts[0].Pt[0], DTPts[YminIndex].Pt[0]); SWAP(RealType, DTPts[0].Pt[1], DTPts[YminIndex].Pt[1]); /* Sort all the other elements but the first obe (with Ymin) according */ /* to angle (0..180) relative to Ymin point: */ qsort(&DTPts[1], (*NumOfPoints) - 1, sizeof(ConvexHullPoint), CompareAngle); /* Eliminate duplicated points. */ for (i = j = 1; i < *NumOfPoints; ) { if (APX_EQ(DTPts[i].Pt[0], DTPts[i - 1].Pt[0]) && APX_EQ(DTPts[i].Pt[1], DTPts[i - 1].Pt[1])) { /* Same point - do not copy. */ i++; } else { DTPts[j].Pt[0] = DTPts[i].Pt[0]; DTPts[j].Pt[1] = DTPts[i].Pt[1]; i++; j++; } } *NumOfPoints = j; CHPts = (ConvexHullPoint *) /* Allocate mem. for CH. */ IritMalloc(sizeof(ConvexHullPoint) * *NumOfPoints); /* Now the points are sorted by angle, so we can scan them and drop non */ /* convex ones (by using 3 "running" points and testing middles angle). */ p1x = DTPts[0].Pt[0]; p1y = DTPts[0].Pt[1]; p2x = DTPts[1].Pt[0]; p2y = DTPts[1].Pt[1]; p3x = DTPts[2].Pt[0]; p3y = DTPts[2].Pt[1]; Index1 = 0; /* Indices of points.*/ Index2 = 1; Index3 = 2; CHIndex = 0; ResetLocalStack(); GlblCHError = FALSE; while (Index1 < *NumOfPoints) { if (IsConvex(p1x, p1y, p2x, p2y, p3x, p3y)) { /* Go forward */ CHPts[CHIndex].Pt[0] = p1x; CHPts[CHIndex++].Pt[1] = p1y; PushIndex(Index1); /* Save that index, we might need that ... */ p1x = p2x; p1y = p2y; Index1 = Index2; p2x = p3x; p2y = p3y; Index2 = Index3++; p3x = DTPts[Index3 % *NumOfPoints].Pt[0]; p3y = DTPts[Index3 % *NumOfPoints].Pt[1]; } else { /* Go backward until convex corner, using last pushed pos. */ p2x = p1x; p2y = p1y; Index2 = Index1; Index1 = PopIndex(); CHIndex--; p1x = DTPts[Index1].Pt[0]; p1y = DTPts[Index1].Pt[1]; } if (GlblCHError) return FALSE; } GEN_COPY(DTPts, CHPts, sizeof(ConvexHullPoint) * CHIndex); *NumOfPoints = CHIndex; IritFree((VoidPtr) CHPts); return TRUE; } /***************************************************************************** * DESCRIPTION: * * Routine to determine if given 3 points form convex angle using cross * * prod. * *****************************************************************************/ static int IsConvex(RealType p1x, RealType p1y, RealType p2x, RealType p2y, RealType p3x, RealType p3y) { return (p1x - p2x) * (p2y - p3y) - (p2x - p3x) * (p1y - p2y) < 0; } int LocalStack[MAX_STACK], StackPointer; /***************************************************************************** * DESCRIPTION: * * Routine to reset the local stack so ConvexHull routine might use it. * *****************************************************************************/ static void ResetLocalStack(void) { StackPointer = 0; } /***************************************************************************** * DESCRIPTION: * * Routine to push one integer on local stack, assuming there is no * * overflow. * *****************************************************************************/ static void PushIndex(int i) { if (StackPointer >= MAX_STACK) { fprintf(stderr, "Max stack for convex hull computation exceeded.\n"); GlblCHError = TRUE; return; } LocalStack[StackPointer++] = i; } /***************************************************************************** * DESCRIPTION: * * Routine to pop one integer from local stack, assuming the stack never * * empty. * *****************************************************************************/ static int PopIndex(void) { if (StackPointer <= 0) { fprintf(stderr, "underflow in convex hull stack.\n"); GlblCHError = TRUE; return 0; } return LocalStack[--StackPointer]; } /***************************************************************************** * DESCRIPTION: * * Routine to Compare to angles of two points relative to Xmin,Ymin point. * * This routine is used by the qsort routine in ConvexHull routine. * * Algorithm: -SIGN(cos(teta)) * cos(teta)^2, see second reference, page 54. * *****************************************************************************/ #if defined(ultrix) && defined(mips) static int CompareAngle(VoidPtr Pt1, VoidPtr Pt2) #else static int CompareAngle(const VoidPtr Pt1, const VoidPtr Pt2) #endif /* ultrix && mips (no const support) */ { RealType Dx, Value1, Value2, *Point1 = (RealType *) Pt1, *Point2 = (RealType *) Pt2; if ((ABS(Point1[XAXIS] - Xmin) < IRIT_EPS) && (ABS(Point1[YAXIS] - Ymin) < IRIT_EPS)) return -1; /* Make Min first */ if ((ABS(Point2[XAXIS] - Xmin) < IRIT_EPS) && (ABS(Point2[YAXIS] - Ymin) < IRIT_EPS)) return 1; /* Make Min first */ Dx = Point1[XAXIS] - Xmin; Value1 = SQR(Dx) * (-SIGN(Dx)) / (SQR(Dx) + SQR(Point1[YAXIS] - Ymin)); Dx = Point2[XAXIS] - Xmin; Value2 = SQR(Dx) * (-SIGN(Dx)) / (SQR(Dx) + SQR(Point2[YAXIS] - Ymin)); if (Value1 > Value2) return -1; else if (Value1 < Value2) return 1; else return 0; }