/****************************************************************************** * SymbPoly.c - Generic Curve/Surface to polygon/polylines conversion routines.* ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, July. 90. * ******************************************************************************/ #include "symb_loc.h" #include "geom_lib.h" #define CURVE_PTS_DIST_RESOLUTION 10 #define SYMB_REL_TOL_BOUND 1e-7 #define SYMB_CRV_POLY_EPS 1e-4 #define SYMB_MAX_CRV_SUBDIV_LEN 1e3 /* No subdiv. of larger crvs. */ #define SYMB_CRV_OPTI_TOL_PT_COORD(WPoints, Points, Axis, Inx) \ (WPoints == NULL ? Points[Axis][Inx] : Points[Axis][Inx] / WPoints[Inx]); STATIC_DATA CagdCrvStruct *GlblDeriv1MagSqrCrv = NULL, *GlblCrvtrCrv = NULL; STATIC_DATA CagdRType GlblCrvtrCrvTMin, GlblCrvtrCrvTMax; static CagdPolylineStruct *SymbCrv2OptiTlrncPolyline(CagdCrvStruct *Crv, CagdRType Tolerance); static CagdPtStruct *SymbCrv2OptiTlrncPolyAux(CagdCrvStruct *Crv, CagdRType Tolerance, CagdRType TMin, CagdRType TMax, CagdBType AddFirstPt); static CagdPolylineStruct *SymbCrv2OptiCrvtrPolyline(CagdCrvStruct *Crv, int SamplesPerCurve); static CagdRType CrvCurvatureEvalFunc(CagdRType t); /***************************************************************************** * DESCRIPTION: M * Routine to convert a single surface to set of triangles approximating it. M * FineNess is a fineness control on result and the larger it is more M * triangles may result. M * A value of 10 is a good starting value. M * NULL is returned in case of an error, otherwise list of CagdPolygonStruct. M * * * PARAMETERS: M * Srf: To approximate into triangles. M * FineNess: Control on accuracy, the higher the finer. M * ComputeNormals: If TRUE, normal information is also computed. M * FourPerFlat: If TRUE, four triangles are created per flat surface. M * If FALSE, only 2 triangles are created. M * ComputeUV: If TRUE, UV values are stored and returned as well. M * * * RETURN VALUE: M * CagdPolygonStruct *: A list of polygons with optional normal and/or M * UV parametric information. M * NULL is returned in case of an error. M * * * SEE ALSO: M * BzrSrf2Polygons, IritSurface2Polygons, IritTrimSrf2Polygons, M * BspSrf2Polygons, TrimSrf2Polygons M * * * KEYWORDS: M * SymbSrf2Polygons, polygonization, surface approximation M *****************************************************************************/ CagdPolygonStruct *SymbSrf2Polygons(CagdSrfStruct *Srf, int FineNess, CagdBType ComputeNormals, CagdBType FourPerFlat, CagdBType ComputeUV) { /* Make sure we do not deal with constant surfaces. */ if (Srf -> UOrder < 2 || Srf -> VOrder < 2) { SYMB_FATAL_ERROR(SYMB_ERR_POLY_CONST_SRF); return NULL; } switch (Srf -> GType) { case CAGD_SBEZIER_TYPE: return BzrSrf2Polygons(Srf, FineNess, ComputeNormals, FourPerFlat, ComputeUV); case CAGD_SBSPLINE_TYPE: return BspSrf2Polygons(Srf, FineNess, ComputeNormals, FourPerFlat, ComputeUV); case CAGD_SPOWER_TYPE: SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT); return NULL; default: SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_SRF); return NULL; } } /***************************************************************************** * DESCRIPTION: M * Routine to convert a single surface to NumOfIsolines polylines in each M * parametric direction with SamplesPerCurve in each isoparametric curve. M * Polyline are always E3 of CagdPolylineStruct type. M * NULL is returned in case of an error, otherwise list of M * CagdPolylineStruct. Attempt is made to extract isolines along C1 M * discontinuities first. M * * * PARAMETERS: M * Srf: Srf to extract isoparametric curves from. M * NumOfIsocurves: To extarct from Srf in each (U or V) direction. M * TolSamples: Tolerance of approximation error (Method = 2) or M * Number of samples to compute on polyline (Method = 0, 1). M * Method: 0 - TolSamples are set uniformly in parametric space, M * 1 - TolSamples are set optimally, considering the M * isocurve's curvature. M * 2 - TolSamples sets the maximum error allowed between the M * piecewise linear approximation and original curve. M * * * RETURN VALUE: M * CagdPolylineStruct *: List of polylines representing a piecewise linear M * approximation of the extracted isoparamteric M * curves or NULL is case of an error. M * * * SEE ALSO: M * BspSrf2Polylines, BzrSrf2Polylines, IritSurface2Polylines, M * IritTrimSrf2Polylines, TrimSrf2Polylines M * * * KEYWORDS: M * SymbSrf2Polylines, polylines, isoparametric curves M *****************************************************************************/ CagdPolylineStruct *SymbSrf2Polylines(CagdSrfStruct *Srf, int NumOfIsocurves[2], CagdRType TolSamples, SymbCrvApproxMethodType Method) { CagdCrvStruct *Crv, *Crvs; CagdPolylineStruct *Poly, *Polys; switch (Method) { case SYMB_CRV_APPROX_TOLERANCE: case SYMB_CRV_APPROX_CURVATURE: Crvs = SymbSrf2Curves(Srf, NumOfIsocurves); Crvs = SymbSubdivCrvsAtPoles(Crvs, SYMB_CRV_POLY_EPS); Polys = NULL; for (Crv = Crvs; Crv != NULL; Crv = Crv -> Pnext) { if (Method == SYMB_CRV_APPROX_TOLERANCE) { Poly = SymbCrv2OptiTlrncPolyline(Crv, TolSamples); } else { Poly = SymbCrv2OptiCrvtrPolyline(Crv, MAX((int) TolSamples, 2)); } Poly -> Pnext = Polys; Polys = Poly; } CagdCrvFreeList(Crvs); return Polys; default: case SYMB_CRV_APPROX_UNIFORM: switch (Srf -> GType) { case CAGD_SBEZIER_TYPE: return BzrSrf2Polylines(Srf, NumOfIsocurves, (int) TolSamples); case CAGD_SBSPLINE_TYPE: return BspSrf2Polylines(Srf, NumOfIsocurves, (int) TolSamples); case CAGD_SPOWER_TYPE: SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT); return NULL; default: SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_SRF); return NULL; } } } /***************************************************************************** * DESCRIPTION: M * Split the given list of curves at the poles of the (rational) curve, M * if any, in place. M * * * PARAMETERS: M * Crvs: List of curves to split, in place. M * Eps: Tolerance to clip the curve away from the poles as computed. M * * * RETURN VALUE: M * CagdCrvStruct *: List of curves with no poles (but with possible M * all negative weights). M * * * SEE ALSO: M * CagdCrvSubdivAtParams M * * * KEYWORDS: M * SymbSubdivCrvsAtPoles M *****************************************************************************/ CagdCrvStruct *SymbSubdivCrvsAtPoles(CagdCrvStruct *Crvs, CagdRType Eps) { int ClipPoles = Cagd2PolyClipPolysAtPoles(CAGD_QUERY_VALUE); CagdCrvStruct *Crv, *RetCrvs = NULL; if (!ClipPoles) return Crvs; while (Crvs != NULL) { LIST_POP(Crv, Crvs); if (!CAGD_IS_RATIONAL_CRV(Crv)) { LIST_PUSH(Crv, RetCrvs); } else { CagdPtStruct *Pts = SymbCrvZeroSet(Crv, 0, Eps * 0.1); if (Pts == NULL) { LIST_PUSH(Crv, RetCrvs); } else { CagdCrvStruct *SplitCrvs = CagdCrvSubdivAtParams(Crv, Pts, 0.0); CagdPtFreeList(Pts); CagdCrvFree(Crv); while (SplitCrvs != NULL) { CagdRType TMin, TMax; /* Trim the end points of the curve segement to get rid */ /* of the poles. */ CagdCrvDomain(SplitCrvs, &TMin, &TMax); if (TMin + Eps < TMax - Eps) { Crv = CagdCrvRegionFromCrv(SplitCrvs, TMin + Eps, TMax - Eps); LIST_PUSH(Crv, RetCrvs); } LIST_POP(Crv, SplitCrvs); CagdCrvFree(Crv); } } } } return RetCrvs; } /***************************************************************************** * DESCRIPTION: M * Routine to extract from a surface NumOfIsoline isocurve list M * in each param. direction. M * Iso parametric curves are sampled equally spaced in parametric space. M * NULL is returned in case of an error, otherwise list of CagdCrvStruct. M * * * PARAMETERS: M * Srf: To extract isoparametric curves from. M * NumOfIsocurves: In each (U or V) direction. M * * * RETURN VALUE: M * CagdCrvStruct *: List of extracted isoparametric curves. These curves M * inherit the order and continuity of the original Srf. M * NULL is returned in case of an error. M * * * SEE ALSO: M * BspSrf2PCurves, BzrSrf2Curves, CagdSrf2Curves M * * * KEYWORDS: M * SymbSrf2Curves, curves, isoparametric curves M *****************************************************************************/ CagdCrvStruct *SymbSrf2Curves(CagdSrfStruct *Srf, int NumOfIsocurves[2]) { switch (Srf -> GType) { case CAGD_SBEZIER_TYPE: return BzrSrf2Curves(Srf, NumOfIsocurves); case CAGD_SBSPLINE_TYPE: return BspSrf2Curves(Srf, NumOfIsocurves); case CAGD_SPOWER_TYPE: SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT); return NULL; default: SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_SRF); return NULL; } } /***************************************************************************** * DESCRIPTION: M * Routine to approx. a single curve as a polyline with TolSamples M * samples/tolerance. Polyline is always E3 CagdPolylineStruct type. M * NULL is returned in case of an error, otherwise CagdPolylineStruct. M * * * PARAMETERS: M * Crv: To approximate as a polyline. M * TolSamples: Tolerance of approximation error (Method = 2) or M * Number of samples to compute on polyline (Method = 0, 1). M * Method: 0 - TolSamples are set uniformly in parametric space, M * 1 - TolSamples are set optimally, considering the M * isocurve's curvature. M * 2 - TolSamples sets the maximum error allowed between the M * piecewise linear approximation and original curve. M * OptiLin: If TRUE, optimize linear curves. M * * * RETURN VALUE: M * CagdPolylineStruct *: A polyline representing the piecewise linear M * approximation from, or NULL in case of an error. M * * * SEE ALSO: M * BspCrv2Polyline, BzrCrv2Polyline, IritCurve2Polylines M * * * KEYWORDS: M * SymbCrv2Polyline, piecewise linear approximation, polyline M *****************************************************************************/ CagdPolylineStruct *SymbCrv2Polyline(CagdCrvStruct *Crv, CagdRType TolSamples, SymbCrvApproxMethodType Method, CagdBType OptiLin) { int i, n, NewCrvs, ClipPoles = Cagd2PolyClipPolysAtPoles(CAGD_QUERY_VALUE); CagdCrvStruct *Crvs, *TCrv; CagdPolylineStruct *P, *PList = NULL; if (ClipPoles) { TCrv = CagdCrvCopy(Crv); TCrv -> Pnext = NULL; Crvs = SymbSubdivCrvsAtPoles(TCrv, SYMB_CRV_POLY_EPS); NewCrvs = TRUE; n = CagdListLength(Crvs); } else { Crvs = Crv; NewCrvs = FALSE; n = 1; } for (i = 0, Crv = Crvs; i < n; i++, Crv = Crv -> Pnext) { switch (Method) { case SYMB_CRV_APPROX_TOLERANCE: P = SymbCrv2OptiTlrncPolyline(Crv, TolSamples); break; case SYMB_CRV_APPROX_CURVATURE: if (Crv -> Order > 2) { P = SymbCrv2OptiCrvtrPolyline(Crv, (int) TolSamples); /* else do uniform sampling on linear curves. */ break; } default: case SYMB_CRV_APPROX_UNIFORM: switch (Crv -> GType) { case CAGD_CBEZIER_TYPE: P = BzrCrv2Polyline(Crv, (int) TolSamples); break; case CAGD_CBSPLINE_TYPE: P = BspCrv2Polyline(Crv, (int) TolSamples, NULL, OptiLin); break; case CAGD_CPOWER_TYPE: TCrv = CnvrtPower2BezierCrv(Crv); P = BzrCrv2Polyline(TCrv, (int) TolSamples); CagdCrvFree(TCrv); break; default: SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV); return NULL; } } LIST_PUSH(P, PList); } if (NewCrvs) CagdCrvFreeList(Crvs); return CagdListReverse(PList); } /***************************************************************************** * DESCRIPTION: * * Routine to convert a single curve to piecewise linear polyline. * * Polyline is always E3 of CagdPolylineStruct type. * * Curve is adaptively sampled and result is guaranteed to be within * * Tolerance from the real curve. * * NULL is returned in case of an error, otherwise CagdPolylineStruct. * * If Crv has an integer "SaveParamVals" attribute, a polyline is returned * * as a ptr attribute "SaveParamVals" with the parameter values of the point * * of the piecewise linear approximation. * * * * PARAMETERS: * * Crv: To approximate as a polyline. * * Tolerance: Maximal distance between the curve and its piecewise * * linear approximation. * * * * RETURN VALUE: * * CagdPolylineStruct *: A polyline close within Tolerance to Crv. * *****************************************************************************/ static CagdPolylineStruct *SymbCrv2OptiTlrncPolyline(CagdCrvStruct *Crv, CagdRType Tolerance) { CagdBType BBoxPosWeights, SaveParamVals = !IP_ATTR_IS_BAD_INT(AttrGetIntAttrib(Crv -> Attr, "SaveParamVals")); int PtListLen; CagdRType TMin, TMax, m; CagdCrvStruct *OpenEndCrv = CAGD_IS_BSPLINE_CRV(Crv) && !BspCrvHasOpenEC(Crv) ? BspCrvOpenEnd(Crv) : Crv; CagdPtStruct *Pt, *PtList; CagdPolylineStruct *Params = NULL, *Poly = NULL; CagdBBoxStruct BBox; /* Ignore zero weights while deriving the bbox. */ BBoxPosWeights = CagdIgnoreNonPosWeightBBox(TRUE); CagdCrvBBox(OpenEndCrv, &BBox); CagdIgnoreNonPosWeightBBox(BBoxPosWeights); m = MAX(BBox.Max[0] - BBox.Min[0], BBox.Max[1] - BBox.Min[1]); if (Tolerance < SYMB_REL_TOL_BOUND * m) Tolerance = SYMB_REL_TOL_BOUND * m; if (CAGD_IS_BEZIER_CRV(OpenEndCrv)) { CagdCrvStruct *TCrv = CnvrtBezier2BsplineCrv(OpenEndCrv); if (OpenEndCrv != Crv) CagdCrvFree(OpenEndCrv); OpenEndCrv = TCrv; } if (SaveParamVals) { CagdCrvDomain(OpenEndCrv, &TMin, &TMax); } else { TMin = 1.0; TMax = -1.0; } PtList = SymbCrv2OptiTlrncPolyAux(OpenEndCrv, Tolerance, TMin, TMax, TRUE); PtListLen = CagdListLength(PtList); if (PtListLen > 0) { int i; Poly = CagdPolylineNew(PtListLen + 1); if (SaveParamVals) Params = CagdPolylineNew(PtListLen + 1); for (Pt = PtList, i = 0; Pt != NULL; Pt = Pt -> Pnext, i++) { PT_COPY(Poly -> Polyline[i].Pt, Pt -> Pt); if (SaveParamVals) { PT_RESET(Params -> Polyline[i].Pt); Params -> Polyline[i].Pt[0] = AttrGetRealAttrib(Pt -> Attr, "SaveParamVals"); } } CagdPtFreeList(PtList); /* Add the last point */ CagdCoerceToE3(Poly -> Polyline[PtListLen].Pt, OpenEndCrv -> Points, OpenEndCrv -> Length - 1, OpenEndCrv -> PType); if (SaveParamVals) { PT_RESET(Params -> Polyline[PtListLen].Pt); Params -> Polyline[PtListLen].Pt[0] = TMax; AttrSetPtrAttrib(&Poly -> Attr, "SaveParamVals", Params); } } if (OpenEndCrv != Crv) CagdCrvFree(OpenEndCrv); return Poly; } /***************************************************************************** * DESCRIPTION: * * Auxiliary function of SymbCrv2OptiTlrncPolyline. * * Gets a curve and test if a line from first point to last point is close * * enough (less than Tolerance) to all other control points. Otherwise, * * the curve is subdivided and this function recurses on the two sub curves. * * Last point of curve is never concatenated onto the returned list. * * if TMin < TMax parameter values are saved as real attributes with the * * name of "SaveParamVals". * * * * PARAMETERS: * * Crv: To approximate as a polyline. * * Tolerance: Maximal distance between the curve and its piecewise * * linear approximation. * * TMin, TMax: Domain of curve to save parameter values. If, however, * * TMin > TMax, no parameter values are saved. * * AddFirstPt: If TRUE, first point on curve is also added. * * * * RETURN VALUE: * * CagdPtStruct *: A list of points approximating Crv to within Tolerance. * *****************************************************************************/ static CagdPtStruct *SymbCrv2OptiTlrncPolyAux(CagdCrvStruct *Crv, CagdRType Tolerance, CagdRType TMin, CagdRType TMax, CagdBType AddFirstPt) { CagdPType LinePt1, LinePt2; CagdVType LineDir, Vec; CagdRType t1, t2, **Points = Crv -> Points, *WPoints = Points[0]; int i, MaxDim = CAGD_NUM_OF_PT_COORD(Crv -> PType), Len = Crv -> Length; LineType Line; PlaneType Plane1, Plane2; CagdPtStruct *RetPtList = AddFirstPt ? CagdPtNew() : NULL; if (MaxDim != 2 && MaxDim != 3) { if (MaxDim > 3) MaxDim = 3; else { SYMB_FATAL_ERROR(SYMB_ERR_ONLY_2D_OR_3D); return NULL; } } LinePt1[0] = SYMB_CRV_OPTI_TOL_PT_COORD(WPoints, Points, 1, 0); LinePt1[1] = SYMB_CRV_OPTI_TOL_PT_COORD(WPoints, Points, 2, 0); LinePt1[2] = MaxDim == 2 ? 0.0 : SYMB_CRV_OPTI_TOL_PT_COORD(WPoints, Points, 3, 0); LinePt2[0] = SYMB_CRV_OPTI_TOL_PT_COORD(WPoints, Points, 1, Len - 1); LinePt2[1] = SYMB_CRV_OPTI_TOL_PT_COORD(WPoints, Points, 2, Len - 1); LinePt2[2] = MaxDim == 2 ? 0.0 : SYMB_CRV_OPTI_TOL_PT_COORD(WPoints, Points, 3, Len - 1); LineDir[0] = LinePt2[0] - LinePt1[0]; LineDir[1] = LinePt2[1] - LinePt1[1]; LineDir[2] = MaxDim == 2 ? 0.0 : LinePt2[2] - LinePt1[2]; if (VEC_SQR_LENGTH(LineDir) < IRIT_UEPS) { /* The two points are identical - select the Y line as alterntive. */ LinePt2[1] += 1.0; LineDir[1] += 1.0; } if (MaxDim == 2) { /* Construct a normalized line equation as "Ax + By + C = 0". */ GMLineFrom2Points(Line, LinePt1, LinePt2); } else { /* Construct two normalized orthogonal plane equations through the */ /* 3-space line as "Ax + By + Cz + D = 0". */ VEC_RESET(Vec); /* Create Vec to be orthogonal to LineDir. */ if (FABS(LineDir[0]) <= FABS(LineDir[1]) && FABS(LineDir[0]) <= FABS(LineDir[2])) Vec[0] = 1.0; else if (FABS(LineDir[1]) <= FABS(LineDir[0]) && FABS(LineDir[1]) <= FABS(LineDir[2])) Vec[1] = 1.0; else Vec[2] = 1.0; CROSS_PROD(Plane1, Vec, LineDir); VEC_NORMALIZE(Plane1); Plane1[3] = -DOT_PROD(Plane1, LinePt1); CROSS_PROD(Plane2, Plane1, LineDir); VEC_NORMALIZE(Plane2); Plane2[3] = -DOT_PROD(Plane2, LinePt1); } if (AddFirstPt) { PT_COPY(RetPtList -> Pt, LinePt1); if (TMin < TMax) { /* Add parameter value as real attribute! */ AttrSetRealAttrib(&RetPtList -> Attr, "SaveParamVals", TMin); } } if (MaxDim == 2) { /* Check rigorously against all interior ctl pts. */ if (WPoints == NULL) { for (i = 1; i < Len - 1; i++) { if (FABS(Line[0] * Points[1][i] + Line[1] * Points[2][i] + Line[2]) > Tolerance) break; } } else { for (i = 1; i < Len - 1; i++) { if (FABS(Line[0] * Points[1][i] + Line[1] * Points[2][i] + Line[2] * WPoints[i]) > Tolerance * FABS(WPoints[i])) break; } } } else { /* Check rigorously against all interior ctl pts. */ if (WPoints == NULL) { for (i = 1; i < Len - 1; i++) { if (FABS(Plane1[0] * Points[1][i] + Plane1[1] * Points[2][i] + Plane1[2] * Points[3][i] + Plane1[3]) > Tolerance || FABS(Plane2[0] * Points[1][i] + Plane2[1] * Points[2][i] + Plane2[2] * Points[3][i] + Plane2[3]) > Tolerance) break; } } else { for (i = 1; i < Len - 1; i++) { if (FABS(Plane1[0] * Points[1][i] + Plane1[1] * Points[2][i] + Plane1[2] * Points[3][i] + Plane1[3] * WPoints[i]) > Tolerance * FABS(WPoints[i]) || FABS(Plane2[0] * Points[1][i] + Plane2[1] * Points[2][i] + Plane2[2] * Points[3][i] + Plane2[3] * WPoints[i]) > Tolerance * FABS(WPoints[i])) break; } } } CagdCrvDomain(Crv, &t1, &t2); if (t2 - t1 > IRIT_EPS && i < Len - 1) { /* Tolerance has not been met. */ CagdCrvStruct *Crv1, *Crv2; CagdRType t; CagdPtStruct *RetPtList1, *RetPtList2; if (Len >= SYMB_MAX_CRV_SUBDIV_LEN) { int i, l, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdRType **Points = Crv -> Points, *WPts = Points[0], *NodeVals = CagdCrvNodes(Crv); /* Return the control polygon instead. */ for (i = AddFirstPt ? 1 : 0; i < Len; i++) { CagdPtStruct *PtNew = CagdPtNew(); LIST_PUSH(PtNew, RetPtList); for (l = 0; l < MIN(MaxCoord, 3); l++) PtNew -> Pt[l] = WPts != NULL ? Points[l + 1][i] / WPts[i] : Points[l + 1][i]; for (l = MIN(MaxCoord, 3); l < 3; l++) PtNew -> Pt[l] = 0.0; if (TMin < TMax) { /* Add parameter value as real attribute! */ AttrSetRealAttrib(&RetPtList -> Attr, "SaveParamVals", NodeVals[i]); } } IritFree(NodeVals); return CagdListReverse(RetPtList); } if (CAGD_IS_BSPLINE_CRV(Crv) && Crv -> Length > Crv -> Order) { t = Crv -> KnotVector[(Crv -> Length + Crv -> Order) >> 1]; if (TMin > TMax) { /* No parameter values are needed */ TMin = t + 1.0; TMax = t - 1.0; } } else if (TMin > TMax) { /* No parameter values are needed */ CagdCrvDomain(Crv, &t1, &t2); t = (t1 + t2) * 0.5; TMin = t + 1.0; TMax = t - 1.0; } else { CagdCrvDomain(Crv, &TMin, &TMax); t = (TMin + TMax) * 0.5; } Crv1 = CagdCrvSubdivAtParam(Crv, t); Crv2 = Crv1 -> Pnext; Crv1 -> Pnext = NULL; RetPtList1 = SymbCrv2OptiTlrncPolyAux(Crv1, Tolerance, TMin, t, FALSE); RetPtList2 = SymbCrv2OptiTlrncPolyAux(Crv2, Tolerance, t, TMax, TRUE); if (RetPtList == NULL) RetPtList = CagdListAppend(RetPtList1, RetPtList2); else RetPtList -> Pnext = CagdListAppend(RetPtList1, RetPtList2); CagdCrvFree(Crv1); CagdCrvFree(Crv2); } return RetPtList; } /***************************************************************************** * DESCRIPTION: * * Routine to convert a single curve to polyline with SamplesPerCurve * * samples, using the curvature field of the curve. * * Polyline is always E3 of CagdPolylineStruct type. * * Curve is sampled at optimal locations as to minimize the L-infinity * * error between the curve and its polyline approximation. * * NULL is returned in case of an error, otherwise CagdPolylineStruct. * * * * PARAMETERS: * * Crv: To approiximate as a polyline. * * SamplesPerCurve: Number of samples one can take off the curve. * * * * RETURN VALUE: * * CagdPolylineStruct *: A polyline with SamplesPerCurve samples that * * approixmates Crv. * *****************************************************************************/ static CagdPolylineStruct *SymbCrv2OptiCrvtrPolyline(CagdCrvStruct *Crv, int SamplesPerCurve) { int i, OldMultInterp = BspMultInterpFlag(FALSE); CagdRType *TVals; CagdPolylineStruct *P = CagdPolylineNew(SamplesPerCurve); CagdCrvStruct *CTmp; CagdPolylnStruct *NewPolyline = P -> Polyline; GlblCrvtrCrv = SymbCrv3DCurvatureSqr(Crv); CTmp = CagdCrvDerive(Crv); GlblDeriv1MagSqrCrv = SymbCrvDotProd(CTmp, CTmp); CagdCrvFree(CTmp); CagdCrvDomain(Crv, &GlblCrvtrCrvTMin, &GlblCrvtrCrvTMax); TVals = GMDistPoint1DWithEnergy(SamplesPerCurve, GlblCrvtrCrvTMin, GlblCrvtrCrvTMax, CURVE_PTS_DIST_RESOLUTION, CrvCurvatureEvalFunc); for (i = 0; i < SamplesPerCurve; i++) { /* Convert to E3 polyline. */ CagdRType *R = CagdCrvEval(Crv, TVals[i]); CagdCoerceToE3(NewPolyline[i].Pt, &R, -1, Crv -> PType); } CagdCrvFree(GlblCrvtrCrv); CagdCrvFree(GlblDeriv1MagSqrCrv); IritFree(TVals); BspMultInterpFlag(OldMultInterp); return P; } /***************************************************************************** * DESCRIPTION: * * Evaluation function of curvature square. This function should return the * * square root of the curvature, scaled by the arclength using GlblCrvtrCrv * * GlblDeriv1MagSqrCrv that contain curvature square and arclength sqaure. * * * * PARAMETERS: * * t: Parameter at which to evalate the global curvature field. * * * * RETURN VALUE: * * CagdRType: Rstimated curvature square. * *****************************************************************************/ static CagdRType CrvCurvatureEvalFunc(CagdRType t) { CagdRType CrvtrSqr, MagSqr, *R; if (t < GlblCrvtrCrvTMin) t = GlblCrvtrCrvTMin; if (t > GlblCrvtrCrvTMax) t = GlblCrvtrCrvTMax; R = CagdCrvEval(GlblCrvtrCrv, t); CrvtrSqr = FABS( R[1] / R[0] ); R = CagdCrvEval(GlblDeriv1MagSqrCrv, t); MagSqr = GlblDeriv1MagSqrCrv -> PType == CAGD_PT_E1_TYPE ? R[1] : R[1] / R[0]; return sqrt( sqrt( CrvtrSqr ) * MagSqr ); }