0028196: Modeling Data - Algorithm 'Extrema_GenLocateExtPS' failed to find the extremum in a case

Extrema_GenLocateExtPS.cxx, Extrema_GenLocateExtPS.hxx:

Adjusting tolerances according to surface sizes is added.
Additional methods for searching solution are added for cases if basic method fails.

Extrema_FuncPSDist.cxx - small bug fixing.

BRepFill_TrimShellCorner.cxx - fixing regression

ChFi3d_Builder_CnCrn.cxx setting parameters for Plate algorithm to improve stability of solution

ProjLib_ComputeApproxOnPolarSurface.cxx - code optimization

Some test cases are modified according to current state of Extrema algorithm

src/BRepFill/BRepFill_TrimShellCorner.cxx

@@ -281,11 +281,30 @@ void BRepFill_TrimShellCorner::Perform()
     }
   }
   
+  Standard_Real aMaxTol = 0.;
+  TopExp_Explorer anExp(myShape1, TopAbs_VERTEX);
+  for (; anExp.More(); anExp.Next())
+  {
+    aMaxTol = Max(aMaxTol, BRep_Tool::Tolerance(TopoDS::Vertex(anExp.Current())));
+  }
+
+  anExp.Init(myShape2, TopAbs_VERTEX);
+  for (; anExp.More(); anExp.Next())
+  {
+    aMaxTol = Max(aMaxTol, BRep_Tool::Tolerance(TopoDS::Vertex(anExp.Current())));
+  }
+
+  Standard_Real aFuzzy = 4.*Precision::Confusion();
   BOPAlgo_PaveFiller aPF;
   TopTools_ListOfShape aLS;
   aLS.Append(myShape1);
   aLS.Append(myShape2);
   aPF.SetArguments(aLS);
