/****************************************************************************** * Crv_Tans.c - computation of curve and point tangents and relations. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, January 96 * ******************************************************************************/ #include "symb_loc.h" #include "user_lib.h" #include "bool_lib.h" #include "iritprsr.h" #include "allocate.h" #define TAN_TAN_DIAGONAL_EPS 3e-3 /***************************************************************************** * DESCRIPTION: M * Computes the points on a C1 freeform planar Bspline curve, Crv, that a M * line tangent to Crv there goes through point Pt. That is, M * M * (C(t) - P) || C'(t), V * M * where || denotes a parallel constraint. M * * * PARAMETERS: M * Crv: To compute its tangent lines through Pt. M * Pt: Point of origin, all tangents to Crv goes through. M * Tolerance: Accuracy of computation. M * * * RETURN VALUE: M * CagdPtStruct *: A list of parameter location on Crv with tangent lines M * through Pt. Parameters are save in the X coordinate. M * * * SEE ALSO: M * SymbCrvCnvxHull, SymbTangentToCrvAtTwoPts, SymbCrvDiameter M * * * KEYWORDS: M * SymbCrvPtTangents M *****************************************************************************/ CagdPtStruct *SymbCrvPtTangents(CagdCrvStruct *Crv, CagdPType Pt, RealType Tolerance) { CagdPtStruct *RetList; CagdCrvStruct *DCrv, *TCrv, *Crv1, *Crv2, *DCrvW, *DCrvX, *DCrvY, *DCrvZ, *CrvW, *CrvX, *CrvY, *CrvZ; CagdRType Trans[3]; /* Make sure the given curve is open end conditioned curve. */ if (CAGD_IS_BEZIER_CRV(Crv)) Crv = CnvrtBezier2BsplineCrv(Crv); else Crv = CagdCrvCopy(Crv); if (CAGD_IS_PERIODIC_CRV(Crv)) { TCrv = CnvrtPeriodic2FloatCrv(Crv); CagdCrvFree(Crv); Crv = TCrv; } if (CAGD_IS_BSPLINE_CRV(Crv) && !BspCrvHasOpenEC(Crv)) { TCrv = BspCrvOpenEnd(Crv); CagdCrvFree(Crv); Crv = TCrv; } DCrv = CagdCrvDerive(Crv); /* Compute 'C(t) - Pt' into 'C(t)': */ PT_COPY(Trans, Pt); PT_SCALE(Trans, -1.0); CagdCrvTransform(Crv, Trans, 1.0); SymbCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ); CagdCrvFree(Crv); if (CrvW != NULL) CagdCrvFree(CrvW); if (CrvZ != NULL) CagdCrvFree(CrvZ); SymbCrvSplitScalar(DCrv, &DCrvW, &DCrvX, &DCrvY, &DCrvZ); CagdCrvFree(DCrv); if (DCrvW != NULL) CagdCrvFree(DCrvW); if (DCrvZ != NULL) CagdCrvFree(DCrvZ); /* Make sure "C(t) - Pt" is parallel to "C'(t)" by making sure that */ /* "C(t) - Pt" is orthogonal to a vertical to "C'(t)": */ Crv1 = SymbCrvMult(CrvX, DCrvY); CagdCrvFree(CrvX); CagdCrvFree(DCrvY); Crv2 = SymbCrvMult(CrvY, DCrvX); CagdCrvFree(CrvY); CagdCrvFree(DCrvX); Crv = SymbCrvSub(Crv1, Crv2); CagdCrvFree(Crv1); CagdCrvFree(Crv2); #ifdef DUMP_PT_CRV_TAN_SYMBOLIC_FUNC IritPrsrStdoutObject(GenCRVObject(Crv)); /* Also a memory leak... */ #endif /* DUMP_PT_CRV_TAN_SYMBOLIC_FUNC */ RetList = SymbCrvZeroSet(Crv, 1, Tolerance); CagdCrvFree(Crv); return RetList; } /***************************************************************************** * DESCRIPTION: M * Computes all the lines that are tangent to Crv at two locations. M * Returned is a list of parameter locations' pairs where the tangent is M * tangent to the curve. M * The tangents are computed as two sets of contours of the solution to M * the two equations of: M * 1. [ C(t) - C(r) ] || C'(t) V * 2. [ C(t) - C(r) ] || C'(r) V * and computing all the intersection points between these two sets of M * contours. M * Note that since equations 1 and 2 are symmetric, one only needs to solve M * for once and flip the notation of r and t. M * * * PARAMETERS: M * Crv: The curve to compute all tangent lines at two locations. M * FineNess: Of numeric search for the zero set (for surface subdivision). M * A positive value (10 is a good start). M * * * RETURN VALUE: M * CagdPtStruct *: A list of parameter location on Crv with tangent lines M * through Pt. Parameters are save in the X & Y coordinate. M * * * SEE ALSO: M * SymbCrvPtTangents, SymbCrvCnvxHull, SymbCrvDiameter M * * * KEYWORDS: M * SymbTangentToCrvAtTwoPts M *****************************************************************************/ CagdPtStruct *SymbTangentToCrvAtTwoPts(CagdCrvStruct *Crv, CagdRType FineNess) { int OldCirc; CagdRType TMin, TMax; static PlaneType Plane = { 1.0, 0.0, 0.0, 0.0 }; /* it is a scalar surface - only X. */ CagdCrvStruct *DCrv; CagdSrfStruct *Srf1, *Srf2, *Srf, *DSrf, *SrfT, *SrfR, *SrfW, *SrfX, *SrfY, *SrfZ, *DSrfW, *DSrfX, *DSrfY, *DSrfZ; IPPolygonStruct *Cntrs1, *Cntrs2, *CntrA, *CntrB; CagdPtStruct *RetList = NULL; /* Make sure the given curve is open end conditioned curve. */ if (CAGD_IS_BEZIER_CRV(Crv)) Crv = CnvrtBezier2BsplineCrv(Crv); else Crv = CagdCrvCopy(Crv); if (CAGD_IS_PERIODIC_CRV(Crv)) { CagdCrvStruct *TCrv = CnvrtPeriodic2FloatCrv(Crv); CagdCrvFree(Crv); Crv = TCrv; } if (CAGD_IS_BSPLINE_CRV(Crv) && !BspCrvHasOpenEC(Crv)) { CagdCrvStruct *TCrv = BspCrvOpenEnd(Crv); CagdCrvFree(Crv); Crv = TCrv; } CagdCrvDomain(Crv, &TMin, &TMax); /* Make sure the domain is zero to one. */ BspKnotAffineTrans(Crv -> KnotVector, Crv -> Length + Crv -> Order, -TMin, 1.0 / (TMax - TMin)); DCrv = CagdCrvDerive(Crv); DSrf = CagdPromoteCrvToSrf(DCrv, CAGD_CONST_U_DIR); SrfT = CagdPromoteCrvToSrf(Crv, CAGD_CONST_U_DIR); SrfR = CagdPromoteCrvToSrf(Crv, CAGD_CONST_V_DIR); CagdCrvFree(Crv); CagdCrvFree(DCrv); Srf1 = SymbSrfSub(SrfR, SrfT); CagdSrfFree(SrfR); CagdSrfFree(SrfT); SymbSrfSplitScalar(Srf1, &SrfW, &SrfX, &SrfY, &SrfZ); CagdSrfFree(Srf1); if (SrfW != NULL) CagdSrfFree(SrfW); if (SrfZ != NULL) CagdSrfFree(SrfZ); SymbSrfSplitScalar(DSrf, &DSrfW, &DSrfX, &DSrfY, &DSrfZ); CagdSrfFree(DSrf); if (DSrfW != NULL) CagdSrfFree(DSrfW); if (DSrfZ != NULL) CagdSrfFree(DSrfZ); /* Compute "parallel" constraint as perpendicular cross product zero. */ Srf1 = SymbSrfMult(SrfX, DSrfY); CagdSrfFree(SrfX); CagdSrfFree(DSrfY); Srf2 = SymbSrfMult(SrfY, DSrfX); CagdSrfFree(SrfY); CagdSrfFree(DSrfX); Srf = SymbSrfSub(Srf1, Srf2); CagdSrfFree(Srf1); CagdSrfFree(Srf2); /* Compute the zero set of the equation as contours. */ OldCirc = IritPrsrSetPolyListCirc(TRUE); Cntrs1 = UserCntrSrfWithPlane(Srf, Plane, FineNess); IritPrsrSetPolyListCirc(OldCirc); /* Move the rt parametric domain from YZ to XY. */ for (CntrA = Cntrs1; CntrA != NULL; CntrA = CntrA -> Pnext) { IPVertexStruct *V; for (V = CntrA -> PVertex; V != NULL; V = V -> Pnext) { V -> Coord[0] = V -> Coord[1]; V -> Coord[1] = V -> Coord[2]; V -> Coord[2] = 0.0; } } #ifdef DUMP_2TAN_SRF IritPrsrStdoutObject(GenSRFObject(Srf)); /* Also a memory leak! */ IritPrsrStdoutObject(GenPOLYObject(Cntrs1)); #endif /* DUMP_2TAN_SRF */ /* Create the flipped r <-> t copy of the contours. */ Cntrs2 = CopyPolygonList(Cntrs1); for (CntrB = Cntrs2; CntrB != NULL; CntrB = CntrB -> Pnext) { IPVertexStruct *V; for (V = CntrB -> PVertex; V != NULL; V = V -> Pnext) SWAP(CagdRType, V -> Coord[0], V -> Coord[1]); } #ifdef DUMP_2TAN_CRV_CNTRS IritPrsrStdoutObject(GenPOLYObject(Cntrs1)); /* Also a memory leak! */ IritPrsrStdoutObject(GenPOLYObject(Cntrs2)); #endif /* DUMP_2TAN_CRV_CNTRS */ /* Compute all the intersection points of both sets of contours. */ for (CntrB = Cntrs2; CntrB != NULL; CntrB = CntrB -> Pnext) { for (CntrA = Cntrs1; CntrA != NULL; CntrA = CntrA -> Pnext) { Bool2DInterStruct *Inter, *Inters = Boolean2DComputeInters(CntrA, CntrB, FALSE); for (Inter = Inters; Inter != NULL; ) { Bool2DInterStruct *NextInter = Inter -> Pnext; /* Eliminate the intersections along the diagonal (t == r) */ /* and due to symmetry - select only one out of each pair. */ if (!APX_EQ_EPS(Inter -> InterPt[0], Inter -> InterPt[1], TAN_TAN_DIAGONAL_EPS) && Inter -> InterPt[0] < Inter -> InterPt[1]) { CagdPtStruct *NewInterPt = CagdPtNew(); #ifdef DUMP_2TAN_CRV_INTER_PTS printf("[OBJECT [COLOR 13] [width 0.04] NONE\n [POINT %f %f 0]\n]\n", Inter -> InterPt[0], Inter -> InterPt[1]); #endif /* DUMP_2TAN_CRV_INTER_PTS */ NewInterPt -> Pnext = RetList; RetList = NewInterPt; NewInterPt -> Pt[0] = TMin + Inter -> InterPt[0] * (TMax - TMin); NewInterPt -> Pt[1] = TMin + Inter -> InterPt[1] * (TMax - TMin); NewInterPt -> Pt[2] = 0.0; } IritFree(Inter); Inter = NextInter; } } } return RetList; }