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