/****************************************************************************** * CnvxHull.c - convex hull computations. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Dec 1987. * ******************************************************************************/ #include "irit_sm.h" #include "misc_lib.h" #include "geom_loc.h" #define GM_CH_EXTRA_STACK 100 /* Spare number of points. */ STATIC_DATA RealType Ymin, Xmin; /* Global min Y point used by CompareAngle. */ STATIC_DATA int *LocalStack, GlblCHError = FALSE, StackPointer = 0; 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 * GMConvexHull M *****************************************************************************/ int GMConvexHull(GMR2PtStruct *DTPts, int *NumOfPoints) { int i, j, YminIndex, Index1, Index2, Index3, CHIndex; RealType p1x, p1y, p2x, p2y, p3x, p3y; GMR2PtStruct *CHPts; if (*NumOfPoints < 4) /* 3 or less points - must all be on convex hull. */ return TRUE; LocalStack = IritMalloc((GM_CH_EXTRA_STACK + *NumOfPoints) * sizeof(int)); 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 one (with Ymin) according */ /* to angle (0..180) relative to Ymin point: */ qsort(&DTPts[1], (*NumOfPoints) - 1, sizeof(GMR2PtStruct), 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 = (GMR2PtStruct *) /* Allocate mem. for CH. */ IritMalloc(sizeof(GMR2PtStruct) * *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 if (CHIndex == 0) { /* First location - cannot pop! */ 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) { IritFree(LocalStack); return FALSE; } } GEN_COPY(DTPts, CHPts, sizeof(GMR2PtStruct) * CHIndex); *NumOfPoints = CHIndex; IritFree(CHPts); IritFree(LocalStack); return TRUE; } /***************************************************************************** * DESCRIPTION: * * Routine to determine if given 3 points form convex angle using cross * * prod. * * * * PARAMETERS: * * p1x, p1y: First point. * * p2x, p2y: Second point. * * p3x, p3y: Third point. * * * * RETURN VALUE: * * int: TRUE if convex, FALSE otherwise. * *****************************************************************************/ static int IsConvex(RealType p1x, RealType p1y, RealType p2x, RealType p2y, RealType p3x, RealType p3y) { return (p1x - p2x) * (p2y - p3y) - (p2x - p3x) * (p1y - p2y) <= SQR(IRIT_UEPS); } /***************************************************************************** * DESCRIPTION: * * Resetsthe local stack so ConvexHull routine might use it. * * * * PARAMETERS: * * None * * * * RETURN VALUE: * * void * *****************************************************************************/ static void ResetLocalStack(void) { StackPointer = 0; } /***************************************************************************** * DESCRIPTION: * * Push one integer on local stack, assuming there is no overflow. * * * * PARAMETERS: * * i: To push. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void PushIndex(int i) { LocalStack[StackPointer++] = i; } /***************************************************************************** * DESCRIPTION: * * Pop one integer from local stack, assuming the stack never empty. * * * * PARAMETERS: * * None * * * * RETURN VALUE: * * int: Poped value * *****************************************************************************/ static int PopIndex(void) { if (StackPointer <= 0) { GEOM_FATAL_ERROR(GEOM_ERR_CH_STACK_UNDERFLOW); GlblCHError = TRUE; return 0; } return LocalStack[--StackPointer]; } /***************************************************************************** * DESCRIPTION: * * Compares two 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. * * * * PARAMETERS: * * Pt1, Pt2: Two points to compare their angles. * * * * RETURN VALUE: * * int: -1, 0 or 1 based on the angle ratios of the two points. * *****************************************************************************/ #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 ((FABS(Point1[0] - Xmin) < IRIT_EPS) && (FABS(Point1[1] - Ymin) < IRIT_EPS)) return -1; /* Make Min first */ if ((FABS(Point2[0] - Xmin) < IRIT_EPS) && (FABS(Point2[1] - Ymin) < IRIT_EPS)) return 1; /* Make Min first */ Dx = Point1[0] - Xmin; Value1 = SQR(Dx) * (-SIGN(Dx)) / (SQR(Dx) + SQR(Point1[1] - Ymin)); Dx = Point2[0] - Xmin; Value2 = SQR(Dx) * (-SIGN(Dx)) / (SQR(Dx) + SQR(Point2[1] - Ymin)); if (Value1 > Value2) return -1; else if (Value1 < Value2) return 1; else return 0; }