0030428: Modeling Algorithms - BRepBndLib::AddOBB() throws exception on empty Compound
[occt.git] / src / BRepBndLib / BRepBndLib_1.cxx
CommitLineData
1a0339b4 1// Copyright (c) 1999-2017 OPEN CASCADE SAS
2//
3// This file is part of Open CASCADE Technology software library.
4//
5// This library is free software; you can redistribute it and/or modify it under
6// the terms of the GNU Lesser General Public License version 2.1 as published
7// by the Free Software Foundation, with special exception defined in the file
8// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
9// distribution for complete text of the license and disclaimer of any warranty.
10//
11// Alternatively, this file may be used under the terms of Open CASCADE
12// commercial license or contractual agreement.
13
14#include <Adaptor3d_HCurve.hxx>
15#include <Adaptor3d_HSurface.hxx>
16#include <GeomAdaptor_HCurve.hxx>
17#include <BRepBndLib.hxx>
18#include <GProp_GProps.hxx>
19#include <TopoDS_Shape.hxx>
20#include <BRep_Tool.hxx>
21#include <TopoDS.hxx>
22#include <Bnd_OBB.hxx>
23#include <BRepGProp.hxx>
24#include <TopExp_Explorer.hxx>
25#include <GProp_PrincipalProps.hxx>
26#include <gp_Ax3.hxx>
27#include <BRepBuilderAPI_Transform.hxx>
28#include <Bnd_Box.hxx>
29#include <NCollection_List.hxx>
30#include <TColgp_Array1OfPnt.hxx>
31#include <TColStd_Array1OfReal.hxx>
32#include <Geom_Plane.hxx>
33#include <Geom_Line.hxx>
34#include <TColStd_Array1OfInteger.hxx>
35#include <BRepAdaptor_Curve.hxx>
36#include <BRepAdaptor_HSurface.hxx>
37
38#include <Geom_OffsetCurve.hxx>
39#include <Geom_BSplineCurve.hxx>
40#include <Geom_BezierCurve.hxx>
41#include <Geom_BSplineSurface.hxx>
42#include <Geom_BezierSurface.hxx>
43
44//=======================================================================
45// Function : IsLinear
46// purpose : Returns TRUE if theC is line-like.
47//=======================================================================
48static Standard_Boolean IsLinear(const Adaptor3d_Curve& theC)
49{
50 const GeomAbs_CurveType aCT = theC.GetType();
51 if(aCT == GeomAbs_OffsetCurve)
52 {
53 return IsLinear(GeomAdaptor_Curve(theC.OffsetCurve()->BasisCurve()));
54 }
55
56 if((aCT == GeomAbs_BSplineCurve) || (aCT == GeomAbs_BezierCurve))
57 {
58 // Indeed, curves with C0-continuity and degree==1, may be
59 // represented with set of points. It will be possible made
60 // in the future.
61
62 return ((theC.Degree() == 1) &&
63 (theC.Continuity() != GeomAbs_C0));
64 }
65
66 if(aCT == GeomAbs_Line)
67 {
68 return Standard_True;
69 }
70
71 return Standard_False;
72}
73
74//=======================================================================
75// Function : IsPlanar
76// purpose : Returns TRUE if theS is plane-like.
77//=======================================================================
78static Standard_Boolean IsPlanar(const Adaptor3d_Surface& theS)
79{
80 const GeomAbs_SurfaceType aST = theS.GetType();
81 if(aST == GeomAbs_OffsetSurface)
82 {
83 return IsPlanar(theS.BasisSurface()->Surface());
84 }
85
86 if(aST == GeomAbs_SurfaceOfExtrusion)
87 {
88 return IsLinear(theS.BasisCurve()->Curve());
89 }
90
91 if((aST == GeomAbs_BSplineSurface) || (aST == GeomAbs_BezierSurface))
92 {
93 if((theS.UDegree() != 1) || (theS.VDegree() != 1))
94 return Standard_False;
95
96 // Indeed, surfaces with C0-continuity and degree==1, may be
97 // represented with set of points. It will be possible made
98 // in the future.
99
100 return ((theS.UContinuity() != GeomAbs_C0) && (theS.VContinuity() != GeomAbs_C0));
101 }
102
103 if(aST == GeomAbs_Plane)
104 {
105 return Standard_True;
106 }
107
108 return Standard_False;
109}
110
111//=======================================================================
112// Function : PointsForOBB
113// purpose : Returns number of points for array.
114//
115// Attention!!!
116// 1. Start index for thePts must be 0 strictly.
117// 2. Currently, infinite edges/faces (e.g. half-space) are not
118// processed correctly because computation of UV-bounds is a costly operation.
119//=======================================================================
120static Standard_Integer PointsForOBB(const TopoDS_Shape& theS,
121 const Standard_Boolean theIsTriangulationUsed,
122 TColgp_Array1OfPnt* thePts = 0,
123 TColStd_Array1OfReal* theArrOfToler = 0)
124{
125 Standard_Integer aRetVal = 0;
126 TopExp_Explorer anExpF, anExpE;
127
128 // get all vertices from the shape
129 for(anExpF.Init(theS, TopAbs_VERTEX); anExpF.More(); anExpF.Next())
130 {
131 const TopoDS_Vertex &aVert = TopoDS::Vertex(anExpF.Current());
132 if(thePts)
133 {
134 const gp_Pnt aP = BRep_Tool::Pnt(aVert);
135 (*thePts)(aRetVal) = aP;
136 }
137
138 if(theArrOfToler)
139 {
140 (*theArrOfToler) (aRetVal) = BRep_Tool::Tolerance(aVert);
141 }
142
143 ++aRetVal;
144 }
145
146 if(aRetVal == 0)
147 return 0;
148
149 // analyze the faces of the shape on planarity and existence of triangulation
150 TopLoc_Location aLoc;
151 for(anExpF.Init(theS, TopAbs_FACE); anExpF.More(); anExpF.Next())
152 {
153 const TopoDS_Face &aF = TopoDS::Face(anExpF.Current());
154 const BRepAdaptor_Surface anAS(aF, Standard_False);
155
156 if (!IsPlanar(anAS.Surface()))
157 {
158 if (!theIsTriangulationUsed)
159 // not planar and triangulation usage disabled
160 return 0;
161 }
162 else
163 {
164 // planar face
165 for(anExpE.Init(aF, TopAbs_EDGE); anExpE.More(); anExpE.Next())
166 {
167 const TopoDS_Edge &anE = TopoDS::Edge(anExpE.Current());
87a64d53 168 if (BRep_Tool::IsGeometric (anE))
1a0339b4 169 {
87a64d53 170 const BRepAdaptor_Curve anAC(anE);
171 if (!IsLinear(anAC))
172 {
173 if (!theIsTriangulationUsed)
174 // not linear and triangulation usage disabled
175 return 0;
176
177 break;
178 }
1a0339b4 179 }
180 }
181
182 if (!anExpE.More())
183 // skip planar face with linear edges as its vertices have already been added
184 continue;
185 }
186
187 // Use triangulation of the face
188 const Handle(Poly_Triangulation) &aTrng = BRep_Tool::Triangulation(aF, aLoc);
189 if (aTrng.IsNull())
190 // no triangulation on the face
191 return 0;
192
193 const Standard_Integer aCNode = aTrng->NbNodes();
194 const TColgp_Array1OfPnt& aNodesArr = aTrng->Nodes();
195 for (Standard_Integer i = 1; i <= aCNode; i++)
196 {
197 if (thePts)
198 {
199 const gp_Pnt aP = aLoc.IsIdentity() ? aNodesArr(i) :
200 aNodesArr(i).Transformed(aLoc);
201 (*thePts)(aRetVal) = aP;
202 }
203
204 if (theArrOfToler)
205 {
206 (*theArrOfToler) (aRetVal) = aTrng->Deflection();
207 }
208
209 ++aRetVal;
210 }
211 }
212
213 // Consider edges without faces
214
215 for(anExpE.Init(theS, TopAbs_EDGE, TopAbs_FACE); anExpE.More(); anExpE.Next())
216 {
217 const TopoDS_Edge &anE = TopoDS::Edge(anExpE.Current());
87a64d53 218 if (BRep_Tool::IsGeometric (anE))
219 {
220 const BRepAdaptor_Curve anAC(anE);
221 if (IsLinear(anAC))
222 {
223 // skip linear edge as its vertices have already been added
224 continue;
225 }
226 }
1a0339b4 227
228 if (!theIsTriangulationUsed)
229 // not linear and triangulation usage disabled
230 return 0;
231
232 const Handle(Poly_Polygon3D) &aPolygon = BRep_Tool::Polygon3D(anE, aLoc);
233 if (aPolygon.IsNull())
234 return 0;
235
236 const Standard_Integer aCNode = aPolygon->NbNodes();
237 const TColgp_Array1OfPnt& aNodesArr = aPolygon->Nodes();
238 for (Standard_Integer i = 1; i <= aCNode; i++)
239 {
240 if (thePts)
241 {
242 const gp_Pnt aP = aLoc.IsIdentity() ? aNodesArr(i) :
243 aNodesArr(i).Transformed(aLoc);
244 (*thePts)(aRetVal) = aP;
245 }
246
247 if (theArrOfToler)
248 {
249 (*theArrOfToler) (aRetVal) = aPolygon->Deflection();
250 }
251
252 ++aRetVal;
253 }
254 }
255
256 return aRetVal;
257}
258
259//=======================================================================
260// Function : IsWCS
261// purpose : Returns 0 if the theDir does not match any axis of WCS.
262// Otherwise, returns the index of correspond axis.
263//=======================================================================
264static Standard_Integer IsWCS(const gp_Dir& theDir)
265{
266 const Standard_Real aToler = Precision::Angular()*Precision::Angular();
267
268 const Standard_Real aX = theDir.X(),
269 aY = theDir.Y(),
270 aZ = theDir.Z();
271
272 const Standard_Real aVx = aY*aY + aZ*aZ,
273 aVy = aX*aX + aZ*aZ,
274 aVz = aX*aX + aY*aY;
275
276 if(aVz < aToler)
277 return 3; // Z-axis
278
279 if(aVy < aToler)
280 return 2; // Y-axis
281
282 if(aVx < aToler)
283 return 1; // X-axis
284
285 return 0;
286}
287
288//=======================================================================
289// Function : CheckPoints
290// purpose : Collects points for DiTO algorithm for OBB construction on
291// linear/planar shapes and shapes having triangulation
292// (http://www.idt.mdh.se/~tla/publ/FastOBBs.pdf).
293//=======================================================================
294static Standard_Boolean CheckPoints(const TopoDS_Shape& theS,
295 const Standard_Boolean theIsTriangulationUsed,
296 const Standard_Boolean theIsShapeToleranceUsed,
297 Bnd_OBB& theOBB)
298{
299 const Standard_Integer aNbPnts = PointsForOBB(theS, theIsTriangulationUsed);
300
301 if(aNbPnts < 1)
302 return Standard_False;
303
304 TColgp_Array1OfPnt anArrPnts(0, theOBB.IsVoid() ? aNbPnts - 1 : aNbPnts + 7);
305 TColStd_Array1OfReal anArrOfTolerances;
306 if(theIsShapeToleranceUsed)
307 {
308 anArrOfTolerances.Resize(anArrPnts.Lower(), anArrPnts.Upper(), Standard_False);
309 anArrOfTolerances.Init(0.0);
310 }
311
312 TColStd_Array1OfReal *aPtrArrTol = theIsShapeToleranceUsed ? &anArrOfTolerances : 0;
313
314 PointsForOBB(theS, theIsTriangulationUsed, &anArrPnts, aPtrArrTol);
315
316 if(!theOBB.IsVoid())
317 {
318 // All points of old OBB have zero-tolerance
319 theOBB.GetVertex(&anArrPnts(aNbPnts));
320 }
321
322#if 0
323 for(Standard_Integer i = anArrPnts.Lower(); i <= anArrPnts.Upper(); i++)
324 {
325 const gp_Pnt &aP = anArrPnts(i);
326 std::cout << "point p" << i << " " << aP.X() << ", " <<
327 aP.Y() << ", " <<
328 aP.Z() << ", "<< std::endl;
329 }
330#endif
331
332 theOBB.ReBuild(anArrPnts, aPtrArrTol);
333
334 return (!theOBB.IsVoid());
335}
336
337//=======================================================================
338// Function : ComputeProperties
339// purpose : Computes properties of theS.
340//=======================================================================
341static void ComputeProperties(const TopoDS_Shape& theS,
342 GProp_GProps& theGCommon)
343{
344 TopExp_Explorer anExp;
345 for(anExp.Init(theS, TopAbs_SOLID); anExp.More(); anExp.Next())
346 {
347 GProp_GProps aG;
348 BRepGProp::VolumeProperties(anExp.Current(), aG, Standard_True);
349 theGCommon.Add(aG);
350 }
351
352 for(anExp.Init(theS, TopAbs_FACE, TopAbs_SOLID); anExp.More(); anExp.Next())
353 {
354 GProp_GProps aG;
355 BRepGProp::SurfaceProperties(anExp.Current(), aG, Standard_True);
356 theGCommon.Add(aG);
357 }
358
359 for(anExp.Init(theS, TopAbs_EDGE, TopAbs_FACE); anExp.More(); anExp.Next())
360 {
361 GProp_GProps aG;
362 BRepGProp::LinearProperties(anExp.Current(), aG, Standard_True);
363 theGCommon.Add(aG);
364 }
365
366 for(anExp.Init(theS, TopAbs_VERTEX, TopAbs_EDGE); anExp.More(); anExp.Next())
367 {
368 GProp_GProps aG(BRep_Tool::Pnt(TopoDS::Vertex(anExp.Current())));
369 theGCommon.Add(aG);
370 }
371}
372
373//=======================================================================
374// Function : ComputePCA
375// purpose : Creates OBB with axes of inertia.
376//=======================================================================
377static void ComputePCA(const TopoDS_Shape& theS,
378 Bnd_OBB& theOBB,
379 const Standard_Boolean theIsTriangulationUsed,
380 const Standard_Boolean theIsOptimal,
381 const Standard_Boolean theIsShapeToleranceUsed)
382{
383 // Compute the transformation matrix to obtain more tight bounding box
384 GProp_GProps aGCommon;
385 ComputeProperties(theS, aGCommon);
386
387 // Transform the shape to the local coordinate system
388 gp_Trsf aTrsf;
389
390 const Standard_Integer anIdx1 =
391 IsWCS(aGCommon.PrincipalProperties().FirstAxisOfInertia());
392 const Standard_Integer anIdx2 =
393 IsWCS(aGCommon.PrincipalProperties().SecondAxisOfInertia());
394
395 if((anIdx1 == 0) || (anIdx2 == 0))
396 {
397 // Coordinate system in which the shape will have the optimal bounding box
398 gp_Ax3 aLocCoordSys(aGCommon.CentreOfMass(),
399 aGCommon.PrincipalProperties().ThirdAxisOfInertia(),
400 aGCommon.PrincipalProperties().FirstAxisOfInertia());
401 aTrsf.SetTransformation(aLocCoordSys);
402 }
403
404 const TopoDS_Shape aST = (aTrsf.Form() == gp_Identity) ? theS :
405 theS.Moved(TopLoc_Location(aTrsf));
406
407 // Initial axis-aligned BndBox
408 Bnd_Box aShapeBox;
409 if(theIsOptimal)
410 {
411 BRepBndLib::AddOptimal(aST, aShapeBox, theIsTriangulationUsed, theIsShapeToleranceUsed);
412 }
413 else
414 {
415 BRepBndLib::Add(aST, aShapeBox);
416 }
0939d4cf 417 if (aShapeBox.IsVoid())
418 {
419 return;
420 }
1a0339b4 421
422 gp_Pnt aPMin = aShapeBox.CornerMin();
423 gp_Pnt aPMax = aShapeBox.CornerMax();
424
425 gp_XYZ aXDir(1, 0, 0);
426 gp_XYZ aYDir(0, 1, 0);
427 gp_XYZ aZDir(0, 0, 1);
428
429 // Compute the center of the box
430 gp_XYZ aCenter = (aPMin.XYZ() + aPMax.XYZ()) / 2.;
431
432 // Compute the half diagonal size of the box.
433 // It takes into account the gap.
434 gp_XYZ anOBBHSize = (aPMax.XYZ() - aPMin.XYZ()) / 2.;
435
436 // Apply transformation if necessary
437 if(aTrsf.Form() != gp_Identity)
438 {
439 aTrsf.Invert();
440 aTrsf.Transforms(aCenter);
441
442 // Make transformation
443 const Standard_Real * aMat = &aTrsf.HVectorialPart().Value(1, 1);
444 // Compute axes directions of the box
445 aXDir = gp_XYZ(aMat[0], aMat[3], aMat[6]);
446 aYDir = gp_XYZ(aMat[1], aMat[4], aMat[7]);
447 aZDir = gp_XYZ(aMat[2], aMat[5], aMat[8]);
448 }
449
450 if(theOBB.IsVoid())
451 {
452 // Create the OBB box
453
454 // Set parameters to the OBB
455 theOBB.SetCenter(aCenter);
456
457 theOBB.SetXComponent(aXDir, anOBBHSize.X());
458 theOBB.SetYComponent(aYDir, anOBBHSize.Y());
459 theOBB.SetZComponent(aZDir, anOBBHSize.Z());
460 theOBB.SetAABox(aTrsf.Form() == gp_Identity);
461 }
462 else
463 {
464 // Recreate the OBB box
465
466 TColgp_Array1OfPnt aListOfPnts(0, 15);
467 theOBB.GetVertex(&aListOfPnts(0));
468
469 const Standard_Real aX = anOBBHSize.X();
470 const Standard_Real aY = anOBBHSize.Y();
471 const Standard_Real aZ = anOBBHSize.Z();
472
473 const gp_XYZ aXext = aX*aXDir,
474 aYext = aY*aYDir,
475 aZext = aZ*aZDir;
476
477 Standard_Integer aPntIdx = 8;
478 aListOfPnts(aPntIdx++) = aCenter - aXext - aYext - aZext;
479 aListOfPnts(aPntIdx++) = aCenter + aXext - aYext - aZext;
480 aListOfPnts(aPntIdx++) = aCenter - aXext + aYext - aZext;
481 aListOfPnts(aPntIdx++) = aCenter + aXext + aYext - aZext;
482 aListOfPnts(aPntIdx++) = aCenter - aXext - aYext + aZext;
483 aListOfPnts(aPntIdx++) = aCenter + aXext - aYext + aZext;
484 aListOfPnts(aPntIdx++) = aCenter - aXext + aYext + aZext;
485 aListOfPnts(aPntIdx++) = aCenter + aXext + aYext + aZext;
486
487 theOBB.ReBuild(aListOfPnts);
488 }
489}
490
491//=======================================================================
492// Function : AddOBB
493// purpose :
494//=======================================================================
495void BRepBndLib::AddOBB(const TopoDS_Shape& theS,
496 Bnd_OBB& theOBB,
497 const Standard_Boolean theIsTriangulationUsed,
498 const Standard_Boolean theIsOptimal,
499 const Standard_Boolean theIsShapeToleranceUsed)
500{
501 if(CheckPoints(theS, theIsTriangulationUsed, theIsShapeToleranceUsed, theOBB))
502 return;
503
504 ComputePCA(theS, theOBB, theIsTriangulationUsed, theIsOptimal, theIsShapeToleranceUsed);
505}