+  if (aMaxTol < 1.005 * Precision::Confusion())
+  {
+    aFuzzy = Max(aPF.FuzzyValue(), aFuzzy);
+    aPF.SetFuzzyValue(aFuzzy);
+  }
   //
   aPF.Perform();
   if (aPF.HasErrors()) {
@@ -856,6 +875,9 @@ Standard_Boolean BRepFill_TrimShellCorner::ChooseSection(const TopoDS_Shape& Com
     TopoDS_Edge FirstEdge = FindEdgeCloseToBisectorPlane(theFirstVertex,
                                                           OldComp,
                                                           myAxeOfBisPlane.Axis());
+    if (FirstEdge.IsNull())
+      return Standard_False;
+
     iter.Initialize(OldComp);
     if (!iter.More())
     {
@@ -865,6 +887,8 @@ Standard_Boolean BRepFill_TrimShellCorner::ChooseSection(const TopoDS_Shape& Com
     TopoDS_Edge LastEdge  = FindEdgeCloseToBisectorPlane(theLastVertex,
                                                           OldComp,
                                                           myAxeOfBisPlane.Axis());
+    if (LastEdge.IsNull())
+      return Standard_False;
 
     BB.Add(NewWire, FirstEdge);
 

src/ChFi3d/ChFi3d_Builder_CnCrn.cxx

@@ -2062,8 +2062,18 @@ void  ChFi3d_Builder::PerformMoreThreeCorner(const Standard_Integer Jndex,
   } 
 
 // declaration for plate 
-  GeomPlate_BuildPlateSurface PSurf(3,10,3,tol2d,tolesp,angular);
+  //GeomPlate_BuildPlateSurface PSurf(3,10,3,tol2d,tolesp,angular);
-
+  //
+  //Sence of Plate parameters and their preferable values :
+  // degree is total order of ordinary or mixed derivatives:
+  // dS/dU, dS/dV have degree 1, d2S/dU2, d2S/dV2, d2S/(dUdV) have degree 2
+  // nbiter - number of iterations, when surface from previous iteration uses as initial surface for next one
+  // practically this process does not converge, using "bad" initial surface leads to much more "bad" solution.
+  // constr is order of constraint: 0 - G0, 1 - G1 ...
+  // Using constraint order > 0 very often causes  unpredicable undulations of solution
+  Standard_Integer degree = 3, nbcurvpnt = 10, nbiter = 1;
+  Standard_Integer constr = 1; //G1
+  GeomPlate_BuildPlateSurface PSurf(degree, nbcurvpnt, nbiter, tol2d, tolesp, angular);
 // calculation of curves on surface for each stripe 
   for (ic=0;ic<nedge;ic++) {
     gp_Pnt2d p2d1, p2d2;
@@ -2088,9 +2098,10 @@ void  ChFi3d_Builder::PerformMoreThreeCorner(const Standard_Integer Jndex,
       Adaptor3d_CurveOnSurface  CurvOnS (Acurv,Asurf);
       Handle(Adaptor3d_CurveOnSurface) HCons =
 	new Adaptor3d_CurveOnSurface(CurvOnS);
-      Order.SetValue(ic,1);
+      //Order.SetValue(ic,1);
+      Order.SetValue(ic, constr);
       Handle(GeomPlate_CurveConstraint) Cont =
-	new GeomPlate_CurveConstraint(HCons,Order.Value(ic),10,tolesp,angular,0.1);
+	new GeomPlate_CurveConstraint(HCons,Order.Value(ic), nbcurvpnt,tolesp,angular,0.1);
       PSurf.Add(Cont);
       
       // calculate indexes of points and of the curve for the DS         

src/Extrema/Extrema_FuncPSDist.cxx

@@ -107,8 +107,8 @@ Standard_Boolean Extrema_FuncPSDist::IsInside(const math_Vector& X)
 {
     if (X(1) < mySurf.FirstUParameter() ||
         X(1) > mySurf.LastUParameter() ||
-        X(2) < mySurf.FirstUParameter() ||
+        X(2) < mySurf.FirstVParameter() ||
-        X(2) > mySurf.LastUParameter() )
+        X(2) > mySurf.LastVParameter() )
     {
       // Point out of borders.
       return Standard_False;

src/Extrema/Extrema_GenLocateExtPS.cxx

@@ -23,10 +23,84 @@
 #include <Extrema_POnSurf.hxx>
 #include <gp_Pnt.hxx>
 #include <math_FunctionSetRoot.hxx>
+#include <math_NewtonFunctionSetRoot.hxx>
 #include <math_BFGS.hxx>
-#include <math_Vector.hxx>
+#include <math_BFGS.hxx>
+#include <math_FRPR.hxx>
 #include <StdFail_NotDone.hxx>
 
+static void CorrectTol(const Standard_Real theU0, const Standard_Real theV0,
+  math_Vector& theTol)
+{
+  //Correct tolerance for large values of UV parameters
+  Standard_Real aTolRef = Precision::PConfusion();
+  Standard_Real anEpsRef = Epsilon(1.);
+  Standard_Real epsu = Epsilon(theU0);
+  const Standard_Real tolog10 = 0.43429;
+  if (epsu > anEpsRef)
+  {
+    Standard_Integer n = RealToInt(tolog10 * Log(epsu / anEpsRef) + 1) + 1;
+    Standard_Integer i;
+    Standard_Real tol = aTolRef;
+    for (i = 1; i <= n; ++i)
+    {
+      tol *= 10.;
+    }
+    theTol(1) = Max(theTol(1), tol);
+  }
+  Standard_Real epsv = Epsilon(theV0);
+  if (epsv > anEpsRef)
+  {
+    Standard_Integer n = RealToInt(tolog10 * Log(epsv / anEpsRef) + 1) + 1;
+    Standard_Integer i;
+    Standard_Real tol = aTolRef;
+    for (i = 1; i <= n; ++i)
+    {
+      tol *= 10.;
+    }
+    theTol(2) = Max(theTol(2), tol);
+  }
+}
+//=======================================================================
+//function : IsMinDist
+//purpose  : 
+//=======================================================================
+
+Standard_Boolean Extrema_GenLocateExtPS::IsMinDist(const gp_Pnt& theP, const Adaptor3d_Surface& theS,
+  const Standard_Real theU0, const Standard_Real theV0)
+{
+  Standard_Real du = Max(theS.UResolution(10.*Precision::Confusion()), 10.*Precision::PConfusion());
+  Standard_Real dv = Max(theS.VResolution(10.*Precision::Confusion()), 10.*Precision::PConfusion());
+  Standard_Real u, v;
+  gp_Pnt aP0 = theS.Value(theU0, theV0);
+  Standard_Real d0 = theP.SquareDistance(aP0);
+  Standard_Integer iu, iv;
+  for (iu = -1; iu <= 1; ++iu)
+  {
+    u = theU0 + iu * du;
+    if (!theS.IsUPeriodic())
+    {
+      u = Max(u, theS.FirstUParameter());
+      u = Min(u, theS.LastUParameter());
+    }
+    for (iv = -1; iv <= 1; ++iv)
+    {
+      if (iu == 0 && iv == 0)
+        continue;
+
+      v = theV0 + iv * dv;
+      if (!theS.IsVPeriodic())
+      {
+        v = Max(v, theS.FirstVParameter());
+        v = Min(v, theS.LastVParameter());
+      }
+      Standard_Real d = theP.SquareDistance(theS.Value(u, v));
+      if (d < d0)
+        return Standard_False;
+    }
+  }
+  return Standard_True;
+}
 //=======================================================================
 //function : Extrema_GenLocateExtPS
 //purpose  : 
@@ -69,36 +143,84 @@ void Extrema_GenLocateExtPS::Perform(const gp_Pnt& theP,
   aBoundSup(1) = mySurf.LastUParameter();
   aBoundSup(2) = mySurf.LastVParameter();
 
-  if (isDistanceCriteria == Standard_False)
+  if (isDistanceCriteria)
   {
-    // Normal projection criteria.
+    // Distance criteria.
-    Extrema_FuncPSNorm F(theP,mySurf);
+    Standard_Real aRelTol = 1.e-8;
+    math_Vector aResPnt(1, 2);  
 
-    math_FunctionSetRoot SR (F, aTol);
+    Extrema_FuncPSDist F(mySurf, theP);
-    SR.Perform(F, aStart, aBoundInf, aBoundSup);
+
-    if (!SR.IsDone()) 
+    math_BFGS aSolver(2, aRelTol);
+    aSolver.Perform(F, aStart);
+
+    if (!aSolver.IsDone())
+    {
+      //Try another method
+      math_FRPR aSolver1(F, aRelTol);
+      aSolver1.Perform(F, aStart);
+      if(!aSolver1.IsDone())
         return;
+      aSolver1.Location(aResPnt);
+      mySqDist = aSolver1.Minimum();
+    }
+    else
+    {
+      aSolver.Location(aResPnt);    
+      mySqDist = aSolver.Minimum();
+    }
 
-    mySqDist = F.SquareDistance(1);
+    myPoint.SetParameters(aResPnt(1), aResPnt(2), mySurf.Value(aResPnt(1), aResPnt(2)));
-    myPoint = F.Point(1);
     myDone = Standard_True;
   }
   else
   {
-    // Distance criteria.
+    // Normal projection criteria.
-    Extrema_FuncPSDist F(mySurf, theP);
+    Extrema_FuncPSNorm F(theP, mySurf);
-    math_BFGS aSolver(2);
-    aSolver.Perform(F, aStart);
 
-    if (!aSolver.IsDone()) 
+    if (mySurf.GetType() == GeomAbs_BSplineSurface)
+    {
+      aTol(1) = myTolU;
+      aTol(2) = myTolV;
+      CorrectTol(theU0, theV0, aTol);
+    }
+
+    Standard_Boolean isCorrectTol = (Abs(aTol(1) - myTolU) > Precision::PConfusion() ||
+      Abs(aTol(2) - myTolV) > Precision::PConfusion());
+
+    math_FunctionSetRoot aSR(F, aTol);
+    aSR.Perform(F, aStart, aBoundInf, aBoundSup);
+
+    if (!aSR.IsDone() || isCorrectTol)
+    {
+      if (isCorrectTol)
+      {
+        aTol(1) = myTolU;
+        aTol(2) = myTolV;
+      }
+      math_NewtonFunctionSetRoot aNSR(F, aTol, Precision::Confusion());
+      aNSR.Perform(F, aStart, aBoundInf, aBoundSup);
+      if (!aSR.IsDone() && !aNSR.IsDone())
+      {
         return;
+      }
+    }
     
-    math_Vector aResPnt(1,2);
+    Standard_Real aNbExt = F.NbExt();
-    aSolver.Location(aResPnt);
+    mySqDist = F.SquareDistance(1);
-    mySqDist = aSolver.Minimum();
+    myPoint = F.Point(1);
-    myPoint.SetParameters(aResPnt(1), aResPnt(2), mySurf.Value(aResPnt(1), aResPnt(2)));
+    Standard_Integer i;
+    for (i = 2; i <= aNbExt; ++i)
+    {
+      if (F.SquareDistance(i) < mySqDist)
+      {
+        mySqDist = F.SquareDistance(i);
+        myPoint = F.Point(i);
+      }
+    }
     myDone = Standard_True;
   }
+
 }
 
 //=======================================================================

src/Extrema/Extrema_GenLocateExtPS.hxx

@@ -63,6 +63,11 @@ public:
   //! Returns the point of the extremum distance.
   Standard_EXPORT const Extrema_POnSurf& Point() const;
 
+  //! Returns True if UV point theU0, theV0 is point of local minimum of square distance between
+  //! point theP and points theS(U, V), U, V are in small area around theU0, theV0
+  Standard_EXPORT static Standard_Boolean IsMinDist(const gp_Pnt& theP, const Adaptor3d_Surface& theS,
+    const Standard_Real theU0, const Standard_Real theV0);
+
 private:
 
   const Extrema_GenLocateExtPS& operator=(const Extrema_GenLocateExtPS&);

src/GeometryTest/GeometryTest_APICommands.cxx

@@ -136,7 +136,7 @@ static Standard_Integer proj (Draw_Interpretor& di, Standard_Integer n, const ch
       const gp_XY aP2d(Draw::Atof(a[5]), Draw::Atof(a[6]));
       GeomAdaptor_Surface aGAS(GS);
       Extrema_GenLocateExtPS aProjector(aGAS, Precision::PConfusion(), Precision::PConfusion());
-      aProjector.Perform(P, aP2d.X(), aP2d.Y());
+      aProjector.Perform(P, aP2d.X(), aP2d.Y(), Standard_False);
       if (!aProjector.IsDone())
       {
         di << "projection failed.";

src/ProjLib/ProjLib_ComputeApproxOnPolarSurface.cxx

@@ -1097,7 +1097,7 @@ Handle(Adaptor2d_Curve2d)
     TColgp_SequenceOfPnt2d Sols;
     Standard_Boolean areManyZeros = Standard_False;
     
-    Curve->D0(Param.Value(1), pntproj) ;
+    pntproj = Pts(1);
     Extrema_ExtPS  aExtPS(pntproj, *Surf, TolU, TolV) ;
     Standard_Real aMinSqDist = RealLast();
     if (aExtPS.IsDone())
@@ -1172,7 +1172,7 @@ Handle(Adaptor2d_Curve2d)
 	      Standard_Integer nbExtOk = 0;
 	      Standard_Integer indExt = 0;
 	      Standard_Integer iT = 1 + (NbOfPnts - 1)/5*i;
-	      Curve->D0( Param.Value(iT), pntproj );
+              pntproj = Pts(iT);
 	      Extrema_ExtPS aTPS( pntproj, *Surf, TolU, TolV );
 	      Dist2Min = 1.e+200;
 	      if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
@@ -1199,7 +1199,7 @@ Handle(Adaptor2d_Curve2d)
 	      Standard_Boolean isFound = Standard_False;
               for (j = tPp + 1; j <= NbOfPnts; ++j)
               {
-		Curve->D0( Param.Value(j), pntproj );
+                pntproj = Pts(j);
 		Extrema_ExtPS aTPS( pntproj, *Surf, TolU, TolV );
 		Dist2Min = RealLast();
 		if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
@@ -1220,14 +1220,16 @@ Handle(Adaptor2d_Curve2d)
 	      }
 
 	      if( isFound ) {
-		gp_Vec2d atV(aPp,aPn);
+		gp_Dir2d atV(gp_Vec2d(aPp,aPn));
 		Standard_Boolean isChosen = Standard_False;
 		for( i = 1; i <= nbSols; i++ ) {
 		  const gp_Pnt2d& aP1 = Sols.Value(i);
-		  gp_Vec2d asV(aP1,aPp);
+		  gp_Dir2d asV(gp_Vec2d(aP1,aPp));
 		  if( asV.Dot(atV) > 0. ) {
 		    isChosen = Standard_True;
-		    Pts2d(1).SetCoord(aP1.X(),aP1.Y());
+                    u = aP1.X();
+                    v = aP1.Y();
+		    Pts2d(1).SetCoord(u, v);
 		    myProjIsDone = Standard_True;
 		    break;
 		  }
@@ -1271,7 +1273,7 @@ Handle(Adaptor2d_Curve2d)
 	if(myProjIsDone) {
 	  myProjIsDone = Standard_False;
 	  Dist2Min = RealLast();
-	  Curve->D0(Param.Value(i), pntproj);
+          pntproj = Pts(i);
           Extrema_GenLocateExtPS  aLocateExtPS (*Surf, TolU, TolV);
           aLocateExtPS.Perform(pntproj, U0, V0);
 

tests/blend/complex/A3

@@ -9,4 +9,4 @@ tscale s 0 0 0 1000
 explode s e
 blend result s 5 s_8
 
-checkprops result -s 3824.84
+checkprops result -s 3740

tests/blend/simple/W4

@@ -8,4 +8,4 @@ restore [locate_data_file CCV_1_c12gsf.rle] s
 explode s E
 blend result s 0.2 s_2 0.2 s_1
 
-checkprops result -s 11.5154
+checkprops result -s 11.5416

tests/bugs/moddata_3/bug28196

@@ -1,12 +1,7 @@
 puts "============"
-puts "CR28196"
+puts "OCC28196: Modeling Data - Algorithm 'Extrema_GenLocateExtPS' finds wrong extremum in a case"
 puts "==========="
 puts ""
-###############################################################################
-# Algorithm 'Extrema_GenLocateExtPS' failed to find the extremum in a case
-###############################################################################
-
-puts "TODO OCC28196 ALL: Error: projection with starting point is worse"
 
 pload MODELING
 
@@ -18,7 +13,7 @@ dset y 0.0187851531995338
 dset z 0.433545417254429
 dset u 13140.3030987283
 dset v 28.5494495724281
-set max_error 1e-5
+set max_tol 1e-5
 
 proj s x y z
 set dist_proj [lindex [length ext_1] end]
@@ -28,7 +23,7 @@ set dist_proj_uv [lindex [length ext_1] end]
 puts "Distance of standard projection $dist_proj"
 puts "Distance of projection with starting point $dist_proj_uv"
 
-if {[expr $dist_proj_uv - $dist_proj] > $max_error} {
+if {[expr $dist_proj_uv - $dist_proj] > $max_tol} {
   puts "Error: projection with starting point is worse than standard projection"
 } else {
   puts "OK: well projection with starting point"

tests/perf/modalg/bug453_2

@@ -22,6 +22,6 @@ tcopy result_1 result
 
 dchrono h2 stop counter blend
 
-checkprops result -s 3.65777e+06 
+checkprops result -s 3.5496e+06 
 checkshape result
 checkview -display result -2d -path ${imagedir}/${test_image}.png

tests/perf/moddata/bug453_3

@@ -1,7 +1,7 @@
 puts "TODO OCC24156 MacOS: Tcl Exception: tolerance ang"
 puts "TODO OCC24156 MacOS: TEST INCOMPLETE"
 puts "TODO OCC27203 ALL: Error: Max tolerance" 
-puts "TODO OCC27203 All: Error : The area of result shape is"
+##puts "TODO OCC27203 All: Error : The area of result shape is"
 
 puts "========"
 puts "OCC453"