1 // Created by: Eugeny MALTCHIKOV
2 // Copyright (c) 2017 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
15 #include <Bnd_OBB.hxx>
17 #include <Bnd_Tools.hxx>
18 #include <Bnd_Range.hxx>
20 #include <BVH_BoxSet.hxx>
21 #include <BVH_LinearBuilder.hxx>
23 #include <BVH_Traverse.hxx>
25 #include <NCollection_Array1.hxx>
26 #include <Precision.hxx>
27 #include <Standard_Dump.hxx>
28 #include <TColStd_Array1OfReal.hxx>
30 //! Auxiliary class to select from the points stored in
31 //! BVH tree the two points giving the extreme projection
32 //! parameters on the axis
33 class OBB_ExtremePointsSelector :
34 public BVH_Traverse <Standard_Real, 3, BVH_BoxSet <Standard_Real, 3, gp_XYZ>, Bnd_Range>
39 OBB_ExtremePointsSelector() :
40 BVH_Traverse <Standard_Real, 3, BVH_BoxSet <Standard_Real, 3, gp_XYZ>, Bnd_Range>(),
41 myPrmMin (RealLast()),
42 myPrmMax (RealFirst())
45 public: //! @name Set axis for projection
48 void SetAxis (const gp_XYZ& theAxis) { myAxis = theAxis; }
50 public: //! @name Clears the points from previous runs
55 myPrmMin = RealLast();
56 myPrmMax = RealFirst();
59 public: //! @name Getting the results
61 //! Returns the minimal projection parameter
62 Standard_Real MinPrm() const { return myPrmMin; }
64 //! Returns the maximal projection parameter
65 Standard_Real MaxPrm() const { return myPrmMax; }
67 //! Returns the minimal projection point
68 const gp_XYZ& MinPnt() const { return myPntMin; }
70 //! Returns the maximal projection point
71 const gp_XYZ& MaxPnt() const { return myPntMax; }
73 public: //! @name Definition of rejection/acceptance rules
75 //! Defines the rules for node rejection
76 virtual Standard_Boolean RejectNode (const BVH_Vec3d& theCMin,
77 const BVH_Vec3d& theCMax,
78 Bnd_Range& theMetric) const Standard_OVERRIDE
80 if (myPrmMin > myPrmMax)
81 // No parameters computed yet
82 return Standard_False;
84 Standard_Real aPrmMin = myPrmMin, aPrmMax = myPrmMax;
85 Standard_Boolean isToReject = Standard_True;
87 // Check if the current node is between already found parameters
88 for (Standard_Integer i = 0; i < 2; ++i)
90 Standard_Real x = !i ? theCMin.x() : theCMax.x();
91 for (Standard_Integer j = 0; j < 2; ++j)
93 Standard_Real y = !j ? theCMin.y() : theCMax.y();
94 for (Standard_Integer k = 0; k < 2; ++k)
96 Standard_Real z = !k ? theCMin.z() : theCMax.z();
98 Standard_Real aPrm = myAxis.Dot (gp_XYZ (x, y, z));
102 isToReject = Standard_False;
104 else if (aPrm > aPrmMax)
107 isToReject = Standard_False;
113 theMetric = Bnd_Range (aPrmMin, aPrmMax);
118 //! Rules for node rejection by the metric
119 virtual Standard_Boolean RejectMetric (const Bnd_Range& theMetric) const Standard_OVERRIDE
121 if (myPrmMin > myPrmMax)
122 // no parameters computed
123 return Standard_False;
125 Standard_Real aMin, aMax;
126 if (!theMetric.GetBounds (aMin, aMax))
128 return Standard_False;
130 // Check if the box of the branch is inside of the already computed parameters
131 return aMin > myPrmMin && aMax < myPrmMax;
134 //! Defines the rules for leaf acceptance
135 virtual Standard_Boolean Accept (const Standard_Integer theIndex,
136 const Bnd_Range&) Standard_OVERRIDE
138 const gp_XYZ& theLeaf = myBVHSet->Element (theIndex);
139 Standard_Real aPrm = myAxis.Dot (theLeaf);
150 return Standard_True;
153 public: //! @name Choosing the best branch
155 //! Returns true if the metric of the left branch is better than the metric of the right
156 virtual Standard_Boolean IsMetricBetter (const Bnd_Range& theLeft,
157 const Bnd_Range& theRight) const Standard_OVERRIDE
159 if (myPrmMin > myPrmMax)
160 // no parameters computed
161 return Standard_True;
163 Standard_Real aMin[2], aMax[2];
164 if (!theLeft.GetBounds (aMin[0], aMax[0]) ||
165 !theRight.GetBounds (aMin[1], aMax[1]))
167 return Standard_True;
169 // Choose branch with larger extension over computed parameters
170 Standard_Real anExt[2] = {0.0, 0.0};
171 for (int i = 0; i < 2; ++i)
173 if (aMin[i] < myPrmMin) anExt[i] += myPrmMin - aMin[i];
174 if (aMax[i] > myPrmMax) anExt[i] += aMax[i] - myPrmMax;
176 return anExt[0] > anExt[1];
179 protected: //! @name Fields
181 gp_XYZ myAxis; //!< Axis to project the points to
182 Standard_Real myPrmMin; //!< Minimal projection parameter
183 Standard_Real myPrmMax; //!< Maximal projection parameter
184 gp_XYZ myPntMin; //!< Minimal projection point
185 gp_XYZ myPntMax; //!< Maximal projection point
188 //! Tool for OBB construction
192 //! Constructor. theL - list of points.
193 //! theLT is a pointer to the list of tolerances
194 //! (i-th element of this array is a tolerance
195 //! of i-th point in theL). If theLT is empty
196 //! then the tolerance of every point is equal to 0.
197 //! Attention! The objects, which theL and theLT links on,
198 //! must be available during all time of OBB creation
199 //! (i.e. while the object of OBBTool exists).
200 OBBTool(const TColgp_Array1OfPnt& theL,
201 const TColStd_Array1OfReal *theLT = 0,
202 Standard_Boolean theIsOptimal = Standard_False);
204 //! DiTO algorithm for OBB construction
205 //! (http://www.idt.mdh.se/~tla/publ/FastOBBs.pdf)
206 void ProcessDiTetrahedron();
208 //! Creates OBB with already computed parameters
209 void BuildBox(Bnd_OBB& theBox);
213 // Computes the extreme points on the set of Initial axes
214 void ComputeExtremePoints ();
216 //! Works with the triangle set by the points in myTriIdx.
217 //! If theIsBuiltTrg == TRUE, new set of triangles will be
219 void ProcessTriangle(const Standard_Integer theIdx1,
220 const Standard_Integer theIdx2,
221 const Standard_Integer theIdx3,
222 const Standard_Boolean theIsBuiltTrg);
224 //! Computes myTriIdx[2]
225 void FillToTriangle3();
227 //! Computes myTriIdx[3] and myTriIdx[4]
228 void FillToTriangle5(const gp_XYZ& theNormal,
229 const gp_XYZ& theBarryCenter);
231 //! Returns half of the Surface area of the box
232 static Standard_Real ComputeQuality(const Standard_Real* const thePrmArr)
234 const Standard_Real aDX = thePrmArr[1] - thePrmArr[0],
235 aDY = thePrmArr[3] - thePrmArr[2],
236 aDZ = thePrmArr[5] - thePrmArr[4];
238 return (aDX*aDY + aDY*aDZ + aDX*aDZ);
242 //! Assignment operator is forbidden
243 OBBTool& operator=(const OBBTool&);
246 //! Params structure stores the two values meaning
247 //! min and max parameters on the axis
251 _ParamMin(RealLast()), _ParamMax(RealFirst())
254 Params(Standard_Real theMin, Standard_Real theMax)
255 : _ParamMin(theMin), _ParamMax(theMax)
258 Standard_Real _ParamMin;
259 Standard_Real _ParamMax;
262 //! Computes the Minimal and maximal parameters on the vector
263 //! connecting the points myLExtremalPoints[theId1] and myLExtremalPoints[theId2]
264 void ComputeParams (const Standard_Integer theId1,
265 const Standard_Integer theId2,
266 Standard_Real &theMin,
267 Standard_Real &theMax)
269 theMin = myParams[theId1][theId2]._ParamMin;
270 theMax = myParams[theId1][theId2]._ParamMax;
274 FindMinMax ((myLExtremalPoints[theId1] - myLExtremalPoints[theId2]).Normalized(), theMin, theMax);
275 myParams[theId1][theId2]._ParamMin = myParams[theId2][theId1]._ParamMin = theMin;
276 myParams[theId1][theId2]._ParamMax = myParams[theId2][theId1]._ParamMax = theMax;
280 //! Looks for the min-max parameters on the axis.
281 //! For optimal case projects all the points on the axis,
282 //! for not optimal - only the set of extreme points.
283 void FindMinMax (const gp_XYZ& theAxis,
284 Standard_Real &theMin,
285 Standard_Real &theMax)
287 theMin = RealLast(), theMax = RealFirst();
290 Project (theAxis, theMin, theMax);
293 for (Standard_Integer i = 0; i < myNbExtremalPoints; ++i)
295 Standard_Real aPrm = theAxis.Dot (myLExtremalPoints[i]);
296 if (aPrm < theMin) theMin = aPrm;
297 if (aPrm > theMax) theMax = aPrm;
302 //! Projects the set of points on the axis
303 void Project (const gp_XYZ& theAxis,
304 Standard_Real& theMin, Standard_Real& theMax,
305 gp_XYZ* thePntMin = 0, gp_XYZ* thePntMax = 0)
307 theMin = RealLast(), theMax = RealFirst();
312 OBB_ExtremePointsSelector anExtremePointsSelector;
313 anExtremePointsSelector.SetBVHSet (myPointBoxSet.get());
314 anExtremePointsSelector.SetAxis (theAxis);
315 anExtremePointsSelector.Select();
316 theMin = anExtremePointsSelector.MinPrm();
317 theMax = anExtremePointsSelector.MaxPrm();
318 if (thePntMin) *thePntMin = anExtremePointsSelector.MinPnt();
319 if (thePntMax) *thePntMax = anExtremePointsSelector.MaxPnt();
323 // Project all points
324 for (Standard_Integer iP = myPntsList.Lower(); iP <= myPntsList.Upper(); ++iP)
326 const gp_XYZ& aPoint = myPntsList(iP).XYZ();
327 const Standard_Real aPrm = theAxis.Dot (aPoint);
346 //! Number of the initial axes.
347 static const Standard_Integer myNbInitAxes = 7;
349 //! Number of extremal points
350 static const Standard_Integer myNbExtremalPoints = 2 * myNbInitAxes;
352 //! The source list of points
353 const TColgp_Array1OfPnt& myPntsList;
355 //! Pointer to the array of tolerances
356 const TColStd_Array1OfReal *myListOfTolers;
358 //! Points of ditetrahedron
359 //! given by their indices in myLExtremalPoints.
360 Standard_Integer myTriIdx[5];
362 //! List of extremal points
363 gp_XYZ myLExtremalPoints[myNbExtremalPoints];
365 //! The axes of the box (always normalized or
366 //! can be null-vector)
369 //! The surface area of the OBB
370 Standard_Real myQualityCriterion;
372 //! Defines if the OBB should be computed more tight.
373 //! Takes more time, but the volume is less.
374 Standard_Boolean myOptimal;
376 //! Point box set organized with BVH
377 opencascade::handle<BVH_BoxSet <Standard_Real, 3, gp_XYZ>> myPointBoxSet;
379 //! Stored min/max parameters for the axes between extremal points
380 Params myParams[myNbExtremalPoints][myNbExtremalPoints];
383 //=======================================================================
384 // Function : SetMinMax
386 // ATTENTION!!! thePrmArr must be initialized before this method calling.
387 //=======================================================================
388 static inline void SetMinMax(Standard_Real* const thePrmArr,
389 const Standard_Real theNewParam)
391 if(theNewParam < thePrmArr[0])
393 thePrmArr[0] = theNewParam;
395 if(theNewParam > thePrmArr[1])
397 thePrmArr[1] = theNewParam;
401 //=======================================================================
402 // Function : Constructor
404 //=======================================================================
406 OBBTool(const TColgp_Array1OfPnt& theL,
407 const TColStd_Array1OfReal *theLT,
408 const Standard_Boolean theIsOptimal) : myPntsList(theL),
409 myListOfTolers(theLT),
410 myQualityCriterion(RealLast()),
411 myOptimal (theIsOptimal)
415 // Use linear builder for BVH construction with 30 elements in the leaf
416 opencascade::handle<BVH_LinearBuilder<Standard_Real, 3> > aLBuilder =
417 new BVH_LinearBuilder<Standard_Real, 3> (30);
418 myPointBoxSet = new BVH_BoxSet <Standard_Real, 3, gp_XYZ> (aLBuilder);
419 myPointBoxSet->SetSize(myPntsList.Length());
421 // Add the points into Set
422 for (Standard_Integer iP = 0; iP < theL.Length(); ++iP)
424 const gp_Pnt& aP = theL (iP);
425 Standard_Real aTol = theLT ? theLT->Value(iP) : Precision::Confusion();
426 BVH_Box <Standard_Real, 3> aBox (BVH_Vec3d (aP.X() - aTol, aP.Y() - aTol, aP.Z() - aTol),
427 BVH_Vec3d (aP.X() + aTol, aP.Y() + aTol, aP.Z() + aTol));
428 myPointBoxSet->Add (aP.XYZ(), aBox);
431 myPointBoxSet->Build();
434 ComputeExtremePoints();
437 //=======================================================================
438 // Function : ComputeExtremePoints
440 //=======================================================================
441 void OBBTool::ComputeExtremePoints()
443 // Six initial axes show great quality on the Optimal OBB, plus
444 // the performance is better (due to the less number of operations).
445 // But they show worse quality for the not optimal approach.
446 //const Standard_Real a = (sqrt(5) - 1) / 2.;
447 //const gp_XYZ anInitialAxes6[myNbInitAxes] = { gp_XYZ (0, 1, a),
448 // gp_XYZ (0, 1, -a),
450 // gp_XYZ (1, -a, 0),
452 // gp_XYZ (a, 0, -1) };
453 const Standard_Real aSqrt3 = Sqrt(3);
454 const gp_XYZ anInitialAxes7[myNbInitAxes] = { gp_XYZ (1.0, 0.0, 0.0),
455 gp_XYZ (0.0, 1.0, 0.0),
456 gp_XYZ (0.0, 0.0, 1.0),
457 gp_XYZ (1.0, 1.0, 1.0) / aSqrt3,
458 gp_XYZ (1.0, 1.0, -1.0) / aSqrt3,
459 gp_XYZ (1.0, -1.0, 1.0) / aSqrt3,
460 gp_XYZ (1.0, -1.0, -1.0) / aSqrt3 };
462 // Set of initial axes
463 const gp_XYZ *anInitialAxesArray = anInitialAxes7;
465 // Min and Max parameter
466 Standard_Real aParams[myNbExtremalPoints];
467 // Look for the extremal points (myLExtremalPoints)
468 for (Standard_Integer anAxeInd = 0, aPrmInd = -1; anAxeInd < myNbInitAxes; ++anAxeInd)
470 Standard_Integer aMinInd = ++aPrmInd, aMaxInd = ++aPrmInd;
471 aParams[aMinInd] = RealLast();
472 aParams[aMaxInd] = -RealLast();
473 Project (anInitialAxesArray[anAxeInd],
474 aParams[aMinInd], aParams[aMaxInd],
475 &myLExtremalPoints[aMinInd], &myLExtremalPoints[aMaxInd]);
478 // For not optimal box it is necessary to compute the max axis
479 // created by the maximally distant extreme points
482 for(Standard_Integer i = 0; i < 5; i++)
483 myTriIdx[i] = INT_MAX;
485 // Compute myTriIdx[0] and myTriIdx[1].
486 Standard_Real aMaxSqDist = -1.0;
487 for (Standard_Integer aPrmInd = 0; aPrmInd < myNbExtremalPoints; aPrmInd += 2)
489 const gp_Pnt &aP1 = myLExtremalPoints[aPrmInd],
490 &aP2 = myLExtremalPoints[aPrmInd + 1];
491 const Standard_Real aSqDist = aP1.SquareDistance(aP2);
492 if (aSqDist > aMaxSqDist)
494 aMaxSqDist = aSqDist;
495 myTriIdx[0] = aPrmInd;
496 myTriIdx[1] = aPrmInd + 1;
500 // Compute the maximal axis orthogonal to the found one
505 //=======================================================================
506 // Function : FillToTriangle3
507 // purpose : Two value of myTriIdx array is known. Let us find myTriIdx[2].
508 // It must be in maximal distance from the infinite axis going
509 // through the points with indexes myTriIdx[0] and myTriIdx[1].
510 //=======================================================================
511 void OBBTool::FillToTriangle3()
513 const gp_XYZ &aP0 = myLExtremalPoints[myTriIdx[0]];
514 const gp_XYZ anAxis = myLExtremalPoints[myTriIdx[1]] - aP0;
515 Standard_Real aMaxSqDist = -1.0;
516 for(Standard_Integer i = 0; i < myNbExtremalPoints; i++)
518 if((i == myTriIdx[0]) || (i == myTriIdx[1]))
521 const gp_XYZ &aP = myLExtremalPoints[i];
522 const Standard_Real aDistToAxe = anAxis.CrossSquareMagnitude(aP - aP0);
523 if(aDistToAxe > aMaxSqDist)
526 aMaxSqDist = aDistToAxe;
531 //=======================================================================
532 // Function : FillToTriangle5
533 // purpose : Three value of myTriIdx array is known.
534 // Let us find myTriIdx[3] and myTriIdx[4].
535 // They must be in the different sides of the plane of
536 // triangle set by points myTriIdx[0], myTriIdx[1] and
537 // myTriIdx[2]. Moreover, the distance from these points
538 // to the triangle plane must be maximal.
539 //=======================================================================
540 void OBBTool::FillToTriangle5(const gp_XYZ& theNormal,
541 const gp_XYZ& theBarryCenter)
543 Standard_Real aParams[2] = {0.0, 0.0};
544 Standard_Integer id3 = -1, id4 = -1;
546 for(Standard_Integer aPtIdx = 0; aPtIdx < myNbExtremalPoints; aPtIdx++)
548 if((aPtIdx == myTriIdx[0]) || (aPtIdx == myTriIdx[1]) || (aPtIdx == myTriIdx[2]))
551 const gp_XYZ &aCurrPoint = myLExtremalPoints[aPtIdx];
552 const Standard_Real aParam = theNormal.Dot(aCurrPoint - theBarryCenter);
554 if (aParam < aParams[0])
559 else if (aParam > aParams[1])
566 // The points must be in the different sides of the triangle plane.
567 if (id3 >= 0 && aParams[0] < -Precision::Confusion())
570 if (id4 >= 0 && aParams[1] > Precision::Confusion())
574 //=======================================================================
575 // Function : ProcessTriangle
576 // purpose : Choose the optimal box with triple axes containing normal
577 // to the triangle and some edge of the triangle (3rd axis is
578 // computed from these two ones).
579 //=======================================================================
580 void OBBTool::ProcessTriangle(const Standard_Integer theIdx1,
581 const Standard_Integer theIdx2,
582 const Standard_Integer theIdx3,
583 const Standard_Boolean theIsBuiltTrg)
585 const Standard_Integer aNbAxes = 3;
587 // All axes must be normalized in order to provide correct area computation
588 // (see ComputeQuality(...) method).
589 int ID1[3] = { theIdx2, theIdx3, theIdx1 },
590 ID2[3] = { theIdx1, theIdx2, theIdx3 };
591 gp_XYZ aYAxis[aNbAxes] = {(myLExtremalPoints[ID1[0]] - myLExtremalPoints[ID2[0]]),
592 (myLExtremalPoints[ID1[1]] - myLExtremalPoints[ID2[1]]),
593 (myLExtremalPoints[ID1[2]] - myLExtremalPoints[ID2[2]])};
595 // Normal to the triangle plane
596 gp_XYZ aZAxis = aYAxis[0].Crossed(aYAxis[1]);
598 Standard_Real aSqMod = aZAxis.SquareModulus();
600 if (aSqMod < Precision::SquareConfusion())
603 aZAxis /= Sqrt(aSqMod);
605 gp_XYZ aXAxis[aNbAxes];
606 for (Standard_Integer i = 0; i < aNbAxes; i++)
607 aXAxis[i] = aYAxis[i].Crossed(aZAxis).Normalized();
610 FillToTriangle5 (aZAxis, myLExtremalPoints[theIdx1]);
612 // Min and Max parameter
613 const Standard_Integer aNbPoints = 2 * aNbAxes;
615 // Compute Min/Max params for ZAxis
616 Standard_Real aParams[aNbPoints];
617 FindMinMax (aZAxis, aParams[4], aParams[5]); // Compute params on ZAxis once
619 Standard_Integer aMinIdx = -1;
620 for(Standard_Integer anAxeInd = 0; anAxeInd < aNbAxes; anAxeInd++)
622 const gp_XYZ &aAX = aXAxis[anAxeInd];
623 // Compute params on XAxis
624 FindMinMax (aAX, aParams[0], aParams[1]);
625 // Compute params on YAxis checking for stored values
626 ComputeParams (ID1[anAxeInd], ID2[anAxeInd], aParams[2], aParams[3]);
628 const Standard_Real anArea = ComputeQuality(aParams);
629 if (anArea < myQualityCriterion)
631 myQualityCriterion = anArea;
639 myAxes[0] = aXAxis[aMinIdx];
640 myAxes[1] = aYAxis[aMinIdx].Normalized();
643 //=======================================================================
644 // Function : ProcessDiTetrahedron
645 // purpose : DiTo-algorithm (http://www.idt.mdh.se/~tla/publ/FastOBBs.pdf)
646 //=======================================================================
647 void OBBTool::ProcessDiTetrahedron()
649 // To compute the optimal OBB it is necessary to check all possible
650 // axes created by the extremal points. It is also necessary to project
651 // all the points on the axis, as for each different axis there will be
652 // different extremal points.
655 for (Standard_Integer i = 0; i < myNbExtremalPoints - 2; i++)
657 for (Standard_Integer j = i + 1; j < myNbExtremalPoints - 1; j++)
659 for (Standard_Integer k = j + 1; k < myNbExtremalPoints; k++)
661 ProcessTriangle (i, j, k, Standard_False);
668 // Use the standard DiTo approach
669 ProcessTriangle(myTriIdx[0], myTriIdx[1], myTriIdx[2], Standard_True);
671 if (myTriIdx[3] <= myNbExtremalPoints)
673 ProcessTriangle(myTriIdx[0], myTriIdx[1], myTriIdx[3], Standard_False);
674 ProcessTriangle(myTriIdx[1], myTriIdx[2], myTriIdx[3], Standard_False);
675 ProcessTriangle(myTriIdx[0], myTriIdx[2], myTriIdx[3], Standard_False);
678 if (myTriIdx[4] <= myNbExtremalPoints)
680 ProcessTriangle(myTriIdx[0], myTriIdx[1], myTriIdx[4], Standard_False);
681 ProcessTriangle(myTriIdx[1], myTriIdx[2], myTriIdx[4], Standard_False);
682 ProcessTriangle(myTriIdx[0], myTriIdx[2], myTriIdx[4], Standard_False);
687 //=======================================================================
688 // Function : BuildBox
690 //=======================================================================
691 void OBBTool::BuildBox(Bnd_OBB& theBox)
695 // In fact, use Precision::SquareConfusion().
696 const Standard_Boolean isOBB = myAxes[0].SquareModulus()*
697 myAxes[1].SquareModulus()*
698 myAxes[2].SquareModulus() > 1.0e-14;
700 const gp_Dir aXDir = isOBB ? myAxes[0] : gp_Dir(1, 0, 0);
701 const gp_Dir aYDir = isOBB ? myAxes[1] : gp_Dir(0, 1, 0);
702 const gp_Dir aZDir = isOBB ? myAxes[2] : gp_Dir(0, 0, 1);
704 const Standard_Integer aNbPoints = 6;
705 Standard_Real aParams[aNbPoints];
707 gp_XYZ aFCurrPoint = myPntsList.First().XYZ();
709 aParams[0] = aParams[1] = aFCurrPoint.Dot(aXDir.XYZ());
710 aParams[2] = aParams[3] = aFCurrPoint.Dot(aYDir.XYZ());
711 aParams[4] = aParams[5] = aFCurrPoint.Dot(aZDir.XYZ());
713 if(myListOfTolers != 0)
715 const Standard_Real aTol = myListOfTolers->First();
724 for(Standard_Integer i = myPntsList.Lower() + 1; i <= myPntsList.Upper(); i++)
726 const gp_XYZ &aCurrPoint = myPntsList(i).XYZ();
727 const Standard_Real aDx = aCurrPoint.Dot(aXDir.XYZ()),
728 aDy = aCurrPoint.Dot(aYDir.XYZ()),
729 aDz = aCurrPoint.Dot(aZDir.XYZ());
731 if(myListOfTolers == 0)
733 SetMinMax(&aParams[0], aDx);
734 SetMinMax(&aParams[2], aDy);
735 SetMinMax(&aParams[4], aDz);
739 const Standard_Real aTol = myListOfTolers->Value(i);
740 aParams[0] = Min(aParams[0], aDx - aTol);
741 aParams[1] = Max(aParams[1], aDx + aTol);
742 aParams[2] = Min(aParams[2], aDy - aTol);
743 aParams[3] = Max(aParams[3], aDy + aTol);
744 aParams[4] = Min(aParams[4], aDz - aTol);
745 aParams[5] = Max(aParams[5], aDz + aTol);
750 const Standard_Real aHX = 0.5*(aParams[1] - aParams[0]);
751 const Standard_Real aHY = 0.5*(aParams[3] - aParams[2]);
752 const Standard_Real aHZ = 0.5*(aParams[5] - aParams[4]);
754 const gp_XYZ aCenter = 0.5*((aParams[1] + aParams[0])*aXDir.XYZ() +
755 (aParams[3] + aParams[2])*aYDir.XYZ() +
756 (aParams[5] + aParams[4])*aZDir.XYZ());
758 theBox.SetCenter(aCenter);
759 theBox.SetXComponent(aXDir, aHX);
760 theBox.SetYComponent(aYDir, aHY);
761 theBox.SetZComponent(aZDir, aHZ);
762 theBox.SetAABox(!isOBB);
765 // =======================================================================
766 // function : ReBuild
767 // purpose : http://www.idt.mdh.se/~tla/publ/
768 // =======================================================================
769 void Bnd_OBB::ReBuild(const TColgp_Array1OfPnt& theListOfPoints,
770 const TColStd_Array1OfReal *theListOfTolerances,
771 const Standard_Boolean theIsOptimal)
773 switch(theListOfPoints.Length())
776 ProcessOnePoint(theListOfPoints.First());
777 if(theListOfTolerances)
778 Enlarge(theListOfTolerances->First());
782 const Standard_Real aTol1 = (theListOfTolerances == 0) ? 0.0 :
783 theListOfTolerances->First();
785 const Standard_Real aTol2 = (theListOfTolerances == 0) ? 0.0 :
786 theListOfTolerances->Last();
788 const gp_XYZ &aP1 = theListOfPoints.First().XYZ(),
789 &aP2 = theListOfPoints.Last().XYZ();
790 const gp_XYZ aDP = aP2 - aP1;
791 const Standard_Real aDPm = aDP.Modulus();
792 myIsAABox = Standard_False;
793 myHDims[1] = myHDims[2] = Max(aTol1, aTol2);
795 if(aDPm < Precision::Confusion())
797 ProcessOnePoint(aP1);
798 Enlarge(myHDims[1] + Precision::Confusion());
802 myHDims[0] = 0.5*(aDPm+aTol1+aTol2);
803 myAxes[0] = aDP/aDPm;
804 if(Abs(myAxes[0].X()) > Abs(myAxes[0].Y()))
806 // Z-coord. is maximal or X-coord. is maximal
807 myAxes[1].SetCoord(-myAxes[0].Z(), 0.0, myAxes[0].X());
811 // Z-coord. is maximal or Y-coord. is maximal
812 myAxes[1].SetCoord(0.0, -myAxes[0].Z(), myAxes[0].Y());
815 myAxes[2] = myAxes[0].Crossed(myAxes[1]).Normalized();
816 myCenter = aP1 + 0.5*(aDPm - aTol1 + aTol2)*myAxes[0];
823 OBBTool aTool(theListOfPoints, theListOfTolerances, theIsOptimal);
824 aTool.ProcessDiTetrahedron();
825 aTool.BuildBox(*this);
828 // =======================================================================
831 // =======================================================================
832 Standard_Boolean Bnd_OBB::IsOut(const Bnd_OBB& theOther) const
834 if (IsVoid() || theOther.IsVoid())
835 return Standard_True;
837 if (myIsAABox && theOther.myIsAABox)
839 return ((Abs(theOther.myCenter.X() - myCenter.X()) > theOther.myHDims[0] + myHDims[0]) ||
840 (Abs(theOther.myCenter.Y() - myCenter.Y()) > theOther.myHDims[1] + myHDims[1]) ||
841 (Abs(theOther.myCenter.Z() - myCenter.Z()) > theOther.myHDims[2] + myHDims[2]));
844 // According to the Separating Axis Theorem for Oriented Bounding Boxes
845 // it is necessary to check the 15 separating axes (Ls):
846 // - 6 axes of the boxes;
847 // - 9 cross products of the axes of the boxes.
848 // If any of these axes is valid, the boxes do not interfere.
850 // The algorithm is following:
851 // 1. Compute the "length" for j-th BndBox (j=1...2) according to the formula:
852 // L(j)=Sum(myHDims[i]*Abs(myAxes[i].Dot(Ls)))
853 // 2. If (theCenter2 - theCenter1).Dot(Ls) > (L(1) + L(2))
854 // then the considered OBBs are not interfered in terms of the axis Ls.
856 // If OBBs are not interfered in terms of at least one axis (of 15) then
857 // they are not interfered at all.
859 // Precomputed difference between centers
860 gp_XYZ D = theOther.myCenter - myCenter;
862 // Check the axes of the this box, i.e. L is one of myAxes
863 // Since the Dot product of two of these directions is null, it could be skipped:
864 // myXDirection.Dot(myYDirection) = 0
866 for(Standard_Integer i = 0; i < 3; ++i)
868 // Length of the second segment
869 Standard_Real aLSegm2 = 0;
870 for(Standard_Integer j = 0; j < 3; ++j)
871 aLSegm2 += theOther.myHDims[j] * Abs(theOther.myAxes[j].Dot(myAxes[i]));
873 // Distance between projected centers
874 Standard_Real aDistCC = Abs(D.Dot(myAxes[i]));
876 if(aDistCC > myHDims[i] + aLSegm2)
877 return Standard_True;
880 // Check the axes of the Other box, i.e. L is one of theOther.myAxes
882 for(Standard_Integer i = 0; i < 3; ++i)
884 // Length of the first segment
885 Standard_Real aLSegm1 = 0.;
886 for(Standard_Integer j = 0; j < 3; ++j)
887 aLSegm1 += myHDims[j] * Abs(myAxes[j].Dot(theOther.myAxes[i]));
889 // Distance between projected centers
890 Standard_Real aDistCC = Abs(D.Dot(theOther.myAxes[i]));
892 if(aDistCC > aLSegm1 + theOther.myHDims[i])
893 return Standard_True;
896 const Standard_Real aTolNull = Epsilon(1.0);
898 // Check the axes produced by the cross products
899 for(Standard_Integer i = 0; i < 3; ++i)
901 for(Standard_Integer j = 0; j < 3; ++j)
904 gp_XYZ aLAxe = myAxes[i].Crossed(theOther.myAxes[j]);
906 const Standard_Real aNorm = aLAxe.Modulus();
912 // Length of the first segment
913 Standard_Real aLSegm1 = 0.;
914 for(Standard_Integer k = 0; k < 3; ++k)
915 aLSegm1 += myHDims[k] * Abs(myAxes[k].Dot(aLAxe));
917 // Length of the second segment
918 Standard_Real aLSegm2 = 0.;
919 for(Standard_Integer k = 0; k < 3; ++k)
920 aLSegm2 += theOther.myHDims[k] * Abs(theOther.myAxes[k].Dot(aLAxe));
922 // Distance between projected centers
923 Standard_Real aDistCC = Abs(D.Dot(aLAxe));
925 if(aDistCC > aLSegm1 + aLSegm2)
926 return Standard_True;
930 return Standard_False;
933 // =======================================================================
936 // =======================================================================
937 Standard_Boolean Bnd_OBB::IsOut(const gp_Pnt& theP) const
939 // 1. Project the point to myAxes[i] (i=0...2).
940 // 2. Check, whether the absolute value of the correspond
941 // projection parameter is greater than myHDims[i].
942 // In this case, IsOut method will return TRUE.
944 const gp_XYZ aRV = theP.XYZ() - myCenter;
946 return ((Abs(myAxes[0].Dot(aRV)) > myHDims[0]) ||
947 (Abs(myAxes[1].Dot(aRV)) > myHDims[1]) ||
948 (Abs(myAxes[2].Dot(aRV)) > myHDims[2]));
951 // =======================================================================
952 // function : IsCompletelyInside
953 // purpose : Checks if every vertex of theOther is completely inside *this
954 // =======================================================================
955 Standard_Boolean Bnd_OBB::IsCompletelyInside(const Bnd_OBB& theOther) const
957 if(IsVoid() || theOther.IsVoid())
958 return Standard_False;
961 theOther.GetVertex(aVert);
962 for(Standard_Integer i = 0; i < 8; i++)
965 return Standard_False;
968 return Standard_True;
971 // =======================================================================
974 // =======================================================================
975 void Bnd_OBB::Add(const gp_Pnt& theP)
979 myCenter = theP.XYZ();
980 myAxes[0] = gp::DX().XYZ();
981 myAxes[1] = gp::DY().XYZ();
982 myAxes[2] = gp::DZ().XYZ();
986 myIsAABox = Standard_True;
993 ReBuild(TColgp_Array1OfPnt(aList[0], 0, 8));
997 // =======================================================================
1000 // =======================================================================
1001 void Bnd_OBB::Add(const Bnd_OBB& theOther)
1003 if (!theOther.IsVoid())
1007 myCenter = theOther.myCenter;
1008 myAxes[0] = theOther.myAxes[0];
1009 myAxes[1] = theOther.myAxes[1];
1010 myAxes[2] = theOther.myAxes[2];
1011 myHDims[0] = theOther.myHDims[0];
1012 myHDims[1] = theOther.myHDims[1];
1013 myHDims[2] = theOther.myHDims[2];
1014 myIsAABox = theOther.myIsAABox;
1019 GetVertex(&aList[0]);
1020 theOther.GetVertex(&aList[8]);
1021 ReBuild(TColgp_Array1OfPnt(aList[0], 0, 15));
1026 //=======================================================================
1027 //function : DumpJson
1029 //=======================================================================
1030 void Bnd_OBB::DumpJson (Standard_OStream& theOStream, const Standard_Integer theDepth) const
1032 DUMP_CLASS_BEGIN (theOStream, Bnd_OBB);
1034 DUMP_FIELD_VALUES_DUMPED (theOStream, theDepth, &myCenter);
1035 DUMP_FIELD_VALUES_DUMPED (theOStream, theDepth, &myAxes[0]);
1036 DUMP_FIELD_VALUES_DUMPED (theOStream, theDepth, &myAxes[1]);
1037 DUMP_FIELD_VALUES_DUMPED (theOStream, theDepth, &myAxes[2]);
1039 DUMP_FIELD_VALUE_NUMERICAL (theOStream, myHDims[0]);
1040 DUMP_FIELD_VALUE_NUMERICAL (theOStream, myHDims[1]);
1041 DUMP_FIELD_VALUE_NUMERICAL (theOStream, myHDims[2]);
1042 DUMP_FIELD_VALUE_NUMERICAL (theOStream, myIsAABox);