d30349d6164e0e81c61dac896a3a62163480cf7e
[occt.git] / src / Extrema / Extrema_ExtPExtS.cxx
1 // Created on: 1999-09-16
2 // Created by: Edward AGAPOV
3 // Copyright (c) 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 #include <Standard_NotImplemented.hxx>
18 #include <Standard_OutOfRange.hxx>
19 #include <StdFail_NotDone.hxx>
20 #include <Adaptor3d_HCurve.hxx>
21 #include <ElCLib.hxx>
22 #include <Extrema_ExtPElC.hxx>
23 #include <Extrema_ExtPExtS.hxx>
24 #include <Extrema_POnCurv.hxx>
25 #include <Extrema_POnSurf.hxx>
26 #include <Precision.hxx>
27 #include <gp.hxx>
28 #include <gp_Ax2.hxx>
29 #include <gp_Dir.hxx>
30 #include <gp_Lin.hxx>
31 #include <gp_Pln.hxx>
32 #include <gp_Pnt.hxx>
33 #include <gp_Vec.hxx>
34 #include <math_FunctionSetRoot.hxx>
35 #include <math_Vector.hxx>
36 #include <Adaptor3d_HSurfaceOfLinearExtrusion.hxx>
37
38 IMPLEMENT_STANDARD_HANDLE (Extrema_ExtPExtS, Standard_Transient)
39 IMPLEMENT_STANDARD_RTTIEXT(Extrema_ExtPExtS, Standard_Transient)
40
41 static gp_Ax2 GetPosition (const Handle(Adaptor3d_HCurve)& C);
42      
43 static void PerformExtPElC (Extrema_ExtPElC& E,
44                             const gp_Pnt& P,
45                             const Handle(Adaptor3d_HCurve)& C,
46                             const Standard_Real Tol);
47      
48 static Standard_Boolean
49   IsCaseAnalyticallyComputable (const GeomAbs_CurveType& theType,
50                                 const gp_Ax2& theCurvePos,
51                                 const gp_Dir& theSurfaceDirection) ;
52
53 static gp_Pnt GetValue(const Standard_Real U,
54                        const Handle(Adaptor3d_HCurve)& C);
55 //=======================================================================
56 //function : Project
57 //purpose  : Returns the projection of a point <Point> on a plane 
58 //           <ThePlane>  along a direction <TheDir>.
59 //=======================================================================
60 static gp_Pnt ProjectPnt(const gp_Ax2& ThePlane,
61                          const gp_Dir& TheDir,
62                          const gp_Pnt& Point) 
63 {
64   gp_Vec PO(Point,ThePlane.Location());
65   Standard_Real Alpha = PO * gp_Vec(ThePlane.Direction());
66   Alpha /= TheDir * ThePlane.Direction();
67   gp_Pnt P;
68   P.SetXYZ(Point.XYZ() + Alpha * TheDir.XYZ());
69   return P;
70 }
71
72 //=======================================================================
73 //function : IsOriginalPnt
74 //purpose  : 
75 //=======================================================================
76
77 static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
78                                        const Extrema_POnSurf* Points,
79                                        const Standard_Integer NbPoints)
80 {
81   for (Standard_Integer i=1; i<=NbPoints; i++) {
82     if (Points[i-1].Value().IsEqual(P, Precision::Confusion())) {
83       return Standard_False;
84     }
85   }
86   return Standard_True;
87 }
88
89 //=======================================================================
90 //function : MakePreciser
91 //purpose  : 
92 //=======================================================================
93
94 void Extrema_ExtPExtS::MakePreciser (Standard_Real& U,
95                                      const gp_Pnt& P,
96                                      const Standard_Boolean isMin,
97                                      const gp_Ax2& OrtogSection) const
98 {
99   if (U>myusup) {
100     U = myusup;
101   } else if (U<myuinf) {
102     U = myuinf;
103   } else {
104     
105     Standard_Real step = (myusup - myuinf) / 30, D2e, D2next ,D2prev;
106     gp_Pnt
107       Pe = ProjectPnt (OrtogSection, myDirection, GetValue(U,myC)),
108       Pprev = ProjectPnt (OrtogSection, myDirection, GetValue(U-step, myC)),
109       Pnext = ProjectPnt (OrtogSection, myDirection, GetValue(U+step, myC));
110     D2e = P.SquareDistance(Pe),
111     D2next = P.SquareDistance(Pnext),
112     D2prev = P.SquareDistance(Pprev);
113     Standard_Boolean notFound;
114     if (isMin) 
115       notFound =  (D2e > D2prev || D2e > D2next);
116     else 
117       notFound = (D2e < D2prev || D2e < D2next);
118     
119     if (notFound && (D2e < D2next && isMin)) {
120       step = -step;
121       D2next = D2prev;
122       Pnext = Pprev;
123     }
124     while (notFound) {
125       U = U + step;
126       if (U > myusup) {
127               U = myusup;
128               break;
129       }
130       if (U < myuinf) {
131               U = myuinf;
132               break;
133       }
134       D2e = D2next;
135       Pe = Pnext;
136       Pnext = ProjectPnt (OrtogSection, myDirection, GetValue(U+step, myC));
137       D2next = P.SquareDistance(Pnext);
138       if (isMin) 
139               notFound = D2e > D2next;
140       else 
141               notFound = D2e < D2next;
142     }
143   }
144 }
145 //=============================================================================
146
147 Extrema_ExtPExtS::Extrema_ExtPExtS ()
148 {
149   myDone = Standard_False;
150 }
151
152 //=============================================================================
153
154 Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt&                                       theP,
155                                     const Handle(Adaptor3d_HSurfaceOfLinearExtrusion)&  theS,
156                                     const Standard_Real                                 theUmin,
157                                     const Standard_Real                                 theUsup,
158                                     const Standard_Real                                 theVmin,
159                                     const Standard_Real                                 theVsup,
160                                     const Standard_Real                                 theTolU,
161                                     const Standard_Real                                 theTolV)
162 {
163   Initialize (theS,
164               theUmin,
165               theUsup,
166               theVmin,
167               theVsup,
168               theTolU,
169               theTolV);
170
171   Perform (theP);
172 }
173 //=============================================================================
174
175 Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt&                                       theP,
176                                     const Handle(Adaptor3d_HSurfaceOfLinearExtrusion)&  theS,
177                                     const Standard_Real                                 theTolU, 
178                                     const Standard_Real                                 theTolV)
179 {
180   Initialize (theS,
181               theS->FirstUParameter(),
182               theS->LastUParameter(),
183               theS->FirstVParameter(),
184               theS->LastVParameter(),
185               theTolU,
186               theTolV);
187
188   Perform (theP);
189 }
190 //=======================================================================
191 //function : Initialize
192 //purpose  : 
193 //=======================================================================
194
195 void Extrema_ExtPExtS::Initialize (const Handle(Adaptor3d_HSurfaceOfLinearExtrusion)&  theS,
196                                    const Standard_Real                                 theUinf,
197                                    const Standard_Real                                 theUsup,
198                                    const Standard_Real                                 theVinf,
199                                    const Standard_Real                                 theVsup,
200                                    const Standard_Real                                 theTolU,
201                                    const Standard_Real                                 theTolV)
202 {
203   myuinf = theUinf;
204   myusup = theUsup;
205   mytolu = theTolU;
206   
207   myvinf = theVinf;
208   myvsup = theVsup;
209   mytolv = theTolV;
210   
211   Handle(Adaptor3d_HCurve) anACurve = theS->BasisCurve();
212
213   myF.Initialize (theS->ChangeSurface());
214   myC = anACurve;
215   myS = theS;
216   myPosition = GetPosition(myC);
217   myDirection = theS->Direction();
218   myIsAnalyticallyComputable = //Standard_False;
219     IsCaseAnalyticallyComputable (myC->GetType(), myPosition, myDirection);
220   
221   if (!myIsAnalyticallyComputable)
222   {
223     myExtPS.Initialize (theS->ChangeSurface(),
224                         32,
225                         32,
226                         theUinf,
227                         theUsup,
228                         theVinf,
229                         theVsup,
230                         theTolU,
231                         theTolV);
232   }
233 }
234
235
236 //=======================================================================
237 //function : Perform
238 //purpose  : 
239 //=======================================================================
240
241 void Extrema_ExtPExtS::Perform (const gp_Pnt& P)
242 {
243   const Standard_Integer NbExtMax = 4; //dimension of arrays
244                                        //myPoint[] and mySqDist[]
245                                        //For "analytical" case
246   myDone = Standard_False;
247   myNbExt = 0;
248   
249   if (!myIsAnalyticallyComputable) {
250     myExtPS.Perform(P);
251     myDone = myExtPS.IsDone();
252 //  modified by NIZHNY-EAP Wed Nov 17 12:59:08 1999 ___BEGIN___
253     myNbExt = myExtPS.NbExt();
254 //  modified by NIZHNY-EAP Wed Nov 17 12:59:09 1999 ___END___
255     return;
256   }
257   
258   gp_Pnt Pe, Pp = ProjectPnt(myPosition,myDirection,P);
259   Extrema_ExtPElC anExt;
260   PerformExtPElC(anExt, Pp, myC, mytolu);
261   if (!anExt.IsDone()) return;
262   
263   gp_Ax2 anOrtogSection (P, myDirection);
264   Standard_Real U,V;
265   Standard_Boolean
266     isMin,
267     isSimpleCase =
268       myDirection.IsParallel(myPosition.Direction(),Precision::Angular());
269   Standard_Integer i, aNbExt = anExt.NbExt();
270   math_Vector UV(1,2), Tol(1,2), UVinf(1,2), UVsup(1,2);
271   Tol(1) = mytolu;   Tol(2) = mytolv;
272   UVinf(1) = myuinf; UVinf(2) = myvinf;
273   UVsup(1) = myusup; UVsup(2) = myvsup;
274   
275   
276   for (i=1; i<=aNbExt; i++) {
277     Extrema_POnCurv POC=anExt.Point(i);
278     U = POC.Parameter();
279     //// modified by jgv, 23.12.2008 for OCC17194 ////
280     if (myC->IsPeriodic())
281     {
282       Standard_Real U2 = U;
283       ElCLib::AdjustPeriodic(myuinf, myuinf + 2.*M_PI, Precision::PConfusion(), U, U2);
284     }
285     //////////////////////////////////////////////////
286     gp_Pnt E = POC.Value();
287     Pe = ProjectPnt(anOrtogSection, myDirection, E);
288     
289     if (isSimpleCase) {
290       V = gp_Vec(E,Pe) * gp_Vec(myDirection);
291       // modified by NIZHNY-MKK  Thu Sep 18 14:46:14 2003.BEGIN
292       //       myPoint[++myNbExt] = Extrema_POnSurf(U, V, Pe);
293       //       myValue[myNbExt] = anExt.Value(i);
294       myPoint[myNbExt] = Extrema_POnSurf(U, V, Pe);
295       mySqDist[myNbExt] = anExt.SquareDistance(i);
296       myNbExt++;
297       if(myNbExt == NbExtMax)
298       {
299         break;
300       }
301       // modified by NIZHNY-MKK  Thu Sep 18 14:46:18 2003.END
302     }
303     else {
304       myF.SetPoint(P);
305       isMin = anExt.IsMin(i);//( Pp.Distance(GetValue(U+10,myC)) > anExt.Value(i) );
306       
307       MakePreciser(U, P, isMin, anOrtogSection);
308       E = GetValue(U, myC);
309       Pe = ProjectPnt (anOrtogSection, myDirection, E),
310       V = gp_Vec(E,Pe) * gp_Vec(myDirection);
311       UV(1) = U;         UV(2) = V;
312       math_FunctionSetRoot aFSR (myF,UV,Tol,UVinf,UVsup);
313 //      for (Standard_Integer k=1 ; k <= myF.NbExt(); 
314       Standard_Integer k;
315       for ( k=1 ; k <= myF.NbExt(); k++) {
316               if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
317                 // modified by NIZHNY-MKK  Thu Sep 18 14:46:41 2003.BEGIN
318                 //        myPoint[++myNbExt] = myF.Point(k);
319                 //        myValue[myNbExt] = myF.Value(k);
320                 myPoint[myNbExt] = myF.Point(k);
321                 mySqDist[myNbExt] = myF.SquareDistance(k);
322                 myNbExt++;
323           if(myNbExt == NbExtMax)
324           {
325             break;
326           }
327                 // modified by NIZHNY-MKK  Thu Sep 18 14:46:43 2003.END
328               }
329       }
330       if(myNbExt == NbExtMax)
331       {
332         break;
333       }
334       // try symmetric point
335       myF.SetPoint(P); //To clear previous solutions
336       U *= -1;
337       MakePreciser(U, P, isMin, anOrtogSection);
338       E = GetValue(U, myC);
339       Pe = ProjectPnt (anOrtogSection, myDirection, E),
340       V = gp_Vec(E,Pe) * gp_Vec(myDirection); 
341       UV(1) = U; UV(2) = V;
342       
343       aFSR.Perform (myF,UV,UVinf,UVsup);
344       
345       for (k=1 ; k <= myF.NbExt(); k++) {
346         if(myF.SquareDistance(k) > Precision::Confusion()*Precision::Confusion())
347         {
348           //Additional checking solution: FSR sometimes is wrong
349           //when starting point is far from solution.
350           Standard_Real dist = Sqrt(myF.SquareDistance(k));
351           math_Vector Vals(1, 2);
352           const Extrema_POnSurf& PonS=myF.Point(k);
353           Standard_Real u, v;
354           PonS.Parameter(u, v);
355           UV(1) = u;
356           UV(2) = v;
357           myF.Value(UV, Vals);
358           gp_Vec du, dv;
359           myS->D1(u, v, Pe, du, dv);
360           Standard_Real mdu = du.Magnitude();
361           Standard_Real mdv = dv.Magnitude();
362           u = Abs(Vals(1));
363           v = Abs(Vals(2));
364           if(mdu > Precision::PConfusion())
365           {
366             if(u / dist / mdu > Precision::PConfusion())
367             {
368               continue;
369             }
370           }
371           if(mdv > Precision::PConfusion())
372           {
373             if(v / dist / mdv > Precision::PConfusion())
374             {
375               continue;
376             }
377           }
378
379         }
380               if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
381                 // modified by NIZHNY-MKK  Thu Sep 18 14:46:59 2003.BEGIN
382                 //        myPoint[++myNbExt] = myF.Point(k);
383                 //        myValue[myNbExt] = myF.Value(k);
384                 myPoint[myNbExt] = myF.Point(k);
385                 mySqDist[myNbExt] = myF.SquareDistance(k);
386                 myNbExt++;
387           if(myNbExt == NbExtMax)
388           {
389             break;
390           }
391                 // modified by NIZHNY-MKK  Thu Sep 18 14:47:04 2003.END
392               }
393       }
394       if(myNbExt == NbExtMax)
395       {
396         break;
397       }
398     }
399   }
400   myDone = Standard_True;
401   return;
402 }
403
404 //=============================================================================
405
406 Standard_Boolean Extrema_ExtPExtS::IsDone () const { return myDone; }
407 //=============================================================================
408
409 Standard_Integer Extrema_ExtPExtS::NbExt () const
410 {
411   if (!IsDone()) { StdFail_NotDone::Raise(); }
412   if (myIsAnalyticallyComputable)
413     return myNbExt;
414   else
415     return myExtPS.NbExt();
416 }
417 //=============================================================================
418
419 Standard_Real Extrema_ExtPExtS::SquareDistance (const Standard_Integer N) const
420 {
421   if (!IsDone()) { StdFail_NotDone::Raise(); }
422   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
423   if (myIsAnalyticallyComputable)
424     // modified by NIZHNY-MKK  Thu Sep 18 14:48:39 2003.BEGIN
425     //     return myValue[N];
426     return mySqDist[N-1];
427   // modified by NIZHNY-MKK  Thu Sep 18 14:48:42 2003.END
428   else
429     return myExtPS.SquareDistance(N);
430 }
431 //=============================================================================
432
433 const Extrema_POnSurf& Extrema_ExtPExtS::Point (const Standard_Integer N) const
434 {
435   if (!IsDone()) { StdFail_NotDone::Raise(); }
436   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
437   if (myIsAnalyticallyComputable) {
438     // modified by NIZHNY-MKK  Thu Sep 18 14:47:40 2003.BEGIN
439     //     return myPoint[N];
440     return myPoint[N-1];
441   }
442   // modified by NIZHNY-MKK  Thu Sep 18 14:47:43 2003.END
443   else
444     return myExtPS.Point(N);
445 }
446 //=============================================================================
447
448
449 static gp_Ax2 GetPosition (const Handle(Adaptor3d_HCurve)& C)
450 {
451   switch (C->GetType()) {
452   case GeomAbs_Line: {
453     gp_Lin L = C->Line();
454     gp_Pln Pln = gp_Pln(L.Location(), L.Direction());
455     //:abv 30.05.02: OCC  - use constructor instead of Set...s() to avoid exception
456     gp_Ax2 Pos ( Pln.Location(), Pln.Position().Direction(), Pln.Position().XDirection() );
457 //     Pos.SetAxis(Pln.XAxis());
458 //     Pos.SetXDirection(Pln.YAxis().Direction());
459 //     Pos.SetYDirection(Pln.Position().Direction());
460     return Pos;
461   }
462   case GeomAbs_Circle:
463     return  C->Circle().Position();
464   case GeomAbs_Ellipse:
465     return  C->Ellipse().Position();
466   case GeomAbs_Hyperbola:
467     return  C->Hyperbola().Position();
468   case GeomAbs_Parabola:
469     return  C->Parabola().Position();
470   default: 
471     return gp_Ax2 ();
472   }
473 }
474 //=============================================================================
475
476 static void PerformExtPElC (Extrema_ExtPElC& E,
477                             const gp_Pnt& P,
478                             const Handle(Adaptor3d_HCurve)& C,
479                             const Standard_Real Tol)
480 {
481   switch (C->GetType()) {
482   case GeomAbs_Hyperbola:
483     E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
484     return;
485   case GeomAbs_Line:
486     E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
487     return;
488   case GeomAbs_Circle:
489     E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
490     return;
491   case GeomAbs_Ellipse:
492     E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
493     return;
494   case GeomAbs_Parabola:
495     E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
496     return;
497   default:
498     return;
499   }
500 }
501
502 //=======================================================================
503 //function : IsCaseAnalyticallyComputable
504 //purpose  : 
505 //=======================================================================
506
507 static Standard_Boolean IsCaseAnalyticallyComputable
508   (const GeomAbs_CurveType& theType,
509    const gp_Ax2& theCurvePos,
510    const gp_Dir& theSurfaceDirection) 
511 {
512   // check type
513   switch (theType) {
514   case GeomAbs_Line:
515   case GeomAbs_Circle:
516   case GeomAbs_Ellipse:
517   case GeomAbs_Hyperbola:
518   case GeomAbs_Parabola:
519     break;
520   default:
521     return  Standard_False;
522   }
523   // check if it is a plane
524   if (Abs(theCurvePos.Direction() * theSurfaceDirection) <= gp::Resolution())
525     return Standard_False;
526   else
527     return Standard_True;
528 }
529 //=======================================================================
530 //function : GetValue
531 //purpose  : 
532 //=======================================================================
533
534 static gp_Pnt GetValue(const Standard_Real U,
535                        const Handle(Adaptor3d_HCurve)& C)
536 {
537   switch (C->GetType()) {
538   case GeomAbs_Line:
539     return ElCLib::Value(U, C->Line());
540   case GeomAbs_Circle:
541     return ElCLib::Value(U, C->Circle());
542   case GeomAbs_Ellipse:
543     return ElCLib::Value(U, C->Ellipse());
544   case GeomAbs_Hyperbola:
545     return ElCLib::Value(U, C->Hyperbola());
546   case GeomAbs_Parabola:
547     return ElCLib::Value(U, C->Parabola());
548   default:
549     return gp_Pnt();
550   }
551 }
552 //=======================================================================
553 //function : GetU
554 //purpose  : 
555 //=======================================================================
556 //#ifdef OCCT_DEBUG
557 //static Standard_Real GetU(const gp_Vec& vec,
558 //                        const gp_Pnt& P,
559 //                        const Handle(Adaptor3d_HCurve)& C)
560 //{
561 //  switch (C->GetType()) {
562 //  case GeomAbs_Line:
563 //    return ElCLib::Parameter(C->Line().Translated(vec), P);
564 //  case GeomAbs_Circle:
565 //    return ElCLib::Parameter(C->Circle().Translated(vec), P);
566 //  case GeomAbs_Ellipse:
567 //    return ElCLib::Parameter(C->Ellipse().Translated(vec), P);
568 //  case GeomAbs_Hyperbola:
569 //    return ElCLib::Parameter(C->Hyperbola().Translated(vec), P);
570 //  case GeomAbs_Parabola:
571 //    return ElCLib::Parameter(C->Parabola().Translated(vec), P);
572 //  default:
573 //    return 0;
574 //  }
575 //}
576 //#endif