/****************************************************************************** * Crv_Skel.c - computation of curve/surface skeleton approximation. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, August 96 * ******************************************************************************/ #include "symb_loc.h" #include "user_lib.h" #include "iritprsr.h" #include "allocate.h" #include "geomat3d.h" #define CONTOUR_EPS 1e-8 /* Level above zero to actually contour. */ #define DIAG_POINT_EPS 1e-3 /* Off diagonal to be considered diagonal. */ #define MAX_DOMAIN 10 /* 10 times the bounding box. */ #define MAX_NUMER_ITER 3 /* Number of numerical iterations. */ #define MAX_ERROR_ALLOWED 1e-6 /* In mid point computation. */ #define SRF_SCALE_FACTOR 1 /* Scaling factor for better contouring. */ static int NumerMarchMidPoint(CagdCrvStruct *Crv1, CagdRType *t1, CagdCrvStruct *Crv2, CagdRType *t2); static CagdRType *ComputeInterMidPoint(CagdCrvStruct *Crv1, CagdRType t1, CagdCrvStruct *Crv2, CagdRType t2); static CagdSrfStruct *SymbCrvBisectorsSrf3D(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2, CagdSrfStruct *DSrf1, CagdSrfStruct *DSrf2); /***************************************************************************** * DESCRIPTION: M * Computes the skeleton curves (bisectors) of a given curve or two. M * If Crv contains a list of two curves the bisector between the two curves M * is computed. Otherwise, Crv self bisectors are computed. M * * * PARAMETERS: M * Crv: Either one or two curves to compute bisectors for. M * Assumes E2 curves. M * UseNrmlTan: If 0, tangent fields are used to compute bisector surface. M * If 1, normal fields instead of tangents are used to M * compute the bisector surface function. If 2, a blend of M * the normal field and tangent field computation is used. M * Tolerance: Accuracy of computation. M * NumerImprove: If TRUE, a numerical improvment stage is applied. M * SameNormal: If TRUE, the bisector should be oriented for inner or M * outer side of the curves, with respect to their normals. M * * * RETURN VALUE: M * CagdCrvStruct *: A list of piecewise linear curves approximating the M * bisectors of Crv. M * * * SEE ALSO: M * SymbCrvCnvxHull, SymbCrvDiameter, SymbCrvBisectorsSrf M * * * KEYWORDS: M * SymbCrvBisectors, bisectors, skeleton M *****************************************************************************/ CagdCrvStruct *SymbCrvBisectors(CagdCrvStruct *Crv, int UseNrmlTan, RealType Tolerance, CagdBType NumerImprove, CagdBType SameNormal) { static PlaneType Plane = { 1.0, 0.0, 0.0, CONTOUR_EPS }; /* A scalar srf - only X. */ CagdBType SelfBisect = Crv -> Pnext == NULL; CagdRType TMin1, TMax1, TMin2, TMax2, XMax, YMax; CagdPType Trans; CagdSrfStruct *BisectSrf = SymbCrvBisectorsSrf(Crv, UseNrmlTan); IPPolygonStruct *Cntrs, *Cntr; CagdCrvStruct *PLCrvList = NULL, *Crv1 = Crv, *Crv2 = Crv -> Pnext ? Crv -> Pnext : Crv; CagdBBoxStruct BBox; CagdCrvDomain(Crv1, &TMin1, &TMax1); CagdCrvDomain(Crv2, &TMin2, &TMax2); PT_RESET(Trans); CagdSrfTransform(BisectSrf, Trans, SRF_SCALE_FACTOR); /* Computes the zero set of the equation as contours. */ Cntrs = UserCntrSrfWithPlane(BisectSrf, Plane, Tolerance); CagdSrfFree(BisectSrf); CagdCrvBBox(Crv, &BBox); XMax = MAX_DOMAIN * MAX(fabs(BBox.Min[0]), fabs(BBox.Max[0])); YMax = MAX_DOMAIN * MAX(fabs(BBox.Min[1]), fabs(BBox.Max[1])); /* Get the st parameters and convert to the bisector's position. */ for (Cntr = Cntrs; Cntr != NULL; Cntr = Cntr -> Pnext) { IPVertexStruct *V, *VPrev = NULL; int i, Len; CagdCrvStruct *PLCrv; CagdRType **Points; /* Filters out the data - no diagonals or below diagonals. */ for (V = Cntr -> PVertex; V != NULL; ) { CagdRType *R, Err; CagdPType Pt1, Pt2, V1, V2, Nrml1, Nrml2; CagdVecStruct *Vec; V -> Coord[1] = BOUND(V -> Coord[1], TMin1, TMax1); V -> Coord[2] = BOUND(V -> Coord[2], TMin2, TMax2); /* Numerically improve the data. */ if (NumerImprove) { NumerMarchMidPoint(Crv1, &V -> Coord[1], Crv2, &V -> Coord[2]); NumerMarchMidPoint(Crv2, &V -> Coord[2], Crv1, &V -> Coord[1]); } Vec = CagdCrvNormalXY(Crv1, V -> Coord[1], TRUE); PT_COPY(Nrml1, Vec -> Vec); Vec = CagdCrvNormalXY(Crv2, V -> Coord[2], TRUE); PT_COPY(Nrml2, Vec -> Vec); R = CagdCrvEval(Crv1, V -> Coord[1]); CagdCoerceToE3(Pt1, &R, -1, Crv1 -> PType); R = CagdCrvEval(Crv2, V -> Coord[2]); CagdCoerceToE3(Pt2, &R, -1, Crv2 -> PType); R = ComputeInterMidPoint(Crv1, V -> Coord[1], Crv2, V -> Coord[2]); PT_SUB(V1, Pt1, R); PT_SUB(V2, Pt2, R); Err = fabs(PT_LENGTH(V1) - PT_LENGTH(V2)); if ((NumerImprove && Err > MAX_ERROR_ALLOWED) || fabs(R[0]) > XMax || fabs(R[1]) > YMax || (SameNormal && DOT_PROD(Nrml1, Nrml2) > 0.0) || (VPrev != NULL && PT_APX_EQ(V -> Coord, VPrev -> Coord)) || (SelfBisect && (APX_EQ_EPS(V -> Coord[1], V -> Coord[2], DIAG_POINT_EPS) || V -> Coord[1] < V -> Coord[2]))) { /* Remove this vertex. */ if (VPrev != NULL) { VPrev -> Pnext = V -> Pnext; IPFreeVertex(V); V = VPrev -> Pnext; } else { /* First vertex in contour. */ Cntr -> PVertex = V -> Pnext; IPFreeVertex(V); V = Cntr -> PVertex; } } else { VPrev = V; V = V -> Pnext; } } Len = IritPrsrVrtxListLen(Cntr -> PVertex); if (Len < 2) continue; PLCrv = BspCrvNew(Len, 2, CAGD_PT_E2_TYPE); Points = PLCrv -> Points; BspKnotUniformOpen(Len, 2, PLCrv -> KnotVector); for (V = Cntr -> PVertex, i = 0; V != NULL; V = V -> Pnext, i++) { CagdRType *R = ComputeInterMidPoint(Crv1, V -> Coord[1], Crv2, V -> Coord[2]); Points[1][i] = R[0]; Points[2][i] = R[1]; } if (PLCrv != NULL) LIST_PUSH(PLCrv, PLCrvList); } IPFreePolygonList(Cntrs); return PLCrvList; } /***************************************************************************** * DESCRIPTION: * * Numerically march along curve 1 in order to improve the result of the * * mid point computation, using Newton Raphson. * * * * PARAMETERS: * * Crv1: First curve of the matching mid point. * * t1: Parameter value of first curve's mid point. * * Crv2 Second curve of the matching mid point. * * t2: Parameter value of second curve's mid point. * * * * RETURN VALUE: * * void TRUE if succesfully improved the mid point, FALSE otherwise. * *****************************************************************************/ static int NumerMarchMidPoint(CagdCrvStruct *Crv1, CagdRType *t1, CagdCrvStruct *Crv2, CagdRType *t2) { int i = 0, Improve = FALSE; CagdPType Pt1, Pt2; CagdRType Error, *R, TMin1, TMax1, NewError, DErrDt1, t1n; CagdCrvDomain(Crv1, &TMin1, &TMax1); R = CagdCrvEval(Crv1, *t1); CagdCoerceToE3(Pt1, &R, -1, Crv1 -> PType); R = CagdCrvEval(Crv2, *t2); CagdCoerceToE3(Pt2, &R, -1, Crv2 -> PType); R = ComputeInterMidPoint(Crv1, *t1, Crv2, *t2); Error = fabs(PT_PT_DIST_SQR(Pt1, R) - PT_PT_DIST_SQR(Pt2, R)); do { t1n = *t1 + IRIT_EPS; if (t1n < TMin1 || t1n > TMax1) break; R = CagdCrvEval(Crv1, t1n); CagdCoerceToE3(Pt1, &R, -1, Crv1 -> PType); R = ComputeInterMidPoint(Crv1, t1n, Crv2, *t2), NewError = fabs(PT_PT_DIST_SQR(Pt1, R) - PT_PT_DIST_SQR(Pt2, R)); /* Compute the change in error as a function of parameter change. */ DErrDt1 = (Error - NewError) / IRIT_EPS; if (DErrDt1 == 0) break; t1n = *t1 + Error / DErrDt1; if (t1n < TMin1 || t1n > TMax1) break; R = CagdCrvEval(Crv1, t1n); CagdCoerceToE3(Pt1, &R, -1, Crv1 -> PType); R = ComputeInterMidPoint(Crv1, t1n, Crv2, *t2); NewError = fabs(PT_PT_DIST_SQR(Pt1, R) - PT_PT_DIST_SQR(Pt2, R)); if (NewError < Error) { Error = NewError; Improve = TRUE; *t1 = t1n; } } while (i++ < MAX_NUMER_ITER && Improve); return Improve; } /***************************************************************************** * DESCRIPTION: * * Computes the intersection point of the normals of the given two points * * on the given two curves. * * * * PARAMETERS: * * Crv1: First curve of the matching mid point. * * t1: Parameter value of first curve's mid point. * * Crv2 Second curve of the matching mid point. * * t2: Parameter value of second curve's mid point. * * * * RETURN VALUE: * * CagdRType *: Point of intersection, statically allocated. * *****************************************************************************/ static CagdRType *ComputeInterMidPoint(CagdCrvStruct *Crv1, CagdRType t1, CagdCrvStruct *Crv2, CagdRType t2) { static CagdPType Inter1; CagdPType Pt1, Pt2, Nrml1, Nrml2, Inter2; CagdRType *R; CagdVecStruct *Vec; R = CagdCrvEval(Crv1, t1); CagdCoerceToE3(Pt1, &R, -1, Crv1 -> PType); R = CagdCrvEval(Crv2, t2); CagdCoerceToE3(Pt2, &R, -1, Crv2 -> PType); Vec = CagdCrvNormalXY(Crv1, t1, TRUE); PT_COPY(Nrml1, Vec -> Vec); Vec = CagdCrvNormalXY(Crv2, t2, TRUE); PT_COPY(Nrml2, Vec -> Vec); CG2PointsFromLineLine(Pt1, Nrml1, Pt2, Nrml2, Inter1, &t1, Inter2, &t2); return Inter1; } /***************************************************************************** * DESCRIPTION: M * Computes the bisector surface definition of a given curve or two. M * If Crv contains a list of two curves the bisector between the two curves M * is computed. Otherwise, Crv self--bisectors are sought. The result is a M * scalar surface whose zero set is the set of bisector(s) of the curves. M * * * PARAMETERS: M * Crv: Either one or two curves to compute bisectors for. Assumes M * E2 curves. M * UseNrmlTan: If 0, tangent fields are used to compute bisector surface. M * If 1, normal fields instead of tangents are used to compute M * the bisector surface function. If 2, a blend of the normal M * field and tangent field computation is used. M * if -1, the curves are assumed three space curve and the real M * bisector surface is returned in such a case. M * * * RETURN VALUE: M * CagdSrfStruct *: A scalar surface whose zero set provides matching M * bisecting points on Crv, if UseNrmlTan >= 0. The real M * bisector surface for three space curves, if UseNrmlTan = -1. M * * * SEE ALSO: M * SymbCrvCnvxHull, SymbCrvDiameter, SymbCrvBisectors, M * SymbCrvBisectorsSrf2, SymbCrvPtBisectorsSrf3D, SymbCrvCrvBisectorSrf3D M * * * KEYWORDS: M * SymbCrvBisectorsSrf, bisectors, skeleton M *****************************************************************************/ CagdSrfStruct *SymbCrvBisectorsSrf(CagdCrvStruct *Crv, int UseNrmlTan) { CagdRType UMin1, UMax1, VMin1, VMax1, UMin2, UMax2, VMin2, VMax2; CagdCrvStruct *Crv1, *Crv2, *DCrv1, *DCrv2; CagdSrfStruct *Res, *Srf1X, *Srf1Y, *Srf2X, *Srf2Y, *DSrf1X, *DSrf1Y, *DSrf2X, *DSrf2Y, *Alpha, *Beta, *TSrf1, *TSrf2, *TSrf3, *TSrf4, *Srf1, *Srf2, *DSrf1, *DSrf2; if (UseNrmlTan == -1) { Crv1 = CagdCoerceCrvTo(Crv, CAGD_PT_E3_TYPE); Crv2 = CagdCoerceCrvTo(Crv -> Pnext ? Crv -> Pnext : Crv, CAGD_PT_E3_TYPE); } else { Crv1 = CagdCoerceCrvTo(Crv, CAGD_PT_E2_TYPE); Crv2 = CagdCoerceCrvTo(Crv -> Pnext ? Crv -> Pnext : Crv, CAGD_PT_E2_TYPE); } if (CAGD_IS_BEZIER_CRV(Crv1)) { CagdCrvStruct *TCrv = CnvrtBezier2BsplineCrv(Crv1); CagdCrvFree(Crv1); Crv1 = TCrv; } if (CAGD_IS_BEZIER_CRV(Crv2)) { CagdCrvStruct *TCrv = CnvrtBezier2BsplineCrv(Crv2); CagdCrvFree(Crv2); Crv2 = TCrv; } DCrv1 = CagdCrvDerive(Crv1); DCrv2 = CagdCrvDerive(Crv2); Srf1 = CagdPromoteCrvToSrf(Crv1, CAGD_CONST_U_DIR); Srf2 = CagdPromoteCrvToSrf(Crv2, CAGD_CONST_V_DIR); CagdSrfDomain(Srf1, &UMin1, &UMax1, &VMin1, &VMax1); CagdSrfDomain(Srf2, &UMin2, &UMax2, &VMin2, &VMax2); BspKnotAffineTrans2(Srf1 -> VKnotVector, Srf1 -> VLength + Srf1 -> VOrder, VMin2, VMax2); BspKnotAffineTrans2(Srf2 -> UKnotVector, Srf2 -> ULength + Srf2 -> UOrder, UMin1, UMax1); DSrf1 = CagdPromoteCrvToSrf(DCrv1, CAGD_CONST_U_DIR); DSrf2 = CagdPromoteCrvToSrf(DCrv2, CAGD_CONST_V_DIR); BspKnotAffineTrans2(DSrf1 -> VKnotVector, DSrf1 -> VLength + DSrf1 -> VOrder, VMin2, VMax2); BspKnotAffineTrans2(DSrf2 -> UKnotVector, DSrf2 -> ULength + DSrf2 -> UOrder, UMin1, UMax1); CagdCrvFree(Crv1); CagdCrvFree(Crv2); CagdCrvFree(DCrv1); CagdCrvFree(DCrv2); if (UseNrmlTan == -1) return SymbCrvBisectorsSrf3D(Srf1, Srf2, DSrf1, DSrf2); SymbSrfSplitScalar(Srf1, &TSrf1, &Srf1X, &Srf1Y, &TSrf1); SymbSrfSplitScalar(Srf2, &TSrf1, &Srf2X, &Srf2Y, &TSrf1); SymbSrfSplitScalar(DSrf1, &TSrf1, &DSrf1X, &DSrf1Y, &TSrf1); SymbSrfSplitScalar(DSrf2, &TSrf1, &DSrf2X, &DSrf2Y, &TSrf1); CagdSrfFree(Srf1); CagdSrfFree(Srf2); CagdSrfFree(DSrf1); CagdSrfFree(DSrf2); TSrf1 = SymbSrfSub(Srf2X, Srf1X); TSrf2 = SymbSrfSub(Srf2Y, Srf1Y); switch (UseNrmlTan) { case 2: /* Use a blend of the two above methods. */ fprintf(stderr, "Blend bisector method is not implemented.\n"); default: case 0: /* Use tangents to compute angular bound. */ TSrf3 = SymbSrfMult(TSrf1, DSrf2X); TSrf4 = SymbSrfMult(TSrf2, DSrf2Y); Alpha = SymbSrfAdd(TSrf3, TSrf4); CagdSrfFree(TSrf3); CagdSrfFree(TSrf4); TSrf3 = SymbSrfMult(DSrf1Y, TSrf2); TSrf4 = SymbSrfMult(DSrf1X, TSrf1); Beta = SymbSrfAdd(TSrf3, TSrf4); CagdSrfFree(TSrf3); CagdSrfFree(TSrf4); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); break; case 1: /* Use normal to compute angular bound. */ TSrf3 = SymbSrfMult(TSrf1, DSrf2Y); TSrf4 = SymbSrfMult(TSrf2, DSrf2X); Alpha = SymbSrfSub(TSrf3, TSrf4); CagdSrfFree(TSrf3); CagdSrfFree(TSrf4); TSrf3 = SymbSrfMult(DSrf1X, TSrf2); TSrf4 = SymbSrfMult(DSrf1Y, TSrf1); Beta = SymbSrfSub(TSrf3, TSrf4); CagdSrfFree(TSrf3); CagdSrfFree(TSrf4); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); break; } /* Adds up the components of "|| Alpha N(t) ||^2 - || Beta N(s) ||^2" */ TSrf1 = SymbSrfMult(Alpha, DSrf1X); Res = SymbSrfMult(TSrf1, TSrf1); CagdSrfFree(TSrf1); TSrf1 = SymbSrfMult(Alpha, DSrf1Y); TSrf2 = SymbSrfMult(TSrf1, TSrf1); TSrf3 = SymbSrfAdd(Res, TSrf2); CagdSrfFree(Res); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); Res = TSrf3; TSrf1 = SymbSrfMult(Beta, DSrf2X); TSrf2 = SymbSrfMult(TSrf1, TSrf1); TSrf3 = SymbSrfSub(Res, TSrf2); CagdSrfFree(Res); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); Res = TSrf3; TSrf1 = SymbSrfMult(Beta, DSrf2Y); TSrf2 = SymbSrfMult(TSrf1, TSrf1); TSrf3 = SymbSrfSub(Res, TSrf2); CagdSrfFree(Res); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); Res = TSrf3; CagdSrfFree(Alpha); CagdSrfFree(Beta); CagdSrfFree(Srf1X); CagdSrfFree(Srf1Y); CagdSrfFree(Srf2X); CagdSrfFree(Srf2Y); CagdSrfFree(DSrf1X); CagdSrfFree(DSrf1Y); CagdSrfFree(DSrf2X); CagdSrfFree(DSrf2Y); return Res; } /***************************************************************************** * DESCRIPTION: M * Computes the bisector surface definition of a given curve or two. M * If Crv contains a list of two curves the bisector between the two curves M * is computed. Otherwise, Crv self--bisectors are sought. The result is a M * scalar surface whose zero set is the set of bisector(s) of the curves. M * Solve for the planar bisector surface in the plane and then elevate in Z M * using the rational function of M * ||P - C1(s)||^2 - ||P - C2(t)||^2. V * * * PARAMETERS: M * Crv: Either one or two curves to compute bisectors for. Assumes M * E2 curves. M * * * RETURN VALUE: M * CagdSrfStruct *: A scalar surface whose zero set provides matching M * bisecting points on Crv, if UseNrmlTan >= 0. The real M * bisector surface for three space curves, if UseNrmlTan = -1. M * * * SEE ALSO: M * SymbCrvCnvxHull, SymbCrvDiameter, SymbCrvBisectors, M * SymbCrvBisectorsSrf, SymbCrvPtBisectorsSrf3D, SymbCrvCrvBisectorSrf3D M * * * KEYWORDS: M * SymbCrvBisectorsSrf, bisectors, skeleton M *****************************************************************************/ CagdSrfStruct *SymbCrvBisectorsSrf2(CagdCrvStruct *Crv) { CagdRType UMin1, UMax1, VMin1, VMax1, UMin2, UMax2, VMin2, VMax2; CagdCrvStruct *Crv1, *Crv2, *DCrv1, *DCrv2; CagdSrfStruct *Srf1, *Srf2, *DSrf1, *DSrf2, *Res, *Denom, *NumerX, *NumerY, *NumerZ, *Srf1X, *Srf1Y, *Srf2X, *Srf2Y, *DSrf1X, *DSrf1Y, *DSrf2X, *DSrf2Y, *TSrf1, *TSrf2, *TSrf3, *TSrf4, *TSrf5, *TSrf6; Crv1 = CagdCoerceCrvTo(Crv, CAGD_PT_E2_TYPE); Crv2 = CagdCoerceCrvTo(Crv -> Pnext ? Crv -> Pnext : Crv, CAGD_PT_E2_TYPE); if (CAGD_IS_BEZIER_CRV(Crv1)) { CagdCrvStruct *TCrv = CnvrtBezier2BsplineCrv(Crv1); CagdCrvFree(Crv1); Crv1 = TCrv; } if (CAGD_IS_BEZIER_CRV(Crv2)) { CagdCrvStruct *TCrv = CnvrtBezier2BsplineCrv(Crv2); CagdCrvFree(Crv2); Crv2 = TCrv; } DCrv1 = CagdCrvDerive(Crv1); DCrv2 = CagdCrvDerive(Crv2); Srf1 = CagdPromoteCrvToSrf(Crv1, CAGD_CONST_U_DIR); Srf2 = CagdPromoteCrvToSrf(Crv2, CAGD_CONST_V_DIR); CagdSrfDomain(Srf1, &UMin1, &UMax1, &VMin1, &VMax1); CagdSrfDomain(Srf2, &UMin2, &UMax2, &VMin2, &VMax2); BspKnotAffineTrans2(Srf1 -> VKnotVector, Srf1 -> VLength + Srf1 -> VOrder, VMin2, VMax2); BspKnotAffineTrans2(Srf2 -> UKnotVector, Srf2 -> ULength + Srf2 -> UOrder, UMin1, UMax1); DSrf1 = CagdPromoteCrvToSrf(DCrv1, CAGD_CONST_U_DIR); DSrf2 = CagdPromoteCrvToSrf(DCrv2, CAGD_CONST_V_DIR); BspKnotAffineTrans2(DSrf1 -> VKnotVector, DSrf1 -> VLength + DSrf1 -> VOrder, VMin2, VMax2); BspKnotAffineTrans2(DSrf2 -> UKnotVector, DSrf2 -> ULength + DSrf2 -> UOrder, UMin1, UMax1); CagdCrvFree(Crv1); CagdCrvFree(Crv2); CagdCrvFree(DCrv1); CagdCrvFree(DCrv2); SymbSrfSplitScalar(Srf1, &TSrf1, &Srf1X, &Srf1Y, &TSrf1); SymbSrfSplitScalar(Srf2, &TSrf1, &Srf2X, &Srf2Y, &TSrf1); SymbSrfSplitScalar(DSrf1, &TSrf1, &DSrf1X, &DSrf1Y, &TSrf1); SymbSrfSplitScalar(DSrf2, &TSrf1, &DSrf2X, &DSrf2Y, &TSrf1); Denom = SymbSrfDeterminant2(DSrf1X, DSrf1Y, DSrf2X, DSrf2Y); TSrf1 = SymbSrfDotProd(Srf1, DSrf1); TSrf2 = SymbSrfDotProd(Srf2, DSrf2); NumerX = SymbSrfDeterminant2(TSrf1, DSrf1Y, TSrf2, DSrf2Y); NumerY = SymbSrfDeterminant2(DSrf1X, TSrf1, DSrf2X, TSrf2); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); CagdSrfFree(Srf1); CagdSrfFree(Srf2); CagdSrfFree(DSrf1); CagdSrfFree(DSrf2); CagdSrfFree(DSrf1X); CagdSrfFree(DSrf1Y); CagdSrfFree(DSrf2X); CagdSrfFree(DSrf2Y); /* Compute || P - C1(s) ||^2 - || P - C2(t) ||^2 */ TSrf1 = SymbSrfMult(Srf1X, Denom); TSrf2 = SymbSrfSub(TSrf1, NumerX); TSrf3 = SymbSrfMult(TSrf2, TSrf2); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); TSrf1 = SymbSrfMult(Srf1Y, Denom); TSrf2 = SymbSrfSub(TSrf1, NumerY); TSrf4 = SymbSrfMult(TSrf2, TSrf2); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); TSrf5 = SymbSrfAdd(TSrf3, TSrf4); CagdSrfFree(TSrf3); CagdSrfFree(TSrf4); TSrf1 = SymbSrfMult(Srf2X, Denom); TSrf2 = SymbSrfSub(TSrf1, NumerX); TSrf3 = SymbSrfMult(TSrf2, TSrf2); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); TSrf1 = SymbSrfMult(Srf2Y, Denom); TSrf2 = SymbSrfSub(TSrf1, NumerY); TSrf4 = SymbSrfMult(TSrf2, TSrf2); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); TSrf6 = SymbSrfAdd(TSrf3, TSrf4); CagdSrfFree(TSrf3); CagdSrfFree(TSrf4); NumerZ = SymbSrfSub(TSrf5, TSrf6); CagdSrfFree(TSrf5); CagdSrfFree(TSrf6); CagdSrfFree(Srf1X); CagdSrfFree(Srf1Y); CagdSrfFree(Srf2X); CagdSrfFree(Srf2Y); CagdMakeSrfsCompatible(&NumerZ, &NumerX, TRUE, TRUE, TRUE, TRUE); CagdMakeSrfsCompatible(&NumerZ, &NumerY, TRUE, TRUE, TRUE, TRUE); CagdMakeSrfsCompatible(&NumerZ, &Denom, TRUE, TRUE, TRUE, TRUE); CagdMakeSrfsCompatible(&NumerY, &NumerX, TRUE, TRUE, TRUE, TRUE); CagdMakeSrfsCompatible(&NumerY, &Denom, TRUE, TRUE, TRUE, TRUE); CagdMakeSrfsCompatible(&NumerX, &Denom, TRUE, TRUE, TRUE, TRUE); Res = SymbSrfMergeScalar(Denom, NumerX, NumerY, NumerZ); CagdSrfFree(NumerX); CagdSrfFree(NumerY); CagdSrfFree(NumerZ); CagdSrfFree(Denom); return Res; } /***************************************************************************** * DESCRIPTION: M * Computes the bisector surface of two curve in arbitrary general three M * space position. M * * * PARAMETERS: M * Crv1, Crv2: Two three space curves to compute their bisector surface. M * * * RETURN VALUE: M * CagdSrfStruct *: The bisector surface. M * * * SEE ALSO: M * SymbCrvDiameter, SymbCrvCnvxHull, SymbCrvBisectorsSrf, M * SymbCrvPtBisectorSrf3D M * * * KEYWORDS: M * SymbCrvCrvBisectorSrf3D, bisectors, skeleton M *****************************************************************************/ CagdSrfStruct *SymbCrvCrvBisectorSrf3D(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2) { CagdSrfStruct *BisectSrf; Crv1 = CagdCrvCopy(Crv1); Crv1 -> Pnext = Crv2; BisectSrf = SymbCrvBisectorsSrf(Crv1, -1); CagdCrvFree(Crv1); return BisectSrf; } /***************************************************************************** * DESCRIPTION: M * Computes the expression of a 3 by 3 determinants. M * * * PARAMETERS: M * Srf11, Srf12, ..., Srf33: The nine factors of the determinant. M * * * RETURN VALUE: M * CagdSrfStruct *: A scalar field representing the determinant computation.M * * * SEE ALSO: M * SymbSrfFff, SymbSrfSff, SymbSrfGaussCurvature, SymbSrfMeanEvolute, M * SymbSrfMeanCurvatureSqr, SymbSrfIsoFocalSrf, SymbSrfCurvatureUpperBound, M * SymbSrfIsoDirNormalCurvatureBound, SymbSrfDeterminant2, M * SymbCrvDeterminant3 M * * * KEYWORDS: M * SymbSrfDeterminant3, determinant M *****************************************************************************/ CagdSrfStruct *SymbSrfDeterminant3(CagdSrfStruct *Srf11, CagdSrfStruct *Srf12, CagdSrfStruct *Srf13, CagdSrfStruct *Srf21, CagdSrfStruct *Srf22, CagdSrfStruct *Srf23, CagdSrfStruct *Srf31, CagdSrfStruct *Srf32, CagdSrfStruct *Srf33) { CagdSrfStruct *Prod1 = SymbSrfDeterminant2(Srf22, Srf23, Srf32, Srf33), *Prod1a = SymbSrfMult(Srf11, Prod1), *Prod2 = SymbSrfDeterminant2(Srf21, Srf23, Srf31, Srf33), *Prod2a = SymbSrfMult(Srf12, Prod2), *Prod3 = SymbSrfDeterminant2(Srf21, Srf22, Srf31, Srf32), *Prod3a = SymbSrfMult(Srf13, Prod3), *Sub12 = SymbSrfSub(Prod1a, Prod2a), *Add123 = SymbSrfAdd(Sub12, Prod3a); CagdSrfFree(Prod1); CagdSrfFree(Prod1a); CagdSrfFree(Prod2); CagdSrfFree(Prod2a); CagdSrfFree(Prod3); CagdSrfFree(Prod3a); CagdSrfFree(Sub12); return Add123; } /***************************************************************************** * DESCRIPTION: * * Computes the bisector surface of two curves in arbitrary general three * * space position. Given are the curves and their derivatives as curves * * already promoted to surfaces. * * * * PARAMETERS: * * Srf1: The first curve promoted to a surface. * * Srf2: The second curve promoted to a surface. * * DSrf1: The derivative of the first curve, promoted to a surface. * * DSrf2: The derivative of the second curve, promoted to a surface. * * * * RETURN VALUE: * * CagdSrfStruct *: The bisector surface between these two curves. * *****************************************************************************/ static CagdSrfStruct *SymbCrvBisectorsSrf3D(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2, CagdSrfStruct *DSrf1, CagdSrfStruct *DSrf2) { static CagdPType Trans = { 0.0, 0.0, 0.0 }; CagdSrfStruct *Ret, *Srf1X, *Srf1Y, *Srf1Z, *Srf2X, *Srf2Y, *Srf2Z, *DSrf1X, *DSrf1Y, *DSrf1Z, *DSrf2X, *DSrf2Y, *DSrf2Z, *D1Srf, *D2Srf, *E12, *E12X, *E12Y, *E12Z, *M12Srf, *TSrf1, *TSrf2, *TSrf3, *TSrf4; SymbSrfSplitScalar(Srf1, &TSrf1, &Srf1X, &Srf1Y, &Srf1Z); SymbSrfSplitScalar(Srf2, &TSrf1, &Srf2X, &Srf2Y, &Srf2Z); SymbSrfSplitScalar(DSrf1, &TSrf1, &DSrf1X, &DSrf1Y, &DSrf1Z); SymbSrfSplitScalar(DSrf2, &TSrf1, &DSrf2X, &DSrf2Y, &DSrf2Z); D1Srf = SymbSrfDotProd(Srf1, DSrf1); D2Srf = SymbSrfDotProd(Srf2, DSrf2); TSrf1 = SymbSrfDotProd(Srf1, Srf1); TSrf2 = SymbSrfDotProd(Srf2, Srf2); M12Srf = SymbSrfSub(TSrf2, TSrf1); CagdSrfTransform(M12Srf, Trans, 0.5); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); E12 = SymbSrfSub(Srf2, Srf1); SymbSrfSplitScalar(E12, &TSrf1, &E12X, &E12Y, &E12Z); CagdSrfFree(E12); TSrf1 = SymbSrfDeterminant3(DSrf1X, DSrf1Y, DSrf1Z, DSrf2X, DSrf2Y, DSrf2Z, E12X, E12Y, E12Z); TSrf2 = SymbSrfDeterminant3(D1Srf, DSrf1Y, DSrf1Z, D2Srf, DSrf2Y, DSrf2Z, M12Srf, E12Y, E12Z); TSrf3 = SymbSrfDeterminant3(DSrf1X, D1Srf, DSrf1Z, DSrf2X, D2Srf, DSrf2Z, E12X, M12Srf, E12Z); TSrf4 = SymbSrfDeterminant3(DSrf1X, DSrf1Y, D1Srf, DSrf2X, DSrf2Y, D2Srf, E12X, E12Y, M12Srf); CagdMakeSrfsCompatible(&TSrf1, &TSrf2, TRUE, TRUE, TRUE, TRUE); CagdMakeSrfsCompatible(&TSrf1, &TSrf3, TRUE, TRUE, TRUE, TRUE); CagdMakeSrfsCompatible(&TSrf1, &TSrf4, TRUE, TRUE, TRUE, TRUE); Ret = SymbSrfMergeScalar(TSrf1, TSrf2, TSrf3, TSrf4); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); CagdSrfFree(TSrf3); CagdSrfFree(TSrf4); CagdSrfFree(D1Srf); CagdSrfFree(D2Srf); CagdSrfFree(M12Srf); CagdSrfFree(E12X); CagdSrfFree(E12Y); CagdSrfFree(E12Z); CagdSrfFree(Srf1); CagdSrfFree(Srf2); CagdSrfFree(DSrf1); CagdSrfFree(DSrf2); CagdSrfFree(Srf1X); CagdSrfFree(Srf1Y); CagdSrfFree(Srf1Z); CagdSrfFree(Srf2X); CagdSrfFree(Srf2Y); CagdSrfFree(Srf2Z); CagdSrfFree(DSrf1X); CagdSrfFree(DSrf1Y); CagdSrfFree(DSrf1Z); CagdSrfFree(DSrf2X); CagdSrfFree(DSrf2Y); CagdSrfFree(DSrf2Z); return Ret; } /***************************************************************************** * DESCRIPTION: M * Computes the bisector surface of a curve in arbitrary general three M * space position and a point in three space. M * * * PARAMETERS: M * Crv: Three space curve to compute its bisector surface with Pt. M * Pt: A point in three space to compute its bisector with Crv. M * RulingScale: The scaling factor for the ruling direction. M * * * RETURN VALUE: M * CagdSrfStruct *: The bisector surface. M * * * SEE ALSO: M * SymbCrvDiameter, SymbCrvCnvxHull, SymbCrvBisectorsSrf, M * SymbCrvCrvBisectorSrf3D, SymbSrfPtBisectorSrf3D M * * * KEYWORDS: M * SymbCrvPtBisectorSrf3D, bisector M *****************************************************************************/ CagdSrfStruct *SymbCrvPtBisectorSrf3D(CagdCrvStruct *Crv, CagdPType Pt, CagdRType RulingScale) { int i; CagdPType Trans; CagdPtStruct *Pt1, *Pt2; CagdCrvStruct *DCrv, *ECrv, *FCrv, *GCrv, *TCrv1, *TCrv2, *TCrv3, *TCrv4, *TCrv5, *TCrv6, *DCrvX, *DCrvY, *DCrvZ, *ECrvX, *ECrvY, *ECrvZ, *GCrvX, *GCrvY, *GCrvZ; CagdSrfStruct *Ret, *TSrf1, *TSrf2, *TSrf3, *TSrf4; if (CAGD_IS_RATIONAL_CRV(Crv)) { SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_NO_SUPPORT); return NULL; } Crv = CagdCoerceCrvTo(Crv, CAGD_PT_E3_TYPE); DCrv = CagdCrvDerive(Crv); GCrv = CagdCrvCopy(Crv); PT_SCALE2(Trans, Pt, -1); CagdCrvTransform(GCrv, Trans, 1.0); ECrv = SymbCrvCrossProd(GCrv, DCrv); FCrv = SymbCrvDotProd(Crv, Crv); PT_CLEAR(Trans); CagdCrvTransform(FCrv, Trans, 0.5); Trans[0] = -DOT_PROD(Pt, Pt) / 2.0; CagdCrvTransform(FCrv, Trans, 1.0); TCrv1 = SymbCrvDotProd(Crv, DCrv); TCrv2 = SymbCrvDotProd(Crv, ECrv); SymbCrvSplitScalar(DCrv, &TCrv3, &DCrvX, &DCrvY, &DCrvZ); SymbCrvSplitScalar(ECrv, &TCrv3, &ECrvX, &ECrvY, &ECrvZ); SymbCrvSplitScalar(GCrv, &TCrv3, &GCrvX, &GCrvY, &GCrvZ); TCrv3 = SymbCrvDeterminant3(GCrvX, GCrvY, GCrvZ, DCrvX, DCrvY, DCrvZ, ECrvX, ECrvY, ECrvZ); TCrv4 = SymbCrvDeterminant3(FCrv, GCrvY, GCrvZ, TCrv1, DCrvY, DCrvZ, TCrv2, ECrvY, ECrvZ); TCrv5 = SymbCrvDeterminant3(GCrvX, FCrv, GCrvZ, DCrvX, TCrv1, DCrvZ, ECrvX, TCrv2, ECrvZ); TCrv6 = SymbCrvDeterminant3(GCrvX, GCrvY, FCrv, DCrvX, DCrvY, TCrv1, ECrvX, ECrvY, TCrv2); CagdMakeCrvsCompatible(&TCrv3, &TCrv4, TRUE, TRUE); CagdMakeCrvsCompatible(&TCrv3, &TCrv5, TRUE, TRUE); CagdMakeCrvsCompatible(&TCrv3, &TCrv6, TRUE, TRUE); CagdCrvFree(TCrv1); CagdCrvFree(TCrv2); TCrv1 = SymbCrvMergeScalar(TCrv3, TCrv4, TCrv5, TCrv6); CagdCrvFree(TCrv3); CagdCrvFree(TCrv4); CagdCrvFree(TCrv5); CagdCrvFree(TCrv6); CagdCrvFree(Crv); CagdCrvFree(DCrv); CagdCrvFree(FCrv); CagdCrvFree(GCrv); CagdCrvFree(DCrvX); CagdCrvFree(DCrvY); CagdCrvFree(DCrvZ); CagdCrvFree(ECrvX); CagdCrvFree(ECrvY); CagdCrvFree(ECrvZ); CagdCrvFree(GCrvX); CagdCrvFree(GCrvY); CagdCrvFree(GCrvZ); /* Convert the ruled surface directrix (TSrf1) into a surface in the U */ /* direction and add the ruling direction as V direction. */ TSrf1 = CagdPromoteCrvToSrf(TCrv1, CAGD_CONST_U_DIR); CagdCrvFree(TCrv1); TSrf2 = CagdPromoteCrvToSrf(ECrv, CAGD_CONST_U_DIR); CagdCrvFree(ECrv); Pt1 = CagdPtNew(); Pt2 = CagdPtNew(); for (i = 0; i < 3; i++) { Pt1 -> Pt[i] = -RulingScale; Pt2 -> Pt[i] = RulingScale; } TCrv2 = CagdMergePtPt(Pt1, Pt2); CagdPtFree(Pt1); CagdPtFree(Pt2); TSrf3 = CagdPromoteCrvToSrf(TCrv2, CAGD_CONST_V_DIR); if (CAGD_IS_BSPLINE_SRF(TSrf2)) { CagdRType UMin, UMax, VMin, VMax; TSrf4 = CnvrtBezier2BsplineSrf(TSrf3); CagdSrfFree(TSrf3); CagdSrfDomain(TSrf2, &UMin, &UMax, &VMin, &VMax); BspKnotAffineTrans2(TSrf4 -> UKnotVector, TSrf4 -> ULength + TSrf4 -> UOrder, UMin, UMax); TSrf3 = TSrf4; } TSrf4 = SymbSrfMult(TSrf2, TSrf3); CagdSrfFree(TSrf2); CagdSrfFree(TSrf3); Ret = SymbSrfAdd(TSrf1, TSrf4); CagdSrfFree(TSrf1); CagdSrfFree(TSrf4); return Ret; } /***************************************************************************** * DESCRIPTION: M * Computes the bisector surface of a surface in arbitrary general three M * space position and a point in three space. M * Solution bisector surface R is derived by solving the three linear M * equations of (S for Srf, P for Pt): M * < dS/du, R > = < dS/du, S > V * < dS/dv, R > = < dS/dv, S > V * < S - P, R > = (< S, S > - < P, P >) / 2 V * * * PARAMETERS: M * Srf: Three space surface too compute its bisector surface with Pt. M * Pt: A point in three space to compute its bisector with Srf. M * * * RETURN VALUE: M * CagdSrfStruct *: The bisector surface. M * * * SEE ALSO: M * SymbCrvDiameter, SymbCrvCnvxHull, SymbCrvBisectorsSrf, M * SymbCrvCrvBisectorSrf3D M * * * KEYWORDS: M * SymbSrfPtBisectorSrf3D, bisector M *****************************************************************************/ CagdSrfStruct *SymbSrfPtBisectorSrf3D(CagdSrfStruct *Srf, CagdPType Pt) { CagdPType Trans; CagdSrfStruct *DuSrfX, *DuSrfY, *DuSrfZ, *DvSrfX, *DvSrfY, *DvSrfZ, *DuSrf, *DvSrf, *B1Srf, *B2Srf, *B3Srf, *SrfXPt, *SrfYPt, *SrfZPt, *TSrf1, *TSrf2, *TSrf3, *TSrf4, *TSrf5, *TSrf6; if (CAGD_IS_RATIONAL_SRF(Srf)) { SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_NO_SUPPORT); return NULL; } if (CAGD_NUM_OF_PT_COORD(Srf -> PType) < 3) Srf = CagdCoerceSrfTo(Srf, CAGD_PT_E3_TYPE); else Srf = CagdSrfCopy(Srf); DuSrf = CagdSrfDerive(Srf, CAGD_CONST_U_DIR); DvSrf = CagdSrfDerive(Srf, CAGD_CONST_V_DIR); /* Forming the Ax = B system. */ B1Srf = SymbSrfDotProd(Srf, DuSrf); B2Srf = SymbSrfDotProd(Srf, DvSrf); B3Srf = SymbSrfDotProd(Srf, Srf); PT_RESET(Trans); Trans[0] = -DOT_PROD(Pt, Pt); CagdSrfTransform(B3Srf, Trans, 0.5); SymbSrfSplitScalar(DuSrf, &TSrf1, &DuSrfX, &DuSrfY, &DuSrfZ); SymbSrfSplitScalar(DvSrf, &TSrf1, &DvSrfX, &DvSrfY, &DvSrfZ); CagdSrfFree(DuSrf); CagdSrfFree(DvSrf); TSrf1 = CagdSrfCopy(Srf); PT_COPY(Trans, Pt); PT_SCALE(Trans, -1.0) CagdSrfTransform(TSrf1, Trans, 1.0); SymbSrfSplitScalar(TSrf1, &TSrf2, &SrfXPt, &SrfYPt, &SrfZPt); CagdSrfFree(TSrf1); /* And solving the lienar system. */ TSrf3 = SymbSrfDeterminant3(DuSrfX, DuSrfY, DuSrfZ, DvSrfX, DvSrfY, DvSrfZ, SrfXPt, SrfYPt, SrfZPt); TSrf4 = SymbSrfDeterminant3(B1Srf, DuSrfY, DuSrfZ, B2Srf, DvSrfY, DvSrfZ, B3Srf, SrfYPt, SrfZPt); TSrf5 = SymbSrfDeterminant3(DuSrfX, B1Srf, DuSrfZ, DvSrfX, B2Srf, DvSrfZ, SrfXPt, B3Srf, SrfZPt); TSrf6 = SymbSrfDeterminant3(DuSrfX, DuSrfY, B1Srf, DvSrfX, DvSrfY, B2Srf, SrfXPt, SrfYPt, B3Srf); CagdSrfFree(DuSrfX); CagdSrfFree(DuSrfY); CagdSrfFree(DuSrfZ); CagdSrfFree(DvSrfX); CagdSrfFree(DvSrfY); CagdSrfFree(DvSrfZ); CagdSrfFree(SrfXPt); CagdSrfFree(SrfYPt); CagdSrfFree(SrfZPt); CagdSrfFree(B1Srf); CagdSrfFree(B2Srf); CagdSrfFree(B3Srf); CagdMakeSrfsCompatible(&TSrf3, &TSrf4, TRUE, TRUE, TRUE, TRUE); CagdMakeSrfsCompatible(&TSrf3, &TSrf5, TRUE, TRUE, TRUE, TRUE); CagdMakeSrfsCompatible(&TSrf3, &TSrf6, TRUE, TRUE, TRUE, TRUE); TSrf1 = SymbSrfMergeScalar(TSrf3, TSrf4, TSrf5, TSrf6); CagdSrfFree(TSrf3); CagdSrfFree(TSrf4); CagdSrfFree(TSrf5); CagdSrfFree(TSrf6); CagdSrfFree(Srf); return TSrf1; } /***************************************************************************** * DESCRIPTION: M * Computes the expression of a 3 by 3 determinants. M * * * PARAMETERS: M * Crv11, Crv12, ..., Crv33: The nine factors of the determinant. M * * * RETURN VALUE: M * CagdCrvStruct *: A scalar field representing the determinant computation.M * * * SEE ALSO: M * SymbCrvDeterminant2, SymbSrfDeterminant3 M * * * KEYWORDS: M * SymbCrvDeterminant3, determinant M *****************************************************************************/ CagdCrvStruct *SymbCrvDeterminant3(CagdCrvStruct *Crv11, CagdCrvStruct *Crv12, CagdCrvStruct *Crv13, CagdCrvStruct *Crv21, CagdCrvStruct *Crv22, CagdCrvStruct *Crv23, CagdCrvStruct *Crv31, CagdCrvStruct *Crv32, CagdCrvStruct *Crv33) { CagdCrvStruct *Prod1 = SymbCrvDeterminant2(Crv22, Crv23, Crv32, Crv33), *Prod1a = SymbCrvMult(Crv11, Prod1), *Prod2 = SymbCrvDeterminant2(Crv21, Crv23, Crv31, Crv33), *Prod2a = SymbCrvMult(Crv12, Prod2), *Prod3 = SymbCrvDeterminant2(Crv21, Crv22, Crv31, Crv32), *Prod3a = SymbCrvMult(Crv13, Prod3), *Sub12 = SymbCrvSub(Prod1a, Prod2a), *Add123 = SymbCrvAdd(Sub12, Prod3a); CagdCrvFree(Prod1); CagdCrvFree(Prod1a); CagdCrvFree(Prod2); CagdCrvFree(Prod2a); CagdCrvFree(Prod3); CagdCrvFree(Prod3a); CagdCrvFree(Sub12); return Add123; } /***************************************************************************** * DESCRIPTION: M * Computes the expression of Crv11 * Crv22 - Crv12 * Crv21, which is a M * determinant of a 2 by 2 matrix. M * * * PARAMETERS: M * Crv11, Crv12, Crv21, Crv22: The four factors of the determinant. M * * * RETURN VALUE: M * CagdCrvStruct *: A scalar field representing the determinant computation.M * * * SEE ALSO: M * SymbCrvDeterminant3, SymbSrfDeterminant2 M * * * KEYWORDS: M * SymbCrvDeterminant2, determinant M *****************************************************************************/ CagdCrvStruct *SymbCrvDeterminant2(CagdCrvStruct *Crv11, CagdCrvStruct *Crv12, CagdCrvStruct *Crv21, CagdCrvStruct *Crv22) { CagdCrvStruct *Prod1 = SymbCrvMult(Crv11, Crv22), *Prod2 = SymbCrvMult(Crv21, Crv12), *Add12 = SymbCrvSub(Prod1, Prod2); CagdCrvFree(Prod1); CagdCrvFree(Prod2); return Add12; }