diff --git a/src/GCPnts/FILES b/src/GCPnts/FILES index 29ead1dda0..f990131585 100755 --- a/src/GCPnts/FILES +++ b/src/GCPnts/FILES @@ -10,7 +10,6 @@ GCPnts_QuasiUniformDeflection.cxx GCPnts_QuasiUniformDeflection.pxx GCPnts_QuasiUniformDeflection.hxx GCPnts_TangentialDeflection.cxx -GCPnts_TangentialDeflection.pxx GCPnts_TangentialDeflection.hxx GCPnts_UniformAbscissa.cxx GCPnts_UniformAbscissa.pxx diff --git a/src/GCPnts/GCPnts_TangentialDeflection.cxx b/src/GCPnts/GCPnts_TangentialDeflection.cxx index ac90aab9ee..7b74e20f2f 100644 --- a/src/GCPnts/GCPnts_TangentialDeflection.cxx +++ b/src/GCPnts/GCPnts_TangentialDeflection.cxx @@ -14,111 +14,451 @@ // Alternatively, this file may be used under the terms of Open CASCADE // commercial license or contractual agreement. +#include #include #include -#include #include #include #include #include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include #include #include -inline static void D0 (const Adaptor3d_Curve& C, const Standard_Real U, gp_Pnt& P) +namespace { - C.D0 (U, P); -} + static const Standard_Real Us3 = 0.3333333333333333333333333333; -inline static void D2 (const Adaptor3d_Curve& C, const Standard_Real U, - gp_Pnt& P, gp_Vec& V1, gp_Vec& V2) -{ - C.D2 (U, P, V1, V2); -} + //! Auxiliary tool to resolve 2D/3D curve classes. + template struct CurveTypes {}; - -static void D0 (const Adaptor2d_Curve2d& C, const Standard_Real U, gp_Pnt& PP) -{ - Standard_Real X, Y; - gp_Pnt2d P; - C.D0 (U, P); - P.Coord (X, Y); - PP.SetCoord (X, Y, 0.0); -} - -static void D2 (const Adaptor2d_Curve2d& C, const Standard_Real U, - gp_Pnt& PP, gp_Vec& VV1, gp_Vec& VV2) -{ - Standard_Real X, Y; - gp_Pnt2d P; - gp_Vec2d V1,V2; - C.D2 (U, P, V1, V2); - P.Coord (X, Y); - PP.SetCoord (X, Y, 0.0); - V1.Coord (X, Y); - VV1.SetCoord (X, Y, 0.0); - V2.Coord (X, Y); - VV2.SetCoord (X, Y, 0.0); -} - -static Standard_Real EstimAngl(const gp_Pnt& P1, const gp_Pnt& Pm, const gp_Pnt& P2) -{ - gp_Vec V1(P1, Pm), V2(Pm, P2); - Standard_Real L = V1.Magnitude() * V2.Magnitude(); - // - if(L > gp::Resolution()) + //! Auxiliary tool to resolve 3D curve classes. + template<> struct CurveTypes { - return V1.CrossMagnitude(V2)/L; + typedef Geom_BezierCurve BezierCurve; + typedef Geom_BSplineCurve BSplineCurve; + typedef GCPnts_DistFunction DistFunction; + typedef GCPnts_DistFunctionMV DistFunctionMV; + }; + + //! Auxiliary tool to resolve 2D curve classes. + template<> struct CurveTypes + { + typedef Geom2d_BezierCurve BezierCurve; + typedef Geom2d_BSplineCurve BSplineCurve; + typedef GCPnts_DistFunction2d DistFunction; + typedef GCPnts_DistFunction2dMV DistFunctionMV; + }; + + inline static void D0 (const Adaptor3d_Curve& C, const Standard_Real U, gp_Pnt& P) + { + C.D0 (U, P); + } + + inline static void D2 (const Adaptor3d_Curve& C, const Standard_Real U, + gp_Pnt& P, gp_Vec& V1, gp_Vec& V2) + { + C.D2 (U, P, V1, V2); + } + + static void D0 (const Adaptor2d_Curve2d& C, const Standard_Real U, gp_Pnt& PP) + { + Standard_Real X, Y; + gp_Pnt2d P; + C.D0 (U, P); + P.Coord (X, Y); + PP.SetCoord (X, Y, 0.0); + } + + static void D2 (const Adaptor2d_Curve2d& C, const Standard_Real U, + gp_Pnt& PP, gp_Vec& VV1, gp_Vec& VV2) + { + Standard_Real X, Y; + gp_Pnt2d P; + gp_Vec2d V1,V2; + C.D2 (U, P, V1, V2); + P.Coord (X, Y); + PP.SetCoord (X, Y, 0.0); + V1.Coord (X, Y); + VV1.SetCoord (X, Y, 0.0); + V2.Coord (X, Y); + VV2.SetCoord (X, Y, 0.0); + } + + static Standard_Real EstimAngl (const gp_Pnt& P1, const gp_Pnt& Pm, const gp_Pnt& P2) + { + gp_Vec V1(P1, Pm), V2(Pm, P2); + Standard_Real L = V1.Magnitude() * V2.Magnitude(); + if (L > gp::Resolution()) + { + return V1.CrossMagnitude(V2)/L; + } + else + { + return 0.; + } + } + + // Return number of interval of continuity on which theParam is located. + // Last parameter is used to increase search speed. + static Standard_Integer getIntervalIdx(const Standard_Real theParam, + TColStd_Array1OfReal& theIntervs, + const Standard_Integer thePreviousIdx) + { + Standard_Integer anIdx; + for(anIdx = thePreviousIdx; anIdx < theIntervs.Upper(); anIdx++) + { + if (theParam >= theIntervs(anIdx) && + theParam <= theIntervs(anIdx + 1)) // Inside of anIdx interval. + { + break; + } + } + return anIdx; + } +} + +//======================================================================= +//function : GCPnts_TangentialDeflection +//purpose : +//======================================================================= +GCPnts_TangentialDeflection::GCPnts_TangentialDeflection() +: myAngularDeflection (0.0), + myCurvatureDeflection (0.0), + myUTol (0.0), + myMinNbPnts (0), + myMinLen(0.0), + myLastU (0.0), + myFirstu (0.0) +{ +} + +//======================================================================= +//function : GCPnts_TangentialDeflection +//purpose : +//======================================================================= +GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (const Adaptor3d_Curve& theC, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints, + const Standard_Real theUTol, + const Standard_Real theMinLen) +: myAngularDeflection (0.0), + myCurvatureDeflection (0.0), + myUTol (0.0), + myMinNbPnts (0), + myMinLen(0.0), + myLastU (0.0), + myFirstu (0.0) +{ + Initialize (theC, theAngularDeflection, theCurvatureDeflection, theMinimumOfPoints, theUTol, theMinLen); +} + +//======================================================================= +//function : GCPnts_TangentialDeflection +//purpose : +//======================================================================= +GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (const Adaptor3d_Curve& theC, + const Standard_Real theFirstParameter, const Standard_Real theLastParameter, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints, + const Standard_Real theUTol, + const Standard_Real theMinLen) +: myAngularDeflection (0.0), + myCurvatureDeflection (0.0), + myUTol (0.0), + myMinNbPnts (0), + myMinLen(0.0), + myLastU (0.0), + myFirstu (0.0) +{ + Initialize (theC, theFirstParameter, theLastParameter, + theAngularDeflection, theCurvatureDeflection, + theMinimumOfPoints, + theUTol, theMinLen); +} + +//======================================================================= +//function : GCPnts_TangentialDeflection +//purpose : +//======================================================================= +GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (const Adaptor2d_Curve2d& theC, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints, + const Standard_Real theUTol, + const Standard_Real theMinLen) +: myAngularDeflection (0.0), + myCurvatureDeflection (0.0), + myUTol (0.0), + myMinNbPnts (0), + myMinLen(0.0), + myLastU (0.0), + myFirstu (0.0) +{ + Initialize (theC, theAngularDeflection, theCurvatureDeflection, theMinimumOfPoints, theUTol, theMinLen); +} + +//======================================================================= +//function : GCPnts_TangentialDeflection +//purpose : +//======================================================================= +GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (const Adaptor2d_Curve2d& theC, + const Standard_Real theFirstParameter, const Standard_Real theLastParameter, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints, + const Standard_Real theUTol, + const Standard_Real theMinLen) +: myAngularDeflection (0.0), + myCurvatureDeflection (0.0), + myUTol (0.0), + myMinNbPnts (0), + myMinLen(0.0), + myLastU (0.0), + myFirstu (0.0) +{ + Initialize (theC, theFirstParameter, theLastParameter, + theAngularDeflection, theCurvatureDeflection, + theMinimumOfPoints, + theUTol, theMinLen); +} + +//======================================================================= +//function : Initialize +//purpose : +//======================================================================= +void GCPnts_TangentialDeflection::Initialize (const Adaptor3d_Curve& theC, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints, + const Standard_Real theUTol, + const Standard_Real theMinLen) +{ + Initialize (theC, theC.FirstParameter(), theC.LastParameter(), + theAngularDeflection, theCurvatureDeflection, + theMinimumOfPoints, + theUTol, theMinLen); +} + +//======================================================================= +//function : Initialize +//purpose : +//======================================================================= +void GCPnts_TangentialDeflection::Initialize (const Adaptor2d_Curve2d& theC, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints, + const Standard_Real theUTol, + const Standard_Real theMinLen) +{ + Initialize (theC, theC.FirstParameter(), theC.LastParameter(), + theAngularDeflection, theCurvatureDeflection, + theMinimumOfPoints, + theUTol, theMinLen); +} + +//======================================================================= +//function : Initialize +//purpose : +//======================================================================= +void GCPnts_TangentialDeflection::Initialize (const Adaptor3d_Curve& theC, + const Standard_Real theFirstParameter, const Standard_Real theLastParameter, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints, + const Standard_Real theUTol, + const Standard_Real theMinLen) +{ + initialize (theC, theFirstParameter, theLastParameter, + theAngularDeflection, theCurvatureDeflection, + theMinimumOfPoints, + theUTol, + theMinLen); +} + +//======================================================================= +//function : Initialize +//purpose : +//======================================================================= +void GCPnts_TangentialDeflection::Initialize (const Adaptor2d_Curve2d& theC, + const Standard_Real theFirstParameter, const Standard_Real theLastParameter, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints, + const Standard_Real theUTol, + const Standard_Real theMinLen) +{ + initialize (theC, theFirstParameter, theLastParameter, + theAngularDeflection, theCurvatureDeflection, + theMinimumOfPoints, + theUTol, + theMinLen); +} + +//======================================================================= +//function : EvaluateDu +//purpose : +//======================================================================= +template +void GCPnts_TangentialDeflection::EvaluateDu (const TheCurve& theC, + const Standard_Real theU, + gp_Pnt& theP, + Standard_Real& theDu, + Standard_Boolean& theNotDone) const +{ + gp_Vec T, N; + D2 (theC, theU, theP, T, N); + Standard_Real Lt = T.Magnitude(); + Standard_Real LTol = Precision::Confusion(); + if (Lt > LTol && N.Magnitude () > LTol) + { + Standard_Real Lc = N.CrossMagnitude (T); + Standard_Real Ln = Lc/Lt; + if (Ln > LTol) + { + theDu = sqrt (8.0 * Max (myCurvatureDeflection, myMinLen) / Ln); + theNotDone = Standard_False; + } + } +} + +//======================================================================= +//function : PerformLinear +//purpose : +//======================================================================= +template +void GCPnts_TangentialDeflection::PerformLinear (const TheCurve& theC) +{ + gp_Pnt P; + D0 (theC, myFirstu, P); + myParameters.Append (myFirstu); + myPoints .Append (P); + if (myMinNbPnts > 2) + { + Standard_Real Du = (myLastU - myFirstu) / myMinNbPnts; + Standard_Real U = myFirstu + Du; + for (Standard_Integer i = 2; i < myMinNbPnts; i++) + { + D0 (theC, U, P); + myParameters.Append (U); + myPoints .Append (P); + U += Du; + } + } + D0 (theC, myLastU, P); + myParameters.Append (myLastU); + myPoints .Append (P); +} + +//======================================================================= +//function : PerformCircular +//purpose : +//======================================================================= +template +void GCPnts_TangentialDeflection::PerformCircular (const TheCurve& theC) +{ + // akm 8/01/02 : check the radius before divide by it + Standard_Real dfR = theC.Circle().Radius(); + Standard_Real Du = GCPnts_TangentialDeflection::ArcAngularStep (dfR, myCurvatureDeflection, myAngularDeflection, myMinLen); + + const Standard_Real aDiff = myLastU - myFirstu; + // Round up number of points to satisfy curvatureDeflection more precisely + Standard_Integer NbPoints = (Standard_Integer)Min(Ceiling(aDiff / Du), 1.0e+6); + NbPoints = Max(NbPoints, myMinNbPnts - 1); + Du = aDiff / NbPoints; + + gp_Pnt P; + Standard_Real U = myFirstu; + for (Standard_Integer i = 1; i <= NbPoints; i++) + { + D0 (theC, U, P); + myParameters.Append (U); + myPoints .Append (P); + U += Du; + } + + D0 (theC, myLastU, P); + myParameters.Append (myLastU); + myPoints .Append (P); +} + +//======================================================================= +//function : initialize +//purpose : +//======================================================================= +template +void GCPnts_TangentialDeflection::initialize (const TheCurve& theC, + const Standard_Real theFirstParameter, const Standard_Real theLastParameter, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints, + const Standard_Real theUTol, + const Standard_Real theMinLen) +{ + Standard_ConstructionError_Raise_if (theCurvatureDeflection < Precision::Confusion() || theAngularDeflection < Precision::Angular(), + "GCPnts_TangentialDeflection::Initialize - Zero Deflection") + myParameters.Clear(); + myPoints .Clear(); + if (theFirstParameter < theLastParameter) + { + myFirstu = theFirstParameter; + myLastU = theLastParameter; } else { - return 0.; + myLastU = theFirstParameter; + myFirstu = theLastParameter; } -} + myUTol = theUTol; + myAngularDeflection = theAngularDeflection; + myCurvatureDeflection = theCurvatureDeflection; + myMinNbPnts = Max (theMinimumOfPoints, 2); + myMinLen = Max (theMinLen, Precision::Confusion()); - -// Return number of interval of continuity on which theParam is located. -// Last parameter is used to increase search speed. -static Standard_Integer getIntervalIdx(const Standard_Real theParam, - TColStd_Array1OfReal& theIntervs, - const Standard_Integer thePreviousIdx) -{ - Standard_Integer anIdx; - for(anIdx = thePreviousIdx; anIdx < theIntervs.Upper(); anIdx++) + switch (theC.GetType()) { - if (theParam >= theIntervs(anIdx) && - theParam <= theIntervs(anIdx + 1)) // Inside of anIdx interval. + case GeomAbs_Line: { + PerformLinear (theC); + break; + } + case GeomAbs_Circle: + { + PerformCircular (theC); + break; + } + case GeomAbs_BSplineCurve: + { + Handle(typename CurveTypes::BSplineCurve) aBS = theC.BSpline(); + if (aBS->NbPoles() == 2) PerformLinear (theC); + else PerformCurve (theC); + break; + } + case GeomAbs_BezierCurve: + { + Handle(typename CurveTypes::BezierCurve) aBZ = theC.Bezier(); + if (aBZ->NbPoles() == 2) PerformLinear (theC); + else PerformCurve (theC); + break; + } + default: + { + PerformCurve (theC); break; } } - return anIdx; -} -// -//======================================================================= -//function : CPnts_TangentialDeflection -//purpose : -//======================================================================= - -GCPnts_TangentialDeflection::GCPnts_TangentialDeflection () -: angularDeflection(0.0), - curvatureDeflection(0.0), - uTol(0.0), - minNbPnts(0), - myMinLen(0.0), - lastu(0.0), - firstu(0.0) -{ } //======================================================================= //function : AddPoint -//purpose : +//purpose : //======================================================================= - Standard_Integer GCPnts_TangentialDeflection::AddPoint (const gp_Pnt& thePnt, const Standard_Real theParam, @@ -126,38 +466,38 @@ Standard_Integer GCPnts_TangentialDeflection::AddPoint { const Standard_Real tol = Precision::PConfusion(); Standard_Integer index = -1; - const Standard_Integer nb = parameters.Length(); + const Standard_Integer nb = myParameters.Length(); for ( Standard_Integer i = 1; index == -1 && i <= nb; i++ ) { - Standard_Real dist = parameters.Value( i ) - theParam; + Standard_Real dist = myParameters.Value( i ) - theParam; if ( fabs( dist ) <= tol ) { index = i; if ( theIsReplace ) { - points.ChangeValue(i) = thePnt; - parameters.ChangeValue(i) = theParam; + myPoints .ChangeValue (i) = thePnt; + myParameters.ChangeValue (i) = theParam; } } else if ( dist > tol ) { - points.InsertBefore( i, thePnt ); - parameters.InsertBefore( i, theParam ); + myPoints .InsertBefore (i, thePnt); + myParameters.InsertBefore (i, theParam); index = i; } } - if ( index == -1 ) + if (index == -1) { - points.Append( thePnt ); - parameters.Append( theParam ); - index = parameters.Length(); + myPoints .Append (thePnt); + myParameters.Append (theParam); + index = myParameters.Length(); } return index; } //======================================================================= //function : ArcAngularStep -//purpose : +//purpose : //======================================================================= Standard_Real GCPnts_TangentialDeflection::ArcAngularStep( const Standard_Real theRadius, @@ -183,36 +523,485 @@ Standard_Real GCPnts_TangentialDeflection::ArcAngularStep( return Du; } -#include -#include -#include -#include +//======================================================================= +//function : PerformCurve +//purpose : +//======================================================================= +template +void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& theC) +{ + Standard_Integer i, j; + gp_XYZ V1, V2; + gp_Pnt MiddlePoint, CurrentPoint, LastPoint; + Standard_Real Du, Dusave, MiddleU, L1, L2; -#define TheCurve Adaptor3d_Curve -#define Handle_TheBezierCurve Handle(Geom_BezierCurve) -#define Handle_TheBSplineCurve Handle(Geom_BSplineCurve) -#define TheMaxCurvLinDist GCPnts_DistFunction -#define TheMaxCurvLinDistMV GCPnts_DistFunctionMV -#include "GCPnts_TangentialDeflection.pxx" -#undef Handle_TheBezierCurve -#undef Handle_TheBSplineCurve -#undef TheCurve -#undef TheMaxCurvLinDist -#undef TheMaxCurvLinDistMV + Standard_Real U1 = myFirstu; + Standard_Real LTol = Precision::Confusion(); // protection longueur nulle + Standard_Real ATol = 1.e-2 * myAngularDeflection; + if (ATol > 1.e-2) + { + ATol = 1.e-2; + } + else if (ATol < 1.e-7) + { + ATol = 1.e-7; + } + D0 (theC, myLastU, LastPoint); -#include -#include -#include -#include -#define TheCurve Adaptor2d_Curve2d -#define Handle_TheBezierCurve Handle(Geom2d_BezierCurve) -#define Handle_TheBSplineCurve Handle(Geom2d_BSplineCurve) -#define TheMaxCurvLinDist GCPnts_DistFunction2d -#define TheMaxCurvLinDistMV GCPnts_DistFunction2dMV -#include "GCPnts_TangentialDeflection.pxx" -#undef Handle_TheBezierCurve -#undef Handle_TheBSplineCurve -#undef TheCurve -#undef TheMaxCurvLinDist -#undef TheMaxCurvLinDistMV + // Initialization du calcul + + Standard_Boolean NotDone = Standard_True; + Dusave = (myLastU - myFirstu) * Us3; + Du = Dusave; + EvaluateDu (theC, U1, CurrentPoint, Du, NotDone); + myParameters.Append (U1); + myPoints .Append (CurrentPoint); + + // Used to detect "isLine" current bspline and in Du computation in general handling. + const Standard_Integer NbInterv = theC.NbIntervals (GeomAbs_CN); + TColStd_Array1OfReal Intervs (1, NbInterv + 1); + theC.Intervals (Intervs, GeomAbs_CN); + + if (NotDone || Du > 5. * Dusave) + { + //C'est soit une droite, soit une singularite : + V1 = (LastPoint.XYZ() - CurrentPoint.XYZ()); + L1 = V1.Modulus (); + if (L1 > LTol) + { + // Si c'est une droite on verifie en calculant minNbPoints : + Standard_Boolean IsLine = Standard_True; + Standard_Integer NbPoints = (myMinNbPnts > 3) ? myMinNbPnts : 3; + switch (theC.GetType()) + { + case GeomAbs_BSplineCurve: + { + Handle(typename CurveTypes::BSplineCurve) BS = theC.BSpline(); + NbPoints = Max(BS->Degree() + 1, NbPoints); + break; + } + case GeomAbs_BezierCurve: + { + Handle(typename CurveTypes::BezierCurve) BZ = theC.Bezier(); + NbPoints = Max(BZ->Degree() + 1, NbPoints); + break; + } + default: + { + break; + } + } + //// + Standard_Real param = 0.; + for (i = 1; i <= NbInterv && IsLine; ++i) + { + // Avoid usage intervals out of [myFirstu, myLastU]. + if ((Intervs(i + 1) < myFirstu) + || (Intervs(i) > myLastU)) + { + continue; + } + + // Fix border points in applicable intervals, to avoid be out of target interval. + if ((Intervs(i) < myFirstu) + && (Intervs(i+1) > myFirstu)) + { + Intervs(i) = myFirstu; + } + if ((Intervs(i) < myLastU) + && (Intervs(i+1) > myLastU)) + { + Intervs(i + 1) = myLastU; + } + + const Standard_Real delta = (Intervs(i+1) - Intervs(i))/NbPoints; + for (j = 1; j <= NbPoints && IsLine; ++j) + { + param = Intervs(i) + j*delta; + D0 (theC, param, MiddlePoint); + V2 = (MiddlePoint.XYZ() - CurrentPoint.XYZ()); + L2 = V2.Modulus (); + if (L2 > LTol) + { + const Standard_Real aAngle = V2.CrossMagnitude(V1)/(L1*L2); + IsLine = (aAngle < ATol); + } + } + } + + if (IsLine) + { + myParameters.Clear(); + myPoints .Clear(); + + PerformLinear (theC); + return; + } + else + { + // c'etait une singularite on continue: + //Du = Dusave; + EvaluateDu (theC, param, MiddlePoint, Du, NotDone); + } + } + else + { + Du = (myLastU - myFirstu) / 2.1; + MiddleU = myFirstu + Du; + D0 (theC, MiddleU, MiddlePoint); + V1 = (MiddlePoint.XYZ() - CurrentPoint.XYZ()); + L1 = V1.Modulus (); + if (L1 < LTol) + { + // L1 < LTol C'est une courbe de longueur nulle, calcul termine : + // on renvoi un segment de 2 points (protection) + myParameters.Append (myLastU); + myPoints .Append (LastPoint); + return; + } + } + } + + if (Du > Dusave) Du = Dusave; + else Dusave = Du; + + if (Du < myUTol) + { + Du = myLastU - myFirstu; + if (Du < myUTol) + { + myParameters.Append (myLastU); + myPoints .Append (LastPoint); + return; + } + } + + // Traitement normal pour une courbe + Standard_Boolean MorePoints = Standard_True; + Standard_Real U2 = myFirstu; + Standard_Real AngleMax = myAngularDeflection * 0.5; // car on prend le point milieu + // Indexes of intervals of U1 and U2, used to handle non-uniform case. + Standard_Integer aIdx[2] = {Intervs.Lower(), Intervs.Lower()}; + Standard_Boolean isNeedToCheck = Standard_False; + gp_Pnt aPrevPoint = myPoints.Last(); + + while (MorePoints) + { + aIdx[0] = getIntervalIdx(U1, Intervs, aIdx[0]); + U2 += Du; + + if (U2 >= myLastU) // Bout de courbe + { + U2 = myLastU; + CurrentPoint = LastPoint; + Du = U2-U1; + Dusave = Du; + } + else + { + D0 (theC, U2, CurrentPoint); // Point suivant + } + + Standard_Real Coef = 0.0, ACoef = 0., FCoef = 0.; + Standard_Boolean Correction, TooLarge, TooSmall; + TooLarge = Standard_False; + Correction = Standard_True; + TooSmall = Standard_False; + + while (Correction) // Ajustement Du + { + if (isNeedToCheck) + { + aIdx[1] = getIntervalIdx(U2, Intervs, aIdx[0]); + if (aIdx[1] > aIdx[0]) // Jump to another polynom. + { + // Set Du to the smallest value and check deflection on it. + if (Du > (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3) + { + Du = (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3; + U2 = U1 + Du; + if (U2 > myLastU) + { + U2 = myLastU; + } + D0 (theC, U2, CurrentPoint); + } + } + } + MiddleU = (U1+U2)*0.5; // Verif / au point milieu + D0 (theC, MiddleU, MiddlePoint); + + V1 = (CurrentPoint.XYZ() - aPrevPoint.XYZ()); // Critere de fleche + V2 = (MiddlePoint.XYZ() - aPrevPoint.XYZ()); + L1 = V1.Modulus (); + + FCoef = (L1 > myMinLen) ? V1.CrossMagnitude(V2) / (L1 * myCurvatureDeflection) : 0.0; + + V1 = (CurrentPoint.XYZ() - MiddlePoint.XYZ()); // Critere d'angle + L1 = V1.Modulus (); + L2 = V2.Modulus (); + if (L1 > myMinLen && L2 > myMinLen) + { + Standard_Real angg = V1.CrossMagnitude(V2) / (L1 * L2); + ACoef = angg / AngleMax; + } + else + { + ACoef = 0.0; + } + + // On retient le plus penalisant + Coef = Max (ACoef, FCoef); + + if (isNeedToCheck && Coef < 0.55) + { + isNeedToCheck = Standard_False; + Du = Dusave; + U2 = U1 + Du; + if (U2 > myLastU) + { + U2 = myLastU; + } + D0 (theC, U2, CurrentPoint); + continue; + } + + if (Coef <= 1.0) + { + if (Abs (myLastU - U2) < myUTol) + { + myParameters.Append (myLastU); + myPoints .Append (LastPoint); + MorePoints = Standard_False; + Correction = Standard_False; + } + else + { + if (Coef >= 0.55 || TooLarge) + { + myParameters.Append (U2); + myPoints .Append (CurrentPoint); + aPrevPoint = CurrentPoint; + Correction = Standard_False; + isNeedToCheck = Standard_True; + } + else if (TooSmall) + { + Correction = Standard_False; + aPrevPoint = CurrentPoint; + } + else + { + TooSmall = Standard_True; + //Standard_Real UUU2 = U2; + Du += Min((U2-U1)*(1.-Coef), Du*Us3); + + U2 = U1 + Du; + if (U2 > myLastU) + { + U2 = myLastU; + } + D0 (theC, U2, CurrentPoint); + } + } + } + else + { + if (Coef >= 1.5) + { + if (!aPrevPoint.IsEqual (myPoints.Last(), Precision::Confusion())) + { + myParameters.Append (U1); + myPoints .Append (aPrevPoint); + } + U2 = MiddleU; + Du = U2-U1; + CurrentPoint = MiddlePoint; + } + else + { + Du*=0.9; + U2 = U1 + Du; + D0 (theC, U2, CurrentPoint); + TooLarge = Standard_True; + } + } + } + + Du = U2-U1; + + if (MorePoints) + { + if (U1 > myFirstu) + { + if (FCoef > ACoef) + { + // La fleche est critere de decoupage + EvaluateDu (theC, U2, CurrentPoint, Du, NotDone); + if (NotDone) + { + Du += (Du-Dusave)*(Du/Dusave); + if (Du > 1.5 * Dusave) Du = 1.5 * Dusave; + if (Du < 0.75* Dusave) Du = 0.75 * Dusave; + } + } + else + { + //L'angle est le critere de decoupage + Du += (Du-Dusave)*(Du/Dusave); + if (Du > 1.5 * Dusave) Du = 1.5 * Dusave; + if (Du < 0.75* Dusave) Du = 0.75 * Dusave; + } + } + + if (Du < myUTol) + { + Du = myLastU - U2; + if (Du < myUTol) + { + myParameters.Append (myLastU); + myPoints .Append (LastPoint); + MorePoints = Standard_False; + } + else if (Du*Us3 > myUTol) + { + Du*=Us3; + } + } + U1 = U2; + Dusave = Du; + } + } + // Recalage avant dernier point : + i = myPoints.Length() - 1; + // Real d = myPoints (i).Distance (myPoints (i+1)); + // if (Abs(myParameters (i) - myParameters (i+1))<= 0.000001 || d < Precision::Confusion()) { + // cout<<"deux points confondus"<= 2) + { + MiddleU = myParameters (i-1); + MiddleU = (myLastU + MiddleU)*0.5; + D0 (theC, MiddleU, MiddlePoint); + myParameters.SetValue (i, MiddleU); + myPoints .SetValue (i, MiddlePoint); + } + + //-- On rajoute des points aux milieux des segments si le nombre + //-- mini de points n'est pas atteint + //-- + Standard_Integer Nbp = myPoints.Length(); + + //std::cout << "GCPnts_TangentialDeflection: Number of Points (" << Nbp << " " << myMinNbPnts << " )" << std::endl; + + while (Nbp < myMinNbPnts) + { + for (i = 2; i <= Nbp; i += 2) + { + MiddleU = (myParameters.Value(i-1) + myParameters.Value(i)) * 0.5; + D0 (theC, MiddleU, MiddlePoint); + myParameters.InsertBefore (i, MiddleU); + myPoints .InsertBefore (i, MiddlePoint); + Nbp++; + } + } + // Additional check for intervals + const Standard_Real MinLen2 = myMinLen * myMinLen; + const Standard_Integer MaxNbp = 10 * Nbp; + for (i = 1; i < Nbp; ++i) + { + U1 = myParameters (i); + U2 = myParameters (i + 1); + + if (U2 - U1 <= myUTol) + { + continue; + } + + // Check maximal deflection on interval; + Standard_Real dmax = 0.; + Standard_Real umax = 0.; + Standard_Real amax = 0.; + EstimDefl (theC, U1, U2, dmax, umax); + const gp_Pnt& P1 = myPoints (i); + const gp_Pnt& P2 = myPoints (i + 1); + D0 (theC, umax, MiddlePoint); + amax = EstimAngl(P1, MiddlePoint, P2); + if (dmax > myCurvatureDeflection || amax > AngleMax) + { + if (umax - U1 > myUTol && U2 - umax > myUTol) + { + if (P1.SquareDistance(MiddlePoint) > MinLen2 + && P2.SquareDistance(MiddlePoint) > MinLen2) + { + myParameters.InsertAfter (i, umax); + myPoints .InsertAfter (i, MiddlePoint); + ++Nbp; + --i; //To compensate ++i in loop header: i must point to first part of split interval + if (Nbp > MaxNbp) + { + break; + } + } + } + } + } +} + +//======================================================================= +//function : EstimDefl +//purpose : +//======================================================================= +template +void GCPnts_TangentialDeflection::EstimDefl (const TheCurve& theC, + const Standard_Real theU1, const Standard_Real theU2, + Standard_Real& theMaxDefl, Standard_Real& theUMax) +{ + const Standard_Real Du = (myLastU - myFirstu); + // + typename CurveTypes::DistFunction aFunc (theC, theU1, theU2); + // + const Standard_Integer aNbIter = 100; + const Standard_Real aRelTol = Max (1.e-3, 2. * myUTol / (Abs(theU1) + Abs(theU2))); + // + math_BrentMinimum anOptLoc (aRelTol, aNbIter, myUTol); + anOptLoc.Perform (aFunc, theU1, (theU1 + theU2) / 2., theU2); + if (anOptLoc.IsDone()) + { + theMaxDefl = Sqrt(-anOptLoc.Minimum()); + theUMax = anOptLoc.Location(); + return; + } + // + math_Vector aLowBorder (1, 1), aUppBorder (1, 1), aSteps (1, 1); + aSteps (1) = Max (0.1 * Du, 100. * myUTol); + const Standard_Integer aNbParticles = Max(8, RealToInt(32 * (theU2 - theU1) / Du)); + aLowBorder (1) = theU1; + aUppBorder (1) = theU2; + // + // + Standard_Real aValue = 0.0; + math_Vector aT (1, 1); + typename CurveTypes::DistFunctionMV aFuncMV(aFunc); + + math_PSO aFinder (&aFuncMV, aLowBorder, aUppBorder, aSteps, aNbParticles); + aFinder.Perform (aSteps, aValue, aT); + // + anOptLoc.Perform (aFunc, + Max (aT(1) - aSteps(1), theU1), + aT (1), + Min (aT(1) + aSteps(1), theU2)); + if (anOptLoc.IsDone()) + { + theMaxDefl = Sqrt(-anOptLoc.Minimum()); + theUMax = anOptLoc.Location(); + return; + } + + theMaxDefl = Sqrt(-aValue); + theUMax = aT (1); +} diff --git a/src/GCPnts/GCPnts_TangentialDeflection.hxx b/src/GCPnts/GCPnts_TangentialDeflection.hxx index 0e66213137..b84917835c 100644 --- a/src/GCPnts/GCPnts_TangentialDeflection.hxx +++ b/src/GCPnts/GCPnts_TangentialDeflection.hxx @@ -26,16 +26,13 @@ #include #include -class Standard_ConstructionError; -class Standard_OutOfRange; - //! Computes a set of points on a curve from package //! Adaptor3d such as between two successive points //! P1(u1)and P2(u2) : -//! +//! @code //! . ||P1P3^P3P2||/||P1P3||*||P3P2|| gp::Resolution() must be //! satisfied at the construction time. //! -//! A minimum number of points can be fixed for a -//! linear or circular element. +//! A minimum number of points can be fixed for a linear or circular element. //! Example: -//! Handle(Geom_BezierCurve) C = new Geom_BezierCurve (Poles); -//! GeomAdaptor_Curve Curve (C); -//! Real CDeflect = 0.01; //Curvature deflection -//! Real ADeflect = 0.09; //Angular deflection +//! @code +//! Handle(Geom_BezierCurve) aCurve = new Geom_BezierCurve (thePoles); +//! GeomAdaptor_Curve aCurveAdaptor (aCurve); +//! double aCDeflect = 0.01; // Curvature deflection +//! double anADeflect = 0.09; // Angular deflection //! -//! GCPnts_TangentialDeflection PointsOnCurve; -//! PointsOnCurve.Initialize (Curve, ADeflect, CDeflect); -//! -//! Real U; -//! gp_Pnt P; -//! for (Integer i=1; i<=PointsOnCurve.NbPoints();i++) { -//! U = PointsOnCurve.Parameter (i); -//! P = PointsOnCurve.Value (i); +//! GCPnts_TangentialDeflection aPointsOnCurve; +//! aPointsOnCurve.Initialize (aCurveAdaptor, anADeflect, aCDeflect); +//! for (int i = 1; i <= aPointsOnCurve.NbPoints(); ++i) +//! { +//! double aU = aPointsOnCurve.Parameter (i); +//! gp_Pnt aPnt = aPointsOnCurve.Value (i); //! } - +//! @endcode class GCPnts_TangentialDeflection { public: -// DEFINE_STANDARD_ALLOC - + //! Empty constructor. + //! @sa Initialize() Standard_EXPORT GCPnts_TangentialDeflection(); - - Standard_EXPORT GCPnts_TangentialDeflection(const Adaptor3d_Curve& C, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7); - - Standard_EXPORT GCPnts_TangentialDeflection(const Adaptor3d_Curve& C, const Standard_Real FirstParameter, const Standard_Real LastParameter, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7); - - Standard_EXPORT GCPnts_TangentialDeflection(const Adaptor2d_Curve2d& C, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7); - - Standard_EXPORT GCPnts_TangentialDeflection(const Adaptor2d_Curve2d& C, const Standard_Real FirstParameter, const Standard_Real LastParameter, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7); - - Standard_EXPORT void Initialize (const Adaptor3d_Curve& C, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7); - - Standard_EXPORT void Initialize (const Adaptor3d_Curve& C, const Standard_Real FirstParameter, const Standard_Real LastParameter, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7); - - Standard_EXPORT void Initialize (const Adaptor2d_Curve2d& C, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7); - - Standard_EXPORT void Initialize (const Adaptor2d_Curve2d& C, const Standard_Real FirstParameter, const Standard_Real LastParameter, const Standard_Real AngularDeflection, const Standard_Real CurvatureDeflection, const Standard_Integer MinimumOfPoints = 2, const Standard_Real UTol = 1.0e-9, const Standard_Real theMinLen = 1.0e-7); - + + //! Constructor for 3D curve. + //! @param theC [in] 3d curve + //! @param theAngularDeflection [in] angular deflection in radians + //! @param theCurvatureDeflection [in] linear deflection + //! @param theMinimumOfPoints [in] minimum number of points + //! @param theUTol [in] tolerance in curve parametric scope + //! @param theMinLen [in] minimal length + Standard_EXPORT GCPnts_TangentialDeflection (const Adaptor3d_Curve& theC, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints = 2, + const Standard_Real theUTol = 1.0e-9, + const Standard_Real theMinLen = 1.0e-7); + + //! Constructor for 3D curve with restricted range. + //! @param theC [in] 3d curve + //! @param theFirstParameter [in] first parameter on curve + //! @param theLastParameter [in] last parameter on curve + //! @param theAngularDeflection [in] angular deflection in radians + //! @param theCurvatureDeflection [in] linear deflection + //! @param theMinimumOfPoints [in] minimum number of points + //! @param theUTo l [in] tolerance in curve parametric scope + //! @param theMinLen [in] minimal length + Standard_EXPORT GCPnts_TangentialDeflection (const Adaptor3d_Curve& theC, + const Standard_Real theFirstParameter, const Standard_Real theLastParameter, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints = 2, + const Standard_Real theUTol = 1.0e-9, + const Standard_Real theMinLen = 1.0e-7); + + //! Constructor for 2D curve. + //! @param theC [in] 2d curve + //! @param theAngularDeflection [in] angular deflection in radians + //! @param theCurvatureDeflection [in] linear deflection + //! @param theMinimumOfPoints [in] minimum number of points + //! @param theUTol [in] tolerance in curve parametric scope + //! @param theMinLen [in] minimal length + Standard_EXPORT GCPnts_TangentialDeflection (const Adaptor2d_Curve2d& theC, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints = 2, + const Standard_Real theUTol = 1.0e-9, + const Standard_Real theMinLen = 1.0e-7); + + //! Constructor for 2D curve with restricted range. + //! @param theC [in] 2d curve + //! @param theFirstParameter [in] first parameter on curve + //! @param theLastParameter [in] last parameter on curve + //! @param theAngularDeflection [in] angular deflection in radians + //! @param theCurvatureDeflection [in] linear deflection + //! @param theMinimumOfPoints [in] minimum number of points + //! @param theUTol [in] tolerance in curve parametric scope + //! @param theMinLen [in] minimal length + Standard_EXPORT GCPnts_TangentialDeflection (const Adaptor2d_Curve2d& theC, + const Standard_Real theFirstParameter, const Standard_Real theLastParameter, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints = 2, + const Standard_Real theUTol = 1.0e-9, + const Standard_Real theMinLen = 1.0e-7); + + //! Initialize algorithm for 3D curve. + //! @param theC [in] 3d curve + //! @param theAngularDeflection [in] angular deflection in radians + //! @param theCurvatureDeflection [in] linear deflection + //! @param theMinimumOfPoints [in] minimum number of points + //! @param theUTol [in] tolerance in curve parametric scope + //! @param theMinLen [in] minimal length + Standard_EXPORT void Initialize (const Adaptor3d_Curve& theC, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints = 2, + const Standard_Real theUTol = 1.0e-9, + const Standard_Real theMinLen = 1.0e-7); + + //! Initialize algorithm for 3D curve with restricted range. + //! @param theC [in] 3d curve + //! @param theFirstParameter [in] first parameter on curve + //! @param theLastParameter [in] last parameter on curve + //! @param theAngularDeflection [in] angular deflection in radians + //! @param theCurvatureDeflection [in] linear deflection + //! @param theMinimumOfPoints [in] minimum number of points + //! @param theUTol [in] tolerance in curve parametric scope + //! @param theMinLen [in] minimal length + Standard_EXPORT void Initialize (const Adaptor3d_Curve& theC, + const Standard_Real theFirstParameter, const Standard_Real theLastParameter, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints = 2, + const Standard_Real theUTol = 1.0e-9, + const Standard_Real theMinLen = 1.0e-7); + + //! Initialize algorithm for 2D curve. + //! @param theC [in] 2d curve + //! @param theAngularDeflection [in] angular deflection in radians + //! @param theCurvatureDeflection [in] linear deflection + //! @param theMinimumOfPoints [in] minimum number of points + //! @param theUTol [in] tolerance in curve parametric scope + //! @param theMinLen [in] minimal length + Standard_EXPORT void Initialize (const Adaptor2d_Curve2d& theC, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints = 2, + const Standard_Real theUTol = 1.0e-9, + const Standard_Real theMinLen = 1.0e-7); + + //! Initialize algorithm for 2D curve with restricted range. + //! @param theC [in] 2d curve + //! @param theFirstParameter [in] first parameter on curve + //! @param theLastParameter [in] last parameter on curve + //! @param theAngularDeflection [in] angular deflection in radians + //! @param theCurvatureDeflection [in] linear deflection + //! @param theMinimumOfPoints [in] minimum number of points + //! @param theUTol [in] tolerance in curve parametric scope + //! @param theMinLen [in] minimal length + Standard_EXPORT void Initialize (const Adaptor2d_Curve2d& theC, + const Standard_Real theFirstParameter, const Standard_Real theLastParameter, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints = 2, + const Standard_Real theUTol = 1.0e-9, + const Standard_Real theMinLen = 1.0e-7); + //! Add point to already calculated points (or replace existing) //! Returns index of new added point //! or founded with parametric tolerance (replaced if theIsReplace is true) - Standard_EXPORT Standard_Integer AddPoint (const gp_Pnt& thePnt, const Standard_Real theParam, const Standard_Boolean theIsReplace = Standard_True); - + Standard_EXPORT Standard_Integer AddPoint (const gp_Pnt& thePnt, + const Standard_Real theParam, + const Standard_Boolean theIsReplace = Standard_True); + Standard_Integer NbPoints () const { - return parameters.Length (); + return myParameters.Length (); } - + Standard_Real Parameter (const Standard_Integer I) const { - return parameters.Value (I); + return myParameters.Value (I); } - + gp_Pnt Value (const Standard_Integer I) const { - return points.Value (I); + return myPoints.Value (I); } - + //! Computes angular step for the arc using the given parameters. - Standard_EXPORT static Standard_Real ArcAngularStep (const Standard_Real theRadius, const Standard_Real theLinearDeflection, const Standard_Real theAngularDeflection, const Standard_Real theMinLength); + Standard_EXPORT static Standard_Real ArcAngularStep (const Standard_Real theRadius, + const Standard_Real theLinearDeflection, + const Standard_Real theAngularDeflection, + const Standard_Real theMinLength); private: - - Standard_EXPORT void PerformLinear (const Adaptor3d_Curve& C); - - Standard_EXPORT void PerformLinear (const Adaptor2d_Curve2d& C); - - Standard_EXPORT void PerformCircular (const Adaptor3d_Curve& C); - - Standard_EXPORT void PerformCircular (const Adaptor2d_Curve2d& C); - - Standard_EXPORT void PerformCurve (const Adaptor3d_Curve& C); - - Standard_EXPORT void PerformCurve (const Adaptor2d_Curve2d& C); - - Standard_EXPORT void EvaluateDu (const Adaptor3d_Curve& C, const Standard_Real U, gp_Pnt& P, Standard_Real& Du, Standard_Boolean& NotDone) const; - - Standard_EXPORT void EvaluateDu (const Adaptor2d_Curve2d& C, const Standard_Real U, gp_Pnt& P, Standard_Real& Du, Standard_Boolean& NotDone) const; - Standard_EXPORT void EstimDefl (const Adaptor3d_Curve& C, const Standard_Real U1, const Standard_Real U2, - Standard_Real& MaxDefl, Standard_Real& UMax); - - Standard_EXPORT void EstimDefl (const Adaptor2d_Curve2d& C, const Standard_Real U1, const Standard_Real U2, - Standard_Real& MaxDefl, Standard_Real& UMax); + template + void initialize (const TheCurve& theC, + const Standard_Real theFirstParameter, const Standard_Real theLastParameter, + const Standard_Real theAngularDeflection, const Standard_Real theCurvatureDeflection, + const Standard_Integer theMinimumOfPoints, + const Standard_Real theUTol, + const Standard_Real theMinLen); + + template + void PerformLinear (const TheCurve& theC); + + template + void PerformCircular (const TheCurve& theC); + + //! Respecting the angle and the deflection, + //! we impose a minimum number of points on a curve. + template + void PerformCurve (const TheCurve& theC); + + template + void EvaluateDu (const TheCurve& theC, + const Standard_Real theU, + gp_Pnt& theP, + Standard_Real& theDu, + Standard_Boolean& theNotDone) const; + + //! Estimation of maximal deflection for interval [theU1, theU2] + template + void EstimDefl (const TheCurve& theC, + const Standard_Real theU1, const Standard_Real theU2, + Standard_Real& theMaxDefl, Standard_Real& theUMax); private: - Standard_Real angularDeflection; - Standard_Real curvatureDeflection; - Standard_Real uTol; - Standard_Integer minNbPnts; + Standard_Real myAngularDeflection; + Standard_Real myCurvatureDeflection; + Standard_Real myUTol; + Standard_Integer myMinNbPnts; Standard_Real myMinLen; - Standard_Real lastu; - Standard_Real firstu; - TColgp_SequenceOfPnt points; - TColStd_SequenceOfReal parameters; + Standard_Real myLastU; + Standard_Real myFirstu; + TColgp_SequenceOfPnt myPoints; + TColStd_SequenceOfReal myParameters; }; #endif // _GCPnts_TangentialDeflection_HeaderFile diff --git a/src/GCPnts/GCPnts_TangentialDeflection.pxx b/src/GCPnts/GCPnts_TangentialDeflection.pxx deleted file mode 100644 index 5cc7038773..0000000000 --- a/src/GCPnts/GCPnts_TangentialDeflection.pxx +++ /dev/null @@ -1,687 +0,0 @@ -// Created on: 1996-11-08 -// Created by: Jean Claude VAUTHIER -// Copyright (c) 1996-1999 Matra Datavision -// Copyright (c) 1999-2014 OPEN CASCADE SAS -// -// This file is part of Open CASCADE Technology software library. -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License version 2.1 as published -// by the Free Software Foundation, with special exception defined in the file -// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT -// distribution for complete text of the license and disclaimer of any warranty. -// -// Alternatively, this file may be used under the terms of Open CASCADE -// commercial license or contractual agreement. - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#define Us3 0.3333333333333333333333333333 - -void GCPnts_TangentialDeflection::EvaluateDu ( - const TheCurve& C, - const Standard_Real U, - gp_Pnt& P, - Standard_Real& Du, - Standard_Boolean& NotDone) const { - - gp_Vec T, N; - D2 (C, U, P, T, N); - Standard_Real Lt = T.Magnitude (); - Standard_Real LTol = Precision::Confusion (); - if (Lt > LTol && N.Magnitude () > LTol) { - Standard_Real Lc = N.CrossMagnitude (T); - Standard_Real Ln = Lc/Lt; - if (Ln > LTol) { - Du = sqrt (8.0 * Max(curvatureDeflection, myMinLen) / Ln); - NotDone = Standard_False; - } - } -} - - -//======================================================================= -//function : GCPnts_TangentialDeflection -//purpose : -//======================================================================= - -GCPnts_TangentialDeflection::GCPnts_TangentialDeflection ( - const TheCurve& C, - const Standard_Real AngularDeflection, - const Standard_Real CurvatureDeflection, - const Standard_Integer MinimumOfPoints, - const Standard_Real UTol, - const Standard_Real theMinLen) - -{ - Initialize (C,AngularDeflection,CurvatureDeflection,MinimumOfPoints,UTol,theMinLen); -} - - -//======================================================================= -//function : GCPnts_TangentialDeflection -//purpose : -//======================================================================= - -GCPnts_TangentialDeflection::GCPnts_TangentialDeflection ( - const TheCurve& C, - const Standard_Real FirstParameter, - const Standard_Real LastParameter, - const Standard_Real AngularDeflection, - const Standard_Real CurvatureDeflection, - const Standard_Integer MinimumOfPoints, - const Standard_Real UTol, - const Standard_Real theMinLen) - -{ - Initialize (C, - FirstParameter, - LastParameter, - AngularDeflection, - CurvatureDeflection, - MinimumOfPoints, - UTol, theMinLen); -} - - - -//======================================================================= -//function : Initialize -//purpose : -//======================================================================= - -void GCPnts_TangentialDeflection::Initialize ( - const TheCurve& C, - const Standard_Real AngularDeflection, - const Standard_Real CurvatureDeflection, - const Standard_Integer MinimumOfPoints, - const Standard_Real UTol, - const Standard_Real theMinLen) - -{ - Initialize (C, - C.FirstParameter (), - C.LastParameter (), - AngularDeflection, - CurvatureDeflection, - MinimumOfPoints, - UTol, theMinLen); -} - - -//======================================================================= -//function : Initialize -//purpose : -//======================================================================= - -void GCPnts_TangentialDeflection::Initialize ( - const TheCurve& C, - const Standard_Real FirstParameter, - const Standard_Real LastParameter, - const Standard_Real AngularDeflection, - const Standard_Real CurvatureDeflection, - const Standard_Integer MinimumOfPoints, - const Standard_Real UTol, - const Standard_Real theMinLen) - -{ - - Standard_ConstructionError_Raise_if (CurvatureDeflection < Precision::Confusion () || - AngularDeflection < Precision::Angular (), - "GCPnts_TangentialDeflection::Initialize - Zero Deflection") - - parameters.Clear (); - points .Clear (); - if (FirstParameter < LastParameter) { - firstu = FirstParameter; - lastu = LastParameter; - } - else { - lastu = FirstParameter; - firstu = LastParameter; - } - uTol = UTol; - angularDeflection = AngularDeflection; - curvatureDeflection = CurvatureDeflection; - minNbPnts = Max (MinimumOfPoints, 2); - myMinLen = Max(theMinLen, Precision::Confusion()); - - switch (C.GetType()) { - - case GeomAbs_Line: - PerformLinear (C); - break; - - case GeomAbs_Circle: - PerformCircular (C); - break; - - case GeomAbs_BSplineCurve: - { - Handle_TheBSplineCurve BS = C.BSpline() ; - if (BS->NbPoles() == 2 ) PerformLinear (C); - else PerformCurve (C); - break; - } - case GeomAbs_BezierCurve: - { - Handle_TheBezierCurve BZ = C.Bezier(); - if (BZ->NbPoles() == 2) PerformLinear (C); - else PerformCurve (C); - break; - } - default: PerformCurve (C); - - } -} - - -//======================================================================= -//function : PerformLinear -//purpose : -//======================================================================= - -void GCPnts_TangentialDeflection::PerformLinear (const TheCurve& C) { - - gp_Pnt P; - D0 (C, firstu, P); - parameters.Append (firstu); - points .Append (P); - if (minNbPnts > 2) - { - Standard_Real Du = (lastu - firstu) / minNbPnts; - Standard_Real U = firstu + Du; - for (Standard_Integer i = 2; i < minNbPnts; i++) - { - D0 (C, U, P); - parameters.Append (U); - points .Append (P); - U += Du; - } - } - D0 (C, lastu, P); - parameters.Append (lastu); - points .Append (P); -} - -//======================================================================= -//function : PerformCircular -//purpose : -//======================================================================= - -void GCPnts_TangentialDeflection::PerformCircular (const TheCurve& C) -{ - // akm 8/01/02 : check the radius before divide by it - Standard_Real dfR = C.Circle().Radius(); - Standard_Real Du = GCPnts_TangentialDeflection::ArcAngularStep( - dfR, curvatureDeflection, angularDeflection, myMinLen); - - const Standard_Real aDiff = lastu - firstu; - // Round up number of points to satisfy curvatureDeflection more precisely - Standard_Integer NbPoints = (Standard_Integer)Min(Ceiling(aDiff / Du), 1.0e+6); - NbPoints = Max(NbPoints, minNbPnts - 1); - Du = aDiff / NbPoints; - - gp_Pnt P; - Standard_Real U = firstu; - for (Standard_Integer i = 1; i <= NbPoints; i++) - { - D0 (C, U, P); - parameters.Append (U); - points .Append (P); - U += Du; - } - D0 (C, lastu, P); - parameters.Append (lastu); - points .Append (P); -} - - -//======================================================================= -//function : PerformCurve -//purpose : On respecte ll'angle et la fleche, on peut imposer un nombre -// minimum de points sur un element lineaire -//======================================================================= -void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C) - -{ - Standard_Integer i, j; - gp_XYZ V1, V2; - gp_Pnt MiddlePoint, CurrentPoint, LastPoint; - Standard_Real Du, Dusave, MiddleU, L1, L2; - - Standard_Real U1 = firstu; - Standard_Real LTol = Precision::Confusion (); //protection longueur nulle - Standard_Real ATol = 1.e-2 * angularDeflection; - if(ATol > 1.e-2) - ATol = 1.e-2; - else if(ATol < 1.e-7) - ATol = 1.e-7; - - D0 (C, lastu, LastPoint); - - //Initialization du calcul - - Standard_Boolean NotDone = Standard_True; - Dusave = (lastu - firstu)*Us3; - Du = Dusave; - EvaluateDu (C, U1, CurrentPoint, Du, NotDone); - parameters.Append (U1); - points .Append (CurrentPoint); - - // Used to detect "isLine" current bspline and in Du computation in general handling. - Standard_Integer NbInterv = C.NbIntervals(GeomAbs_CN); - TColStd_Array1OfReal Intervs(1, NbInterv+1); - C.Intervals(Intervs, GeomAbs_CN); - - if (NotDone || Du > 5. * Dusave) { - //C'est soit une droite, soit une singularite : - V1 = (LastPoint.XYZ() - CurrentPoint.XYZ()); - L1 = V1.Modulus (); - if (L1 > LTol) - { - //Si c'est une droite on verifie en calculant minNbPoints : - Standard_Boolean IsLine = Standard_True; - Standard_Integer NbPoints = (minNbPnts > 3) ? minNbPnts : 3; - switch (C.GetType()) { - case GeomAbs_BSplineCurve: - { - Handle_TheBSplineCurve BS = C.BSpline() ; - NbPoints = Max(BS->Degree() + 1, NbPoints); - break; - } - case GeomAbs_BezierCurve: - { - Handle_TheBezierCurve BZ = C.Bezier(); - NbPoints = Max(BZ->Degree() + 1, NbPoints); - break; - } - default: - ;} - //// - Standard_Real param = 0.; - for (i = 1; i <= NbInterv && IsLine; ++i) - { - // Avoid usage intervals out of [firstu, lastu]. - if ((Intervs(i+1) < firstu) || - (Intervs(i) > lastu)) - { - continue; - } - // Fix border points in applicable intervals, to avoid be out of target interval. - if ((Intervs(i) < firstu) && - (Intervs(i+1) > firstu)) - { - Intervs(i) = firstu; - } - if ((Intervs(i) < lastu) && - (Intervs(i+1) > lastu)) - { - Intervs(i + 1) = lastu; - } - - Standard_Real delta = (Intervs(i+1) - Intervs(i))/NbPoints; - for (j = 1; j <= NbPoints && IsLine; ++j) - { - param = Intervs(i) + j*delta; - D0 (C, param, MiddlePoint); - V2 = (MiddlePoint.XYZ() - CurrentPoint.XYZ()); - L2 = V2.Modulus (); - if (L2 > LTol) - { - const Standard_Real aAngle = V2.CrossMagnitude(V1)/(L1*L2); - IsLine = (aAngle < ATol); - } - } - } - - if (IsLine) - { - parameters.Clear(); - points .Clear(); - - PerformLinear(C); - return; - } - else - { - //c'etait une singularite on continue : - //Du = Dusave; - EvaluateDu (C, param, MiddlePoint, Du, NotDone); - } - } - else - { - Du = (lastu-firstu)/2.1; - MiddleU = firstu + Du; - D0 (C, MiddleU, MiddlePoint); - V1 = (MiddlePoint.XYZ() - CurrentPoint.XYZ()); - L1 = V1.Modulus (); - if (L1 < LTol) - { - // L1 < LTol C'est une courbe de longueur nulle, calcul termine : - // on renvoi un segment de 2 points (protection) - parameters.Append (lastu); - points .Append (LastPoint); - return; - } - } - } - - if (Du > Dusave) Du = Dusave; - else Dusave = Du; - - if (Du < uTol) { - Du = lastu - firstu; - if (Du < uTol) { - parameters.Append (lastu); - points .Append (LastPoint); - return; - } - } - - //Traitement normal pour une courbe - Standard_Boolean MorePoints = Standard_True; - Standard_Real U2 = firstu; - Standard_Real AngleMax = angularDeflection * 0.5; //car on prend le point milieu - Standard_Integer aIdx[2] = {Intervs.Lower(), Intervs.Lower()}; // Indexes of intervals of U1 and U2, used to handle non-uniform case. - Standard_Boolean isNeedToCheck = Standard_False; - gp_Pnt aPrevPoint = points.Last(); - - while (MorePoints) { - aIdx[0] = getIntervalIdx(U1, Intervs, aIdx[0]); - U2 += Du; - - if (U2 >= lastu) { //Bout de courbe - U2 = lastu; - CurrentPoint = LastPoint; - Du = U2-U1; - Dusave = Du; - } - else D0 (C, U2, CurrentPoint); //Point suivant - - Standard_Real Coef = 0.0, ACoef = 0., FCoef = 0.; - Standard_Boolean Correction, TooLarge, TooSmall; - TooLarge = Standard_False; - Correction = Standard_True; - TooSmall = Standard_False; - - while (Correction) { //Ajustement Du - if (isNeedToCheck) - { - aIdx[1] = getIntervalIdx(U2, Intervs, aIdx[0]); - if (aIdx[1] > aIdx[0]) // Jump to another polynom. - { - if (Du > (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3) // Set Du to the smallest value and check deflection on it. - { - Du = (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3; - U2 = U1 + Du; - if (U2 > lastu) - U2 = lastu; - D0 (C, U2, CurrentPoint); - } - } - } - MiddleU = (U1+U2)*0.5; //Verif / au point milieu - D0 (C, MiddleU, MiddlePoint); - - V1 = (CurrentPoint.XYZ() - aPrevPoint.XYZ()); //Critere de fleche - V2 = (MiddlePoint.XYZ() - aPrevPoint.XYZ()); - L1 = V1.Modulus (); - - FCoef = (L1 > myMinLen) ? - V1.CrossMagnitude(V2)/(L1*curvatureDeflection) : 0.0; - - V1 = (CurrentPoint.XYZ() - MiddlePoint.XYZ()); //Critere d'angle - L1 = V1.Modulus (); - L2 = V2.Modulus (); - if (L1 > myMinLen && L2 > myMinLen) - { - Standard_Real angg = V1.CrossMagnitude(V2) / (L1 * L2); - ACoef = angg / AngleMax; - } - else - ACoef = 0.0; - - //On retient le plus penalisant - Coef = Max(ACoef, FCoef); - - if (isNeedToCheck && Coef < 0.55) - { - isNeedToCheck = Standard_False; - Du = Dusave; - U2 = U1 + Du; - if (U2 > lastu) - U2 = lastu; - D0 (C, U2, CurrentPoint); - continue; - } - - if (Coef <= 1.0) { - if (Abs (lastu-U2) < uTol) { - parameters.Append (lastu); - points .Append (LastPoint); - MorePoints = Standard_False; - Correction = Standard_False; - } - else { - if (Coef >= 0.55 || TooLarge) { - parameters.Append (U2); - points .Append (CurrentPoint); - aPrevPoint = CurrentPoint; - Correction = Standard_False; - isNeedToCheck = Standard_True; - } - else if (TooSmall) { - Correction = Standard_False; - aPrevPoint = CurrentPoint; - } - else { - TooSmall = Standard_True; - //Standard_Real UUU2 = U2; - Du += Min((U2-U1)*(1.-Coef), Du*Us3); - - U2 = U1 + Du; - if (U2 > lastu) - U2 = lastu; - D0 (C, U2, CurrentPoint); - } - } - } - else { - - if (Coef >= 1.5) { - if (!aPrevPoint.IsEqual(points.Last(), Precision::Confusion())) - { - parameters.Append (U1); - points .Append (aPrevPoint); - } - U2 = MiddleU; - Du = U2-U1; - CurrentPoint = MiddlePoint; - } - else { - Du*=0.9; - U2 = U1 + Du; - D0 (C, U2, CurrentPoint); - TooLarge = Standard_True; - } - - } - } - - Du = U2-U1; - - if (MorePoints) { - if (U1 > firstu) { - if (FCoef > ACoef) { - //La fleche est critere de decoupage - EvaluateDu (C, U2, CurrentPoint, Du, NotDone); - if (NotDone) { - Du += (Du-Dusave)*(Du/Dusave); - if (Du > 1.5 * Dusave) Du = 1.5 * Dusave; - if (Du < 0.75* Dusave) Du = 0.75 * Dusave; - } - } - else { - //L'angle est le critere de decoupage - Du += (Du-Dusave)*(Du/Dusave); - if (Du > 1.5 * Dusave) Du = 1.5 * Dusave; - if (Du < 0.75* Dusave) Du = 0.75 * Dusave; - } - } - - if (Du < uTol) { - Du = lastu - U2; - if (Du < uTol) { - parameters.Append (lastu); - points .Append (LastPoint); - MorePoints = Standard_False; - } - else if (Du*Us3 > uTol) Du*=Us3; - } - U1 = U2; - Dusave = Du; - } - } - //Recalage avant dernier point : - i = points.Length()-1; -// Real d = points (i).Distance (points (i+1)); -// if (Abs(parameters (i) - parameters (i+1))<= 0.000001 || d < Precision::Confusion()) { -// cout<<"deux points confondus"<= 2) { - MiddleU = parameters (i-1); - MiddleU = (lastu + MiddleU)*0.5; - D0 (C, MiddleU, MiddlePoint); - parameters.SetValue (i, MiddleU); - points .SetValue (i, MiddlePoint); - } - - //-- On rajoute des points aux milieux des segments si le nombre - //-- mini de points n'est pas atteint - //-- - Standard_Integer Nbp = points.Length(); - - //std::cout << "GCPnts_TangentialDeflection: Number of Points (" << Nbp << " " << minNbPnts << " )" << std::endl; - - while(Nbp < minNbPnts) - { - for (i = 2; i <= Nbp; i += 2) - { - MiddleU = (parameters.Value(i-1)+parameters.Value(i))*0.5; - D0 (C, MiddleU, MiddlePoint); - parameters.InsertBefore(i,MiddleU); - points.InsertBefore(i,MiddlePoint); - Nbp++; - } - } - //Additional check for intervals - Standard_Real MinLen2 = myMinLen * myMinLen; - Standard_Integer MaxNbp = 10 * Nbp; - for(i = 1; i < Nbp; ++i) - { - U1 = parameters(i); - U2 = parameters(i + 1); - - if (U2 - U1 <= uTol) - { - continue; - } - - // Check maximal deflection on interval; - Standard_Real dmax = 0.; - Standard_Real umax = 0.; - Standard_Real amax = 0.; - EstimDefl(C, U1, U2, dmax, umax); - const gp_Pnt& P1 = points(i); - const gp_Pnt& P2 = points(i+1); - D0(C, umax, MiddlePoint); - amax = EstimAngl(P1, MiddlePoint, P2); - if(dmax > curvatureDeflection || amax > AngleMax) - { - if(umax - U1 > uTol && U2 - umax > uTol) - { - if (P1.SquareDistance(MiddlePoint) > MinLen2 && P2.SquareDistance(MiddlePoint) > MinLen2) - { - parameters.InsertAfter(i, umax); - points.InsertAfter(i, MiddlePoint); - ++Nbp; - --i; //To compensate ++i in loop header: i must point to first part of split interval - if(Nbp > MaxNbp) - { - break; - } - } - } - } - } - // -} - -//======================================================================= -//function : EstimDefl -//purpose : Estimation of maximal deflection for interval [U1, U2] -// -//======================================================================= -void GCPnts_TangentialDeflection::EstimDefl (const TheCurve& C, - const Standard_Real U1, const Standard_Real U2, - Standard_Real& MaxDefl, Standard_Real& UMax) -{ - Standard_Real Du = (lastu - firstu); - // - TheMaxCurvLinDist aFunc(C, U1, U2); - // - const Standard_Integer aNbIter = 100; - Standard_Real reltol = Max(1.e-3, 2.*uTol/((Abs(U1) + Abs(U2)))); - // - math_BrentMinimum anOptLoc(reltol, aNbIter, uTol); - anOptLoc.Perform(aFunc, U1, (U1+U2)/2., U2); - if(anOptLoc.IsDone()) - { - MaxDefl = Sqrt(-anOptLoc.Minimum()); - UMax = anOptLoc.Location(); - return; - } - // - math_Vector aLowBorder(1,1); - math_Vector aUppBorder(1,1); - math_Vector aSteps(1,1); - // - aSteps(1) = Max(0.1 * Du, 100. * uTol); - Standard_Integer aNbParticles = Max(8, RealToInt(32 * (U2 - U1) / Du)); - // - aLowBorder(1) = U1; - aUppBorder(1) = U2; - // - // - Standard_Real aValue; - math_Vector aT(1,1); - TheMaxCurvLinDistMV aFuncMV(aFunc); - - math_PSO aFinder(&aFuncMV, aLowBorder, aUppBorder, aSteps, aNbParticles); - aFinder.Perform(aSteps, aValue, aT); - // - anOptLoc.Perform(aFunc, Max(aT(1) - aSteps(1), U1) , aT(1), Min(aT(1) + aSteps(1),U2)); - if(anOptLoc.IsDone()) - { - MaxDefl = Sqrt(-anOptLoc.Minimum()); - UMax = anOptLoc.Location(); - return; - } - MaxDefl = Sqrt(-aValue); - UMax = aT(1); - // -}