0026075: Make Extrema_GenExtCC return IsParallel flag in case of parallel curves

1) Added check for parallel curves.
2) Changed unefficient o(n^2) duplicates deleting algorithm to o(n) algorithm.
3) Deleted useless upper level duplicates deleting algorithm.

Test-case for issue #26075

src/Extrema/Extrema_ExtCC.cxx

@@ -678,42 +678,37 @@ void Extrema_ExtCC::Results(const Extrema_ECC&   AlgExt,
 			     const Standard_Real  Ut21,
 			     const Standard_Real  Ut22)
 {
-  Standard_Integer i, j,NbExt;
+  Standard_Integer i, NbExt;
-  Standard_Real Val, U, U2,Uj,U2j;
+  Standard_Real Val, U, U2;
-  Extrema_POnCurv P1, P2,P1j,P2j;
+  Extrema_POnCurv P1, P2;
-  Standard_Boolean IsExtrema;
 
   myDone = AlgExt.IsDone();
-  if (myDone) {
+  if (myDone)
+  {
+    myIsPar = AlgExt.IsParallel();
     NbExt = AlgExt.NbExt();
-    for (i = 1; i <= NbExt; i++) {
+    for (i = 1; i <= NbExt; i++)
+    {
       AlgExt.Points(i, P1, P2);
       U = P1.Parameter();
       U2 = P2.Parameter();
-      IsExtrema=Standard_True;
-      for  (j=1;j<=mynbext;j++)
-      { P1j=mypoints.Value(2*j-1);
-        P2j=mypoints.Value(2*j);
-        Uj=P1j.Parameter();
-        U2j=P2j.Parameter();
-        if ((Abs(Uj-U)<=myTol[0]) && (Abs(U2j-U2)<=myTol[1]))     
-        IsExtrema=Standard_False;}
 
-      if (IsExtrema) 
+      // Check points to be into param space.
+      if (Extrema_CurveTool::IsPeriodic(*((Adaptor3d_Curve*)myC[0])))
       {
-//  Verification de la validite des parametres
-	if (Extrema_CurveTool::IsPeriodic(*((Adaptor3d_Curve*)myC[0]))) {
         U = ElCLib::InPeriod(U, Ut11, Ut11+Extrema_CurveTool::Period(*((Adaptor3d_Curve*)myC[0])));
       }
-	if (Extrema_CurveTool::IsPeriodic(*((Adaptor3d_Curve*)myC[1]))) {
+      if (Extrema_CurveTool::IsPeriodic(*((Adaptor3d_Curve*)myC[1])))
+      {
         U2 = ElCLib::InPeriod(U2, Ut21, Ut21+Extrema_CurveTool::Period(*((Adaptor3d_Curve*)myC[1])));
       }
 
       if ((U  >= Ut11 - RealEpsilon())  &&
           (U  <= Ut12 + RealEpsilon())  &&
           (U2 >= Ut21 - RealEpsilon())  &&
-	   (U2 <= Ut22 + RealEpsilon()))   
+          (U2 <= Ut22 + RealEpsilon())   )
-         { mynbext++;
+      {
+        mynbext++;
         Val = AlgExt.SquareDistance(i);
         mySqDist.Append(Val);
         P1.SetValues(U, P1.Value());
@@ -723,5 +718,4 @@ void Extrema_ExtCC::Results(const Extrema_ECC&   AlgExt,
       }
     }
   }
-  }  
 }

src/Extrema/Extrema_ExtCC2d.cxx

@@ -512,6 +512,7 @@ void Extrema_ExtCC2d::Results(const Extrema_ECC2d& AlgExt,
   myDone = AlgExt.IsDone();
   if (myDone)
   {
+    myIsPar = AlgExt.IsParallel();
     if (!myIsPar)
     {
       NbExt = AlgExt.NbExt();

src/Extrema/Extrema_GenExtCC.cdl

@@ -66,6 +66,10 @@ is
     	---Purpose: Returns True if the distances are found.
     	is static;
 
+    IsParallel (me) returns Boolean
+      ---Purpose: Returns state of myParallel flag.
+      is static;
+
     NbExt (me) returns Integer
     	---Purpose: Returns the number of extremum distances.
     	raises  NotDone from StdFail,
@@ -96,6 +100,7 @@ is
     	is static;
 
 fields
+  myParallel : Boolean;
   myCurveMinTol : Real from Standard;
   myLowBorder : Vector from math;
   myUppBorder : Vector from math;

src/Extrema/Extrema_GenExtCC.gxx

@@ -14,6 +14,8 @@
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
+#include <algorithm>
+
 #include <Extrema_GlobOptFuncCC.hxx>
 #include <math_GlobOptMin.hxx>
 #include <Standard_NullObject.hxx>
@@ -21,13 +23,82 @@
 #include <StdFail_NotDone.hxx>
 #include <TColStd_Array1OfReal.hxx>
 #include <Precision.hxx>
+#include <NCollection_Vector.hxx>
+#include <NCollection_CellFilter.hxx>
+
+// Comparator, used in std::sort.
+static Standard_Boolean comp(const gp_XY& theA,
+                             const gp_XY& theB)
+{
+  if (theA.X() < theB.X())
+  {
+    return Standard_True;
+  }
+  else
+  {
+    if (theA.X() == theB.X())
+    {
+      if (theA.Y() <= theB.Y())
+        return Standard_True;
+    }
+  }
+
+  return Standard_False;
+}
+
+class Extrema_CCPointsInspector : public NCollection_CellFilter_InspectorXY
+{
+public:
+  typedef gp_XY Target;
+  //! Constructor; remembers the tolerance
+  Extrema_CCPointsInspector (const Standard_Real theTol)
+  {
+    myTol = theTol * theTol;
+    myIsFind = Standard_False;
+  }
+
+  void ClearFind()
+  {
+    myIsFind = Standard_False;
+  }
+
+  Standard_Boolean isFind()
+  {
+    return myIsFind;
+  }
+
+  //! Set current point to search for coincidence
+  void SetCurrent (const gp_XY& theCurPnt) 
+  { 
+    myCurrent = theCurPnt;
+  }
+
+  //! Implementation of inspection method
+  NCollection_CellFilter_Action Inspect (const Target& theObject)
+  {
+    gp_XY aPt = myCurrent.Subtracted(theObject);
+    const Standard_Real aSQDist = aPt.SquareModulus();
+    if(aSQDist < myTol)
+    {
+      myIsFind = Standard_True;
+    }
+
+    return CellFilter_Keep;
+  }
+
+private:
+  Standard_Real myTol;
+  gp_XY myCurrent;
+  Standard_Boolean myIsFind;
+};
 
 //=======================================================================
 //function : Extrema_GenExtCC
 //purpose  : 
 //=======================================================================
 Extrema_GenExtCC::Extrema_GenExtCC()
-: myCurveMinTol(Precision::PConfusion()),
+: myParallel(Standard_False),
+  myCurveMinTol(Precision::PConfusion()),
   myLowBorder(1,2),
   myUppBorder(1,2),
   myDone(Standard_False)
@@ -41,7 +112,8 @@ Extrema_GenExtCC::Extrema_GenExtCC()
 //=======================================================================
 Extrema_GenExtCC::Extrema_GenExtCC(const Curve1& C1,
                                    const Curve2& C2)
-: myCurveMinTol(Precision::PConfusion()),
+: myParallel(Standard_False),
+  myCurveMinTol(Precision::PConfusion()),
   myLowBorder(1,2),
   myUppBorder(1,2),
   myDone(Standard_False)
@@ -64,7 +136,8 @@ Extrema_GenExtCC::Extrema_GenExtCC(const Curve1& C1,
                                    const Standard_Real Usup,
                                    const Standard_Real Vinf,
                                    const Standard_Real Vsup)
-: myCurveMinTol(Precision::PConfusion()),
+: myParallel(Standard_False),
+  myCurveMinTol(Precision::PConfusion()),
   myLowBorder(1,2),
   myUppBorder(1,2),
   myDone(Standard_False)
@@ -112,6 +185,7 @@ void Extrema_GenExtCC::SetTolerance(Standard_Real theTol)
 void Extrema_GenExtCC::Perform()
 {
   myDone = Standard_False;
+  myParallel = Standard_False;
 
   Curve1 &C1 = *(Curve1*)myC[0];
   Curve2 &C2 = *(Curve2*)myC[1];
@@ -131,14 +205,19 @@ void Extrema_GenExtCC::Perform()
   Standard_Real aSameTol = myCurveMinTol / (aDiscTol);
   aFinder.SetTol(aDiscTol, aSameTol);
 
-  Standard_Real anEps1 = (myUppBorder(1) - myLowBorder(1)) * Precision::Confusion();
+  // Size computed to have cell index inside of int32 value.
-  Standard_Real anEps2 = (myUppBorder(2) - myLowBorder(2)) * Precision::Confusion();
+  const Standard_Real aCellSize = Max(anIntervals1.Upper() - anIntervals1.Lower(),
+                                      anIntervals2.Upper() - anIntervals2.Lower())
+                                  * Precision::PConfusion() / (2.0 * Sqrt(2.0));
+  Extrema_CCPointsInspector anInspector(Precision::PConfusion());
+  NCollection_CellFilter<Extrema_CCPointsInspector> aFilter(aCellSize);
+  NCollection_Vector<gp_XY> aPnts;
 
   Standard_Integer i,j,k;
   math_Vector aFirstBorderInterval(1,2);
   math_Vector aSecondBorderInterval(1,2);
-  Standard_Real aF = RealLast(); // best functional value
+  Standard_Real aF = RealLast(); // Best functional value.
-  Standard_Real aCurrF = RealLast(); // current functional value computed on current interval
+  Standard_Real aCurrF = RealLast(); // Current functional value computed on current interval.
   for(i = 1; i <= aNbInter[0]; i++)
   {
     for(j = 1; j <= aNbInter[1]; j++)
@@ -151,14 +230,14 @@ void Extrema_GenExtCC::Perform()
       aFinder.SetLocalParams(aFirstBorderInterval, aSecondBorderInterval);
       aFinder.Perform();
 
-      // check that solution found on current interval is not worse than previous
+      // Check that solution found on current interval is not worse than previous.
       aCurrF = aFinder.GetF();
       if (aCurrF >= aF + aSameTol * aValueTol)
       {
         continue;
       }
 
-      // clean previously computed solution if current one is better
+      // Clean previously computed solution if current one is better.
       if (aCurrF > aF - aSameTol * aValueTol)
       {
         if (aCurrF < aF)
@@ -167,31 +246,88 @@ void Extrema_GenExtCC::Perform()
       else
       {
         aF = aCurrF;
-        myPoints1.Clear();
+        aFilter.Reset(aCellSize);
-        myPoints2.Clear();
+        aPnts.Clear();
       }
 
-      // save found solutions avoiding repetitions
+      // Save found solutions avoiding repetitions.
       math_Vector sol(1,2);
       for(k = 1; k <= aFinder.NbExtrema(); k++)
       {
         aFinder.Points(k, sol);
+        gp_XY aPnt2d(sol(1), sol(2));
 
-        // avoid duplicated points
+        gp_XY aXYmin = anInspector.Shift(aPnt2d, -aCellSize);
-        Standard_Boolean isNew = Standard_True;
+        gp_XY aXYmax = anInspector.Shift(aPnt2d,  aCellSize);
-        for (Standard_Integer iSol = 1; isNew && iSol <= myPoints1.Length(); iSol++)
+
+        anInspector.ClearFind();
+        anInspector.SetCurrent(aPnt2d);
+        aFilter.Inspect(aXYmin, aXYmax, anInspector);
+        if (!anInspector.isFind())
         {
-          if (Abs(myPoints1(iSol) - sol(1)) < anEps1 &&
+          // Point is out of close cells, add new one.
-              Abs(myPoints2(iSol) - sol(2)) < anEps2)
+          aFilter.Add(aPnt2d, aPnt2d);
-            isNew = Standard_False;
+          aPnts.Append(gp_XY(sol(1), sol(2)));
         }
-        if (isNew)
+      }
+    }
+  }
+
+  if (aPnts.Size() == 0)
   {
-          myPoints1.Append(sol(1));
+    // No solutions.
-          myPoints2.Append(sol(2));
+    myDone = Standard_False;
+    return;
+  }
+
+  // Check for infinity solutions case, for this:
+  // Sort points lexicographically and check midpoint between each two neighboring points.
+  // If all midpoints functional value is acceptable
+  // then set myParallel flag to true and return one soulution.
+  std::sort(aPnts.begin(), aPnts.end(), comp);
+  Standard_Boolean isParallel = Standard_False;
+  Standard_Real aVal = 0.0;
+  math_Vector aVec(1,2, 0.0);
+
+  // Avoid mark parallel case when have duplicates out of tolerance.
+  // Bad conditioned task: bug25635_1, bug23706_10, bug23706_13.
+  const Standard_Integer aMinNbInfSol = 100;
+  if (aPnts.Size() >= aMinNbInfSol)
+  {
+    isParallel = Standard_True;
+    for(Standard_Integer anIdx = aPnts.Lower(); anIdx <= aPnts.Upper() - 1; anIdx++)
+    {
+      const gp_XY& aCurrent = aPnts(anIdx);
+      const gp_XY& aNext    = aPnts(anIdx + 1);
+
+      aVec(1) = (aCurrent.X() + aNext.X()) * 0.5;
+      aVec(2) = (aCurrent.Y() + aNext.Y()) * 0.5;
+
+      aFunc.Value(aVec, aVal);
+
+      if (Abs(aVal - aF) > Precision::Confusion())
+      {
+        isParallel = Standard_False;
+        break;
       }
     }
   }
+
+  if (isParallel)
+  {
+    const gp_XY& aCurrent = aPnts.First();
+    myPoints1.Append(aCurrent.X());
+    myPoints2.Append(aCurrent.Y());
+    myParallel = Standard_True;
+  }
+  else
+  {
+    for(Standard_Integer anIdx = aPnts.Lower(); anIdx <= aPnts.Upper(); anIdx++)
+    {
+      const gp_XY& aCurrent = aPnts(anIdx);
+      myPoints1.Append(aCurrent.X());
+      myPoints2.Append(aCurrent.Y());
+    }
   }
 
   myDone = Standard_True;
@@ -206,6 +342,15 @@ Standard_Boolean Extrema_GenExtCC::IsDone() const
   return myDone; 
 }
 
+//=======================================================================
+//function : IsParallel
+//purpose  : 
+//=======================================================================
+Standard_Boolean Extrema_GenExtCC::IsParallel() const 
+{
+  return myParallel; 
+}
+
 //=======================================================================
 //function : NbExt
 //purpose  : 

src/math/math_GlobOptMin.cxx

@@ -497,13 +497,15 @@ Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
 {
   Standard_Integer i,j;
   Standard_Boolean isSame = Standard_True;
+  math_Vector aTol(1, myN);
+  aTol = (myB -  myA) * mySameTol;
 
   for(i = 0; i < mySolCount; i++)
   {
     isSame = Standard_True;
     for(j = 1; j <= myN; j++)
     {
-      if ((Abs(thePnt(j) - myY(i * myN + j))) > (myB(j) -  myA(j)) * mySameTol)
+      if ((Abs(thePnt(j) - myY(i * myN + j))) > aTol(j))
       {
         isSame = Standard_False;
         break;

tests/bugs/modalg_6/bug26075

@@ -0,0 +1,18 @@
+puts "========"
+puts "OCC26075"
+puts "========"
+puts ""
+###########################################################################
+# Make Extrema_GenExtCC return IsParallel flag in case of parallel curves
+###########################################################################
+
+restore [locate_data_file dist1-s1.brep] s1
+restore [locate_data_file dist1-s2.brep] s2
+mkcurve c1 s1
+mkcurve c2 s2
+set bug_info [extrema c1 c2]
+
+set bug_info [string range $bug_info [expr {[string first "\n" $bug_info] + 1}] [expr {[string last "\n" $bug_info] - 1}]]
+if {$bug_info != "No solutions!"} {
+  puts "ERROR: OCC26075 is reproduced. Flag IsParallel is not returned."
+}