/****************************************************************************** * Cagd_cci.c - curve curve intersection routines. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Nov. 96. * ******************************************************************************/ #include "cagd_loc.h" #include "geomat3d.h" #define X_COORD(i) (WPoints == NULL ? XPoints[i] : XPoints[i] / WPoints[i]) #define Y_COORD(i) (WPoints == NULL ? YPoints[i] : YPoints[i] / WPoints[i]) #define ZERO_APX_EQ(x, y) (ABS(x - y) < Eps * 10) static CagdPtStruct *GlblInterList = NULL; static CagdCrvStruct *GlblTanCrv1 = NULL, *GlblTanCrv2 = NULL; static void CagdCrvCrvInterAux(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2, CagdRType Eps); static CagdBType CagdCrvCrvInterNumer(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2, CagdRType Eps); static void InsertInterPoints(CagdRType t1, CagdRType t2, CagdRType Eps); /***************************************************************************** * DESCRIPTION: M * Given a curve, computes the angular span of its tangent field, in XY M * plane. M * * * PARAMETERS: M * Crv: Curve to consider M * Dir: General, median, direction of tangent field, in XY plane. M * AngularSpan: Maximal deviation of tangent field from Dir, in radians. M * * * RETURN VALUE: M * CagdBType: TRUE if angular span of curve is less than 180 degrees, M * FALSE otherwise. In the later case, Dir and Angle are invalid. M * * * SEE ALSO: M * CagdCrvCrvInter M * * * KEYWORDS: M * CagdCrvTanAngularSpan M *****************************************************************************/ CagdBType CagdCrvTanAngularSpan(CagdCrvStruct *Crv, CagdVType Dir, CagdRType *AngularSpan) { int i, Length = Crv -> Length; CagdRType R, *WPoints = Crv -> Points[0], *XPoints = Crv -> Points[1], *YPoints = Crv -> Points[2]; CagdVType V; *AngularSpan = 1.0; Dir[0] = X_COORD(Length - 1) - X_COORD(0); Dir[1] = Y_COORD(Length - 1) - Y_COORD(0); Dir[2] = 0.0; if ((R = sqrt(SQR(Dir[0]) + SQR(Dir[1]))) < IRIT_EPS) return FALSE; R = 1.0 / R; PT_SCALE(Dir, R); for (i = 1; i < Length; i++) { V[0] = X_COORD(i) - X_COORD(i - 1); V[1] = Y_COORD(i) - Y_COORD(i - 1); V[2] = 0.0; if ((R = sqrt(SQR(V[0]) + SQR(V[1]))) < IRIT_EPS) continue; R = 1.0 / R; PT_SCALE(V, R); if ((R = DOT_PROD(V, Dir)) < 0.0) return FALSE; if (*AngularSpan > R) *AngularSpan = R; } *AngularSpan = BOUND(*AngularSpan, 0.0, 1.0); *AngularSpan = acos(*AngularSpan); return TRUE; } /***************************************************************************** * DESCRIPTION: M * Computes the intersection points, if any, of the two given curves. M * * * PARAMETERS: M * Crv1, Crv2: Two curves to compute their intersection points. M * Eps: Accuracy of computation. M * * * RETURN VALUE: M * CagdPtStruct *: List of intersection points. Each points would M * contain (u1, u, 0.0). M * * * SEE ALSO: M * CagdCrvTanAngularSpan M * * * KEYWORDS: M * CagdCrvCrvInter M *****************************************************************************/ CagdPtStruct *CagdCrvCrvInter(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2, CagdRType Eps) { CagdBType IsBezier1 = FALSE, IsBezier2 = FALSE; GlblInterList = NULL; if (CAGD_IS_BEZIER_CRV(Crv1)) { IsBezier1 = TRUE; Crv1 = CnvrtBezier2BsplineCrv(Crv1); } if (CAGD_IS_BEZIER_CRV(Crv2)) { IsBezier2 = TRUE; Crv2 = CnvrtBezier2BsplineCrv(Crv2); } GlblTanCrv1 = CagdCrvDerive(Crv1); GlblTanCrv2 = CagdCrvDerive(Crv2); CagdCrvCrvInterAux(Crv1, Crv2, Eps); CagdCrvFree(GlblTanCrv1); CagdCrvFree(GlblTanCrv2); if (IsBezier1) CagdCrvFree(Crv1); if (IsBezier2) CagdCrvFree(Crv2); return GlblInterList; } /***************************************************************************** * DESCRIPTION: * * Auxliary function of CagdCrvCrvInter. * * * * PARAMETERS: * * Crv1, Crv2: Two curves to compute their intersection points. * * Eps: Accuracy of computation. * * * * RETURN VALUE: * * CagdPtStruct *: List of intersection points. Each points would * * contain (u1, u, 0.0). * *****************************************************************************/ static void CagdCrvCrvInterAux(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2, CagdRType Eps) { CagdBBoxStruct BBox1, BBox2; CagdVType Dir1, Dir2; CagdRType Angle1, Angle2, TMin1, TMax1, TMin2, TMax2; CagdCrvStruct *Crv1a, *Crv1b, *Crv2a, *Crv2b; CagdCrvBBox(Crv1, &BBox1); CagdCrvBBox(Crv2, &BBox2); if (BBox1.Max[0] < BBox2.Min[0] || BBox1.Max[1] < BBox2.Min[1] || BBox2.Max[0] < BBox1.Min[0] || BBox2.Max[1] < BBox1.Min[1]) { /* The bounding boxes do not overlap - return. */ return; } if (CagdCrvTanAngularSpan(Crv1, Dir1, &Angle1) && CagdCrvTanAngularSpan(Crv2, Dir2, &Angle2) && acos(DOT_PROD(Dir1, Dir2)) > Angle1 + Angle2) { /* The tangent field's cones do not overlap - only one intersection */ /* occur between the two curves - try to find it numerically. */ if (CagdCrvCrvInterNumer(Crv1, Crv2, Eps)) return; } /* Subdivide the two curves and recurse. */ CagdCrvDomain(Crv1, &TMin1, &TMax1); CagdCrvDomain(Crv2, &TMin2, &TMax2); if (TMax1 - TMin1 < IRIT_EPS ||TMax2 - TMin2 < IRIT_EPS) { /* Failed in the numerical approach - stop the subdivision! */ InsertInterPoints((TMin1 + TMax1) / 2.0, (TMin2 + TMax2) / 2.0, Eps); return; } Crv1a = CagdCrvSubdivAtParam(Crv1, (TMin1 + TMax1) / 2.0); Crv1b = Crv1a -> Pnext; Crv1a -> Pnext = NULL; Crv2a = CagdCrvSubdivAtParam(Crv2, (TMin2 + TMax2) / 2.0); Crv2b = Crv2a -> Pnext; Crv2a -> Pnext = NULL; CagdCrvCrvInterAux(Crv1a, Crv2a, Eps); CagdCrvCrvInterAux(Crv1a, Crv2b, Eps); CagdCrvCrvInterAux(Crv1b, Crv2a, Eps); CagdCrvCrvInterAux(Crv1b, Crv2b, Eps); CagdCrvFree(Crv1a); CagdCrvFree(Crv1b); CagdCrvFree(Crv2a); CagdCrvFree(Crv2b); } /***************************************************************************** * DESCRIPTION: * * Uses Newton raphson to try and converge on an intersection points of * * two curves. The two curves are assumed to have at most one intersection. * * * * PARAMETERS: * * Crv1, Crv2: Two curves to compute their intersection points. * * Eps: Accuracy of computation. * * * * RETURN VALUE: * * CagdBType: TRUE if found intersection, FALSE otherwise. * *****************************************************************************/ static CagdBType CagdCrvCrvInterNumer(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2, CagdRType Eps) { CagdRType t1, t2, TMin1, TMax1, TMin2, TMax2, CrntDist = IRIT_INFNTY; CagdPType Pt1, Pt2; CagdVType Tan1, Tan2; CagdCrvDomain(Crv1, &TMin1, &TMax1); CagdCrvDomain(Crv2, &TMin2, &TMax2); t1 = (TMin1 + TMax1) / 2.0; t2 = (TMin2 + TMax2) / 2.0; Pt1[2] = Pt2[2] = Tan1[2] = Tan2[2] = 0.0; while (TRUE) { CagdRType Dist, *R, Inter1Param, Inter2Param; CagdPType Inter1, Inter2; R = CagdCrvEval(Crv1, t1); CagdCoerceToE2(Pt1, &R, -1, Crv1 -> PType); R = CagdCrvEval(GlblTanCrv1, t1); CagdCoerceToE2(Tan1, &R, -1, GlblTanCrv1 -> PType); R = CagdCrvEval(Crv2, t2); CagdCoerceToE2(Pt2, &R, -1, Crv2 -> PType); R = CagdCrvEval(GlblTanCrv2, t2); CagdCoerceToE2(Tan2, &R, -1, GlblTanCrv2 -> PType); if ((Dist = PT_PT_DIST(Pt1, Pt2)) < Eps) { /* Done - found the intersection points. */ InsertInterPoints(t1, t2, Eps); return TRUE; } else if (Dist * 1.1 > CrntDist) {/* Failed to significantly improve. */ return FALSE; } else { CrntDist = Dist; } if (!CG2PointsFromLineLine(Pt1, Tan1, Pt2, Tan2, Inter1, &Inter1Param, Inter2, &Inter2Param)) return FALSE; t1 += Inter1Param; t1 = BOUND(t1, TMin1, TMax1); t2 += Inter2Param; t2 = BOUND(t2, TMin2, TMax2); } } /***************************************************************************** * DESCRIPTION: * * Insert t1/t2 values into GlblInterList, provided no equal t1 value exists * * already in the list. List is in ascending order with respect to t1. * * * * PARAMETERS: * * t1, t2: New parameter values to insert to global GlblInterList list. * * Eps: Accuracy of computation. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void InsertInterPoints(CagdRType t1, CagdRType t2, CagdRType Eps) { CagdPtStruct *PtTmp, *PtLast, *Pt; Pt = CagdPtNew(); Pt -> Pt[0] = t1; Pt -> Pt[1] = t2; Pt -> Pt[2] = 0.0; if (GlblInterList) { for (PtTmp = GlblInterList, PtLast = NULL; PtTmp != NULL; PtLast = PtTmp, PtTmp = PtTmp -> Pnext) { if (ZERO_APX_EQ(PtTmp -> Pt[0], t1)) { IritFree(Pt); return; } if (PtTmp -> Pt[0] > t1) break; } if (PtTmp) { /* Insert the new point in the middle of the list. */ Pt -> Pnext = PtTmp; if (PtLast) PtLast -> Pnext = Pt; else GlblInterList = Pt; } else { /* Insert the new point as the last point in the list. */ PtLast -> Pnext = Pt; } } else GlblInterList = Pt; }