OpenCascade/occt
1
0
Fork 0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-04-04 18:06:22 +03:00
Code Packages Releases Activity

0030595: Oriented Bounding Box seems not optimal for some shapes

Browse Source 
Add possibility of construction of the Optimal Oriented Bounding Box from set of points (the case of shape with triangulation).

The interface of the BRepBndLib::AddOBB method is not changed, but the option <theIsOptimal> now controls also the construction of the OBB from Set of points.
The slightly modified DiTo algorithm will be used, checking all possible axes created by the extreme points.
The performance of the construction of the Optimal OBB is lower but the quality is usually much higher (can't be worse by definition).

Test cases for the issue.
This commit is contained in:
emv 2019-04-18 11:17:18 +03:00
parent 6b41f0f335
commit 1bb67d3844
12 changed files with 579 additions and 129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
dox/user_guides/modeling_data/images/modeling_data_obb_125K.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 18 KiB

BIN
dox/user_guides/modeling_data/images/modeling_data_opt_obb_125K.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 18 KiB

BIN
dox/user_guides/modeling_data/images/modeling_data_pca_obb_125K.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 15 KiB

21
dox/user_guides/modeling_data/modeling_data.md
View File

@ -1341,7 +1341,26 @@ Further, let us consider the triangle \f$ T_{0}\left \langle p_{0}, p_{1}, p_{2}
<span>10.</span> Compute the center of OBB and its half dimensions.<br>
<span>11.</span> Create OBB using the center, axes and half dimensions.<br>
This algorithm is implemented in the *Bnd_OBB::ReBuild(...)* method.
@subsubsection occt_modat_6_1_1_opt Creation of Optimal OBB from set of points
For creation of the optimal OBB from set of points the same algorithm as described above is used but with some simplifications in logic and increased computation time.
For the optimal OBB it is necessary to check all possible axes which can be created by the extremal points. And since the extremal points are only valid for the initial axes it is necessary to project the whole set of points on each axis.
This approach usually provides much tighter OBB but the performance is lower. The complexity of the algorithm is still linear and with use of BVH for the set of points it is O(N + C*log(N)).
Here is the example of optimal and not optimal OBB for the model using the set of 125K nodes:
<table align="center">
<tr>
<td>@figure{/user_guides/modeling_data/images/modeling_data_obb_125K.png,"Not optimal OBB by DiTo-14",160}</td>
<td>@figure{/user_guides/modeling_data/images/modeling_data_opt_obb_125K.png,"Optimal OBB by DiTo-14",160}</td>
<td>@figure{/user_guides/modeling_data/images/modeling_data_pca_obb_125K.png,"Not optimal OBB by PCA",160}</td>
</tr>
</table>
Computation of the not optimal OBB in this case took 0.007 sec, optimal - 0.1 sec, which is about 14 times slower. Such performance is comparable to creation of the OBB for this shape by PCA approach (see below) which takes about 0.17 sec.
The computation of optimal OBB is controlled by the same *theIsOptimal* flag in the BRepBndLib::AddOBB method as for PCA algorithm.
These algorithms are implemented in the *Bnd_OBB::ReBuild(...)* method.
@subsubsection occt_modat_6_1_2 Creation of OBB based on Axes of inertia

5
src/BRepBndLib/BRepBndLib.hxx
View File

@ -93,9 +93,8 @@ public:
//! be ignored at all.
//! If theIsShapeToleranceUsed == TRUE then resulting box will be
//! extended on the tolerance of the shape.
//! theIsOptimal flag defines the algorithm for construction of initial
//! Bnd_Box for the second method (if theIsOptimal == TRUE then
//! this box will be created by AddOptimal(...) method).
//! theIsOptimal flag defines whether to look for the more tight
//! OBB for the cost of performance or not.
Standard_EXPORT static
void AddOBB(const TopoDS_Shape& theS,
Bnd_OBB& theOBB,

5
src/BRepBndLib/BRepBndLib_1.cxx
View File

@ -293,6 +293,7 @@ static Standard_Integer IsWCS(const gp_Dir& theDir)
//=======================================================================
static Standard_Boolean CheckPoints(const TopoDS_Shape& theS,
const Standard_Boolean theIsTriangulationUsed,
const Standard_Boolean theIsOptimal,
const Standard_Boolean theIsShapeToleranceUsed,
Bnd_OBB& theOBB)
{
@ -329,7 +330,7 @@ static Standard_Boolean CheckPoints(const TopoDS_Shape& theS,
}
#endif
theOBB.ReBuild(anArrPnts, aPtrArrTol);
theOBB.ReBuild(anArrPnts, aPtrArrTol, theIsOptimal);
return (!theOBB.IsVoid());
}
@ -498,7 +499,7 @@ void BRepBndLib::AddOBB(const TopoDS_Shape& theS,
const Standard_Boolean theIsOptimal,
const Standard_Boolean theIsShapeToleranceUsed)
{
if(CheckPoints(theS, theIsTriangulationUsed, theIsShapeToleranceUsed, theOBB))
if (CheckPoints(theS, theIsTriangulationUsed, theIsOptimal, theIsShapeToleranceUsed, theOBB))
return;
ComputePCA(theS, theOBB, theIsTriangulationUsed, theIsOptimal, theIsShapeToleranceUsed);

4
src/BRepTest/BRepTest_BasicCommands.cxx
View File

@ -1486,7 +1486,9 @@ void BRepTest::BasicCommands(Draw_Interpretor& theCommands)
"\n\t\t: -noTriangulation Force use of exact geometry for calculation"
"\n\t\t: even if triangulation is present."
"\n\t\t: -optimal Force calculation of optimal (more tight) AABB."
"\n\t\t: In case of OBB, applies to initial AABB used in OBB calculation."
"\n\t\t: In case of OBB:"
"\n\t\t: - for PCA approach applies to initial AABB used in OBB calculation"
"\n\t\t: - for DiTo approach modifies the DiTo algorithm to check more axes."
"\n\t\t: -extToler Include tolerance of the shape in the resulting box."
"\n\t\t:"
"\n\t\t: Output options:"

555
src/Bnd/Bnd_OBB.cxx
View File

@ -14,11 +14,177 @@
#include <Bnd_OBB.hxx>
#include <Bnd_B3d.hxx>
#include <Bnd_Tools.hxx>
#include <Bnd_Range.hxx>
#include <BVH_BoxSet.hxx>
#include <BVH_LinearBuilder.hxx>
#include <BVH_Traverse.hxx>
#include <NCollection_Array1.hxx>
#include <Precision.hxx>
#include <TColStd_Array1OfReal.hxx>
//! Auxiliary class to select from the points stored in
//! BVH tree the two points giving the extreme projection
//! parameters on the axis
class OBB_ExtremePointsSelector :
public BVH_Traverse <Standard_Real, 3, BVH_BoxSet <Standard_Real, 3, gp_XYZ>, Bnd_Range>
{
public:
//! Constructor
OBB_ExtremePointsSelector() :
BVH_Traverse <Standard_Real, 3, BVH_BoxSet <Standard_Real, 3, gp_XYZ>, Bnd_Range>(),
myPrmMin (RealLast()),
myPrmMax (RealFirst())
{}
public: //! @name Set axis for projection
//! Sets the axis
void SetAxis (const gp_XYZ& theAxis) { myAxis = theAxis; }
public: //! @name Clears the points from previous runs
//! Clear
void Clear()
{
myPrmMin = RealLast();
myPrmMax = RealFirst();
}
public: //! @name Getting the results
//! Returns the minimal projection parameter
Standard_Real MinPrm() const { return myPrmMin; }
//! Returns the maximal projection parameter
Standard_Real MaxPrm() const { return myPrmMax; }
//! Returns the minimal projection point
const gp_XYZ& MinPnt() const { return myPntMin; }
//! Returns the maximal projection point
const gp_XYZ& MaxPnt() const { return myPntMax; }
public: //! @name Definition of rejection/acceptance rules
//! Defines the rules for node rejection
virtual Standard_Boolean RejectNode (const BVH_Vec3d& theCMin,
const BVH_Vec3d& theCMax,
Bnd_Range& theMetric) const Standard_OVERRIDE
{
if (myPrmMin > myPrmMax)
// No parameters computed yet
return Standard_False;
Standard_Real aPrmMin = myPrmMin, aPrmMax = myPrmMax;
Standard_Boolean isToReject = Standard_True;
// Check if the current node is between already found parameters
for (Standard_Integer i = 0; i < 2; ++i)
{
Standard_Real x = !i ? theCMin.x() : theCMax.x();
for (Standard_Integer j = 0; j < 2; ++j)
{
Standard_Real y = !j ? theCMin.y() : theCMax.y();
for (Standard_Integer k = 0; k < 2; ++k)
{
Standard_Real z = !k ? theCMin.z() : theCMax.z();
Standard_Real aPrm = myAxis.Dot (gp_XYZ (x, y, z));
if (aPrm < aPrmMin)
{
aPrmMin = aPrm;
isToReject = Standard_False;
}
else if (aPrm > aPrmMax)
{
aPrmMax = aPrm;
isToReject = Standard_False;
}
}
}
}
theMetric = Bnd_Range (aPrmMin, aPrmMax);
return isToReject;
}
//! Rules for node rejection by the metric
virtual Standard_Boolean RejectMetric (const Bnd_Range& theMetric) const Standard_OVERRIDE
{
if (myPrmMin > myPrmMax)
// no parameters computed
return Standard_False;
Standard_Real aMin, aMax;
if (!theMetric.GetBounds (aMin, aMax))
// void metric
return Standard_False;
// Check if the box of the branch is inside of the already computed parameters
return aMin > myPrmMin && aMax < myPrmMax;
}
//! Defines the rules for leaf acceptance
virtual Standard_Boolean Accept (const Standard_Integer theIndex,
const Bnd_Range&) Standard_OVERRIDE
{
const gp_XYZ& theLeaf = myBVHSet->Element (theIndex);
Standard_Real aPrm = myAxis.Dot (theLeaf);
if (aPrm < myPrmMin)
{
myPrmMin = aPrm;
myPntMin = theLeaf;
}
if (aPrm > myPrmMax)
{
myPrmMax = aPrm;
myPntMax = theLeaf;
}
return Standard_True;
}
public: //! @name Choosing the best branch
//! Returns true if the metric of the left branch is better than the metric of the right
virtual Standard_Boolean IsMetricBetter (const Bnd_Range& theLeft,
const Bnd_Range& theRight) const Standard_OVERRIDE
{
if (myPrmMin > myPrmMax)
// no parameters computed
return Standard_True;
Standard_Real aMin[2], aMax[2];
if (!theLeft.GetBounds (aMin[0], aMax[0]) ||
!theRight.GetBounds (aMin[1], aMax[1]))
// void metrics
return Standard_True;
// Choose branch with larger extension over computed parameters
Standard_Real anExt[2] = {0.0, 0.0};
for (int i = 0; i < 2; ++i)
{
if (aMin[i] < myPrmMin) anExt[i] += myPrmMin - aMin[i];
if (aMax[i] > myPrmMax) anExt[i] += aMax[i] - myPrmMax;
}
return anExt[0] > anExt[1];
}
protected: //! @name Fields
gp_XYZ myAxis; //!< Axis to project the points to
Standard_Real myPrmMin; //!< Minimal projection parameter
Standard_Real myPrmMax; //!< Maximal projection parameter
gp_XYZ myPntMin; //!< Minimal projection point
gp_XYZ myPntMax; //!< Maximal projection point
};
//! Tool for OBB construction
class OBBTool
{
public:
@ -31,7 +197,8 @@ public:
//! must be available during all time of OBB creation
//! (i.e. while the object of OBBTool exists).
OBBTool(const TColgp_Array1OfPnt& theL,
const TColStd_Array1OfReal *theLT = 0);
const TColStd_Array1OfReal *theLT = 0,
Standard_Boolean theIsOptimal = Standard_False);
//! DiTO algorithm for OBB construction
//! (http://www.idt.mdh.se/~tla/publ/FastOBBs.pdf)
@ -41,6 +208,10 @@ public:
void BuildBox(Bnd_OBB& theBox);
protected:
// Computes the extreme points on the set of Initial axes
void ComputeExtremePoints ();
//! Works with the triangle set by the points in myTriIdx.
//! If theIsBuiltTrg == TRUE, new set of triangles will be
//! recomputed.
@ -71,6 +242,106 @@ protected:
OBBTool& operator=(const OBBTool&);
private:
//! Params structure stores the two values meaning
//! min and max parameters on the axis
struct Params
{
Params() :
_ParamMin(RealLast()), _ParamMax(RealFirst())
{}
Params(Standard_Real theMin, Standard_Real theMax)
: _ParamMin(theMin), _ParamMax(theMax)
{}
Standard_Real _ParamMin;
Standard_Real _ParamMax;
};
//! Computes the Minimal and maximal parameters on the vector
//! connecting the points myLExtremalPoints[theId1] and myLExtremalPoints[theId2]
void ComputeParams (const Standard_Integer theId1,
const Standard_Integer theId2,
Standard_Real &theMin,
Standard_Real &theMax)
{
theMin = myParams[theId1][theId2]._ParamMin;
theMax = myParams[theId1][theId2]._ParamMax;
if (theMin > theMax)
{
FindMinMax ((myLExtremalPoints[theId1] - myLExtremalPoints[theId2]).Normalized(), theMin, theMax);
myParams[theId1][theId2]._ParamMin = myParams[theId2][theId1]._ParamMin = theMin;
myParams[theId1][theId2]._ParamMax = myParams[theId2][theId1]._ParamMax = theMax;
}
}
//! Looks for the min-max parameters on the axis.
//! For optimal case projects all the points on the axis,
//! for not optimal - only the set of extreme points.
void FindMinMax (const gp_XYZ& theAxis,
Standard_Real &theMin,
Standard_Real &theMax)
{
theMin = RealLast(), theMax = RealFirst();
if (myOptimal)
Project (theAxis, theMin, theMax);
else
{
for (Standard_Integer i = 0; i < myNbExtremalPoints; ++i)
{
Standard_Real aPrm = theAxis.Dot (myLExtremalPoints[i]);
if (aPrm < theMin) theMin = aPrm;
if (aPrm > theMax) theMax = aPrm;
}
}
}
//! Projects the set of points on the axis
void Project (const gp_XYZ& theAxis,
Standard_Real& theMin, Standard_Real& theMax,
gp_XYZ* thePntMin = 0, gp_XYZ* thePntMax = 0)
{
theMin = RealLast(), theMax = RealFirst();
if (myOptimal)
{
// Project BVH
OBB_ExtremePointsSelector anExtremePointsSelector;
anExtremePointsSelector.SetBVHSet (myPointBoxSet.get());
anExtremePointsSelector.SetAxis (theAxis);
anExtremePointsSelector.Select();
theMin = anExtremePointsSelector.MinPrm();
theMax = anExtremePointsSelector.MaxPrm();
if (thePntMin) *thePntMin = anExtremePointsSelector.MinPnt();
if (thePntMax) *thePntMax = anExtremePointsSelector.MaxPnt();
}
else
{
// Project all points
for (Standard_Integer iP = myPntsList.Lower(); iP <= myPntsList.Upper(); ++iP)
{
const gp_XYZ& aPoint = myPntsList(iP).XYZ();
const Standard_Real aPrm = theAxis.Dot (aPoint);
if (aPrm < theMin)
{
theMin = aPrm;
if (thePntMin)
*thePntMin = aPoint;
}
if (aPrm > theMax)
{
theMax = aPrm;
if (thePntMax)
*thePntMax = aPoint;
}
}
}
}
private:
//! Number of the initial axes.
static const Standard_Integer myNbInitAxes = 7;
@ -96,6 +367,16 @@ private:
//! The surface area of the OBB
Standard_Real myQualityCriterion;
//! Defines if the OBB should be computed more tight.
//! Takes more time, but the volume is less.
Standard_Boolean myOptimal;
//! Point box set organized with BVH
opencascade::handle<BVH_BoxSet <Standard_Real, 3, gp_XYZ>> myPointBoxSet;
//! Stored min/max parameters for the axes between extremal points
Params myParams[myNbExtremalPoints][myNbExtremalPoints];
};
//=======================================================================
@ -110,7 +391,7 @@ static inline void SetMinMax(Standard_Real* const thePrmArr,
{
thePrmArr[0] = theNewParam;
}
else if(theNewParam > thePrmArr[1])
if(theNewParam > thePrmArr[1])
{
thePrmArr[1] = theNewParam;
}
@ -122,76 +403,102 @@ static inline void SetMinMax(Standard_Real* const thePrmArr,
//=======================================================================
OBBTool::
OBBTool(const TColgp_Array1OfPnt& theL,
const TColStd_Array1OfReal *theLT) :myPntsList(theL),
myListOfTolers(theLT),
myQualityCriterion(RealLast())
const TColStd_Array1OfReal *theLT,
const Standard_Boolean theIsOptimal) : myPntsList(theL),
myListOfTolers(theLT),
myQualityCriterion(RealLast()),
myOptimal (theIsOptimal)
{
if (myOptimal)
{
// Use linear builder for BVH construction with 30 elements in the leaf
opencascade::handle<BVH_LinearBuilder<Standard_Real, 3> > aLBuilder =
new BVH_LinearBuilder<Standard_Real, 3> (30);
myPointBoxSet = new BVH_BoxSet <Standard_Real, 3, gp_XYZ> (aLBuilder);
myPointBoxSet->SetSize(myPntsList.Length());
// Add the points into Set
for (Standard_Integer iP = 0; iP < theL.Length(); ++iP)
{
const gp_Pnt& aP = theL (iP);
Standard_Real aTol = theLT ? theLT->Value(iP) : Precision::Confusion();
BVH_Box <Standard_Real, 3> aBox (BVH_Vec3d (aP.X() - aTol, aP.Y() - aTol, aP.Z() - aTol),
BVH_Vec3d (aP.X() + aTol, aP.Y() + aTol, aP.Z() + aTol));
myPointBoxSet->Add (aP.XYZ(), aBox);
}
myPointBoxSet->Build();
}
ComputeExtremePoints();
}
//=======================================================================
// Function : ComputeExtremePoints
// purpose :
//=======================================================================
void OBBTool::ComputeExtremePoints()
{
// Six initial axes show great quality on the Optimal OBB, plus
// the performance is better (due to the less number of operations).
// But they show worse quality for the not optimal approach.
//const Standard_Real a = (sqrt(5) - 1) / 2.;
//const gp_XYZ anInitialAxes6[myNbInitAxes] = { gp_XYZ (0, 1, a),
// gp_XYZ (0, 1, -a),
// gp_XYZ (1, a, 0),
// gp_XYZ (1, -a, 0),
// gp_XYZ (a, 0, 1),
// gp_XYZ (a, 0, -1) };
const Standard_Real aSqrt3 = Sqrt(3);
// Origin of all initial axis is (0,0,0).
// All axes must be normalized.
const gp_XYZ anInitialAxesArray[myNbInitAxes] = {gp_XYZ(1.0, 0.0, 0.0),
gp_XYZ(0.0, 1.0, 0.0),
gp_XYZ(0.0, 0.0, 1.0),
gp_XYZ(1.0, 1.0, 1.0) / aSqrt3,
gp_XYZ(1.0, 1.0, -1.0) / aSqrt3,
gp_XYZ(1.0, -1.0, 1.0) / aSqrt3,
gp_XYZ(1.0, -1.0, -1.0) / aSqrt3};
const gp_XYZ anInitialAxes7[myNbInitAxes] = { gp_XYZ (1.0, 0.0, 0.0),
gp_XYZ (0.0, 1.0, 0.0),
gp_XYZ (0.0, 0.0, 1.0),
gp_XYZ (1.0, 1.0, 1.0) / aSqrt3,
gp_XYZ (1.0, 1.0, -1.0) / aSqrt3,
gp_XYZ (1.0, -1.0, 1.0) / aSqrt3,
gp_XYZ (1.0, -1.0, -1.0) / aSqrt3 };
// Minimal and maximal point on every axis
const Standard_Integer aNbPoints = 2 * myNbInitAxes;
// Set of initial axes
const gp_XYZ *anInitialAxesArray = anInitialAxes7;
for(Standard_Integer i = 0; i < 5; i++)
{
myTriIdx[i] = INT_MAX;
}
// Min and Max parameter
Standard_Real aParams[aNbPoints];
for(Standard_Integer i = 0; i < aNbPoints; i += 2)
{
aParams[i] = RealLast();
aParams[i + 1] = RealFirst();
}
Standard_Real aParams[myNbExtremalPoints];
// Look for the extremal points (myLExtremalPoints)
for(Standard_Integer i = myPntsList.Lower() ; i <= myPntsList.Upper(); i++)
for (Standard_Integer anAxeInd = 0, aPrmInd = -1; anAxeInd < myNbInitAxes; ++anAxeInd)
{
const gp_XYZ &aCurrPoint = myPntsList(i).XYZ();
for(Standard_Integer anAxeInd = 0, aPrmInd = 0; anAxeInd < myNbInitAxes; anAxeInd++, aPrmInd++)
{
const Standard_Real aParam = aCurrPoint.Dot(anInitialAxesArray[anAxeInd]);
if(aParam < aParams[aPrmInd])
{
myLExtremalPoints[aPrmInd] = aCurrPoint;
aParams[aPrmInd] = aParam;
}
aPrmInd++;
if(aParam > aParams[aPrmInd])
{
myLExtremalPoints[aPrmInd] = aCurrPoint;
aParams[aPrmInd] = aParam;
}
}
Standard_Integer aMinInd = ++aPrmInd, aMaxInd = ++aPrmInd;
aParams[aMinInd] = RealLast();
aParams[aMaxInd] = -RealLast();
Project (anInitialAxesArray[anAxeInd],
aParams[aMinInd], aParams[aMaxInd],
&myLExtremalPoints[aMinInd], &myLExtremalPoints[aMaxInd]);
}
// Compute myTriIdx[0] and myTriIdx[1].
Standard_Real aMaxSqDist = -1.0;
for(Standard_Integer aPrmInd = 0; aPrmInd < aNbPoints; aPrmInd += 2)
// For not optimal box it is necessary to compute the max axis
// created by the maximally distant extreme points
if (!myOptimal)
{
const gp_Pnt &aP1 = myLExtremalPoints[aPrmInd],
&aP2 = myLExtremalPoints[aPrmInd + 1];
const Standard_Real aSqDist = aP1.SquareDistance(aP2);
if(aSqDist > aMaxSqDist)
{
aMaxSqDist = aSqDist;
myTriIdx[0] = aPrmInd;
myTriIdx[1] = aPrmInd + 1;
}
}
for(Standard_Integer i = 0; i < 5; i++)
myTriIdx[i] = INT_MAX;
FillToTriangle3();
// Compute myTriIdx[0] and myTriIdx[1].
Standard_Real aMaxSqDist = -1.0;
for (Standard_Integer aPrmInd = 0; aPrmInd < myNbExtremalPoints; aPrmInd += 2)
{
const gp_Pnt &aP1 = myLExtremalPoints[aPrmInd],
&aP2 = myLExtremalPoints[aPrmInd + 1];
const Standard_Real aSqDist = aP1.SquareDistance(aP2);
if (aSqDist > aMaxSqDist)
{
aMaxSqDist = aSqDist;
myTriIdx[0] = aPrmInd;
myTriIdx[1] = aPrmInd + 1;
}
}
// Compute the maximal axis orthogonal to the found one
FillToTriangle3();
}
}
//=======================================================================
@ -233,6 +540,7 @@ void OBBTool::FillToTriangle5(const gp_XYZ& theNormal,
const gp_XYZ& theBarryCenter)
{
Standard_Real aParams[2] = {0.0, 0.0};
Standard_Integer id3 = -1, id4 = -1;
for(Standard_Integer aPtIdx = 0; aPtIdx < myNbExtremalPoints; aPtIdx++)
{
@ -242,28 +550,24 @@ void OBBTool::FillToTriangle5(const gp_XYZ& theNormal,
const gp_XYZ &aCurrPoint = myLExtremalPoints[aPtIdx];
const Standard_Real aParam = theNormal.Dot(aCurrPoint - theBarryCenter);
if(aParam < aParams[0])
if (aParam < aParams[0])
{
myTriIdx[3] = aPtIdx;
id3 = aPtIdx;
aParams[0] = aParam;
}
else if(aParam > aParams[1])
else if (aParam > aParams[1])
{
myTriIdx[4] = aPtIdx;
id4 = aPtIdx;
aParams[1] = aParam;
}
}
// The points must be in the different sides of the triangle plane.
if(aParams[0] > -Precision::Confusion())
{
myTriIdx[3] = INT_MAX;
}
if (id3 >= 0 && aParams[0] < -Precision::Confusion())
myTriIdx[3] = id3;
if(aParams[1] < Precision::Confusion())
{
myTriIdx[4] = INT_MAX;
}
if (id4 >= 0 && aParams[1] > Precision::Confusion())
myTriIdx[4] = id4;
}
//=======================================================================
@ -279,71 +583,60 @@ void OBBTool::ProcessTriangle(const Standard_Integer theIdx1,
{
const Standard_Integer aNbAxes = 3;
//Some vertex of the triangle
const gp_XYZ aP0 = myLExtremalPoints[theIdx1];
// All axes must be normalized in order to provide correct area computation
// (see ComputeQuality(...) method).
gp_XYZ aYAxis[aNbAxes] = {(myLExtremalPoints[theIdx2] - myLExtremalPoints[theIdx1]),
(myLExtremalPoints[theIdx3] - myLExtremalPoints[theIdx2]),
(myLExtremalPoints[theIdx1] - myLExtremalPoints[theIdx3])};
int ID1[3] = { theIdx2, theIdx3, theIdx1 },
ID2[3] = { theIdx1, theIdx2, theIdx3 };
gp_XYZ aYAxis[aNbAxes] = {(myLExtremalPoints[ID1[0]] - myLExtremalPoints[ID2[0]]),
(myLExtremalPoints[ID1[1]] - myLExtremalPoints[ID2[1]]),
(myLExtremalPoints[ID1[2]] - myLExtremalPoints[ID2[2]])};
// Normal to the triangle plane
gp_XYZ aZAxis = aYAxis[0].Crossed(aYAxis[1]);
Standard_Real aSqMod = aZAxis.SquareModulus();
if(aSqMod < Precision::SquareConfusion())
if (aSqMod < Precision::SquareConfusion())
return;
aZAxis /= Sqrt(aSqMod);
gp_XYZ aXAxis[aNbAxes];
for(Standard_Integer i = 0; i < aNbAxes; i++)
{
for (Standard_Integer i = 0; i < aNbAxes; i++)
aXAxis[i] = aYAxis[i].Crossed(aZAxis).Normalized();
aYAxis[i].Normalize();
}
if(theIsBuiltTrg)
FillToTriangle5(aZAxis, aP0);
if (theIsBuiltTrg)
FillToTriangle5 (aZAxis, myLExtremalPoints[theIdx1]);
// Min and Max parameter
const Standard_Integer aNbPoints = 2 * aNbAxes;
// Compute Min/Max params for ZAxis
Standard_Real aParams[aNbPoints];
FindMinMax (aZAxis, aParams[4], aParams[5]); // Compute params on ZAxis once
Standard_Integer aMinIdx = -1;
for(Standard_Integer anAxeInd = 0; anAxeInd < aNbAxes; anAxeInd++)
{
const gp_XYZ &aAX = aXAxis[anAxeInd],
&aAY = aYAxis[anAxeInd];
Standard_Real aParams[aNbPoints] = {0.0, 0.0, 0.0,
0.0, 0.0, 0.0};
for(Standard_Integer aPtIdx = 0; aPtIdx < myNbExtremalPoints; aPtIdx++)
{
if(aPtIdx == theIdx1)
continue;
const gp_XYZ aCurrPoint = myLExtremalPoints[aPtIdx] - aP0;
SetMinMax(&aParams[0], aAX.Dot(aCurrPoint));
SetMinMax(&aParams[2], aAY.Dot(aCurrPoint));
SetMinMax(&aParams[4], aZAxis.Dot(aCurrPoint));
}
const gp_XYZ &aAX = aXAxis[anAxeInd];
// Compute params on XAxis
FindMinMax (aAX, aParams[0], aParams[1]);
// Compute params on YAxis checking for stored values
ComputeParams (ID1[anAxeInd], ID2[anAxeInd], aParams[2], aParams[3]);
const Standard_Real anArea = ComputeQuality(aParams);
if(anArea < myQualityCriterion)
if (anArea < myQualityCriterion)
{
myQualityCriterion = anArea;
aMinIdx = anAxeInd;
}
}
if(aMinIdx < 0)
if (aMinIdx < 0)
return;
myAxes[0] = aXAxis[aMinIdx];
myAxes[1] = aYAxis[aMinIdx];
myAxes[1] = aYAxis[aMinIdx].Normalized();
myAxes[2] = aZAxis;
}
//=======================================================================
@ -352,20 +645,41 @@ void OBBTool::ProcessTriangle(const Standard_Integer theIdx1,
//=======================================================================
void OBBTool::ProcessDiTetrahedron()
{
ProcessTriangle(myTriIdx[0], myTriIdx[1], myTriIdx[2], Standard_True);
if(myTriIdx[3] <= myNbExtremalPoints)
// To compute the optimal OBB it is necessary to check all possible
// axes created by the extremal points. It is also necessary to project
// all the points on the axis, as for each different axis there will be
// different extremal points.
if (myOptimal)
{
ProcessTriangle(myTriIdx[0], myTriIdx[1], myTriIdx[3], Standard_False);
ProcessTriangle(myTriIdx[1], myTriIdx[2], myTriIdx[3], Standard_False);
ProcessTriangle(myTriIdx[0], myTriIdx[2], myTriIdx[3], Standard_False);
for (Standard_Integer i = 0; i < myNbExtremalPoints - 2; i++)
{
for (Standard_Integer j = i + 1; j < myNbExtremalPoints - 1; j++)
{
for (Standard_Integer k = j + 1; k < myNbExtremalPoints; k++)
{
ProcessTriangle (i, j, k, Standard_False);
}
}
}
}
if(myTriIdx[4] <= myNbExtremalPoints)
else
{
ProcessTriangle(myTriIdx[0], myTriIdx[1], myTriIdx[4], Standard_False);
ProcessTriangle(myTriIdx[1], myTriIdx[2], myTriIdx[4], Standard_False);
ProcessTriangle(myTriIdx[0], myTriIdx[2], myTriIdx[4], Standard_False);
// Use the standard DiTo approach
ProcessTriangle(myTriIdx[0], myTriIdx[1], myTriIdx[2], Standard_True);
if (myTriIdx[3] <= myNbExtremalPoints)
{
ProcessTriangle(myTriIdx[0], myTriIdx[1], myTriIdx[3], Standard_False);
ProcessTriangle(myTriIdx[1], myTriIdx[2], myTriIdx[3], Standard_False);
ProcessTriangle(myTriIdx[0], myTriIdx[2], myTriIdx[3], Standard_False);
}
if (myTriIdx[4] <= myNbExtremalPoints)
{
ProcessTriangle(myTriIdx[0], myTriIdx[1], myTriIdx[4], Standard_False);
ProcessTriangle(myTriIdx[1], myTriIdx[2], myTriIdx[4], Standard_False);
ProcessTriangle(myTriIdx[0], myTriIdx[2], myTriIdx[4], Standard_False);
}
}
}
@ -452,7 +766,8 @@ void OBBTool::BuildBox(Bnd_OBB& theBox)
// purpose : http://www.idt.mdh.se/~tla/publ/
// =======================================================================
void Bnd_OBB::ReBuild(const TColgp_Array1OfPnt& theListOfPoints,
const TColStd_Array1OfReal *theListOfTolerances)
const TColStd_Array1OfReal *theListOfTolerances,
const Standard_Boolean theIsOptimal)
{
switch(theListOfPoints.Length())
{
@ -504,7 +819,7 @@ void Bnd_OBB::ReBuild(const TColgp_Array1OfPnt& theListOfPoints,
break;
}
OBBTool aTool(theListOfPoints, theListOfTolerances);
OBBTool aTool(theListOfPoints, theListOfTolerances, theIsOptimal);
aTool.ProcessDiTetrahedron();
aTool.BuildBox(*this);
}

9
src/Bnd/Bnd_OBB.hxx
View File

@ -95,12 +95,17 @@ public:
myCenter.SetCoord(0.5*(aX2 + aX1), 0.5*(aY2 + aY1), 0.5*(aZ2 + aZ1));
}
//! Created new OBB covering every point in theListOfPoints.
//! Creates new OBB covering every point in theListOfPoints.
//! Tolerance of every such point is set by *theListOfTolerances array.
//! If this array is not void (not null-pointer) then the resulted Bnd_OBB
//! will be enlarged using tolerances of points lying on the box surface.
//! <theIsOptimal> flag defines the mode in which the OBB will be built.
//! Constructing Optimal box takes more time, but the resulting box is usually
//! more tight. In case of construction of Optimal OBB more possible
//! axes are checked.
Standard_EXPORT void ReBuild(const TColgp_Array1OfPnt& theListOfPoints,
const TColStd_Array1OfReal *theListOfTolerances = 0);
const TColStd_Array1OfReal *theListOfTolerances = 0,
const Standard_Boolean theIsOptimal = Standard_False);
//! Sets the center of OBB
void SetCenter(const gp_Pnt& theCenter)

41
tests/bugs/modalg_7/bug30595_1 Normal file
View File

@ -0,0 +1,41 @@
puts "==============================================================="
puts "0030595: Oriented Bounding Box seems not optimal for some shapes"
puts "==============================================================="
puts ""
# average volumes of OBBs on different sets
set OBB_Vol_Exp 615000
# set relative error to 1%
set eps 0.01
restore [locate_data_file bug30595_UC4_P_2K.brep] s1
# build OBB
dchrono s1_time start
bounding s1 -obb -shape bs1 -optimal
dchrono s1_time stop counter s1_OBB
# check volume
checkprops bs1 -v $OBB_Vol_Exp -deps $eps
restore [locate_data_file bug30595_UC4_P_13K.brep] s2
# build OBB
dchrono s2_time start
bounding s2 -obb -shape bs2 -optimal
dchrono s2_time stop counter s2_OBB
# check volume
checkprops bs2 -v $OBB_Vol_Exp -deps $eps
restore [locate_data_file bug30595_UC4_P_125K.brep] s3
# build OBB
dchrono s3_time start
bounding s3 -obb -shape bs3 -optimal
dchrono s3_time stop counter s3_OBB
# check volume
checkprops bs3 -v $OBB_Vol_Exp -deps $eps
smallview +X+Z
donly s3 bs3; fit
checkview -screenshot -2d -path ${imagedir}/${test_image}_xz.png
smallview +X+Y
donly s3 bs3; fit
checkview -screenshot -2d -path ${imagedir}/${test_image}_xy.png

18
tests/bugs/modalg_7/bug30595_2 Normal file
View File

@ -0,0 +1,18 @@
puts "==============================================================="
puts "0030595: Oriented Bounding Box seems not optimal for some shapes"
puts "==============================================================="
puts ""
pload XSDRAW
stepread [locate_data_file bug30595_UC1.stp] s *
incmesh s_1 0.1
dchrono stime start
bounding s_1 -obb -shape bs -optimal
dchrono stime stop counter s_OBB
checkprops bs -v 1.32656e+07 -deps 1.e-5
smallview
donly s_1 bs; fit
checkview -screenshot -2d -path ${imagedir}/${test_image}.png

50
tests/bugs/modalg_7/bug30595_3 Normal file
View File

@ -0,0 +1,50 @@
puts "==============================================================="
puts "0030595: Oriented Bounding Box seems not optimal for some shapes"
puts "==============================================================="
puts ""
# test is the copy of the test case bug29311_2
# but computing the optimal OBB comparing to tight AABB
# with 1.e-6% precision
set NbIters 101
set step [expr 360.0/($NbIters-1) ]
restore [locate_data_file bug29237_no_overlap.rhs.brep] a
# Create AABB for a and put it into "r1" variable
# Draw[]> bounding a -shape r1
# The volume of one AABB is
# Draw[]> vprops r1 1.0e-12 -full
# 32736000.184203226
set Vexp 32736000.184203226
set VMax 0
set MaxIteration 0
for {set i 1} { $i <= $NbIters} { incr i } {
bounding a -obb -shape rr$i -optimal
regexp {Mass +: +([-0-9.+eE]+)} [vprops rr$i 1.0e-12 -full] full Vreal
if { $Vreal > $VMax } {
set VMax $Vreal
set MaxIteration $i
copy a amax
}
if { $i != $NbIters } { trotate a 283 162 317 2 7 9 $step }
}
set aDeltaMax [ expr 100.0*abs($VMax/$Vexp - 1.0) ]
puts "Delta of computation not greater than $aDeltaMax %. Maximal delta is achieved in $MaxIteration iteration. See \"amax\" shape."
if { $aDeltaMax > 1.e-6 } {
puts "Error: The obtained OBB(s) is not precise."
}
smallview
donly amax rr${MaxIteration}
fit
checkview -screenshot -2d -path ${imagedir}/${test_image}.png