0027371: Regression: BRepExtrema works too much slower in 691 (from 670)
[occt.git] / src / Extrema / Extrema_GenExtCC.gxx
CommitLineData
b311480e 1// Created on: 1995-07-18
2// Created by: Modelistation
3// Copyright (c) 1995-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 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
973c2be1 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
20a216fe 17#include <algorithm>
18
4bbaf12b 19#include <Extrema_GlobOptFuncCC.hxx>
20#include <math_GlobOptMin.hxx>
7fd59977 21#include <Standard_NullObject.hxx>
22#include <Standard_OutOfRange.hxx>
4bbaf12b 23#include <StdFail_NotDone.hxx>
24#include <TColStd_Array1OfReal.hxx>
7fd59977 25#include <Precision.hxx>
20a216fe 26#include <NCollection_Vector.hxx>
27#include <NCollection_CellFilter.hxx>
28
29// Comparator, used in std::sort.
30static 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 {
afb27815 41 if (theA.Y() < theB.Y())
20a216fe 42 return Standard_True;
43 }
44 }
45
46 return Standard_False;
47}
48
49class Extrema_CCPointsInspector : public NCollection_CellFilter_InspectorXY
50{
51public:
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
89private:
90 Standard_Real myTol;
91 gp_XY myCurrent;
92 Standard_Boolean myIsFind;
93};
7fd59977 94
5493d334 95//=======================================================================
96//function : Extrema_GenExtCC
97//purpose :
98//=======================================================================
4bbaf12b 99Extrema_GenExtCC::Extrema_GenExtCC()
836d7b64 100: myIsFindSingleSolution(Standard_False),
101 myParallel(Standard_False),
20a216fe 102 myCurveMinTol(Precision::PConfusion()),
6ca1fa70 103 myLowBorder(1,2),
4bbaf12b 104 myUppBorder(1,2),
105 myDone(Standard_False)
7fd59977 106{
6ca1fa70 107 myC[0] = myC[1] = 0;
7fd59977 108}
109
5493d334 110//=======================================================================
111//function : Extrema_GenExtCC
112//purpose :
113//=======================================================================
114Extrema_GenExtCC::Extrema_GenExtCC(const Curve1& C1,
115 const Curve2& C2)
836d7b64 116: myIsFindSingleSolution(Standard_False),
117 myParallel(Standard_False),
20a216fe 118 myCurveMinTol(Precision::PConfusion()),
6ca1fa70 119 myLowBorder(1,2),
4bbaf12b 120 myUppBorder(1,2),
121 myDone(Standard_False)
7fd59977 122{
4bbaf12b 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();
7fd59977 129}
130
5493d334 131//=======================================================================
132//function : Extrema_GenExtCC
133//purpose :
134//=======================================================================
135Extrema_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)
836d7b64 141: myIsFindSingleSolution(Standard_False),
142 myParallel(Standard_False),
20a216fe 143 myCurveMinTol(Precision::PConfusion()),
6ca1fa70 144 myLowBorder(1,2),
4bbaf12b 145 myUppBorder(1,2),
146 myDone(Standard_False)
7fd59977 147{
4bbaf12b 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;
7fd59977 154}
155
5493d334 156//=======================================================================
157//function : SetParams
158//purpose :
159//=======================================================================
160void 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)
7fd59977 166{
4bbaf12b 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;
7fd59977 173}
174
5493d334 175//=======================================================================
176//function : SetTolerance
177//purpose :
178//=======================================================================
179void Extrema_GenExtCC::SetTolerance(Standard_Real theTol)
180{
181 myCurveMinTol = theTol;
182}
183
184//=======================================================================
185//function : Perform
186//purpose :
187//=======================================================================
188void Extrema_GenExtCC::Perform()
7fd59977 189{
190 myDone = Standard_False;
20a216fe 191 myParallel = Standard_False;
7fd59977 192
4bbaf12b 193 Curve1 &C1 = *(Curve1*)myC[0];
194 Curve2 &C2 = *(Curve2*)myC[1];
195
196 Standard_Integer aNbInter[2];
5333268d 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
4bbaf12b 208 TColStd_Array1OfReal anIntervals1(1, aNbInter[0] + 1);
209 TColStd_Array1OfReal anIntervals2(1, aNbInter[1] + 1);
5333268d 210 C1.Intervals(anIntervals1, aContinuity);
211 C2.Intervals(anIntervals2, aContinuity);
4bbaf12b 212
246c7a75 213 // Lipchitz constant computation.
1907fb9a 214 Standard_Real aLC = 9.0; // Default value.
246c7a75 215 const Standard_Real aMaxDer1 = 1.0 / C1.Resolution(1.0);
216 const Standard_Real aMaxDer2 = 1.0 / C2.Resolution(1.0);
217 const Standard_Real aMaxDer = Max(aMaxDer1, aMaxDer2) * Sqrt(2.0);
218 if (aLC > aMaxDer)
219 aLC = aMaxDer;
220
221 // Change constant value according to the concrete curve types.
1907fb9a 222 Standard_Boolean isConstLockedFlag = Standard_False;
223 if (C1.GetType() == GeomAbs_Line)
224 {
225 Standard_Real aMaxDer = 1.0 / C2.Resolution(1.0);
226 if (aLC > aMaxDer)
227 {
228 isConstLockedFlag = Standard_True;
229 aLC = aMaxDer;
230 }
231 }
232 if (C2.GetType() == GeomAbs_Line)
233 {
234 Standard_Real aMaxDer = 1.0 / C1.Resolution(1.0);
235 if (aLC > aMaxDer)
236 {
237 isConstLockedFlag = Standard_True;
238 aLC = aMaxDer;
239 }
240 }
241
6ca1fa70 242 Extrema_GlobOptFuncCCC2 aFunc (C1, C2);
1907fb9a 243 math_GlobOptMin aFinder(&aFunc, myLowBorder, myUppBorder, aLC);
244 aFinder.SetLipConstState(isConstLockedFlag);
5333268d 245 aFinder.SetContinuity(aContinuity == GeomAbs_C2 ? 2 : 1);
5493d334 246 Standard_Real aDiscTol = 1.0e-2;
247 Standard_Real aValueTol = 1.0e-2;
248 Standard_Real aSameTol = myCurveMinTol / (aDiscTol);
249 aFinder.SetTol(aDiscTol, aSameTol);
1907fb9a 250 aFinder.SetFunctionalMinimalValue(0.0); // Best distance cannot be lower than 0.0.
4bbaf12b 251
20a216fe 252 // Size computed to have cell index inside of int32 value.
253 const Standard_Real aCellSize = Max(anIntervals1.Upper() - anIntervals1.Lower(),
254 anIntervals2.Upper() - anIntervals2.Lower())
255 * Precision::PConfusion() / (2.0 * Sqrt(2.0));
256 Extrema_CCPointsInspector anInspector(Precision::PConfusion());
a7653f4f 257 NCollection_CellFilter<Extrema_CCPointsInspector> aFilter(aCellSize);
20a216fe 258 NCollection_Vector<gp_XY> aPnts;
6ca1fa70 259
4bbaf12b 260 Standard_Integer i,j,k;
261 math_Vector aFirstBorderInterval(1,2);
262 math_Vector aSecondBorderInterval(1,2);
20a216fe 263 Standard_Real aF = RealLast(); // Best functional value.
264 Standard_Real aCurrF = RealLast(); // Current functional value computed on current interval.
4bbaf12b 265 for(i = 1; i <= aNbInter[0]; i++)
266 {
267 for(j = 1; j <= aNbInter[1]; j++)
268 {
269 aFirstBorderInterval(1) = anIntervals1(i);
270 aFirstBorderInterval(2) = anIntervals2(j);
271 aSecondBorderInterval(1) = anIntervals1(i + 1);
272 aSecondBorderInterval(2) = anIntervals2(j + 1);
273
274 aFinder.SetLocalParams(aFirstBorderInterval, aSecondBorderInterval);
836d7b64 275 aFinder.Perform(GetSingleSolutionFlag());
4bbaf12b 276
20a216fe 277 // Check that solution found on current interval is not worse than previous.
4bbaf12b 278 aCurrF = aFinder.GetF();
6ca1fa70 279 if (aCurrF >= aF + aSameTol * aValueTol)
4bbaf12b 280 {
6ca1fa70 281 continue;
282 }
283
20a216fe 284 // Clean previously computed solution if current one is better.
6ca1fa70 285 if (aCurrF > aF - aSameTol * aValueTol)
286 {
287 if (aCurrF < aF)
4bbaf12b 288 aF = aCurrF;
6ca1fa70 289 }
290 else
291 {
292 aF = aCurrF;
20a216fe 293 aFilter.Reset(aCellSize);
294 aPnts.Clear();
6ca1fa70 295 }
7fd59977 296
20a216fe 297 // Save found solutions avoiding repetitions.
6ca1fa70 298 math_Vector sol(1,2);
299 for(k = 1; k <= aFinder.NbExtrema(); k++)
4bbaf12b 300 {
6ca1fa70 301 aFinder.Points(k, sol);
20a216fe 302 gp_XY aPnt2d(sol(1), sol(2));
6ca1fa70 303
20a216fe 304 gp_XY aXYmin = anInspector.Shift(aPnt2d, -aCellSize);
305 gp_XY aXYmax = anInspector.Shift(aPnt2d, aCellSize);
306
307 anInspector.ClearFind();
308 anInspector.SetCurrent(aPnt2d);
309 aFilter.Inspect(aXYmin, aXYmax, anInspector);
310 if (!anInspector.isFind())
6ca1fa70 311 {
20a216fe 312 // Point is out of close cells, add new one.
313 aFilter.Add(aPnt2d, aPnt2d);
314 aPnts.Append(gp_XY(sol(1), sol(2)));
6ca1fa70 315 }
4bbaf12b 316 }
7fd59977 317 }
4bbaf12b 318 }
7fd59977 319
20a216fe 320 if (aPnts.Size() == 0)
321 {
322 // No solutions.
323 myDone = Standard_False;
324 return;
325 }
326
327 // Check for infinity solutions case, for this:
328 // Sort points lexicographically and check midpoint between each two neighboring points.
329 // If all midpoints functional value is acceptable
1907fb9a 330 // then set myParallel flag to true and return one solution.
20a216fe 331 std::sort(aPnts.begin(), aPnts.end(), comp);
332 Standard_Boolean isParallel = Standard_False;
333 Standard_Real aVal = 0.0;
334 math_Vector aVec(1,2, 0.0);
335
336 // Avoid mark parallel case when have duplicates out of tolerance.
337 // Bad conditioned task: bug25635_1, bug23706_10, bug23706_13.
338 const Standard_Integer aMinNbInfSol = 100;
339 if (aPnts.Size() >= aMinNbInfSol)
340 {
341 isParallel = Standard_True;
342 for(Standard_Integer anIdx = aPnts.Lower(); anIdx <= aPnts.Upper() - 1; anIdx++)
343 {
344 const gp_XY& aCurrent = aPnts(anIdx);
345 const gp_XY& aNext = aPnts(anIdx + 1);
346
347 aVec(1) = (aCurrent.X() + aNext.X()) * 0.5;
348 aVec(2) = (aCurrent.Y() + aNext.Y()) * 0.5;
349
350 aFunc.Value(aVec, aVal);
351
352 if (Abs(aVal - aF) > Precision::Confusion())
353 {
354 isParallel = Standard_False;
355 break;
356 }
357 }
358 }
359
360 if (isParallel)
361 {
362 const gp_XY& aCurrent = aPnts.First();
363 myPoints1.Append(aCurrent.X());
364 myPoints2.Append(aCurrent.Y());
365 myParallel = Standard_True;
366 }
367 else
368 {
369 for(Standard_Integer anIdx = aPnts.Lower(); anIdx <= aPnts.Upper(); anIdx++)
370 {
371 const gp_XY& aCurrent = aPnts(anIdx);
372 myPoints1.Append(aCurrent.X());
373 myPoints2.Append(aCurrent.Y());
374 }
375 }
376
7fd59977 377 myDone = Standard_True;
378}
7fd59977 379
5493d334 380//=======================================================================
381//function : IsDone
382//purpose :
383//=======================================================================
384Standard_Boolean Extrema_GenExtCC::IsDone() const
4bbaf12b 385{
386 return myDone;
387}
7fd59977 388
5493d334 389//=======================================================================
20a216fe 390//function : IsParallel
391//purpose :
392//=======================================================================
393Standard_Boolean Extrema_GenExtCC::IsParallel() const
394{
395 return myParallel;
396}
397
398//=======================================================================
5493d334 399//function : NbExt
400//purpose :
401//=======================================================================
402Standard_Integer Extrema_GenExtCC::NbExt() const
7fd59977 403{
404 StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::NbExt()")
4bbaf12b 405
6ca1fa70 406 return myPoints1.Length();
7fd59977 407}
7fd59977 408
5493d334 409//=======================================================================
410//function : SquareDistance
411//purpose :
412//=======================================================================
413Standard_Real Extrema_GenExtCC::SquareDistance(const Standard_Integer N) const
7fd59977 414{
415 StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::SquareDistance()")
416 Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::SquareDistance()")
4bbaf12b 417
418 return Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)).SquareDistance(Tool2::Value(*((Curve2*)myC[1]), myPoints2(N)));
7fd59977 419}
7fd59977 420
5493d334 421//=======================================================================
422//function : Points
423//purpose :
424//=======================================================================
425void Extrema_GenExtCC::Points(const Standard_Integer N,
426 POnC& P1,
427 POnC& P2) const
7fd59977 428{
4bbaf12b 429 StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::Points()")
430 Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::Points()")
431
432 P1.SetValues(myPoints1(N), Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)));
433 P2.SetValues(myPoints2(N), Tool2::Value(*((Curve2*)myC[1]), myPoints2(N)));
20a216fe 434}
836d7b64 435
436//=======================================================================
437//function : SetSingleSolutionFlag
438//purpose :
439//=======================================================================
440void Extrema_GenExtCC::SetSingleSolutionFlag(const Standard_Boolean theFlag)
441{
442 myIsFindSingleSolution = theFlag;
443}
444
445//=======================================================================
446//function : GetSingleSolutionFlag
447//purpose :
448//=======================================================================
449Standard_Boolean Extrema_GenExtCC::GetSingleSolutionFlag() const
450{
451 return myIsFindSingleSolution;
452}