src/Extrema/Extrema_GenExtCC.gxx

   1 // Created on: 1995-07-18
   2 // Created by: Modelistation
   3 // Copyright (c) 1995-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <algorithm>
  18
  19 #include <Extrema_GlobOptFuncCC.hxx>
  20 #include <math_GlobOptMin.hxx>
  21 #include <Standard_NullObject.hxx>
  22 #include <Standard_OutOfRange.hxx>
  23 #include <StdFail_NotDone.hxx>
  24 #include <TColStd_Array1OfReal.hxx>
  25 #include <Precision.hxx>
  26 #include <NCollection_Vector.hxx>
  27 #include <NCollection_CellFilter.hxx>
  28
  29 // Comparator, used in std::sort.
  30 static Standard_Boolean comp(const gp_XY& theA,
  31                              const gp_XY& theB)
  32 {
  33   if (theA.X() < theB.X())
  34   {
  35     return Standard_True;
  36   }
  37   else
  38   {
  39     if (theA.X() == theB.X())
  40     {
  41       if (theA.Y() <= theB.Y())
  42         return Standard_True;
  43     }
  44   }
  45
  46   return Standard_False;
  47 }
  48
  49 class Extrema_CCPointsInspector : public NCollection_CellFilter_InspectorXY
  50 {
  51 public:
  52   typedef gp_XY Target;
  53   //! Constructor; remembers the tolerance
  54   Extrema_CCPointsInspector (const Standard_Real theTol)
  55   {
  56     myTol = theTol * theTol;
  57     myIsFind = Standard_False;
  58   }
  59
  60   void ClearFind()
  61   {
  62     myIsFind = Standard_False;
  63   }
  64
  65   Standard_Boolean isFind()
  66   {
  67     return myIsFind;
  68   }
  69
  70   //! Set current point to search for coincidence
  71   void SetCurrent (const gp_XY& theCurPnt)
  72   {
  73     myCurrent = theCurPnt;
  74   }
  75
  76   //! Implementation of inspection method
  77   NCollection_CellFilter_Action Inspect (const Target& theObject)
  78   {
  79     gp_XY aPt = myCurrent.Subtracted(theObject);
  80     const Standard_Real aSQDist = aPt.SquareModulus();
  81     if(aSQDist < myTol)
  82     {
  83       myIsFind = Standard_True;
  84     }
  85
  86     return CellFilter_Keep;
  87   }
  88
  89 private:
  90   Standard_Real myTol;
  91   gp_XY myCurrent;
  92   Standard_Boolean myIsFind;
  93 };
  94
  95 //=======================================================================
  96 //function : Extrema_GenExtCC
  97 //purpose  :
  98 //=======================================================================
  99 Extrema_GenExtCC::Extrema_GenExtCC()
 100 : myParallel(Standard_False),
 101   myCurveMinTol(Precision::PConfusion()),
 102   myLowBorder(1,2),
 103   myUppBorder(1,2),
 104   myDone(Standard_False)
 105 {
 106   myC[0] = myC[1] = 0;
 107 }
 108
 109 //=======================================================================
 110 //function : Extrema_GenExtCC
 111 //purpose  :
 112 //=======================================================================
 113 Extrema_GenExtCC::Extrema_GenExtCC(const Curve1& C1,
 114                                    const Curve2& C2)
 115 : myParallel(Standard_False),
 116   myCurveMinTol(Precision::PConfusion()),
 117   myLowBorder(1,2),
 118   myUppBorder(1,2),
 119   myDone(Standard_False)
 120 {
 121   myC[0] = (Standard_Address)&C1;
 122   myC[1] = (Standard_Address)&C2;
 123   myLowBorder(1) = C1.FirstParameter();
 124   myLowBorder(2) = C2.FirstParameter();
 125   myUppBorder(1) = C1.LastParameter();
 126   myUppBorder(2) = C2.LastParameter();
 127 }
 128
 129 //=======================================================================
 130 //function : Extrema_GenExtCC
 131 //purpose  :
 132 //=======================================================================
 133 Extrema_GenExtCC::Extrema_GenExtCC(const Curve1& C1,
 134                                    const Curve2& C2,
 135                                    const Standard_Real Uinf,
 136                                    const Standard_Real Usup,
 137                                    const Standard_Real Vinf,
 138                                    const Standard_Real Vsup)
 139 : myParallel(Standard_False),
 140   myCurveMinTol(Precision::PConfusion()),
 141   myLowBorder(1,2),
 142   myUppBorder(1,2),
 143   myDone(Standard_False)
 144 {
 145   myC[0] = (Standard_Address)&C1;
 146   myC[1] = (Standard_Address)&C2;
 147   myLowBorder(1) = Uinf;
 148   myLowBorder(2) = Vinf;
 149   myUppBorder(1) = Usup;
 150   myUppBorder(2) = Vsup;
 151 }
 152
 153 //=======================================================================
 154 //function : SetParams
 155 //purpose  :
 156 //=======================================================================
 157 void Extrema_GenExtCC::SetParams(const Curve1& C1,
 158                                  const Curve2& C2,
 159                                  const Standard_Real Uinf,
 160                                  const Standard_Real Usup,
 161                                  const Standard_Real Vinf,
 162                                  const Standard_Real Vsup)
 163 {
 164   myC[0] = (Standard_Address)&C1;
 165   myC[1] = (Standard_Address)&C2;
 166   myLowBorder(1) = Uinf;
 167   myLowBorder(2) = Vinf;
 168   myUppBorder(1) = Usup;
 169   myUppBorder(2) = Vsup;
 170 }
 171
 172 //=======================================================================
 173 //function : SetTolerance
 174 //purpose  :
 175 //=======================================================================
 176 void Extrema_GenExtCC::SetTolerance(Standard_Real theTol)
 177 {
 178   myCurveMinTol = theTol;
 179 }
 180
 181 //=======================================================================
 182 //function : Perform
 183 //purpose  :
 184 //=======================================================================
 185 void Extrema_GenExtCC::Perform()
 186 {
 187   myDone = Standard_False;
 188   myParallel = Standard_False;
 189
 190   Curve1 &C1 = *(Curve1*)myC[0];
 191   Curve2 &C2 = *(Curve2*)myC[1];
 192
 193   Standard_Integer aNbInter[2];
 194   aNbInter[0] = C1.NbIntervals(GeomAbs_C2);
 195   aNbInter[1] = C2.NbIntervals(GeomAbs_C2);
 196   TColStd_Array1OfReal anIntervals1(1, aNbInter[0] + 1);
 197   TColStd_Array1OfReal anIntervals2(1, aNbInter[1] + 1);
 198   C1.Intervals(anIntervals1, GeomAbs_C2);
 199   C2.Intervals(anIntervals2, GeomAbs_C2);
 200
 201   Extrema_GlobOptFuncCCC2 aFunc (C1, C2);
 202   math_GlobOptMin aFinder(&aFunc, myLowBorder, myUppBorder);
 203   Standard_Real aDiscTol = 1.0e-2;
 204   Standard_Real aValueTol = 1.0e-2;
 205   Standard_Real aSameTol = myCurveMinTol / (aDiscTol);
 206   aFinder.SetTol(aDiscTol, aSameTol);
 207
 208   // Size computed to have cell index inside of int32 value.
 209   const Standard_Real aCellSize = Max(anIntervals1.Upper() - anIntervals1.Lower(),
 210                                       anIntervals2.Upper() - anIntervals2.Lower())
 211                                   * Precision::PConfusion() / (2.0 * Sqrt(2.0));
 212   Extrema_CCPointsInspector anInspector(Precision::PConfusion());
 213   NCollection_CellFilter<Extrema_CCPointsInspector> aFilter(aCellSize);
 214   NCollection_Vector<gp_XY> aPnts;
 215
 216   Standard_Integer i,j,k;
 217   math_Vector aFirstBorderInterval(1,2);
 218   math_Vector aSecondBorderInterval(1,2);
 219   Standard_Real aF = RealLast(); // Best functional value.
 220   Standard_Real aCurrF = RealLast(); // Current functional value computed on current interval.
 221   for(i = 1; i <= aNbInter[0]; i++)
 222   {
 223     for(j = 1; j <= aNbInter[1]; j++)
 224     {
 225       aFirstBorderInterval(1) = anIntervals1(i);
 226       aFirstBorderInterval(2) = anIntervals2(j);
 227       aSecondBorderInterval(1) = anIntervals1(i + 1);
 228       aSecondBorderInterval(2) = anIntervals2(j + 1);
 229
 230       aFinder.SetLocalParams(aFirstBorderInterval, aSecondBorderInterval);
 231       aFinder.Perform();
 232
 233       // Check that solution found on current interval is not worse than previous.
 234       aCurrF = aFinder.GetF();
 235       if (aCurrF >= aF + aSameTol * aValueTol)
 236       {
 237         continue;
 238       }
 239
 240       // Clean previously computed solution if current one is better.
 241       if (aCurrF > aF - aSameTol * aValueTol)
 242       {
 243         if (aCurrF < aF)
 244           aF = aCurrF;
 245       }
 246       else
 247       {
 248         aF = aCurrF;
 249         aFilter.Reset(aCellSize);
 250         aPnts.Clear();
 251       }
 252
 253       // Save found solutions avoiding repetitions.
 254       math_Vector sol(1,2);
 255       for(k = 1; k <= aFinder.NbExtrema(); k++)
 256       {
 257         aFinder.Points(k, sol);
 258         gp_XY aPnt2d(sol(1), sol(2));
 259
 260         gp_XY aXYmin = anInspector.Shift(aPnt2d, -aCellSize);
 261         gp_XY aXYmax = anInspector.Shift(aPnt2d,  aCellSize);
 262
 263         anInspector.ClearFind();
 264         anInspector.SetCurrent(aPnt2d);
 265         aFilter.Inspect(aXYmin, aXYmax, anInspector);
 266         if (!anInspector.isFind())
 267         {
 268           // Point is out of close cells, add new one.
 269           aFilter.Add(aPnt2d, aPnt2d);
 270           aPnts.Append(gp_XY(sol(1), sol(2)));
 271         }
 272       }
 273     }
 274   }
 275
 276   if (aPnts.Size() == 0)
 277   {
 278     // No solutions.
 279     myDone = Standard_False;
 280     return;
 281   }
 282
 283   // Check for infinity solutions case, for this:
 284   // Sort points lexicographically and check midpoint between each two neighboring points.
 285   // If all midpoints functional value is acceptable
 286   // then set myParallel flag to true and return one soulution.
 287   std::sort(aPnts.begin(), aPnts.end(), comp);
 288   Standard_Boolean isParallel = Standard_False;
 289   Standard_Real aVal = 0.0;
 290   math_Vector aVec(1,2, 0.0);
 291
 292   // Avoid mark parallel case when have duplicates out of tolerance.
 293   // Bad conditioned task: bug25635_1, bug23706_10, bug23706_13.
 294   const Standard_Integer aMinNbInfSol = 100;
 295   if (aPnts.Size() >= aMinNbInfSol)
 296   {
 297     isParallel = Standard_True;
 298     for(Standard_Integer anIdx = aPnts.Lower(); anIdx <= aPnts.Upper() - 1; anIdx++)
 299     {
 300       const gp_XY& aCurrent = aPnts(anIdx);
 301       const gp_XY& aNext    = aPnts(anIdx + 1);
 302
 303       aVec(1) = (aCurrent.X() + aNext.X()) * 0.5;
 304       aVec(2) = (aCurrent.Y() + aNext.Y()) * 0.5;
 305
 306       aFunc.Value(aVec, aVal);
 307
 308       if (Abs(aVal - aF) > Precision::Confusion())
 309       {
 310         isParallel = Standard_False;
 311         break;
 312       }
 313     }
 314   }
 315
 316   if (isParallel)
 317   {
 318     const gp_XY& aCurrent = aPnts.First();
 319     myPoints1.Append(aCurrent.X());
 320     myPoints2.Append(aCurrent.Y());
 321     myParallel = Standard_True;
 322   }
 323   else
 324   {
 325     for(Standard_Integer anIdx = aPnts.Lower(); anIdx <= aPnts.Upper(); anIdx++)
 326     {
 327       const gp_XY& aCurrent = aPnts(anIdx);
 328       myPoints1.Append(aCurrent.X());
 329       myPoints2.Append(aCurrent.Y());
 330     }
 331   }
 332
 333   myDone = Standard_True;
 334 }
 335
 336 //=======================================================================
 337 //function : IsDone
 338 //purpose  :
 339 //=======================================================================
 340 Standard_Boolean Extrema_GenExtCC::IsDone() const
 341 {
 342   return myDone;
 343 }
 344
 345 //=======================================================================
 346 //function : IsParallel
 347 //purpose  :
 348 //=======================================================================
 349 Standard_Boolean Extrema_GenExtCC::IsParallel() const
 350 {
 351   return myParallel;
 352 }
 353
 354 //=======================================================================
 355 //function : NbExt
 356 //purpose  :
 357 //=======================================================================
 358 Standard_Integer Extrema_GenExtCC::NbExt() const
 359 {
 360   StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::NbExt()")
 361
 362   return myPoints1.Length();
 363 }
 364
 365 //=======================================================================
 366 //function : SquareDistance
 367 //purpose  :
 368 //=======================================================================
 369 Standard_Real Extrema_GenExtCC::SquareDistance(const Standard_Integer N) const
 370 {
 371   StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::SquareDistance()")
 372   Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::SquareDistance()")
 373
 374   return Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)).SquareDistance(Tool2::Value(*((Curve2*)myC[1]), myPoints2(N)));
 375 }
 376
 377 //=======================================================================
 378 //function : Points
 379 //purpose  :
 380 //=======================================================================
 381 void Extrema_GenExtCC::Points(const Standard_Integer N,
 382                               POnC& P1,
 383                               POnC& P2) const
 384 {
 385   StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::Points()")
 386   Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::Points()")
 387
 388   P1.SetValues(myPoints1(N), Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)));
 389   P2.SetValues(myPoints2(N), Tool2::Value(*((Curve2*)myC[1]), myPoints2(N)));
 390 }