/****************************************************************************** * Offset.c - computes offset approximation to curves and surfaces. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, March. 93. * ******************************************************************************/ #include "symb_loc.h" #include "extra_fn.h" #define MAX_OFFSET_IMPROVE_ITERS 20 #define NORMAL_PERTURB 1e-4 #define OFFSET_TRIM_EPS 1.25 #define OFFSET_INFLECTION_IRIT_EPS 1e-5 #define MAX_OFFSET_MATCH_ITERS 5 static CagdCrvStruct *SymbCrvCrvConvAux(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2, CagdRType OffsetDist, CagdRType Tolerance, CagdVType T1, CagdVType T2); /***************************************************************************** * DESCRIPTION: M * Given a curve and an offset amount OffsetDist, returns an approximation to M * the offset curve by offseting the control polygon in the normal direction. M * * * PARAMETERS: M * Crv: To approximate its offset curve with distance OffsetDist. M * OffsetDist: Amount of offset. Negative denotes other offset direction. M * BezInterp: If TRUE, control points are interpolated when the curve is M * reduced to a Bezier form. Otherwise, control points are M * translated OffsetDist amount only, under estimating the M * Offset. M * * * RETURN VALUE: M * CagdCrvStruct *: An approximation to the offset curve. M * * * SEE ALSO: M * SymbCrvSubdivOffset, SymbSrfOffset, SymbSrfSubdivOffset, M * SymbCrvAdapOffset, SymbCrvAdapOffsetTrim, SymbCrvLeastSquarOffset, M * SymbCrvMatchingOffset M * * * KEYWORDS: M * SymbCrvOffset, offset M *****************************************************************************/ CagdCrvStruct *SymbCrvOffset(CagdCrvStruct *Crv, CagdRType OffsetDist, CagdBType BezInterp) { CagdBType IsBezierCurve = FALSE, IsRational = CAGD_IS_RATIONAL_CRV(Crv); int i, j, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType), Order = Crv -> Order, Length = Crv -> Length; CagdBType HasNewKV = FALSE; CagdVecStruct *N; CagdCrvStruct *OffsetCrv = CagdCrvCopy(Crv); CagdRType *Nodes, *NodePtr, **Points = OffsetCrv -> Points, *KV = NULL; switch (Crv -> GType) { case CAGD_CBEZIER_TYPE: HasNewKV = TRUE; KV = BspKnotUniformOpen(Length, Order, NULL); IsBezierCurve = TRUE; break; case CAGD_CBSPLINE_TYPE: KV = Crv -> KnotVector; IsBezierCurve = Crv -> Length == Crv -> Order; break; case CAGD_CPOWER_TYPE: SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT); return NULL; default: SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV); return NULL; } Nodes = BspKnotNodes(KV, Length + Order, Order); NodePtr = Nodes; /* Interpolate the compute control points instead of under estimating */ /* This offset. */ if (BezInterp && IsBezierCurve) { CagdCrvStruct *TCrv; if (IsRational) { TCrv = CagdCoerceCrvTo(OffsetCrv, CAGD_MAKE_PT_TYPE(FALSE, CAGD_NUM_OF_PT_COORD(OffsetCrv -> PType))); CagdCrvFree(OffsetCrv); OffsetCrv = TCrv; Points = OffsetCrv -> Points; } for (j = 0; j < Length; j++, NodePtr++) { CagdRType *R = CagdCrvEval(Crv, *NodePtr); if ((N = CagdCrvNormalXY(Crv, *NodePtr, TRUE)) == NULL && (N = CagdCrvNormalXY(Crv, NodePtr == Nodes ? *NodePtr + NORMAL_PERTURB : *NodePtr - NORMAL_PERTURB, TRUE)) == NULL) { SYMB_FATAL_ERROR(SYMB_ERR_CANNOT_COMP_NORMAL); return NULL; } for (i = 1; i <= MaxCoord; i++) Points[i][j] = R[i] / (IsRational ? R[0] : 1.0) + N -> Vec[i - 1] * OffsetDist; } TCrv = CagdCrvCopy(OffsetCrv); for (i = 1; i <= MaxCoord; i++) BzrCrvInterp(TCrv -> Points[i], OffsetCrv -> Points[i], Length); CagdCrvFree(OffsetCrv); OffsetCrv = TCrv; } else { for (j = 0; j < Length; j++, NodePtr++) { if ((N = CagdCrvNormalXY(Crv, *NodePtr, TRUE)) == NULL && (N = CagdCrvNormalXY(Crv, NodePtr == Nodes ? *NodePtr + NORMAL_PERTURB : *NodePtr - NORMAL_PERTURB, TRUE)) == NULL) { SYMB_FATAL_ERROR(SYMB_ERR_CANNOT_COMP_NORMAL); return NULL; } for (i = 1; i <= MaxCoord; i++) Points[i][j] += N -> Vec[i - 1] * OffsetDist * (IsRational ? Points[W][j] : 1.0); } } if (HasNewKV) IritFree((VoidPtr) KV); IritFree((VoidPtr) Nodes); return OffsetCrv; } /***************************************************************************** * DESCRIPTION: M * Given a curve and an offset amount OffsetDist, returns an approximation to M * the offset curve by offseting the control polygon in the normal direction. M * If resulting offset is not satisfying the required tolerance the curve M * is subdivided and the algorithm recurses on both parts. M * * * PARAMETERS: M * Crv: To approximate its offset curve with distance OffsetDist. M * OffsetDist: Amount of offset. Negative denotes other offset direction. M * Tolerance: Accuracy control. M * BezInterp: If TRUE, control points are interpolated when the curve is M * reduced to a Bezier form. Otherwise, control points are M * translated OffsetDist amount only, under estimating the M * Offset. M * * * RETURN VALUE: M * CagdCrvStruct *: An approximation to the offset curve, to within M * Tolerance. M * * * SEE ALSO: M * SymbCrvOffset, SymbSrfOffset, SymbSrfSubdivOffset, SymbCrvAdapOffset, M * SymbCrvAdapOffsetTrim, SymbCrvLeastSquarOffset, SymbCrvMatchingOffset M * * * KEYWORDS: M * SymbCrvSubdivOffset, offset M *****************************************************************************/ CagdCrvStruct *SymbCrvSubdivOffset(CagdCrvStruct *Crv, CagdRType OffsetDist, CagdRType Tolerance, CagdBType BezInterp) { CagdCrvStruct *OffCrv = SymbCrvOffset(Crv, OffsetDist, BezInterp), *Dist = SymbCrvSub(Crv, OffCrv), *DistSqr = SymbCrvDotProd(Dist, Dist); CagdRType *R, MinVal, MaxVal, TMin, TMax; CagdCrvFree(Dist); R = SymbExtremumCntPtVals(DistSqr -> Points, DistSqr -> Length, TRUE); MinVal = R[1] < 0.0 ? 0.0 : sqrt(R[1]); R = SymbExtremumCntPtVals(DistSqr -> Points, DistSqr -> Length, FALSE); MaxVal = R[1] < 0.0 ? 0.0 : sqrt(R[1]); CagdCrvFree(DistSqr); CagdCrvDomain(Crv, &TMin, &TMax); if ((ABS(MinVal - ABS(OffsetDist)) > Tolerance || ABS(MaxVal - ABS(OffsetDist)) > Tolerance) && TMax - TMin > NORMAL_PERTURB * 10.0) { CagdBType NewCrv; CagdCrvStruct *Crv1, *Crv2, *OffCrv1, *OffCrv2; if (CAGD_IS_BEZIER_CRV(Crv)) { /* Promote to a Bspline so we can keep track of parametrization. */ Crv = CnvrtBezier2BsplineCrv(Crv); NewCrv = TRUE; } else NewCrv = FALSE; if (TMax - TMin > NORMAL_PERTURB * 10.0) { CagdCrvFree(OffCrv); Crv1 = CagdCrvSubdivAtParam(Crv, (TMin + TMax) / 2.0); Crv2 = Crv1 -> Pnext; Crv1 -> Pnext = NULL; OffCrv1 = SymbCrvSubdivOffset(Crv1, OffsetDist, Tolerance, BezInterp); OffCrv2 = SymbCrvSubdivOffset(Crv2, OffsetDist, Tolerance, BezInterp); CagdCrvFree(Crv1); CagdCrvFree(Crv2); OffCrv = CagdMergeCrvCrv(OffCrv1, OffCrv2, TRUE); CagdCrvFree(OffCrv1); CagdCrvFree(OffCrv2); } if (NewCrv) CagdCrvFree(Crv); } return OffCrv; } /***************************************************************************** * DESCRIPTION: M * Given a surface and an offset amount OffsetDist, returns an approximation M * to the offset surface by offseting the control mesh in the normal M * direction. M * * * PARAMETERS: M * Srf: To approximate its offset surface with distance OffsetDist. M * OffsetDist: Amount of offset. Negative denotes other offset direction. M * * * RETURN VALUE: M * CagdSrfStruct *: An approximation to the offset surface. M * * * SEE ALSO: M * SymbCrvOffset, SymbCrvSubdivOffset, SymbSrfSubdivOffset, M * SymbCrvAdapOffset, SymbCrvAdapOffsetTrim, SymbCrvLeastSquarOffset, M * SymbCrvMatchingOffset M * * * KEYWORDS: M * SymbSrfOffset, offset M *****************************************************************************/ CagdSrfStruct *SymbSrfOffset(CagdSrfStruct *Srf, CagdRType OffsetDist) { CagdBType IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf); int i, Row, Col, MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType), UOrder = Srf -> UOrder, VOrder = Srf -> VOrder, ULength = Srf -> ULength, VLength = Srf -> VLength; CagdBType HasNewKV = FALSE; CagdVecStruct *N; CagdSrfStruct *OffsetSrf = CagdSrfCopy(Srf); CagdRType *UNodes, *VNodes, *UNodePtr, *VNodePtr, **Points = OffsetSrf -> Points, *UKV = NULL, *VKV = NULL; switch (Srf -> GType) { case CAGD_SBEZIER_TYPE: HasNewKV = TRUE; UKV = BspKnotUniformOpen(ULength, UOrder, NULL); VKV = BspKnotUniformOpen(VLength, VOrder, NULL); break; case CAGD_SBSPLINE_TYPE: UKV = Srf -> UKnotVector; VKV = Srf -> VKnotVector; break; case CAGD_SPOWER_TYPE: SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT); return NULL; default: SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV); return NULL; } UNodes = BspKnotNodes(UKV, ULength + UOrder, UOrder); VNodes = BspKnotNodes(VKV, VLength + VOrder, VOrder); if (IsNotRational) for (Row = 0, VNodePtr = VNodes; Row < VLength; Row++, VNodePtr++) for (Col = 0, UNodePtr = UNodes; Col < ULength; Col++, UNodePtr++) { N = CagdSrfNormal(Srf, *UNodePtr, *VNodePtr, TRUE); for (i = 1; i <= MaxCoord; i++) Points[i][CAGD_MESH_UV(OffsetSrf, Col, Row)] += N -> Vec[i - 1] * OffsetDist; } else for (Row = 0, VNodePtr = VNodes; Row < VLength; Row++, VNodePtr++) for (Col = 0, UNodePtr = UNodes; Col < ULength; Col++, UNodePtr++) { N = CagdSrfNormal(Srf, *UNodePtr, *VNodePtr, TRUE); for (i = 1; i <= MaxCoord; i++) Points[i][CAGD_MESH_UV(OffsetSrf, Col, Row)] += N -> Vec[i - 1] * OffsetDist * Points[W][CAGD_MESH_UV(OffsetSrf, Col, Row)]; } if (HasNewKV) { IritFree((VoidPtr) UKV); IritFree((VoidPtr) VKV); } IritFree((VoidPtr) UNodes); IritFree((VoidPtr) VNodes); return OffsetSrf; } /***************************************************************************** * DESCRIPTION: M * Given a surface and an offset amount OffsetDist, returns an approximation M * to the offset surface by offseting the control mesh in the normal M * direction. M * If resulting offset is not satisfying the required tolerance the surface M * is subdivided and the algorithm recurses on both parts. M * * * PARAMETERS: M * Srf: To approximate its offset surface with distance OffsetDist. M * OffsetDist: Amount of offset. Negative denotes other offset direction. M * Tolerance: Accuracy control. M * * * RETURN VALUE: M * CagdSrfStruct *: An approximation to the offset surface, to within M * Tolerance. M * * * SEE ALSO: M * SymbCrvOffset, SymbCrvSubdivOffset, SymbSrfOffset, SymbCrvAdapOffset, M * SymbCrvAdapOffsetTrim, SymbCrvLeastSquarOffset, SymbCrvMatchingOffset M * * * KEYWORDS: M * SymbSrfSubdivOffset, offset M *****************************************************************************/ CagdSrfStruct *SymbSrfSubdivOffset(CagdSrfStruct *Srf, CagdRType OffsetDist, CagdRType Tolerance) { CagdSrfStruct *OffSrf = SymbSrfOffset(Srf, OffsetDist), *Dist = SymbSrfSub(Srf, OffSrf), *DistSqr = SymbSrfDotProd(Dist, Dist); CagdRType *R, MinVal, MaxVal; CagdSrfDirType Dir; CagdSrfFree(Dist); R = SymbExtremumCntPtVals(DistSqr -> Points, DistSqr -> ULength * DistSqr -> VLength, TRUE); MinVal = R[1] < 0.0 ? 0.0 : sqrt(R[1]); R = SymbExtremumCntPtVals(DistSqr -> Points, DistSqr -> ULength * DistSqr -> VLength, FALSE); MaxVal = R[1] < 0.0 ? 0.0 : sqrt(R[1]); CagdSrfFree(DistSqr); if (ABS(MinVal - ABS(OffsetDist)) > Tolerance || ABS(MaxVal - ABS(OffsetDist)) > Tolerance) { CagdSrfStruct *Srf1, *Srf2, *OffSrf1, *OffSrf2; CagdRType UMin, UMax, VMin, VMax, t; CagdBType NewSrf; /* Make sure it is a Bspline surface. */ if (CAGD_IS_BEZIER_SRF(Srf)) { /* Promote to a Bspline so we can keep track of parametrization. */ Srf = CnvrtBezier2BsplineSrf(Srf); NewSrf = TRUE; } else NewSrf = FALSE; CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax); if (MAX(UMax - UMin, VMax - VMin) > NORMAL_PERTURB * 10.0) { CagdSrfFree(OffSrf); if (UMax - UMin > VMax - VMin) { Dir = CAGD_CONST_U_DIR; t = (UMin + UMax) / 2.0; } else { Dir = CAGD_CONST_V_DIR; t = (VMin + VMax) / 2.0; } Srf1 = CagdSrfSubdivAtParam(Srf, t, Dir); Srf2 = Srf1 -> Pnext; Srf1 -> Pnext = NULL; OffSrf1 = SymbSrfSubdivOffset(Srf1, OffsetDist, Tolerance); OffSrf2 = SymbSrfSubdivOffset(Srf2, OffsetDist, Tolerance); CagdSrfFree(Srf1); CagdSrfFree(Srf2); OffSrf = CagdMergeSrfSrf(OffSrf1, OffSrf2, Dir, TRUE, TRUE); CagdSrfFree(OffSrf1); CagdSrfFree(OffSrf2); } if (NewSrf) CagdSrfFree(Srf); } return OffSrf; } /***************************************************************************** * DESCRIPTION: M * Given a curve and an offset amount OffsetDist, returns an approximation to M * the offset curve by offseting the control polygon in the normal direction. M * This function computes an approximation to the offset using M * OffsetAprxFunc, measure the error and use it to refine and decrease the M * error adaptively. M * Bezier curves are promoted to Bsplines curves. M * See also: Gershon Elber and Elaine Cohen, "Error Bounded Variable M * Distance Offset Operator for Free Form Curves and Surfaces". International M * Journal of Computational Geometry & Applications, Vol. 1, Num. 1, March M * 1991, pp 67-78. M * * * PARAMETERS: M * OrigCrv: To approximate its offset curve with distance OffsetDist. M * OffsetDist: Amount of offset. Negative denotes other offset direction. M * OffsetError: Tolerance control. M * OffsetAprxFunc: A function that can be used to approximate an offset M * of a curve. If NULL SymbCrvOffset function is selected. M * BezInterp: If TRUE, control points are interpolated when the curve is M * reduced to a Bezier form. Otherwise, control points are M * translated OffsetDist amount only, under estimating the M * Offset. M * * * RETURN VALUE: M * CagdCrvStruct *: An approximation to the offset curve, to within M * OffsetError. M * * * SEE ALSO: M * SymbCrvOffset, SymbCrvSubdivOffset, SymbSrfOffset, SymbSrfSubdivOffset, M * SymbCrvAdapOffsetTrim, SymbCrvLeastSquarOffset, SymbCrvMatchingOffset M * * * KEYWORDS: M * SymbCrvAdapOffset, offset M *****************************************************************************/ CagdCrvStruct *SymbCrvAdapOffset(CagdCrvStruct *OrigCrv, CagdRType OffsetDist, CagdRType OffsetError, SymbOffCrvFuncType OffsetAprxFunc, CagdBType BezInterp) { int i; CagdBType IsRational = CAGD_IS_RATIONAL_CRV(OrigCrv); CagdRType Min, Max, TMin, TMax, OffsetDistSqr = SQR(OffsetDist); CagdCrvStruct *OffsetCrv, *Crv; switch (OrigCrv -> GType) { case CAGD_CBEZIER_TYPE: Crv = CnvrtBezier2BsplineCrv(OrigCrv); break; case CAGD_CBSPLINE_TYPE: Crv = CagdCrvCopy(OrigCrv); break; case CAGD_CPOWER_TYPE: default: SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV); Crv = NULL; break; } if (OffsetAprxFunc == NULL) OffsetAprxFunc = SymbCrvOffset; CagdCrvDomain(Crv, &TMin, &TMax); for (i = 0; i < MAX_OFFSET_IMPROVE_ITERS; i++) { CagdCrvStruct *DiffCrv, *DistSqrCrv; OffsetCrv = OffsetAprxFunc(Crv, OffsetDist, BezInterp); DiffCrv = SymbCrvSub(OffsetCrv, Crv); DistSqrCrv = SymbCrvDotProd(DiffCrv, DiffCrv); CagdCrvFree(DiffCrv); CagdCrvMinMax(DistSqrCrv, 1, &Min, &Max); if (OffsetDistSqr - Min < OffsetError && Max - OffsetDistSqr < OffsetError) { /* Error is within bounds - returns this offset approximation. */ CagdCrvFree(DistSqrCrv); break; } else { /* Refine in regions where the error is too large. */ int j, k, Length = DistSqrCrv -> Length, Order = DistSqrCrv -> Order, KVLen = Length + Order; CagdRType *KV = DistSqrCrv -> KnotVector, *Nodes = BspKnotNodes(KV, KVLen, Order), *RefKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * 2 * Length); for (j = k = 0; j < Length; j++) { CagdRType *Pt = CagdCrvEval(DistSqrCrv, Nodes[j]), V = OffsetDistSqr - (IsRational ? Pt[1] / Pt[0] : Pt[1]); if (ABS(V) > OffsetError) { int Index = BspKnotLastIndexLE(KV, KVLen, Nodes[j]); if (APX_EQ(KV[Index], Nodes[j])) { if (j > 0) RefKV[k++] = (Nodes[j] + Nodes[j - 1]) / 2.0; if (j < Length - 1) RefKV[k++] = (Nodes[j] + Nodes[j + 1]) / 2.0; } else RefKV[k++] = Nodes[j]; } } CagdCrvFree(DistSqrCrv); IritFree((VoidPtr) Nodes); if (k == 0) { /* No refinement was found needed - return current curve. */ IritFree((VoidPtr) RefKV); break; } else { CagdCrvStruct *CTmp = CagdCrvRefineAtParams(Crv, FALSE, RefKV, k); IritFree((VoidPtr) RefKV); CagdCrvFree(Crv); Crv = CTmp; } } } CagdCrvFree(Crv); return OffsetCrv; } /***************************************************************************** * DESCRIPTION: M * Same function as CagdCrvAdapOffset, but trims the self intersection loops. M * See also: Gershon Elber and Elaine Cohen, "Error Bounded Variable M * Distance Offset Operator for Free Form Curves and Surfaces". International M * Journal of Computational Geometry & Applications, Vol. 1, Num. 1, March M * 1991, pp 67-78. M * * * PARAMETERS: M * OrigCrv: To approximate its offset curve with distance OffsetDist. M * OffsetDist: Amount of offset. Negative denotes other offset direction. M * OffsetError: Tolerance control. M * OffsetAprxFunc: A function that can be used to approximate an offset M * of a curve. If NULL SymbCrvOffset function is selected. M * Third parameter of SymbOffCrvFuncType is optional. M * BezInterp: If TRUE, control points are interpolated when the curve is M * reduced to a Bezier form. Otherwise, control points are M * translated OffsetDist amount only, under estimating the M * Offset. M * * * RETURN VALUE: M * CagdCrvStruct *: An approximation to the offset curve, to within M * OffsetError. M * * * SEE ALSO: M * SymbCrvOffset, SymbCrvSubdivOffset, SymbSrfOffset, SymbSrfSubdivOffset, M * SymbCrvAdapOffset, SymbCrvLeastSquarOffset, SymbCrvMatchingOffset M * * * KEYWORDS: M * SymbCrvAdapOffsetTrim, offset M *****************************************************************************/ CagdCrvStruct *SymbCrvAdapOffsetTrim(CagdCrvStruct *OrigCrv, CagdRType OffsetDist, CagdRType OffsetError, SymbOffCrvFuncType OffsetAprxFunc, CagdBType BezInterp) { CagdBType IsRational = CAGD_IS_RATIONAL_CRV(OrigCrv); CagdCrvStruct *CrvtrSqr, *CrvtrSign, *Crv, *Crvs, *LastCrv, *CrvOffsetList = NULL; CagdPtStruct *Pts, *PtsHead, *LastPts; CagdRType OffsetDistSqr1 = 1.0 / SQR(OffsetDist); switch (OrigCrv -> GType) { case CAGD_CBEZIER_TYPE: Crv = CnvrtBezier2BsplineCrv(OrigCrv); break; case CAGD_CBSPLINE_TYPE: Crv = CagdCrvCopy(OrigCrv); break; default: SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV); Crv = NULL; break; } if (OffsetDist == 0.0) return Crv; CrvtrSqr = SymbCrv2DCurvatureSqr(Crv); CrvtrSign = SymbCrv2DCurvatureSign(Crv); PtsHead = SymbCrvConstSet(CrvtrSqr, 1, SQR(OffsetError), OffsetDistSqr1 / SQR(OFFSET_TRIM_EPS)); for (Pts = PtsHead, LastPts = NULL; Pts != NULL;) { CagdRType *CrvtrSignPt, CrvtrSignVal; CrvtrSignPt = CagdCrvEval(CrvtrSign, Pts -> Pt[0]); CrvtrSignVal = IsRational ? CrvtrSignPt[1] / CrvtrSignPt[0] : CrvtrSignPt[1]; if (CrvtrSignVal * OffsetDist > 0.0) { /* Remove point from list. */ if (LastPts != NULL) { LastPts -> Pnext = Pts -> Pnext; CagdPtFree(Pts); Pts = LastPts -> Pnext; } else { PtsHead = PtsHead -> Pnext; CagdPtFree(Pts); Pts = PtsHead; } } else { LastPts = Pts; Pts = Pts -> Pnext; } } for (Pts = PtsHead; Pts != NULL || Crv != NULL; Pts = (Pts != NULL ? Pts -> Pnext : NULL)) { CagdCrvStruct *Crv1 = Pts ? CagdCrvSubdivAtParam(Crv, Pts -> Pt[0]) : CagdCrvCopy(Crv), *Crv2 = Pts ? Crv1 -> Pnext : NULL; CagdRType Min, Max, *CrvtrSqrPt, CrvtrSqrVal, *CrvtrSignPt, CrvtrSignVal; CagdCrvDomain(Crv1, &Min, &Max); CrvtrSqrPt = CagdCrvEval(CrvtrSqr, (Min + Max) / 2); CrvtrSqrVal = CrvtrSqrPt[1] / CrvtrSqrPt[0]; CrvtrSignPt = CagdCrvEval(CrvtrSign, (Min + Max) / 2); CrvtrSignVal = IsRational ? CrvtrSignPt[1] / CrvtrSignPt[0] : CrvtrSignPt[1]; if (CrvtrSqrVal < OffsetDistSqr1 / SQR(OFFSET_TRIM_EPS) || CrvtrSignVal * OffsetDist > 0.0) { CagdCrvStruct *Crv1Off = SymbCrvAdapOffset(Crv1, OffsetDist, OffsetError, OffsetAprxFunc, BezInterp); LIST_PUSH(Crv1Off, CrvOffsetList); } CagdCrvFree(Crv1); CagdCrvFree(Crv); Crv = Crv2; } Crvs = CagdListReverse(CrvOffsetList); CrvOffsetList = NULL; LastCrv = Crvs; Crvs = Crvs -> Pnext; for (Crv = Crvs; Crv != NULL; Crv = Crv -> Pnext) { CagdCrvStruct *CTmp, *CTmp2; CagdSrfStruct *DistSrf = SymbSrfDistCrvCrv(LastCrv, Crv); CagdPtStruct *IPts = SymbSrfDistFindPoints(DistSrf, OffsetError, FALSE); CagdSrfFree(DistSrf); if (IPts != NULL) { if (IPts -> Pnext) { SYMB_FATAL_ERROR(SYMB_ERR_TOO_COMPLEX); return NULL; } else { CagdRType TMin, TMax; CagdCrvDomain(LastCrv, &TMin, &TMax); CTmp = CagdCrvRegionFromCrv(LastCrv, TMin, IPts -> Pt[0]); CagdCrvFree(LastCrv); LastCrv = CTmp; CagdCrvDomain(Crv, &TMin, &TMax); CTmp = CagdCrvRegionFromCrv(Crv, IPts -> Pt[1], TMax); CTmp2 = CagdMergeCrvCrv(LastCrv, CTmp, FALSE); CagdCrvFree(CTmp); CagdCrvFree(LastCrv); LastCrv = CTmp2; } } else { /* Simply chain the pieces together. */ CTmp = CagdMergeCrvCrv(LastCrv, Crv, FALSE); CagdCrvFree(LastCrv); LastCrv = CTmp; } } CagdPtFreeList(PtsHead); CagdCrvFree(CrvtrSqr); CagdCrvFree(CrvtrSign); return LastCrv; } /***************************************************************************** * DESCRIPTION: M * Given a curve and an offset amount OffsetDist, returns an approximation to M * the offset curve by least square fitting a curve to samples taken on the M * offset curve. M * Resulting curve of order Order (degree of Crv if Order == 0) will have M * NumOfDOF control points that least sqaure fit NumOfSamples samples on the M * offset curve. M * Tolerance will be updated to hold an error distance measure. M * * * PARAMETERS: M * Crv: To approximate its offset curve with distance OffsetDist. M * OffsetDist: Amount of offset. Negative denotes other offset direction. M * NumOfSamples: Number of samples to sample the offset curve at. M * NumOfDOF: Number of degrees of freedom on the newly computed offset M * approximation. This is thesame as the number of control M * points the new curve will have. M * Order: Of the newly constructed offset approximation. If equal to M * zero, the order of Crv will be used. M * Tolerance: To return an error estimate in the L-infinity norm. M * * * RETURN VALUE: M * CagdCrvStruct *: An approximation to the offset curve. M * * * SEE ALSO: M * SymbCrvOffset, SymbCrvSubdivOffset, SymbSrfOffset, SymbSrfSubdivOffset, M * SymbCrvAdapOffset, SymbCrvAdapOffsetTrim, SymbCrvMatchingOffset M * * * KEYWORDS: M * SymbCrvLeastSquarOffset, offset M *****************************************************************************/ CagdCrvStruct *SymbCrvLeastSquarOffset(CagdCrvStruct *Crv, CagdRType OffsetDist, int NumOfSamples, int NumOfDOF, int Order, CagdRType *Tolerance) { int i; CagdRType TMin, TMax, t, dt; CagdVType Tangent; CagdCrvStruct *OffApprox, *TangentCrv = CagdCrvDerive(Crv); CagdPtStruct *Pt = NULL, *PtList = NULL; /* Create NumOfSamples samples on the offset curve. */ CagdCrvDomain(Crv, &TMin, &TMax); dt = (TMax - TMin) / (NumOfSamples - 1); for (i = 0, t = TMin; i < NumOfSamples; i++, t += dt) { CagdRType *R; if (t > TMax) /* Take care of round off errors. */ t = TMax; if (PtList == NULL) PtList = Pt = CagdPtNew(); else { Pt -> Pnext = CagdPtNew(); Pt = Pt -> Pnext; } R = CagdCrvEval(Crv, t); CagdCoerceToE3(Pt -> Pt, &R, -1, Crv -> PType); R = CagdCrvEval(TangentCrv, t); CagdCoerceToE2(Tangent, &R, -1, TangentCrv -> PType); Tangent[2] = 0.0; PT_NORMALIZE(Tangent); Pt -> Pt[0] += Tangent[1] * OffsetDist; Pt -> Pt[1] -= Tangent[0] * OffsetDist; } OffApprox = BspCrvInterpPts(PtList, Order == 0 ? Crv -> Order : Order, MIN(NumOfDOF, NumOfSamples), CAGD_UNIFORM_PARAM, Crv -> Periodic); *Tolerance = BspCrvInterpPtsError(OffApprox, PtList, CAGD_UNIFORM_PARAM, Crv -> Periodic); CagdPtFreeList(PtList); CagdCrvFree(TangentCrv); return OffApprox; } /***************************************************************************** * DESCRIPTION: M * Computes an offset to a freeform curve using matching of tangent fields. M * The given curve is split at all its inflection points, made sure it spans M * less than 90 degrees, and then is matched against an arc of the proper M * angular span of tangents. M * Unlike other offset methods, this method allways preserves the distance M * between the original curve ans its offset. The error in this methods can M * surface only in the non orthogonality of the offset direction. M * * * PARAMETERS: M * Crv: To approximate its offset curve with distance OffsetDist. M * OffsetDist: Amount of offset. Negative denotes other offset direction. M * Tolerance: Of angular discrepancy that is allowed. M * * * RETURN VALUE: M * CagdCrvStruct *: The offset curve approximation. M * * * SEE ALSO: M * SymbCrvOffset, SymbCrvSubdivOffset, SymbSrfOffset, SymbSrfSubdivOffset, M * SymbCrvAdapOffset, SymbCrvAdapOffsetTrim, SymbCrvLeastSquarOffset, M * SymbCrvCrvConvolution M * * * KEYWORDS: M * SymbCrvMatchingOffset M *****************************************************************************/ CagdCrvStruct *SymbCrvMatchingOffset(CagdCrvStruct *Crv, CagdRType OffsetDist, CagdRType Tolerance) { CagdCrvStruct *OffCrv, *OffCrvList = NULL, *OrigCrv = Crv; CagdPtStruct *Pt, *Pts = SymbCrv2DInflectionPts(Crv, OFFSET_INFLECTION_IRIT_EPS); for (Pt = Pts; Pt != NULL; Pt = Pt -> Pnext) { CagdCrvStruct *Crv1 = CagdCrvSubdivAtParam(Crv, Pt -> Pt[0]), *Crv2 = Crv1 -> Pnext, *Off1Crv = SymbCrvCrvConvolution(Crv1, NULL, OffsetDist, Tolerance); LIST_PUSH(Off1Crv, OffCrvList); Crv = Crv2; CagdCrvFree(Crv1); } CagdPtFreeList(Pts); OffCrv = SymbCrvCrvConvolution(Crv, NULL, OffsetDist, Tolerance); LIST_PUSH(OffCrv, OffCrvList); if (Crv != OrigCrv) CagdCrvFree(Crv); OffCrvList = CagdListReverse(OffCrvList); OffCrv = CagdMergeCrvList(OffCrvList, TRUE); CagdCrvFreeList(OffCrvList); return OffCrv; } /***************************************************************************** * DESCRIPTION: M * Computes the convolution of the given two curves by matching their M * tangents and reparametrizing Crv2. M * If Crv2 is NULL, an Arc of radius OffsetDist is used, resulting in an M * offset operation of Crv1. M * Both Crv1 and Crv2 are assumed to have no inflection points and to M * span the same angular domain. That is Crv1'(0) || Crv2'(0) and similarly M * Crv1'(1) || Crv2'(1), where || denotes parallel. M * * * PARAMETERS: M * Crv1, Crv2: The two curves to convolve. M * OffsetDist: Amount of offset, if Crv2 == NULL. Negative value denotes M * other offset/convolution direction. M * Tolerance: Of angular discrepancy that is allowed. M * * * RETURN VALUE: M * CagdCrvStruct *: The offset curve approximation. M * * * SEE ALSO: M * SymbCrvOffset, SymbCrvSubdivOffset, SymbSrfOffset, SymbSrfSubdivOffset, M * SymbCrvAdapOffset, SymbCrvAdapOffsetTrim, SymbCrvLeastSquarOffset, M * SymbCrvMatchingOffset M * * * KEYWORDS: M * SymbCrvCrvConvolution M *****************************************************************************/ CagdCrvStruct *SymbCrvCrvConvolution(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2, CagdRType OffsetDist, CagdRType Tolerance) { int i; CagdCrvStruct *Crv1E2 = CagdCoerceCrvTo(Crv1, CAGD_PT_E2_TYPE), *DCrv1E2 = CagdCrvDerive(Crv1E2); CagdRType **Points = DCrv1E2 -> Points, *XPts = Points[1], *YPts = Points[2]; CagdCrvFree(Crv1E2); /* Scans the curve's Hodograph for a vector that is more than 90 degrees */ /* away from the initial vector, and if found one, must subdivide Crv. */ for (i = 1; i < DCrv1E2 -> Length; i++) if (XPts[0] * XPts[i] + YPts[0] * YPts[i] < 0) break; if (i < DCrv1E2 -> Length) { CagdRType TMin, TMax; CagdCrvStruct *ConvCrv, *Crv1Subdiv, *ConvCrv1Sub1, *Crv2Subdiv, *ConvCrv1Sub2; CagdCrvDomain(Crv1, &TMin, &TMax); Crv1Subdiv = CagdCrvSubdivAtParam(Crv1, (TMin + TMax) / 2.0); if (Crv2 != NULL) Crv2Subdiv = CagdCrvSubdivAtParam(Crv2, (TMin + TMax) / 2.0); ConvCrv1Sub1 = SymbCrvCrvConvolution(Crv1Subdiv, Crv2 != NULL ? Crv2Subdiv : NULL, OffsetDist, Tolerance); ConvCrv1Sub2 = SymbCrvCrvConvolution(Crv1Subdiv -> Pnext, Crv2 != NULL ? Crv2Subdiv -> Pnext : NULL, OffsetDist, Tolerance); CagdCrvFreeList(Crv1Subdiv); if (Crv2 != NULL) CagdCrvFreeList(Crv2Subdiv); ConvCrv = CagdMergeCrvCrv(ConvCrv1Sub1, ConvCrv1Sub2, TRUE); CagdCrvFree(ConvCrv1Sub1); CagdCrvFree(ConvCrv1Sub2); CagdCrvFree(DCrv1E2); return ConvCrv; } else { /* Curve does not span more than 90 degs - approximate by matching. */ CagdVType T1, T2; T1[0] = -YPts[0]; T1[1] = XPts[0]; T1[2] = 0.0; T2[0] = -YPts[DCrv1E2 -> Length - 1]; T2[1] = XPts[DCrv1E2 -> Length - 1]; T2[2] = 0.0; CagdCrvFree(DCrv1E2); return SymbCrvCrvConvAux(Crv1, Crv2, OffsetDist, Tolerance, T1, T2); } } /***************************************************************************** * DESCRIPTION: * * An Auxiliary function of SymbCrvCrvConvolution. Matches a given curve(s) * * with no inflection points and angular tangent span between T1 and T2 of * * less than 90 degrees to an arc of same angular span and return the * * approximated convolved curve. * *****************************************************************************/ static CagdCrvStruct *SymbCrvCrvConvAux(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2, CagdRType OffsetDist, CagdRType Tolerance, CagdVType T1, CagdVType T2) { CagdBType NewCrv2; int NumIters = 0, Reduce = 3, SampleSet = 15; CagdRType MaxError = 0.0; CagdPtStruct Start, Center, End; CagdCrvStruct *CTmp, *Crv1Match, *ConvCrv; if (Crv2 == NULL) { /* Constructs the proper arc with the right orientation as detected */ /* by the cross product of the two angular tangents T1, and T2. */ NewCrv2 = TRUE; PT_NORMALIZE(T1); PT_SCALE(T1, fabs(OffsetDist)); PT_NORMALIZE(T2); PT_SCALE(T2, fabs(OffsetDist)); PT_COPY(Start.Pt, T1); PT_COPY(End.Pt, T2); PT_RESET(Center.Pt); CTmp = BzrCrvCreateArc(&Start, &Center, &End); Crv2 = CagdCoerceCrvTo(CTmp, CAGD_PT_P2_TYPE); CagdCrvFree(CTmp); if (T1[0] * T2[1] - T1[1] * T2[0] > 0) { /* Reverse the generated arc by 180 rotation along Z. */ CagdVType Translate; PT_RESET(Translate); CagdCrvTransform(Crv2, Translate, -1); } else OffsetDist *= -1.0; } else NewCrv2 = FALSE; do { int i; CagdCrvStruct *DCrv1Match, *ErrorCrv, *ErrorCrvSqr, *ErrorCrvE1, *DCrv1MatchLenSqr, *RecipDCrv1MatchLenSqr; CagdRType *Pts; if ((Crv1Match = CagdMatchingTwoCurves(Crv2, Crv1, Reduce, SampleSet, 2, FALSE, FALSE, NULL)) != NULL) { DCrv1Match = CagdCrvDerive(Crv1Match); /* Compute the error functional as */ /* ^2 */ /* ErrorCrv(t) = ---------------------------- */ /* */ DCrv1MatchLenSqr = SymbCrvDotProd(DCrv1Match, DCrv1Match); RecipDCrv1MatchLenSqr = SymbCrvInvert(DCrv1MatchLenSqr); CagdCrvFree(DCrv1MatchLenSqr); ErrorCrv = SymbCrvDotProd(Crv2, DCrv1Match); CagdCrvFree(DCrv1Match); ErrorCrvSqr = SymbCrvMult(ErrorCrv, ErrorCrv); CagdCrvFree(ErrorCrv); ErrorCrv = SymbCrvMult(ErrorCrvSqr, RecipDCrv1MatchLenSqr); CagdCrvFree(ErrorCrvSqr); CagdCrvFree(RecipDCrv1MatchLenSqr); ErrorCrvE1 = CagdCoerceCrvTo(ErrorCrv, CAGD_PT_E1_TYPE); CagdCrvFree(ErrorCrv); MaxError = 0.0; /* Get a bound on the maximal angular error. */ for (i = 0, Pts = ErrorCrvE1 -> Points[1]; i < ErrorCrvE1 -> Length; i++, Pts++) { if (MaxError < *Pts) MaxError = *Pts; } CagdCrvFree(ErrorCrvE1); MaxError = sqrt(MaxError); MaxError /= fabs(OffsetDist); /* Convert to degrees relative to the 90 orthogonal. */ MaxError = fabs((M_PI / 2 - acos(MaxError)) * 180 / M_PI); } else MaxError = IRIT_INFNTY; if (MaxError > Tolerance && ++NumIters < MAX_OFFSET_MATCH_ITERS) { /* Going for another round - free current match/set new params. */ if (Crv1Match != NULL) CagdCrvFree(Crv1Match); Reduce += Reduce; SampleSet += SampleSet; } } while (MaxError > Tolerance && NumIters < MAX_OFFSET_MATCH_ITERS); if (OffsetDist > 0) ConvCrv = SymbCrvAdd(Crv1Match, Crv2); else ConvCrv = SymbCrvSub(Crv1Match, Crv2); if (NewCrv2) CagdCrvFree(Crv2); CagdCrvFree(Crv1Match); return ConvCrv; } /***************************************************************************** * DESCRIPTION: M * Computes an elevated surface emenating from the given C^1 continuous M * curve in all directions like a fire front. The surface gets away from Crv M * in a slope of 45 degrees. This elevated surface is an approximation of M * the real envelope only, as prescribed by Tolerance. M * If the given curve is closed, it is assume to be C^1 at the end point as M * well. For a close curve, two surfaces are actually returned - one for the M * inside and one for the outside firefront. M * This function employs SymbCrvSubdivOffset for the offset computations. M * * * PARAMETERS: M * Crv: The curve to process. M * Height: The height of the elevated surface (also the width of the M * offset operation. M * Tolerance: Accuracy of the elevated surface approximation. M * * * RETURN VALUE: M * CagdSrfStruct *: A freeform surface approximating the elevated surface M * for open curve Crv, or two surfaces for the case of closed M * curve Crv. M * * * KEYWORDS: M * SymbEnvOffsetFromCrv M *****************************************************************************/ CagdSrfStruct *SymbEnvOffsetFromCrv(CagdCrvStruct *Crv, CagdRType Height, CagdRType Tolerance) { CagdBType ClosedCurve, IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); int MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdPType Trans, PtStart, PtEnd; CagdRType TMin, TMax, *R; CagdCrvStruct *TCrv, *Crv1, *Crv2, *CrvOff1, *CrvOff2; CagdSrfStruct *Srf1, *Srf2, *TSrf, *TSrf2; /* Check if the curve is closed or not. */ CagdCrvDomain(Crv, &TMin, &TMax); R = CagdCrvEval(Crv, TMin); CagdCoerceToE3(PtStart, &R, -1, Crv -> PType); R = CagdCrvEval(Crv, TMax); CagdCoerceToE3(PtEnd, &R, -1, Crv -> PType); ClosedCurve = PT_APX_EQ(PtStart, PtEnd); /* Make sure we have a three dimensional curve. */ if (MaxCoord < 3) { Crv1 = CagdCoerceCrvTo(Crv, IsNotRational ? CAGD_PT_E3_TYPE : CAGD_PT_P3_TYPE); } else Crv1 = CagdCrvCopy(Crv); Crv2 = CagdCrvReverse(Crv); /* Compute offsets in the plane and translate to proper heights. */ TCrv = SymbCrvSubdivOffset(Crv1, Height, Tolerance, FALSE), CrvOff1 = CagdCoerceCrvTo(TCrv, IsNotRational ? CAGD_PT_E3_TYPE : CAGD_PT_P3_TYPE); CagdCrvFree(TCrv); TCrv = SymbCrvSubdivOffset(Crv2, Height, Tolerance, FALSE); CrvOff2 = CagdCoerceCrvTo(TCrv, IsNotRational ? CAGD_PT_E3_TYPE : CAGD_PT_P3_TYPE); CagdCrvFree(TCrv); Trans[0] = Trans[1] = 0.0; Trans[2] = Height; CagdCrvTransform(CrvOff1, Trans, 1.0); CagdCrvTransform(CrvOff2, Trans, 1.0); /* Create the end closing half cones, for open curve segments. */ if (ClosedCurve) { int Len1, i; CagdRType **Points; /* Make sure end points of CrvOff1, CrvOff2 are identical as well. */ Points = CrvOff1 -> Points; Len1 = CrvOff1 -> Length - 1; for (i = IsNotRational; i <= MaxCoord; i++) Points[i][0] = Points[i][Len1] = (Points[i][0] + Points[i][Len1]) / 2.0; Points = CrvOff2 -> Points; Len1 = CrvOff2 -> Length - 1; for (i = IsNotRational; i <= MaxCoord; i++) Points[i][0] = Points[i][Len1] = (Points[i][0] + Points[i][Len1]) / 2.0; Srf1 = CagdRuledSrf(Crv1, CrvOff1, 2, 2); Srf2 = CagdRuledSrf(Crv2, CrvOff2, 2, 2); Srf1 -> Pnext = Srf2; } else { CagdPtStruct Pt1, Pt2; CagdCrvStruct *DiagLine; MatrixType Mat; CagdVecStruct *Tan; CagdSrfStruct *HalfCone; /* Create half a cone with the proper dimensions. */ PT_RESET(Pt1.Pt); Pt2.Pt[0] = Height; Pt2.Pt[1] = 0.0; Pt2.Pt[2] = Height; DiagLine = CagdMergePtPt(&Pt1, &Pt2); if (IsNotRational) { TSrf = CagdSurfaceRevPolynomialApprox(DiagLine); HalfCone = CagdSrfRegionFromSrf(TSrf, 0.0, 2.0, CAGD_CONST_U_DIR); } else { TSrf = CagdSurfaceRev(DiagLine); HalfCone = CagdSrfRegionFromSrf(TSrf, 0.0, 2.0, CAGD_CONST_U_DIR); } /* Compute the two side walls. */ Srf1 = CagdRuledSrf(Crv1, CrvOff1, 2, 2); Srf2 = CagdRuledSrf(Crv2, CrvOff2, 2, 2); /* Position the first half cone at the end of the first side wall */ /* and merge with the first side wall. */ TSrf = CagdSrfCopy(HalfCone); Tan = CagdCrvTangent(Crv, TMax, TRUE); MatGenMatRotZ1(atan2(Tan -> Vec[1], Tan -> Vec[0]) - M_PI / 2, Mat); CagdSrfMatTransform(TSrf, Mat); CagdSrfTransform(TSrf, PtEnd, 1.0); TSrf2 = CagdMergeSrfSrf(Srf1, TSrf, CAGD_CONST_U_DIR, TRUE, FALSE); CagdSrfFree(Srf1); CagdSrfFree(TSrf); /* Merge in the second side wall. */ Srf1 = CagdMergeSrfSrf(TSrf2, Srf2, CAGD_CONST_U_DIR, TRUE, FALSE); CagdSrfFree(TSrf2); CagdSrfFree(Srf2); /* Position the second half cone at the beginning of the first side */ /* wall and merge at the end of the second side wall. */ TSrf = HalfCone; Tan = CagdCrvTangent(Crv, TMin, TRUE); MatGenMatRotZ1(atan2(Tan -> Vec[1], Tan -> Vec[0]) + M_PI / 2, Mat); CagdSrfMatTransform(TSrf, Mat); CagdSrfTransform(TSrf, PtStart, 1.0); TSrf2 = CagdMergeSrfSrf(Srf1, TSrf, CAGD_CONST_U_DIR, TRUE, FALSE); CagdSrfFree(Srf1); CagdSrfFree(TSrf); Srf1 = TSrf2; } CagdCrvFree(Crv1); CagdCrvFree(Crv2); CagdCrvFree(CrvOff1); CagdCrvFree(CrvOff2); return Srf1; }