0031016: Projection of an ellipse is a B-spline in some cases

Improve projection of ellipse and circle on a plane in case of the same parametrization of the original curve and the projected one is not necessary. Now the projection is a canonical curve instead on B-spline.

src/ProjLib/ProjLib_ProjectOnPlane.cxx

@@ -44,6 +44,9 @@
 #include <Geom_Parabola.hxx>
 #include <Geom_Hyperbola.hxx>
 #include <Geom_Ellipse.hxx>
+#include <GeomLib_Tool.hxx>
+#include <math_Jacobi.hxx>
+#include <math_Matrix.hxx>
 
 
 
@@ -516,9 +519,6 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_HCurve)&    C,
   gp_Ax2 Axis;
   Standard_Real R1 =0., R2 =0.;
 
-  if ( Type != GeomAbs_Line)  // on garde le parametrage
-    myKeepParam = Standard_True;
-  else                        // on prend le choix utilisateur.
   myKeepParam = KeepParametrization;
 
   switch ( Type) {
@@ -648,12 +648,13 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_HCurve)&    C,
 
       Standard_Real Tol2 = myTolerance*myTolerance;
       if (VDx.SquareMagnitude() < Tol2 ||
-	  VDy.SquareMagnitude() < Tol2    ) {
+          VDy.SquareMagnitude() < Tol2 ||
+          VDx.CrossSquareMagnitude(VDy) < Tol2) {
         myIsApprox = Standard_True;
       }
 
-      if (!myIsApprox && 
+      if (!myIsApprox)
-	  gp_Dir(VDx).IsNormal(gp_Dir(VDy),Precision::Angular())) {
+      {
 	Dx = gp_Dir(VDx);
 	Dy = gp_Dir(VDy);
 	gp_Pnt O    = Axis.Location();
@@ -662,6 +663,59 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_HCurve)&    C,
 	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 = !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)
+        {
           gp_Ax2 Axe(P, Dx^Dy, Dx);
 
           if (Abs(Major - Minor) < Precision::Confusion()) {
@@ -685,6 +739,13 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_HCurve)&    C,
           }
           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
@@ -692,20 +753,26 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_HCurve)&    C,
 	myResult = new GeomAdaptor_HCurve(aGACurve);
 //  Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
       }
-      }
+      else if (GeomCirclePtr || GeomEllipsePtr)
-      else {
+      {
-	myIsApprox = Standard_True;
+        Handle(Geom_Curve) aResultCurve = GeomCirclePtr;
-	myType     = GeomAbs_BSplineCurve;
+        if (aResultCurve.IsNull())
-        PerformApprox(myCurve,myPlane,myDirection,ApproxCurve);
+          aResultCurve = GeomEllipsePtr;
-//  Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
+        // start and end parameters of the projected curve
-	GeomAdaptor_Curve aGACurve(ApproxCurve);
+        Standard_Real aParFirst = myCurve->FirstParameter();
-	myResult = new GeomAdaptor_HCurve(aGACurve);
+        Standard_Real aParLast = myCurve->LastParameter();
-//  Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
+        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:
     {
+      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)
@@ -757,6 +824,7 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_HCurve)&    C,
     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)
