/****************************************************************************** * SymbZero.c - computes the zeros and extremes of a given object. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, March 93. * ******************************************************************************/ #include "symb_loc.h" #define SET_DEFAULT_IRIT_EPS 1e-3 #define SET_DEFAULT_IRIT_EPS 1e-3 #define ZERO_APX_EQ(x, y) (ABS(x - y) < GlblSetEpsilon * 10) #define MID_RATIO 0.5123456789 static CagdPtStruct *GlblPtList = NULL; static CagdRType CrvTMin, CrvTMax, GlblSetEpsilon = SET_DEFAULT_IRIT_EPS; static void SymbScalarCrvConstSet(CagdCrvStruct *Crv, CagdRType ConstVal); static void SymbScalarCrvExtremSet(CagdCrvStruct *Crv); static void SymbInsertNewParam(CagdRType t); /***************************************************************************** * DESCRIPTION: M * Computes the zero set of a given curve, in the given axis (1-3 for X-Z). M * Returned is a list of the zero set points holding the parameter values M * att Pt[0] of each point. M * * * PARAMETERS: M * Crv: To compute its zero set. M * Axis: The axis of Crv to compute zero set for. M * Epsilon: Tolerance control. M * * * RETURN VALUE: M * CagdPtStruct *: List of parameter values form which Crv is zero in M * axis Axis. M * * * KEYWORDS: M * SymbCrvZeroSet, zero set, symbolic computation M *****************************************************************************/ CagdPtStruct *SymbCrvZeroSet(CagdCrvStruct *Crv, int Axis, CagdRType Epsilon) { return SymbCrvConstSet(Crv, Axis, Epsilon, 0.0); } /***************************************************************************** * DESCRIPTION: M * Computes the extremum set of a given curve, in given axis (1-3 for X-Z). M * Returned is a list of the extreme set points holding the parameter M * values at Pt[0] of each point. M * One could compute the derivative of the curve and find its zero set. M * However, for rational curves, this will double the degree and slow down M * the computation considerably. M * * * PARAMETERS: M * Crv: To compute its extremum set. M * Axis: The axis of Crv to compute extremum set for. M * Epsilon: Tolerance control. M * * * RETURN VALUE: M * CagdPtStruct *: List of parameter values form which Crv has an M * extremum value in axis Axis. M * * * KEYWORDS: M * SymbCrvExtremSet, extremum set, symbolic computation M *****************************************************************************/ CagdPtStruct *SymbCrvExtremSet(CagdCrvStruct *Crv, int Axis, CagdRType Epsilon) { CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *TCrv, *ScalarCrv = NULL; GlblSetEpsilon = Epsilon; SymbCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ); switch (Axis) { case 1: if (CrvX) ScalarCrv = CrvX; else SYMB_FATAL_ERROR(SYMB_ERR_OUT_OF_RANGE); break; case 2: if (CrvY) ScalarCrv = CrvY; else SYMB_FATAL_ERROR(SYMB_ERR_OUT_OF_RANGE); break; case 3: if (CrvZ) ScalarCrv = CrvZ; else SYMB_FATAL_ERROR(SYMB_ERR_OUT_OF_RANGE); break; default: SYMB_FATAL_ERROR(SYMB_ERR_OUT_OF_RANGE); } Crv = SymbCrvMergeScalar(CrvW, ScalarCrv, NULL, NULL); if (CrvW) CagdCrvFree(CrvW); if (CrvX) CagdCrvFree(CrvX); if (CrvY) CagdCrvFree(CrvY); if (CrvZ) CagdCrvFree(CrvZ); GlblPtList = NULL; if (CAGD_IS_BEZIER_CRV(Crv)) { /* We need to save domains. */ TCrv = CnvrtBezier2BsplineCrv(Crv); CagdCrvFree(Crv); Crv = TCrv; } CagdCrvDomain(Crv, &CrvTMin, &CrvTMax); SymbScalarCrvExtremSet(Crv); CagdCrvFree(Crv); return GlblPtList; } /***************************************************************************** * DESCRIPTION: M * Computes the constant set of a given curve, in the given axis M * (1-3 for X-Z). M * Returned is a list of the constant set points holding the parameter M * values at Pt[0] of each point. M * * * PARAMETERS: M * Crv: To compute its constant set. M * Axis: The axis of Crv to compute constant set for. M * Epsilon: Tolerance control. M * ConstVal: The value at which to compute the constant set. M * * * RETURN VALUE: M * CagdPtStruct *: List of parameter values form which Crv has an M * value of ConstVal in axis Axis. M * * * KEYWORDS: M * SymbCrvConstSet, constant set, zero set, symbolic computation M *****************************************************************************/ CagdPtStruct *SymbCrvConstSet(CagdCrvStruct *Crv, int Axis, CagdRType Epsilon, CagdRType ConstVal) { CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *TCrv, *ScalarCrv = NULL; GlblSetEpsilon = Epsilon; SymbCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ); switch (Axis) { case 1: if (CrvX) ScalarCrv = CrvX; else SYMB_FATAL_ERROR(SYMB_ERR_OUT_OF_RANGE); break; case 2: if (CrvY) ScalarCrv = CrvY; else SYMB_FATAL_ERROR(SYMB_ERR_OUT_OF_RANGE); break; case 3: if (CrvZ) ScalarCrv = CrvZ; else SYMB_FATAL_ERROR(SYMB_ERR_OUT_OF_RANGE); break; default: SYMB_FATAL_ERROR(SYMB_ERR_OUT_OF_RANGE); } Crv = SymbCrvMergeScalar(CrvW, ScalarCrv, NULL, NULL); if (CrvW) CagdCrvFree(CrvW); if (CrvX) CagdCrvFree(CrvX); if (CrvY) CagdCrvFree(CrvY); if (CrvZ) CagdCrvFree(CrvZ); GlblPtList = NULL; if (CAGD_IS_BEZIER_CRV(Crv)) { /* We need to save domains. */ TCrv = CnvrtBezier2BsplineCrv(Crv); CagdCrvFree(Crv); Crv = TCrv; } CagdCrvDomain(Crv, &CrvTMin, &CrvTMax); SymbScalarCrvConstSet(Crv, ConstVal); CagdCrvFree(Crv); return GlblPtList; } /***************************************************************************** * DESCRIPTION: * * Computes the zero set of a scalar curve. Curve might be rational. * * Assumes open end condition on the curve parametric domain. * * Zero set points are append to the GlblPtList point list. * * * * PARAMETERS: * * Crv: To compute its constant set. * * ConstVal: Thep value at which to compute the constant set. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void SymbScalarCrvConstSet(CagdCrvStruct *Crv, CagdRType ConstVal) { int i, Len = Crv -> Length, Sign = 0; CagdRType *ScalarPts = Crv -> Points[1], *ScalarWPts = Crv -> Points[0]; if (SymbCrvPosNegWeights(Crv)) { i = 0; /* Force subdivision of the curve. */ } else { for (i = 0; i < Len; i++) { CagdRType V = (ScalarWPts ? ScalarPts[i] / ScalarWPts[i] : ScalarPts[i]) - ConstVal; int NewSign = SIGN(V); if (Sign * NewSign < 0) break; else if (NewSign) Sign = NewSign; } } if (i < Len) { CagdRType TMin, TMax, TMid; /* Control polygon is both positive and negative. */ CagdCrvDomain(Crv, &TMin, &TMax); TMid = TMin * MID_RATIO + TMax * (1.0 - MID_RATIO); if (TMax - TMin < GlblSetEpsilon / 5.0) { /* Small enough - stop. */ SymbInsertNewParam((TMax + TMin) / 2.0); } else { /* Needs to subdiv. */ CagdCrvStruct *Crv1 = CagdCrvSubdivAtParam(Crv, TMid), *Crv2 = Crv1 -> Pnext; SymbScalarCrvConstSet(Crv1, ConstVal); SymbScalarCrvConstSet(Crv2, ConstVal); CagdCrvFree(Crv1); CagdCrvFree(Crv2); } } } /***************************************************************************** * DESCRIPTION: * * Computes the extrem set of a scalar curve. Curve might be rational. * * Assumes open end condition on the curve parametric domain. * * Extreme set points are append to the GlblPtList point list. * * * * PARAMETERS: * * Crv: To compute its extremum set. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void SymbScalarCrvExtremSet(CagdCrvStruct *Crv) { int i, Len = Crv -> Length, Sign = 0; CagdRType *ScalarPts = Crv -> Points[1], *ScalarWPts = Crv -> Points[0]; if (SymbCrvPosNegWeights(Crv)) { i = 0; /* Force subdivision of the curve. */ } else { CagdRType PrevV = ScalarWPts ? ScalarPts[0] / ScalarWPts[0] : ScalarPts[0]; for (i = 0; i < Len; i++) { CagdRType V = (ScalarWPts ? ScalarPts[i] / ScalarWPts[i] : ScalarPts[i]), DV = V - PrevV; int NewSign = SIGN(DV); if (Sign * NewSign < 0) break; else if (NewSign) Sign = NewSign; PrevV = V; } } if (i < Len) { CagdRType TMin, TMax, TMid; /* Control polygon is both increasing and decreasing. */ CagdCrvDomain(Crv, &TMin, &TMax); TMid = TMin * MID_RATIO + TMax * (1.0 - MID_RATIO); if (TMax - TMin < GlblSetEpsilon / 5.0) { /* Small enough - stop. */ SymbInsertNewParam((TMax + TMin) / 2.0); } else { /* Needs to subdiv. */ CagdCrvStruct *Crv1 = CagdCrvSubdivAtParam(Crv, TMid), *Crv2 = Crv1 -> Pnext; SymbScalarCrvExtremSet(Crv1); SymbScalarCrvExtremSet(Crv2); CagdCrvFree(Crv1); CagdCrvFree(Crv2); } } else { CagdRType V1, V2, TMin, TMax; CagdCrvDomain(Crv, &TMin, &TMax); /* This segment is monotone. Test if end points are extremes. */ V1 = ScalarWPts ? ScalarPts[0] / ScalarWPts[0] : ScalarPts[0]; V2 = ScalarWPts ? ScalarPts[1] / ScalarWPts[1] : ScalarPts[1]; if (APX_EQ(V1, V2)) SymbInsertNewParam(TMin); V1 = ScalarWPts ? ScalarPts[Len - 2] / ScalarWPts[Len - 2] : ScalarPts[Len - 2]; V2 = ScalarWPts ? ScalarPts[Len - 1] / ScalarWPts[Len - 1] : ScalarPts[Len - 1]; if (APX_EQ(V1, V2)) SymbInsertNewParam(TMax); } } /***************************************************************************** * DESCRIPTION: M * Returns TRUE iff the Crv is not rational or rational with weights that are M * entirely positive or entirely negative. M * * * PARAMETERS: M * Crv: To examine for same sign weights, if any. M * * * RETURN VALUE: M * CagdBType: TRUE if no weights or all of same sign. M * * * KEYWORDS: M * SymbCrvPosNegWeights, symbolic computation M *****************************************************************************/ CagdBType SymbCrvPosNegWeights(CagdCrvStruct *Crv) { int i; CagdBType HasNeg, HasPos; CagdRType *Weights = Crv -> Points[0]; if (Weights == NULL) return FALSE; /* Curve is not rational. */ for (HasNeg = HasPos = FALSE, i = Crv -> Length - 1; i >= 0; i--) { HasNeg |= *Weights < 0.0; HasPos |= *Weights++ > 0.0; } return HasNeg && HasPos; } /***************************************************************************** * DESCRIPTION: * * Insert a single t value into existing GlblPtList, provided no equal t * * value exists already in the list. List is ordered incrementally. * * * * PARAMETERS: * * t: New value to insert to global GlblPtList list. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void SymbInsertNewParam(CagdRType t) { CagdPtStruct *PtTmp, *PtLast, *Pt; if (APX_EQ(t, CrvTMin) || APX_EQ(t, CrvTMax)) return; Pt = CagdPtNew(); Pt -> Pt[0] = t; if (GlblPtList) { for (PtTmp = GlblPtList, PtLast = NULL; PtTmp != NULL; PtLast = PtTmp, PtTmp = PtTmp -> Pnext) { if (ZERO_APX_EQ(PtTmp -> Pt[0], t)) { IritFree(Pt); return; } if (PtTmp -> Pt[0] > t) break; } if (PtTmp) { /* Insert the new point in the middle of the list. */ Pt -> Pnext = PtTmp; if (PtLast) PtLast -> Pnext = Pt; else GlblPtList = Pt; } else { /* Insert the new point as the last point in the list. */ PtLast -> Pnext = Pt; } } else GlblPtList = Pt; }