/****************************************************************************** * Bzr_Sym.c - Bezier symbolic computation. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Sep. 92. * ******************************************************************************/ #include "symb_loc.h" #define PINDEX_IJ(i,j) P[(i) + (j) * (m+1)] #define QINDEX_IJ(i,j) Q[(i) + (j) * m] static void BzrSrfFactorBilinearAux0(CagdRType *P, CagdRType *Q, CagdRType *A, int m, int n); static void BzrSrfFactorBilinearAux3(CagdRType *P, CagdRType *Q, CagdRType *A, int m, int n); static void BzrSrfFactorBilinearAux1(CagdRType *P, CagdRType *Q, CagdRType *A, int m, int n); static void BzrSrfFactorBilinearAux2(CagdRType *P, CagdRType *Q, CagdRType *A, int m, int n); /***************************************************************************** * 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 * * * SEE ALSO: M * BzrCrvMultPtsVecs M * * * KEYWORDS: M * BzrCrvMult, product M *****************************************************************************/ CagdCrvStruct *BzrCrvMult(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2) { STATIC_DATA int CpPtsLen = 0; STATIC_DATA CagdRType *CpPts1 = NULL, *CpPts2 = NULL; CagdBType IsNotRational; int i, j, ij, l, MaxCoord, IsScalar, ProdOrder, DupCrvs, Order1 = Crv1 -> Order, Order2 = Crv2 -> Order, Degree1 = Order1 - 1, Degree2 = Order2 - 1, ProdDegree = Degree1 + Degree2; CagdCrvStruct *ProdCrv; CagdRType **PPoints, **Points1, **Points2; if (!CAGD_IS_BEZIER_CRV(Crv1) || !CAGD_IS_BEZIER_CRV(Crv2)) { SYMB_FATAL_ERROR(SYMB_ERR_BZR_CRV_EXPECT); return NULL; } if (Crv1 -> PType != Crv2 -> PType) { Crv1 = CagdCrvCopy(Crv1); Crv2 = CagdCrvCopy(Crv2); if (!CagdMakeCrvsCompatible(&Crv1, &Crv2, FALSE, FALSE)) { SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT); return NULL; } DupCrvs = TRUE; } else DupCrvs = FALSE; ProdCrv = BzrCrvNew(ProdOrder = Order1 + Order2 - 1, Crv1 -> PType); IsNotRational = !CAGD_IS_RATIONAL_CRV(ProdCrv); MaxCoord = CAGD_NUM_OF_PT_COORD(ProdCrv -> PType); IsScalar = IsNotRational && MaxCoord == 1; PPoints = ProdCrv -> Points; Points1 = Crv1 -> Points; Points2 = Crv2 -> Points; for (l = IsNotRational; l <= MaxCoord; l++) ZAP_MEM(PPoints[l], sizeof(CagdRType) * ProdOrder); /* Allocate temporary data, if necessary. */ if (CpPtsLen < MAX(Order1, Order2)) { CpPtsLen = MAX(Order1, Order2) * 2; if (CpPts1 != NULL) IritFree(CpPts1); if (CpPts2 != NULL) IritFree(CpPts2); CpPts1 = (CagdRType *) IritMalloc(sizeof(CagdRType) * CpPtsLen); CpPts2 = (CagdRType *) IritMalloc(sizeof(CagdRType) * CpPtsLen); } if (ProdOrder < CAGD_MAX_BEZIER_CACHE_ORDER) { /* Use cache. */ for (l = IsNotRational; l <= MaxCoord; l++) { CagdRType *PPts = PPoints[l]; /* Place the combinatorial coefficients with the control points. */ for (i = 0; i < Order1; i++) CpPts1[i] = Points1[l][i] * CagdIChooseKTable[Degree1][i]; for (i = 0; i < Order2; i++) CpPts2[i] = Points2[l][i] * CagdIChooseKTable[Degree2][i]; /* Do the convolution. */ for (i = 0; i < Order1; i++) { for (j = 0, ij = i; j < Order2; j++, ij++) PPts[ij] += CpPts1[i] * CpPts2[j]; } /* Update the denominator combinatorial coefficient. */ for (i = 0; i < ProdOrder; i++) PPts[i] /= CagdIChooseKTable[ProdDegree][i]; } } else { for (l = IsNotRational; l <= MaxCoord; l++) { CagdRType *PPts = PPoints[l]; /* Place the combinatorial coefficients with the control points. */ for (i = 0; i < Order1; i++) CpPts1[i] = Points1[l][i] * CagdIChooseK(i, Degree1); for (i = 0; i < Order2; i++) CpPts2[i] = Points2[l][i] * CagdIChooseK(i, Degree2); /* Do the convolution. */ for (i = 0; i < Order1; i++) { for (j = 0, ij = i; j < Order2; j++, ij++) PPts[ij] += CpPts1[i] * CpPts2[j]; } /* Update the denominator combinatorial coefficient. */ for (i = 0; i < ProdOrder; i++) PPts[i] /= CagdIChooseK(i, ProdDegree); } } if (DupCrvs) { CagdCrvFree(Crv1); CagdCrvFree(Crv2); } return ProdCrv; } /***************************************************************************** * DESCRIPTION: M * Given two Bezier scalar curves as vectors Ptsi of orders Orderi, multiply M * them. M * * * PARAMETERS: M * Pts1: First vector of scalars of first Bezier curve. M * Order1: Order of first Bezier curve. M * Pts2: Second vector of scalars of second Bezier curve. M * Order2: Order of second Bezier curve. M * ProdPts: Result vector of scalars of product Bezier curve. Result M * vector is of length Order1+Order2-1. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * BzrCrvMult M * * * KEYWORDS: M * BzrCrvMultPtsVecs, product M *****************************************************************************/ void BzrCrvMultPtsVecs(CagdRType *Pts1, int Order1, CagdRType *Pts2, int Order2, CagdRType *ProdPts) { STATIC_DATA int CpPtsLen = 0; STATIC_DATA CagdRType *CpPts1 = NULL, *CpPts2 = NULL; int i, j, ij, Degree1 = Order1 - 1, Degree2 = Order2 - 1, ProdOrder = Order1 + Order2 - 1, ProdDegree = ProdOrder - 1; ZAP_MEM(ProdPts, sizeof(CagdRType) * ProdOrder); /* Allocate temporary data, if necessary. */ if (CpPtsLen < MAX(Order1, Order2)) { CpPtsLen = MAX(Order1, Order2) * 2; if (CpPts1 != NULL) IritFree(CpPts1); if (CpPts2 != NULL) IritFree(CpPts2); CpPts1 = (CagdRType *) IritMalloc(sizeof(CagdRType) * CpPtsLen); CpPts2 = (CagdRType *) IritMalloc(sizeof(CagdRType) * CpPtsLen); } if (ProdOrder < CAGD_MAX_BEZIER_CACHE_ORDER) { /* Use cache. */ /* Place the combinatorial coefficients with the control points. */ for (i = 0; i < Order1; i++) CpPts1[i] = Pts1[i] * CagdIChooseKTable[Degree1][i]; for (i = 0; i < Order2; i++) CpPts2[i] = Pts2[i] * CagdIChooseKTable[Degree2][i]; /* Do the convolution. */ for (i = 0; i < Order1; i++) { for (j = 0, ij = i; j < Order2; j++, ij++) ProdPts[ij] += CpPts1[i] * CpPts2[j]; } /* Update the denominator combinatorial coefficient. */ for (i = 0; i < ProdOrder; i++) ProdPts[i] /= CagdIChooseKTable[ProdDegree][i]; } else { /* Place the combinatorial coefficients with the control points. */ for (i = 0; i < Order1; i++) CpPts1[i] = Pts1[i] * CagdIChooseK(i, Degree1); for (i = 0; i < Order2; i++) CpPts2[i] = Pts2[i] * CagdIChooseK(i, Degree2); /* Do the convolution. */ for (i = 0; i < Order1; i++) { for (j = 0, ij = i; j < Order2; j++, ij++) ProdPts[ij] += CpPts1[i] * CpPts2[j]; } /* Update the denominator combinatorial coefficient. */ for (i = 0; i < ProdOrder; i++) ProdPts[i] /= CagdIChooseK(i, ProdDegree); } } /***************************************************************************** * 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) { STATIC_DATA int CpPtsLen = 0; STATIC_DATA CagdRType *CpPts1 = NULL, *CpPts2 = NULL; CagdBType IsNotRational; int i, j, k, l, m, MaxCoord, UDegree, VDegree, UOrder, VOrder, Size, DupSrfs, 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 **PPoints, **Points1, **Points2; if (!CAGD_IS_BEZIER_SRF(Srf1) || !CAGD_IS_BEZIER_SRF(Srf2)) { SYMB_FATAL_ERROR(SYMB_ERR_BZR_SRF_EXPECT); return NULL; } /* Allocate temporary data, if necessary. */ if (CpPtsLen < MAX(UOrder1, UOrder2) * MAX(VOrder1, VOrder2)) { CpPtsLen = MAX(UOrder1, UOrder2) * MAX(VOrder1, VOrder2) * 2; if (CpPts1 != NULL) IritFree(CpPts1); if (CpPts2 != NULL) IritFree(CpPts2); CpPts1 = (CagdRType *) IritMalloc(sizeof(CagdRType) * CpPtsLen); CpPts2 = (CagdRType *) IritMalloc(sizeof(CagdRType) * CpPtsLen); } if (Srf1 -> PType != Srf2 -> PType) { Srf1 = CagdSrfCopy(Srf1); Srf2 = CagdSrfCopy(Srf2); if (!CagdMakeSrfsCompatible(&Srf1, &Srf2, FALSE, FALSE, FALSE, FALSE)) { SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT); return NULL; } DupSrfs = TRUE; } else DupSrfs = FALSE; 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); Size = UOrder * VOrder; UDegree = UOrder - 1; VDegree = VOrder - 1; PPoints = ProdSrf -> Points; Points1 = Srf1 -> Points; Points2 = Srf2 -> Points; for (k = IsNotRational; k <= MaxCoord; k++) ZAP_MEM(PPoints[k], sizeof(CagdRType) * Size); if (UOrder < CAGD_MAX_BEZIER_CACHE_ORDER && VOrder < CAGD_MAX_BEZIER_CACHE_ORDER) { for (k = IsNotRational; k <= MaxCoord; k++) { CagdRType *p, *p1, *p2, *PPts = PPoints[k]; /* Place the combinatorial coefficients with the control points. */ for (p1 = Points1[k], p = CpPts1, j = l = 0; j < VOrder1; j++) { for (i = 0; i < UOrder1; i++, l++) *p++ = *p1++ * CagdIChooseKTable[VDegree1][j] * CagdIChooseKTable[UDegree1][i]; } for (p2 = Points2[k], p = CpPts2, j = l = 0; j < VOrder2; j++) { for (i = 0; i < UOrder2; i++, l++) *p++ = *p2++ * CagdIChooseKTable[VDegree2][j] * CagdIChooseKTable[UDegree2][i]; } /* Do the convolution. */ for (j = 0; j < VOrder1; j++) { for (i = 0; i < UOrder1; i++) { CagdRType r1 = CpPts1[CAGD_MESH_UV(Srf1, i, j)]; p2 = CpPts2; for (m = 0; m < VOrder2; m++) { p = &PPts[CAGD_MESH_UV(ProdSrf, i, j + m)]; for (l = 0; l++ < UOrder2; ) *p++ += r1 * *p2++; } } } /* Update the denominator combinatorial coefficient. */ for (j = 0; j < VOrder; j++) { p = &PPts[CAGD_MESH_UV(ProdSrf, 0, j)]; for (i = 0; i < UOrder; i++) *p++ /= CagdIChooseKTable[UDegree][i] * CagdIChooseKTable[VDegree][j]; } } } else { /* The original product - easier to follow but not so optimized. */ /* */ /* int il, jm; */ /* */ /* for (i = 0; i < UOrder1; i++) { */ /* for (j = 0; j < VOrder1; j++) { */ /* for (l = 0, il = i; l < UOrder2; l++, il++) { */ /* for (m = 0, jm = j; m < VOrder2; m++, jm++) { */ /* int Index = CAGD_MESH_UV(ProdSrf, il, jm), */ /* Index1 = CAGD_MESH_UV(Srf1, i, j), */ /* Index2 = CAGD_MESH_UV(Srf2, l, m); */ /* CagdRType */ /* Coef = CagdIChooseK(i, UDegree1) * */ /* CagdIChooseK(l, UDegree2) * */ /* CagdIChooseK(j, VDegree1) * */ /* CagdIChooseK(m, VDegree2) / */ /* (CagdIChooseK(il, UDegree) * */ /* CagdIChooseK(jm, VDegree)); */ /* */ /* for (k = IsNotRational; k <= MaxCoord; k++) */ /* PPoints[k][Index] += Coef * */ /* Points1[k][Index1] * Points2[k][Index2]; */ /* } */ /* } */ /* } */ /* } */ for (k = IsNotRational; k <= MaxCoord; k++) { CagdRType *p, *p1, *p2, *PPts = PPoints[k]; /* Place the combinatorial coefficients with the control points. */ for (p1 = Points1[k], p = CpPts1, j = l = 0; j < VOrder1; j++) { for (i = 0; i < UOrder1; i++, l++) *p++ = *p1++ * CagdIChooseK(j, VDegree1) * CagdIChooseK(i, UDegree1); } for (p2 = Points2[k], p = CpPts2, j = l = 0; j < VOrder2; j++) { for (i = 0; i < UOrder2; i++, l++) *p++ = *p2++ * CagdIChooseK(j, VDegree2) * CagdIChooseK(i, UDegree2); } /* Do the convolution. */ for (j = 0; j < VOrder1; j++) { for (i = 0; i < UOrder1; i++) { CagdRType r1 = CpPts1[CAGD_MESH_UV(Srf1, i, j)]; p2 = CpPts2; for (m = 0; m < VOrder2; m++) { p = &PPts[CAGD_MESH_UV(ProdSrf, i, j + m)]; for (l = 0; l++ < UOrder2; ) *p++ += r1 * *p2++; } } } /* Update the denominator combinatorial coefficient. */ for (j = 0; j < VOrder; j++) { p = &PPts[CAGD_MESH_UV(ProdSrf, 0, j)]; for (i = 0; i < UOrder; i++) *p++ /= CagdIChooseK(i, UDegree) * CagdIChooseK(j, VDegree); } } } if (DupSrfs) { CagdSrfFree(Srf1); CagdSrfFree(Srf2); } return ProdSrf; } /***************************************************************************** * DESCRIPTION: M * Factors out a (u - v) term from a scalar surface, assuming it has one. M * S(u,v) / {A[0] (1-u)(1-v) + A[1] u(1-v) + A[2] (1-u)v + A[3] uv} V * V * If S(P) = Bilinear(A) * S(Q), then V * A[0] (m-i) (n-j) Q[i][j] + A[1] i (n-j) Q[i-1][j] + V * A[2] (m-i) j Q[i][j-1] + A[3] i j Q[i-1][j-1] = m n P[i][j] V * * * PARAMETERS: M * Srf: To factor out a bilinear term from. M * A: Four coeficients of the scalar bilinear. M * * * RETURN VALUE: M * CagdSrfStruct *: Srf / (u - v). M * * * SEE ALSO: M * BspSrfFactorUMinusV M * * * KEYWORDS: M * BzrSrfFactorBilinear M *****************************************************************************/ CagdSrfStruct *BzrSrfFactorBilinear(CagdSrfStruct *Srf, CagdRType *A) { int m = Srf -> ULength - 1, n = Srf -> VLength - 1; CagdSrfStruct *FctrSrf = BzrSrfNew(m, n, Srf -> PType); CagdRType *P = Srf -> Points[1], *Q = FctrSrf -> Points[1]; int Selected; /* select the largest absolute value of A[i] */ if (FABS(A[0]) >= FABS(A[1])) { if (FABS(A[0]) >= FABS(A[2])) { if (FABS(A[0]) >= FABS(A[3])) Selected = 0; else Selected = 3; } else { if (FABS(A[2]) >= FABS(A[3])) Selected = 2; else Selected = 3; } } else { if (FABS(A[1]) >= FABS(A[2])) { if (FABS(A[1]) >= FABS(A[3])) Selected = 1; else Selected = 3; } else { if (FABS(A[2]) >= FABS(A[3])) Selected = 2; else Selected = 3; } } if (FABS(A[Selected]) < IRIT_UEPS) { SYMB_FATAL_ERROR(SYMB_ERR_DIV_ZERO); } switch(Selected) { case 0: BzrSrfFactorBilinearAux0(P, Q, A, m, n); break; case 1: BzrSrfFactorBilinearAux1(P, Q, A, m, n); break; case 2: BzrSrfFactorBilinearAux2(P, Q, A, m, n); break; case 3: BzrSrfFactorBilinearAux3(P, Q, A, m, n); break; default: break; } return FctrSrf; } /***************************************************************************** * AUXILIARY: * * Auxiliary function of BzrSrfFactorBilinear. * *****************************************************************************/ static void BzrSrfFactorBilinearAux0(CagdRType *P, CagdRType *Q, CagdRType *A, int m, int n) { int i, j; QINDEX_IJ(0,0) = PINDEX_IJ(0,0) / A[0]; for (i = 1; i < m; i++) QINDEX_IJ(i,0) = (m * PINDEX_IJ(i,0) - A[1] * i * QINDEX_IJ(i - 1,0)) / (A[0] * (m - i)); for (j = 1; j < n; j++) QINDEX_IJ(0,j) = (n * PINDEX_IJ(0,j) - A[2] * j * QINDEX_IJ(0,j - 1)) / (A[0] * (n - j)); for (i = 1; i < m; i++) for (j = 1; j < n; j++) QINDEX_IJ(i,j) = ((m * n) * PINDEX_IJ(i,j) - (A[1] * i * (n - j)) * QINDEX_IJ(i - 1, j) - (A[2] * (m - i) * j) * QINDEX_IJ(i, j - 1) - (A[3] * i * j) * QINDEX_IJ(i - 1, j - 1)) / (A[0] * (m - i) * (n - j)); } /***************************************************************************** * AUXILIARY: * * Auxiliary function of BzrSrfFactorBilinear. * *****************************************************************************/ static void BzrSrfFactorBilinearAux3(CagdRType *P, CagdRType *Q, CagdRType *A, int m, int n) { int i, j; QINDEX_IJ(m-1,n-1) = PINDEX_IJ(m,n) / A[3]; for (i = m-1; i; i--) QINDEX_IJ(i - 1,n-1) = (m * PINDEX_IJ(i,n) - A[2] * (m - i) * QINDEX_IJ(i,n-1)) / (A[3] * i); for (j = n-1; j; j--) QINDEX_IJ(m-1,j - 1) = (n * PINDEX_IJ(m,j) - A[1] * (n - j) * QINDEX_IJ(m-1,j)) / (A[3] * j); for (i = m-1; i; i--) for (j = n-1; j; j--) QINDEX_IJ(i - 1,j - 1) = ((m * n) * PINDEX_IJ(i,j) - (A[0] * (m - i) * (n - j)) * QINDEX_IJ(i,j) - (A[1] * i * (n - j)) * QINDEX_IJ(i - 1,j) - (A[2] * (m - i) * j) * QINDEX_IJ(i,j - 1)) / (A[3] * i * j); } /***************************************************************************** * AUXILIARY: * * Auxiliary function of BzrSrfFactorBilinear. * *****************************************************************************/ static void BzrSrfFactorBilinearAux1(CagdRType *P, CagdRType *Q, CagdRType *A, int m, int n) { int i, j; QINDEX_IJ(m-1, 0) = PINDEX_IJ(m,0) / A[1]; for (i = m-1; i; i--) QINDEX_IJ(i - 1,0) = (m * PINDEX_IJ(i,0) - A[0] * (m - i) * QINDEX_IJ(i,0)) / (A[1] * i); for (j = 1; j < n; j++) QINDEX_IJ(m-1,j) = (n * PINDEX_IJ(m,j) - A[3] * j * QINDEX_IJ(m-1,j - 1)) / (A[1] * (n - j)); for (i = m-1; i; i--) for (j = 1; j < n; j++) QINDEX_IJ(i - 1,j) = ((m * n) * PINDEX_IJ(i,j) - (A[0] * (m - i) * (n - j)) * QINDEX_IJ(i, j) - (A[2] * (m - i) * j) * QINDEX_IJ(i, j - 1) - (A[3] * i * j) * QINDEX_IJ(i - 1, j - 1)) / (A[1] * i * (n - j)); } /***************************************************************************** * AUXILIARY: * * Auxiliary function of BzrSrfFactorBilinear. * *****************************************************************************/ static void BzrSrfFactorBilinearAux2(CagdRType *P, CagdRType *Q, CagdRType *A, int m, int n) { int i, j; QINDEX_IJ(0, n-1) = PINDEX_IJ(0,n) / A[2]; for (i = 1; i < m; i++) QINDEX_IJ(i,n-1) = (m * PINDEX_IJ(i,n) - A[3] * i * QINDEX_IJ(i - 1,n-1)) / (A[2] * (m - i)); for (j = n-1; j; j--) QINDEX_IJ(0,j - 1) = (n * PINDEX_IJ(0,j) - A[0] * (n - j) * QINDEX_IJ(0,j)) / (A[2] * j); for (i = 1; i < m; i++) for (j = n-1; j; j--) QINDEX_IJ(i,j - 1) = ((m * n) * PINDEX_IJ(i,j) - (A[0] * (m - i) * (n - j)) * QINDEX_IJ(i, j) - (A[1] * i * (n - j)) * QINDEX_IJ(i - 1, j) - (A[3] * i * j) * QINDEX_IJ(i - 1, j - 1)) / (A[2] * (m - i) * j); } /***************************************************************************** * DESCRIPTION: M * Factors out a (u - v) term from a scalar surface, assuming it has one. M * * * PARAMETERS: M * Srf: To factor out a (u - v) term from. M * * * RETURN VALUE: M * CagdSrfStruct *: Srf / (u - v). M * * * SEE ALSO: M * BspSrfFactorUMinusV, BzrSrfFactorBilinear M * * * KEYWORDS: M * BzrSrfFactorUMinusV M *****************************************************************************/ CagdSrfStruct *BzrSrfFactorUMinusV(CagdSrfStruct *Srf) { int i, j, m = Srf -> ULength - 1, n = Srf -> VLength - 1; CagdSrfStruct *FctrSrf = BzrSrfNew(m, n, Srf -> PType); CagdRType *P = Srf -> Points[1], *Q = FctrSrf -> Points[1]; for (j = 0; j < n; j++) Q[j * m] = -n * P[(j + 1) * (m + 1)] / (j + 1); for (i = 1; i < m; i++) for (j = 0; j < n; j++) Q[i + j * m] = ((i * (n - j - 1)) * Q[i - 1 + (j + 1) * m] - (m * n) * P[i +(j + 1) * (m + 1)]) / ((m - i) * (j + 1)); return FctrSrf; } /***************************************************************************** * 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 { SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_EXPECTED); return NULL; } 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); return NULL; } 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 { SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_EXPECTED); return NULL; } 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); return NULL; } 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 * polynomials. 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 polynomials M * approximation is computed for Crv. M * M * In case 3 a list of polynomial cubic curves is returned. Tol is then M * used for the distance tolerance error measure for the approximation. M * If the total length of control polygon is (approximately) more than M * MaxLen, 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 * * * RETURN VALUE: M * CagdCrvStruct *: A list of cubic Bezier polynomials approximating Crv. M * * * SEE ALSO: M * BzrApproxBzrCrvAsQuadraticPoly, BzrApproxBzrCrvAsQuadratics, M * BzrApproxBzrCrvAsCubicPoly M * * * KEYWORDS: M * BzrApproxBzrCrvAsCubics, conversion, approximation M *****************************************************************************/ CagdCrvStruct *BzrApproxBzrCrvAsCubics(CagdCrvStruct *Crv, CagdRType Tol, CagdRType MaxLen) { CagdBType NewCrv = FALSE; CagdCrvStruct *TCrv, *CubicCrvs, *CubicCrvsMaxLen = NULL; if (CAGD_IS_PERIODIC_CRV(Crv)) { NewCrv = TRUE; Crv = CnvrtPeriodic2FloatCrv(Crv); } if (CAGD_IS_BSPLINE_CRV(Crv)) { if (!BspCrvHasOpenEC(Crv)) { TCrv = BspCrvOpenEnd(Crv); if (NewCrv) CagdCrvFree(Crv); Crv = TCrv; NewCrv = TRUE; } TCrv = CnvrtBspline2BezierCrv(Crv); if (NewCrv) CagdCrvFree(Crv); for (Crv = TCrv, CubicCrvs = NULL; Crv != NULL; Crv = Crv -> Pnext) { CubicCrvs = CagdListAppend(CubicCrvs, BzrApproxBzrCrvAsCubics(Crv, Tol, MaxLen)); } CagdCrvFreeList(TCrv); return CubicCrvs; } /* Only Bezier Curves should get to here. */ if (!CAGD_IS_BEZIER_CRV(Crv)) { SYMB_FATAL_ERROR(SYMB_ERR_BZR_CRV_EXPECT); return NULL; } 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); return NULL; } CubicCrvs = BzrApproxBzrCrvAsCubicPoly(Crv, Tol * Tol); } for (TCrv = CubicCrvs; TCrv != NULL; TCrv = TCrv -> Pnext) { CagdCrvStruct *MaxLenCrv = CagdLimitCrvArcLen(TCrv, MaxLen); CubicCrvsMaxLen = CagdListAppend(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 (squared). M * * * RETURN VALUE: M * CagdCrvStruct *: List of cubic Bezier curves approximating Crv. M * * * SEE ALSO: M * BzrApproxBzrCrvAsQuadraticPoly, BzrApproxBzrCrvAsQuadratics, M * BzrApproxBzrCrvAsCubics 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 = MIN(CAGD_NUM_OF_PT_COORD(Crv -> PType), 3); 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 computing 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); CagdCrvFree(CubicBzr); if (NewCrv) CagdCrvFree(Crv); return CagdListAppend(Crv1Approx, Crv2Approx);; } } /***************************************************************************** * DESCRIPTION: M * Given a Bezier curve - convert it to (possibly) piecewise quadratic M * polynomials. M * If the curve is M * 1. A quadratic - a copy if it is returned. M * 2. Lower than quadratic - a degree raised (to a quadratic) curve is M * returned. M * 3. Higher than quadratic - a C^1 continuous piecewise quadratic M * polynomial approximation is computed for Crv. M * M * In case 3, a list of polynomial quadratic curves is returned. Tol is M * then used for the distance tolerance error measure for the approximation. M * If the total length of control polygon is (approximately) more than M * MaxLen, the curve is subdivided until this is not the case. M * * * PARAMETERS: M * Crv: To approximate using quadratic Bezier polynomials. M * Tol: Accuracy control. M * MaxLen: Maximum arc length of curve. M * * * RETURN VALUE: M * CagdCrvStruct *: A list of quadratic Bezier polynomials approximating M * Crv. M * * * SEE ALSO: M * BzrApproxBzrCrvAsQuadraticPoly, BzrApproxBzrCrvAsCubics, M * BzrApproxBzrCrvAsCubicPoly M * * * KEYWORDS: M * BzrApproxBzrCrvAsQuadratics, conversion, approximation M *****************************************************************************/ CagdCrvStruct *BzrApproxBzrCrvAsQuadratics(CagdCrvStruct *Crv, CagdRType Tol, CagdRType MaxLen) { CagdBType NewCrv = FALSE; CagdCrvStruct *TCrv, *QuadraticCrvs, *QuadraticCrvsMaxLen = NULL; if (CAGD_IS_PERIODIC_CRV(Crv)) { NewCrv = TRUE; Crv = CnvrtPeriodic2FloatCrv(Crv); } if (CAGD_IS_BSPLINE_CRV(Crv)) { if (!BspCrvHasOpenEC(Crv)) { TCrv = BspCrvOpenEnd(Crv); if (NewCrv) CagdCrvFree(Crv); Crv = TCrv; NewCrv = TRUE; } TCrv = CnvrtBspline2BezierCrv(Crv); if (NewCrv) CagdCrvFree(Crv); for (Crv = TCrv, QuadraticCrvs = NULL; Crv != NULL; Crv = Crv -> Pnext) { QuadraticCrvs = CagdListAppend(QuadraticCrvs, BzrApproxBzrCrvAsQuadratics(Crv, Tol, MaxLen)); } CagdCrvFreeList(TCrv); return QuadraticCrvs; } /* Only Bezier Curves should get to here. */ if (!CAGD_IS_BEZIER_CRV(Crv)) { SYMB_FATAL_ERROR(SYMB_ERR_BZR_CRV_EXPECT); return NULL; } switch (Crv -> Order) { case 2: QuadraticCrvs = BzrCrvDegreeRaiseN(Crv, 3); break; case 3: QuadraticCrvs = CagdCrvCopy(Crv); break; default: if (Crv -> Order < 3) { SYMB_FATAL_ERROR(SYMB_ERR_OUT_OF_RANGE); return NULL; } QuadraticCrvs = BzrApproxBzrCrvAsQuadraticPoly(Crv, Tol * Tol); } for (TCrv = QuadraticCrvs; TCrv != NULL; TCrv = TCrv -> Pnext) { CagdCrvStruct *MaxLenCrv = CagdLimitCrvArcLen(TCrv, MaxLen); QuadraticCrvsMaxLen = CagdListAppend(QuadraticCrvsMaxLen, MaxLenCrv); } if (NewCrv) CagdCrvFree(Crv); CagdCrvFreeList(QuadraticCrvs); return QuadraticCrvsMaxLen; } /***************************************************************************** * 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, using pairs of quadratics. 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 (squared). M * * * RETURN VALUE: M * CagdCrvStruct *: List of cubic Bezier curves approximating Crv. M * * * SEE ALSO: M * BzrApproxBzrCrvAsQuadratics, BzrApproxBzrCrvAsCubics, M * BzrApproxBzrCrvAsCubicPoly M * * * KEYWORDS: M * BzrApproxBzrCrvAsQuadraticPoly, approximation, conversion M *****************************************************************************/ CagdCrvStruct *BzrApproxBzrCrvAsQuadraticPoly(CagdCrvStruct *Crv, CagdRType Tol2) { CagdBType NewCrv = FALSE, ApproxOK = TRUE, IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); int i, Order = Crv -> Order, MaxCoord = MIN(CAGD_NUM_OF_PT_COORD(Crv -> PType), 3); CagdPointType PType = Crv -> PType, QuadraticPType = CAGD_MAKE_PT_TYPE(FALSE, MaxCoord); CagdCrvStruct *DistSqrCrv, *DiffCrv, *QuadBsp, *QuadraticBzr1 = BzrCrvNew(3, QuadraticPType), *QuadraticBzr2 = BzrCrvNew(3, QuadraticPType); CagdRType E3Pt1[3], E3Pt2[3], E3Pt3[3], E3Pt4[3], Tan1[3], Tan2[3], TMin, TMax, **QuadraticPoints1 = QuadraticBzr1 -> Points, **QuadraticPoints2 = QuadraticBzr2 -> 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); /* Tangents at the end of the two points must match. */ PT_SUB(Tan1, E3Pt2, E3Pt1); PT_SUB(Tan2, E3Pt4, E3Pt3); /* Divide by 2 (quadratic degree) and by 2 since we have two pieces. */ PT_SCALE(Tan1, (Order - 1.0) / 4.0); PT_SCALE(Tan2, (Order - 1.0) / 4.0); for (i = 0; i < MaxCoord; i++) { /* End of the two points must match. */ QuadraticPoints1[i + 1][0] = E3Pt1[i]; QuadraticPoints2[i + 1][2] = E3Pt4[i]; /* Middle points are fully governed from the C^1 continuity. */ QuadraticPoints1[i + 1][1] = E3Pt1[i] + Tan1[i]; QuadraticPoints2[i + 1][1] = E3Pt4[i] - Tan2[i]; /* End point of first cyrve == start point of second curve - middle */ /* point due to the C^1 continuity. constraint. */ QuadraticPoints1[i + 1][2] = QuadraticPoints2[i + 1][0] = (QuadraticPoints1[i + 1][1] + QuadraticPoints2[i + 1][1]) * 0.5; } /* Compare the two curves by computing the distance square between */ /* corresponding parameter values. */ QuadBsp = CagdMergeCrvCrv(QuadraticBzr1, QuadraticBzr2, FALSE); CagdCrvDomain(Crv, &TMin, &TMax); BspKnotAffineTransOrder2(QuadBsp -> KnotVector, QuadBsp -> Order, QuadBsp -> Length + QuadBsp -> Order, TMin, TMax); DiffCrv = SymbCrvSub(Crv, QuadBsp); 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); CagdCrvFree(QuadBsp); if (ApproxOK) { if (NewCrv) CagdCrvFree(Crv); QuadraticBzr1 -> Pnext = QuadraticBzr2; return QuadraticBzr1; } else { CagdCrvStruct *Crv1 = BzrCrvSubdivAtParam(Crv, 0.5), *Crv2 = Crv1 -> Pnext, *Crv1Approx = BzrApproxBzrCrvAsQuadraticPoly(Crv1, Tol2), *Crv2Approx = BzrApproxBzrCrvAsQuadraticPoly(Crv2, Tol2); CagdCrvFree(Crv1); CagdCrvFree(Crv2); CagdCrvFree(QuadraticBzr1); CagdCrvFree(QuadraticBzr2); if (NewCrv) CagdCrvFree(Crv); return CagdListAppend(Crv1Approx, Crv2Approx); } }