0027261: Incorrect bounding boxes computed for the b-spline faces
[occt.git] / src / BndLib / BndLib_AddSurface.cxx
1 // Created on: 1995-07-24
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 - Fri Aug 27 12:29:04 2004 OCC6503
18
19 #include <Adaptor3d_HSurface.hxx>
20 #include <Adaptor3d_Surface.hxx>
21 #include <Bnd_Box.hxx>
22 #include <BndLib.hxx>
23 #include <BndLib_AddSurface.hxx>
24 #include <BSplCLib.hxx>
25 #include <ElSLib.hxx>
26 #include <Geom_BezierSurface.hxx>
27 #include <Geom_BSplineSurface.hxx>
28 #include <GeomAbs_SurfaceType.hxx>
29 #include <gp_Pln.hxx>
30 #include <gp_Pnt.hxx>
31 #include <Precision.hxx>
32 #include <TColgp_Array2OfPnt.hxx>
33 #include <TColStd_Array1OfInteger.hxx>
34 #include <TColStd_Array1OfReal.hxx>
35
36 //=======================================================================
37 //function : Add
38 //purpose  : 
39 //=======================================================================
40 void BndLib_AddSurface::Add(const Adaptor3d_Surface& S,
41                             const Standard_Real Tol,
42                             Bnd_Box& B ) 
43 {
44
45   BndLib_AddSurface::Add(S,
46                          S.FirstUParameter(),
47                          S.LastUParameter (),
48                          S.FirstVParameter(),
49                          S.LastVParameter (),Tol,B);
50 }
51 //=======================================================================
52 //function : NbUSamples
53 //purpose  : 
54 //=======================================================================
55
56 static Standard_Integer NbUSamples(const Adaptor3d_Surface& S) 
57 {
58   Standard_Integer N;
59   GeomAbs_SurfaceType Type = S.GetType();
60   switch (Type) {
61   case GeomAbs_BezierSurface: 
62     {
63       N = 2*S.NbUPoles();
64       break;
65     }
66   case GeomAbs_BSplineSurface: 
67     {
68       const Handle(Geom_BSplineSurface)& BS = S.BSpline();
69       N = 2*(BS->UDegree() + 1)*(BS->NbUKnots() -1);
70       break;
71     }
72   default:
73     N = 33;
74   }
75   return Min (50,N);
76 }
77
78 //=======================================================================
79 //function : NbVSamples
80 //purpose  : 
81 //=======================================================================
82
83 static Standard_Integer NbVSamples(const Adaptor3d_Surface& S) 
84 {
85   Standard_Integer N;
86   GeomAbs_SurfaceType Type = S.GetType();
87   switch (Type) {
88   case GeomAbs_BezierSurface:
89     {
90       N = 2*S.NbVPoles();
91       break;
92     }
93   case GeomAbs_BSplineSurface:
94     {
95       const Handle(Geom_BSplineSurface)& BS = S.BSpline();
96       N = 2*(BS->VDegree() + 1)*(BS->NbVKnots() - 1) ;
97       break;
98     }
99   default:
100     N = 33;
101   }
102   return Min(50,N);
103 }
104
105 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 Begin
106 static gp_Pnt BaryCenter(const gp_Pln        &aPlane,
107                          const Standard_Real  aUMin,
108                          const Standard_Real  aUMax,
109                          const Standard_Real  aVMin,
110                          const Standard_Real  aVMax)
111 {
112   Standard_Real aU, aV;
113   Standard_Boolean isU1Inf = Precision::IsInfinite(aUMin);
114   Standard_Boolean isU2Inf = Precision::IsInfinite(aUMax);
115   Standard_Boolean isV1Inf = Precision::IsInfinite(aVMin);
116   Standard_Boolean isV2Inf = Precision::IsInfinite(aVMax);
117
118   if (isU1Inf && isU2Inf)
119     aU = 0;
120   else if (isU1Inf)
121     aU = aUMax - 10.;
122   else if (isU2Inf)
123     aU = aUMin + 10.;
124   else
125     aU = (aUMin + aUMax)/2.;
126     
127   if (isV1Inf && isV2Inf)
128     aV = 0;
129   else if (isV1Inf)
130     aV = aVMax - 10.;
131   else if (isV2Inf)
132     aV = aVMin + 10.;
133   else
134     aV = (aVMin + aVMax)/2.;
135
136   gp_Pnt aCenter = ElSLib::Value(aU, aV, aPlane);
137
138   return aCenter;
139 }
140
141 static void TreatInfinitePlane(const gp_Pln        &aPlane,
142                                const Standard_Real  aUMin,
143                                const Standard_Real  aUMax,
144                                const Standard_Real  aVMin,
145                                const Standard_Real  aVMax,
146                                const Standard_Real  aTol,
147                                      Bnd_Box       &aB)
148 {
149   // Get 3 coordinate axes of the plane.
150   const gp_Dir        &aNorm        = aPlane.Axis().Direction();
151   const Standard_Real  anAngularTol = RealEpsilon();
152
153   // Get location of the plane as its barycenter
154   gp_Pnt aLocation = BaryCenter(aPlane, aUMin, aUMax, aVMin, aVMax);
155
156   if (aNorm.IsParallel(gp::DX(), anAngularTol)) {
157     aB.Add(aLocation);
158     aB.OpenYmin();
159     aB.OpenYmax();
160     aB.OpenZmin();
161     aB.OpenZmax();
162   } else if (aNorm.IsParallel(gp::DY(), anAngularTol)) {
163     aB.Add(aLocation);
164     aB.OpenXmin();
165     aB.OpenXmax();
166     aB.OpenZmin();
167     aB.OpenZmax();
168   } else if (aNorm.IsParallel(gp::DZ(), anAngularTol)) {
169     aB.Add(aLocation);
170     aB.OpenXmin();
171     aB.OpenXmax();
172     aB.OpenYmin();
173     aB.OpenYmax();
174   } else {
175     aB.SetWhole();
176     return;
177   }
178
179   aB.Enlarge(aTol);
180 }
181
182 // Compute start and finish indexes used in convex hull.
183 // theMinIdx - minimum poles index, that can be used.
184 // theMaxIdx - maximum poles index, that can be used.
185 // theShiftCoeff - shift between flatknots array and poles array.
186 // This vaule should be equal to 1 in case of non periodic BSpline,
187 // and (degree + 1) - mults(the lowest index).
188 void ComputePolesIndexes(const TColStd_Array1OfReal &theFlatKnots,
189                          const Standard_Integer theDegree,
190                          const Standard_Real theMin,
191                          const Standard_Real theMax,
192                          const Standard_Integer theMinIdx,
193                          const Standard_Integer theMaxIdx,
194                          const Standard_Integer theShiftCoeff,
195                          Standard_Integer &theOutMinIdx,
196                          Standard_Integer &theOutMaxIdx)
197 {
198   // Set initial values for the result indexes to handle situation when requested parameter space
199   // is slightly greater than B-spline parameter space.
200   theOutMinIdx = theFlatKnots.Lower();
201   theOutMaxIdx = theFlatKnots.Upper();
202
203   // Compute first and last used flat knots.
204   for(Standard_Integer aKnotIdx = theFlatKnots.Lower();
205       aKnotIdx < theFlatKnots.Upper();
206       aKnotIdx++)
207   {
208     if (theFlatKnots(aKnotIdx) <= theMin)
209       theOutMinIdx = aKnotIdx;
210
211     if (theFlatKnots(theFlatKnots.Upper() - aKnotIdx + theFlatKnots.Lower()) >= theMax)
212       theOutMaxIdx = theFlatKnots.Upper() - aKnotIdx + theFlatKnots.Lower();
213   }
214
215   theOutMinIdx = Max(theOutMinIdx - 2 * theDegree + 2 - theShiftCoeff, theMinIdx);
216   theOutMaxIdx = Min(theOutMaxIdx - 2 + theDegree + 1 - theShiftCoeff, theMaxIdx);
217 }
218
219 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 End
220 //=======================================================================
221 //function : Add
222 //purpose  : 
223 //=======================================================================
224 void BndLib_AddSurface::Add(const Adaptor3d_Surface& S,
225                             const Standard_Real UMin,
226                             const Standard_Real UMax,
227                             const Standard_Real VMin,
228                             const Standard_Real VMax,
229                             const Standard_Real Tol,
230                             Bnd_Box& B ) 
231 {
232   GeomAbs_SurfaceType Type = S.GetType(); // skv OCC6503
233
234   if (Precision::IsInfinite(VMin) ||
235       Precision::IsInfinite(VMax) ||
236       Precision::IsInfinite(UMin) ||
237       Precision::IsInfinite(UMax)   ) {
238 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 Begin
239 //     B.SetWhole();
240 //     return;
241     switch (Type) {
242     case GeomAbs_Plane: 
243       {
244         TreatInfinitePlane(S.Plane(), UMin, UMax, VMin, VMax, Tol, B);
245         return;
246       }
247     default: 
248       {
249         B.SetWhole();
250         return;
251       }
252     }
253 //  Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503 End
254   }
255
256 //   GeomAbs_SurfaceType Type = S.GetType(); // skv OCC6503
257
258   switch (Type) {
259
260   case GeomAbs_Plane: 
261     {
262       gp_Pln Plan = S.Plane();
263       B.Add(ElSLib::Value(UMin,VMin,Plan)); 
264       B.Add(ElSLib::Value(UMin,VMax,Plan)); 
265       B.Add(ElSLib::Value(UMax,VMin,Plan)); 
266       B.Add(ElSLib::Value(UMax,VMax,Plan)); 
267       B.Enlarge(Tol);
268       break;
269     }
270   case GeomAbs_Cylinder: 
271     {
272       BndLib::Add(S.Cylinder(),UMin,UMax,VMin,VMax,Tol,B);
273       break;
274     }
275   case GeomAbs_Cone: 
276     {
277       BndLib::Add(S.Cone(),UMin,UMax,VMin,VMax,Tol,B);
278       break;
279     }
280   case GeomAbs_Torus: 
281     {
282       BndLib::Add(S.Torus(),UMin,UMax,VMin,VMax,Tol,B);
283       break;
284     }
285   case GeomAbs_Sphere: 
286     {
287       if (Abs(UMin) < Precision::Angular() &&
288           Abs(UMax - 2.*M_PI) < Precision::Angular() &&
289           Abs(VMin + M_PI/2.) < Precision::Angular() &&
290           Abs(VMax - M_PI/2.) < Precision::Angular()) // a whole sphere
291         BndLib::Add(S.Sphere(),Tol,B);
292       else
293         BndLib::Add(S.Sphere(),UMin,UMax,VMin,VMax,Tol,B);
294       break;
295     }
296   case GeomAbs_OffsetSurface: 
297     {
298       Handle(Adaptor3d_HSurface) HS = S.BasisSurface();
299       Add (HS->Surface(),UMin,UMax,VMin,VMax,Tol,B);
300       B.Enlarge(S.OffsetValue());
301       B.Enlarge(Tol);
302       break;
303     } 
304   case GeomAbs_BezierSurface:
305   case GeomAbs_BSplineSurface: 
306     {
307       Standard_Boolean isUseConvexHullAlgorithm = Standard_True;
308       Standard_Real PTol = Precision::Parametric(Precision::Confusion());
309       // Borders of underlying geometry.
310       Standard_Real anUMinParam = UMin, anUMaxParam = UMax,// BSpline case.
311                      aVMinParam = VMin,  aVMaxParam = VMax;
312       if (Type == GeomAbs_BezierSurface)
313       {
314         // Bezier surface:
315         // All of poles used for any parameter,
316         // thats why in case of trimmed parameters handled by grid algorithm.
317
318         if (Abs(UMin-S.FirstUParameter()) > PTol ||
319             Abs(VMin-S.FirstVParameter()) > PTol ||
320             Abs(UMax-S.LastUParameter ()) > PTol ||
321             Abs(VMax-S.LastVParameter ()) > PTol )
322         {
323           // Borders not equal to topology borders.
324           isUseConvexHullAlgorithm = Standard_False;
325         }
326       }
327       else
328       {
329         // BSpline:
330         // If Umin, Vmin, Umax, Vmax lies inside geometry bounds then:
331         // use convex hull algorithm,
332         // if Umin, VMin, Umax, Vmax lies outside then:
333         // use grid algorithm on analytic continuation (default case).
334         S.BSpline()->Bounds(anUMinParam, anUMaxParam, aVMinParam, aVMaxParam);
335
336         if ( (UMin - anUMinParam) < -PTol ||
337              (VMin -  aVMinParam) < -PTol ||
338              (UMax - anUMaxParam) >  PTol ||
339              (VMax -  aVMaxParam) >  PTol )
340         {
341           // Out of geometry borders.
342           isUseConvexHullAlgorithm = Standard_False;
343         }
344       }
345
346       if (isUseConvexHullAlgorithm)
347       {
348           TColgp_Array2OfPnt Tp(1,S.NbUPoles(),1,S.NbVPoles());
349           Standard_Integer UMinIdx = 0, UMaxIdx = 0;
350           Standard_Integer VMinIdx = 0, VMaxIdx = 0;
351           if (Type == GeomAbs_BezierSurface)
352           {
353             S.Bezier()->Poles(Tp);
354
355             UMinIdx = Tp.LowerRow();
356             UMaxIdx = Tp.UpperRow();
357             VMinIdx = Tp.LowerCol();
358             VMaxIdx = Tp.UpperCol();
359           }
360           else
361           {
362             S.BSpline()->Poles(Tp);
363
364             UMinIdx = Tp.LowerRow();
365             UMaxIdx = Tp.UpperRow();
366             VMinIdx = Tp.LowerCol();
367             VMaxIdx = Tp.UpperCol();
368
369             if (UMin > anUMinParam ||
370                 UMax < anUMaxParam)
371             {
372               Standard_Integer anUFlatKnotsCount = S.BSpline()->NbUPoles() + S.BSpline()->UDegree() + 1;
373               Standard_Integer aShift = 1;
374
375               if (S.BSpline()->IsUPeriodic())
376               {
377                 TColStd_Array1OfInteger aMults(1, S.BSpline()->NbUKnots());
378                 S.BSpline()->UMultiplicities(aMults);
379                 anUFlatKnotsCount = BSplCLib::KnotSequenceLength(aMults, S.BSpline()->UDegree(), Standard_True);
380
381                 aShift = S.BSpline()->UDegree() + 1 - S.BSpline()->UMultiplicity(1);
382               }
383
384               TColStd_Array1OfReal anUFlatKnots(1, anUFlatKnotsCount);
385               S.BSpline()->UKnotSequence(anUFlatKnots);
386
387               ComputePolesIndexes(anUFlatKnots,
388                                   S.BSpline()->UDegree(),
389                                   UMin, UMax,
390                                   UMinIdx, UMaxIdx,  // Min and Max Indexes
391                                   aShift,
392                                   UMinIdx, UMaxIdx); // the Output indexes
393             }
394
395             if (VMin > aVMinParam ||
396                 VMax < aVMaxParam)
397             {
398               Standard_Integer anVFlatKnotsCount = S.BSpline()->NbVPoles() + S.BSpline()->VDegree() + 1;
399               Standard_Integer aShift = 1;
400
401               if (S.BSpline()->IsVPeriodic())
402               {
403                 TColStd_Array1OfInteger aMults(1, S.BSpline()->NbVKnots());
404                 S.BSpline()->VMultiplicities(aMults);
405                 anVFlatKnotsCount = BSplCLib::KnotSequenceLength(aMults, S.BSpline()->VDegree(), Standard_True);
406
407                 aShift = S.BSpline()->VDegree() + 1 - S.BSpline()->VMultiplicity(1);
408               }
409
410               TColStd_Array1OfReal anVFlatKnots(1, anVFlatKnotsCount);
411               S.BSpline()->VKnotSequence(anVFlatKnots);
412
413               ComputePolesIndexes(anVFlatKnots,
414                                   S.BSpline()->VDegree(),
415                                   VMin, VMax,
416                                   VMinIdx, VMaxIdx,  // Min and Max Indexes
417                                   aShift,
418                                   VMinIdx, VMaxIdx); // the Output indexes
419             }
420
421           }
422
423           // Use poles to build convex hull.
424           for (Standard_Integer i = UMinIdx; i <= UMaxIdx; i++)
425           {
426             for (Standard_Integer j = VMinIdx; j <= VMaxIdx; j++)
427             {
428               B.Add(Tp(i,j));
429             }
430           }
431
432           B.Enlarge(Tol);
433           break;
434       }
435     }
436   default: 
437     {
438       Standard_Integer Nu = NbUSamples(S);
439       Standard_Integer Nv = NbVSamples(S);
440       gp_Pnt P;
441       for (Standard_Integer i =1 ;i<=Nu;i++){
442         Standard_Real U = UMin + ((UMax-UMin)*(i-1)/(Nu-1));
443         for (Standard_Integer j=1 ;j<=Nv;j++){
444           Standard_Real V = VMin + ((VMax-VMin)*(j-1)/(Nv-1));
445           S.D0(U,V,P);
446           B.Add(P);
447         }
448       }
449       B.Enlarge(Tol);
450     }
451   }
452 }