/****************************************************************************** * MS_Circ.c - minimum spanning circle/cone. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, May 2002. * ******************************************************************************/ #include "irit_sm.h" #include "misc_lib.h" #include "geom_loc.h" static void GMMinSpanCircWithPt(GMR2PtStruct *DTPts, int NumOfPoints, GMR2PtStruct *Q, GMR2PtStruct *Center, RealType *Radius); static void GMMinSpanCircWith2Pts(GMR2PtStruct *DTPts, int NumOfPoints, GMR2PtStruct *Q1, GMR2PtStruct *Q2, GMR2PtStruct *Center, RealType *Radius); static void GMMinSpanConeWithPt(VectorType *DTPts, int NumOfPoints, VectorType Q, VectorType Center, RealType *Radius); static void GMMinSpanConeWith2Pts(VectorType *DTPts, int NumOfPoints, VectorType Q1, VectorType Q2, VectorType Center, RealType *Radius); /***************************************************************************** * DESCRIPTION: M * Minimum spanning circle (MSC) computation of a set of points. M * Algorithm is based on Section 4.7 of "Computational Geometry, Algorithms M * and Applications" by M. de Berg et. al. M * * * PARAMETERS: M * DTPts: The set of point to compute their MSC. M * NumOfPoints: Number of points in set DTPts. M * Center: Of computed MSC. M * Radius: Of computed MSC. M * * * RETURN VALUE: M * int: TRUE if successful, FALSE otherwise. M * * * SEE ALSO: M * GMMinSpanCone M * * * KEYWORDS: M * GMMinSpanCirc, minimum spanning circle M *****************************************************************************/ int GMMinSpanCirc(GMR2PtStruct *DTPts, int NumOfPoints, GMR2PtStruct *Center, RealType *Radius) { int i; RealType RadiusSqr; if (NumOfPoints < 2) { *Center = *DTPts; *Radius = 0.0; return TRUE; } /* Compute a trivial bound to first two point. */ Center -> Pt[0] = (DTPts[0].Pt[0] + DTPts[1].Pt[0]) * 0.5; Center -> Pt[1] = (DTPts[0].Pt[1] + DTPts[1].Pt[1]) * 0.5; RadiusSqr = (SQR(DTPts[0].Pt[0] - DTPts[1].Pt[0]) + SQR(DTPts[0].Pt[1] - DTPts[1].Pt[1])) * 0.25; /* And examine the rest if inside. */ for (i = 2; i < NumOfPoints; i++) { RealType RadSqr = SQR(DTPts[i].Pt[0] - Center -> Pt[0]) + SQR(DTPts[i].Pt[1] - Center -> Pt[1]); if (RadSqr > RadiusSqr) GMMinSpanCircWithPt(DTPts, i, &DTPts[i], Center, &RadiusSqr); } *Radius = sqrt(RadiusSqr); return TRUE; } /***************************************************************************** * DESCRIPTION: * * Auxiliary function of GMMinSpanCirc. * * * * PARAMETERS: * * DTPts: The set of point to compute their MSC. * * NumOfPoints: Number of points in set DTPts. * * Q: Extra point that must be on the boundary of the MSC. * * Center: Of computed MSC. * * RadiusSqr: Of computed MSC. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void GMMinSpanCircWithPt(GMR2PtStruct *DTPts, int NumOfPoints, GMR2PtStruct *Q, GMR2PtStruct *Center, RealType *RadiusSqr) { int i; if (NumOfPoints < 1) { GEOM_FATAL_ERROR(GEOM_ERR_MSC_TOO_FEW_PTS); return; } /* Compute a trivial bound to first two point. */ Center -> Pt[0] = (DTPts[0].Pt[0] + Q -> Pt[0]) * 0.5; Center -> Pt[1] = (DTPts[0].Pt[1] + Q -> Pt[1]) * 0.5; *RadiusSqr = (SQR(DTPts[0].Pt[0] - Q -> Pt[0]) + SQR(DTPts[0].Pt[1] - Q -> Pt[1])) * 0.25; /* And examine the rest if inside. */ for (i = 1; i < NumOfPoints; i++) { RealType RadSqr = SQR(DTPts[i].Pt[0] - Center -> Pt[0]) + SQR(DTPts[i].Pt[1] - Center -> Pt[1]); if (RadSqr > *RadiusSqr) GMMinSpanCircWith2Pts(DTPts, i, &DTPts[i], Q, Center, RadiusSqr); } } /***************************************************************************** * DESCRIPTION: * * Auxiliary function of GMMinSpanCirc. * * * * PARAMETERS: * * DTPts: The set of point to compute their MSC. * * NumOfPoints: Number of points in set DTPts. * * Q1, Q2: Extra points that must be on the boundary of the MSC. * * Center: Of computed MSC. * * RadiusSqr: Of computed MSC. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void GMMinSpanCircWith2Pts(GMR2PtStruct *DTPts, int NumOfPoints, GMR2PtStruct *Q1, GMR2PtStruct *Q2, GMR2PtStruct *Center, RealType *RadiusSqr) { int i; /* Compute a trivial bound to first two point. */ Center -> Pt[0] = (Q1 -> Pt[0] + Q2 -> Pt[0]) * 0.5; Center -> Pt[1] = (Q1 -> Pt[1] + Q2 -> Pt[1]) * 0.5; *RadiusSqr = (SQR(Q2 -> Pt[0] - Q1 -> Pt[0]) + SQR(Q2 -> Pt[1] - Q1 -> Pt[1])) * 0.25; /* And examine the rest if inside. */ for (i = 0; i < NumOfPoints; i++) { RealType RadSqr = SQR(DTPts[i].Pt[0] - Center -> Pt[0]) + SQR(DTPts[i].Pt[1] - Center -> Pt[1]); if (RadSqr > *RadiusSqr) { PointType C, Pt1, Pt2, Pt3; /* Compute the circle through Q1, Q2, and DTPts[i]. */ PT2D_COPY(Pt1, Q1 -> Pt); PT2D_COPY(Pt2, Q2 -> Pt); PT2D_COPY(Pt3, DTPts[i].Pt); Pt1[2] = Pt2[2] = Pt3[2] = 0.0; if (!GMCircleFrom3Points(C, RadiusSqr, Pt1, Pt2, Pt3)) { /* Three points are collinear. */ GEOM_FATAL_ERROR(GEOM_ERR_MSC_COLIN_CIRC); return; } *RadiusSqr = SQR(*RadiusSqr); Center -> Pt[0] = C[0]; Center -> Pt[1] = C[1]; } } } /***************************************************************************** * DESCRIPTION: M * Minimum spanning cone (MSN) computation of a set of vectors. M * Find a central vector as the average of all given vectors and find the M * vector with maximal angular distance from it. M * * * PARAMETERS: M * DTVecs: The set of vectors to compute their MSN. M * VecsNormalized: TRUE if vectors are normalized, FALSE otherwise. M * NumOfVecs: Number of vectors in set DTVecs. M * Center: Of computed MSN. M * Angle: Of computed MSN, in radians. M * * * RETURN VALUE: M * int: TRUE if successful, FALSE otherwise. M * * * SEE ALSO: M * GMMinSpanCirc, GMMinSpanCone M * * * KEYWORDS: M * GMMinSpanConeAvg M *****************************************************************************/ int GMMinSpanConeAvg(VectorType *DTVecs, int VecsNormalized, int NumOfVecs, VectorType Center, RealType *Angle) { int i; RealType IProd, MinIProd = 1.0; VectorType *NrmlDTVecs; if (NumOfVecs < 2) { GEOM_FATAL_ERROR(GEOM_ERR_MSC_TOO_FEW_PTS); return FALSE; } if (!VecsNormalized) { NrmlDTVecs = (VectorType *) IritMalloc(NumOfVecs * sizeof(VectorType)); GEN_COPY(NrmlDTVecs, DTVecs, NumOfVecs * sizeof(VectorType)); } else NrmlDTVecs = DTVecs; /* Compute the center of the cone. */ VEC_RESET(Center); for (i = 0; i < NumOfVecs; i++) { if (!VecsNormalized) VEC_NORMALIZE(NrmlDTVecs[i]); VEC_ADD(Center, Center, NrmlDTVecs[i]); } VEC_NORMALIZE(Center); /* Compute the angular openning of the cone. */ for (i = 0; i < NumOfVecs; i++) { IProd = DOT_PROD(Center, NrmlDTVecs[i]); if (MinIProd > IProd) MinIProd = IProd; } *Angle = acos(MinIProd); if (!VecsNormalized) IritFree(NrmlDTVecs); return TRUE; } /***************************************************************************** * DESCRIPTION: M * Minimum spanning cone (MSN) computation of a set of vectors. M * Algorithm is based on the Minimum Spanning Circle in Section 4.7 of M * "Computational Geometry, Algorithms and Applications" by M. de Berg et. al.M * * * PARAMETERS: M * DTVecs: The set of vectors to compute their MSN. M * VecsNormalized: TRUE if vectors are normalized, FALSE otherwise. M * NumOfVecs: Number of vectors in set DTVecs. M * Center: Of computed MSN. M * Angle: Of computed MSN, in radians. M * * * RETURN VALUE: M * int: TRUE if successful, FALSE otherwise. M * * * SEE ALSO: M * GMMinSpanCirc, GMMinSpanConeAvg M * * * KEYWORDS: M * GMMinSpanCone M *****************************************************************************/ int GMMinSpanCone(VectorType *DTVecs, int VecsNormalized, int NumOfVecs, VectorType Center, RealType *Angle) { int i; VectorType *NrmlDTVecs; if (NumOfVecs < 2) { VEC_COPY(Center, DTVecs[0]); *Angle = 0.0; return TRUE; } if (!VecsNormalized) { NrmlDTVecs = (VectorType *) IritMalloc(NumOfVecs * sizeof(VectorType)); GEN_COPY(NrmlDTVecs, DTVecs, NumOfVecs * sizeof(VectorType)); for (i = 0; i < NumOfVecs; i++) { VEC_NORMALIZE(NrmlDTVecs[i]); } } else NrmlDTVecs = DTVecs; /* Compute a trivial bound to first two vectors. */ PT_BLEND(Center, NrmlDTVecs[0], NrmlDTVecs[1], 0.5); PT_NORMALIZE(Center); *Angle = acos(DOT_PROD(NrmlDTVecs[0], NrmlDTVecs[1])) * 0.5; /* And examine the rest if inside. */ for (i = 2; i < NumOfVecs; i++) { RealType Ang = acos(DOT_PROD(NrmlDTVecs[i], Center)); if (Ang > *Angle) GMMinSpanConeWithPt(NrmlDTVecs, i, NrmlDTVecs[i], Center, Angle); } if (!VecsNormalized) IritFree(NrmlDTVecs); return TRUE; } /***************************************************************************** * DESCRIPTION: * * Auxiliary function of GMMinSpanCone. * * * * PARAMETERS: * * DTVecs: The set of vector to compute their MSC. * * NumOfVecs: Number of vectors in set DTVecs. * * Q: Extra vector that must be on the boundary of the MSC. * * Center: Of computed MSC. * * Angle: Of computed MSC. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void GMMinSpanConeWithPt(VectorType *DTVecs, int NumOfVecs, VectorType Q, VectorType Center, RealType *Angle) { int i; if (NumOfVecs < 1) { GEOM_FATAL_ERROR(GEOM_ERR_MSC_TOO_FEW_PTS); return; } /* Compute a trivial bound to first two vector. */ PT_BLEND(Center, DTVecs[0], Q, 0.5); PT_NORMALIZE(Center); *Angle = acos(DOT_PROD(DTVecs[0], Q)) * 0.5; /* And examine the rest if inside. */ for (i = 1; i < NumOfVecs; i++) { RealType Ang = acos(DOT_PROD(DTVecs[i], Center)); if (Ang > *Angle) GMMinSpanConeWith2Pts(DTVecs, i, DTVecs[i], Q, Center, Angle); } } /***************************************************************************** * DESCRIPTION: * * Auxiliary function of GMMinSpanCone. * * * * PARAMETERS: * * DTVecs: The set of vector to compute their MSC. * * NumOfVecs: Number of vectors in set DTVecs. * * Q1, Q2: Extra vectors that must be on the boundary of the MSC. * * Center: Of computed MSC. * * Angle: Of computed MSC. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void GMMinSpanConeWith2Pts(VectorType *DTVecs, int NumOfVecs, VectorType Q1, VectorType Q2, VectorType Center, RealType *Angle) { int i; /* Compute a trivial bound to first two vector. */ PT_BLEND(Center, Q1, Q2, 0.5); PT_NORMALIZE(Center); *Angle = acos(DOT_PROD(Q1, Q2)) * 0.5; /* And examine the rest if inside. */ for (i = 0; i < NumOfVecs; i++) { RealType Ang = acos(DOT_PROD(DTVecs[i], Center)); if (Ang > *Angle) { RealType R; PlaneType Plane; /* Compute the cone through Q1, Q2, and DTVecs[i]. */ GMPlaneFrom3Points(Plane, Q1, Q2, DTVecs[i]); VEC_COPY(Center, Plane); R = DOT_PROD(Center, Q1); if (R < 0.0) { VEC_SCALE(Center, -1); R = -R; } *Angle = acos(R); } } }