0025621: CAST analysis - Avoid constructors not supplying an initial value for all...
[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 : myuinf(0.0),
149   myusup(0.0),
150   mytolu(0.0),
151   myvinf(0.0),
152   myvsup(0.0),
153   mytolv(0.0),
154   myIsAnalyticallyComputable(Standard_False),
155   myDone(Standard_False),
156   myNbExt(Standard_False)
157 {
158 }
159
160 //=============================================================================
161
162 Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt&                                       theP,
163                                     const Handle(Adaptor3d_HSurfaceOfLinearExtrusion)&  theS,
164                                     const Standard_Real                                 theUmin,
165                                     const Standard_Real                                 theUsup,
166                                     const Standard_Real                                 theVmin,
167                                     const Standard_Real                                 theVsup,
168                                     const Standard_Real                                 theTolU,
169                                     const Standard_Real                                 theTolV)
170 : myuinf(theUmin),
171   myusup(theUsup),
172   mytolu(theTolU),
173   myvinf(theVmin),
174   myvsup(theVsup),
175   mytolv(theTolV),
176   myS   (theS),
177   myIsAnalyticallyComputable(Standard_False),
178   myDone(Standard_False),
179   myNbExt(Standard_False)
180 {
181   Initialize (theS,
182               theUmin,
183               theUsup,
184               theVmin,
185               theVsup,
186               theTolU,
187               theTolV);
188
189   Perform (theP);
190 }
191 //=============================================================================
192
193 Extrema_ExtPExtS::Extrema_ExtPExtS (const gp_Pnt&                                       theP,
194                                     const Handle(Adaptor3d_HSurfaceOfLinearExtrusion)&  theS,
195                                     const Standard_Real                                 theTolU, 
196                                     const Standard_Real                                 theTolV)
197 : myuinf(theS->FirstUParameter()),
198   myusup(theS->LastUParameter()),
199   mytolu(theTolU),
200   myvinf(theS->FirstVParameter()),
201   myvsup(theS->LastVParameter()),
202   mytolv(theTolV),
203   myS   (theS),
204   myIsAnalyticallyComputable(Standard_False),
205   myDone(Standard_False),
206   myNbExt(Standard_False)
207 {
208   Initialize (theS,
209               theS->FirstUParameter(),
210               theS->LastUParameter(),
211               theS->FirstVParameter(),
212               theS->LastVParameter(),
213               theTolU,
214               theTolV);
215
216   Perform (theP);
217 }
218
219 //=======================================================================
220 //function : Initialize
221 //purpose  : 
222 //=======================================================================
223
224 void Extrema_ExtPExtS::Initialize (const Handle(Adaptor3d_HSurfaceOfLinearExtrusion)&  theS,
225                                    const Standard_Real                                 theUinf,
226                                    const Standard_Real                                 theUsup,
227                                    const Standard_Real                                 theVinf,
228                                    const Standard_Real                                 theVsup,
229                                    const Standard_Real                                 theTolU,
230                                    const Standard_Real                                 theTolV)
231 {
232   myuinf = theUinf;
233   myusup = theUsup;
234   mytolu = theTolU;
235   
236   myvinf = theVinf;
237   myvsup = theVsup;
238   mytolv = theTolV;
239   
240   Handle(Adaptor3d_HCurve) anACurve = theS->BasisCurve();
241
242   myF.Initialize (theS->ChangeSurface());
243   myC = anACurve;
244   myS = theS;
245   myPosition = GetPosition(myC);
246   myDirection = theS->Direction();
247   myIsAnalyticallyComputable = //Standard_False;
248     IsCaseAnalyticallyComputable (myC->GetType(), myPosition, myDirection);
249   
250   if (!myIsAnalyticallyComputable)
251   {
252     myExtPS.Initialize (theS->ChangeSurface(),
253                         32,
254                         32,
255                         theUinf,
256                         theUsup,
257                         theVinf,
258                         theVsup,
259                         theTolU,
260                         theTolV);
261   }
262 }
263
264
265 //=======================================================================
266 //function : Perform
267 //purpose  : 
268 //=======================================================================
269
270 void Extrema_ExtPExtS::Perform (const gp_Pnt& P)
271 {
272   const Standard_Integer NbExtMax = 4; //dimension of arrays
273                                        //myPoint[] and mySqDist[]
274                                        //For "analytical" case
275   myDone = Standard_False;
276   myNbExt = 0;
277   
278   if (!myIsAnalyticallyComputable) {
279     myExtPS.Perform(P);
280     myDone = myExtPS.IsDone();
281 //  modified by NIZHNY-EAP Wed Nov 17 12:59:08 1999 ___BEGIN___
282     myNbExt = myExtPS.NbExt();
283 //  modified by NIZHNY-EAP Wed Nov 17 12:59:09 1999 ___END___
284     return;
285   }
286   
287   gp_Pnt Pe, Pp = ProjectPnt(myPosition,myDirection,P);
288   Extrema_ExtPElC anExt;
289   PerformExtPElC(anExt, Pp, myC, mytolu);
290   if (!anExt.IsDone()) return;
291   
292   gp_Ax2 anOrtogSection (P, myDirection);
293   Standard_Real U,V;
294   Standard_Boolean
295     isMin,
296     isSimpleCase =
297       myDirection.IsParallel(myPosition.Direction(),Precision::Angular());
298   Standard_Integer i, aNbExt = anExt.NbExt();
299   math_Vector UV(1,2), Tol(1,2), UVinf(1,2), UVsup(1,2);
300   Tol(1) = mytolu;   Tol(2) = mytolv;
301   UVinf(1) = myuinf; UVinf(2) = myvinf;
302   UVsup(1) = myusup; UVsup(2) = myvsup;
303   
304   
305   for (i=1; i<=aNbExt; i++) {
306     Extrema_POnCurv POC=anExt.Point(i);
307     U = POC.Parameter();
308     //// modified by jgv, 23.12.2008 for OCC17194 ////
309     if (myC->IsPeriodic())
310     {
311       Standard_Real U2 = U;
312       ElCLib::AdjustPeriodic(myuinf, myuinf + 2.*M_PI, Precision::PConfusion(), U, U2);
313     }
314     //////////////////////////////////////////////////
315     gp_Pnt E = POC.Value();
316     Pe = ProjectPnt(anOrtogSection, myDirection, E);
317     
318     if (isSimpleCase) {
319       V = gp_Vec(E,Pe) * gp_Vec(myDirection);
320       // modified by NIZHNY-MKK  Thu Sep 18 14:46:14 2003.BEGIN
321       //       myPoint[++myNbExt] = Extrema_POnSurf(U, V, Pe);
322       //       myValue[myNbExt] = anExt.Value(i);
323       myPoint[myNbExt] = Extrema_POnSurf(U, V, Pe);
324       mySqDist[myNbExt] = anExt.SquareDistance(i);
325       myNbExt++;
326       if(myNbExt == NbExtMax)
327       {
328         break;
329       }
330       // modified by NIZHNY-MKK  Thu Sep 18 14:46:18 2003.END
331     }
332     else {
333       myF.SetPoint(P);
334       isMin = anExt.IsMin(i);//( Pp.Distance(GetValue(U+10,myC)) > anExt.Value(i) );
335       
336       MakePreciser(U, P, isMin, anOrtogSection);
337       E = GetValue(U, myC);
338       Pe = ProjectPnt (anOrtogSection, myDirection, E),
339       V = gp_Vec(E,Pe) * gp_Vec(myDirection);
340       UV(1) = U;         UV(2) = V;
341       math_FunctionSetRoot aFSR (myF,UV,Tol,UVinf,UVsup);
342 //      for (Standard_Integer k=1 ; k <= myF.NbExt(); 
343       Standard_Integer k;
344       for ( k=1 ; k <= myF.NbExt(); k++) {
345               if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
346                 // modified by NIZHNY-MKK  Thu Sep 18 14:46:41 2003.BEGIN
347                 //        myPoint[++myNbExt] = myF.Point(k);
348                 //        myValue[myNbExt] = myF.Value(k);
349                 myPoint[myNbExt] = myF.Point(k);
350                 mySqDist[myNbExt] = myF.SquareDistance(k);
351                 myNbExt++;
352           if(myNbExt == NbExtMax)
353           {
354             break;
355           }
356                 // modified by NIZHNY-MKK  Thu Sep 18 14:46:43 2003.END
357               }
358       }
359       if(myNbExt == NbExtMax)
360       {
361         break;
362       }
363       // try symmetric point
364       myF.SetPoint(P); //To clear previous solutions
365       U *= -1;
366       MakePreciser(U, P, isMin, anOrtogSection);
367       E = GetValue(U, myC);
368       Pe = ProjectPnt (anOrtogSection, myDirection, E),
369       V = gp_Vec(E,Pe) * gp_Vec(myDirection); 
370       UV(1) = U; UV(2) = V;
371       
372       aFSR.Perform (myF,UV,UVinf,UVsup);
373       
374       for (k=1 ; k <= myF.NbExt(); k++) {
375         if(myF.SquareDistance(k) > Precision::Confusion()*Precision::Confusion())
376         {
377           //Additional checking solution: FSR sometimes is wrong
378           //when starting point is far from solution.
379           Standard_Real dist = Sqrt(myF.SquareDistance(k));
380           math_Vector Vals(1, 2);
381           const Extrema_POnSurf& PonS=myF.Point(k);
382           Standard_Real u, v;
383           PonS.Parameter(u, v);
384           UV(1) = u;
385           UV(2) = v;
386           myF.Value(UV, Vals);
387           gp_Vec du, dv;
388           myS->D1(u, v, Pe, du, dv);
389           Standard_Real mdu = du.Magnitude();
390           Standard_Real mdv = dv.Magnitude();
391           u = Abs(Vals(1));
392           v = Abs(Vals(2));
393           if(mdu > Precision::PConfusion())
394           {
395             if(u / dist / mdu > Precision::PConfusion())
396             {
397               continue;
398             }
399           }
400           if(mdv > Precision::PConfusion())
401           {
402             if(v / dist / mdv > Precision::PConfusion())
403             {
404               continue;
405             }
406           }
407
408         }
409               if (IsOriginalPnt(myF.Point(k).Value(), myPoint, myNbExt)) {
410                 // modified by NIZHNY-MKK  Thu Sep 18 14:46:59 2003.BEGIN
411                 //        myPoint[++myNbExt] = myF.Point(k);
412                 //        myValue[myNbExt] = myF.Value(k);
413                 myPoint[myNbExt] = myF.Point(k);
414                 mySqDist[myNbExt] = myF.SquareDistance(k);
415                 myNbExt++;
416           if(myNbExt == NbExtMax)
417           {
418             break;
419           }
420                 // modified by NIZHNY-MKK  Thu Sep 18 14:47:04 2003.END
421               }
422       }
423       if(myNbExt == NbExtMax)
424       {
425         break;
426       }
427     }
428   }
429   myDone = Standard_True;
430   return;
431 }
432
433 //=============================================================================
434
435 Standard_Boolean Extrema_ExtPExtS::IsDone () const { return myDone; }
436 //=============================================================================
437
438 Standard_Integer Extrema_ExtPExtS::NbExt () const
439 {
440   if (!IsDone()) { StdFail_NotDone::Raise(); }
441   if (myIsAnalyticallyComputable)
442     return myNbExt;
443   else
444     return myExtPS.NbExt();
445 }
446 //=============================================================================
447
448 Standard_Real Extrema_ExtPExtS::SquareDistance (const Standard_Integer N) const
449 {
450   if (!IsDone()) { StdFail_NotDone::Raise(); }
451   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
452   if (myIsAnalyticallyComputable)
453     // modified by NIZHNY-MKK  Thu Sep 18 14:48:39 2003.BEGIN
454     //     return myValue[N];
455     return mySqDist[N-1];
456   // modified by NIZHNY-MKK  Thu Sep 18 14:48:42 2003.END
457   else
458     return myExtPS.SquareDistance(N);
459 }
460 //=============================================================================
461
462 const Extrema_POnSurf& Extrema_ExtPExtS::Point (const Standard_Integer N) const
463 {
464   if (!IsDone()) { StdFail_NotDone::Raise(); }
465   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
466   if (myIsAnalyticallyComputable) {
467     // modified by NIZHNY-MKK  Thu Sep 18 14:47:40 2003.BEGIN
468     //     return myPoint[N];
469     return myPoint[N-1];
470   }
471   // modified by NIZHNY-MKK  Thu Sep 18 14:47:43 2003.END
472   else
473     return myExtPS.Point(N);
474 }
475 //=============================================================================
476
477
478 static gp_Ax2 GetPosition (const Handle(Adaptor3d_HCurve)& C)
479 {
480   switch (C->GetType()) {
481   case GeomAbs_Line: {
482     gp_Lin L = C->Line();
483     gp_Pln Pln = gp_Pln(L.Location(), L.Direction());
484     //:abv 30.05.02: OCC  - use constructor instead of Set...s() to avoid exception
485     gp_Ax2 Pos ( Pln.Location(), Pln.Position().Direction(), Pln.Position().XDirection() );
486 //     Pos.SetAxis(Pln.XAxis());
487 //     Pos.SetXDirection(Pln.YAxis().Direction());
488 //     Pos.SetYDirection(Pln.Position().Direction());
489     return Pos;
490   }
491   case GeomAbs_Circle:
492     return  C->Circle().Position();
493   case GeomAbs_Ellipse:
494     return  C->Ellipse().Position();
495   case GeomAbs_Hyperbola:
496     return  C->Hyperbola().Position();
497   case GeomAbs_Parabola:
498     return  C->Parabola().Position();
499   default: 
500     return gp_Ax2 ();
501   }
502 }
503 //=============================================================================
504
505 static void PerformExtPElC (Extrema_ExtPElC& E,
506                             const gp_Pnt& P,
507                             const Handle(Adaptor3d_HCurve)& C,
508                             const Standard_Real Tol)
509 {
510   switch (C->GetType()) {
511   case GeomAbs_Hyperbola:
512     E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
513     return;
514   case GeomAbs_Line:
515     E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
516     return;
517   case GeomAbs_Circle:
518     E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
519     return;
520   case GeomAbs_Ellipse:
521     E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
522     return;
523   case GeomAbs_Parabola:
524     E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
525     return;
526   default:
527     return;
528   }
529 }
530
531 //=======================================================================
532 //function : IsCaseAnalyticallyComputable
533 //purpose  : 
534 //=======================================================================
535
536 static Standard_Boolean IsCaseAnalyticallyComputable
537   (const GeomAbs_CurveType& theType,
538    const gp_Ax2& theCurvePos,
539    const gp_Dir& theSurfaceDirection) 
540 {
541   // check type
542   switch (theType) {
543   case GeomAbs_Line:
544   case GeomAbs_Circle:
545   case GeomAbs_Ellipse:
546   case GeomAbs_Hyperbola:
547   case GeomAbs_Parabola:
548     break;
549   default:
550     return  Standard_False;
551   }
552   // check if it is a plane
553   if (Abs(theCurvePos.Direction() * theSurfaceDirection) <= gp::Resolution())
554     return Standard_False;
555   else
556     return Standard_True;
557 }
558 //=======================================================================
559 //function : GetValue
560 //purpose  : 
561 //=======================================================================
562
563 static gp_Pnt GetValue(const Standard_Real U,
564                        const Handle(Adaptor3d_HCurve)& C)
565 {
566   switch (C->GetType()) {
567   case GeomAbs_Line:
568     return ElCLib::Value(U, C->Line());
569   case GeomAbs_Circle:
570     return ElCLib::Value(U, C->Circle());
571   case GeomAbs_Ellipse:
572     return ElCLib::Value(U, C->Ellipse());
573   case GeomAbs_Hyperbola:
574     return ElCLib::Value(U, C->Hyperbola());
575   case GeomAbs_Parabola:
576     return ElCLib::Value(U, C->Parabola());
577   default:
578     return gp_Pnt();
579   }
580 }
581 //=======================================================================
582 //function : GetU
583 //purpose  : 
584 //=======================================================================
585 //#ifdef OCCT_DEBUG
586 //static Standard_Real GetU(const gp_Vec& vec,
587 //                        const gp_Pnt& P,
588 //                        const Handle(Adaptor3d_HCurve)& C)
589 //{
590 //  switch (C->GetType()) {
591 //  case GeomAbs_Line:
592 //    return ElCLib::Parameter(C->Line().Translated(vec), P);
593 //  case GeomAbs_Circle:
594 //    return ElCLib::Parameter(C->Circle().Translated(vec), P);
595 //  case GeomAbs_Ellipse:
596 //    return ElCLib::Parameter(C->Ellipse().Translated(vec), P);
597 //  case GeomAbs_Hyperbola:
598 //    return ElCLib::Parameter(C->Hyperbola().Translated(vec), P);
599 //  case GeomAbs_Parabola:
600 //    return ElCLib::Parameter(C->Parabola().Translated(vec), P);
601 //  default:
602 //    return 0;
603 //  }
604 //}
605 //#endif