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