0026252: GeomAdaptor_Surface should use inner adaptor to calculate values of complex...
[occt.git] / src / Extrema / Extrema_ExtPRevS.cxx
1 // Created on: 1999-09-21
2 // Created by: Edward AGAPOV
3 // Copyright (c) 1999-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
18 #include <Adaptor3d_HCurve.hxx>
19 #include <GeomAdaptor_HSurfaceOfRevolution.hxx>
20 #include <ElCLib.hxx>
21 #include <Extrema_ExtPElC.hxx>
22 #include <Extrema_ExtPRevS.hxx>
23 #include <Extrema_GenExtPS.hxx>
24 #include <Extrema_POnCurv.hxx>
25 #include <Extrema_POnSurf.hxx>
26 #include <gp_Ax1.hxx>
27 #include <gp_Ax2.hxx>
28 #include <gp_Lin.hxx>
29 #include <gp_Pln.hxx>
30 #include <gp_Pnt.hxx>
31 #include <gp_Trsf.hxx>
32 #include <gp_Vec.hxx>
33 #include <Precision.hxx>
34 #include <Standard_OutOfRange.hxx>
35 #include <Standard_Type.hxx>
36 #include <StdFail_NotDone.hxx>
37
38 static gp_Ax2 GetPosition (const GeomAdaptor_SurfaceOfRevolution& S)//const Handle(Adaptor_HCurve)& C)
39 {
40   Handle(Adaptor3d_HCurve) C = S.BasisCurve();
41   
42   switch (C->GetType()) {
43     
44   case GeomAbs_Line: {
45     gp_Lin L = C->Line();
46     gp_Dir N = S.AxeOfRevolution().Direction();
47     if (N.IsParallel(L.Direction(), Precision::Angular())) {
48       gp_Vec OO (L.Location(), S.AxeOfRevolution().Location());
49       if (OO.Magnitude() <= gp::Resolution()) {
50         OO = gp_Vec(L.Location(), ElCLib::Value(100,L));
51         if (N.IsParallel(OO, Precision::Angular()))
52           return gp_Ax2(); // Line and axe of revolution coinside
53       }
54       N ^= OO;
55     }
56     else {
57       N ^= L.Direction();
58     }
59     return gp_Ax2 (L.Location(), N, L.Direction());
60   }
61   case GeomAbs_Circle:
62     return C->Circle().Position();
63   case GeomAbs_Ellipse:
64     return C->Ellipse().Position();
65   case GeomAbs_Hyperbola:
66     return C->Hyperbola().Position();
67   case GeomAbs_Parabola:
68     return C->Parabola().Position();
69   default:
70     return gp_Ax2();
71   }
72 }
73
74 //=======================================================================
75 //function : HasSingularity
76 //purpose  : 
77 //=======================================================================
78
79 static Standard_Boolean HasSingularity(const GeomAdaptor_SurfaceOfRevolution& S) 
80 {
81
82   const Handle(Adaptor3d_HCurve) C = S.BasisCurve();
83   gp_Dir N = S.AxeOfRevolution().Direction();
84   gp_Pnt P = S.AxeOfRevolution().Location();
85
86   gp_Lin L(P, N);
87
88   P = C->Value(C->FirstParameter());
89
90   if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;
91
92   P = C->Value(C->LastParameter());
93
94   if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;
95   
96   return Standard_False;
97 }
98
99 //=============================================================================
100
101 static void PerformExtPElC (Extrema_ExtPElC& E,
102                             const gp_Pnt& P,
103                             const Handle(Adaptor3d_HCurve)& C,
104                             const Standard_Real Tol)
105 {
106   switch (C->GetType()) {
107   case GeomAbs_Hyperbola:
108     E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
109     return;
110   case GeomAbs_Line:
111     E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
112     return;
113   case GeomAbs_Circle:
114     E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
115     return;
116   case GeomAbs_Ellipse:
117     E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
118     return;
119   case GeomAbs_Parabola:
120     E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
121     return;
122   default:
123     return ;
124   }
125 }
126
127 //=======================================================================
128 //function : IsCaseAnalyticallyComputable
129 //purpose  : 
130 //=======================================================================
131
132 static Standard_Boolean IsCaseAnalyticallyComputable
133   (const GeomAbs_CurveType& theType,
134    const gp_Ax2& theCurvePos,
135    const gp_Ax1& AxeOfRevolution) 
136 {
137   // check type
138   switch (theType) {
139   case GeomAbs_Line:
140   case GeomAbs_Circle:
141   case GeomAbs_Ellipse:
142   case GeomAbs_Hyperbola:
143   case GeomAbs_Parabola:
144     break;
145   default:
146     return  Standard_False;
147   }
148 //  the axe of revolution must be in the plane of the curve.
149   gp_Pln pl(theCurvePos.Location(), theCurvePos.Direction());
150   gp_Pnt p1 = AxeOfRevolution.Location();
151   Standard_Real dist = 100., dist2 = dist * dist;
152   Standard_Real aThreshold = Precision::Angular() * Precision::Angular() * dist2;
153   gp_Pnt p2 = AxeOfRevolution.Location().XYZ() + dist * AxeOfRevolution.Direction().XYZ();
154
155   if((pl.SquareDistance(p1) < aThreshold) && 
156      (pl.SquareDistance(p2) < aThreshold))
157     return Standard_True;
158   return Standard_False;
159   //   gp_Vec V (AxeOfRevolution.Location(),theCurvePos.Location());
160   //   if (Abs( V * theCurvePos.Direction()) <= gp::Resolution())
161   //     return Standard_True;
162   //   else
163   //     return Standard_False;
164 }
165
166 //=======================================================================
167 //function : IsOriginalPnt
168 //purpose  : 
169 //=======================================================================
170
171 static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
172                                        const Extrema_POnSurf* Points,
173                                        const Standard_Integer NbPoints)
174 {
175   for (Standard_Integer i=1; i<=NbPoints; i++) {
176     if (Points[i-1].Value().IsEqual(P, Precision::Confusion())) {
177       return Standard_False;
178     }
179   }
180   return Standard_True;
181 }
182 //=======================================================================
183 //function : IsExtremum
184 //purpose  : 
185 //=======================================================================
186
187 static Standard_Boolean IsExtremum (const Standard_Real U, const Standard_Real V,
188                                     const gp_Pnt& P, const Adaptor3d_SurfacePtr& S,
189                                     gp_Pnt& E,        Standard_Real& Dist2,
190                                     const Standard_Boolean IsVSup,
191                                     const Standard_Boolean IsMin)
192 {
193   E = S->Value(U,V);
194   Dist2 = P.SquareDistance(E);
195   if (IsMin) 
196     return (Dist2 < P.SquareDistance(S->Value(U+1,V)) &&
197             Dist2 < P.SquareDistance(S->Value(U-1,V)) &&
198             Dist2 < P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
199   else
200     return (Dist2 > P.SquareDistance(S->Value(U+1,V)) &&
201             Dist2 > P.SquareDistance(S->Value(U-1,V)) &&
202             Dist2 > P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
203 }  
204 //=======================================================================
205 //function : Extrema_ExtPRevS
206 //purpose  : 
207 //=======================================================================
208
209 Extrema_ExtPRevS::Extrema_ExtPRevS() 
210 {
211   myDone = Standard_False;
212 }
213 //=======================================================================
214 //function : Extrema_ExtPRevS
215 //purpose  : 
216 //=======================================================================
217
218 Extrema_ExtPRevS::Extrema_ExtPRevS (const gp_Pnt&                                 theP,
219                                     const Handle(GeomAdaptor_HSurfaceOfRevolution)& theS,
220                                     const Standard_Real                           theUmin,
221                                     const Standard_Real                           theUsup,
222                                     const Standard_Real                           theVmin,
223                                     const Standard_Real                           theVsup,
224                                     const Standard_Real                           theTolU,
225                                     const Standard_Real                           theTolV)
226 {
227   Initialize (theS,
228               theUmin,
229               theUsup,
230               theVmin,
231               theVsup,
232               theTolU,
233               theTolV);
234
235   Perform (theP);
236 }
237 //=======================================================================
238 //function : Extrema_ExtPRevS
239 //purpose  : 
240 //=======================================================================
241
242 Extrema_ExtPRevS::Extrema_ExtPRevS (const gp_Pnt&                                 theP,
243                                     const Handle(GeomAdaptor_HSurfaceOfRevolution)& theS,
244                                     const Standard_Real                           theTolU,
245                                     const Standard_Real                           theTolV)
246 {
247   Initialize (theS,
248               theS->FirstUParameter(),
249               theS->LastUParameter(),
250               theS->FirstVParameter(),
251               theS->LastVParameter(),
252               theTolU,
253               theTolV);
254
255   Perform (theP);
256 }
257 //=======================================================================
258 //function : Initialize
259 //purpose  : 
260 //=======================================================================
261
262 void Extrema_ExtPRevS::Initialize (const Handle(GeomAdaptor_HSurfaceOfRevolution)& theS,
263                                    const Standard_Real                           theUmin,
264                                    const Standard_Real                           theUsup,
265                                    const Standard_Real                           theVmin,
266                                    const Standard_Real                           theVsup,
267                                    const Standard_Real                           theTolU,
268                                    const Standard_Real                           theTolV)
269 {
270   myvinf = theVmin;
271   myvsup = theVsup;
272   mytolv = theTolV;
273
274   Handle(Adaptor3d_HCurve) anACurve = theS->BasisCurve();
275   
276   if (myS != theS)
277   {
278     myS = theS;
279     myPosition = GetPosition (theS->ChangeSurface());
280     myIsAnalyticallyComputable =
281       IsCaseAnalyticallyComputable (anACurve->GetType(), myPosition, theS->AxeOfRevolution());
282   }
283
284   if (!myIsAnalyticallyComputable)
285   {
286     Standard_Integer aNbu = 32, aNbv = 32;
287
288     if (HasSingularity (theS->ChangeSurface()))
289     {
290       aNbv = 100;
291     }
292
293     myExtPS.Initialize (theS->ChangeSurface(),
294                         aNbu,
295                         aNbv,
296                         theUmin,
297                         theUsup,
298                         theVmin,
299                         theVsup,
300                         theTolU,
301                         theTolV);
302   }
303 }
304 //=======================================================================
305 //function : Perform
306 //purpose  : 
307 //=======================================================================
308
309 void Extrema_ExtPRevS::Perform(const gp_Pnt& P)
310 {
311   myDone = Standard_False;
312   myNbExt = 0;
313   
314   if (!myIsAnalyticallyComputable) {
315     
316     myExtPS.Perform(P);
317     myDone = myExtPS.IsDone();
318     myNbExt = myExtPS.NbExt();
319     return;
320   }
321   
322   Handle(Adaptor3d_HCurve) anACurve = myS->BasisCurve();
323
324   gp_Ax1 Ax = myS->AxeOfRevolution();
325   gp_Vec Dir = Ax.Direction(), Z = myPosition.Direction();
326   gp_Pnt O = Ax.Location();
327
328   Standard_Real OPdir = gp_Vec(O,P).Dot(Dir);
329   gp_Pnt Pp = P.Translated(Dir.Multiplied(-OPdir));
330   if (O.IsEqual(Pp,Precision::Confusion())) // P is on the AxeOfRevolution
331     return;
332   
333   Standard_Real U,V;
334   gp_Pnt P1, Ppp;
335   Standard_Real OPpz = gp_Vec(O,Pp).Dot(Z);
336   if (Abs(OPpz) <= gp::Resolution()) {
337     Ppp = Pp;
338     U = 0;
339   }
340   else {
341     Ppp = Pp.Translated(Z.Multiplied(-OPpz));
342     if (O.IsEqual(Ppp,Precision::Confusion())) 
343       U = M_PI/2;
344     else {
345       U = gp_Vec(O,Ppp).AngleWithRef(gp_Vec(O,Pp),Dir);
346     }
347   }
348
349   gp_Vec OPpp (O,Ppp), OPq (O, myS->Value(M_PI/2,0));
350   if (U != M_PI/2) {
351     if (Abs(OPq.Magnitude()) <= gp::Resolution()) 
352       OPq = gp_Vec(O, myS->Value(M_PI/2,anACurve->LastParameter()/10));
353     if (OPpp.AngleWithRef(OPq,Dir) < 0)
354       U += M_PI;
355   }
356   
357   gp_Trsf T;
358   T.SetRotation(Ax, -U);
359   P1 = P.Transformed(T);
360   
361   gp_Pnt E;
362   Standard_Real Dist2;
363   Standard_Integer i;
364   
365   Extrema_ExtPElC anExt;
366   PerformExtPElC(anExt, P1, anACurve, mytolv);
367   
368   if (anExt.IsDone()) {
369     myDone = Standard_True;
370     for (i=1; i<=anExt.NbExt(); i++) {
371       Extrema_POnCurv POC=anExt.Point(i);
372       V = POC.Parameter();
373       if (V > myvsup) {
374         //       if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
375         //                         Standard_True, anExt.IsMin(i))) continue;
376         Standard_Real newV = myvsup;
377
378         if((anACurve->GetType() == GeomAbs_Circle) || 
379            (anACurve->GetType() == GeomAbs_Ellipse)) {
380           newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
381
382           if (newV > myvsup) {
383             newV -= 2. * M_PI;
384
385             if (newV + mytolv < myvinf) {
386               newV = myvsup;
387             } else if (newV < myvinf) {
388               newV = myvinf;
389             }
390           }
391         }
392         V = newV;
393
394         if ( !IsExtremum (U, V, P, &(myS->ChangeSurface()), E, Dist2,
395                            Standard_True, anExt.IsMin(i))) {
396           continue;
397         }
398       } else if (V < myvinf) {
399         //      if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
400         //                        Standard_False, anExt.IsMin(i))) continue;
401
402         Standard_Real newV = myvinf;
403
404         if((anACurve->GetType() == GeomAbs_Circle) || 
405            (anACurve->GetType() == GeomAbs_Ellipse)) {
406           newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
407           
408           if(newV < myvinf) {
409             newV += 2. * M_PI;
410  
411             if (newV - mytolv > myvsup) {
412               newV = myvinf;
413             } else if (newV > myvsup) {
414               newV = myvsup;
415             }
416           }
417         }
418         V = newV;
419
420         if ( !IsExtremum (U, V, P, &(myS->ChangeSurface()), E, Dist2,
421                           Standard_False, anExt.IsMin(i))) continue;
422       } else {
423         E = myS->Value(U,V);
424         Dist2 = P.SquareDistance(E);
425       }
426       if (IsOriginalPnt(E, myPoint, myNbExt)) {
427         myPoint[myNbExt] = Extrema_POnSurf(U,V,E);
428         mySqDist[myNbExt] = Dist2;
429         myNbExt++;
430       }
431     }
432   }
433   T.SetRotation(Ax, M_PI);
434   P1.Transform(T);
435   
436   PerformExtPElC(anExt, P1, anACurve, mytolv);
437   if (anExt.IsDone()) {
438     myDone = Standard_True;
439
440     U += M_PI;
441     
442     for (i=1; i<=anExt.NbExt(); i++) {
443       Extrema_POnCurv POC=anExt.Point(i);
444       V = POC.Parameter();
445       if (V > myvsup) {
446         //      if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
447         //                         Standard_True, anExt.IsMin(i))) continue;
448
449         Standard_Real newV = myvsup;
450
451         if((anACurve->GetType() == GeomAbs_Circle) || 
452            (anACurve->GetType() == GeomAbs_Ellipse)) {
453           newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
454
455           if (newV > myvsup) {
456             newV -= 2. * M_PI;
457
458             if (newV + mytolv < myvinf) {
459               newV = myvsup;
460             } else if (newV < myvinf) {
461               newV = myvinf;
462             }
463           }
464         }
465         V = newV;
466         
467         if ( !IsExtremum (U, V, P, &(myS->ChangeSurface()), E, Dist2,
468                            Standard_True, anExt.IsMin(i))) continue;
469       } else if (V < myvinf) {
470         //      if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
471         //                        Standard_False, anExt.IsMin(i))) continue;
472         Standard_Real newV = myvinf;
473
474         if((anACurve->GetType() == GeomAbs_Circle) || 
475            (anACurve->GetType() == GeomAbs_Ellipse)) {
476           newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
477           
478           if(newV < myvinf) {
479             newV += 2. * M_PI;
480  
481             if (newV - mytolv > myvsup) {
482               newV = myvinf;
483             } else if (newV > myvsup) {
484               newV = myvsup;
485             }
486           }
487         }
488         V = newV;
489
490         if ( !IsExtremum (U, V, P, &(myS->ChangeSurface()), E, Dist2,
491                           Standard_False, anExt.IsMin(i))) continue;
492       } else {
493         E = myS->Value(U,V);
494         Dist2 = P.SquareDistance(E);
495       }
496       if (IsOriginalPnt(E, myPoint, myNbExt)) {
497         myPoint[myNbExt] = Extrema_POnSurf(U,V,E);
498         mySqDist[myNbExt] = Dist2;
499         myNbExt++;
500       }
501     }
502   }
503 }
504
505
506 //=======================================================================
507 //function : IsDone
508 //purpose  : 
509 //=======================================================================
510
511 Standard_Boolean Extrema_ExtPRevS::IsDone() const
512 {
513   return myDone; 
514 }
515
516
517 //=======================================================================
518 //function : NbExt
519 //purpose  : 
520 //=======================================================================
521
522 Standard_Integer Extrema_ExtPRevS::NbExt() const
523 {
524   if (!IsDone()) { StdFail_NotDone::Raise(); }
525   return myNbExt;
526 }
527
528
529 //=======================================================================
530 //function : Value
531 //purpose  : 
532 //=======================================================================
533
534 Standard_Real Extrema_ExtPRevS::SquareDistance(const Standard_Integer N) const
535 {
536   if (!IsDone()) { StdFail_NotDone::Raise(); }
537   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
538   if (myIsAnalyticallyComputable)
539     return mySqDist[N-1];
540   else
541     return myExtPS.SquareDistance(N);
542 }
543 //=======================================================================
544 //function : Point
545 //purpose  : 
546 //=======================================================================
547
548 const Extrema_POnSurf& Extrema_ExtPRevS::Point(const Standard_Integer N) const
549 {
550   if (!IsDone()) { StdFail_NotDone::Raise(); }
551   if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
552   if (myIsAnalyticallyComputable)
553     return myPoint[N-1];
554   else
555     return myExtPS.Point(N);
556 }
557
558