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 <TColStd_ListOfInteger.hxx>
26 #include <Precision.hxx>
27 #include <NCollection_Vector.hxx>
28 #include <NCollection_CellFilter.hxx>
29 #include <GCPnts_AbscissaPoint.hxx>
31 // Comparator, used in std::sort.
32 static Standard_Boolean comp(const gp_XY& theA,
35 if (theA.X() < theB.X())
41 if (theA.X() == theB.X())
43 if (theA.Y() < theB.Y())
48 return Standard_False;
51 static void ChangeIntervals(Handle(TColStd_HArray1OfReal)& theInts, const Standard_Integer theNbInts)
53 Standard_Integer aNbInts = theInts->Length() - 1;
54 Standard_Integer aNbAdd = theNbInts - aNbInts;
55 Handle(TColStd_HArray1OfReal) aNewInts = new TColStd_HArray1OfReal(1, theNbInts + 1);
56 Standard_Integer aNbLast = theInts->Length();
60 aNewInts->SetValue(1, theInts->First());
61 aNewInts->SetValue(theNbInts + 1, theInts->Last());
62 Standard_Real dt = (theInts->Last() - theInts->First()) / theNbInts;
63 Standard_Real t = theInts->First() + dt;
64 for (i = 2; i <= theNbInts; ++i, t += dt)
66 aNewInts->SetValue(i, t);
72 for (i = 1; i <= aNbLast; ++i)
74 aNewInts->SetValue(i, theInts->Value(i));
79 Standard_Real anLIntMax = -1.;
80 Standard_Integer aMaxInd = -1;
81 for (i = 1; i < aNbLast; ++i)
83 Standard_Real anL = aNewInts->Value(i + 1) - aNewInts->Value(i);
91 Standard_Real t = (aNewInts->Value(aMaxInd + 1) + aNewInts->Value(aMaxInd)) / 2.;
92 for (i = aNbLast; i > aMaxInd; --i)
94 aNewInts->SetValue(i + 1, aNewInts->Value(i));
98 aNewInts->SetValue(aMaxInd + 1, t);
102 class Extrema_CCPointsInspector : public NCollection_CellFilter_InspectorXY
105 typedef gp_XY Target;
106 //! Constructor; remembers the tolerance
107 Extrema_CCPointsInspector (const Standard_Real theTol)
109 myTol = theTol * theTol;
110 myIsFind = Standard_False;
115 myIsFind = Standard_False;
118 Standard_Boolean isFind()
123 //! Set current point to search for coincidence
124 void SetCurrent (const gp_XY& theCurPnt)
126 myCurrent = theCurPnt;
129 //! Implementation of inspection method
130 NCollection_CellFilter_Action Inspect (const Target& theObject)
132 gp_XY aPt = myCurrent.Subtracted(theObject);
133 const Standard_Real aSQDist = aPt.SquareModulus();
136 myIsFind = Standard_True;
139 return CellFilter_Keep;
145 Standard_Boolean myIsFind;
148 //=======================================================================
149 //function : ProjPOnC
150 //purpose : Projects the point on the curve and returns the minimal
151 // projection distance
152 //=======================================================================
153 static Standard_Real ProjPOnC(const Pnt& theP,
154 Extrema_GExtPC& theProjTool)
156 Standard_Real aDist = ::RealLast();
157 theProjTool.Perform(theP);
158 if (theProjTool.IsDone() && theProjTool.NbExt())
160 for (Standard_Integer i = 1; i <= theProjTool.NbExt(); ++i)
162 Standard_Real aD = theProjTool.SquareDistance(i);
171 //=======================================================================
172 //function : Extrema_GenExtCC
174 //=======================================================================
175 Extrema_GenExtCC::Extrema_GenExtCC()
176 : myIsFindSingleSolution(Standard_False),
177 myParallel(Standard_False),
178 myCurveMinTol(Precision::PConfusion()),
181 myDone(Standard_False)
186 //=======================================================================
187 //function : Extrema_GenExtCC
189 //=======================================================================
190 Extrema_GenExtCC::Extrema_GenExtCC(const Curve1& C1,
192 : myIsFindSingleSolution(Standard_False),
193 myParallel(Standard_False),
194 myCurveMinTol(Precision::PConfusion()),
197 myDone(Standard_False)
199 myC[0] = (Standard_Address)&C1;
200 myC[1] = (Standard_Address)&C2;
201 myLowBorder(1) = C1.FirstParameter();
202 myLowBorder(2) = C2.FirstParameter();
203 myUppBorder(1) = C1.LastParameter();
204 myUppBorder(2) = C2.LastParameter();
207 //=======================================================================
208 //function : Extrema_GenExtCC
210 //=======================================================================
211 Extrema_GenExtCC::Extrema_GenExtCC(const Curve1& C1,
213 const Standard_Real Uinf,
214 const Standard_Real Usup,
215 const Standard_Real Vinf,
216 const Standard_Real Vsup)
217 : myIsFindSingleSolution(Standard_False),
218 myParallel(Standard_False),
219 myCurveMinTol(Precision::PConfusion()),
222 myDone(Standard_False)
224 myC[0] = (Standard_Address)&C1;
225 myC[1] = (Standard_Address)&C2;
226 myLowBorder(1) = Uinf;
227 myLowBorder(2) = Vinf;
228 myUppBorder(1) = Usup;
229 myUppBorder(2) = Vsup;
232 //=======================================================================
233 //function : SetParams
235 //=======================================================================
236 void Extrema_GenExtCC::SetParams(const Curve1& C1,
238 const Standard_Real Uinf,
239 const Standard_Real Usup,
240 const Standard_Real Vinf,
241 const Standard_Real Vsup)
243 myC[0] = (Standard_Address)&C1;
244 myC[1] = (Standard_Address)&C2;
245 myLowBorder(1) = Uinf;
246 myLowBorder(2) = Vinf;
247 myUppBorder(1) = Usup;
248 myUppBorder(2) = Vsup;
251 //=======================================================================
252 //function : SetTolerance
254 //=======================================================================
255 void Extrema_GenExtCC::SetTolerance(Standard_Real theTol)
257 myCurveMinTol = theTol;
260 //=======================================================================
263 //=======================================================================
264 void Extrema_GenExtCC::Perform()
266 myDone = Standard_False;
267 myParallel = Standard_False;
269 Curve1 &C1 = *(Curve1*)myC[0];
270 Curve2 &C2 = *(Curve2*)myC[1];
272 Standard_Integer aNbInter[2];
273 GeomAbs_Shape aContinuity = GeomAbs_C2;
274 aNbInter[0] = C1.NbIntervals(aContinuity);
275 aNbInter[1] = C2.NbIntervals(aContinuity);
277 if (aNbInter[0] * aNbInter[1] > 100)
279 aContinuity = GeomAbs_C1;
280 aNbInter[0] = C1.NbIntervals(aContinuity);
281 aNbInter[1] = C2.NbIntervals(aContinuity);
284 Standard_Real anL[2];
285 Standard_Integer indmax = -1, indmin = -1;
286 const Standard_Real mult = 20.;
287 if (!(Precision::IsInfinite(C1.FirstParameter()) || Precision::IsInfinite(C1.LastParameter()) ||
288 Precision::IsInfinite(C2.FirstParameter()) || Precision::IsInfinite(C2.LastParameter())))
290 anL[0] = GCPnts_AbscissaPoint::Length(C1);
291 anL[1] = GCPnts_AbscissaPoint::Length(C2);
292 if (anL[0] / aNbInter[0] > mult * anL[1] / aNbInter[1])
297 else if (anL[1] / aNbInter[1] > mult * anL[0] / aNbInter[0])
303 Standard_Integer aNbIntOpt = 0;
306 aNbIntOpt = RealToInt(anL[indmax] * aNbInter[indmin] / anL[indmin] / (mult / 4.)) + 1;
307 if (aNbIntOpt > 100 || aNbIntOpt < aNbInter[indmax])
313 if (aNbIntOpt * aNbInter[indmin] > 100)
315 aNbIntOpt = 100 / aNbInter[indmin];
316 if (aNbIntOpt < aNbInter[indmax])
324 Handle(TColStd_HArray1OfReal) anIntervals1 = new TColStd_HArray1OfReal(1, aNbInter[0] + 1);
325 Handle(TColStd_HArray1OfReal) anIntervals2 = new TColStd_HArray1OfReal(1, aNbInter[1] + 1);
326 C1.Intervals(anIntervals1->ChangeArray1(), aContinuity);
327 C2.Intervals(anIntervals2->ChangeArray1(), aContinuity);
332 //Change anIntervals1
333 ChangeIntervals(anIntervals1, aNbIntOpt);
334 aNbInter[0] = anIntervals1->Length() - 1;
338 //Change anIntervals2;
339 ChangeIntervals(anIntervals2, aNbIntOpt);
340 aNbInter[1] = anIntervals2->Length() - 1;
344 // Lipchitz constant computation.
345 const Standard_Real aMaxLC = 10000.;
346 Standard_Real aLC = 9.0; // Default value.
347 const Standard_Real aMaxDer1 = 1.0 / C1.Resolution(1.0);
348 const Standard_Real aMaxDer2 = 1.0 / C2.Resolution(1.0);
349 Standard_Real aMaxDer = Max(aMaxDer1, aMaxDer2) * Sqrt(2.0);
353 // Change constant value according to the concrete curve types.
354 Standard_Boolean isConstLockedFlag = Standard_False;
355 //To prevent LipConst to became too small
356 const Standard_Real aCR = 0.001;
357 if (aMaxDer1 / aMaxDer < aCR || aMaxDer2 / aMaxDer < aCR)
359 isConstLockedFlag = Standard_True;
361 if (aMaxDer > aMaxLC)
364 isConstLockedFlag = Standard_True;
366 if (C1.GetType() == GeomAbs_Line)
368 aMaxDer = 1.0 / C2.Resolution(1.0);
371 isConstLockedFlag = Standard_True;
375 if (C2.GetType() == GeomAbs_Line)
377 aMaxDer = 1.0 / C1.Resolution(1.0);
380 isConstLockedFlag = Standard_True;
385 Extrema_GlobOptFuncCCC2 aFunc (C1, C2);
386 math_GlobOptMin aFinder(&aFunc, myLowBorder, myUppBorder, aLC);
387 aFinder.SetLipConstState(isConstLockedFlag);
388 aFinder.SetContinuity(aContinuity == GeomAbs_C2 ? 2 : 1);
389 Standard_Real aDiscTol = 1.0e-2;
390 Standard_Real aValueTol = 1.0e-2;
391 Standard_Real aSameTol = myCurveMinTol / (aDiscTol);
392 aFinder.SetTol(aDiscTol, aSameTol);
393 aFinder.SetFunctionalMinimalValue(0.0); // Best distance cannot be lower than 0.0.
395 // Size computed to have cell index inside of int32 value.
396 const Standard_Real aCellSize = Max(Max(anIntervals1->Last() - anIntervals1->First(),
397 anIntervals2->Last() - anIntervals2->First())
398 * Precision::PConfusion() / (2.0 * Sqrt(2.0)),
399 Precision::PConfusion());
400 Extrema_CCPointsInspector anInspector(aCellSize);
401 NCollection_CellFilter<Extrema_CCPointsInspector> aFilter(aCellSize);
402 NCollection_Vector<gp_XY> aPnts;
404 Standard_Integer i,j,k;
405 math_Vector aFirstBorderInterval(1,2);
406 math_Vector aSecondBorderInterval(1,2);
407 Standard_Real aF = RealLast(); // Best functional value.
408 Standard_Real aCurrF = RealLast(); // Current functional value computed on current interval.
409 for(i = 1; i <= aNbInter[0]; i++)
411 for(j = 1; j <= aNbInter[1]; j++)
413 aFirstBorderInterval(1) = anIntervals1->Value(i);
414 aFirstBorderInterval(2) = anIntervals2->Value(j);
415 aSecondBorderInterval(1) = anIntervals1->Value(i + 1);
416 aSecondBorderInterval(2) = anIntervals2->Value(j + 1);
418 aFinder.SetLocalParams(aFirstBorderInterval, aSecondBorderInterval);
419 aFinder.Perform(GetSingleSolutionFlag());
421 // Check that solution found on current interval is not worse than previous.
422 aCurrF = aFinder.GetF();
423 if (aCurrF >= aF + aSameTol * aValueTol)
428 // Clean previously computed solution if current one is better.
429 if (aCurrF > aF - aSameTol * aValueTol)
437 aFilter.Reset(aCellSize);
441 // Save found solutions avoiding repetitions.
442 math_Vector sol(1,2);
443 for(k = 1; k <= aFinder.NbExtrema(); k++)
445 aFinder.Points(k, sol);
446 gp_XY aPnt2d(sol(1), sol(2));
448 gp_XY aXYmin = anInspector.Shift(aPnt2d, -aCellSize);
449 gp_XY aXYmax = anInspector.Shift(aPnt2d, aCellSize);
451 anInspector.ClearFind();
452 anInspector.SetCurrent(aPnt2d);
453 aFilter.Inspect(aXYmin, aXYmax, anInspector);
454 if (!anInspector.isFind())
456 // Point is out of close cells, add new one.
457 aFilter.Add(aPnt2d, aPnt2d);
458 aPnts.Append(gp_XY(sol(1), sol(2)));
464 const Standard_Integer aNbSol = aPnts.Length();
468 myDone = Standard_False;
472 myDone = Standard_True;
477 const gp_XY& aSol = aPnts.First();
478 myPoints1.Append(aSol.X());
479 myPoints2.Append(aSol.Y());
483 // More than one solution is found.
484 // Check for infinity solutions case, for this:
485 // Sort points lexicographically and check midpoint between each two neighboring points.
486 // If all midpoints functional value is acceptable then check the projection distances
487 // of the bounding points of the curves onto the opposite curves.
488 // If these distances are also acceptable set myParallel flag to true and return one solution.
489 std::sort(aPnts.begin(), aPnts.end(), comp);
491 // Solutions to pass into result.
492 // If the parallel segment is found, save only extreme solutions on that segment.
493 // The first and last solutions will always be the extreme ones, thus save them unconditionally.
494 TColStd_ListOfInteger aSolutions;
496 // Manages the addition of the solution into result.
497 // Set it to TRUE to add the first solution.
498 Standard_Boolean bSaveSolution = Standard_True;
500 // Define direction of the second curve relatively the first one
501 // (it will be needed for projection).
502 Standard_Boolean bDirsCoinside = Standard_True;
503 // Check also if the found solutions are not concentrated in one point
504 // on any of the curves. And if they are, avoid marking the curves as parallel.
505 Standard_Boolean bDifferentSolutions = Standard_False;
507 Standard_Boolean isParallel = Standard_True;
508 Standard_Real aVal = 0.0;
509 math_Vector aVec(1, 2, 0.0);
511 // Iterate on all solutions and collect the extreme solutions on all parallel segments.
512 for (Standard_Integer anIdx = 0; anIdx < aNbSol - 1; anIdx++)
514 const gp_XY& aCurrent = aPnts(anIdx);
515 const gp_XY& aNext = aPnts(anIdx + 1);
517 aVec(1) = (aCurrent.X() + aNext.X()) * 0.5;
518 aVec(2) = (aCurrent.Y() + aNext.Y()) * 0.5;
520 aFunc.Value(aVec, aVal);
522 if (Abs(aVal - aF) < Precision::Confusion())
524 // It seems the parallel segment is found.
525 // Save only extreme solutions on that segment.
528 // Add current solution as the beginning of the parallel segment.
529 aSolutions.Append(anIdx);
530 // Do not keep the next solution in current parallel segment.
531 bSaveSolution = Standard_False;
536 // Mid point does not satisfy the tolerance criteria, curves are not parallel.
537 isParallel = Standard_False;
538 // Add current solution as the last one in previous parallel segment.
539 aSolutions.Append(anIdx);
540 // Save also the next solution as the first one in next parallel segment.
541 bSaveSolution = Standard_True;
544 if (!bDifferentSolutions)
546 if (aNext.X() > aCurrent.X())
548 if (aNext.Y() > aCurrent.Y())
550 bDifferentSolutions = Standard_True;
551 bDirsCoinside = Standard_True;
553 else if (aNext.Y() < aCurrent.Y())
555 bDifferentSolutions = Standard_True;
556 bDirsCoinside = Standard_False;
561 // Save the last solution
562 aSolutions.Append(aNbSol - 1);
564 if (!bDifferentSolutions)
565 isParallel = Standard_False;
569 // For the check on parallel case it is also necessary to check additionally
570 // if the ends of the curves do not diverge. For this, project the bounding
571 // points of the curves on the opposite curves and check the distances.
573 Standard_Real aT1[2] = {myLowBorder(1), myUppBorder(1)};
574 Standard_Real aT2[2] = {bDirsCoinside ? myLowBorder(2) : myUppBorder(2),
575 bDirsCoinside ? myUppBorder(2) : myLowBorder(2)};
577 Extrema_GExtPC anExtPC1, anExtPC2;
578 anExtPC1.Initialize(C1, myLowBorder(1), myUppBorder(1));
579 anExtPC2.Initialize(C2, myLowBorder(2), myUppBorder(2));
581 for (Standard_Integer iT = 0; isParallel && (iT < 2); ++iT)
583 Standard_Real aDist1 = ProjPOnC(C1.Value(aT1[iT]), anExtPC2);
584 Standard_Real aDist2 = ProjPOnC(C2.Value(aT2[iT]), anExtPC1);
585 isParallel = (Abs(Min(aDist1, aDist2) - aF) < Precision::Confusion());
591 // Keep only one solution
592 const gp_XY& aSol = aPnts.First();
593 myPoints1.Append(aSol.X());
594 myPoints2.Append(aSol.Y());
595 myParallel = Standard_True;
599 // Keep all saved solutions
600 TColStd_ListIteratorOfListOfInteger aItSol(aSolutions);
601 for (; aItSol.More(); aItSol.Next())
603 const gp_XY& aSol = aPnts(aItSol.Value());
604 myPoints1.Append(aSol.X());
605 myPoints2.Append(aSol.Y());
610 //=======================================================================
613 //=======================================================================
614 Standard_Boolean Extrema_GenExtCC::IsDone() const
619 //=======================================================================
620 //function : IsParallel
622 //=======================================================================
623 Standard_Boolean Extrema_GenExtCC::IsParallel() const
625 if (!IsDone()) throw StdFail_NotDone();
629 //=======================================================================
632 //=======================================================================
633 Standard_Integer Extrema_GenExtCC::NbExt() const
635 if (!IsDone()) throw StdFail_NotDone();
636 return myPoints1.Length();
639 //=======================================================================
640 //function : SquareDistance
642 //=======================================================================
643 Standard_Real Extrema_GenExtCC::SquareDistance(const Standard_Integer N) const
645 if (N < 1 || N > NbExt())
647 throw Standard_OutOfRange();
650 return Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)).SquareDistance(Tool2::Value(*((Curve2*)myC[1]), myPoints2(N)));
653 //=======================================================================
656 //=======================================================================
657 void Extrema_GenExtCC::Points(const Standard_Integer N,
663 throw StdFail_InfiniteSolutions();
666 if (N < 1 || N > NbExt())
668 throw Standard_OutOfRange();
671 P1.SetValues(myPoints1(N), Tool1::Value(*((Curve1*)myC[0]), myPoints1(N)));
672 P2.SetValues(myPoints2(N), Tool2::Value(*((Curve2*)myC[1]), myPoints2(N)));
675 //=======================================================================
676 //function : SetSingleSolutionFlag
678 //=======================================================================
679 void Extrema_GenExtCC::SetSingleSolutionFlag(const Standard_Boolean theFlag)
681 myIsFindSingleSolution = theFlag;
684 //=======================================================================
685 //function : GetSingleSolutionFlag
687 //=======================================================================
688 Standard_Boolean Extrema_GenExtCC::GetSingleSolutionFlag() const
690 return myIsFindSingleSolution;