0028599: Replacement of old Boolean operations with new ones in BRepProj_Projection...
[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   myDone = Standard_False;
117   // Try analytic computation of extrema
118   Standard_Boolean isComputeAnalytic = Standard_True;
119
120   switch(myCtype) {
121
122   case GeomAbs_Line: 
123     {
124
125       switch(myStype) {
126       case GeomAbs_Sphere:
127         myExtElCS.Perform(C.Line(), myS->Sphere());
128         break;   
129       case GeomAbs_Cylinder:
130         myExtElCS.Perform(C.Line(), myS->Cylinder());
131         break;
132       case GeomAbs_Plane:
133         myExtElCS.Perform(C.Line(), myS->Plane());
134         if (myExtElCS.IsParallel())   break;
135         Standard_FALLTHROUGH
136
137       case GeomAbs_Torus:
138       case GeomAbs_Cone:
139       case GeomAbs_BezierSurface:
140       case GeomAbs_BSplineSurface:
141       case GeomAbs_SurfaceOfRevolution:
142       case GeomAbs_SurfaceOfExtrusion:
143       case GeomAbs_OffsetSurface:
144       case GeomAbs_OtherSurface:
145         {
146           Standard_Real cfirst = myucinf, clast = myucsup;
147           Standard_Real ufirst = myS->FirstUParameter(), ulast = myS->LastUParameter(), 
148             vfirst = myS->FirstVParameter(), vlast = myS->LastVParameter();
149
150           if(Precision::IsInfinite(Abs(cfirst)) || Precision::IsInfinite(Abs(clast))) {
151
152             Bnd_Box aSurfBox;
153             BndLib_AddSurface::Add(*myS, ufirst, ulast, vfirst, vlast, Precision::Confusion(), aSurfBox);
154             Standard_Real xmin, ymin, zmin, xmax, ymax, zmax;
155             aSurfBox.Get(xmin, ymin, zmin, xmax, ymax, zmax);
156             Standard_Real tmin = Precision::Infinite(), tmax = -tmin;
157             gp_Lin aLin = C.Line();
158
159
160             if(!( Precision::IsInfinite(Abs(xmin)) || Precision::IsInfinite(Abs(xmax)) || 
161               Precision::IsInfinite(Abs(ymin)) || Precision::IsInfinite(Abs(ymax)) || 
162               Precision::IsInfinite(Abs(zmin)) || Precision::IsInfinite(Abs(zmax)))  ) {
163
164                 Extrema_ExtPElC anExt;
165                 Extrema_POnCurv aPntOnLin;
166                 Standard_Real aParOnLin;
167                 Standard_Real lim = Precision::Infinite();
168                 gp_Pnt aLimPntArray[8];
169
170                 aLimPntArray[0].SetCoord(xmin, ymin, zmin);
171                 aLimPntArray[1].SetCoord(xmax, ymin, zmin);
172                 aLimPntArray[2].SetCoord(xmin, ymax, zmin);
173                 aLimPntArray[3].SetCoord(xmax, ymax, zmin);
174                 aLimPntArray[4].SetCoord(xmin, ymin, zmax);
175                 aLimPntArray[5].SetCoord(xmax, ymin, zmax);
176                 aLimPntArray[6].SetCoord(xmin, ymax, zmax);
177                 aLimPntArray[7].SetCoord(xmax, ymax, zmax);
178
179                 for(i = 0; i <= 7; i++) {
180                   anExt.Perform(aLimPntArray[i], aLin, Precision::Confusion(), -lim, lim);
181                   aPntOnLin = anExt.Point(1);
182                   aParOnLin = aPntOnLin.Parameter();
183                   tmin = Min(aParOnLin, tmin);
184                   tmax = Max(aParOnLin, tmax);
185                 }
186
187             }
188             else {
189               tmin = -1.e+50;
190               tmax =  1.e+50;
191             }
192
193
194             cfirst = Max(cfirst, tmin);
195             clast  = Min(clast,  tmax);
196
197           }
198
199           if (myS->IsUPeriodic())
200             NbU = 13;
201           if (myS->IsVPeriodic())
202             NbV = 13;
203
204           Extrema_GenExtCS Ext(C, *myS, NbT, NbU, NbV, cfirst, clast, ufirst, ulast,
205             vfirst, vlast, mytolC, mytolS);
206
207           myDone = Ext.IsDone();
208           if (myDone) {
209             Standard_Integer NbExt = Ext.NbExt();
210             Standard_Real T,U,V;
211             Extrema_POnCurv PC;
212             Extrema_POnSurf PS;
213             for (i = 1; i <= NbExt; i++) {
214               PC = Ext.PointOnCurve(i);
215               PS = Ext.PointOnSurface(i);
216               T = PC.Parameter();
217               PS.Parameter(U, V);
218               AddSolution(C, T, U, V, PC.Value(), PS.Value(), Ext.SquareDistance(i));
219             }
220           }
221           return;
222
223         }
224       }
225       break;
226     }
227     //  Modified by skv - Thu Jul  7 12:29:34 2005 OCC9134 Begin
228   case GeomAbs_Circle:
229     {
230       if(myStype == GeomAbs_Cylinder) {
231         myExtElCS.Perform(C.Circle(), myS->Cylinder());
232         break;
233       }
234       else if(myStype == GeomAbs_Plane)
235       {
236         myExtElCS.Perform(C.Circle(), myS->Plane());
237         break;
238       }
239       else if (myStype == GeomAbs_Sphere)
240       {
241         myExtElCS.Perform(C.Circle(), myS->Sphere());
242         break;
243       }
244     }
245     Standard_FALLTHROUGH
246   case GeomAbs_Hyperbola: 
247     {
248       if(myCtype == GeomAbs_Hyperbola && myStype == GeomAbs_Plane) {
249         //  Modified by skv - Thu Jul  7 12:29:34 2005 OCC9134 End
250         myExtElCS.Perform(C.Hyperbola(), myS->Plane());
251         break;
252       }
253     }
254     Standard_FALLTHROUGH
255   default:
256     {
257       isComputeAnalytic = Standard_False;
258       break;
259     }
260   }
261
262   if (isComputeAnalytic)
263   {
264     if (myExtElCS.IsDone())
265     {
266       myDone = Standard_True;
267       myIsPar = myExtElCS.IsParallel();
268       if (myIsPar)
269       {
270         mySqDist.Append(myExtElCS.SquareDistance(1));
271       }
272       else
273       {
274         Standard_Integer NbExt = myExtElCS.NbExt();
275         for (i = 1; i <= NbExt; i++)
276         {
277           Extrema_POnCurv PC;
278           Extrema_POnSurf PS;
279           myExtElCS.Points(i, PC, PS);
280           Standard_Real Ucurve = PC.Parameter();
281           Standard_Real U, V;
282           PS.Parameter(U, V);
283           AddSolution(C, Ucurve, U, V, PC.Value(), PS.Value(), myExtElCS.SquareDistance(i));
284         }
285
286         if (mySqDist.Length() == 0 && NbExt > 0)
287         {
288           // Analytical extrema seem to be out of curve/surface boundaries.
289           // Try extremity points of curve.
290           gp_Pnt aPOnC[2], aPOnS[2];
291           Standard_Real aT[2] = { myucinf, myucsup }, U[2], V[2];
292           Standard_Real aDist[2] = { -1, -1 };
293           for (i = 0; i < 2; ++i)
294           {
295             if (Precision::IsInfinite(aT[i]))
296               continue;
297
298             aPOnC[i] = C.Value(aT[i]);
299             switch (myStype)
300             {
301               case GeomAbs_Plane:
302               {
303                 ElSLib::Parameters(myS->Plane(), aPOnC[i], U[i], V[i]);
304                 aPOnS[i] = ElSLib::Value(U[i], V[i], myS->Plane());
305                 break;
306               }
307               case GeomAbs_Sphere:
308               {
309                 ElSLib::Parameters(myS->Sphere(), aPOnC[i], U[i], V[i]);
310                 aPOnS[i] = ElSLib::Value(U[i], V[i], myS->Sphere());
311                 break;
312               }
313               case GeomAbs_Cylinder:
314               {
315                 ElSLib::Parameters(myS->Cylinder(), aPOnC[i], U[i], V[i]);
316                 aPOnS[i] = ElSLib::Value(U[i], V[i], myS->Cylinder());
317                 break;
318               }
319               case GeomAbs_Torus:
320               {
321                 ElSLib::Parameters(myS->Torus(), aPOnC[i], U[i], V[i]);
322                 aPOnS[i] = ElSLib::Value(U[i], V[i], myS->Torus());
323                 break;
324               }
325               case GeomAbs_Cone:
326               {
327                 ElSLib::Parameters(myS->Cone(), aPOnC[i], U[i], V[i]);
328                 aPOnS[i] = ElSLib::Value(U[i], V[i], myS->Cone());
329                 break;
330               }
331               default:
332                 continue;
333             }
334
335             aDist[i] = aPOnC[i].SquareDistance(aPOnS[i]);
336           }
337
338           Standard_Boolean bAdd[2] = {Standard_False, Standard_False};
339
340           // Choose solution to add
341           if (aDist[0] >= 0. && aDist[1] >= 0.)
342           {
343             Standard_Real aDiff = aDist[0] - aDist[1];
344             // Both computed -> take only minimal
345             if (Abs(aDiff) < Precision::Confusion())
346               // Add both
347               bAdd[0] = bAdd[1] = Standard_True;
348             else if (aDiff < 0)
349               // Add first
350               bAdd[0] = Standard_True;
351             else
352               // Add second
353               bAdd[1] = Standard_True;
354           }
355           else if (aDist[0] >= 0.)
356             // Add first
357             bAdd[0] = Standard_True;
358           else if (aDist[1] >= 0.)
359             // Add second
360             bAdd[1] = Standard_True;
361
362           for (i = 0; i < 2; ++i)
363           {
364             if (bAdd[i])
365               AddSolution(C, aT[i], U[i], V[i], aPOnC[i], aPOnS[i], aDist[i]);
366           }
367         }
368       }
369       return;
370     }
371   }
372
373   // Elementary extrema is not done, try generic solution
374   Extrema_GenExtCS Ext;
375   Ext.Initialize(*myS, NbU, NbV, mytolS);
376   if (myCtype == GeomAbs_Hyperbola) {
377     Standard_Real tmin = Max(-20., C.FirstParameter());
378     Standard_Real tmax = Min(20., C.LastParameter());
379     Ext.Perform(C, NbT, tmin, tmax, mytolC); // to avoid overflow
380   }
381   else {
382     if ((myCtype == GeomAbs_Circle       && NbT < 13) ||
383       (myCtype == GeomAbs_BSplineCurve && NbT < 13))
384     {
385       NbT = 13;
386     }
387     Ext.Perform(C, NbT, mytolC);
388   }
389
390   myDone = Ext.IsDone();
391   if (myDone) {
392     Standard_Integer NbExt = Ext.NbExt();
393     Standard_Real T, U, V;
394     Extrema_POnCurv PC;
395     Extrema_POnSurf PS;
396     for (i = 1; i <= NbExt; i++) {
397       PC = Ext.PointOnCurve(i);
398       PS = Ext.PointOnSurface(i);
399       T = PC.Parameter();
400       PS.Parameter(U, V);
401       AddSolution(C, T, U, V, PC.Value(), PS.Value(), Ext.SquareDistance(i));
402     }
403
404     //Add sharp points
405     Standard_Integer SolNumber = mySqDist.Length();
406     Standard_Address CopyC = (Standard_Address)&C;
407     Adaptor3d_Curve& aC = *(Adaptor3d_Curve*)CopyC;
408     Standard_Integer NbIntervals = aC.NbIntervals(GeomAbs_C1);
409     TColStd_Array1OfReal SharpPoints(1, NbIntervals + 1);
410     aC.Intervals(SharpPoints, GeomAbs_C1);
411
412     Extrema_ExtPS aProjPS;
413     aProjPS.Initialize(*myS,
414       myS->FirstUParameter(),
415       myS->LastUParameter(),
416       myS->FirstVParameter(),
417       myS->LastVParameter(),
418       mytolS,
419       mytolS);
420
421     for (i = 2; i < SharpPoints.Upper(); ++i)
422     {
423       T = SharpPoints(i);
424       gp_Pnt aPnt = C.Value(T);
425       aProjPS.Perform(aPnt);
426       if (!aProjPS.IsDone())
427         continue;
428       Standard_Integer NbProj = aProjPS.NbExt(), jmin = 0;
429       Standard_Real MinSqDist = RealLast();
430       for (j = 1; j <= NbProj; j++)
431       {
432         Standard_Real aSqDist = aProjPS.SquareDistance(j);
433         if (aSqDist < MinSqDist)
434         {
435           MinSqDist = aSqDist;
436           jmin = j;
437         }
438       }
439       if (jmin != 0)
440       {
441         aProjPS.Point(jmin).Parameter(U, V);
442         AddSolution(C, T, U, V,
443           aPnt, aProjPS.Point(jmin).Value(), MinSqDist);
444       }
445     }
446     //Cut sharp solutions to keep only minimum and maximum
447     Standard_Integer imin = SolNumber + 1, imax = mySqDist.Length();
448     for (i = SolNumber + 1; i <= mySqDist.Length(); i++)
449     {
450       if (mySqDist(i) < mySqDist(imin))
451         imin = i;
452       if (mySqDist(i) > mySqDist(imax))
453         imax = i;
454     }
455     if (mySqDist.Length() > SolNumber + 2)
456     {
457       Standard_Real MinSqDist = mySqDist(imin);
458       Standard_Real MaxSqDist = mySqDist(imax);
459       Extrema_POnCurv MinPC = myPOnC(imin);
460       Extrema_POnCurv MaxPC = myPOnC(imax);
461       Extrema_POnSurf MinPS = myPOnS(imin);
462       Extrema_POnSurf MaxPS = myPOnS(imax);
463
464       mySqDist.Remove(SolNumber + 1, mySqDist.Length());
465       myPOnC.Remove(SolNumber + 1, myPOnC.Length());
466       myPOnS.Remove(SolNumber + 1, myPOnS.Length());
467
468       mySqDist.Append(MinSqDist);
469       myPOnC.Append(MinPC);
470       myPOnS.Append(MinPS);
471       mySqDist.Append(MaxSqDist);
472       myPOnC.Append(MaxPC);
473       myPOnS.Append(MaxPS);
474     }
475   }
476 }
477
478
479 Standard_Boolean Extrema_ExtCS::IsDone() const
480 {
481   return myDone;
482 }
483
484 Standard_Boolean Extrema_ExtCS::IsParallel() const
485 {
486   return myIsPar;
487 }
488
489
490 Standard_Real Extrema_ExtCS::SquareDistance(const Standard_Integer N) const
491 {
492   if(!myDone) throw StdFail_NotDone();
493   if (myIsPar && N != 1) throw StdFail_InfiniteSolutions();
494   if ((N < 1) || (N > mySqDist.Length())) throw Standard_OutOfRange();
495   return mySqDist.Value(N);
496 }
497
498
499 Standard_Integer Extrema_ExtCS::NbExt() const
500 {
501   if(!myDone) throw StdFail_NotDone();
502   return myPOnC.Length();
503 }
504
505
506
507 void Extrema_ExtCS::Points(const Standard_Integer N,
508   Extrema_POnCurv&       P1,
509   Extrema_POnSurf&       P2) const
510 {
511   if(!myDone) throw StdFail_NotDone();
512   P1 = myPOnC.Value(N);
513   P2 = myPOnS.Value(N);
514 }
515
516 Standard_Boolean Extrema_ExtCS::AddSolution(const Adaptor3d_Curve& theCurve,
517   const Standard_Real aT,
518   const Standard_Real aU,
519   const Standard_Real aV,
520   const gp_Pnt& PointOnCurve,
521   const gp_Pnt& PointOnSurf,
522   const Standard_Real SquareDist)
523 {
524   Standard_Boolean Added = Standard_False;
525
526   Standard_Real T = aT, U = aU, V = aV;
527
528   if (theCurve.IsPeriodic())
529     T = ElCLib::InPeriod(T, myucinf, myucinf + theCurve.Period());
530   if (myS->IsUPeriodic())
531     U = ElCLib::InPeriod(U, myuinf, myuinf + myS->UPeriod());
532   if (myS->IsVPeriodic())
533     V = ElCLib::InPeriod(V, myvinf, myvinf + myS->VPeriod());
534
535   Extrema_POnCurv aPC;
536   Extrema_POnSurf aPS;
537   if ((myucinf-T) <= mytolC && (T-myucsup) <= mytolC &&
538     (myuinf-U) <= mytolS && (U-myusup) <= mytolS &&
539     (myvinf-V) <= mytolS && (V-myvsup) <= mytolS)
540   {
541     Standard_Boolean IsNewSolution = Standard_True;
542     for (Standard_Integer j = 1; j <= mySqDist.Length(); j++)
543     {
544       aPC = myPOnC(j);
545       aPS = myPOnS(j);
546       Standard_Real Tj = aPC.Parameter();
547       Standard_Real Uj, Vj;
548       aPS.Parameter(Uj, Vj);
549       if (Abs(T - Tj) <= mytolC &&
550         Abs(U - Uj) <= mytolS &&
551         Abs(V - Vj) <= mytolS)
552       {
553         IsNewSolution = Standard_False;
554         break;
555       }
556     }
557     if (IsNewSolution)
558     {
559       mySqDist.Append(SquareDist);
560       aPC.SetValues(T, PointOnCurve);
561       myPOnC.Append(aPC);
562       myPOnS.Append(Extrema_POnSurf(U, V, PointOnSurf));
563       Added = Standard_True;
564     }
565   }
566   return Added;
567 }