From 1f37f1d50a91cda73f8a4445525abb3e60a31f84 Mon Sep 17 00:00:00 2001 From: ifv Date: Thu, 3 Mar 2022 15:05:23 +0300 Subject: [PATCH] 0031661: Modeling Data - Exception when projecting parabola or hyperbola to plane ProjLib/ProjLib_ProjectOnPlane.cxx - formatting 0031661: Modeling Data - Algorithm crashes when projecting parabola or hyperbola to plane ProjLib/ProjLib_ProjectOnPlane.cxx - building of analytical parabola and hyperbola is added bugs/moddata_3/bug31661_* - new test cases are added --- src/ProjLib/ProjLib_ProjectOnPlane.cxx | 1645 ++++++++++++++---------- src/ProjLib/ProjLib_ProjectOnPlane.hxx | 5 +- tests/bugs/moddata_3/bug31661_1 | 26 + tests/bugs/moddata_3/bug31661_2 | 26 + 4 files changed, 1020 insertions(+), 682 deletions(-) create mode 100644 tests/bugs/moddata_3/bug31661_1 create mode 100644 tests/bugs/moddata_3/bug31661_2 diff --git a/src/ProjLib/ProjLib_ProjectOnPlane.cxx b/src/ProjLib/ProjLib_ProjectOnPlane.cxx index 07d3292333..4330a8666c 100644 --- a/src/ProjLib/ProjLib_ProjectOnPlane.cxx +++ b/src/ProjLib/ProjLib_ProjectOnPlane.cxx @@ -47,8 +47,14 @@ #include #include #include +#include +#include +#include +#include +#include - +const Standard_Real aParabolaLimit = 20000.; +const Standard_Real aHyperbolaLimit = 10.; //======================================================================= //function : OnPlane_Value @@ -56,18 +62,18 @@ //======================================================================= static gp_Pnt OnPlane_Value(const Standard_Real U, - const Handle(Adaptor3d_Curve)& aCurvePtr, - const gp_Ax3& Pl, - const gp_Dir& D) + const Handle(Adaptor3d_Curve)& aCurvePtr, + const gp_Ax3& Pl, + const gp_Dir& D) { // PO . Z / Z = Pl.Direction() // Proj(u) = P(u) + ------- * D avec \ O = Pl.Location() // D . Z gp_Pnt Point = aCurvePtr->Value(U); - - gp_Vec PO(Point,Pl.Location()); - + + gp_Vec PO(Point, Pl.Location()); + Standard_Real Alpha = PO * gp_Vec(Pl.Direction()); Alpha /= D * Pl.Direction(); Point.SetXYZ(Point.XYZ() + Alpha * D.XYZ()); @@ -81,25 +87,25 @@ static gp_Pnt OnPlane_Value(const Standard_Real U, //======================================================================= static gp_Vec OnPlane_DN(const Standard_Real U, - const Standard_Integer DerivativeRequest, - const Handle(Adaptor3d_Curve)& aCurvePtr, - const gp_Ax3& Pl, - const gp_Dir& D) + const Standard_Integer DerivativeRequest, + const Handle(Adaptor3d_Curve)& aCurvePtr, + const gp_Ax3& Pl, + const gp_Dir& D) { // PO . Z / Z = Pl.Direction() // Proj(u) = P(u) + ------- * D avec \ O = Pl.Location() // D . Z - gp_Vec Vector = aCurvePtr->DN(U,DerivativeRequest); + gp_Vec Vector = aCurvePtr->DN(U, DerivativeRequest); gp_Dir Z = Pl.Direction(); - - Standard_Real - Alpha = Vector * gp_Vec(Z); - Alpha /= D * Z; - Vector.SetXYZ( Vector.XYZ() - Alpha * D.XYZ()); - return Vector ; + Standard_Real + Alpha = Vector * gp_Vec(Z); + Alpha /= D * Z; + + Vector.SetXYZ(Vector.XYZ() - Alpha * D.XYZ()); + return Vector; } //======================================================================= @@ -108,27 +114,27 @@ static gp_Vec OnPlane_DN(const Standard_Real U, //======================================================================= static Standard_Boolean OnPlane_D1(const Standard_Real U, - gp_Pnt& P, - gp_Vec& V, - const Handle(Adaptor3d_Curve)& aCurvePtr, - const gp_Ax3& Pl, - const gp_Dir& D) + gp_Pnt& P, + gp_Vec& V, + const Handle(Adaptor3d_Curve)& aCurvePtr, + const gp_Ax3& Pl, + const gp_Dir& D) { Standard_Real Alpha; gp_Pnt Point; gp_Vec Vector; - + gp_Dir Z = Pl.Direction(); aCurvePtr->D1(U, Point, Vector); // evaluate the point as in `OnPlane_Value` - gp_Vec PO(Point,Pl.Location()); - Alpha = PO * gp_Vec(Z); + gp_Vec PO(Point, Pl.Location()); + Alpha = PO * gp_Vec(Z); Alpha /= D * Z; P.SetXYZ(Point.XYZ() + Alpha * D.XYZ()); - + // evaluate the derivative. // // d(Proj) d(P) 1 d(P) @@ -136,10 +142,10 @@ static Standard_Boolean OnPlane_D1(const Standard_Real U, // dU dU ( D . Z) dU // - Alpha = Vector * gp_Vec(Z); - Alpha /= D * Z; + Alpha = Vector * gp_Vec(Z); + Alpha /= D * Z; - V.SetXYZ( Vector.XYZ() - Alpha * D.XYZ()); + V.SetXYZ(Vector.XYZ() - Alpha * D.XYZ()); return Standard_True; } @@ -149,25 +155,25 @@ static Standard_Boolean OnPlane_D1(const Standard_Real U, //======================================================================= static Standard_Boolean OnPlane_D2(const Standard_Real U, - gp_Pnt& P, - gp_Vec& V1, - gp_Vec& V2, - const Handle(Adaptor3d_Curve) & aCurvePtr, - const gp_Ax3& Pl, - const gp_Dir& D) + gp_Pnt& P, + gp_Vec& V1, + gp_Vec& V2, + const Handle(Adaptor3d_Curve) & aCurvePtr, + const gp_Ax3& Pl, + const gp_Dir& D) { Standard_Real Alpha; gp_Pnt Point; gp_Vec Vector1, - Vector2; - + Vector2; + gp_Dir Z = Pl.Direction(); aCurvePtr->D2(U, Point, Vector1, Vector2); // evaluate the point as in `OnPlane_Value` - gp_Vec PO(Point,Pl.Location()); - Alpha = PO * gp_Vec(Z); + gp_Vec PO(Point, Pl.Location()); + Alpha = PO * gp_Vec(Z); Alpha /= D * Z; P.SetXYZ(Point.XYZ() + Alpha * D.XYZ()); @@ -178,15 +184,15 @@ static Standard_Boolean OnPlane_D2(const Standard_Real U, // dU dU ( D . Z) dU // - Alpha = Vector1 * gp_Vec(Z); - Alpha /= D * Z; + Alpha = Vector1 * gp_Vec(Z); + Alpha /= D * Z; - V1.SetXYZ( Vector1.XYZ() - Alpha * D.XYZ()); + V1.SetXYZ(Vector1.XYZ() - Alpha * D.XYZ()); - Alpha = Vector2 * gp_Vec(Z); - Alpha /= D * Z; + Alpha = Vector2 * gp_Vec(Z); + Alpha /= D * Z; - V2.SetXYZ( Vector2.XYZ() - Alpha * D.XYZ()); + V2.SetXYZ(Vector2.XYZ() - Alpha * D.XYZ()); return Standard_True; } @@ -196,30 +202,30 @@ static Standard_Boolean OnPlane_D2(const Standard_Real U, //======================================================================= static Standard_Boolean OnPlane_D3(const Standard_Real U, - gp_Pnt& P, - gp_Vec& V1, - gp_Vec& V2, - gp_Vec& V3, - const Handle(Adaptor3d_Curve)& aCurvePtr, - const gp_Ax3& Pl, - const gp_Dir& D) + gp_Pnt& P, + gp_Vec& V1, + gp_Vec& V2, + gp_Vec& V3, + const Handle(Adaptor3d_Curve)& aCurvePtr, + const gp_Ax3& Pl, + const gp_Dir& D) { Standard_Real Alpha; gp_Pnt Point; gp_Vec Vector1, - Vector2, - Vector3; - + Vector2, + Vector3; + gp_Dir Z = Pl.Direction(); aCurvePtr->D3(U, Point, Vector1, Vector2, Vector3); // evaluate the point as in `OnPlane_Value` - gp_Vec PO(Point,Pl.Location()); - Alpha = PO * gp_Vec(Z); + gp_Vec PO(Point, Pl.Location()); + Alpha = PO * gp_Vec(Z); Alpha /= D * Z; P.SetXYZ(Point.XYZ() + Alpha * D.XYZ()); - + // evaluate the derivative. // // d(Proj) d(P) 1 d(P) @@ -227,19 +233,19 @@ static Standard_Boolean OnPlane_D3(const Standard_Real U, // dU dU ( D . Z) dU // - Alpha = Vector1 * gp_Vec(Z); - Alpha /= D * Z; + Alpha = Vector1 * gp_Vec(Z); + Alpha /= D * Z; - V1.SetXYZ( Vector1.XYZ() - Alpha * D.XYZ()); + V1.SetXYZ(Vector1.XYZ() - Alpha * D.XYZ()); - Alpha = Vector2 * gp_Vec(Z); - Alpha /= D * Z; + Alpha = Vector2 * gp_Vec(Z); + Alpha /= D * Z; - V2.SetXYZ( Vector2.XYZ() - Alpha * D.XYZ()); - Alpha = Vector3 * gp_Vec(Z); - Alpha /= D * Z; + V2.SetXYZ(Vector2.XYZ() - Alpha * D.XYZ()); + Alpha = Vector3 * gp_Vec(Z); + Alpha /= D * Z; - V3.SetXYZ( Vector3.XYZ() - Alpha * D.XYZ()); + V3.SetXYZ(Vector3.XYZ() - Alpha * D.XYZ()); return Standard_True; } @@ -255,43 +261,76 @@ class ProjLib_OnPlane : public AppCont_Function gp_Ax3 myPlane; gp_Dir myDirection; -public : +public: - ProjLib_OnPlane(const Handle(Adaptor3d_Curve)& C, - const gp_Ax3& Pl, - const gp_Dir& D) -: myCurve(C), - myPlane(Pl), - myDirection(D) + ProjLib_OnPlane(const Handle(Adaptor3d_Curve)& C, + const gp_Ax3& Pl, + const gp_Dir& D) + : myCurve(C), + myPlane(Pl), + myDirection(D) { myNbPnt = 1; myNbPnt2d = 0; } Standard_Real FirstParameter() const - {return myCurve->FirstParameter();} - - Standard_Real LastParameter() const - {return myCurve->LastParameter();} + { + return myCurve->FirstParameter(); + } + + Standard_Real LastParameter() const + { + return myCurve->LastParameter(); + } + + Standard_Boolean Value(const Standard_Real theT, + NCollection_Array1& /*thePnt2d*/, + NCollection_Array1& thePnt) const + { + thePnt(1) = OnPlane_Value(theT, myCurve, myPlane, myDirection); + return Standard_True; + } - Standard_Boolean Value(const Standard_Real theT, - NCollection_Array1& /*thePnt2d*/, - NCollection_Array1& thePnt) const - { - thePnt(1) = OnPlane_Value(theT, myCurve, myPlane, myDirection); - return Standard_True; - } - Standard_Boolean D1(const Standard_Real theT, - NCollection_Array1& /*theVec2d*/, - NCollection_Array1& theVec) const + NCollection_Array1& /*theVec2d*/, + NCollection_Array1& theVec) const { gp_Pnt aDummyPnt; - return OnPlane_D1(theT, aDummyPnt, theVec(1),myCurve,myPlane,myDirection); + return OnPlane_D1(theT, aDummyPnt, theVec(1), myCurve, myPlane, myDirection); } }; +//======================================================================= +// class : ProjLib_MaxCurvature +//purpose : Use to search apex of parabola or hyperbola, which is its projection +// on a plane. Apex is point with maximal curvature +//======================================================================= + +class ProjLib_MaxCurvature : public math_Function + +{ + +public: + + ProjLib_MaxCurvature(LProp3d_CLProps& theProps): + myProps(&theProps) + { + } + + virtual Standard_Boolean Value(const Standard_Real X, Standard_Real& F) + { + myProps->SetParameter(X); + F = -myProps->Curvature(); + return Standard_True; + } + +private: + + LProp3d_CLProps* myProps; + +}; //=====================================================================// @@ -309,20 +348,30 @@ public : //purpose : //======================================================================= -static void PerformApprox (const Handle(Adaptor3d_Curve)& C, - const gp_Ax3& Pl, - const gp_Dir& D, - Handle(Geom_BSplineCurve) &BSplineCurvePtr) +static void PerformApprox(const Handle(Adaptor3d_Curve)& C, + const gp_Ax3& Pl, + const gp_Dir& D, + Handle(Geom_BSplineCurve) &BSplineCurvePtr) { - ProjLib_OnPlane F(C,Pl,D); + ProjLib_OnPlane F(C, Pl, D); Standard_Integer Deg1, Deg2; Deg1 = 8; Deg2 = 8; - - Approx_FitAndDivide Fit(Deg1,Deg2,Precision::Approximation(), - Precision::PApproximation(),Standard_True); - Fit.SetMaxSegments(100); + if (C->GetType() == GeomAbs_Parabola) + { + Deg1 = 2; Deg2 = 2; + } + Standard_Integer aNbSegm = 100; + if (C->GetType() == GeomAbs_Hyperbola) + { + Deg1 = 14; + Deg2 = 14; + aNbSegm = 1000; + } + Approx_FitAndDivide Fit(Deg1, Deg2, Precision::Approximation(), + Precision::PApproximation(), Standard_True); + Fit.SetMaxSegments(aNbSegm); Fit.Perform(F); if (!Fit.IsAllApproximated()) { @@ -331,64 +380,77 @@ static void PerformApprox (const Handle(Adaptor3d_Curve)& C, Standard_Integer i; Standard_Integer NbCurves = Fit.NbMultiCurves(); Standard_Integer MaxDeg = 0; - + // Pour transformer la MultiCurve en BSpline, il faut que toutes // les Bezier la constituant aient le meme degre -> Calcul de MaxDeg - Standard_Integer NbPoles = 1; + Standard_Integer NbPoles = 1; for (i = 1; i <= NbCurves; i++) { Standard_Integer Deg = Fit.Value(i).Degree(); - MaxDeg = Max ( MaxDeg, Deg); + MaxDeg = Max(MaxDeg, Deg); } NbPoles = MaxDeg * NbCurves + 1; //Poles sur la BSpline - TColgp_Array1OfPnt Poles( 1, NbPoles); - - TColgp_Array1OfPnt TempPoles( 1, MaxDeg + 1); //pour augmentation du degre - - TColStd_Array1OfReal Knots( 1, NbCurves + 1); //Noeuds de la BSpline - + TColgp_Array1OfPnt Poles(1, NbPoles); + + TColgp_Array1OfPnt TempPoles(1, MaxDeg + 1); //pour augmentation du degre + + TColStd_Array1OfReal Knots(1, NbCurves + 1); //Noeuds de la BSpline + Standard_Integer Compt = 1; + Standard_Real anErrMax = 0., anErr3d, anErr2d; for (i = 1; i <= Fit.NbMultiCurves(); i++) { - Fit.Parameters(i, Knots(i), Knots(i+1)); - - AppParCurves_MultiCurve MC = Fit.Value( i); //Charge la Ieme Curve - TColgp_Array1OfPnt LocalPoles( 1, MC.Degree() + 1);//Recupere les poles + Fit.Parameters(i, Knots(i), Knots(i + 1)); + Fit.Error(i, anErr3d, anErr2d); + anErrMax = Max(anErrMax, anErr3d); + AppParCurves_MultiCurve MC = Fit.Value(i); //Charge la Ieme Curve + TColgp_Array1OfPnt LocalPoles(1, MC.Degree() + 1);//Recupere les poles MC.Curve(1, LocalPoles); - + //Augmentation eventuelle du degre - if (MaxDeg > MC.Degree() ) { + if (MaxDeg > MC.Degree()) { BSplCLib::IncreaseDegree(MaxDeg, LocalPoles, BSplCLib::NoWeights(), - TempPoles, BSplCLib::NoWeights()); + TempPoles, BSplCLib::NoWeights()); //mise a jour des poles de la PCurve - for (Standard_Integer j = 1 ; j <= MaxDeg + 1; j++) { - Poles.SetValue( Compt, TempPoles( j)); - Compt++; + for (Standard_Integer j = 1; j <= MaxDeg + 1; j++) { + Poles.SetValue(Compt, TempPoles(j)); + Compt++; } } else { //mise a jour des poles de la PCurve - for (Standard_Integer j = 1 ; j <= MaxDeg + 1; j++) { - Poles.SetValue( Compt, LocalPoles( j)); - Compt++; + for (Standard_Integer j = 1; j <= MaxDeg + 1; j++) { + Poles.SetValue(Compt, LocalPoles(j)); + Compt++; } - } - + } + Compt--; } - + //mise a jour des fields de ProjLib_Approx - Standard_Integer - NbKnots = NbCurves + 1; + Standard_Integer + NbKnots = NbCurves + 1; TColStd_Array1OfInteger Mults(1, NbKnots); - Mults.SetValue( 1, MaxDeg + 1); - for ( i = 2; i <= NbCurves; i++) { - Mults.SetValue( i, MaxDeg); + Mults.SetValue(1, MaxDeg + 1); + for (i = 2; i <= NbCurves; i++) { + Mults.SetValue(i, MaxDeg); } - Mults.SetValue( NbKnots, MaxDeg + 1); - BSplineCurvePtr = - new Geom_BSplineCurve(Poles,Knots,Mults,MaxDeg,Standard_False); + Mults.SetValue(NbKnots, MaxDeg + 1); + BSplineCurvePtr = + new Geom_BSplineCurve(Poles, Knots, Mults, MaxDeg, Standard_False); + + //Try to smooth + Standard_Integer m1 = MaxDeg - 1; + for (i = 2; i < NbKnots; ++i) + { + if (BSplineCurvePtr->Multiplicity(i) == MaxDeg) + { + BSplineCurvePtr->RemoveKnot(i, m1, anErrMax); + } + } + } @@ -398,12 +460,12 @@ static void PerformApprox (const Handle(Adaptor3d_Curve)& C, //======================================================================= ProjLib_ProjectOnPlane::ProjLib_ProjectOnPlane() : -myKeepParam(Standard_False), -myFirstPar(0.), -myLastPar(0.), -myTolerance(0.), -myType (GeomAbs_OtherCurve), -myIsApprox (Standard_False) + myKeepParam(Standard_False), + myFirstPar(0.), + myLastPar(0.), + myTolerance(0.), + myType(GeomAbs_OtherCurve), + myIsApprox(Standard_False) { } @@ -412,15 +474,15 @@ myIsApprox (Standard_False) //purpose : //======================================================================= -ProjLib_ProjectOnPlane::ProjLib_ProjectOnPlane(const gp_Ax3& Pl) : -myPlane (Pl) , -myDirection (Pl.Direction()) , -myKeepParam(Standard_False), -myFirstPar(0.), -myLastPar(0.), -myTolerance(0.), -myType (GeomAbs_OtherCurve), -myIsApprox (Standard_False) +ProjLib_ProjectOnPlane::ProjLib_ProjectOnPlane(const gp_Ax3& Pl) : + myPlane(Pl), + myDirection(Pl.Direction()), + myKeepParam(Standard_False), + myFirstPar(0.), + myLastPar(0.), + myTolerance(0.), + myType(GeomAbs_OtherCurve), + myIsApprox(Standard_False) { } @@ -430,20 +492,20 @@ myIsApprox (Standard_False) //======================================================================= ProjLib_ProjectOnPlane::ProjLib_ProjectOnPlane(const gp_Ax3& Pl, - const gp_Dir& D ) : -myPlane (Pl) , -myDirection (D) , -myKeepParam(Standard_False), -myFirstPar(0.), -myLastPar(0.), -myTolerance(0.), -myType (GeomAbs_OtherCurve), -myIsApprox (Standard_False) + const gp_Dir& D) : + myPlane(Pl), + myDirection(D), + myKeepParam(Standard_False), + myFirstPar(0.), + myLastPar(0.), + myTolerance(0.), + myType(GeomAbs_OtherCurve), + myIsApprox(Standard_False) { -// if ( Abs(D * Pl.Direction()) < Precision::Confusion()) { -// throw Standard_ConstructionError -// ("ProjLib_ProjectOnPlane: The Direction and the Plane are parallel"); -// } + // if ( Abs(D * Pl.Direction()) < Precision::Confusion()) { + // throw Standard_ConstructionError + // ("ProjLib_ProjectOnPlane: The Direction and the Plane are parallel"); + // } } //======================================================================= @@ -482,14 +544,14 @@ Handle(Adaptor3d_Curve) ProjLib_ProjectOnPlane::ShallowCopy() const //======================================================================= static gp_Pnt ProjectPnt(const gp_Ax3& ThePlane, - const gp_Dir& TheDir, - const gp_Pnt& Point) + const gp_Dir& TheDir, + const gp_Pnt& Point) { - gp_Vec PO(Point,ThePlane.Location()); + gp_Vec PO(Point, ThePlane.Location()); Standard_Real Alpha = PO * gp_Vec(ThePlane.Direction()); Alpha /= TheDir * ThePlane.Direction(); - + gp_Pnt P; P.SetXYZ(Point.XYZ() + Alpha * TheDir.XYZ()); @@ -504,13 +566,13 @@ static gp_Pnt ProjectPnt(const gp_Ax3& ThePlane, //======================================================================= static gp_Vec ProjectVec(const gp_Ax3& ThePlane, - const gp_Dir& TheDir, - const gp_Vec& Vec) + const gp_Dir& TheDir, + const gp_Vec& Vec) { gp_Vec D = Vec; gp_Vec Z = ThePlane.Direction(); - D -= ( (Vec * Z) / (TheDir * Z)) * TheDir; + D -= ((Vec * Z) / (TheDir * Z)) * TheDir; return D; } @@ -521,461 +583,463 @@ static gp_Vec ProjectVec(const gp_Ax3& ThePlane, //======================================================================= void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_Curve)& C, - const Standard_Real Tolerance, - const Standard_Boolean KeepParametrization) - + const Standard_Real Tolerance, + const Standard_Boolean KeepParametrization) + { - myCurve = C; - myType = GeomAbs_OtherCurve; - myIsApprox = Standard_False; - myTolerance = Tolerance ; + myCurve = C; + myType = GeomAbs_OtherCurve; + myIsApprox = Standard_False; + myTolerance = Tolerance; Handle(Geom_BSplineCurve) ApproxCurve; Handle(GeomAdaptor_Curve) aGAHCurve; Handle(Geom_Line) GeomLinePtr; - Handle(Geom_Circle) GeomCirclePtr ; - Handle(Geom_Ellipse) GeomEllipsePtr ; - Handle(Geom_Hyperbola) GeomHyperbolaPtr ; + Handle(Geom_Circle) GeomCirclePtr; + Handle(Geom_Ellipse) GeomEllipsePtr; + Handle(Geom_Hyperbola) GeomHyperbolaPtr; + Handle(Geom_Parabola) GeomParabolaPtr; - gp_Lin aLine; + gp_Lin aLine; gp_Elips Elips; -// gp_Hypr Hypr ; + // gp_Hypr Hypr ; - Standard_Integer num_knots ; + Standard_Integer num_knots; GeomAbs_CurveType Type = C->GetType(); gp_Ax2 Axis; - Standard_Real R1 =0., R2 =0.; + Standard_Real R1 = 0., R2 = 0.; myKeepParam = KeepParametrization; - switch ( Type) { - case GeomAbs_Line: - { - // P(u) = O + u * Xc - // ==> Q(u) = f(P(u)) - // = f(O) + u * f(Xc) + switch (Type) { + case GeomAbs_Line: + { + // P(u) = O + u * Xc + // ==> Q(u) = f(P(u)) + // = f(O) + u * f(Xc) - gp_Lin L = myCurve->Line(); - gp_Vec Xc = ProjectVec(myPlane,myDirection,gp_Vec(L.Direction())); + gp_Lin L = myCurve->Line(); + gp_Vec Xc = ProjectVec(myPlane, myDirection, gp_Vec(L.Direction())); - if ( Xc.Magnitude() < Precision::Confusion()) { // line orthog au plan - myType = GeomAbs_BSplineCurve; - gp_Pnt P = ProjectPnt(myPlane,myDirection,L.Location()); - TColStd_Array1OfInteger Mults(1,2); Mults.Init(2); - TColgp_Array1OfPnt Poles(1,2); Poles.Init(P); - TColStd_Array1OfReal Knots(1,2); - Knots(1) = myCurve->FirstParameter(); - Knots(2) = myCurve->LastParameter(); - Handle(Geom_BSplineCurve) BSP = - new Geom_BSplineCurve(Poles,Knots,Mults,1); + if (Xc.Magnitude() < Precision::Confusion()) { // line orthog au plan + myType = GeomAbs_BSplineCurve; + gp_Pnt P = ProjectPnt(myPlane, myDirection, L.Location()); + TColStd_Array1OfInteger Mults(1, 2); Mults.Init(2); + TColgp_Array1OfPnt Poles(1, 2); Poles.Init(P); + TColStd_Array1OfReal Knots(1, 2); + Knots(1) = myCurve->FirstParameter(); + Knots(2) = myCurve->LastParameter(); + Handle(Geom_BSplineCurve) BSP = + new Geom_BSplineCurve(Poles, Knots, Mults, 1); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(BSP); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - } - else if ( Abs( Xc.Magnitude() - 1.) < Precision::Confusion()) { - myType = GeomAbs_Line; - gp_Pnt P = ProjectPnt(myPlane,myDirection,L.Location()); - myFirstPar = myCurve->FirstParameter(); - myLastPar = myCurve->LastParameter(); - aLine = gp_Lin(P,gp_Dir(Xc)); - GeomLinePtr = new Geom_Line(aLine); - -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(GeomLinePtr, - myCurve->FirstParameter(), - myCurve->LastParameter() ); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin + GeomAdaptor_Curve aGACurve(BSP); + myResult = new GeomAdaptor_Curve(aGACurve); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End + } + else if (Abs(Xc.Magnitude() - 1.) < Precision::Confusion()) { + myType = GeomAbs_Line; + gp_Pnt P = ProjectPnt(myPlane, myDirection, L.Location()); + myFirstPar = myCurve->FirstParameter(); + myLastPar = myCurve->LastParameter(); + aLine = gp_Lin(P, gp_Dir(Xc)); + GeomLinePtr = new Geom_Line(aLine); + + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin + GeomAdaptor_Curve aGACurve(GeomLinePtr, + myCurve->FirstParameter(), + myCurve->LastParameter()); + myResult = new GeomAdaptor_Curve(aGACurve); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End + } + else { + myType = GeomAbs_Line; + gp_Pnt P = ProjectPnt(myPlane, myDirection, L.Location()); + aLine = gp_Lin(P, gp_Dir(Xc)); + Standard_Real Udeb, Ufin; + + // eval the first and last parameters of the projected curve + Udeb = myCurve->FirstParameter(); + Ufin = myCurve->LastParameter(); + gp_Pnt P1 = ProjectPnt(myPlane, myDirection, + myCurve->Value(Udeb)); + gp_Pnt P2 = ProjectPnt(myPlane, myDirection, + myCurve->Value(Ufin)); + myFirstPar = gp_Vec(aLine.Direction()).Dot(gp_Vec(P, P1)); + myLastPar = gp_Vec(aLine.Direction()).Dot(gp_Vec(P, P2)); + GeomLinePtr = new Geom_Line(aLine); + if (!myKeepParam) { + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin + GeomAdaptor_Curve aGACurve(GeomLinePtr, + myFirstPar, + myLastPar); + myResult = new GeomAdaptor_Curve(aGACurve); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End } else { - myType = GeomAbs_Line; - gp_Pnt P = ProjectPnt(myPlane,myDirection,L.Location()); - aLine = gp_Lin(P,gp_Dir(Xc)); - Standard_Real Udeb, Ufin; - - // eval the first and last parameters of the projected curve - Udeb = myCurve->FirstParameter(); - Ufin = myCurve->LastParameter(); - gp_Pnt P1 = ProjectPnt(myPlane,myDirection, - myCurve->Value(Udeb)); - gp_Pnt P2 = ProjectPnt(myPlane,myDirection, - myCurve->Value(Ufin)); - myFirstPar = gp_Vec(aLine.Direction()).Dot(gp_Vec(P,P1)); - myLastPar = gp_Vec(aLine.Direction()).Dot(gp_Vec(P,P2)); - GeomLinePtr = new Geom_Line(aLine); - if (!myKeepParam) { -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(GeomLinePtr, - myFirstPar, - myLastPar) ; - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - } - else { - myType = GeomAbs_BSplineCurve; - // - // make a linear BSpline of degree 1 between the end points of - // the projected line - // - Handle(Geom_TrimmedCurve) NewTrimCurvePtr = - new Geom_TrimmedCurve(GeomLinePtr, - myFirstPar, - myLastPar) ; - - Handle(Geom_BSplineCurve) NewCurvePtr = - GeomConvert::CurveToBSplineCurve(NewTrimCurvePtr) ; - num_knots = NewCurvePtr->NbKnots() ; - TColStd_Array1OfReal BsplineKnots(1,num_knots) ; - NewCurvePtr->Knots(BsplineKnots) ; - - BSplCLib::Reparametrize(myCurve->FirstParameter(), - myCurve->LastParameter(), - BsplineKnots) ; - - NewCurvePtr->SetKnots(BsplineKnots) ; -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(NewCurvePtr); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - } + myType = GeomAbs_BSplineCurve; + // + // make a linear BSpline of degree 1 between the end points of + // the projected line + // + Handle(Geom_TrimmedCurve) NewTrimCurvePtr = + new Geom_TrimmedCurve(GeomLinePtr, + myFirstPar, + myLastPar); + + Handle(Geom_BSplineCurve) NewCurvePtr = + GeomConvert::CurveToBSplineCurve(NewTrimCurvePtr); + num_knots = NewCurvePtr->NbKnots(); + TColStd_Array1OfReal BsplineKnots(1, num_knots); + NewCurvePtr->Knots(BsplineKnots); + + BSplCLib::Reparametrize(myCurve->FirstParameter(), + myCurve->LastParameter(), + BsplineKnots); + + NewCurvePtr->SetKnots(BsplineKnots); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin + GeomAdaptor_Curve aGACurve(NewCurvePtr); + myResult = new GeomAdaptor_Curve(aGACurve); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End } - break; } + break; + } case GeomAbs_Circle: - { - // Pour le cercle et l ellipse on a les relations suivantes: - // ( Rem : pour le cercle R1 = R2 = R) - // P(u) = O + R1 * Cos(u) * Xc + R2 * Sin(u) * Yc - // ==> Q(u) = f(P(u)) - // = f(O) + R1 * Cos(u) * f(Xc) + R2 * Sin(u) * f(Yc) - - gp_Circ Circ = myCurve->Circle(); - Axis = Circ.Position(); - R1 = R2 = Circ.Radius(); + { + // Pour le cercle et l ellipse on a les relations suivantes: + // ( Rem : pour le cercle R1 = R2 = R) + // P(u) = O + R1 * Cos(u) * Xc + R2 * Sin(u) * Yc + // ==> Q(u) = f(P(u)) + // = f(O) + R1 * Cos(u) * f(Xc) + R2 * Sin(u) * f(Yc) - } - Standard_FALLTHROUGH + gp_Circ Circ = myCurve->Circle(); + Axis = Circ.Position(); + R1 = R2 = Circ.Radius(); + + } + Standard_FALLTHROUGH case GeomAbs_Ellipse: + { + if (Type == GeomAbs_Ellipse) { + gp_Elips E = myCurve->Ellipse(); + Axis = E.Position(); + R1 = E.MajorRadius(); + R2 = E.MinorRadius(); + } + + // Common Code for CIRCLE & ELLIPSE begin here + gp_Dir X = Axis.XDirection(); + gp_Dir Y = Axis.YDirection(); + gp_Vec VDx = ProjectVec(myPlane, myDirection, X); + gp_Vec VDy = ProjectVec(myPlane, myDirection, Y); + gp_Dir Dx, Dy; + + Standard_Real Tol2 = myTolerance*myTolerance; + if (VDx.SquareMagnitude() < Tol2 || + VDy.SquareMagnitude() < Tol2 || + VDx.CrossSquareMagnitude(VDy) < Tol2) { - if ( Type == GeomAbs_Ellipse) { - gp_Elips E = myCurve->Ellipse(); - Axis = E.Position(); - R1 = E.MajorRadius(); - R2 = E.MinorRadius(); - } + myIsApprox = Standard_True; + } - // Common Code for CIRCLE & ELLIPSE begin here - gp_Dir X = Axis.XDirection(); - gp_Dir Y = Axis.YDirection(); - gp_Vec VDx = ProjectVec(myPlane,myDirection,X); - gp_Vec VDy = ProjectVec(myPlane,myDirection,Y); - gp_Dir Dx,Dy; + if (!myIsApprox) + { + Dx = gp_Dir(VDx); + Dy = gp_Dir(VDy); + gp_Pnt O = Axis.Location(); + gp_Pnt P = ProjectPnt(myPlane, myDirection, O); + gp_Pnt Px = ProjectPnt(myPlane, myDirection, O.Translated(R1*gp_Vec(X))); + gp_Pnt Py = ProjectPnt(myPlane, myDirection, O.Translated(R2*gp_Vec(Y))); + Standard_Real Major = P.Distance(Px); + Standard_Real Minor = P.Distance(Py); - Standard_Real Tol2 = myTolerance*myTolerance; - if (VDx.SquareMagnitude() < Tol2 || - VDy.SquareMagnitude() < Tol2 || - VDx.CrossSquareMagnitude(VDy) < Tol2) + if (myKeepParam) { - myIsApprox = Standard_True; + myIsApprox = !gp_Dir(VDx).IsNormal(gp_Dir(VDy), Precision::Angular()); + } + else + { + // Since it is not necessary to keep the same parameter for the point on the original and on the projected curves, + // we will use the following approach to find axes of the projected ellipse and provide the canonical curve: + // https://www.geometrictools.com/Documentation/ParallelProjectionEllipse.pdf + math_Matrix aMatrA(1, 2, 1, 2); + // A = Jp^T * Pr(Je), where + // Pr(Je) - projection of axes of original ellipse to the target plane + // Jp - X and Y axes of the target plane + aMatrA(1, 1) = myPlane.XDirection().XYZ().Dot(VDx.XYZ()); + aMatrA(1, 2) = myPlane.XDirection().XYZ().Dot(VDy.XYZ()); + aMatrA(2, 1) = myPlane.YDirection().XYZ().Dot(VDx.XYZ()); + aMatrA(2, 2) = myPlane.YDirection().XYZ().Dot(VDy.XYZ()); + + math_Matrix aMatrDelta2(1, 2, 1, 2, 0.0); + // | 1/MajorRad^2 0 | + // Delta^2 = | | + // | 0 1/MajorRad^2 | + aMatrDelta2(1, 1) = 1.0 / (R1 * R1); + aMatrDelta2(2, 2) = 1.0 / (R2 * R2); + + math_Matrix aMatrAInv = aMatrA.Inverse(); + math_Matrix aMatrM = aMatrAInv.Transposed() * aMatrDelta2 * aMatrAInv; + + // perform eigenvalues calculation + math_Jacobi anEigenCalc(aMatrM); + if (anEigenCalc.IsDone()) + { + // radii of the projected ellipse + Minor = 1.0 / Sqrt(anEigenCalc.Value(1)); + Major = 1.0 / Sqrt(anEigenCalc.Value(2)); + + // calculate the rotation angle for the plane axes to meet the correct axes of the projected ellipse + // (swap eigenvectors in respect to major and minor axes) + const math_Matrix& anEigenVec = anEigenCalc.Vectors(); + gp_Trsf2d aTrsfInPlane; + aTrsfInPlane.SetValues(anEigenVec(1, 2), anEigenVec(1, 1), 0.0, + anEigenVec(2, 2), anEigenVec(2, 1), 0.0); + gp_Trsf aRot; + aRot.SetRotation(gp_Ax1(P, myPlane.Direction()), aTrsfInPlane.RotationPart()); + + Dx = myPlane.XDirection().Transformed(aRot); + Dy = myPlane.YDirection().Transformed(aRot); + } + else + { + myIsApprox = Standard_True; + } } if (!myIsApprox) { - Dx = gp_Dir(VDx); - Dy = gp_Dir(VDy); - gp_Pnt O = Axis.Location(); - gp_Pnt P = ProjectPnt(myPlane,myDirection,O); - gp_Pnt Px = ProjectPnt(myPlane,myDirection,O.Translated(R1*gp_Vec(X))); - gp_Pnt Py = ProjectPnt(myPlane,myDirection,O.Translated(R2*gp_Vec(Y))); - Standard_Real Major = P.Distance(Px); - Standard_Real Minor = P.Distance(Py); + gp_Ax2 Axe(P, Dx^Dy, Dx); - if (myKeepParam) - { - myIsApprox = !gp_Dir(VDx).IsNormal(gp_Dir(VDy), Precision::Angular()); + if (Abs(Major - Minor) < Precision::Confusion()) { + myType = GeomAbs_Circle; + gp_Circ Circ(Axe, Major); + GeomCirclePtr = new Geom_Circle(Circ); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin + GeomAdaptor_Curve aGACurve(GeomCirclePtr); + myResult = new GeomAdaptor_Curve(aGACurve); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End } - else - { - // Since it is not necessary to keep the same parameter for the point on the original and on the projected curves, - // we will use the following approach to find axes of the projected ellipse and provide the canonical curve: - // https://www.geometrictools.com/Documentation/ParallelProjectionEllipse.pdf - math_Matrix aMatrA(1, 2, 1, 2); - // A = Jp^T * Pr(Je), where - // Pr(Je) - projection of axes of original ellipse to the target plane - // Jp - X and Y axes of the target plane - aMatrA(1, 1) = myPlane.XDirection().XYZ().Dot(VDx.XYZ()); - aMatrA(1, 2) = myPlane.XDirection().XYZ().Dot(VDy.XYZ()); - aMatrA(2, 1) = myPlane.YDirection().XYZ().Dot(VDx.XYZ()); - aMatrA(2, 2) = myPlane.YDirection().XYZ().Dot(VDy.XYZ()); + else if (Major > Minor) { + myType = GeomAbs_Ellipse; + Elips = gp_Elips(Axe, Major, Minor); - math_Matrix aMatrDelta2(1, 2, 1, 2, 0.0); - // | 1/MajorRad^2 0 | - // Delta^2 = | | - // | 0 1/MajorRad^2 | - aMatrDelta2(1, 1) = 1.0 / (R1 * R1); - aMatrDelta2(2, 2) = 1.0 / (R2 * R2); - - math_Matrix aMatrAInv = aMatrA.Inverse(); - math_Matrix aMatrM = aMatrAInv.Transposed() * aMatrDelta2 * aMatrAInv; - - // perform eigenvalues calculation - math_Jacobi anEigenCalc(aMatrM); - if (anEigenCalc.IsDone()) - { - // radii of the projected ellipse - Minor = 1.0 / Sqrt(anEigenCalc.Value(1)); - Major = 1.0 / Sqrt(anEigenCalc.Value(2)); - - // calculate the rotation angle for the plane axes to meet the correct axes of the projected ellipse - // (swap eigenvectors in respect to major and minor axes) - const math_Matrix& anEigenVec = anEigenCalc.Vectors(); - gp_Trsf2d aTrsfInPlane; - aTrsfInPlane.SetValues(anEigenVec(1, 2), anEigenVec(1, 1), 0.0, - anEigenVec(2, 2), anEigenVec(2, 1), 0.0); - gp_Trsf aRot; - aRot.SetRotation(gp_Ax1(P, myPlane.Direction()), aTrsfInPlane.RotationPart()); - - Dx = myPlane.XDirection().Transformed(aRot); - Dy = myPlane.YDirection().Transformed(aRot); - } - else - { - myIsApprox = Standard_True; - } + GeomEllipsePtr = new Geom_Ellipse(Elips); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin + GeomAdaptor_Curve aGACurve(GeomEllipsePtr); + myResult = new GeomAdaptor_Curve(aGACurve); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End } - - if (!myIsApprox) - { - gp_Ax2 Axe(P, Dx^Dy, Dx); - - if (Abs(Major - Minor) < Precision::Confusion()) { - myType = GeomAbs_Circle; - gp_Circ Circ(Axe, Major); - GeomCirclePtr = new Geom_Circle(Circ); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(GeomCirclePtr); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - } - else if ( Major > Minor) { - myType = GeomAbs_Ellipse; - Elips = gp_Elips( Axe, Major, Minor); - - GeomEllipsePtr = new Geom_Ellipse(Elips); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(GeomEllipsePtr); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - } - else { - myIsApprox = Standard_True; - } + else { + myIsApprox = Standard_True; } } - - // No way to build the canonical curve, approximate as B-spline - if (myIsApprox) - { - myType = GeomAbs_BSplineCurve; - PerformApprox(myCurve,myPlane,myDirection,ApproxCurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(ApproxCurve); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - } - else if (GeomCirclePtr || GeomEllipsePtr) - { - Handle(Geom_Curve) aResultCurve = GeomCirclePtr; - if (aResultCurve.IsNull()) - aResultCurve = GeomEllipsePtr; - // start and end parameters of the projected curve - Standard_Real aParFirst = myCurve->FirstParameter(); - Standard_Real aParLast = myCurve->LastParameter(); - gp_Pnt aPntFirst = ProjectPnt(myPlane, myDirection, myCurve->Value(aParFirst)); - gp_Pnt aPntLast = ProjectPnt(myPlane, myDirection, myCurve->Value(aParLast)); - GeomLib_Tool::Parameter(aResultCurve, aPntFirst, Precision::Confusion(), myFirstPar); - GeomLib_Tool::Parameter(aResultCurve, aPntLast, Precision::Confusion(), myLastPar); - while (myLastPar <= myFirstPar) - myLastPar += myResult->Period(); - } } - break; - case GeomAbs_Parabola: + + // No way to build the canonical curve, approximate as B-spline + if (myIsApprox) { - myKeepParam = Standard_True; - // P(u) = O + (u*u)/(4*f) * Xc + u * Yc - // ==> Q(u) = f(P(u)) - // = f(O) + (u*u)/(4*f) * f(Xc) + u * f(Yc) - - gp_Parab Parab = myCurve->Parabola(); - gp_Ax2 AxeRef = Parab.Position(); - gp_Vec Xc = ProjectVec(myPlane,myDirection,gp_Vec(AxeRef.XDirection())); - gp_Vec Yc = ProjectVec(myPlane,myDirection,gp_Vec(AxeRef.YDirection())); - - // fix for case when no one action is done. 28.03.2002 - Standard_Boolean alocalIsDone = Standard_False; - - if ( Abs( Yc.Magnitude() - 1.) < Precision::Confusion()) { - gp_Pnt P = ProjectPnt(myPlane,myDirection,AxeRef.Location()); - if ( Xc.Magnitude() < Precision::Confusion()) { - myType = GeomAbs_Line; - aLine = gp_Lin( P, gp_Dir(Yc)); - - GeomLinePtr = new Geom_Line(aLine) ; -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(GeomLinePtr); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - alocalIsDone = Standard_True; - } - else if ( Xc.IsNormal(Yc,Precision::Angular())) { - myType = GeomAbs_Parabola; - Standard_Real F = Parab.Focal() / Xc.Magnitude(); - gp_Parab aParab = gp_Parab( gp_Ax2(P,Xc^Yc,Xc), F); - Handle(Geom_Parabola) GeomParabolaPtr = - new Geom_Parabola(aParab) ; -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(GeomParabolaPtr); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - alocalIsDone = Standard_True; - } - } - if (!alocalIsDone)/*else*/ { - myIsApprox = Standard_True; - myType = GeomAbs_BSplineCurve; - PerformApprox(myCurve,myPlane,myDirection,ApproxCurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(ApproxCurve); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - } - } - break; - case GeomAbs_Hyperbola: - { - myKeepParam = Standard_True; - // P(u) = O + R1 * Cosh(u) * Xc + R2 * Sinh(u) * Yc - // ==> Q(u) = f(P(u)) - // = f(O) + R1 * Cosh(u) * f(Xc) + R2 * Sinh(u) * f(Yc) - - gp_Hypr Hypr = myCurve->Hyperbola(); - gp_Ax2 AxeRef = Hypr.Position(); - gp_Vec Xc = ProjectVec(myPlane,myDirection,gp_Vec(AxeRef.XDirection())); - gp_Vec Yc = ProjectVec(myPlane,myDirection,gp_Vec(AxeRef.YDirection())); - gp_Pnt P = ProjectPnt(myPlane,myDirection,AxeRef.Location()); - Standard_Real aR1 = Hypr.MajorRadius(); - Standard_Real aR2 = Hypr.MinorRadius(); - gp_Dir Z = myPlane.Direction(); - - if ( Xc.Magnitude() < Precision::Confusion()) { - myType = GeomAbs_Hyperbola; - gp_Dir X = gp_Dir(Yc) ^ Z; - Hypr = gp_Hypr(gp_Ax2( P, Z, X), 0., aR2 * Yc.Magnitude()); - GeomHyperbolaPtr = - new Geom_Hyperbola(Hypr) ; -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(GeomHyperbolaPtr); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - } - else if ( Yc.Magnitude() < Precision::Confusion()) { - myType = GeomAbs_Hyperbola; - Hypr = - gp_Hypr(gp_Ax2(P, Z, gp_Dir(Xc)), aR1 * Xc.Magnitude(), 0.); - GeomHyperbolaPtr = - new Geom_Hyperbola(Hypr) ; -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(GeomHyperbolaPtr); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - } - else if ( Xc.IsNormal(Yc,Precision::Angular())) { - myType = GeomAbs_Hyperbola; - Hypr = gp_Hypr( gp_Ax2( P, gp_Dir( Xc ^ Yc), gp_Dir( Xc)), - aR1 * Xc.Magnitude(), aR2 * Yc.Magnitude() ); - GeomHyperbolaPtr = - new Geom_Hyperbola(Hypr) ; -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(GeomHyperbolaPtr); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - } - else { - myIsApprox = Standard_True; - myType = GeomAbs_BSplineCurve; - PerformApprox(myCurve,myPlane,myDirection,ApproxCurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(ApproxCurve); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - } - } - break; - case GeomAbs_BezierCurve: - { - Handle(Geom_BezierCurve) BezierCurvePtr = - myCurve->Bezier() ; - Standard_Integer NbPoles = - BezierCurvePtr->NbPoles() ; - - Handle(Geom_BezierCurve) ProjCu = - Handle(Geom_BezierCurve)::DownCast(BezierCurvePtr->Copy()); - - myKeepParam = Standard_True; - myIsApprox = Standard_False; - myType = Type; - for ( Standard_Integer i = 1; i <= NbPoles; i++) { - ProjCu->SetPole - (i,ProjectPnt(myPlane,myDirection,BezierCurvePtr->Pole(i))); - } - -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(ProjCu); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - } - break ; - case GeomAbs_BSplineCurve: - { - Handle(Geom_BSplineCurve) BSplineCurvePtr = - myCurve->BSpline() ; - // - // make a copy of the curve and projects its poles - // - Handle(Geom_BSplineCurve) ProjectedBSplinePtr = - Handle(Geom_BSplineCurve)::DownCast(BSplineCurvePtr->Copy()) ; - - myKeepParam = Standard_True; - myIsApprox = Standard_False; - myType = Type; - for ( Standard_Integer i = 1; i <= BSplineCurvePtr->NbPoles(); i++) { - ProjectedBSplinePtr->SetPole - (i,ProjectPnt(myPlane,myDirection,BSplineCurvePtr->Pole(i))); - } - -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin - GeomAdaptor_Curve aGACurve(ProjectedBSplinePtr); - myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End - } - break; - default: - { - myKeepParam = Standard_True; - myIsApprox = Standard_True; - myType = GeomAbs_BSplineCurve; - PerformApprox(myCurve,myPlane,myDirection,ApproxCurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin + myType = GeomAbs_BSplineCurve; + PerformApprox(myCurve, myPlane, myDirection, ApproxCurve); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin GeomAdaptor_Curve aGACurve(ApproxCurve); myResult = new GeomAdaptor_Curve(aGACurve); -// Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End } - break; + else if (GeomCirclePtr || GeomEllipsePtr) + { + Handle(Geom_Curve) aResultCurve = GeomCirclePtr; + if (aResultCurve.IsNull()) + aResultCurve = GeomEllipsePtr; + // start and end parameters of the projected curve + Standard_Real aParFirst = myCurve->FirstParameter(); + Standard_Real aParLast = myCurve->LastParameter(); + gp_Pnt aPntFirst = ProjectPnt(myPlane, myDirection, myCurve->Value(aParFirst)); + gp_Pnt aPntLast = ProjectPnt(myPlane, myDirection, myCurve->Value(aParLast)); + GeomLib_Tool::Parameter(aResultCurve, aPntFirst, Precision::Confusion(), myFirstPar); + GeomLib_Tool::Parameter(aResultCurve, aPntLast, Precision::Confusion(), myLastPar); + while (myLastPar <= myFirstPar) + myLastPar += myResult->Period(); + } + } + break; + case GeomAbs_Parabola: + { + // P(u) = O + (u*u)/(4*f) * Xc + u * Yc + // ==> Q(u) = f(P(u)) + // = f(O) + (u*u)/(4*f) * f(Xc) + u * f(Yc) + + gp_Parab Parab = myCurve->Parabola(); + gp_Ax2 AxeRef = Parab.Position(); + gp_Vec Xc = ProjectVec(myPlane, myDirection, gp_Vec(AxeRef.XDirection())); + gp_Vec Yc = ProjectVec(myPlane, myDirection, gp_Vec(AxeRef.YDirection())); + gp_Pnt P = ProjectPnt(myPlane, myDirection, AxeRef.Location()); + + myIsApprox = Standard_False; + + if ((Abs(Yc.Magnitude() - 1.) < Precision::Confusion()) && + (Xc.Magnitude() < Precision::Confusion())) + { + myType = GeomAbs_Line; + aLine = gp_Lin(P, gp_Dir(Yc)); + GeomLinePtr = new Geom_Line(aLine); + } + else if (Xc.IsNormal(Yc, Precision::Angular())) { + myType = GeomAbs_Parabola; + Standard_Real F = Parab.Focal() / Xc.Magnitude(); + gp_Parab aProjParab = gp_Parab(gp_Ax2(P, Xc^Yc, Xc), F); + GeomParabolaPtr = + new Geom_Parabola(aProjParab); + } + else if (Yc.Magnitude() < Precision::Confusion() || + Yc.IsParallel(Xc, Precision::Angular())) + { + myIsApprox = Standard_True; + } + else if(!myKeepParam) + { + // Try building parabola with help of apex position + myIsApprox = !BuildParabolaByApex(GeomParabolaPtr); + } + else + { + myIsApprox = Standard_True; + } + + if (!myIsApprox) + { + GetTrimmedResult(GeomParabolaPtr); + } + else + { + BuildByApprox(aParabolaLimit); + } + } + break; + case GeomAbs_Hyperbola: + { + // P(u) = O + R1 * Cosh(u) * Xc + R2 * Sinh(u) * Yc + // ==> Q(u) = f(P(u)) + // = f(O) + R1 * Cosh(u) * f(Xc) + R2 * Sinh(u) * f(Yc) + + gp_Hypr Hypr = myCurve->Hyperbola(); + gp_Ax2 AxeRef = Hypr.Position(); + gp_Vec Xc = ProjectVec(myPlane, myDirection, gp_Vec(AxeRef.XDirection())); + gp_Vec Yc = ProjectVec(myPlane, myDirection, gp_Vec(AxeRef.YDirection())); + gp_Pnt P = ProjectPnt(myPlane, myDirection, AxeRef.Location()); + Standard_Real aR1 = Hypr.MajorRadius(); + Standard_Real aR2 = Hypr.MinorRadius(); + gp_Dir Z = myPlane.Direction(); + myIsApprox = Standard_False; + + if (Xc.Magnitude() < Precision::Confusion()) { + myType = GeomAbs_Hyperbola; + gp_Dir X = gp_Dir(Yc) ^ Z; + Hypr = gp_Hypr(gp_Ax2(P, Z, X), 0., aR2 * Yc.Magnitude()); + GeomHyperbolaPtr = + new Geom_Hyperbola(Hypr); + } + else if (Yc.Magnitude() < Precision::Confusion()) { + myType = GeomAbs_Hyperbola; + Hypr = + gp_Hypr(gp_Ax2(P, Z, gp_Dir(Xc)), aR1 * Xc.Magnitude(), 0.); + GeomHyperbolaPtr = + new Geom_Hyperbola(Hypr); + } + else if (Xc.IsNormal(Yc, Precision::Angular())) { + myType = GeomAbs_Hyperbola; + Hypr = gp_Hypr(gp_Ax2(P, gp_Dir(Xc ^ Yc), gp_Dir(Xc)), + aR1 * Xc.Magnitude(), aR2 * Yc.Magnitude()); + GeomHyperbolaPtr = + new Geom_Hyperbola(Hypr); + } + else if (Yc.Magnitude() < Precision::Confusion() || + Yc.IsParallel(Xc, Precision::Angular())) + { + myIsApprox = Standard_True; + } + else if(!myKeepParam) + { + myIsApprox = !BuildHyperbolaByApex(GeomHyperbolaPtr); + } + else + { + myIsApprox = Standard_True; + } + if ( !myIsApprox ) + { + GetTrimmedResult(GeomHyperbolaPtr); + } + else + { + BuildByApprox(aHyperbolaLimit); + } + } + break; + case GeomAbs_BezierCurve: + { + Handle(Geom_BezierCurve) BezierCurvePtr = + myCurve->Bezier(); + Standard_Integer NbPoles = + BezierCurvePtr->NbPoles(); + + Handle(Geom_BezierCurve) ProjCu = + Handle(Geom_BezierCurve)::DownCast(BezierCurvePtr->Copy()); + + myKeepParam = Standard_True; + myIsApprox = Standard_False; + myType = Type; + for (Standard_Integer i = 1; i <= NbPoles; i++) { + ProjCu->SetPole + (i, ProjectPnt(myPlane, myDirection, BezierCurvePtr->Pole(i))); + } + + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin + GeomAdaptor_Curve aGACurve(ProjCu); + myResult = new GeomAdaptor_Curve(aGACurve); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End + } + break; + case GeomAbs_BSplineCurve: + { + Handle(Geom_BSplineCurve) BSplineCurvePtr = + myCurve->BSpline(); + // + // make a copy of the curve and projects its poles + // + Handle(Geom_BSplineCurve) ProjectedBSplinePtr = + Handle(Geom_BSplineCurve)::DownCast(BSplineCurvePtr->Copy()); + + myKeepParam = Standard_True; + myIsApprox = Standard_False; + myType = Type; + for (Standard_Integer i = 1; i <= BSplineCurvePtr->NbPoles(); i++) { + ProjectedBSplinePtr->SetPole + (i, ProjectPnt(myPlane, myDirection, BSplineCurvePtr->Pole(i))); + } + + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin + GeomAdaptor_Curve aGACurve(ProjectedBSplinePtr); + myResult = new GeomAdaptor_Curve(aGACurve); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End + } + break; + default: + { + myKeepParam = Standard_True; + myIsApprox = Standard_True; + myType = GeomAbs_BSplineCurve; + PerformApprox(myCurve, myPlane, myDirection, ApproxCurve); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin + GeomAdaptor_Curve aGACurve(ApproxCurve); + myResult = new GeomAdaptor_Curve(aGACurve); + // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End + } + break; } } @@ -984,7 +1048,7 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_Curve)& C, //purpose : //======================================================================= -const gp_Ax3& ProjLib_ProjectOnPlane::GetPlane() const +const gp_Ax3& ProjLib_ProjectOnPlane::GetPlane() const { return myPlane; } @@ -994,7 +1058,7 @@ const gp_Ax3& ProjLib_ProjectOnPlane::GetPlane() const //purpose : //======================================================================= -const gp_Dir& ProjLib_ProjectOnPlane::GetDirection() const +const gp_Dir& ProjLib_ProjectOnPlane::GetDirection() const { return myDirection; } @@ -1025,9 +1089,9 @@ const Handle(GeomAdaptor_Curve)& ProjLib_ProjectOnPlane::GetResult() const //purpose : //======================================================================= -Standard_Real ProjLib_ProjectOnPlane::FirstParameter() const +Standard_Real ProjLib_ProjectOnPlane::FirstParameter() const { - if ( myKeepParam || myIsApprox) + if (myKeepParam || myIsApprox) return myCurve->FirstParameter(); else return myFirstPar; @@ -1039,9 +1103,9 @@ Standard_Real ProjLib_ProjectOnPlane::FirstParameter() const //purpose : //======================================================================= -Standard_Real ProjLib_ProjectOnPlane::LastParameter() const +Standard_Real ProjLib_ProjectOnPlane::LastParameter() const { - if ( myKeepParam || myIsApprox) + if (myKeepParam || myIsApprox) return myCurve->LastParameter(); else return myLastPar; @@ -1055,7 +1119,7 @@ Standard_Real ProjLib_ProjectOnPlane::LastParameter() const GeomAbs_Shape ProjLib_ProjectOnPlane::Continuity() const { - return myCurve->Continuity() ; + return myCurve->Continuity(); } @@ -1066,7 +1130,7 @@ GeomAbs_Shape ProjLib_ProjectOnPlane::Continuity() const Standard_Integer ProjLib_ProjectOnPlane::NbIntervals(const GeomAbs_Shape S) const { - return myCurve->NbIntervals(S) ; + return myCurve->NbIntervals(S); } @@ -1075,10 +1139,10 @@ Standard_Integer ProjLib_ProjectOnPlane::NbIntervals(const GeomAbs_Shape S) cons //purpose : //======================================================================= -void ProjLib_ProjectOnPlane::Intervals(TColStd_Array1OfReal& T, - const GeomAbs_Shape S) const +void ProjLib_ProjectOnPlane::Intervals(TColStd_Array1OfReal& T, + const GeomAbs_Shape S) const { - myCurve->Intervals(T,S) ; + myCurve->Intervals(T, S); } //======================================================================= @@ -1086,20 +1150,20 @@ void ProjLib_ProjectOnPlane::Intervals(TColStd_Array1OfReal& T, //purpose : //======================================================================= -Handle(Adaptor3d_Curve) +Handle(Adaptor3d_Curve) ProjLib_ProjectOnPlane::Trim(const Standard_Real First, - const Standard_Real Last, - const Standard_Real Tolerance) const + const Standard_Real Last, + const Standard_Real Tolerance) const { - if (myType != GeomAbs_OtherCurve){ - return myResult->Trim(First,Last,Tolerance) ; + if (myType != GeomAbs_OtherCurve) { + return myResult->Trim(First, Last, Tolerance); } else { - throw Standard_NotImplemented ("ProjLib_ProjectOnPlane::Trim() - curve of unsupported type"); + throw Standard_NotImplemented("ProjLib_ProjectOnPlane::Trim() - curve of unsupported type"); } } - + //======================================================================= //function : IsClosed //purpose : @@ -1107,7 +1171,7 @@ ProjLib_ProjectOnPlane::Trim(const Standard_Real First, Standard_Boolean ProjLib_ProjectOnPlane::IsClosed() const { - return myCurve->IsClosed() ; + return myCurve->IsClosed(); } @@ -1118,9 +1182,9 @@ Standard_Boolean ProjLib_ProjectOnPlane::IsClosed() const Standard_Boolean ProjLib_ProjectOnPlane::IsPeriodic() const { - if ( myIsApprox) + if (myIsApprox) return Standard_False; - else + else return myCurve->IsPeriodic(); } @@ -1132,13 +1196,13 @@ Standard_Boolean ProjLib_ProjectOnPlane::IsPeriodic() const Standard_Real ProjLib_ProjectOnPlane::Period() const { - if ( !IsPeriodic()) { + if (!IsPeriodic()) { throw Standard_NoSuchObject("ProjLib_ProjectOnPlane::Period"); } - - if ( myIsApprox) + + if (myIsApprox) return Standard_False; - else + else return myCurve->Period(); } @@ -1148,19 +1212,19 @@ Standard_Real ProjLib_ProjectOnPlane::Period() const //purpose : //======================================================================= -gp_Pnt ProjLib_ProjectOnPlane::Value(const Standard_Real U) const +gp_Pnt ProjLib_ProjectOnPlane::Value(const Standard_Real U) const { - if (myType != GeomAbs_OtherCurve) { + if (myType != GeomAbs_OtherCurve) { return myResult->Value(U); } else { return OnPlane_Value(U, - myCurve, - myPlane, - myDirection); - + myCurve, + myPlane, + myDirection); + } -} +} //======================================================================= @@ -1168,16 +1232,16 @@ gp_Pnt ProjLib_ProjectOnPlane::Value(const Standard_Real U) const //purpose : //======================================================================= -void ProjLib_ProjectOnPlane::D0(const Standard_Real U , gp_Pnt& P) const +void ProjLib_ProjectOnPlane::D0(const Standard_Real U, gp_Pnt& P) const { if (myType != GeomAbs_OtherCurve) { - myResult->D0(U,P) ; + myResult->D0(U, P); } else { P = OnPlane_Value(U, - myCurve, - myPlane, - myDirection); + myCurve, + myPlane, + myDirection); } } @@ -1188,19 +1252,19 @@ void ProjLib_ProjectOnPlane::D0(const Standard_Real U , gp_Pnt& P) const //======================================================================= void ProjLib_ProjectOnPlane::D1(const Standard_Real U, - gp_Pnt& P , - gp_Vec& V ) const + gp_Pnt& P, + gp_Vec& V) const { if (myType != GeomAbs_OtherCurve) { - myResult->D1(U,P,V) ; + myResult->D1(U, P, V); } else { OnPlane_D1(U, - P, - V, - myCurve, - myPlane, - myDirection); + P, + V, + myCurve, + myPlane, + myDirection); } } @@ -1210,22 +1274,22 @@ void ProjLib_ProjectOnPlane::D1(const Standard_Real U, //purpose : //======================================================================= -void ProjLib_ProjectOnPlane::D2(const Standard_Real U, - gp_Pnt& P, - gp_Vec& V1, - gp_Vec& V2) const +void ProjLib_ProjectOnPlane::D2(const Standard_Real U, + gp_Pnt& P, + gp_Vec& V1, + gp_Vec& V2) const { - if (myType != GeomAbs_OtherCurve) { - myResult->D2(U,P,V1,V2) ; + if (myType != GeomAbs_OtherCurve) { + myResult->D2(U, P, V1, V2); } else { OnPlane_D2(U, - P, - V1, - V2, - myCurve, - myPlane, - myDirection); + P, + V1, + V2, + myCurve, + myPlane, + myDirection); } } @@ -1235,24 +1299,24 @@ void ProjLib_ProjectOnPlane::D2(const Standard_Real U, //purpose : //======================================================================= -void ProjLib_ProjectOnPlane::D3(const Standard_Real U, - gp_Pnt& P, - gp_Vec& V1, - gp_Vec& V2, - gp_Vec& V3) const +void ProjLib_ProjectOnPlane::D3(const Standard_Real U, + gp_Pnt& P, + gp_Vec& V1, + gp_Vec& V2, + gp_Vec& V3) const { - if (myType != GeomAbs_OtherCurve) { - myResult->D3(U,P,V1,V2,V3) ; + if (myType != GeomAbs_OtherCurve) { + myResult->D3(U, P, V1, V2, V3); } - else { + else { OnPlane_D3(U, - P, - V1, - V2, - V3, - myCurve, - myPlane, - myDirection); + P, + V1, + V2, + V3, + myCurve, + myPlane, + myDirection); } } @@ -1262,20 +1326,20 @@ void ProjLib_ProjectOnPlane::D3(const Standard_Real U, //purpose : //======================================================================= -gp_Vec ProjLib_ProjectOnPlane::DN(const Standard_Real U, - const Standard_Integer DerivativeRequest) - const +gp_Vec ProjLib_ProjectOnPlane::DN(const Standard_Real U, + const Standard_Integer DerivativeRequest) + const { if (myType != GeomAbs_OtherCurve) { - return myResult->DN(U,DerivativeRequest) ; + return myResult->DN(U, DerivativeRequest); } else { return OnPlane_DN(U, - DerivativeRequest, - myCurve, - myPlane, - myDirection); - } + DerivativeRequest, + myCurve, + myPlane, + myDirection); + } } @@ -1285,16 +1349,16 @@ gp_Vec ProjLib_ProjectOnPlane::DN(const Standard_Real U, //======================================================================= Standard_Real ProjLib_ProjectOnPlane::Resolution -(const Standard_Real Tolerance) const +(const Standard_Real Tolerance) const { if (myType != GeomAbs_OtherCurve) { - return myResult->Resolution(Tolerance) ; + return myResult->Resolution(Tolerance); } else { return 0; } } - + //======================================================================= //function : GetType @@ -1359,7 +1423,7 @@ gp_Hypr ProjLib_ProjectOnPlane::Hyperbola() const if (myType != GeomAbs_Hyperbola) throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:Hyperbola"); - return myResult->Hyperbola() ; + return myResult->Hyperbola(); } @@ -1372,8 +1436,8 @@ gp_Parab ProjLib_ProjectOnPlane::Parabola() const { if (myType != GeomAbs_Parabola) throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:Parabola"); - - return myResult->Parabola() ; + + return myResult->Parabola(); } //======================================================================= @@ -1384,10 +1448,10 @@ gp_Parab ProjLib_ProjectOnPlane::Parabola() const Standard_Integer ProjLib_ProjectOnPlane::Degree() const { if ((GetType() != GeomAbs_BSplineCurve) && - (GetType() != GeomAbs_BezierCurve)) + (GetType() != GeomAbs_BezierCurve)) throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:Degree"); - if ( myIsApprox) + if (myIsApprox) return myResult->Degree(); else return myCurve->Degree(); @@ -1398,13 +1462,13 @@ Standard_Integer ProjLib_ProjectOnPlane::Degree() const //purpose : //======================================================================= -Standard_Boolean ProjLib_ProjectOnPlane::IsRational() const +Standard_Boolean ProjLib_ProjectOnPlane::IsRational() const { if ((GetType() != GeomAbs_BSplineCurve) && - (GetType() != GeomAbs_BezierCurve)) + (GetType() != GeomAbs_BezierCurve)) throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:IsRational"); - - if ( myIsApprox) + + if (myIsApprox) return myResult->IsRational(); else return myCurve->IsRational(); @@ -1418,10 +1482,10 @@ Standard_Boolean ProjLib_ProjectOnPlane::IsRational() const Standard_Integer ProjLib_ProjectOnPlane::NbPoles() const { if ((GetType() != GeomAbs_BSplineCurve) && - (GetType() != GeomAbs_BezierCurve)) + (GetType() != GeomAbs_BezierCurve)) throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:NbPoles"); - - if ( myIsApprox) + + if (myIsApprox) return myResult->NbPoles(); else return myCurve->NbPoles(); @@ -1432,12 +1496,12 @@ Standard_Integer ProjLib_ProjectOnPlane::NbPoles() const //purpose : //======================================================================= -Standard_Integer ProjLib_ProjectOnPlane::NbKnots() const +Standard_Integer ProjLib_ProjectOnPlane::NbKnots() const { - if ( GetType() != GeomAbs_BSplineCurve) + if (GetType() != GeomAbs_BSplineCurve) throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:NbKnots"); - - if ( myIsApprox) + + if (myIsApprox) return myResult->NbKnots(); else return myCurve->NbKnots(); @@ -1454,11 +1518,11 @@ Handle(Geom_BezierCurve) ProjLib_ProjectOnPlane::Bezier() const if (myType != GeomAbs_BezierCurve) throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:Bezier"); - return myResult->Bezier() ; + return myResult->Bezier(); } //======================================================================= -//function : Bezier +//function : BSpline //purpose : //======================================================================= @@ -1467,6 +1531,225 @@ Handle(Geom_BSplineCurve) ProjLib_ProjectOnPlane::BSpline() const if (myType != GeomAbs_BSplineCurve) throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:BSpline"); - return myResult->BSpline() ; + return myResult->BSpline(); } +//======================================================================= +//function : GetTrimmedResult +//purpose : +//======================================================================= + +void ProjLib_ProjectOnPlane::GetTrimmedResult(const Handle(Geom_Curve)& theProjCurve) +{ + gp_Lin aLin; + gp_Parab aParab; + gp_Hypr aHypr; + if (myType == GeomAbs_Line) + { + aLin = Handle(Geom_Line)::DownCast(theProjCurve)->Lin(); + } + else if (myType == GeomAbs_Parabola) + { + aParab = Handle(Geom_Parabola)::DownCast(theProjCurve)->Parab(); + } + else if (myType == GeomAbs_Hyperbola) + { + aHypr = Handle(Geom_Hyperbola)::DownCast(theProjCurve)->Hypr(); + } + + myFirstPar = theProjCurve->FirstParameter(); + myLastPar = theProjCurve->LastParameter(); + if (!Precision::IsInfinite(myCurve->FirstParameter())) + { + gp_Pnt aP = myCurve->Value(myCurve->FirstParameter()); + aP = ProjectPnt(myPlane, myDirection, aP); + if (myType == GeomAbs_Line) + { + myFirstPar = ElCLib::Parameter(aLin, aP); + } + else if (myType == GeomAbs_Parabola) + { + myFirstPar = ElCLib::Parameter(aParab, aP); + } + else if (myType == GeomAbs_Hyperbola) + { + myFirstPar = ElCLib::Parameter(aHypr, aP); + } + else + { + GeomLib_Tool::Parameter(theProjCurve, aP, Precision::Confusion(), myFirstPar); + } + } + if (!Precision::IsInfinite(myCurve->LastParameter())) + { + gp_Pnt aP = myCurve->Value(myCurve->LastParameter()); + aP = ProjectPnt(myPlane, myDirection, aP); + if (myType == GeomAbs_Line) + { + myLastPar = ElCLib::Parameter(aLin, aP); + } + else if (myType == GeomAbs_Parabola) + { + myLastPar = ElCLib::Parameter(aParab, aP); + } + else if (myType == GeomAbs_Hyperbola) + { + myLastPar = ElCLib::Parameter(aHypr, aP); + } + else + { + GeomLib_Tool::Parameter(theProjCurve, aP, Precision::Confusion(), myLastPar); + } + } + myResult = new GeomAdaptor_Curve(theProjCurve, myFirstPar, myLastPar); + +} + +//======================================================================= +//function : BuildParabolaByApex +//purpose : +//======================================================================= + +Standard_Boolean ProjLib_ProjectOnPlane::BuildParabolaByApex(Handle(Geom_Curve)& theGeomParabolaPtr) +{ + // + //Searching parabola apex as point with maximal curvature + Standard_Real aF = myCurve->Parabola().Focal(); + GeomAbs_CurveType aCurType = myType; + myType = GeomAbs_OtherCurve; //To provide correct calculation of derivativesb by projection for + //copy of instance; + Handle(Adaptor3d_Curve) aProjCrv = ShallowCopy(); + myType = aCurType; + LProp3d_CLProps aProps(aProjCrv, 2, Precision::Confusion()); + ProjLib_MaxCurvature aMaxCur(aProps); + math_BrentMinimum aSolver(Precision::PConfusion()); + aSolver.Perform(aMaxCur, -10.*aF, 0., 10.*aF); + + if (!aSolver.IsDone()) + { + return Standard_False; + } + + Standard_Real aT; + aT = aSolver.Location(); + aProps.SetParameter(aT); + gp_Pnt aP0 = aProps.Value(); + gp_Vec aDY = aProps.D1(); + gp_Dir anYDir(aDY); + gp_Dir anXDir; + Standard_Real aCurv = aProps.Curvature(); + if (Precision::IsInfinite(aCurv) || aCurv < Precision::Confusion()) + { + return Standard_False; + } + aProps.Normal(anXDir); + // + gp_Lin anXLine(aP0, anXDir); + gp_Pnt aP1 = Value(aT + 10.*aF); + // + Standard_Real anX = ElCLib::LineParameter(anXLine.Position(), aP1); + Standard_Real anY = anXLine.Distance(aP1); + Standard_Real aNewF = anY * anY / 4. / anX; + gp_Dir anN = anXDir^anYDir; + gp_Ax2 anA2(aP0, anN, anXDir); + gce_MakeParab aMkParab(anA2, aNewF); + if (!aMkParab.IsDone()) + { + return Standard_False; + } + + gp_Parab aProjParab = aMkParab.Value(); + + myType = GeomAbs_Parabola; + theGeomParabolaPtr = new Geom_Parabola(aProjParab); + //GetTrimmedResult(theGeomParabolaPtr); + + return Standard_True; +} + +//======================================================================= +//function : BuildHyperbolaByApex +//purpose : +//======================================================================= + +Standard_Boolean ProjLib_ProjectOnPlane::BuildHyperbolaByApex(Handle(Geom_Curve)& theGeomHyperbolaPtr) +{ + //Try to build hyperbola with help of apex position + GeomAbs_CurveType aCurType = myType; + myType = GeomAbs_OtherCurve; //To provide correct calculation of derivativesb by projection for + //copy of instance; + Handle(Adaptor3d_Curve) aProjCrv = ShallowCopy(); + myType = aCurType; + //Searching hyperbola apex as point with maximal curvature + LProp3d_CLProps aProps(aProjCrv, 2, Precision::Confusion()); + ProjLib_MaxCurvature aMaxCur(aProps); + math_BrentMinimum aSolver(Precision::PConfusion()); + aSolver.Perform(aMaxCur, -5., 0., 5.); + + if (aSolver.IsDone()) + { + Standard_Real aT; + aT = aSolver.Location(); + aProps.SetParameter(aT); + Standard_Real aCurv = aProps.Curvature(); + if (Precision::IsInfinite(aCurv) || aCurv < Precision::Confusion()) + { + return Standard_False; + } + else + { + gp_Hypr Hypr = myCurve->Hyperbola(); + gp_Ax2 AxeRef = Hypr.Position(); + gp_Pnt P = ProjectPnt(myPlane, myDirection, AxeRef.Location()); + gp_Dir Z = myPlane.Direction(); + gp_Pnt aP0 = aProps.Value(); + gp_Dir anXDir = gce_MakeDir(P, aP0); + gp_Dir anYDir = gce_MakeDir(aProps.D1()); + // + Standard_Real aMajRad = P.Distance(aP0); + gp_Pnt aP1 = Value(aT + 1.); + gp_Vec aV(P, aP1); + Standard_Real anX = aV * anXDir; + Standard_Real anY = aV * anYDir; + Standard_Real aMinRad = anY / Sqrt(anX * anX / aMajRad / aMajRad - 1.); + gp_Ax2 anA2(P, Z, anXDir); + gp_Hypr anHypr(anA2, aMajRad, aMinRad); + theGeomHyperbolaPtr = + new Geom_Hyperbola(anHypr); + myType = GeomAbs_Hyperbola; + } + } + else + { + return Standard_False; + } + return Standard_True; +} + +//======================================================================= +//function : BuilByApprox +//purpose : +//======================================================================= + +void ProjLib_ProjectOnPlane::BuildByApprox(const Standard_Real theLimitParameter) +{ + myType = GeomAbs_BSplineCurve; + Handle(Geom_BSplineCurve) anApproxCurve; + if (Precision::IsInfinite(myCurve->FirstParameter()) || + Precision::IsInfinite(myCurve->LastParameter())) + { + //To avoid exception in approximation + Standard_Real f = Max(-theLimitParameter, myCurve->FirstParameter()); + Standard_Real l = Min(theLimitParameter, myCurve->LastParameter()); + Handle(Adaptor3d_Curve) aTrimCurve = myCurve->Trim(f, l, Precision::Confusion()); + PerformApprox(aTrimCurve, myPlane, myDirection, anApproxCurve); + } + else + { + PerformApprox(myCurve, myPlane, myDirection, anApproxCurve); + } + myFirstPar = anApproxCurve->FirstParameter(); + myLastPar = anApproxCurve->LastParameter(); + GeomAdaptor_Curve aGACurve(anApproxCurve); + myResult = new GeomAdaptor_Curve(aGACurve); +} \ No newline at end of file diff --git a/src/ProjLib/ProjLib_ProjectOnPlane.hxx b/src/ProjLib/ProjLib_ProjectOnPlane.hxx index d885bb5ee8..ec94cbb7b1 100644 --- a/src/ProjLib/ProjLib_ProjectOnPlane.hxx +++ b/src/ProjLib/ProjLib_ProjectOnPlane.hxx @@ -192,9 +192,12 @@ public: protected: + void GetTrimmedResult(const Handle(Geom_Curve)& theProjCurve); + Standard_Boolean BuildParabolaByApex(Handle(Geom_Curve)& theGeomParabolaPtr); + Standard_Boolean BuildHyperbolaByApex(Handle(Geom_Curve)& theGeomParabolaPtr); - + void BuildByApprox(const Standard_Real theLimitParameter); private: diff --git a/tests/bugs/moddata_3/bug31661_1 b/tests/bugs/moddata_3/bug31661_1 new file mode 100644 index 0000000000..b30ef8fe2a --- /dev/null +++ b/tests/bugs/moddata_3/bug31661_1 @@ -0,0 +1,26 @@ +puts "" +puts "==========================================================================" +puts "OCC31661: Modeling Data - Algorithm crashes when projecting parabola or hyperbola to plane" +puts "==========================================================================" +puts "" + +parabola p 0 0 0 1 1 1 2 0 -2 10 +plane pln 0 0 0 0 0 1 +projonplane r p pln 0 + +if {![regexp {Parabola} [dump r]]} { + puts "ERROR: Projected curve is not a parabola" +} + +trim p p -100 100 +projonplane rp p pln 0 + +if {![regexp {Parabola} [dump rp]]} { + puts "ERROR: Projected curve is not a parabola" +} + +checklength rp -l 408.40363195229503 +bounds rp t1 t2 +checkreal t1 [dval t1] -91.077748943768597 1.e-7 1.e-7 +checkreal t2 [dval t2] 72.221567418462357 1.e-7 1.e-7 + diff --git a/tests/bugs/moddata_3/bug31661_2 b/tests/bugs/moddata_3/bug31661_2 new file mode 100644 index 0000000000..0e601af77e --- /dev/null +++ b/tests/bugs/moddata_3/bug31661_2 @@ -0,0 +1,26 @@ +puts "" +puts "==========================================================================" +puts "OCC31661: Modeling Data - Algorithm crashes when projecting parabola or hyperbola to plane" +puts "==========================================================================" +puts "" + +hyperbola h 0 0 0 1 1 1 2 0 -2 10 10 +plane pln 0 0 0 0 0 1 +projonplane rh h pln 0 + +if {![regexp {Hyperbola} [dump rh]]} { + puts "ERROR: Projected curve is not a hyperbola" +} + +trim h h -5 5 +projonplane rh h pln 0 + +if {![regexp {Hyperbola} [dump rh]]} { + puts "ERROR: Projected curve is not a hyperbola" +} + +checklength rh -l 1664.3732976598988 +bounds rh t1 t2 +checkreal t1 [dval t1] -5.23179933356147 1.e-7 1.e-7 +checkreal t2 [dval t2] 4.76820064934972 1.e-7 1.e-7 +