0028346: Function ProjectOnSegments of ShapeAnalysis_Curve returns only single soluti...
[occt.git] / src / Extrema / Extrema_GExtPC.gxx
1 // Created on: 1992-10-19
2 // Created by: Laurent PAINNOT
3 // Copyright (c) 1992-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 // Modified by cma, Fri Dec 10 18:11:56 1993
18
19 #include Extrema_EPC_hxx
20 #include ThePOnC_hxx
21 #include TheSequenceOfPOnC_hxx
22 #include ThePoint_hxx
23 #include TheVector_hxx
24 #include TheExtPElC_hxx
25 #include <StdFail_NotDone.hxx>
26 #include <Standard_Failure.hxx>
27 #include <GeomAbs_CurveType.hxx>
28 #include <Precision.hxx>
29 #include <ElCLib.hxx>
30 #include <TColStd_Array1OfReal.hxx>
31 #include <TColStd_HArray1OfReal.hxx>
32 #include <NCollection_Array1.hxx>
33
34
35 //=======================================================================
36 //function : Perform
37 //purpose  : 
38 //=======================================================================
39 void Extrema_GExtPC::Perform(const ThePoint& P)
40 {
41   mySqDist.Clear();
42   mypoint.Clear();
43   myismin.Clear();
44   Standard_Integer i, NbExt, n;                     
45   Standard_Real U;
46   mysample = 17;
47   Standard_Real t3d = Precision::Confusion();
48   if (Precision::IsInfinite(myuinf)) mydist1 = RealLast();
49   else {
50     Pf = TheCurveTool::Value(*((TheCurve*)myC), myuinf);
51     mydist1 = P.SquareDistance(Pf);
52   }
53   
54   if (Precision::IsInfinite(myusup)) mydist2 = RealLast();
55   else {
56     Pl = TheCurveTool::Value(*((TheCurve*)myC), myusup);
57     mydist2 = P.SquareDistance(Pl);
58   }
59
60   TheCurve & aCurve = *((TheCurve*)myC);
61
62   switch(type) {
63   case GeomAbs_Circle: 
64     {
65       myExtPElC.Perform(P, TheCurveTool::Circle(aCurve), t3d, myuinf, myusup);
66     }
67     break;
68   case GeomAbs_Ellipse: 
69     {
70       myExtPElC.Perform(P, TheCurveTool::Ellipse(aCurve), t3d, myuinf, myusup);
71     }
72     break;
73   case GeomAbs_Parabola: 
74     {
75       myExtPElC.Perform(P, TheCurveTool::Parabola(aCurve), t3d,myuinf,myusup);
76     }
77     break;
78   case GeomAbs_Hyperbola: 
79     {
80       myExtPElC.Perform(P,TheCurveTool::Hyperbola(aCurve),t3d, myuinf, myusup);
81     }
82     break;
83   case GeomAbs_Line: 
84     {
85       myExtPElC.Perform(P, TheCurveTool::Line(aCurve), t3d, myuinf, myusup);
86     }
87     break;
88   case GeomAbs_BezierCurve:
89     {
90       myintuinf = myuinf;
91       myintusup = myusup;
92       mysample = (TheCurveTool::Bezier(aCurve))->NbPoles() * 2;
93       myExtPC.Initialize(aCurve);
94       IntervalPerform(P);
95       return;
96     }
97   case GeomAbs_BSplineCurve: 
98     {
99       const Standard_Integer 
100         aFirstIdx = TheCurveTool::BSpline(aCurve)->FirstUKnotIndex(),
101         aLastIdx  = TheCurveTool::BSpline(aCurve)->LastUKnotIndex();
102       // const reference can not be used due to implementation of BRep_Adaptor.
103       TColStd_Array1OfReal aKnots(aFirstIdx, aLastIdx);
104       TheCurveTool::BSpline(aCurve)->Knots(aKnots); 
105
106       // Workaround to work with:
107       // blend, where knots may be moved from parameter space.
108       Standard_Real aPeriodJump = 0.0;
109       // Avoid problem with too close knots.
110       const Standard_Real aTolCoeff = (myusup - myuinf) * Precision::PConfusion();
111       if (TheCurveTool::IsPeriodic(aCurve))
112       {
113         Standard_Integer aPeriodShift =
114           Standard_Integer ((myuinf - aKnots(aFirstIdx)) / TheCurveTool::Period(aCurve));
115         
116         if (myuinf < aKnots(aFirstIdx) - aTolCoeff)
117           aPeriodShift--;
118
119         aPeriodJump = TheCurveTool::Period(aCurve) * aPeriodShift;
120       }
121
122       Standard_Integer anIdx;
123
124       // Find first and last used knot.
125       Standard_Integer aFirstUsedKnot = aFirstIdx,
126                         aLastUsedKnot = aLastIdx;
127       for(anIdx = aFirstIdx; anIdx <= aLastIdx; anIdx++)
128       {
129         Standard_Real aKnot = aKnots(anIdx) + aPeriodJump;
130         if (myuinf >= aKnot - aTolCoeff)
131           aFirstUsedKnot = anIdx;
132         else
133           break;
134
135       }
136       for(anIdx = aLastIdx; anIdx >= aFirstIdx; anIdx--)
137       {
138         Standard_Real aKnot = aKnots(anIdx) + aPeriodJump;
139         if (myusup <= aKnot + aTolCoeff)
140           aLastUsedKnot = anIdx;
141         else
142           break;
143       }
144
145       if (aFirstUsedKnot == aLastUsedKnot)
146       {
147         // Degenerated case:
148         // Some bounds lies out of curve param space.
149         // In this case build one interval with [myuinf, myusup].
150         // Parameters of these indexes will be redefined.
151         aFirstUsedKnot = aFirstIdx;
152         aLastUsedKnot  = aFirstIdx + 1;
153       }
154
155       mysample = (TheCurveTool::BSpline(aCurve))->Degree() + 1;
156
157       // Fill sample points.
158       Standard_Integer aValIdx = 1;
159       NCollection_Array1<Standard_Real> aVal(1, (mysample) * (aLastUsedKnot - aFirstUsedKnot) + 1);
160       NCollection_Array1<Standard_Real> aParam(1, (mysample) * (aLastUsedKnot - aFirstUsedKnot) + 1);
161       for(anIdx = aFirstUsedKnot; anIdx < aLastUsedKnot; anIdx++)
162       {
163         Standard_Real aF = aKnots(anIdx) + aPeriodJump,
164                       aL = aKnots(anIdx + 1) + aPeriodJump;
165
166         if (anIdx == aFirstUsedKnot)
167           aF = myuinf;
168         if (anIdx == aLastUsedKnot - 1)
169           aL = myusup;
170
171         Standard_Real aStep = (aL - aF) / mysample;
172         for(Standard_Integer aPntIdx = 0; aPntIdx < mysample; aPntIdx++)
173         {
174           Standard_Real aCurrentParam = aF + aStep * aPntIdx;
175           aVal(aValIdx) = TheCurveTool::Value(aCurve, aCurrentParam).SquareDistance(P);
176           aParam(aValIdx) = aCurrentParam;
177           aValIdx++;
178         }
179       }
180       // Fill last point.
181       aVal(aValIdx) = TheCurveTool::Value(aCurve, myusup).SquareDistance(P);
182       aParam(aValIdx) = myusup;
183
184        myExtPC.Initialize(aCurve);
185
186       // Find extremas.
187       for(anIdx = aVal.Lower() + 1; anIdx < aVal.Upper(); anIdx++)
188       {
189         if (aVal(anIdx) <= Precision::SquareConfusion())
190         {
191           mySqDist.Append(aVal(anIdx));
192           myismin.Append(Standard_True);
193           mypoint.Append(ThePOnC(aParam(anIdx), TheCurveTool::Value(aCurve, aParam(anIdx))));
194         }
195         if ((aVal(anIdx) >= aVal(anIdx + 1) &&
196              aVal(anIdx) >= aVal(anIdx - 1)) ||
197             (aVal(anIdx) <= aVal(anIdx + 1) &&
198              aVal(anIdx) <= aVal(anIdx - 1)) )
199         {
200           myintuinf = aParam(anIdx - 1);
201           myintusup = aParam(anIdx + 1);
202
203           IntervalPerform(P);
204         }
205       }
206
207       // Solve on the first and last intervals.
208       if (mydist1 > Precision::SquareConfusion() && !Precision::IsPositiveInfinite(mydist1))
209       {
210         ThePoint aP1, aP2;
211         TheVector aV1, aV2;
212         TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower()),     aP1, aV1);
213         TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower() + 1), aP2, aV2);
214         TheVector aBase1(P, aP1), aBase2(P, aP2);
215         Standard_Real aVal1 = aV1.Dot(aBase1); // Derivative of (C(u) - P)^2
216         Standard_Real aVal2 = aV2.Dot(aBase2); // Derivative of (C(u) - P)^2
217         if(!(Precision::IsInfinite(aVal1) || Precision::IsInfinite(aVal2)))
218         {
219           // Derivatives have opposite signs - min or max inside of interval (sufficient condition).
220           // Necessary condition - when point lies on curve.
221           // Necessary condition - when derivative of point is too small.
222           if(aVal1 * aVal2 <= 0.0 ||
223             aBase1.Dot(aBase2) <= 0.0 ||
224             2.0 * Abs(aVal1) < Precision::Confusion() )
225           {
226             myintuinf = aParam(aVal.Lower());
227             myintusup = aParam(aVal.Lower() + 1);
228             IntervalPerform(P);
229           }
230         }
231       }
232
233       if (mydist2 > Precision::SquareConfusion() && !Precision::IsPositiveInfinite(mydist2))
234       {
235         ThePoint aP1, aP2;
236         TheVector aV1, aV2;
237         TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper() - 1), aP1, aV1);
238         TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper()),     aP2, aV2);
239         TheVector aBase1(P, aP1), aBase2(P, aP2);
240         Standard_Real aVal1 = aV1.Dot(aBase1); // Derivative of (C(u) - P)^2
241         Standard_Real aVal2 = aV2.Dot(aBase2); // Derivative of (C(u) - P)^2
242
243         if(!(Precision::IsInfinite(aVal1) || Precision::IsInfinite(aVal2)))
244         {
245           // Derivatives have opposite signs - min or max inside of interval (sufficient condition).
246           // Necessary condition - when point lies on curve.
247           // Necessary condition - when derivative of point is too small.
248           if(aVal1 * aVal2 <= 0.0 ||
249             aBase1.Dot(aBase2) <= 0.0 ||
250             2.0 * Abs(aVal2) < Precision::Confusion() )
251           {
252             myintuinf = aParam(aVal.Upper() - 1);
253             myintusup = aParam(aVal.Upper());
254             IntervalPerform(P);
255           }
256         }
257       }
258
259       mydone = Standard_True;
260       break;
261     }
262   default:
263     {
264       const Standard_Integer aMaxSample = 17;
265       Standard_Boolean IntExtIsDone = Standard_False;
266       Standard_Boolean IntIsNotValid;
267       Handle(TColStd_HArray1OfReal) theHInter;
268       n = TheCurveTool::NbIntervals(aCurve, GeomAbs_C2);
269       if (n > 1)
270       {
271         theHInter = new TColStd_HArray1OfReal(1, n + 1);
272         TheCurveTool::Intervals(aCurve, theHInter->ChangeArray1(), GeomAbs_C2);
273       }
274       else
275       {
276         theHInter = TheCurveTool::DeflCurvIntervals(aCurve);
277         n = theHInter->Length() - 1;
278       }
279       mysample = Max(mysample / n, aMaxSample);
280       Standard_Real maxint = 0.;
281       for (i = 1; i <= n; ++i)
282       {
283         Standard_Real dt = theHInter->Value(i + 1) - theHInter->Value(i);
284         if (maxint < dt)
285         {
286           maxint = dt;
287         }
288       }
289       Standard_Boolean isPeriodic = TheCurveTool::IsPeriodic(aCurve);
290       TheVector V1;
291       ThePoint PP;
292       Standard_Real s1 = 0.0 ;
293       Standard_Real s2 = 0.0;
294       myExtPC.Initialize(aCurve);
295       for (i = 1; i <= n; i++)
296       {
297         myintuinf = theHInter->Value(i);
298         myintusup = theHInter->Value(i+1);
299         mysample = Max(RealToInt(aMaxSample*(myintusup - myintuinf) / maxint), 3);
300         
301         Standard_Real anInfToCheck = myintuinf;
302         Standard_Real aSupToCheck = myintusup;
303         
304         if (isPeriodic) {
305           Standard_Real aPeriod = TheCurveTool::Period(aCurve);
306           anInfToCheck = ElCLib::InPeriod(myintuinf, myuinf, myuinf+aPeriod);
307           aSupToCheck = myintusup+(anInfToCheck-myintuinf);
308         }    
309         IntIsNotValid = (myuinf > aSupToCheck) || (myusup < anInfToCheck);
310
311         if(IntIsNotValid) continue;
312
313         if (myuinf >= anInfToCheck) anInfToCheck = myuinf;
314         if (myusup <= aSupToCheck) aSupToCheck = myusup;
315         if((aSupToCheck - anInfToCheck) <= mytolu) continue;
316         
317         if (i != 1)
318           {
319           TheCurveTool::D1(aCurve, myintuinf, PP, V1);
320           s1 = (TheVector(P, PP))*V1;
321           if (s1*s2 < 0.0) {
322             mySqDist.Append(PP.SquareDistance(P));
323             myismin.Append((s1 < 0.0));
324             mypoint.Append(ThePOnC(myintuinf, PP));
325           }
326         }
327         if (i != n) {
328           TheCurveTool::D1(aCurve, myintusup, PP, V1);
329           s2 = (TheVector(P, PP))*V1;
330         }
331
332         IntervalPerform(P);
333         IntExtIsDone = IntExtIsDone || mydone;
334       }
335
336       mydone = IntExtIsDone;
337       break;
338     }
339   }
340
341   // Postprocessing.
342   if (type == GeomAbs_BSplineCurve ||
343       type == GeomAbs_OffsetCurve ||
344       type == GeomAbs_OtherCurve)
345   {
346     // Additional checking if the point is on the first or last point of the curve
347     // and does not added yet.
348     if (mydist1 < Precision::SquareConfusion() || 
349         mydist2 < Precision::SquareConfusion())
350     {
351       Standard_Boolean isFirstAdded = Standard_False;
352       Standard_Boolean isLastAdded  = Standard_False;
353       Standard_Integer aNbPoints = mypoint.Length();
354       for (i = 1; i <= aNbPoints; i++)
355       {
356         U = mypoint.Value(i).Parameter();
357         if (Abs(U - myuinf) < mytolu)
358           isFirstAdded = Standard_True;
359         else if (Abs(myusup - U) < mytolu)
360           isLastAdded = Standard_True;
361       }
362       if (!isFirstAdded && mydist1 < Precision::SquareConfusion())
363       {
364         mySqDist.Prepend(mydist1);
365         myismin.Prepend(Standard_True);
366         mypoint.Prepend(ThePOnC(myuinf, Pf));
367       }
368       if (!isLastAdded && mydist2 < Precision::SquareConfusion())
369       {
370         mySqDist.Append(mydist2);
371         myismin.Append(Standard_True);
372         mypoint.Append(ThePOnC(myusup, Pl));
373       }
374       mydone = Standard_True;
375     }
376   }
377   else
378   {
379     // In analytical case
380     mydone = myExtPElC.IsDone();
381     if (mydone)
382     {
383       NbExt = myExtPElC.NbExt();
384       for (i = 1; i <= NbExt; i++)
385       {
386         // Verification de la validite des parametres:
387         ThePOnC PC = myExtPElC.Point(i);
388         U = PC.Parameter();
389         if (TheCurveTool::IsPeriodic(aCurve))
390         {
391           U = ElCLib::InPeriod(U, myuinf, myuinf+TheCurveTool::Period(aCurve));
392         }
393         if ((U >= myuinf-mytolu) && (U <= myusup+mytolu))
394         {
395           PC.SetValues(U, myExtPElC.Point(i).Value());
396           mySqDist.Append(myExtPElC.SquareDistance(i));
397           myismin.Append(myExtPElC.IsMin(i));
398           mypoint.Append(PC);
399         }
400       }
401     }
402   }
403 }
404
405
406 //=======================================================================
407 //function : Initialize
408 //purpose  : 
409 //=======================================================================
410
411 void Extrema_GExtPC::Initialize(const TheCurve&     C,
412                                 const Standard_Real Uinf,
413                                 const Standard_Real Usup,
414                                 const Standard_Real TolF) 
415 {
416   myC = (Standard_Address)&C;
417   myintuinf = myuinf = Uinf;
418   myintusup = myusup = Usup;
419   mytolf = TolF;
420   mytolu = TheCurveTool::Resolution(*((TheCurve*)myC), Precision::Confusion());
421   type = TheCurveTool::GetType(C);
422   mydone = Standard_False;
423   mydist1 = RealLast();
424   mydist2 = RealLast();
425   mysample = 17;
426 }
427
428
429 //=======================================================================
430 //function : IntervalPerform
431 //purpose  : 
432 //=======================================================================
433
434 void Extrema_GExtPC::IntervalPerform(const ThePoint& P)
435 {
436   Standard_Integer i;
437   Standard_Real U;
438   myExtPC.Initialize(mysample, myintuinf, myintusup, mytolu, mytolf);
439   myExtPC.Perform(P);
440   mydone = myExtPC.IsDone();
441   if (mydone)
442   {
443     Standard_Integer NbExt = myExtPC.NbExt();
444     for (i = 1; i <= NbExt; i++)
445     {
446       // Verification de la validite des parametres pour le cas trimme:
447       ThePOnC PC = myExtPC.Point(i);
448       U = PC.Parameter();
449       if (TheCurveTool::IsPeriodic(*((TheCurve*)myC)))
450       {
451         U = ElCLib::InPeriod(U, myuinf, myuinf+TheCurveTool::Period(*((TheCurve*)myC)));
452       }
453       if ((U >= myuinf - mytolu) && (U <= myusup + mytolu))
454       {
455         PC.SetValues(U, PC.Value());
456         mySqDist.Append(myExtPC.SquareDistance(i));
457         myismin.Append(myExtPC.IsMin(i));
458         mypoint.Append(PC);
459       }
460     }
461   }
462 }
463
464
465
466
467 //=======================================================================
468 //function : Extrema_GExtPC
469 //purpose  : 
470 //=======================================================================
471
472 Extrema_GExtPC::Extrema_GExtPC()
473 {
474   myC = 0;
475   mydone = Standard_False;
476   mydist1 = RealLast();
477   mydist2 = RealLast();
478   mytolu = 0.0;
479   mytolf = 0.0;
480   mysample = 17;
481   myintuinf = myintusup = myuinf = myusup = Precision::Infinite();
482   type = GeomAbs_OtherCurve;
483 }
484
485 //=======================================================================
486 //function : Extrema_GExtPC
487 //purpose  : 
488 //=======================================================================
489
490 Extrema_GExtPC::Extrema_GExtPC(const ThePoint&           P, 
491                                const TheCurve&           C,
492                                const Standard_Real       Uinf,
493                                const Standard_Real       Usup,
494                                const Standard_Real       TolF) 
495 {
496   Initialize(C, Uinf, Usup, TolF);
497   Perform(P);
498 }
499
500 //=======================================================================
501 //function : Extrema_GExtPC
502 //purpose  : 
503 //=======================================================================
504
505 Extrema_GExtPC::Extrema_GExtPC(const ThePoint&     P, 
506                                const TheCurve&     C,
507                                const Standard_Real TolF) 
508 {
509   Initialize(C, TheCurveTool::FirstParameter(C), 
510              TheCurveTool::LastParameter(C), TolF);
511   Perform(P);
512 }
513
514
515 //=======================================================================
516 //function : IsDone
517 //purpose  : 
518 //=======================================================================
519
520 Standard_Boolean Extrema_GExtPC::IsDone() const
521 {
522   return mydone;
523 }
524
525
526 //=======================================================================
527 //function : Value
528 //purpose  : 
529 //=======================================================================
530
531 Standard_Real Extrema_GExtPC::SquareDistance(const Standard_Integer N) const 
532 {
533   if(!mydone) throw StdFail_NotDone();
534   if ((N < 1) || (N > mySqDist.Length())) throw Standard_OutOfRange();
535   return mySqDist.Value(N);
536 }
537
538
539 //=======================================================================
540 //function : NbExt
541 //purpose  : 
542 //=======================================================================
543
544 Standard_Integer Extrema_GExtPC::NbExt() const
545 {
546   if(!mydone) throw StdFail_NotDone();
547   return mySqDist.Length();
548 }
549
550
551 //=======================================================================
552 //function : IsMin
553 //purpose  : 
554 //=======================================================================
555
556 Standard_Boolean Extrema_GExtPC::IsMin(const Standard_Integer N) const
557 {
558   if(!mydone) throw StdFail_NotDone();
559   if ((N < 1) || (N > mySqDist.Length())) throw Standard_OutOfRange();
560   return myismin.Value(N);
561 }
562
563
564
565 //=======================================================================
566 //function : Point
567 //purpose  : 
568 //=======================================================================
569
570 const ThePOnC & Extrema_GExtPC::Point(const Standard_Integer N) const
571 {
572   if(!mydone) throw StdFail_NotDone();
573   if ((N < 1) || (N > mySqDist.Length())) throw Standard_OutOfRange();
574   return mypoint.Value(N);
575 }
576
577
578 //=======================================================================
579 //function : TrimmedDistances
580 //purpose  : 
581 //=======================================================================
582
583 void Extrema_GExtPC::TrimmedSquareDistances(Standard_Real& dist1, 
584                                       Standard_Real& dist2,
585                                       ThePoint& P1,
586                                       ThePoint& P2) const
587 {
588   dist1 = mydist1;
589   dist2 = mydist2;
590   P1 = Pf;
591   P2 = Pl;
592 }