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