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
6 // This file is part of Open CASCADE Technology software library.
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.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
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>
29 // Comparator, used in std::sort.
30 static Standard_Boolean comp(const gp_XY& theA,
33 if (theA.X() < theB.X())
39 if (theA.X() == theB.X())
41 if (theA.Y() <= theB.Y())
46 return Standard_False;
49 class Extrema_CCPointsInspector : public NCollection_CellFilter_InspectorXY
53 //! Constructor; remembers the tolerance
54 Extrema_CCPointsInspector (const Standard_Real theTol)
56 myTol = theTol * theTol;
57 myIsFind = Standard_False;
62 myIsFind = Standard_False;
65 Standard_Boolean isFind()
70 //! Set current point to search for coincidence
71 void SetCurrent (const gp_XY& theCurPnt)
73 myCurrent = theCurPnt;
76 //! Implementation of inspection method
77 NCollection_CellFilter_Action Inspect (const Target& theObject)
79 gp_XY aPt = myCurrent.Subtracted(theObject);
80 const Standard_Real aSQDist = aPt.SquareModulus();
83 myIsFind = Standard_True;
86 return CellFilter_Keep;
92 Standard_Boolean myIsFind;
95 //=======================================================================
96 //function : Extrema_GenExtCC
98 //=======================================================================
99 Extrema_GenExtCC::Extrema_GenExtCC()
100 : myParallel(Standard_False),
101 myCurveMinTol(Precision::PConfusion()),
104 myDone(Standard_False)
109 //=======================================================================
110 //function : Extrema_GenExtCC
112 //=======================================================================
113 Extrema_GenExtCC::Extrema_GenExtCC(const Curve1& C1,
115 : myParallel(Standard_False),
116 myCurveMinTol(Precision::PConfusion()),
119 myDone(Standard_False)
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();
129 //=======================================================================
130 //function : Extrema_GenExtCC
132 //=======================================================================
133 Extrema_GenExtCC::Extrema_GenExtCC(const Curve1& C1,
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()),
143 myDone(Standard_False)
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;
153 //=======================================================================
154 //function : SetParams
156 //=======================================================================
157 void Extrema_GenExtCC::SetParams(const Curve1& C1,
159 const Standard_Real Uinf,
160 const Standard_Real Usup,
161 const Standard_Real Vinf,
162 const Standard_Real Vsup)
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;
172 //=======================================================================
173 //function : SetTolerance
175 //=======================================================================
176 void Extrema_GenExtCC::SetTolerance(Standard_Real theTol)
178 myCurveMinTol = theTol;
181 //=======================================================================
184 //=======================================================================
185 void Extrema_GenExtCC::Perform()
187 myDone = Standard_False;
188 myParallel = Standard_False;
190 Curve1 &C1 = *(Curve1*)myC[0];
191 Curve2 &C2 = *(Curve2*)myC[1];
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);
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);
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(Extrema_CCPointsInspector::Dimension, aCellSize);
214 NCollection_Vector<gp_XY> aPnts;
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++)
223 for(j = 1; j <= aNbInter[1]; j++)
225 aFirstBorderInterval(1) = anIntervals1(i);
226 aFirstBorderInterval(2) = anIntervals2(j);
227 aSecondBorderInterval(1) = anIntervals1(i + 1);
228 aSecondBorderInterval(2) = anIntervals2(j + 1);
230 aFinder.SetLocalParams(aFirstBorderInterval, aSecondBorderInterval);
233 // Check that solution found on current interval is not worse than previous.
234 aCurrF = aFinder.GetF();
235 if (aCurrF >= aF + aSameTol * aValueTol)
240 // Clean previously computed solution if current one is better.
241 if (aCurrF > aF - aSameTol * aValueTol)
249 aFilter.Reset(aCellSize);
253 // Save found solutions avoiding repetitions.
254 math_Vector sol(1,2);
255 for(k = 1; k <= aFinder.NbExtrema(); k++)
257 aFinder.Points(k, sol);
258 gp_XY aPnt2d(sol(1), sol(2));
260 gp_XY aXYmin = anInspector.Shift(aPnt2d, -aCellSize);
261 gp_XY aXYmax = anInspector.Shift(aPnt2d, aCellSize);
263 anInspector.ClearFind();
264 anInspector.SetCurrent(aPnt2d);
265 aFilter.Inspect(aXYmin, aXYmax, anInspector);
266 if (!anInspector.isFind())
268 // Point is out of close cells, add new one.
269 aFilter.Add(aPnt2d, aPnt2d);
270 aPnts.Append(gp_XY(sol(1), sol(2)));
276 if (aPnts.Size() == 0)
279 myDone = Standard_False;
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);
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)
297 isParallel = Standard_True;
298 for(Standard_Integer anIdx = aPnts.Lower(); anIdx <= aPnts.Upper() - 1; anIdx++)
300 const gp_XY& aCurrent = aPnts(anIdx);
301 const gp_XY& aNext = aPnts(anIdx + 1);
303 aVec(1) = (aCurrent.X() + aNext.X()) * 0.5;
304 aVec(2) = (aCurrent.Y() + aNext.Y()) * 0.5;
306 aFunc.Value(aVec, aVal);
308 if (Abs(aVal - aF) > Precision::Confusion())
310 isParallel = Standard_False;
318 const gp_XY& aCurrent = aPnts.First();
319 myPoints1.Append(aCurrent.X());
320 myPoints2.Append(aCurrent.Y());
321 myParallel = Standard_True;
325 for(Standard_Integer anIdx = aPnts.Lower(); anIdx <= aPnts.Upper(); anIdx++)
327 const gp_XY& aCurrent = aPnts(anIdx);
328 myPoints1.Append(aCurrent.X());
329 myPoints2.Append(aCurrent.Y());
333 myDone = Standard_True;
336 //=======================================================================
339 //=======================================================================
340 Standard_Boolean Extrema_GenExtCC::IsDone() const
345 //=======================================================================
346 //function : IsParallel
348 //=======================================================================
349 Standard_Boolean Extrema_GenExtCC::IsParallel() const
354 //=======================================================================
357 //=======================================================================
358 Standard_Integer Extrema_GenExtCC::NbExt() const
360 StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::NbExt()")
362 return myPoints1.Length();
365 //=======================================================================
366 //function : SquareDistance
368 //=======================================================================
369 Standard_Real Extrema_GenExtCC::SquareDistance(const Standard_Integer N) const
371 StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::SquareDistance()")
372 Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::SquareDistance()")
374 return Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)).SquareDistance(Tool2::Value(*((Curve2*)myC[1]), myPoints2(N)));
377 //=======================================================================
380 //=======================================================================
381 void Extrema_GenExtCC::Points(const Standard_Integer N,
385 StdFail_NotDone_Raise_if (!myDone, "Extrema_GenExtCC::Points()")
386 Standard_OutOfRange_Raise_if ((N < 1 || N > NbExt()), "Extrema_GenExtCC::Points()")
388 P1.SetValues(myPoints1(N), Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)));
389 P2.SetValues(myPoints2(N), Tool2::Value(*((Curve2*)myC[1]), myPoints2(N)));