From 43dbdb15db3e510d0167b8c1fa43eeac67491fcc Mon Sep 17 00:00:00 2001
From: ifv <ifv@opencascade.com>
Date: Wed, 5 Nov 2014 16:44:37 +0300
Subject: [PATCH] 0025427: Algorithm of building plane from wire hangs

Test-case for issue #25427
---
 src/BRepLib/BRepLib_FindSurface.cxx | 291 +++++++++++++++-------------
 tests/bugs/modalg_5/bug25427        |  17 ++
 2 files changed, 173 insertions(+), 135 deletions(-)
 create mode 100644 tests/bugs/modalg_5/bug25427

diff --git a/src/BRepLib/BRepLib_FindSurface.cxx b/src/BRepLib/BRepLib_FindSurface.cxx
index e010750065..5363865588 100644
--- a/src/BRepLib/BRepLib_FindSurface.cxx
+++ b/src/BRepLib/BRepLib_FindSurface.cxx
@@ -58,7 +58,7 @@
 //purpose  : 
 //=======================================================================
 static Standard_Real Controle(const TColgp_SequenceOfPnt& thePoints,
-			      const Handle(Geom_Plane)& thePlane)
+  const Handle(Geom_Plane)& thePlane)
 {
   Standard_Real dfMaxDist=0.;
   Standard_Real a,b,c,d, dist;
@@ -79,9 +79,9 @@ static Standard_Real Controle(const TColgp_SequenceOfPnt& thePoints,
 //           the first vertex of theEdge2 in parametric space of theFace
 //=======================================================================
 inline static Standard_Boolean Is2DConnected(const TopoDS_Edge& theEdge1,
-					     const TopoDS_Edge& theEdge2,
-					     const Handle(Geom_Surface)& theSurface,
-					     const TopLoc_Location& theLocation)
+  const TopoDS_Edge& theEdge2,
+  const Handle(Geom_Surface)& theSurface,
+  const TopLoc_Location& theLocation)
 {
   Standard_Real f,l;
   //TopLoc_Location aLoc;
@@ -110,8 +110,8 @@ inline static Standard_Boolean Is2DConnected(const TopoDS_Edge& theEdge1,
 //=======================================================================
 
 static Standard_Boolean Is2DClosed(const TopoDS_Shape&         theShape,
-                                   const Handle(Geom_Surface)& theSurface,
-				   const TopLoc_Location& theLocation)
+  const Handle(Geom_Surface)& theSurface,
+  const TopLoc_Location& theLocation)
 {
   try
   {
@@ -164,9 +164,9 @@ BRepLib_FindSurface::BRepLib_FindSurface()
 //purpose  : 
 //=======================================================================
 BRepLib_FindSurface::BRepLib_FindSurface(const TopoDS_Shape&    S, 
-					 const Standard_Real    Tol,
-					 const Standard_Boolean OnlyPlane,
-                                         const Standard_Boolean OnlyClosed)
+  const Standard_Real    Tol,
+  const Standard_Boolean OnlyPlane,
+  const Standard_Boolean OnlyClosed)
 {
   Init(S,Tol,OnlyPlane,OnlyClosed);
 }
@@ -175,9 +175,9 @@ BRepLib_FindSurface::BRepLib_FindSurface(const TopoDS_Shape&    S,
 //purpose  : 
 //=======================================================================
 void BRepLib_FindSurface::Init(const TopoDS_Shape&    S, 
-			       const Standard_Real    Tol,
-			       const Standard_Boolean OnlyPlane,
-                               const Standard_Boolean OnlyClosed)
+  const Standard_Real    Tol,
+  const Standard_Boolean OnlyPlane,
+  const Standard_Boolean OnlyClosed)
 {
   myTolerance = Tol;
   myTolReached = 0.;
@@ -214,24 +214,24 @@ void BRepLib_FindSurface::Init(const TopoDS_Shape&    S,
     // check the other edges
     for (ex.Init(S,TopAbs_EDGE); ex.More(); ex.Next()) {
       if (!E.IsSame(ex.Current())) {
-	j = 0;
-	for(;;) {
-	  j++;
-	  BRep_Tool::CurveOnSurface(TopoDS::Edge(ex.Current()),
-				    PPC,SS,L,ff,ll,j);
-	  if (SS.IsNull()) {
-	    break;
-	  }
-	  if (SS == mySurface) {
-	    break;
-	  }
-	  SS.Nullify();
-	}
+        j = 0;
+        for(;;) {
+          j++;
+          BRep_Tool::CurveOnSurface(TopoDS::Edge(ex.Current()),
+            PPC,SS,L,ff,ll,j);
+          if (SS.IsNull()) {
+            break;
+          }
+          if (SS == mySurface) {
+            break;
+          }
+          SS.Nullify();
+        }
 
-	if (SS.IsNull()) {
-	  mySurface.Nullify();
-	  break;
-	}
+        if (SS.IsNull()) {
+          mySurface.Nullify();
+          break;
+        }
       }
     }
 
@@ -262,7 +262,7 @@ void BRepLib_FindSurface::Init(const TopoDS_Shape&    S,
   //                       distances from neighboring points (_only_ same edge)
   //                    2. Minimizing the weighed sum of squared deviations
   //                       compute coefficients of the sought plane.
-  
+
   TColgp_SequenceOfPnt aPoints;
   TColStd_SequenceOfReal aWeight;
 
@@ -284,64 +284,84 @@ void BRepLib_FindSurface::Init(const TopoDS_Shape&    S,
     {
     case GeomAbs_BezierCurve:
       {
-	// Put all poles for bezier
-	Handle(Geom_BezierCurve) GC = c.Bezier();
-	Standard_Integer iNbPol = GC->NbPoles();
-	if ( iNbPol < 2)
-	  // Degenerate
-	  continue;
-	else
-	{
-	  Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
-	  GC->Poles(aPoles->ChangeArray1());
-	  gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
-	  Standard_Real dfDistPrev = 0., dfDistNext;
-	  for (Standard_Integer iPol=1; iPol<=iNbPol; iPol++)
-	  {
-	    if (iPol<iNbPol)
-	    {
-	      aPoleNext = aPoles->Value(iPol+1);
-	      dfDistNext = aPolePrev.Distance(aPoleNext);
-	    }
-	    else
-	      dfDistNext = 0.;
-	    aPoints.Append (aPolePrev);
-	    aWeight.Append (dfDistPrev+dfDistNext);
-	    dfDistPrev = dfDistNext;
-	    aPolePrev = aPoleNext;
-	  }
-	}
+        // Put all poles for bezier
+        Handle(Geom_BezierCurve) GC = c.Bezier();
+        Standard_Integer iNbPol = GC->NbPoles();
+        Standard_Real tf = GC->FirstParameter();
+        Standard_Real tl = GC->LastParameter();
+        Standard_Real r =  (dfUl - dfUf) / (tl - tf);
+        r *= iNbPol;
+        if ( iNbPol < 2 || r < 1.)
+          // Degenerate
+          continue;
+        else
+        {
+          Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
+          GC->Poles(aPoles->ChangeArray1());
+          gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
+          Standard_Real dfDistPrev = 0., dfDistNext;
+          for (Standard_Integer iPol=1; iPol<=iNbPol; iPol++)
+          {
+            if (iPol<iNbPol)
+            {
+              aPoleNext = aPoles->Value(iPol+1);
+              dfDistNext = aPolePrev.Distance(aPoleNext);
+            }
+            else
+              dfDistNext = 0.;
+            aPoints.Append (aPolePrev);
+            aWeight.Append (dfDistPrev+dfDistNext);
+            dfDistPrev = dfDistNext;
+            aPolePrev = aPoleNext;
+          }
+        }
       }
       break;
     case GeomAbs_BSplineCurve:
       {
-	// Put all poles for bspline
-	Handle(Geom_BSplineCurve) GC = c.BSpline();
-	Standard_Integer iNbPol = GC->NbPoles();
-	if ( iNbPol < 2)
-	  // Degenerate
-	  continue;
-	else
-	{
-	  Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
-	  GC->Poles(aPoles->ChangeArray1());
-	  gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
-	  Standard_Real dfDistPrev = 0., dfDistNext;
-	  for (Standard_Integer iPol=1; iPol<=iNbPol; iPol++)
-	  {
-	    if (iPol<iNbPol)
-	    {
-	      aPoleNext = aPoles->Value(iPol+1);
-	      dfDistNext = aPolePrev.Distance(aPoleNext);
-	    }
-	    else
-	      dfDistNext = 0.;
-	    aPoints.Append (aPolePrev);
-	    aWeight.Append (dfDistPrev+dfDistNext);
-	    dfDistPrev = dfDistNext;
-	    aPolePrev = aPoleNext;
-	  }
-	}
+        // Put all poles for bspline
+        Handle(Geom_BSplineCurve) GC = c.BSpline();
+        Standard_Integer iNbPol = GC->NbPoles();
+        Standard_Real tf = GC->FirstParameter();
+        Standard_Real tl = GC->LastParameter();
+        Standard_Real r =  (dfUl - dfUf) / (tl - tf);
+        r *= iNbPol;
+        if ( iNbPol < 2 || r < 1.)
+          // Degenerate
+          continue;
+        else
+        {
+          const Standard_Integer aNbPolMax = 200;
+          Standard_Integer incr = 1;
+          if(iNbPol > aNbPolMax)
+          {
+            Standard_Integer nb = iNbPol;
+            while(nb > aNbPolMax)
+            {
+              incr++;
+              nb = (iNbPol-1) / incr;
+            }
+          }
+          Handle(TColgp_HArray1OfPnt) aPoles = new (TColgp_HArray1OfPnt) (1, iNbPol);
+          GC->Poles(aPoles->ChangeArray1());
+          gp_Pnt aPolePrev = aPoles->Value(1), aPoleNext;
+          Standard_Real dfDistPrev = 0., dfDistNext;
+          Standard_Integer iPol;
+          for (iPol = 1; iPol <= iNbPol; iPol += incr)
+          {
+            if (iPol <= iNbPol - incr)
+            {
+              aPoleNext = aPoles->Value(iPol+incr);
+              dfDistNext = aPolePrev.Distance(aPoleNext);
+            }
+            else
+              dfDistNext = 0.;
+            aPoints.Append (aPolePrev);
+            aWeight.Append (dfDistPrev+dfDistNext);
+            dfDistPrev = dfDistNext;
+            aPolePrev = aPoleNext;
+          }
+        }
       }
       break;
 
@@ -351,41 +371,41 @@ void BRepLib_FindSurface::Init(const TopoDS_Shape&    S,
     case GeomAbs_Hyperbola:
     case GeomAbs_Parabola:
       if (c.GetType() == GeomAbs_Line)
-	// Two points on straight segment
-	iNbPoints=2;
+        // Two points on straight segment
+        iNbPoints=2;
       else
-	// Four points on otheranalitical curves
-	iNbPoints=4;
+        // Four points on otheranalitical curves
+        iNbPoints=4;
     default:
       { 
-	// Put some points on other curves
-	if (iNbPoints==0)
-	  iNbPoints = 15 + c.NbIntervals(GeomAbs_C3);
-	Standard_Real dfDelta = (dfUl-dfUf)/(iNbPoints-1);
-	Standard_Integer iPoint;
-	Standard_Real dfU;
-	gp_Pnt aPointPrev = c.Value(dfUf), aPointNext;
-	Standard_Real dfDistPrev = 0., dfDistNext;
-	for (iPoint=1, dfU=dfUf+dfDelta; 
-	     iPoint<=iNbPoints; 
-	     iPoint++, dfU+=dfDelta) 
-	{
-	  if (iPoint<iNbPoints)
-	  {
-	    aPointNext = c.Value(dfU);
-	    dfDistNext = aPointPrev.Distance(aPointNext);
-	  }
-	  else
-	    dfDistNext = 0.;
-	  aPoints.Append (aPointPrev);
-	  aWeight.Append (dfDistPrev+dfDistNext);
-	  dfDistPrev = dfDistNext;
-	  aPointPrev = aPointNext;
-	}
+        // Put some points on other curves
+        if (iNbPoints==0)
+          iNbPoints = 15 + c.NbIntervals(GeomAbs_C3);
+        Standard_Real dfDelta = (dfUl-dfUf)/(iNbPoints-1);
+        Standard_Integer iPoint;
+        Standard_Real dfU;
+        gp_Pnt aPointPrev = c.Value(dfUf), aPointNext;
+        Standard_Real dfDistPrev = 0., dfDistNext;
+        for (iPoint=1, dfU=dfUf+dfDelta; 
+          iPoint<=iNbPoints; 
+          iPoint++, dfU+=dfDelta) 
+        {
+          if (iPoint<iNbPoints)
+          {
+            aPointNext = c.Value(dfU);
+            dfDistNext = aPointPrev.Distance(aPointNext);
+          }
+          else
+            dfDistNext = 0.;
+          aPoints.Append (aPointPrev);
+          aWeight.Append (dfDistPrev+dfDistNext);
+          dfDistPrev = dfDistNext;
+          aPointPrev = aPointNext;
+        }
       } // default:
     } // switch (c.GetType()) ...
   } // for (ex.Init(S,TopAbs_EDGE); ex.More() && control; ex.Next()) ...
-  
+
   if (aPoints.Length() < 3) {
     return;
   }
@@ -414,14 +434,14 @@ void BRepLib_FindSurface::Init(const TopoDS_Shape&    S,
     gp_XYZ p=aPoints(iPoint).XYZ()-aBaryCenter;
     Standard_Real w=aWeight(iPoint)/dfMaxWeight;
     aMat(1,1)+=w*p.X()*p.X(); 
-        aMat(1,2)+=w*p.X()*p.Y(); 
-            aMat(1,3)+=w*p.X()*p.Z();
+    aMat(1,2)+=w*p.X()*p.Y(); 
+    aMat(1,3)+=w*p.X()*p.Z();
     aMat(2,1)+=w*p.Y()*p.X();  
-        aMat(2,2)+=w*p.Y()*p.Y();  
-            aMat(2,3)+=w*p.Y()*p.Z();
+    aMat(2,2)+=w*p.Y()*p.Y();  
+    aMat(2,3)+=w*p.Y()*p.Z();
     aMat(3,1)+=w*p.Z()*p.X();  
-        aMat(3,2)+=w*p.Z()*p.Y(); 
-            aMat(3,3)+=w*p.Z()*p.Z();
+    aMat(3,2)+=w*p.Z()*p.Y(); 
+    aMat(3,3)+=w*p.Z()*p.Z();
     aVec(1) -= w*p.X();
     aVec(2) -= w*p.Y();
     aVec(3) -= w*p.Z();
@@ -450,26 +470,27 @@ void BRepLib_FindSurface::Init(const TopoDS_Shape&    S,
   if (!isSolved || myTolerance < dfDist)  {
     gp_Pnt aFirstPnt=aPoints(1);
     for (iPoint=2; iPoint<=aPoints.Length(); iPoint++)  {
-      gp_Vec aDir(aFirstPnt,aPoints(iPoint));
+      const gp_Pnt& aNextPnt = aPoints(iPoint); 
+      gp_Vec aDir(aFirstPnt, aNextPnt);
       Standard_Real dfSide=aDir.Magnitude();
       if (dfSide<myTolerance) {
-	continue; // degeneration
+        continue; // degeneration
       }
       for (Standard_Integer iP1=iPoint+1; iP1<=aPoints.Length(); iP1++) {
-	gp_Vec aCross = gp_Vec(aFirstPnt,aPoints(iP1)) ^ aDir ;
-	if (aCross.Magnitude() > dfSide*myTolerance) {
-	  Handle(Geom_Plane) aPlane2 = new Geom_Plane(aFirstPnt, aCross);
-	  Standard_Real dfDist2 = Controle (aPoints, aPlane2);
-	  if (dfDist2 < myTolerance)  {
-	    myTolReached = dfDist2;
-	    mySurface = aPlane2;
-	    return;
-	  }
-	  if (dfDist2 < dfDist)  {
-	    dfDist = dfDist2;
-	    aPlane = aPlane2;
-	  }
-	}
+        gp_Vec aCross = gp_Vec(aFirstPnt,aPoints(iP1)) ^ aDir ;
+        if (aCross.Magnitude() > dfSide*myTolerance) {
+          Handle(Geom_Plane) aPlane2 = new Geom_Plane(aFirstPnt, aCross);
+          Standard_Real dfDist2 = Controle (aPoints, aPlane2);
+          if (dfDist2 < myTolerance)  {
+            myTolReached = dfDist2;
+            mySurface = aPlane2;
+            return;
+          }
+          if (dfDist2 < dfDist)  {
+            dfDist = dfDist2;
+            aPlane = aPlane2;
+          }
+        }
       }
     }
   }
diff --git a/tests/bugs/modalg_5/bug25427 b/tests/bugs/modalg_5/bug25427
new file mode 100644
index 0000000000..4272b49b2f
--- /dev/null
+++ b/tests/bugs/modalg_5/bug25427
@@ -0,0 +1,17 @@
+puts "========"
+puts "OCC25427"
+puts "========"
+puts ""
+###############################################
+# Algorithm of building plane from wire hangs
+###############################################
+
+smallview
+
+restore [locate_data_file bug25427_w.brep] w
+
+mkplane p w 1
+
+fit
+
+set only_screen_axo 1