| 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 approximation. |
| 214 | Standard_Real aLC = 9.0; // Default value. |
| 215 | Standard_Boolean isConstLockedFlag = Standard_False; |
| 216 | if (C1.GetType() == GeomAbs_Line) |
| 217 | { |
| 218 | Standard_Real aMaxDer = 1.0 / C2.Resolution(1.0); |
| 219 | if (aLC > aMaxDer) |
| 220 | { |
| 221 | isConstLockedFlag = Standard_True; |
| 222 | aLC = aMaxDer; |
| 223 | } |
| 224 | } |
| 225 | if (C2.GetType() == GeomAbs_Line) |
| 226 | { |
| 227 | Standard_Real aMaxDer = 1.0 / C1.Resolution(1.0); |
| 228 | if (aLC > aMaxDer) |
| 229 | { |
| 230 | isConstLockedFlag = Standard_True; |
| 231 | aLC = aMaxDer; |
| 232 | } |
| 233 | } |
| 234 | |
| 235 | Extrema_GlobOptFuncCCC2 aFunc (C1, C2); |
| 236 | math_GlobOptMin aFinder(&aFunc, myLowBorder, myUppBorder, aLC); |
| 237 | aFinder.SetLipConstState(isConstLockedFlag); |
| 238 | aFinder.SetContinuity(aContinuity == GeomAbs_C2 ? 2 : 1); |
| 239 | Standard_Real aDiscTol = 1.0e-2; |
| 240 | Standard_Real aValueTol = 1.0e-2; |
| 241 | Standard_Real aSameTol = myCurveMinTol / (aDiscTol); |
| 242 | aFinder.SetTol(aDiscTol, aSameTol); |
| 243 | aFinder.SetFunctionalMinimalValue(0.0); // Best distance cannot be lower than 0.0. |
| 244 | |
| 245 | // Size computed to have cell index inside of int32 value. |
| 246 | const Standard_Real aCellSize = Max(anIntervals1.Upper() - anIntervals1.Lower(), |
| 247 | anIntervals2.Upper() - anIntervals2.Lower()) |
| 248 | * Precision::PConfusion() / (2.0 * Sqrt(2.0)); |
| 249 | Extrema_CCPointsInspector anInspector(Precision::PConfusion()); |
| 250 | NCollection_CellFilter<Extrema_CCPointsInspector> aFilter(aCellSize); |
| 251 | NCollection_Vector<gp_XY> aPnts; |
| 252 | |
| 253 | Standard_Integer i,j,k; |
| 254 | math_Vector aFirstBorderInterval(1,2); |
| 255 | math_Vector aSecondBorderInterval(1,2); |
| 256 | Standard_Real aF = RealLast(); // Best functional value. |
| 257 | Standard_Real aCurrF = RealLast(); // Current functional value computed on current interval. |
| 258 | for(i = 1; i <= aNbInter[0]; i++) |
| 259 | { |
| 260 | for(j = 1; j <= aNbInter[1]; j++) |
| 261 | { |
| 262 | aFirstBorderInterval(1) = anIntervals1(i); |
| 263 | aFirstBorderInterval(2) = anIntervals2(j); |
| 264 | aSecondBorderInterval(1) = anIntervals1(i + 1); |
| 265 | aSecondBorderInterval(2) = anIntervals2(j + 1); |
| 266 | |
| 267 | aFinder.SetLocalParams(aFirstBorderInterval, aSecondBorderInterval); |
| 268 | aFinder.Perform(GetSingleSolutionFlag()); |
| 269 | |
| 270 | // Check that solution found on current interval is not worse than previous. |
| 271 | aCurrF = aFinder.GetF(); |
| 272 | if (aCurrF >= aF + aSameTol * aValueTol) |
| 273 | { |
| 274 | continue; |
| 275 | } |
| 276 | |
| 277 | // Clean previously computed solution if current one is better. |
| 278 | if (aCurrF > aF - aSameTol * aValueTol) |
| 279 | { |
| 280 | if (aCurrF < aF) |
| 281 | aF = aCurrF; |
| 282 | } |
| 283 | else |
| 284 | { |
| 285 | aF = aCurrF; |
| 286 | aFilter.Reset(aCellSize); |
| 287 | aPnts.Clear(); |
| 288 | } |
| 289 | |
| 290 | // Save found solutions avoiding repetitions. |
| 291 | math_Vector sol(1,2); |
| 292 | for(k = 1; k <= aFinder.NbExtrema(); k++) |
| 293 | { |
| 294 | aFinder.Points(k, sol); |
| 295 | gp_XY aPnt2d(sol(1), sol(2)); |
| 296 | |
| 297 | gp_XY aXYmin = anInspector.Shift(aPnt2d, -aCellSize); |
| 298 | gp_XY aXYmax = anInspector.Shift(aPnt2d, aCellSize); |
| 299 | |
| 300 | anInspector.ClearFind(); |
| 301 | anInspector.SetCurrent(aPnt2d); |
| 302 | aFilter.Inspect(aXYmin, aXYmax, anInspector); |
| 303 | if (!anInspector.isFind()) |
| 304 | { |
| 305 | // Point is out of close cells, add new one. |
| 306 | aFilter.Add(aPnt2d, aPnt2d); |
| 307 | aPnts.Append(gp_XY(sol(1), sol(2))); |
| 308 | } |
| 309 | } |
| 310 | } |
| 311 | } |
| 312 | |
| 313 | if (aPnts.Size() == 0) |
| 314 | { |
| 315 | // No solutions. |
| 316 | myDone = Standard_False; |
| 317 | return; |
| 318 | } |
| 319 | |
| 320 | // Check for infinity solutions case, for this: |
| 321 | // Sort points lexicographically and check midpoint between each two neighboring points. |
| 322 | // If all midpoints functional value is acceptable |
| 323 | // then set myParallel flag to true and return one solution. |
| 324 | std::sort(aPnts.begin(), aPnts.end(), comp); |
| 325 | Standard_Boolean isParallel = Standard_False; |
| 326 | Standard_Real aVal = 0.0; |
| 327 | math_Vector aVec(1,2, 0.0); |
| 328 | |
| 329 | // Avoid mark parallel case when have duplicates out of tolerance. |
| 330 | // Bad conditioned task: bug25635_1, bug23706_10, bug23706_13. |
| 331 | const Standard_Integer aMinNbInfSol = 100; |
| 332 | if (aPnts.Size() >= aMinNbInfSol) |
| 333 | { |
| 334 | isParallel = Standard_True; |
| 335 | for(Standard_Integer anIdx = aPnts.Lower(); anIdx <= aPnts.Upper() - 1; anIdx++) |
| 336 | { |
| 337 | const gp_XY& aCurrent = aPnts(anIdx); |
| 338 | const gp_XY& aNext = aPnts(anIdx + 1); |
| 339 | |
| 340 | aVec(1) = (aCurrent.X() + aNext.X()) * 0.5; |
| 341 | aVec(2) = (aCurrent.Y() + aNext.Y()) * 0.5; |
| 342 | |
| 343 | aFunc.Value(aVec, aVal); |
| 344 | |
| 345 | if (Abs(aVal - aF) > Precision::Confusion()) |
| 346 | { |
| 347 | isParallel = Standard_False; |
| 348 | break; |
| 349 | } |
| 350 | } |
| 351 | } |
| 352 | |
| 353 | if (isParallel) |
| 354 | { |
| 355 | const gp_XY& aCurrent = aPnts.First(); |
| 356 | myPoints1.Append(aCurrent.X()); |
| 357 | myPoints2.Append(aCurrent.Y()); |
| 358 | myParallel = Standard_True; |
| 359 | } |
| 360 | else |
| 361 | { |
| 362 | for(Standard_Integer anIdx = aPnts.Lower(); anIdx <= aPnts.Upper(); anIdx++) |
| 363 | { |
| 364 | const gp_XY& aCurrent = aPnts(anIdx); |
| 365 | myPoints1.Append(aCurrent.X()); |
| 366 | myPoints2.Append(aCurrent.Y()); |
| 367 | } |
| 368 | } |
| 369 | |
| 370 | myDone = Standard_True; |
| 371 | } |
| 372 | |
| 373 | //======================================================================= |
| 374 | //function : IsDone |
| 375 | //purpose : |
| 376 | //======================================================================= |
| 377 | Standard_Boolean Extrema_GenExtCC::IsDone() const |
| 378 | { |
| 379 | return myDone; |
| 380 | } |
| 381 | |
| 382 | //======================================================================= |
| 383 | //function : IsParallel |
| 384 | //purpose : |
| 385 | //======================================================================= |
| 386 | Standard_Boolean Extrema_GenExtCC::IsParallel() const |
| 387 | { |
| 388 | return myParallel; |
| 389 | } |
| 390 | |
| 391 | //======================================================================= |
| 392 | //function : NbExt |
| 393 | //purpose : |
| 394 | //======================================================================= |
| 395 | Standard_Integer Extrema_GenExtCC::NbExt() const |
| 396 | { |
| 397 | StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::NbExt()") |
| 398 | |
| 399 | return myPoints1.Length(); |
| 400 | } |
| 401 | |
| 402 | //======================================================================= |
| 403 | //function : SquareDistance |
| 404 | //purpose : |
| 405 | //======================================================================= |
| 406 | Standard_Real Extrema_GenExtCC::SquareDistance(const Standard_Integer N) const |
| 407 | { |
| 408 | StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::SquareDistance()") |
| 409 | Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::SquareDistance()") |
| 410 | |
| 411 | return Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)).SquareDistance(Tool2::Value(*((Curve2*)myC[1]), myPoints2(N))); |
| 412 | } |
| 413 | |
| 414 | //======================================================================= |
| 415 | //function : Points |
| 416 | //purpose : |
| 417 | //======================================================================= |
| 418 | void Extrema_GenExtCC::Points(const Standard_Integer N, |
| 419 | POnC& P1, |
| 420 | POnC& P2) const |
| 421 | { |
| 422 | StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::Points()") |
| 423 | Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::Points()") |
| 424 | |
| 425 | P1.SetValues(myPoints1(N), Tool1::Value(*((Curve1*)myC[0]), myPoints1(N))); |
| 426 | P2.SetValues(myPoints2(N), Tool2::Value(*((Curve2*)myC[1]), myPoints2(N))); |
| 427 | } |
| 428 | |
| 429 | //======================================================================= |
| 430 | //function : SetSingleSolutionFlag |
| 431 | //purpose : |
| 432 | //======================================================================= |
| 433 | void Extrema_GenExtCC::SetSingleSolutionFlag(const Standard_Boolean theFlag) |
| 434 | { |
| 435 | myIsFindSingleSolution = theFlag; |
| 436 | } |
| 437 | |
| 438 | //======================================================================= |
| 439 | //function : GetSingleSolutionFlag |
| 440 | //purpose : |
| 441 | //======================================================================= |
| 442 | Standard_Boolean Extrema_GenExtCC::GetSingleSolutionFlag() const |
| 443 | { |
| 444 | return myIsFindSingleSolution; |
| 445 | } |