@@ -824,6 +892,7 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_HCurve)&    C,
       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++) {
@@ -847,6 +916,7 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_HCurve)&    C,
       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++) {
@@ -862,6 +932,7 @@ void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_HCurve)&    C,
     break;
   default:
     {
+      myKeepParam = Standard_True;
       myIsApprox = Standard_True;
       myType     = GeomAbs_BSplineCurve;
       PerformApprox(myCurve,myPlane,myDirection,ApproxCurve);

tests/bugs/modalg_7/bug31016_01

@@ -0,0 +1,53 @@
+puts "========"
+puts "OCC31016"
+puts "========"
+puts ""
+###################################################################
+# Projection of an ellipse or a circle is a B-spline in some cases
+###################################################################
+
+ellipse c 0 0 0  0 0 1  2 1 0  20 10
+plane p 0 0 0  0 1 10
+projonplane r c p 0
+
+if {![regexp {Ellipse} [dump r]]} {
+  puts "ERROR: Projected curve is not an ellipse"
+}
+
+proc toPlane {pln curve param} {
+  upvar $pln p
+  upvar $curve crv
+  cvalue crv $param x y z
+  set pr [proj p x y z]
+  regexp {Parameters: ([-0-9.+eE]+) ([-0-9.+eE]+)} $pr full u v
+  svalue p $u $v x y z
+
+  set pntOnPlane {}
+  lappend pntOnPlane [dval x]
+  lappend pntOnPlane [dval y]
+  lappend pntOnPlane [dval z]
+  return $pntOnPlane
+}
+
+# calculate a parametric shift on the projected curve
+set pnt [toPlane p c 0]
+parameters r [lindex $pnt 0] [lindex $pnt 1] [lindex $pnt 2] 0.1 shift
+
+set nbPoints 100
+for {set i 0} {$i <= $nbPoints} {incr i} {
+  set par0 [expr 8*atan(1)*$i/$nbPoints]
+  set par1 [expr $par0 + [dval shift]]
+
+  set pnt [toPlane p c $par0]
+  cvalue r $par1 X Y Z
+
+  set dx [expr [lindex $pnt 0]-[dval X]]
+  set dy [expr [lindex $pnt 1]-[dval Y]]
+  set dz [expr [lindex $pnt 2]-[dval Z]]
+
+  if {[expr $dx*$dx + $dy*$dy + $dz*$dz] < 1.e-14} {
+    puts "OK: Projection correct"
+  } else {
+    puts "ERROR: Projection incorrect"
+  }
+}

tests/bugs/modalg_7/bug31016_02

@@ -0,0 +1,55 @@
+puts "========"
+puts "OCC31016"
+puts "========"
+puts ""
+###################################################################
+# Projection of an ellipse or a circle is a B-spline in some cases
+###################################################################
+
+set projDir { 1 1 1 }
+
+ellipse c 0 0 0  0 0 1  2 1 0  20 10
+plane p 0 0 0  0 1 10
+projonplane r c p [lindex $projDir 0] [lindex $projDir 1] [lindex $projDir 2] 0
+
+if {![regexp {Ellipse} [dump r]]} {
+  puts "ERROR: Projected curve is not an ellipse"
+}
+
+proc toPlane {pln dir curve param} {
+  upvar $pln p
+  upvar $curve crv
+  cvalue crv $param x y z
+  line ln x y z [lindex $dir 0] [lindex $dir 1] [lindex $dir 2]
+  intersect pt ln p
+  regexp {Point : ([-0-9.+eE]+), ([-0-9.+eE]+), ([-0-9.+eE]+)} [dump pt] full x y z
+
+  set pntOnPlane {}
+  lappend pntOnPlane $x
+  lappend pntOnPlane $y
+  lappend pntOnPlane $z
+  return $pntOnPlane
+}
+
+# calculate a parametric shift on the projected curve
+set pnt [toPlane p $projDir c 0]
+parameters r [lindex $pnt 0] [lindex $pnt 1] [lindex $pnt 2] 0.1 shift
+
+set nbPoints 100
+for {set i 0} {$i <= $nbPoints} {incr i} {
+  set par0 [expr 8*atan(1)*$i/$nbPoints]
+  set par1 [expr $par0 + [dval shift]]
+
+  set pnt [toPlane p $projDir c $par0]
+  cvalue r $par1 X Y Z
+
+  set dx [expr [lindex $pnt 0]-[dval X]]
+  set dy [expr [lindex $pnt 1]-[dval Y]]
+  set dz [expr [lindex $pnt 2]-[dval Z]]
+
+  if {[expr $dx*$dx + $dy*$dy + $dz*$dz] < 1.e-14} {
+    puts "OK: Projection correct"
+  } else {
+    puts "ERROR: Projection incorrect"
+  }
+}

tests/bugs/modalg_7/bug31016_03

@@ -0,0 +1,53 @@
+puts "========"
+puts "OCC31016"
+puts "========"
+puts ""
+###################################################################
+# Projection of an ellipse or a circle is a B-spline in some cases
+###################################################################
+
+circle c 0 0 0  10
+plane p 0 0 0  1 0 1
+projonplane r c p 0
+
+if {![regexp {Ellipse} [dump r]]} {
+  puts "ERROR: Projected curve is not an ellipse"
+}
+
+proc toPlane {pln curve param} {
+  upvar $pln p
+  upvar $curve crv
+  cvalue crv $param x y z
+  set pr [proj p x y z]
+  regexp {Parameters: ([-0-9.+eE]+) ([-0-9.+eE]+)} $pr full u v
+  svalue p $u $v x y z
+
+  set pntOnPlane {}
+  lappend pntOnPlane [dval x]
+  lappend pntOnPlane [dval y]
+  lappend pntOnPlane [dval z]
+  return $pntOnPlane
+}
+
+# calculate a parametric shift on the projected curve
+set pnt [toPlane p c 0]
+parameters r [lindex $pnt 0] [lindex $pnt 1] [lindex $pnt 2] 0.1 shift
+
+set nbPoints 100
+for {set i 0} {$i <= $nbPoints} {incr i} {
+  set par0 [expr 8*atan(1)*$i/$nbPoints]
+  set par1 [expr $par0 + [dval shift]]
+
+  set pnt [toPlane p c $par0]
+  cvalue r $par1 X Y Z
+
+  set dx [expr [lindex $pnt 0]-[dval X]]
+  set dy [expr [lindex $pnt 1]-[dval Y]]
+  set dz [expr [lindex $pnt 2]-[dval Z]]
+
+  if {[expr $dx*$dx + $dy*$dy + $dz*$dz] < 1.e-14} {
+    puts "OK: Projection correct"
+  } else {
+    puts "ERROR: Projection incorrect"
+  }
+}

tests/bugs/modalg_7/bug31016_04

@@ -0,0 +1,55 @@
+puts "========"
+puts "OCC31016"
+puts "========"
+puts ""
+###################################################################
+# Projection of an ellipse or a circle is a B-spline in some cases
+###################################################################
+
+set projDir { 1 1 1 }
+
+circle c 0 0 0  10
+plane p 0 0 0  1 0 1
+projonplane r c p [lindex $projDir 0] [lindex $projDir 1] [lindex $projDir 2] 0
+
+if {![regexp {Ellipse} [dump r]]} {
+  puts "ERROR: Projected curve is not an ellipse"
+}
+
+proc toPlane {pln dir curve param} {
+  upvar $pln p
+  upvar $curve crv
+  cvalue crv $param x y z
+  line ln x y z [lindex $dir 0] [lindex $dir 1] [lindex $dir 2]
+  intersect pt ln p
+  regexp {Point : ([-0-9.+eE]+), ([-0-9.+eE]+), ([-0-9.+eE]+)} [dump pt] full x y z
+
+  set pntOnPlane {}
+  lappend pntOnPlane $x
+  lappend pntOnPlane $y
+  lappend pntOnPlane $z
+  return $pntOnPlane
+}
+
+# calculate a parametric shift on the projected curve
+set pnt [toPlane p $projDir c 0]
+parameters r [lindex $pnt 0] [lindex $pnt 1] [lindex $pnt 2] 0.1 shift
+
+set nbPoints 100
+for {set i 0} {$i <= $nbPoints} {incr i} {
+  set par0 [expr 8*atan(1)*$i/$nbPoints]
+  set par1 [expr $par0 + [dval shift]]
+
+  set pnt [toPlane p $projDir c $par0]
+  cvalue r $par1 X Y Z
+
+  set dx [expr [lindex $pnt 0]-[dval X]]
+  set dy [expr [lindex $pnt 1]-[dval Y]]
+  set dz [expr [lindex $pnt 2]-[dval Z]]
+
+  if {[expr $dx*$dx + $dy*$dy + $dz*$dz] < 1.e-14} {
+    puts "OK: Projection correct"
+  } else {
+    puts "ERROR: Projection incorrect"
+  }
+}

tests/bugs/modalg_7/bug31016_05

@@ -0,0 +1,53 @@
+puts "========"
+puts "OCC31016"
+puts "========"
+puts ""
+###################################################################
+# Projection of an ellipse or a circle is a B-spline in some cases
+###################################################################
+
+ellipse c 0 0 0  0.866025403784439 0 0.5  0.5 0 -0.866025403784439  20 10
+plane p 0 0 0  0 0 1
+projonplane r c p 0
+
+if {![regexp {Circle} [dump r]]} {
+  puts "ERROR: Projected curve is not a circle"
+}
+
+proc toPlane {pln curve param} {
+  upvar $pln p
+  upvar $curve crv
+  cvalue crv $param x y z
+  set pr [proj p x y z]
+  regexp {Parameters: ([-0-9.+eE]+) ([-0-9.+eE]+)} $pr full u v
+  svalue p $u $v x y z
+
+  set pntOnPlane {}
+  lappend pntOnPlane [dval x]
+  lappend pntOnPlane [dval y]
+  lappend pntOnPlane [dval z]
+  return $pntOnPlane
+}
+
+# calculate a parametric shift on the projected curve
+set pnt [toPlane p c 0]
+parameters r [lindex $pnt 0] [lindex $pnt 1] [lindex $pnt 2] 0.1 shift
+
+set nbPoints 100
+for {set i 0} {$i <= $nbPoints} {incr i} {
+  set par0 [expr 8*atan(1)*$i/$nbPoints]
+  set par1 [expr $par0 + [dval shift]]
+
+  set pnt [toPlane p c $par0]
+  cvalue r $par1 X Y Z
+
+  set dx [expr [lindex $pnt 0]-[dval X]]
+  set dy [expr [lindex $pnt 1]-[dval Y]]
+  set dz [expr [lindex $pnt 2]-[dval Z]]
+
+  if {[expr $dx*$dx + $dy*$dy + $dz*$dz] < 1.e-14} {
+    puts "OK: Projection correct"
+  } else {
+    puts "ERROR: Projection incorrect"
+  }
+}

tests/bugs/modalg_7/bug31016_06

@@ -0,0 +1,55 @@
+puts "========"
+puts "OCC31016"
+puts "========"
+puts ""
+###################################################################
+# Projection of an ellipse or a circle is a B-spline in some cases
+###################################################################
+
+set projDir { 1 1 1 }
+
+ellipse c 0 0 0  0 0 1  2 1 0  20 10
+plane p 0 0 0  0 0 1
+projonplane r c p [lindex $projDir 0] [lindex $projDir 1] [lindex $projDir 2] 0
+
+if {![regexp {Ellipse} [dump r]]} {
+  puts "ERROR: Projected curve is not an ellipse"
+}
+
+proc toPlane {pln dir curve param} {
+  upvar $pln p
+  upvar $curve crv
+  cvalue crv $param x y z
+  line ln x y z [lindex $dir 0] [lindex $dir 1] [lindex $dir 2]
+  intersect pt ln p
+  regexp {Point : ([-0-9.+eE]+), ([-0-9.+eE]+), ([-0-9.+eE]+)} [dump pt] full x y z
+
+  set pntOnPlane {}
+  lappend pntOnPlane $x
+  lappend pntOnPlane $y
+  lappend pntOnPlane $z
+  return $pntOnPlane
+}
+
+# calculate a parametric shift on the projected curve
+set pnt [toPlane p $projDir c 0]
+parameters r [lindex $pnt 0] [lindex $pnt 1] [lindex $pnt 2] 0.1 shift
+
+set nbPoints 100
+for {set i 0} {$i <= $nbPoints} {incr i} {
+  set par0 [expr 8*atan(1)*$i/$nbPoints]
+  set par1 [expr $par0 + [dval shift]]
+
+  set pnt [toPlane p $projDir c $par0]
+  cvalue r $par1 X Y Z
+
+  set dx [expr [lindex $pnt 0]-[dval X]]
+  set dy [expr [lindex $pnt 1]-[dval Y]]
+  set dz [expr [lindex $pnt 2]-[dval Z]]
+
+  if {[expr $dx*$dx + $dy*$dy + $dz*$dz] < 1.e-14} {
+    puts "OK: Projection correct"
+  } else {
+    puts "ERROR: Projection incorrect"
+  }
+}

tests/bugs/modalg_7/bug31016_07

@@ -0,0 +1,51 @@
+puts "========"
+puts "OCC31016"
+puts "========"
+puts ""
+###################################################################
+# Projection of an ellipse or a circle is a B-spline in some cases
+###################################################################
+
+ellipse c 0 0 0  0 0 1  2 1 0  20 10
+plane p 0 0 0  0 1 10
+projonplane r c p 1
+
+if {![regexp {BSplineCurve} [dump r]]} {
+  puts "ERROR: Projected curve is not a B-spline curve"
+}
+
+proc toPlane {pln curve param} {
+  upvar $pln p
+  upvar $curve crv
+  cvalue crv $param x y z
+  set pr [proj p x y z]
+  regexp {Parameters: ([-0-9.+eE]+) ([-0-9.+eE]+)} $pr full u v
+  svalue p $u $v x y z
+
+  set pntOnPlane {}
+  lappend pntOnPlane [dval x]
+  lappend pntOnPlane [dval y]
+  lappend pntOnPlane [dval z]
+  return $pntOnPlane
+}
+
+# calculate a parametric shift on the projected curve
+set pnt [toPlane p c 0]
+
+set nbPoints 100
+for {set i 0} {$i <= $nbPoints} {incr i} {
+  set par0 [expr 8*atan(1)*$i/$nbPoints]
+
+  set pnt [toPlane p c $par0]
+  cvalue r $par0 X Y Z
+
+  set dx [expr [lindex $pnt 0]-[dval X]]
+  set dy [expr [lindex $pnt 1]-[dval Y]]
+  set dz [expr [lindex $pnt 2]-[dval Z]]
+
+  if {[expr $dx*$dx + $dy*$dy + $dz*$dz] < 1.e-12} {
+    puts "OK: Projection correct"
+  } else {
+    puts "ERROR: Projection incorrect"
+  }
+}

tests/bugs/modalg_7/bug31016_08

@@ -0,0 +1,53 @@
+puts "========"
+puts "OCC31016"
+puts "========"
+puts ""
+###################################################################
+# Projection of an ellipse or a circle is a B-spline in some cases
+###################################################################
+
+set projDir { 1 1 1 }
+
+ellipse c 0 0 0  0 0 1  2 1 0  20 10
+plane p 0 0 0  0 1 10
+projonplane r c p [lindex $projDir 0] [lindex $projDir 1] [lindex $projDir 2] 1
+
+if {![regexp {BSplineCurve} [dump r]]} {
+  puts "ERROR: Projected curve is not a B-spline curve"
+}
+
+proc toPlane {pln dir curve param} {
+  upvar $pln p
+  upvar $curve crv
+  cvalue crv $param x y z
+  line ln x y z [lindex $dir 0] [lindex $dir 1] [lindex $dir 2]
+  intersect pt ln p
+  regexp {Point : ([-0-9.+eE]+), ([-0-9.+eE]+), ([-0-9.+eE]+)} [dump pt] full x y z
+
+  set pntOnPlane {}
+  lappend pntOnPlane $x
+  lappend pntOnPlane $y
+  lappend pntOnPlane $z
+  return $pntOnPlane
+}
+
+# calculate a parametric shift on the projected curve
+set pnt [toPlane p $projDir c 0]
+
+set nbPoints 100
+for {set i 0} {$i <= $nbPoints} {incr i} {
+  set par0 [expr 8*atan(1)*$i/$nbPoints]
+
+  set pnt [toPlane p $projDir c $par0]
+  cvalue r $par0 X Y Z
+
+  set dx [expr [lindex $pnt 0]-[dval X]]
+  set dy [expr [lindex $pnt 1]-[dval Y]]
+  set dz [expr [lindex $pnt 2]-[dval Z]]
+
+  if {[expr $dx*$dx + $dy*$dy + $dz*$dz] < 1.e-12} {
+    puts "OK: Projection correct"
+  } else {
+    puts "ERROR: Projection incorrect"
+  }
+}

tests/bugs/modalg_7/bug31016_09

@@ -0,0 +1,51 @@
+puts "========"
+puts "OCC31016"
+puts "========"
+puts ""
+###################################################################
+# Projection of an ellipse or a circle is a B-spline in some cases
+###################################################################
+
+circle c 0 0 0  10
+plane p 0 0 0  1 0 1
+projonplane r c p 1
+
+if {![regexp {BSplineCurve} [dump r]]} {
+  puts "ERROR: Projected curve is not a B-spline curve"
+}
+
+proc toPlane {pln curve param} {
+  upvar $pln p
+  upvar $curve crv
+  cvalue crv $param x y z
+  set pr [proj p x y z]
+  regexp {Parameters: ([-0-9.+eE]+) ([-0-9.+eE]+)} $pr full u v
+  svalue p $u $v x y z
+
+  set pntOnPlane {}
+  lappend pntOnPlane [dval x]
+  lappend pntOnPlane [dval y]
+  lappend pntOnPlane [dval z]
+  return $pntOnPlane
+}
+
+# calculate a parametric shift on the projected curve
+set pnt [toPlane p c 0]
+
+set nbPoints 100
+for {set i 0} {$i <= $nbPoints} {incr i} {
+  set par0 [expr 8*atan(1)*$i/$nbPoints]
+
+  set pnt [toPlane p c $par0]
+  cvalue r $par0 X Y Z
+
+  set dx [expr [lindex $pnt 0]-[dval X]]
+  set dy [expr [lindex $pnt 1]-[dval Y]]
+  set dz [expr [lindex $pnt 2]-[dval Z]]
+
+  if {[expr $dx*$dx + $dy*$dy + $dz*$dz] < 1.e-12} {
+    puts "OK: Projection correct"
+  } else {
+    puts "ERROR: Projection incorrect"
+  }
+}

tests/bugs/modalg_7/bug31016_10

@@ -0,0 +1,53 @@
+puts "========"
+puts "OCC31016"
+puts "========"
+puts ""
+###################################################################
+# Projection of an ellipse or a circle is a B-spline in some cases
+###################################################################
+
+set projDir { 1 1 1 }
+
+circle c 0 0 0  10
+plane p 0 0 0  1 0 1
+projonplane r c p [lindex $projDir 0] [lindex $projDir 1] [lindex $projDir 2] 1
+
+if {![regexp {BSplineCurve} [dump r]]} {
+  puts "ERROR: Projected curve is not a B-spline curve"
+}
+
+proc toPlane {pln dir curve param} {
+  upvar $pln p
+  upvar $curve crv
+  cvalue crv $param x y z
+  line ln x y z [lindex $dir 0] [lindex $dir 1] [lindex $dir 2]
+  intersect pt ln p
+  regexp {Point : ([-0-9.+eE]+), ([-0-9.+eE]+), ([-0-9.+eE]+)} [dump pt] full x y z
+
+  set pntOnPlane {}
+  lappend pntOnPlane $x
+  lappend pntOnPlane $y
+  lappend pntOnPlane $z
+  return $pntOnPlane
+}
+
+# calculate a parametric shift on the projected curve
+set pnt [toPlane p $projDir c 0]
+
+set nbPoints 100
+for {set i 0} {$i <= $nbPoints} {incr i} {
+  set par0 [expr 8*atan(1)*$i/$nbPoints]
+
+  set pnt [toPlane p $projDir c $par0]
+  cvalue r $par0 X Y Z
+
+  set dx [expr [lindex $pnt 0]-[dval X]]
+  set dy [expr [lindex $pnt 1]-[dval Y]]
+  set dz [expr [lindex $pnt 2]-[dval Z]]
+
+  if {[expr $dx*$dx + $dy*$dy + $dz*$dz] < 1.e-12} {
+    puts "OK: Projection correct"
+  } else {
+    puts "ERROR: Projection incorrect"
+  }
+}

tests/bugs/modalg_7/bug31016_11

@@ -0,0 +1,51 @@
+puts "========"
+puts "OCC31016"
+puts "========"
+puts ""
+###################################################################
+# Projection of an ellipse or a circle is a B-spline in some cases
+###################################################################
+
+ellipse c 0 0 0  0.866025403784439 0 0.5  0.5 0 -0.866025403784439  20 10
+plane p 0 0 0  0 0 1
+projonplane r c p 1
+
+if {![regexp {Circle} [dump r]]} {
+  puts "ERROR: Projected curve is not a circle"
+}
+
+proc toPlane {pln curve param} {
+  upvar $pln p
+  upvar $curve crv
+  cvalue crv $param x y z
+  set pr [proj p x y z]
+  regexp {Parameters: ([-0-9.+eE]+) ([-0-9.+eE]+)} $pr full u v
+  svalue p $u $v x y z
+
+  set pntOnPlane {}
+  lappend pntOnPlane [dval x]
+  lappend pntOnPlane [dval y]
+  lappend pntOnPlane [dval z]
+  return $pntOnPlane
+}
+
+# calculate a parametric shift on the projected curve
+set pnt [toPlane p c 0]
+
+set nbPoints 100
+for {set i 0} {$i <= $nbPoints} {incr i} {
+  set par0 [expr 8*atan(1)*$i/$nbPoints]
+
+  set pnt [toPlane p c $par0]
+  cvalue r $par0 X Y Z
+
+  set dx [expr [lindex $pnt 0]-[dval X]]
+  set dy [expr [lindex $pnt 1]-[dval Y]]
+  set dz [expr [lindex $pnt 2]-[dval Z]]
+
+  if {[expr $dx*$dx + $dy*$dy + $dz*$dz] < 1.e-12} {
+    puts "OK: Projection correct"
+  } else {
+    puts "ERROR: Projection incorrect"
+  }
+}

tests/bugs/modalg_7/bug31016_12

@@ -0,0 +1,53 @@
+puts "========"
+puts "OCC31016"
+puts "========"
+puts ""
+###################################################################
+# Projection of an ellipse or a circle is a B-spline in some cases
+###################################################################
+
+set projDir { 1 1 1 }
+
+ellipse c 0 0 0  0 0 1  2 1 0  20 10
+plane p 0 0 0  0 0 1
+projonplane r c p [lindex $projDir 0] [lindex $projDir 1] [lindex $projDir 2] 1
+
+if {![regexp {Ellipse} [dump r]]} {
+  puts "ERROR: Projected curve is not an ellipse"
+}
+
+proc toPlane {pln dir curve param} {
+  upvar $pln p
+  upvar $curve crv
+  cvalue crv $param x y z
+  line ln x y z [lindex $dir 0] [lindex $dir 1] [lindex $dir 2]
+  intersect pt ln p
+  regexp {Point : ([-0-9.+eE]+), ([-0-9.+eE]+), ([-0-9.+eE]+)} [dump pt] full x y z
+
+  set pntOnPlane {}
+  lappend pntOnPlane $x
+  lappend pntOnPlane $y
+  lappend pntOnPlane $z
+  return $pntOnPlane
+}
+
+# calculate a parametric shift on the projected curve
+set pnt [toPlane p $projDir c 0]
+
+set nbPoints 100
+for {set i 0} {$i <= $nbPoints} {incr i} {
+  set par0 [expr 8*atan(1)*$i/$nbPoints]
+
+  set pnt [toPlane p $projDir c $par0]
+  cvalue r $par0 X Y Z
+
+  set dx [expr [lindex $pnt 0]-[dval X]]
+  set dy [expr [lindex $pnt 1]-[dval Y]]
+  set dz [expr [lindex $pnt 2]-[dval Z]]
+
+  if {[expr $dx*$dx + $dy*$dy + $dz*$dz] < 1.e-14} {
+    puts "OK: Projection correct"
+  } else {
+    puts "ERROR: Projection incorrect"
+  }
+}