]> OCCT Git - occt.git/commitdiff
0031661: Modeling Data - Exception when projecting parabola or hyperbola to plane
authorifv <ifv@opencascade.com>
Thu, 3 Mar 2022 12:05:23 +0000 (15:05 +0300)
committersmoskvin <smoskvin@opencascade.com>
Wed, 6 Apr 2022 00:15:07 +0000 (03:15 +0300)
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
src/ProjLib/ProjLib_ProjectOnPlane.hxx
tests/bugs/moddata_3/bug31661_1 [new file with mode: 0644]
tests/bugs/moddata_3/bug31661_2 [new file with mode: 0644]

index 07d32923336b8e15e4d1a6ccb549684dcfe6fbde..4330a8666c8306df76afb4c243ac2c2d4a438280 100644 (file)
 #include <GeomLib_Tool.hxx>
 #include <math_Jacobi.hxx>
 #include <math_Matrix.hxx>
+#include <gce_MakeParab.hxx>
+#include <gce_MakeDir.hxx>
+#include <LProp3d_CLProps.hxx>
+#include <math_Function.hxx>
+#include <math_BrentMinimum.hxx>
 
-
+const Standard_Real aParabolaLimit = 20000.;
+const Standard_Real aHyperbolaLimit = 10.;
 
 //=======================================================================
 //function : OnPlane_Value
 //=======================================================================
 
 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;
+    Alpha = Vector * gp_Vec(Z);
+  Alpha /= D * Z;
 
-  Vector.SetXYZ( Vector.XYZ() - Alpha * D.XYZ());
-  return Vector ;
+  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();}
-  
+  {
+    return myCurve->FirstParameter();
+  }
+
   Standard_Real LastParameter() const
-    {return myCurve->LastParameter();}
+  {
+    return myCurve->LastParameter();
+  }
+
+  Standard_Boolean Value(const Standard_Real   theT,
+    NCollection_Array1<gp_Pnt2d>& /*thePnt2d*/,
+    NCollection_Array1<gp_Pnt>&   thePnt) const
+  {
+    thePnt(1) = OnPlane_Value(theT, myCurve, myPlane, myDirection);
+    return Standard_True;
+  }
 
-    Standard_Boolean Value(const Standard_Real   theT,
-                           NCollection_Array1<gp_Pnt2d>& /*thePnt2d*/,
-                           NCollection_Array1<gp_Pnt>&   thePnt) const
-    {
-      thePnt(1) = OnPlane_Value(theT, myCurve, myPlane, myDirection);
-      return Standard_True;
-    }
-  
   Standard_Boolean D1(const Standard_Real   theT,
-                      NCollection_Array1<gp_Vec2d>& /*theVec2d*/,
-                      NCollection_Array1<gp_Vec>&   theVec) const
+    NCollection_Array1<gp_Vec2d>& /*theVec2d*/,
+    NCollection_Array1<gp_Vec>&   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);
+
+  //Try to smooth
+  Standard_Integer m1 = MaxDeg - 1;
+  for (i = 2; i < NbKnots; ++i)
+  {
+    if (BSplineCurvePtr->Multiplicity(i) == MaxDeg)
+    {
+      BSplineCurvePtr->RemoveKnot(i, m1, anErrMax);
+    }
   }
-  Mults.SetValue( NbKnots, MaxDeg + 1);
-  BSplineCurvePtr = 
-    new Geom_BSplineCurve(Poles,Knots,Mults,MaxDeg,Standard_False);
+
 }
 
 
@@ -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)
-
-      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);
-
-//  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
+  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()));
+
+    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
+    }
+    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;
-
-      Standard_Real Tol2 = myTolerance*myTolerance;
-      if (VDx.SquareMagnitude() < Tol2 ||
-          VDy.SquareMagnitude() < Tol2 ||
-          VDx.CrossSquareMagnitude(VDy) < Tol2)
+    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);
+
+      if (myKeepParam)
       {
-        myIsApprox = Standard_True;
+        myIsApprox = !gp_Dir(VDx).IsNormal(gp_Dir(VDy), Precision::Angular());
       }
-
-      if (!myIsApprox)
+      else
       {
-       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);
-
-        if (myKeepParam)
+        // 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())
         {
-          myIsApprox = !gp_Dir(VDx).IsNormal(gp_Dir(VDy), Precision::Angular());
+          // 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
         {
-          // 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)
-        {
-          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;
-          }
+          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)
+      if (!myIsApprox)
       {
-        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();
+        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;
+        }
       }
     }
-    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
-      }
+      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:
+    else if (GeomCirclePtr || GeomEllipsePtr)
     {
-      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)
+      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_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();
+    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());
 
-      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
-      }
+    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);
     }
-    break;
-  case GeomAbs_BezierCurve:
+    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()))
     {
-      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
+      myIsApprox = Standard_True;
     }
-    break ;
-  case GeomAbs_BSplineCurve:
+    else if(!myKeepParam)
     {
-      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
+      // Try building parabola with help of apex position
+      myIsApprox = !BuildParabolaByApex(GeomParabolaPtr);
     }
-    break;
-  default:
+    else
     {
-      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;
+
+    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
index d885bb5ee8831afa30c79f9b7f56e077fbf8ba2d..ec94cbb7b1516f8014476df6aac50c77bbc71e01 100644 (file)
@@ -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 (file)
index 0000000..b30ef8f
--- /dev/null
@@ -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 (file)
index 0000000..0e601af
--- /dev/null
@@ -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
+