0029311: Implementation of the Oriented Bounding Boxes (OBB) functionality
[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());
168 const BRepAdaptor_Curve anAC(anE);
169 if (!IsLinear(anAC))
170 {
171 if (!theIsTriangulationUsed)
172 // not linear and triangulation usage disabled
173 return 0;
174
175 break;
176 }
177 }
178
179 if (!anExpE.More())
180 // skip planar face with linear edges as its vertices have already been added
181 continue;
182 }
183
184 // Use triangulation of the face
185 const Handle(Poly_Triangulation) &aTrng = BRep_Tool::Triangulation(aF, aLoc);
186 if (aTrng.IsNull())
187 // no triangulation on the face
188 return 0;
189
190 const Standard_Integer aCNode = aTrng->NbNodes();
191 const TColgp_Array1OfPnt& aNodesArr = aTrng->Nodes();
192 for (Standard_Integer i = 1; i <= aCNode; i++)
193 {
194 if (thePts)
195 {
196 const gp_Pnt aP = aLoc.IsIdentity() ? aNodesArr(i) :
197 aNodesArr(i).Transformed(aLoc);
198 (*thePts)(aRetVal) = aP;
199 }
200
201 if (theArrOfToler)
202 {
203 (*theArrOfToler) (aRetVal) = aTrng->Deflection();
204 }
205
206 ++aRetVal;
207 }
208 }
209
210 // Consider edges without faces
211
212 for(anExpE.Init(theS, TopAbs_EDGE, TopAbs_FACE); anExpE.More(); anExpE.Next())
213 {
214 const TopoDS_Edge &anE = TopoDS::Edge(anExpE.Current());
215 const BRepAdaptor_Curve anAC(anE);
216
217 if (IsLinear(anAC))
218 // skip linear edge as its vertices have already been added
219 continue;
220
221 if (!theIsTriangulationUsed)
222 // not linear and triangulation usage disabled
223 return 0;
224
225 const Handle(Poly_Polygon3D) &aPolygon = BRep_Tool::Polygon3D(anE, aLoc);
226 if (aPolygon.IsNull())
227 return 0;
228
229 const Standard_Integer aCNode = aPolygon->NbNodes();
230 const TColgp_Array1OfPnt& aNodesArr = aPolygon->Nodes();
231 for (Standard_Integer i = 1; i <= aCNode; i++)
232 {
233 if (thePts)
234 {
235 const gp_Pnt aP = aLoc.IsIdentity() ? aNodesArr(i) :
236 aNodesArr(i).Transformed(aLoc);
237 (*thePts)(aRetVal) = aP;
238 }
239
240 if (theArrOfToler)
241 {
242 (*theArrOfToler) (aRetVal) = aPolygon->Deflection();
243 }
244
245 ++aRetVal;
246 }
247 }
248
249 return aRetVal;
250}
251
252//=======================================================================
253// Function : IsWCS
254// purpose : Returns 0 if the theDir does not match any axis of WCS.
255// Otherwise, returns the index of correspond axis.
256//=======================================================================
257static Standard_Integer IsWCS(const gp_Dir& theDir)
258{
259 const Standard_Real aToler = Precision::Angular()*Precision::Angular();
260
261 const Standard_Real aX = theDir.X(),
262 aY = theDir.Y(),
263 aZ = theDir.Z();
264
265 const Standard_Real aVx = aY*aY + aZ*aZ,
266 aVy = aX*aX + aZ*aZ,
267 aVz = aX*aX + aY*aY;
268
269 if(aVz < aToler)
270 return 3; // Z-axis
271
272 if(aVy < aToler)
273 return 2; // Y-axis
274
275 if(aVx < aToler)
276 return 1; // X-axis
277
278 return 0;
279}
280
281//=======================================================================
282// Function : CheckPoints
283// purpose : Collects points for DiTO algorithm for OBB construction on
284// linear/planar shapes and shapes having triangulation
285// (http://www.idt.mdh.se/~tla/publ/FastOBBs.pdf).
286//=======================================================================
287static Standard_Boolean CheckPoints(const TopoDS_Shape& theS,
288 const Standard_Boolean theIsTriangulationUsed,
289 const Standard_Boolean theIsShapeToleranceUsed,
290 Bnd_OBB& theOBB)
291{
292 const Standard_Integer aNbPnts = PointsForOBB(theS, theIsTriangulationUsed);
293
294 if(aNbPnts < 1)
295 return Standard_False;
296
297 TColgp_Array1OfPnt anArrPnts(0, theOBB.IsVoid() ? aNbPnts - 1 : aNbPnts + 7);
298 TColStd_Array1OfReal anArrOfTolerances;
299 if(theIsShapeToleranceUsed)
300 {
301 anArrOfTolerances.Resize(anArrPnts.Lower(), anArrPnts.Upper(), Standard_False);
302 anArrOfTolerances.Init(0.0);
303 }
304
305 TColStd_Array1OfReal *aPtrArrTol = theIsShapeToleranceUsed ? &anArrOfTolerances : 0;
306
307 PointsForOBB(theS, theIsTriangulationUsed, &anArrPnts, aPtrArrTol);
308
309 if(!theOBB.IsVoid())
310 {
311 // All points of old OBB have zero-tolerance
312 theOBB.GetVertex(&anArrPnts(aNbPnts));
313 }
314
315#if 0
316 for(Standard_Integer i = anArrPnts.Lower(); i <= anArrPnts.Upper(); i++)
317 {
318 const gp_Pnt &aP = anArrPnts(i);
319 std::cout << "point p" << i << " " << aP.X() << ", " <<
320 aP.Y() << ", " <<
321 aP.Z() << ", "<< std::endl;
322 }
323#endif
324
325 theOBB.ReBuild(anArrPnts, aPtrArrTol);
326
327 return (!theOBB.IsVoid());
328}
329
330//=======================================================================
331// Function : ComputeProperties
332// purpose : Computes properties of theS.
333//=======================================================================
334static void ComputeProperties(const TopoDS_Shape& theS,
335 GProp_GProps& theGCommon)
336{
337 TopExp_Explorer anExp;
338 for(anExp.Init(theS, TopAbs_SOLID); anExp.More(); anExp.Next())
339 {
340 GProp_GProps aG;
341 BRepGProp::VolumeProperties(anExp.Current(), aG, Standard_True);
342 theGCommon.Add(aG);
343 }
344
345 for(anExp.Init(theS, TopAbs_FACE, TopAbs_SOLID); anExp.More(); anExp.Next())
346 {
347 GProp_GProps aG;
348 BRepGProp::SurfaceProperties(anExp.Current(), aG, Standard_True);
349 theGCommon.Add(aG);
350 }
351
352 for(anExp.Init(theS, TopAbs_EDGE, TopAbs_FACE); anExp.More(); anExp.Next())
353 {
354 GProp_GProps aG;
355 BRepGProp::LinearProperties(anExp.Current(), aG, Standard_True);
356 theGCommon.Add(aG);
357 }
358
359 for(anExp.Init(theS, TopAbs_VERTEX, TopAbs_EDGE); anExp.More(); anExp.Next())
360 {
361 GProp_GProps aG(BRep_Tool::Pnt(TopoDS::Vertex(anExp.Current())));
362 theGCommon.Add(aG);
363 }
364}
365
366//=======================================================================
367// Function : ComputePCA
368// purpose : Creates OBB with axes of inertia.
369//=======================================================================
370static void ComputePCA(const TopoDS_Shape& theS,
371 Bnd_OBB& theOBB,
372 const Standard_Boolean theIsTriangulationUsed,
373 const Standard_Boolean theIsOptimal,
374 const Standard_Boolean theIsShapeToleranceUsed)
375{
376 // Compute the transformation matrix to obtain more tight bounding box
377 GProp_GProps aGCommon;
378 ComputeProperties(theS, aGCommon);
379
380 // Transform the shape to the local coordinate system
381 gp_Trsf aTrsf;
382
383 const Standard_Integer anIdx1 =
384 IsWCS(aGCommon.PrincipalProperties().FirstAxisOfInertia());
385 const Standard_Integer anIdx2 =
386 IsWCS(aGCommon.PrincipalProperties().SecondAxisOfInertia());
387
388 if((anIdx1 == 0) || (anIdx2 == 0))
389 {
390 // Coordinate system in which the shape will have the optimal bounding box
391 gp_Ax3 aLocCoordSys(aGCommon.CentreOfMass(),
392 aGCommon.PrincipalProperties().ThirdAxisOfInertia(),
393 aGCommon.PrincipalProperties().FirstAxisOfInertia());
394 aTrsf.SetTransformation(aLocCoordSys);
395 }
396
397 const TopoDS_Shape aST = (aTrsf.Form() == gp_Identity) ? theS :
398 theS.Moved(TopLoc_Location(aTrsf));
399
400 // Initial axis-aligned BndBox
401 Bnd_Box aShapeBox;
402 if(theIsOptimal)
403 {
404 BRepBndLib::AddOptimal(aST, aShapeBox, theIsTriangulationUsed, theIsShapeToleranceUsed);
405 }
406 else
407 {
408 BRepBndLib::Add(aST, aShapeBox);
409 }
410
411 gp_Pnt aPMin = aShapeBox.CornerMin();
412 gp_Pnt aPMax = aShapeBox.CornerMax();
413
414 gp_XYZ aXDir(1, 0, 0);
415 gp_XYZ aYDir(0, 1, 0);
416 gp_XYZ aZDir(0, 0, 1);
417
418 // Compute the center of the box
419 gp_XYZ aCenter = (aPMin.XYZ() + aPMax.XYZ()) / 2.;
420
421 // Compute the half diagonal size of the box.
422 // It takes into account the gap.
423 gp_XYZ anOBBHSize = (aPMax.XYZ() - aPMin.XYZ()) / 2.;
424
425 // Apply transformation if necessary
426 if(aTrsf.Form() != gp_Identity)
427 {
428 aTrsf.Invert();
429 aTrsf.Transforms(aCenter);
430
431 // Make transformation
432 const Standard_Real * aMat = &aTrsf.HVectorialPart().Value(1, 1);
433 // Compute axes directions of the box
434 aXDir = gp_XYZ(aMat[0], aMat[3], aMat[6]);
435 aYDir = gp_XYZ(aMat[1], aMat[4], aMat[7]);
436 aZDir = gp_XYZ(aMat[2], aMat[5], aMat[8]);
437 }
438
439 if(theOBB.IsVoid())
440 {
441 // Create the OBB box
442
443 // Set parameters to the OBB
444 theOBB.SetCenter(aCenter);
445
446 theOBB.SetXComponent(aXDir, anOBBHSize.X());
447 theOBB.SetYComponent(aYDir, anOBBHSize.Y());
448 theOBB.SetZComponent(aZDir, anOBBHSize.Z());
449 theOBB.SetAABox(aTrsf.Form() == gp_Identity);
450 }
451 else
452 {
453 // Recreate the OBB box
454
455 TColgp_Array1OfPnt aListOfPnts(0, 15);
456 theOBB.GetVertex(&aListOfPnts(0));
457
458 const Standard_Real aX = anOBBHSize.X();
459 const Standard_Real aY = anOBBHSize.Y();
460 const Standard_Real aZ = anOBBHSize.Z();
461
462 const gp_XYZ aXext = aX*aXDir,
463 aYext = aY*aYDir,
464 aZext = aZ*aZDir;
465
466 Standard_Integer aPntIdx = 8;
467 aListOfPnts(aPntIdx++) = aCenter - aXext - aYext - aZext;
468 aListOfPnts(aPntIdx++) = aCenter + aXext - aYext - aZext;
469 aListOfPnts(aPntIdx++) = aCenter - aXext + aYext - aZext;
470 aListOfPnts(aPntIdx++) = aCenter + aXext + aYext - aZext;
471 aListOfPnts(aPntIdx++) = aCenter - aXext - aYext + aZext;
472 aListOfPnts(aPntIdx++) = aCenter + aXext - aYext + aZext;
473 aListOfPnts(aPntIdx++) = aCenter - aXext + aYext + aZext;
474 aListOfPnts(aPntIdx++) = aCenter + aXext + aYext + aZext;
475
476 theOBB.ReBuild(aListOfPnts);
477 }
478}
479
480//=======================================================================
481// Function : AddOBB
482// purpose :
483//=======================================================================
484void BRepBndLib::AddOBB(const TopoDS_Shape& theS,
485 Bnd_OBB& theOBB,
486 const Standard_Boolean theIsTriangulationUsed,
487 const Standard_Boolean theIsOptimal,
488 const Standard_Boolean theIsShapeToleranceUsed)
489{
490 if(CheckPoints(theS, theIsTriangulationUsed, theIsShapeToleranceUsed, theOBB))
491 return;
492
493 ComputePCA(theS, theOBB, theIsTriangulationUsed, theIsOptimal, theIsShapeToleranceUsed);
494}