From 27d82c622ce2bc6d2017cd43536b3e1c10da0423 Mon Sep 17 00:00:00 2001
From: ifv <ifv@opencascade.com>
Date: Fri, 14 Oct 2022 12:59:06 +0300
Subject: [PATCH] 0033170: Modeling Algorithms - Checking for canonical
 geometry: plane detection problems

GeomLib_IsPlanarSurface.cxx - using poles for checking BSpline, Bezier curves and surface changed
                              on checking by curve, surface points.

BRepOffset_MakeOffset.cxx - set normal of plane surface according to normal of initial face surface

tests/cr/bugs/bug33170 - new test case added
---
 src/BRepOffset/BRepOffset_MakeOffset.cxx |  47 ++++++++
 src/GeomLib/GeomLib_IsPlanarSurface.cxx  | 144 ++++++-----------------
 tests/cr/bugs/bug33170                   |  30 +++++
 tests/cr/grids.list                      |   1 +
 4 files changed, 117 insertions(+), 105 deletions(-)
 create mode 100644 tests/cr/bugs/bug33170

diff --git a/src/BRepOffset/BRepOffset_MakeOffset.cxx b/src/BRepOffset/BRepOffset_MakeOffset.cxx
index 7da4639d31..4579874778 100644
--- a/src/BRepOffset/BRepOffset_MakeOffset.cxx
+++ b/src/BRepOffset/BRepOffset_MakeOffset.cxx
@@ -4953,6 +4953,53 @@ Standard_Boolean BRepOffset_MakeOffset::IsPlanar()
       if (aPlanarityChecker.IsPlanar())
       {
         gp_Pln aPln = aPlanarityChecker.Plan();
+        Standard_Real u1, u2, v1, v2, um, vm;
+        aSurf->Bounds(u1, u2, v1, v2);
+        Standard_Boolean isInf1 = Precision::IsInfinite(u1), isInf2 = Precision::IsInfinite(u2);
+        if (!isInf1 && !isInf2)
+        {
+          um = (u1 + u2) / 2.;
+        }
+        else if(isInf1 && !isInf2)
+        {
+          um = u2 - 1.;
+        }
+        else if(!isInf1 && isInf2)
+        {
+          um = u1 + 1.;
+        }
+        else //isInf1 && isInf2
+        {
+          um = 0.;
+        }
+        isInf1 = Precision::IsInfinite(v1), isInf2 = Precision::IsInfinite(v2);
+        if (!isInf1 && !isInf2)
+        {
+          vm = (v1 + v2) / 2.;
+        }
+        else if (isInf1 && !isInf2)
+        {
+          vm = v2 - 1.;
+        }
+        else if(!isInf1 && isInf2)
+        {
+          vm = v1 + 1.;
+        }
+        else //isInf1 && isInf2
+        {
+          vm = 0.;
+        }
+        gp_Pnt aP;
+        gp_Vec aD1, aD2;
+        aBAS.D1(um, vm, aP, aD1, aD2);
+        gp_Vec aNorm = aD1.Crossed(aD2);
+        gp_Dir aPlnNorm = aPln.Position().Direction();
+        if (aNorm.Dot(aPlnNorm) < 0.)
+        {
+          aPlnNorm.Reverse();
+          gp_Ax1 anAx(aPln.Position().Location(), aPlnNorm);
+          aPln.SetAxis(anAx);
+        }
         Handle(Geom_Plane) aPlane = new Geom_Plane(aPln);
         TopoDS_Face aPlanarFace;
         aBB.MakeFace(aPlanarFace, aPlane, aTolForFace);
diff --git a/src/GeomLib/GeomLib_IsPlanarSurface.cxx b/src/GeomLib/GeomLib_IsPlanarSurface.cxx
index b0c04de57c..16ec51a963 100644
--- a/src/GeomLib/GeomLib_IsPlanarSurface.cxx
+++ b/src/GeomLib/GeomLib_IsPlanarSurface.cxx
@@ -30,18 +30,6 @@
 #include <TColgp_Array1OfPnt.hxx>
 #include <TColgp_HArray1OfPnt.hxx>
 
-static Standard_Boolean Controle(const TColgp_Array1OfPnt& P,
-			  const gp_Pln& Plan,
-			  const Standard_Real Tol) 
-{
-  Standard_Integer ii;
-  Standard_Boolean B=Standard_True;
-
-  for (ii=1; ii<=P.Length() && B; ii++)  
-      B = (Plan.Distance(P(ii)) < Tol);
-    
-  return B;
-}
 
 static Standard_Boolean Controle(const TColgp_Array1OfPnt& Poles,
 				 const Standard_Real Tol,
@@ -49,51 +37,36 @@ static Standard_Boolean Controle(const TColgp_Array1OfPnt& Poles,
 				 gp_Pln& Plan) 
 {
   Standard_Boolean IsPlan = Standard_False;
-  Standard_Boolean Essai = Standard_True;
   Standard_Real gx,gy,gz;
-  Standard_Integer Nb = Poles.Length();
-    gp_Pnt Bary;
+  gp_Pnt Bary;
   gp_Dir DX, DY;
+  Standard_Real aTolSingular = Precision::Confusion();
 
-  
-  if (Nb > 10) {
-    // Test allege (pour une rejection rapide)
-    TColgp_Array1OfPnt Aux(1,5);
-    Aux(1) = Poles(1);
-    Aux(2) = Poles(Nb/3);
-    Aux(3) = Poles(Nb/2);
-    Aux(4) = Poles(Nb/2+Nb/3);
-    Aux(5) =  Poles(Nb);
-    GeomLib::Inertia(Aux, Bary, DX, DY, gx, gy, gz);
-    Essai = (gz<Tol);
-  }
-  
-  if (Essai) { // Test Grandeur nature...
-    GeomLib::Inertia(Poles, Bary, DX, DY, gx, gy, gz);
-    if (gz<Tol && gy>Tol) {
-      gp_Pnt P;
-      gp_Vec DU, DV;
-      Standard_Real umin, umax, vmin, vmax;
-      S->Bounds(umin, umax, vmin, vmax);
-      S->D1( (umin+umax)/2, (vmin+vmax)/2, P, DU, DV);
-      // On prend DX le plus proche possible de DU
-      gp_Dir du(DU);
-      Standard_Real Angle1 = du.Angle(DX);
-      Standard_Real Angle2 = du.Angle(DY);
-      if (Angle1 > M_PI/2) Angle1 = M_PI-Angle1;
-      if (Angle2 > M_PI/2) Angle2 = M_PI-Angle2;
-      if (Angle2 < Angle1) {
-	du = DY; DY = DX; DX = du;
-      }
-      if (DX.Angle(DU) > M_PI/2) DX.Reverse();
-      if (DY.Angle(DV) > M_PI/2) DY.Reverse();      
-
-      gp_Ax3 axe(Bary, DX^DY, DX);
-      Plan.SetPosition(axe);
-      Plan.SetLocation(Bary);
-      IsPlan = Standard_True;
+    
+  GeomLib::Inertia(Poles, Bary, DX, DY, gx, gy, gz);
+  if (gz < Tol && gy > aTolSingular) {
+    gp_Pnt P;
+    gp_Vec DU, DV;
+    Standard_Real umin, umax, vmin, vmax;
+    S->Bounds(umin, umax, vmin, vmax);
+    S->D1((umin + umax) / 2, (vmin + vmax) / 2, P, DU, DV);
+    // On prend DX le plus proche possible de DU
+    gp_Dir du(DU);
+    Standard_Real Angle1 = du.Angle(DX);
+    Standard_Real Angle2 = du.Angle(DY);
+    if (Angle1 > M_PI / 2) Angle1 = M_PI - Angle1;
+    if (Angle2 > M_PI / 2) Angle2 = M_PI - Angle2;
+    if (Angle2 < Angle1) {
+      du = DY; DY = DX; DX = du;
     }
-  }   
+    if (DX.Angle(DU) > M_PI / 2) DX.Reverse();
+    if (DY.Angle(DV) > M_PI / 2) DY.Reverse();
+
+    gp_Ax3 axe(Bary, DX^DY, DX);
+    Plan.SetPosition(axe);
+    Plan.SetLocation(Bary);
+    IsPlan = Standard_True;
+  }
   return IsPlan;
 }
 
@@ -106,8 +79,6 @@ static Standard_Boolean Controle(const Handle(Geom_Curve)& C,
   GeomAbs_CurveType Type;
   GeomAdaptor_Curve AC(C);
   Type = AC.GetType();
-  Handle(TColgp_HArray1OfPnt) TabP;
-  TabP.Nullify();
 
   switch (Type) {
   case GeomAbs_Line :
@@ -131,40 +102,27 @@ static Standard_Boolean Controle(const Handle(Geom_Curve)& C,
   case GeomAbs_BezierCurve:
     {
       Nb =  AC.NbPoles();
-      Handle (Geom_BezierCurve) BZ = AC.Bezier();
-      TabP = new (TColgp_HArray1OfPnt) (1, AC.NbPoles());
-      for (ii=1; ii<=Nb; ii++) 
-	TabP->SetValue(ii, BZ->Pole(ii));
       break; 
     }
   case GeomAbs_BSplineCurve:
     {
       Nb =  AC.NbPoles();
-      Handle (Geom_BSplineCurve) BZ = AC.BSpline();
-      TabP = new (TColgp_HArray1OfPnt) (1, AC.NbPoles());
-      for (ii=1; ii<=Nb; ii++) 
-	TabP->SetValue(ii, BZ->Pole(ii));
       break; 
     }
-    default :
-      {
-	Nb = 8 + 3*AC.NbIntervals(GeomAbs_CN);
-      }
-  }
- 
-  if (TabP.IsNull()) {
-    Standard_Real u, du, f, l, d;
-    f = AC.FirstParameter();
-    l = AC.LastParameter();
-    du = (l-f)/(Nb-1);
-    for (ii=1; ii<=Nb && B ; ii++) {
-      u = (ii-1)*du + f;
-      d = Plan.Distance(C->Value(u));
-      B = (d < Tol); 
+  default :
+    {
+	 Nb = 8 + 3*AC.NbIntervals(GeomAbs_CN);
     }
   }
-  else {
-    B = Controle(TabP->Array1(), Plan, Tol);
+ 
+  Standard_Real u, du, f, l, d;
+  f = AC.FirstParameter();
+  l = AC.LastParameter();
+  du = (l - f) / (Nb - 1);
+  for (ii = 1; ii <= Nb && B; ii++) {
+    u = (ii - 1)*du + f;
+    d = Plan.Distance(C->Value(u));
+    B = d < Tol;
   }
 
   return B;
@@ -196,30 +154,6 @@ GeomLib_IsPlanarSurface::GeomLib_IsPlanarSurface(const Handle(Geom_Surface)& S,
      IsPlan = Standard_False;
      break;
     }
-  case GeomAbs_BezierSurface :
-  case GeomAbs_BSplineSurface :
-    {
-      Standard_Integer ii, jj, kk,
-      NbU = AS.NbUPoles(), NbV = AS.NbVPoles(); 
-      TColgp_Array1OfPnt Poles(1, NbU*NbV);
-      if (Type == GeomAbs_BezierSurface) {
-	Handle(Geom_BezierSurface) BZ;
-	BZ = AS.Bezier();
-	for(ii=1, kk=1; ii<=NbU; ii++)
-	  for(jj=1; jj<=NbV; jj++,kk++)
-	    Poles(kk) = BZ->Pole(ii,jj);
-      }
-      else {
-	Handle(Geom_BSplineSurface) BS;
-	BS = AS.BSpline();
-	for(ii=1, kk=1; ii<=NbU; ii++)
-	  for(jj=1; jj<=NbV; jj++,kk++)
-	    Poles(kk) = BS->Pole(ii,jj);
-      }
-
-      IsPlan =  Controle(Poles, Tol, S, myPlan);
-      break;      
-    }
 
   case GeomAbs_SurfaceOfRevolution :
     {
@@ -299,7 +233,7 @@ GeomLib_IsPlanarSurface::GeomLib_IsPlanarSurface(const Handle(Geom_Surface)& S,
       break;
     }
 
-    default :
+  default :
       {
 	Standard_Integer NbU,NbV, ii, jj, kk; 
 	NbU = 8 + 3*AS.NbUIntervals(GeomAbs_CN);
diff --git a/tests/cr/bugs/bug33170 b/tests/cr/bugs/bug33170
new file mode 100644
index 0000000000..9313e54b76
--- /dev/null
+++ b/tests/cr/bugs/bug33170
@@ -0,0 +1,30 @@
+puts "============"
+puts "0033170: Modeling Algorithms - Checking for canonical geometry: plane detection problems" 
+puts "============"  
+puts ""               
+
+set ExpectGap 0.0051495320504590563
+brestore [locate_data_file bug33170.brep] f
+set log1 [getanasurf asurf1 f pln  0.006]
+regexp {Gap = +([-0-9.+eE]+)} $log1 full gap1
+if {[isdraw asurf1]} {
+  set log [dump asurf1]
+  if { [regexp {Plane} $log ] != 1 } {
+     puts "Error: surface is not a plane"
+  }
+} else {
+  puts "Error: required surface is not got"
+}
+checkreal FoundGap1 $gap1 $ExpectGap 1.0e-9 0.0
+# 
+set log2 [getanasurf asurf2 f pln  1.]
+regexp {Gap = +([-0-9.+eE]+)} $log1 full gap2
+if {[isdraw asurf2]} {
+  set log [dump asurf2]
+  if { [regexp {Plane} $log ] != 1 } {
+     puts "Error: surface is not a plane"
+  }
+} else {
+  puts "Error: required surface is not got"
+}
+checkreal FoundGap2 $gap2 $ExpectGap 1.0e-9 0.0
diff --git a/tests/cr/grids.list b/tests/cr/grids.list
index 8eb5c913e5..809bfa6f44 100644
--- a/tests/cr/grids.list
+++ b/tests/cr/grids.list
@@ -1,2 +1,3 @@
 001 base
 002 approx
+003 bugs