/****************************************************************************** * Bzr_Sym.c - Bezier symbolic computation. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Sep. 92. * ******************************************************************************/ #include "symb_loc.h" /***************************************************************************** * DESCRIPTION: M * Given two Bezier curves - multiply them coordinatewise. M * The two curves are promoted to same point type before the multiplication M * can take place. M * * * PARAMETERS: M * Crv1, Crv2: The two curves to multiply. M * * * RETURN VALUE: M * CagdCrvStruct *: The product Crv1 * Crv2 coordinatewise. M * * * KEYWORDS: M * BzrCrvMult, product M *****************************************************************************/ CagdCrvStruct *BzrCrvMult(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2) { CagdBType IsNotRational; int i, j, k, MaxCoord, Order, Order1 = Crv1 -> Order, Order2 = Crv2 -> Order, Degree1 = Order1 - 1, Degree2 = Order2 - 1; CagdCrvStruct *ProdCrv; CagdRType **Points, **Points1, **Points2; if (!CAGD_IS_BEZIER_CRV(Crv1) || !CAGD_IS_BEZIER_CRV(Crv2)) SYMB_FATAL_ERROR(SYMB_ERR_BZR_CRV_EXPECT); Crv1 = CagdCrvCopy(Crv1); Crv2 = CagdCrvCopy(Crv2); if (!CagdMakeCrvsCompatible(&Crv1, &Crv2, FALSE, FALSE)) SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT); ProdCrv = BzrCrvNew(Order = Order1 + Order2 - 1, Crv1 -> PType); IsNotRational = !CAGD_IS_RATIONAL_CRV(ProdCrv); MaxCoord = CAGD_NUM_OF_PT_COORD(ProdCrv -> PType); Points = ProdCrv -> Points; Points1 = Crv1 -> Points; Points2 = Crv2 -> Points; for (i = 0; i < Order; i++) for (k = IsNotRational; k <= MaxCoord; k++) Points[k][i] = 0.0; for (i = 0; i < Order1; i++) { for (j = 0; j < Order2; j++) { CagdRType Coef = CagdIChooseK(i, Degree1) * CagdIChooseK(j, Degree2) / CagdIChooseK(i + j, Degree1 + Degree2); for (k = IsNotRational; k <= MaxCoord; k++) Points[k][i+j] += Coef * Points1[k][i] * Points2[k][j]; } } CagdCrvFree(Crv1); CagdCrvFree(Crv2); return ProdCrv; } /***************************************************************************** * DESCRIPTION: M * Given two Bezier curve lists - multiply them one at a time. M * Return a Bezier curve lists representing their products. M * * * PARAMETERS: M * Crv1Lst: First list of Bezier curves to multiply. M * Crv2Lst: Second list of Bezier curves to multiply. M * * * RETURN VALUE: M * CagdCrvStruct *: A list of product curves M * * * KEYWORDS: M * BzrCrvMultList, product M *****************************************************************************/ CagdCrvStruct *BzrCrvMultList(CagdCrvStruct *Crv1Lst, CagdCrvStruct *Crv2Lst) { CagdCrvStruct *Crv1, *Crv2, *ProdCrvTail = NULL, *ProdCrvList = NULL; for (Crv1 = Crv1Lst, Crv2 = Crv2Lst; Crv1 != NULL && Crv2 != NULL; Crv1 = Crv1 -> Pnext, Crv2 = Crv2 -> Pnext) { CagdCrvStruct *ProdCrv = BzrCrvMult(Crv1, Crv2); if (ProdCrvList == NULL) ProdCrvList = ProdCrvTail = ProdCrv; else { ProdCrvTail -> Pnext = ProdCrv; ProdCrvTail = ProdCrv; } } return ProdCrvList; } /***************************************************************************** * DESCRIPTION: M * Given two Bezier surfaces - multiply them coordinatewise. M * The two surfaces are promoted to same point type before multiplication M * can take place. M * * * PARAMETERS: M * Srf1, Srf2: The two surfaces to multiply. M * * * RETURN VALUE: M * CagdSrfStruct *: The product Srf1 * Srf2 coordinatewise. M * * * KEYWORDS: M * BzrSrfMult, product M *****************************************************************************/ CagdSrfStruct *BzrSrfMult(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2) { CagdBType IsNotRational; int i, j, k, l, m, MaxCoord, UOrder, VOrder, UOrder1 = Srf1 -> UOrder, VOrder1 = Srf1 -> VOrder, UOrder2 = Srf2 -> UOrder, VOrder2 = Srf2 -> VOrder, UDegree1 = UOrder1 - 1, VDegree1 = VOrder1 - 1, UDegree2 = UOrder2 - 1, VDegree2 = VOrder2 - 1; CagdSrfStruct *ProdSrf; CagdRType **Points, **Points1, **Points2; if (!CAGD_IS_BEZIER_SRF(Srf1) || !CAGD_IS_BEZIER_SRF(Srf2)) SYMB_FATAL_ERROR(SYMB_ERR_BZR_SRF_EXPECT); Srf1 = CagdSrfCopy(Srf1); Srf2 = CagdSrfCopy(Srf2); if (!CagdMakeSrfsCompatible(&Srf1, &Srf2, FALSE, FALSE, FALSE, FALSE)) SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT); ProdSrf = BzrSrfNew(UOrder = UOrder1 + UOrder2 - 1, VOrder = VOrder1 + VOrder2 - 1, Srf1 -> PType); IsNotRational = !CAGD_IS_RATIONAL_SRF(ProdSrf); MaxCoord = CAGD_NUM_OF_PT_COORD(ProdSrf -> PType); Points = ProdSrf -> Points; Points1 = Srf1 -> Points; Points2 = Srf2 -> Points; for (k = IsNotRational; k <= MaxCoord; k++) { CagdRType *Pts = Points[k]; for (i = 0; i < UOrder * VOrder; i++) *Pts++ = 0.0; } for (i = 0; i < UOrder1; i++) { for (j = 0; j < VOrder1; j++) { for (l = 0; l < UOrder2; l++) { for (m = 0; m < VOrder2; m++) { int Index = CAGD_MESH_UV(ProdSrf, i + l, j + m), Index1 = CAGD_MESH_UV(Srf1, i, j), Index2 = CAGD_MESH_UV(Srf2, l, m); CagdRType Coef1 = CagdIChooseK(i, UDegree1) * CagdIChooseK(l, UDegree2) / CagdIChooseK(i + l, UDegree1 + UDegree2), Coef2 = CagdIChooseK(j, VDegree1) * CagdIChooseK(m, VDegree2) / CagdIChooseK(j + m, VDegree1 + VDegree2); for (k = IsNotRational; k <= MaxCoord; k++) Points[k][Index] += Coef1 * Coef2 * Points1[k][Index1] * Points2[k][Index2]; } } } } CagdSrfFree(Srf1); CagdSrfFree(Srf2); return ProdSrf; } /***************************************************************************** * DESCRIPTION: M * Given a rational Bezier curve - computes its derivative curve (Hodograph) M * using the quotient rule for differentiation. M * * * PARAMETERS: M * Crv: Bezier curve to differentiate. M * * * RETURN VALUE: M * CagdCrvStruct *: Differentiated rational Bezier curve. M * * * SEE ALSO: M * BzrCrvDerive, BspCrvDerive, BspCrvDeriveRational, CagdCrvDerive M * * * KEYWORDS: M * BzrCrvDeriveRational, derivatives M *****************************************************************************/ CagdCrvStruct *BzrCrvDeriveRational(CagdCrvStruct *Crv) { CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *DCrvW, *DCrvX, *DCrvY, *DCrvZ, *TCrv1, *TCrv2, *DeriveCrv; SymbCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ); if (CrvW) DCrvW = BzrCrvDerive(CrvW); else { DCrvW = NULL; SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_EXPECTED); } if (CrvX) { DCrvX = BzrCrvDerive(CrvX); TCrv1 = BzrCrvMult(DCrvX, CrvW); TCrv2 = BzrCrvMult(CrvX, DCrvW); CagdCrvFree(CrvX); CagdCrvFree(DCrvX); CrvX = SymbCrvSub(TCrv1, TCrv2); CagdCrvFree(TCrv1); CagdCrvFree(TCrv2); } if (CrvY) { DCrvY = BzrCrvDerive(CrvY); TCrv1 = BzrCrvMult(DCrvY, CrvW); TCrv2 = BzrCrvMult(CrvY, DCrvW); CagdCrvFree(CrvY); CagdCrvFree(DCrvY); CrvY = SymbCrvSub(TCrv1, TCrv2); CagdCrvFree(TCrv1); CagdCrvFree(TCrv2); } if (CrvZ) { DCrvZ = BzrCrvDerive(CrvZ); TCrv1 = BzrCrvMult(DCrvZ, CrvW); TCrv2 = BzrCrvMult(CrvZ, DCrvW); CagdCrvFree(CrvZ); CagdCrvFree(DCrvZ); CrvZ = SymbCrvSub(TCrv1, TCrv2); CagdCrvFree(TCrv1); CagdCrvFree(TCrv2); } TCrv1 = BzrCrvMult(CrvW, CrvW); CagdCrvFree(CrvW); CrvW = TCrv1; if (!CagdMakeCrvsCompatible(&CrvW, &CrvX, TRUE, TRUE) || !CagdMakeCrvsCompatible(&CrvW, &CrvY, TRUE, TRUE) || !CagdMakeCrvsCompatible(&CrvW, &CrvZ, TRUE, TRUE)) SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT); DeriveCrv = SymbCrvMergeScalar(CrvW, CrvX, CrvY, CrvZ); if (CrvX) CagdCrvFree(CrvX); if (CrvY) CagdCrvFree(CrvY); if (CrvZ) CagdCrvFree(CrvZ); if (CrvW) { CagdCrvFree(CrvW); CagdCrvFree(DCrvW); } return DeriveCrv; } /***************************************************************************** * DESCRIPTION: M * Given a rational Bezier surface - computes its derivative surface in M * direction Dir, using the quotient rule for differentiation. M * * * PARAMETERS: M * Srf: Bezier surface to differentiate. M * Dir: Direction of Differentiation. Either U or V. M * * * RETURN VALUE: M * CagdSrfStruct *: Differentiated rational Bezier surface. M * * * SEE ALSO: M * CagdSrfDerive, BzrSrfDerive, BspSrfDerive, BspSrfDeriveRational M * * * KEYWORDS: M * BzrSrfDeriveRational, derivatives M *****************************************************************************/ CagdSrfStruct *BzrSrfDeriveRational(CagdSrfStruct *Srf, CagdSrfDirType Dir) { CagdSrfStruct *SrfW, *SrfX, *SrfY, *SrfZ, *DSrfW, *DSrfX, *DSrfY, *DSrfZ, *TSrf1, *TSrf2, *DeriveSrf; SymbSrfSplitScalar(Srf, &SrfW, &SrfX, &SrfY, &SrfZ); if (SrfW) DSrfW = BzrSrfDerive(SrfW, Dir); else { DSrfW = NULL; SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_EXPECTED); } if (SrfX) { DSrfX = BzrSrfDerive(SrfX, Dir); TSrf1 = BzrSrfMult(DSrfX, SrfW); TSrf2 = BzrSrfMult(SrfX, DSrfW); CagdSrfFree(SrfX); CagdSrfFree(DSrfX); SrfX = SymbSrfSub(TSrf1, TSrf2); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); } if (SrfY) { DSrfY = BzrSrfDerive(SrfY, Dir); TSrf1 = BzrSrfMult(DSrfY, SrfW); TSrf2 = BzrSrfMult(SrfY, DSrfW); CagdSrfFree(SrfY); CagdSrfFree(DSrfY); SrfY = SymbSrfSub(TSrf1, TSrf2); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); } if (SrfZ) { DSrfZ = BzrSrfDerive(SrfZ, Dir); TSrf1 = BzrSrfMult(DSrfZ, SrfW); TSrf2 = BzrSrfMult(SrfZ, DSrfW); CagdSrfFree(SrfZ); CagdSrfFree(DSrfZ); SrfZ = SymbSrfSub(TSrf1, TSrf2); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); } TSrf1 = BzrSrfMult(SrfW, SrfW); CagdSrfFree(SrfW); SrfW = TSrf1; if (!CagdMakeSrfsCompatible(&SrfW, &SrfX, TRUE, TRUE, TRUE, TRUE) || !CagdMakeSrfsCompatible(&SrfW, &SrfY, TRUE, TRUE, TRUE, TRUE) || !CagdMakeSrfsCompatible(&SrfW, &SrfZ, TRUE, TRUE, TRUE, TRUE)) SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT); DeriveSrf = SymbSrfMergeScalar(SrfW, SrfX, SrfY, SrfZ); if (SrfX) CagdSrfFree(SrfX); if (SrfY) CagdSrfFree(SrfY); if (SrfZ) CagdSrfFree(SrfZ); if (SrfW) { CagdSrfFree(SrfW); CagdSrfFree(DSrfW); } return DeriveSrf; } /***************************************************************************** * DESCRIPTION: M * Given a Bezier curve - convert it to (possibly) piecewise cubic. M * If the curve is M * 1. A cubic - a copy if it is returned. M * 2. Lower than cubic - a degree raised (to a cubic) curve is returned. M * 3. Higher than cubic - a C^1 continuous piecewise cubic approximation M * is computed for Crv. M * M * In case 3 a list of polynomial cubic curves is returned. Tol is then used M * for the distance tolerance error measure for the approximation. M * If, however, NoRational is set, rational curves of any order will also M * be approximated using cubic polynomials. M * Furthermore if the total length of control polygon is more than MaxLen, M * the curve is subdivided until this is not the case. M * * * PARAMETERS: M * Crv: To approximate using cubic Bezier polynomials. M * Tol: Accuracy control. M * MaxLen: Maximum arc length of curve. M * NoRational: Do we want to approximate rational curves as well? M * * * RETURN VALUE: M * CagdCrvStruct *: A list of cubic Bezier polynomials approximating Crv. M * * * KEYWORDS: M * BzrApproxBzrCrvAsCubics, conversion, approximation M *****************************************************************************/ CagdCrvStruct *BzrApproxBzrCrvAsCubics(CagdCrvStruct *Crv, CagdRType Tol, CagdRType MaxLen, CagdBType NoRational) { CagdBType NewCrv = FALSE; CagdCrvStruct *TCrv, *CubicCrvs, *CubicCrvsMaxLen = NULL; if (CAGD_IS_RATIONAL_CRV(Crv) && NoRational) return BzrApproxBzrCrvAsCubicPoly(Crv, Tol * Tol); if (CAGD_IS_PERIODIC_CRV(Crv)) { NewCrv = TRUE; Crv = CnvrtPeriodic2FloatCrv(Crv); } if (CAGD_IS_BSPLINE_CRV(Crv) && !BspCrvHasOpenEC(Crv)) { TCrv = BspCrvOpenEnd(Crv); if (NewCrv) CagdCrvFree(Crv); Crv = TCrv; NewCrv = TRUE; } switch (Crv -> Order) { case 2: CubicCrvs = BzrCrvDegreeRaiseN(Crv, 4); break; case 3: CubicCrvs = BzrCrvDegreeRaise(Crv); break; case 4: CubicCrvs = CagdCrvCopy(Crv); break; default: if (Crv -> Order < 4) SYMB_FATAL_ERROR(SYMB_ERR_OUT_OF_RANGE); CubicCrvs = BzrApproxBzrCrvAsCubicPoly(Crv, Tol * Tol); } for (TCrv = CubicCrvs; TCrv != NULL; TCrv = TCrv -> Pnext) { CagdCrvStruct *TCrv2, *MaxLenCrv = SymbLimitCrvArcLen(TCrv, MaxLen); for (TCrv2 = MaxLenCrv; TCrv2 -> Pnext != NULL; TCrv2 = TCrv2 -> Pnext); TCrv2 -> Pnext = CubicCrvsMaxLen; CubicCrvsMaxLen = MaxLenCrv; } if (NewCrv) CagdCrvFree(Crv); CagdCrvFreeList(CubicCrvs); return CubicCrvsMaxLen; } /***************************************************************************** * DESCRIPTION: M * Given a Bezier curve with order larger than cubic, approximate it using M * piecewise cubic curves. M * A C^1 continuous approximation is computed so that the approximation is M * at most sqrt(Tol2) away from the given curve. M * Input curve can be rational, although output is always polynomial. M * * * PARAMETERS: M * Crv: To approximate using cubic Bezier curves. M * Tol2: Accuracy control. M * * * RETURN VALUE: M * CagdCrvStruct *: List of cubic Bezier curves approximating Crv. M * * * KEYWORDS: M * BzrApproxBzrCrvAsCubicPoly, approximation, conversion M *****************************************************************************/ CagdCrvStruct *BzrApproxBzrCrvAsCubicPoly(CagdCrvStruct *Crv, CagdRType Tol2) { CagdBType NewCrv = FALSE, ApproxOK = TRUE, IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); int i, Order = Crv -> Order, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdPointType PType = Crv -> PType, CubicPType = CAGD_MAKE_PT_TYPE(FALSE, MaxCoord); CagdCrvStruct *DistSqrCrv, *DiffCrv, *CubicBzr = BzrCrvNew(4, CubicPType); CagdRType E3Pt1[3], E3Pt2[3], E3Pt3[3], E3Pt4[3], Tan1[3], Tan2[3], **CubicPoints = CubicBzr -> Points, **Points = Crv -> Points; if (CAGD_IS_PERIODIC_CRV(Crv)) { NewCrv = TRUE; Crv = CnvrtPeriodic2FloatCrv(Crv); } if (CAGD_IS_BSPLINE_CRV(Crv) && !BspCrvHasOpenEC(Crv)) { CagdCrvStruct *TCrv = BspCrvOpenEnd(Crv); if (NewCrv) CagdCrvFree(Crv); Crv = TCrv; NewCrv = TRUE; } CagdCoerceToE3(E3Pt1, Points, 0, PType); CagdCoerceToE3(E3Pt2, Points, 1, PType); CagdCoerceToE3(E3Pt3, Points, Order - 2, PType); CagdCoerceToE3(E3Pt4, Points, Order - 1, PType); /* End of the two points must match. */ for (i = 0; i < MaxCoord; i++) { CubicPoints[i + 1][0] = E3Pt1[i]; CubicPoints[i + 1][3] = E3Pt4[i]; } /* Tangents at the end of the two points must match. */ PT_SUB(Tan1, E3Pt2, E3Pt1); PT_SUB(Tan2, E3Pt4, E3Pt3); PT_SCALE(Tan1, (Order - 1.0) / 3.0); PT_SCALE(Tan2, (Order - 1.0) / 3.0); for (i = 0; i < MaxCoord; i++) { CubicPoints[i + 1][1] = E3Pt1[i] + Tan1[i]; CubicPoints[i + 1][2] = E3Pt4[i] - Tan2[i]; } /* Compare the two curves by computed the distance square between */ /* corresponding parameter values. */ DiffCrv = SymbCrvSub(Crv, CubicBzr); DistSqrCrv = SymbCrvDotProd(DiffCrv, DiffCrv); CagdCrvFree(DiffCrv); Points = DistSqrCrv -> Points; if (IsNotRational) { CagdRType *Dist = Points[1]; for (i = DistSqrCrv -> Order - 1; i >= 0; i--) { if (*Dist++ > Tol2) { ApproxOK = FALSE; break; } } } else { CagdRType *Dist0 = Points[0], *Dist1 = Points[1]; for (i = DistSqrCrv -> Order - 1; i >= 0; i--) { if (*Dist1++ / *Dist0++ > Tol2) { ApproxOK = FALSE; break; } } } CagdCrvFree(DistSqrCrv); if (ApproxOK) { if (NewCrv) CagdCrvFree(Crv); return CubicBzr; } else { CagdCrvStruct *Crv1 = BzrCrvSubdivAtParam(Crv, 0.5), *Crv2 = Crv1 -> Pnext, *Crv1Approx = BzrApproxBzrCrvAsCubicPoly(Crv1, Tol2), *Crv2Approx = BzrApproxBzrCrvAsCubicPoly(Crv2, Tol2); CagdCrvFree(Crv1); CagdCrvFree(Crv2); for (Crv1 = Crv1Approx; Crv1 -> Pnext != NULL; Crv1 = Crv1 ->Pnext); Crv1 -> Pnext = Crv2Approx; if (NewCrv) CagdCrvFree(Crv); return Crv1Approx; } }