0027261: Incorrect bounding boxes computed for the b-spline faces
[occt.git] / src / BndLib / BndLib_AddSurface.cxx
CommitLineData
b311480e 1// Created on: 1995-07-24
2// Created by: Modelistation
3// Copyright (c) 1995-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 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
973c2be1 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
17// Modified by skv - Fri Aug 27 12:29:04 2004 OCC6503
18
7fd59977 19#include <Adaptor3d_HSurface.hxx>
42cf5bc1 20#include <Adaptor3d_Surface.hxx>
21#include <Bnd_Box.hxx>
7fd59977 22#include <BndLib.hxx>
42cf5bc1 23#include <BndLib_AddSurface.hxx>
9bf0740b 24#include <BSplCLib.hxx>
7fd59977 25#include <ElSLib.hxx>
7fd59977 26#include <Geom_BezierSurface.hxx>
42cf5bc1 27#include <Geom_BSplineSurface.hxx>
28#include <GeomAbs_SurfaceType.hxx>
29#include <gp_Pln.hxx>
30#include <gp_Pnt.hxx>
7fd59977 31#include <Precision.hxx>
42cf5bc1 32#include <TColgp_Array2OfPnt.hxx>
33#include <TColStd_Array1OfInteger.hxx>
34#include <TColStd_Array1OfReal.hxx>
7fd59977 35
36//=======================================================================
37//function : Add
38//purpose :
39//=======================================================================
7fd59977 40void 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
56static 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
83static 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
106static 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
141static 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();
7fd59977 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}
7fd59977 181
9bf0740b 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).
188void 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{
447c7e54 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
9bf0740b 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
7fd59977 220//=======================================================================
221//function : Add
222//purpose :
223//=======================================================================
7fd59977 224void 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) {
9bf0740b 259
7fd59977 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() &&
9bf0740b 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);
7fd59977 292 else
9bf0740b 293 BndLib::Add(S.Sphere(),UMin,UMax,VMin,VMax,Tol,B);
7fd59977 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 {
9bf0740b 307 Standard_Boolean isUseConvexHullAlgorithm = Standard_True;
7fd59977 308 Standard_Real PTol = Precision::Parametric(Precision::Confusion());
9bf0740b 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;
7fd59977 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++){
9bf0740b 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 }
7fd59977 448 }
449 B.Enlarge(Tol);
450 }
451 }
452}