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 : myIsFindSingleSolution(Standard_False),
 101   myParallel(Standard_False),
 102   myCurveMinTol(Precision::PConfusion()),
 103   myLowBorder(1,2),
 104   myUppBorder(1,2),
 105   myDone(Standard_False)
 106 {
 107   myC[0] = myC[1] = 0;
 108 }
 109
 110 //=======================================================================
 111 //function : Extrema_GenExtCC
 112 //purpose  :
 113 //=======================================================================
 114 Extrema_GenExtCC::Extrema_GenExtCC(const Curve1& C1,
 115                                    const Curve2& C2)
 116 : myIsFindSingleSolution(Standard_False),
 117   myParallel(Standard_False),
 118   myCurveMinTol(Precision::PConfusion()),
 119   myLowBorder(1,2),
 120   myUppBorder(1,2),
 121   myDone(Standard_False)
 122 {
 123   myC[0] = (Standard_Address)&C1;
 124   myC[1] = (Standard_Address)&C2;
 125   myLowBorder(1) = C1.FirstParameter();
 126   myLowBorder(2) = C2.FirstParameter();
 127   myUppBorder(1) = C1.LastParameter();
 128   myUppBorder(2) = C2.LastParameter();
 129 }
 130
 131 //=======================================================================
 132 //function : Extrema_GenExtCC
 133 //purpose  :
 134 //=======================================================================
 135 Extrema_GenExtCC::Extrema_GenExtCC(const Curve1& C1,
 136                                    const Curve2& C2,
 137                                    const Standard_Real Uinf,
 138                                    const Standard_Real Usup,
 139                                    const Standard_Real Vinf,
 140                                    const Standard_Real Vsup)
 141 : myIsFindSingleSolution(Standard_False),
 142   myParallel(Standard_False),
 143   myCurveMinTol(Precision::PConfusion()),
 144   myLowBorder(1,2),
 145   myUppBorder(1,2),
 146   myDone(Standard_False)
 147 {
 148   myC[0] = (Standard_Address)&C1;
 149   myC[1] = (Standard_Address)&C2;
 150   myLowBorder(1) = Uinf;
 151   myLowBorder(2) = Vinf;
 152   myUppBorder(1) = Usup;
 153   myUppBorder(2) = Vsup;
 154 }
 155
 156 //=======================================================================
 157 //function : SetParams
 158 //purpose  :
 159 //=======================================================================
 160 void Extrema_GenExtCC::SetParams(const Curve1& C1,
 161                                  const Curve2& C2,
 162                                  const Standard_Real Uinf,
 163                                  const Standard_Real Usup,
 164                                  const Standard_Real Vinf,
 165                                  const Standard_Real Vsup)
 166 {
 167   myC[0] = (Standard_Address)&C1;
 168   myC[1] = (Standard_Address)&C2;
 169   myLowBorder(1) = Uinf;
 170   myLowBorder(2) = Vinf;
 171   myUppBorder(1) = Usup;
 172   myUppBorder(2) = Vsup;
 173 }
 174
 175 //=======================================================================
 176 //function : SetTolerance
 177 //purpose  :
 178 //=======================================================================
 179 void Extrema_GenExtCC::SetTolerance(Standard_Real theTol)
 180 {
 181   myCurveMinTol = theTol;
 182 }
 183
 184 //=======================================================================
 185 //function : Perform
 186 //purpose  :
 187 //=======================================================================
 188 void Extrema_GenExtCC::Perform()
 189 {
 190   myDone = Standard_False;
 191   myParallel = Standard_False;
 192
 193   Curve1 &C1 = *(Curve1*)myC[0];
 194   Curve2 &C2 = *(Curve2*)myC[1];
 195
 196   Standard_Integer aNbInter[2];
 197   GeomAbs_Shape aContinuity = GeomAbs_C2;
 198   aNbInter[0] = C1.NbIntervals(aContinuity);
 199   aNbInter[1] = C2.NbIntervals(aContinuity);
 200
 201   if (aNbInter[0] * aNbInter[1] > 100)
 202   {
 203     aContinuity = GeomAbs_C1;
 204     aNbInter[0] = C1.NbIntervals(aContinuity);
 205     aNbInter[1] = C2.NbIntervals(aContinuity);
 206   }
 207
 208   TColStd_Array1OfReal anIntervals1(1, aNbInter[0] + 1);
 209   TColStd_Array1OfReal anIntervals2(1, aNbInter[1] + 1);
 210   C1.Intervals(anIntervals1, aContinuity);
 211   C2.Intervals(anIntervals2, aContinuity);
 212
 213   // Lipchitz constant computation.
 214   const Standard_Real aMaxLC = 10000.;
 215   Standard_Real aLC = 9.0; // Default value.
 216   const Standard_Real aMaxDer1 = 1.0 / C1.Resolution(1.0);
 217   const Standard_Real aMaxDer2 = 1.0 / C2.Resolution(1.0);
 218   Standard_Real aMaxDer = Max(aMaxDer1, aMaxDer2) * Sqrt(2.0);
 219   if (aLC > aMaxDer)
 220     aLC = aMaxDer;
 221
 222   // Change constant value according to the concrete curve types.
 223   Standard_Boolean isConstLockedFlag = Standard_False;
 224   //To prevent LipConst to became too small
 225   const Standard_Real aCR = 0.001;
 226   if (aMaxDer1 / aMaxDer < aCR || aMaxDer2 / aMaxDer < aCR)
 227   {
 228     isConstLockedFlag = Standard_True;
 229   }
 230   if (aMaxDer > aMaxLC)
 231   {
 232     aLC = aMaxLC;
 233     isConstLockedFlag = Standard_True;
 234   }
 235   if (C1.GetType() == GeomAbs_Line)
 236   {
 237     aMaxDer = 1.0 / C2.Resolution(1.0);
 238     if (aLC > aMaxDer)
 239     {
 240       isConstLockedFlag = Standard_True;
 241       aLC = aMaxDer;
 242     }
 243   }
 244   if (C2.GetType() == GeomAbs_Line)
 245   {
 246     aMaxDer = 1.0 / C1.Resolution(1.0);
 247     if (aLC > aMaxDer)
 248     {
 249       isConstLockedFlag = Standard_True;
 250       aLC = aMaxDer;
 251     }
 252   }
 253
 254   Extrema_GlobOptFuncCCC2 aFunc (C1, C2);
 255   math_GlobOptMin aFinder(&aFunc, myLowBorder, myUppBorder, aLC);
 256   aFinder.SetLipConstState(isConstLockedFlag);
 257   aFinder.SetContinuity(aContinuity == GeomAbs_C2 ? 2 : 1);
 258   Standard_Real aDiscTol = 1.0e-2;
 259   Standard_Real aValueTol = 1.0e-2;
 260   Standard_Real aSameTol = myCurveMinTol / (aDiscTol);
 261   aFinder.SetTol(aDiscTol, aSameTol);
 262   aFinder.SetFunctionalMinimalValue(0.0); // Best distance cannot be lower than 0.0.
 263
 264   // Size computed to have cell index inside of int32 value.
 265   const Standard_Real aCellSize = Max(anIntervals1.Last() - anIntervals1.First(),
 266                                       anIntervals2.Last() - anIntervals2.First())
 267                                   * Precision::PConfusion() / (2.0 * Sqrt(2.0));
 268   Extrema_CCPointsInspector anInspector(aCellSize);
 269   NCollection_CellFilter<Extrema_CCPointsInspector> aFilter(aCellSize);
 270   NCollection_Vector<gp_XY> aPnts;
 271
 272   Standard_Integer i,j,k;
 273   math_Vector aFirstBorderInterval(1,2);
 274   math_Vector aSecondBorderInterval(1,2);
 275   Standard_Real aF = RealLast(); // Best functional value.
 276   Standard_Real aCurrF = RealLast(); // Current functional value computed on current interval.
 277   for(i = 1; i <= aNbInter[0]; i++)
 278   {
 279     for(j = 1; j <= aNbInter[1]; j++)
 280     {
 281       aFirstBorderInterval(1) = anIntervals1(i);
 282       aFirstBorderInterval(2) = anIntervals2(j);
 283       aSecondBorderInterval(1) = anIntervals1(i + 1);
 284       aSecondBorderInterval(2) = anIntervals2(j + 1);
 285
 286       aFinder.SetLocalParams(aFirstBorderInterval, aSecondBorderInterval);
 287       aFinder.Perform(GetSingleSolutionFlag());
 288
 289       // Check that solution found on current interval is not worse than previous.
 290       aCurrF = aFinder.GetF();
 291       if (aCurrF >= aF + aSameTol * aValueTol)
 292       {
 293         continue;
 294       }
 295
 296       // Clean previously computed solution if current one is better.
 297       if (aCurrF > aF - aSameTol * aValueTol)
 298       {
 299         if (aCurrF < aF)
 300           aF = aCurrF;
 301       }
 302       else
 303       {
 304         aF = aCurrF;
 305         aFilter.Reset(aCellSize);
 306         aPnts.Clear();
 307       }
 308
 309       // Save found solutions avoiding repetitions.
 310       math_Vector sol(1,2);
 311       for(k = 1; k <= aFinder.NbExtrema(); k++)
 312       {
 313         aFinder.Points(k, sol);
 314         gp_XY aPnt2d(sol(1), sol(2));
 315
 316         gp_XY aXYmin = anInspector.Shift(aPnt2d, -aCellSize);
 317         gp_XY aXYmax = anInspector.Shift(aPnt2d,  aCellSize);
 318
 319         anInspector.ClearFind();
 320         anInspector.SetCurrent(aPnt2d);
 321         aFilter.Inspect(aXYmin, aXYmax, anInspector);
 322         if (!anInspector.isFind())
 323         {
 324           // Point is out of close cells, add new one.
 325           aFilter.Add(aPnt2d, aPnt2d);
 326           aPnts.Append(gp_XY(sol(1), sol(2)));
 327         }
 328       }
 329     }
 330   }
 331
 332   if (aPnts.Size() == 0)
 333   {
 334     // No solutions.
 335     myDone = Standard_False;
 336     return;
 337   }
 338
 339   // Check for infinity solutions case, for this:
 340   // Sort points lexicographically and check midpoint between each two neighboring points.
 341   // If all midpoints functional value is acceptable
 342   // then set myParallel flag to true and return one solution.
 343   std::sort(aPnts.begin(), aPnts.end(), comp);
 344   Standard_Boolean isParallel = Standard_False;
 345   Standard_Real aVal = 0.0;
 346   math_Vector aVec(1,2, 0.0);
 347
 348   // Avoid mark parallel case when have duplicates out of tolerance.
 349   // Bad conditioned task: bug25635_1, bug23706_10, bug23706_13.
 350   if (aPnts.Size() >= 2)
 351   {
 352     isParallel = Standard_True;
 353     for(Standard_Integer anIdx = aPnts.Lower(); anIdx <= aPnts.Upper() - 1; anIdx++)
 354     {
 355       const gp_XY& aCurrent = aPnts(anIdx);
 356       const gp_XY& aNext    = aPnts(anIdx + 1);
 357
 358       aVec(1) = (aCurrent.X() + aNext.X()) * 0.5;
 359       aVec(2) = (aCurrent.Y() + aNext.Y()) * 0.5;
 360
 361       aFunc.Value(aVec, aVal);
 362
 363       if (Abs(aVal - aF) > Precision::Confusion())
 364       {
 365         isParallel = Standard_False;
 366         break;
 367       }
 368     }
 369   }
 370
 371   if (isParallel)
 372   {
 373     const gp_XY& aCurrent = aPnts.First();
 374     myPoints1.Append(aCurrent.X());
 375     myPoints2.Append(aCurrent.Y());
 376     myParallel = Standard_True;
 377   }
 378   else
 379   {
 380     for(Standard_Integer anIdx = aPnts.Lower(); anIdx <= aPnts.Upper(); anIdx++)
 381     {
 382       const gp_XY& aCurrent = aPnts(anIdx);
 383       myPoints1.Append(aCurrent.X());
 384       myPoints2.Append(aCurrent.Y());
 385     }
 386   }
 387
 388   myDone = Standard_True;
 389 }
 390
 391 //=======================================================================
 392 //function : IsDone
 393 //purpose  :
 394 //=======================================================================
 395 Standard_Boolean Extrema_GenExtCC::IsDone() const
 396 {
 397   return myDone;
 398 }
 399
 400 //=======================================================================
 401 //function : IsParallel
 402 //purpose  :
 403 //=======================================================================
 404 Standard_Boolean Extrema_GenExtCC::IsParallel() const
 405 {
 406   return myParallel;
 407 }
 408
 409 //=======================================================================
 410 //function : NbExt
 411 //purpose  :
 412 //=======================================================================
 413 Standard_Integer Extrema_GenExtCC::NbExt() const
 414 {
 415   StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::NbExt()")
 416
 417   return myPoints1.Length();
 418 }
 419
 420 //=======================================================================
 421 //function : SquareDistance
 422 //purpose  :
 423 //=======================================================================
 424 Standard_Real Extrema_GenExtCC::SquareDistance(const Standard_Integer N) const
 425 {
 426   StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::SquareDistance()")
 427   Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::SquareDistance()")
 428
 429   return Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)).SquareDistance(Tool2::Value(*((Curve2*)myC[1]), myPoints2(N)));
 430 }
 431
 432 //=======================================================================
 433 //function : Points
 434 //purpose  :
 435 //=======================================================================
 436 void Extrema_GenExtCC::Points(const Standard_Integer N,
 437                               POnC& P1,
 438                               POnC& P2) const
 439 {
 440   StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::Points()")
 441   Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::Points()")
 442
 443   P1.SetValues(myPoints1(N), Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)));
 444   P2.SetValues(myPoints2(N), Tool2::Value(*((Curve2*)myC[1]), myPoints2(N)));
 445 }
 446
 447 //=======================================================================
 448 //function : SetSingleSolutionFlag
 449 //purpose  :
 450 //=======================================================================
 451 void Extrema_GenExtCC::SetSingleSolutionFlag(const Standard_Boolean theFlag)
 452 {
 453   myIsFindSingleSolution = theFlag;
 454 }
 455
 456 //=======================================================================
 457 //function : GetSingleSolutionFlag
 458 //purpose  :
 459 //=======================================================================
 460 Standard_Boolean Extrema_GenExtCC::GetSingleSolutionFlag() const
 461 {
 462   return myIsFindSingleSolution;
 463 }