1b8ae22f282830517e0f32a9d1b71c9e906a14bb
[occt.git] / src / Extrema / Extrema_ExtCS.cxx
1 // Created on: 1995-07-18
2 // Created by: Modelistation
3 // Copyright (c) 1995-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 skv - Thu Jul  7 12:29:34 2005 OCC9134
18
19 #include <Adaptor3d_Curve.hxx>
20 #include <Adaptor3d_Surface.hxx>
21 #include <Bnd_Box.hxx>
22 #include <BndLib_AddSurface.hxx>
23 #include <ElCLib.hxx>
24 #include <ElSLib.hxx>
25 #include <Extrema_ExtCS.hxx>
26 #include <Extrema_ExtPElC.hxx>
27 #include <Extrema_ExtPElS.hxx>
28 #include <Extrema_ExtPS.hxx>
29 #include <Extrema_GenExtCS.hxx>
30 #include <Extrema_POnCurv.hxx>
31 #include <Extrema_POnSurf.hxx>
32 #include <GeomAbs_CurveType.hxx>
33 #include <gp_Cone.hxx>
34 #include <gp_Cylinder.hxx>
35 #include <gp_Lin.hxx>
36 #include <gp_Pln.hxx>
37 #include <gp_Pnt.hxx>
38 #include <gp_Sphere.hxx>
39 #include <gp_Torus.hxx>
40 #include <Precision.hxx>
41 #include <Standard_NotImplemented.hxx>
42 #include <Standard_OutOfRange.hxx>
43 #include <Standard_TypeMismatch.hxx>
44 #include <StdFail_InfiniteSolutions.hxx>
45 #include <StdFail_NotDone.hxx>
46 #include <TColStd_Array1OfReal.hxx>
47
48 Extrema_ExtCS::Extrema_ExtCS() 
49 {
50   myDone = Standard_False;
51 }
52
53 Extrema_ExtCS::Extrema_ExtCS(const Adaptor3d_Curve&   C,
54   const Adaptor3d_Surface& S,
55   const Standard_Real    TolC,
56   const Standard_Real    TolS)
57
58 {
59   Initialize(S, S.FirstUParameter(), S.LastUParameter(), 
60     S.FirstVParameter(), S.LastVParameter(), 
61     TolC, TolS);
62   Perform(C, C.FirstParameter(), C.LastParameter());
63 }
64
65 Extrema_ExtCS::Extrema_ExtCS(const Adaptor3d_Curve&   C,
66   const Adaptor3d_Surface& S,
67   const Standard_Real    UCinf,
68   const Standard_Real    UCsup,
69   const Standard_Real    Uinf,  
70   const Standard_Real    Usup,
71   const Standard_Real    Vinf,  
72   const Standard_Real    Vsup,
73   const Standard_Real    TolC,
74   const Standard_Real    TolS)
75
76 {
77   Initialize(S, Uinf, Usup, Vinf, Vsup, TolC, TolS);
78   Perform(C, UCinf, UCsup);
79 }
80
81
82 void Extrema_ExtCS::Initialize(const Adaptor3d_Surface& S,
83   const Standard_Real    Uinf,  
84   const Standard_Real    Usup,
85   const Standard_Real    Vinf,  
86   const Standard_Real    Vsup,
87   const Standard_Real    TolC,
88   const Standard_Real    TolS)
89 {
90   myS = (Adaptor3d_SurfacePtr)&S;
91   myIsPar = Standard_False;
92   myuinf  = Uinf;
93   myusup  = Usup;
94   myvinf  = Vinf;
95   myvsup  = Vsup;
96   mytolC  = TolC;
97   mytolS  = TolS;
98   myStype  = myS->GetType();
99 }
100
101
102 void Extrema_ExtCS::Perform(const Adaptor3d_Curve& C,
103   const Standard_Real  Uinf,
104   const Standard_Real  Usup)
105 {
106   myucinf = Uinf;
107   myucsup = Usup;
108   myPOnS.Clear();
109   myPOnC.Clear();
110   mySqDist.Clear();
111   Standard_Integer i, j;
112   Standard_Integer NbT, NbU, NbV;
113   NbT = NbU = NbV = 10;
114   GeomAbs_CurveType myCtype  = C.GetType();
115
116
117   switch(myCtype) {
118
119   case GeomAbs_Line: 
120     {
121
122       switch(myStype) {
123       case GeomAbs_Sphere:
124         myExtElCS.Perform(C.Line(), myS->Sphere());
125         break;   
126       case GeomAbs_Cylinder:
127         myExtElCS.Perform(C.Line(), myS->Cylinder());
128         break;
129       case GeomAbs_Plane:
130         myExtElCS.Perform(C.Line(), myS->Plane());
131         if (myExtElCS.IsParallel())   break;
132
133       case GeomAbs_Torus:
134       case GeomAbs_Cone:
135       case GeomAbs_BezierSurface:
136       case GeomAbs_BSplineSurface:
137       case GeomAbs_SurfaceOfRevolution:
138       case GeomAbs_SurfaceOfExtrusion:
139       case GeomAbs_OffsetSurface:
140       case GeomAbs_OtherSurface:
141         {
142           Standard_Real cfirst = myucinf, clast = myucsup;
143           Standard_Real ufirst = myS->FirstUParameter(), ulast = myS->LastUParameter(), 
144             vfirst = myS->FirstVParameter(), vlast = myS->LastVParameter();
145
146           if(Precision::IsInfinite(Abs(cfirst)) || Precision::IsInfinite(Abs(clast))) {
147
148             Bnd_Box aSurfBox;
149             BndLib_AddSurface::Add(*myS, ufirst, ulast, vfirst, vlast, Precision::Confusion(), aSurfBox);
150             Standard_Real xmin, ymin, zmin, xmax, ymax, zmax;
151             aSurfBox.Get(xmin, ymin, zmin, xmax, ymax, zmax);
152             Standard_Real tmin = Precision::Infinite(), tmax = -tmin;
153             gp_Lin aLin = C.Line();
154
155
156             if(!( Precision::IsInfinite(Abs(xmin)) || Precision::IsInfinite(Abs(xmax)) || 
157               Precision::IsInfinite(Abs(ymin)) || Precision::IsInfinite(Abs(ymax)) || 
158               Precision::IsInfinite(Abs(zmin)) || Precision::IsInfinite(Abs(zmax)))  ) {
159
160                 Extrema_ExtPElC anExt;
161                 Extrema_POnCurv aPntOnLin;
162                 Standard_Real aParOnLin;
163                 Standard_Real lim = Precision::Infinite();
164                 gp_Pnt aLimPntArray[8];
165
166                 aLimPntArray[0].SetCoord(xmin, ymin, zmin);
167                 aLimPntArray[1].SetCoord(xmax, ymin, zmin);
168                 aLimPntArray[2].SetCoord(xmin, ymax, zmin);
169                 aLimPntArray[3].SetCoord(xmax, ymax, zmin);
170                 aLimPntArray[4].SetCoord(xmin, ymin, zmax);
171                 aLimPntArray[5].SetCoord(xmax, ymin, zmax);
172                 aLimPntArray[6].SetCoord(xmin, ymax, zmax);
173                 aLimPntArray[7].SetCoord(xmax, ymax, zmax);
174
175                 for(i = 0; i <= 7; i++) {
176                   anExt.Perform(aLimPntArray[i], aLin, Precision::Confusion(), -lim, lim);
177                   aPntOnLin = anExt.Point(1);
178                   aParOnLin = aPntOnLin.Parameter();
179                   tmin = Min(aParOnLin, tmin);
180                   tmax = Max(aParOnLin, tmax);
181                 }
182
183             }
184             else {
185               tmin = -1.e+50;
186               tmax =  1.e+50;
187             }
188
189
190             cfirst = Max(cfirst, tmin);
191             clast  = Min(clast,  tmax);
192
193           }
194
195           if (myS->IsUPeriodic())
196             NbU = 13;
197           if (myS->IsVPeriodic())
198             NbV = 13;
199
200           Extrema_GenExtCS Ext(C, *myS, NbT, NbU, NbV, cfirst, clast, ufirst, ulast,
201             vfirst, vlast, mytolC, mytolS);
202
203           myDone = Ext.IsDone();
204           if (myDone) {
205             Standard_Integer NbExt = Ext.NbExt();
206             Standard_Real T,U,V;
207             Extrema_POnCurv PC;
208             Extrema_POnSurf PS;
209             for (i = 1; i <= NbExt; i++) {
210               PC = Ext.PointOnCurve(i);
211               PS = Ext.PointOnSurface(i);
212               T = PC.Parameter();
213               PS.Parameter(U, V);
214               AddSolution(C, T, U, V, PC.Value(), PS.Value(), Ext.SquareDistance(i));
215             }
216           }
217           return;
218
219         }
220       }
221       break;
222     }
223     //  Modified by skv - Thu Jul  7 12:29:34 2005 OCC9134 Begin
224   case GeomAbs_Circle:
225     {
226       if(myStype == GeomAbs_Cylinder) {
227         myExtElCS.Perform(C.Circle(), myS->Cylinder());
228         break;
229       }
230       else if(myStype == GeomAbs_Plane)
231       {
232         myExtElCS.Perform(C.Circle(), myS->Plane());
233         break;
234       }
235     }
236   case GeomAbs_Hyperbola: 
237     {
238       if(myCtype == GeomAbs_Hyperbola && myStype == GeomAbs_Plane) {
239         //  Modified by skv - Thu Jul  7 12:29:34 2005 OCC9134 End
240         myExtElCS.Perform(C.Hyperbola(), myS->Plane());
241         break;
242       }
243     }
244   default:
245     {
246       Extrema_GenExtCS Ext;
247       Ext.Initialize(*myS, NbU, NbV, mytolS);
248       if(myCtype == GeomAbs_Hyperbola) {
249         Standard_Real tmin = Max(-20., C.FirstParameter());
250         Standard_Real tmax = Min(20., C.LastParameter());
251         Ext.Perform(C, NbT, tmin, tmax, mytolC); // to avoid overflow
252       }
253       else {
254         if((myCtype == GeomAbs_Circle       && NbT < 13) ||
255            (myCtype == GeomAbs_BSplineCurve && NbT < 13) )
256         {
257           NbT = 13;
258         }
259         Ext.Perform(C, NbT, mytolC);
260       }
261
262       myDone = Ext.IsDone();
263       if (myDone) {
264         Standard_Integer NbExt = Ext.NbExt();
265         Standard_Real T,U,V;
266         Extrema_POnCurv PC;
267         Extrema_POnSurf PS;
268         for (i = 1; i <= NbExt; i++) {
269           PC = Ext.PointOnCurve(i);
270           PS = Ext.PointOnSurface(i);
271           T = PC.Parameter();
272           PS.Parameter(U, V);
273           AddSolution(C, T, U, V, PC.Value(), PS.Value(), Ext.SquareDistance(i));
274         }
275
276         //Add sharp points
277         Standard_Integer SolNumber = mySqDist.Length();
278         Standard_Address CopyC = (Standard_Address)&C;
279         Adaptor3d_Curve& aC = *(Adaptor3d_Curve*)CopyC;
280         Standard_Integer NbIntervals = aC.NbIntervals(GeomAbs_C1);
281         TColStd_Array1OfReal SharpPoints(1, NbIntervals+1);
282         aC.Intervals(SharpPoints, GeomAbs_C1);
283         
284         Extrema_ExtPS aProjPS;
285         aProjPS.Initialize (*myS,
286                             myS->FirstUParameter(),
287                             myS->LastUParameter(),
288                             myS->FirstVParameter(),
289                             myS->LastVParameter(),
290                             mytolS,
291                             mytolS);
292
293         for (i = 2; i < SharpPoints.Upper(); ++i)
294         {
295           T = SharpPoints(i);
296           gp_Pnt aPnt = C.Value(T);
297           aProjPS.Perform (aPnt);
298           if (!aProjPS.IsDone())
299             continue;
300           Standard_Integer NbProj = aProjPS.NbExt(), jmin = 0;
301           Standard_Real MinSqDist = RealLast();
302           for (j = 1; j <= NbProj; j++)
303           {
304             Standard_Real aSqDist = aProjPS.SquareDistance(j);
305             if (aSqDist < MinSqDist)
306             {
307               MinSqDist = aSqDist;
308               jmin = j;
309             }
310           }
311           if (jmin != 0)
312           {
313             aProjPS.Point(jmin).Parameter(U,V);
314             AddSolution(C, T, U, V,
315                         aPnt, aProjPS.Point(jmin).Value(), MinSqDist);
316           }
317         }
318         //Cut sharp solutions to keep only minimum and maximum
319         Standard_Integer imin = SolNumber + 1, imax = mySqDist.Length();
320         for (i = SolNumber + 1; i <= mySqDist.Length(); i++)
321         {
322           if (mySqDist(i) < mySqDist(imin))
323             imin = i;
324           if (mySqDist(i) > mySqDist(imax))
325             imax = i;
326         }
327         if (mySqDist.Length() > SolNumber + 2)
328         {
329           Standard_Real MinSqDist = mySqDist(imin);
330           Standard_Real MaxSqDist = mySqDist(imax);
331           Extrema_POnCurv MinPC = myPOnC(imin);
332           Extrema_POnCurv MaxPC = myPOnC(imax);
333           Extrema_POnSurf MinPS = myPOnS(imin);
334           Extrema_POnSurf MaxPS = myPOnS(imax);
335
336           mySqDist.Remove(SolNumber + 1, mySqDist.Length());
337           myPOnC.Remove(SolNumber + 1, myPOnC.Length());
338           myPOnS.Remove(SolNumber + 1, myPOnS.Length());
339
340           mySqDist.Append(MinSqDist);
341           myPOnC.Append(MinPC);
342           myPOnS.Append(MinPS);
343           mySqDist.Append(MaxSqDist);
344           myPOnC.Append(MaxPC);
345           myPOnS.Append(MaxPS);
346         }
347       }
348       return;
349     }
350     break;
351   }
352
353   myDone = myExtElCS.IsDone();
354   if (myDone) {
355     myIsPar = myExtElCS.IsParallel();
356     if (myIsPar) {
357       mySqDist.Append(myExtElCS.SquareDistance(1));
358     }
359     else {
360       Standard_Integer NbExt = myExtElCS.NbExt();
361       Standard_Real U, V;
362       for (i = 1; i <= NbExt; i++) {
363         Extrema_POnCurv PC;
364         Extrema_POnSurf PS;
365         myExtElCS.Points(i, PC, PS);
366         Standard_Real Ucurve = PC.Parameter();
367         PS.Parameter(U, V);
368         AddSolution(C, Ucurve, U, V, PC.Value(), PS.Value(), myExtElCS.SquareDistance(i));
369       }
370       if(mySqDist.Length() == 0 && NbExt > 0)
371       {
372         //Analytical extremas seem to be out of curve/surface boundaries.
373         //For plane it is possible to add extremity points of curve
374         if(myStype == GeomAbs_Plane)
375         {
376           gp_Pln aPln = myS->Plane();
377           gp_Pnt PC, PP;
378           if(!Precision::IsInfinite(myucinf))
379           {
380             PC = C.Value(myucinf);
381             ElSLib::PlaneParameters(aPln.Position(), PC, U, V);
382             PP = ElSLib::PlaneValue(U, V, aPln.Position());
383             AddSolution(C, myucinf, U, V, PC, PP, PC.SquareDistance(PP));
384           }
385           if(!Precision::IsInfinite(myucsup))
386           {
387             PC = C.Value(myucsup);
388             ElSLib::PlaneParameters(aPln.Position(), PC, U, V);
389             PP = ElSLib::PlaneValue(U, V, aPln.Position());
390             AddSolution(C, myucsup, U, V, PC, PP, PC.SquareDistance(PP));
391           }
392         }
393       }
394     }
395   }
396
397 }
398
399
400 Standard_Boolean Extrema_ExtCS::IsDone() const
401 {
402   return myDone;
403 }
404
405 Standard_Boolean Extrema_ExtCS::IsParallel() const
406 {
407   return myIsPar;
408 }
409
410
411 Standard_Real Extrema_ExtCS::SquareDistance(const Standard_Integer N) const
412 {
413   if(!myDone) throw StdFail_NotDone();
414   if (myIsPar && N != 1) throw StdFail_InfiniteSolutions();
415   if ((N < 1) || (N > mySqDist.Length())) throw Standard_OutOfRange();
416   return mySqDist.Value(N);
417 }
418
419
420 Standard_Integer Extrema_ExtCS::NbExt() const
421 {
422   if(!myDone) throw StdFail_NotDone();
423   return myPOnC.Length();
424 }
425
426
427
428 void Extrema_ExtCS::Points(const Standard_Integer N,
429   Extrema_POnCurv&       P1,
430   Extrema_POnSurf&       P2) const
431 {
432   if(!myDone) throw StdFail_NotDone();
433   P1 = myPOnC.Value(N);
434   P2 = myPOnS.Value(N);
435 }
436
437 Standard_Boolean Extrema_ExtCS::AddSolution(const Adaptor3d_Curve& theCurve,
438   const Standard_Real aT,
439   const Standard_Real aU,
440   const Standard_Real aV,
441   const gp_Pnt& PointOnCurve,
442   const gp_Pnt& PointOnSurf,
443   const Standard_Real SquareDist)
444 {
445   Standard_Boolean Added = Standard_False;
446
447   Standard_Real T = aT, U = aU, V = aV;
448
449   if (theCurve.IsPeriodic())
450     T = ElCLib::InPeriod(T, myucinf, myucinf + theCurve.Period());
451   if (myS->IsUPeriodic())
452     U = ElCLib::InPeriod(U, myuinf, myuinf + myS->UPeriod());
453   if (myS->IsVPeriodic())
454     V = ElCLib::InPeriod(V, myvinf, myvinf + myS->VPeriod());
455
456   Extrema_POnCurv aPC;
457   Extrema_POnSurf aPS;
458   if ((myucinf-T) <= mytolC && (T-myucsup) <= mytolC &&
459     (myuinf-U) <= mytolS && (U-myusup) <= mytolS &&
460     (myvinf-V) <= mytolS && (V-myvsup) <= mytolS)
461   {
462     Standard_Boolean IsNewSolution = Standard_True;
463     for (Standard_Integer j = 1; j <= mySqDist.Length(); j++)
464     {
465       aPC = myPOnC(j);
466       aPS = myPOnS(j);
467       Standard_Real Tj = aPC.Parameter();
468       Standard_Real Uj, Vj;
469       aPS.Parameter(Uj, Vj);
470       if (Abs(T - Tj) <= mytolC &&
471         Abs(U - Uj) <= mytolS &&
472         Abs(V - Vj) <= mytolS)
473       {
474         IsNewSolution = Standard_False;
475         break;
476       }
477     }
478     if (IsNewSolution)
479     {
480       mySqDist.Append(SquareDist);
481       aPC.SetValues(T, PointOnCurve);
482       myPOnC.Append(aPC);
483       myPOnS.Append(Extrema_POnSurf(U, V, PointOnSurf));
484       Added = Standard_True;
485     }
486   }
487   return Added;
488 }