0030595: Oriented Bounding Box seems not optimal for some shapes
authoremv <emv@opencascade.com>
Thu, 18 Apr 2019 08:17:18 +0000 (11:17 +0300)
committeremv <emv@opencascade.com>
Mon, 17 Jun 2019 06:26:54 +0000 (09:26 +0300)
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.

12 files changed:
dox/user_guides/modeling_data/images/modeling_data_obb_125K.png [new file with mode: 0644]
dox/user_guides/modeling_data/images/modeling_data_opt_obb_125K.png [new file with mode: 0644]
dox/user_guides/modeling_data/images/modeling_data_pca_obb_125K.png [new file with mode: 0644]
dox/user_guides/modeling_data/modeling_data.md
src/BRepBndLib/BRepBndLib.hxx
src/BRepBndLib/BRepBndLib_1.cxx
src/BRepTest/BRepTest_BasicCommands.cxx
src/Bnd/Bnd_OBB.cxx
src/Bnd/Bnd_OBB.hxx
tests/bugs/modalg_7/bug30595_1 [new file with mode: 0644]
tests/bugs/modalg_7/bug30595_2 [new file with mode: 0644]
tests/bugs/modalg_7/bug30595_3 [new file with mode: 0644]

diff --git a/dox/user_guides/modeling_data/images/modeling_data_obb_125K.png b/dox/user_guides/modeling_data/images/modeling_data_obb_125K.png
new file mode 100644 (file)
index 0000000..439458a
Binary files /dev/null and b/dox/user_guides/modeling_data/images/modeling_data_obb_125K.png differ
diff --git a/dox/user_guides/modeling_data/images/modeling_data_opt_obb_125K.png b/dox/user_guides/modeling_data/images/modeling_data_opt_obb_125K.png
new file mode 100644 (file)
index 0000000..b245a56
Binary files /dev/null and b/dox/user_guides/modeling_data/images/modeling_data_opt_obb_125K.png differ
diff --git a/dox/user_guides/modeling_data/images/modeling_data_pca_obb_125K.png b/dox/user_guides/modeling_data/images/modeling_data_pca_obb_125K.png
new file mode 100644 (file)
index 0000000..70f8d18
Binary files /dev/null and b/dox/user_guides/modeling_data/images/modeling_data_pca_obb_125K.png differ
index 5c31123..912e0a1 100644 (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
 
index 6f34387..c643cca 100644 (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,
index c70c2df..1dbf954 100644 (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);
index 7e7b36a..c128542 100644 (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:"
index e9b7f36..0120cfe 100644 (file)
 
 #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)
 {
-  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};
-
-  // Minimal and maximal point on every axis
-  const Standard_Integer aNbPoints = 2 * myNbInitAxes;
-
-  for(Standard_Integer i = 0; i < 5; i++)
+  if (myOptimal)
   {
-    myTriIdx[i] = INT_MAX;
+    // 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);
+  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 };
+
+  // Set of initial axes
+  const gp_XYZ *anInitialAxesArray = anInitialAxes7;
+
   // Min and Max parameter
-  Standard_Real aParams[aNbPoints];
-  for(Standard_Integer i = 0; i < aNbPoints; i += 2)
+  Standard_Real aParams[myNbExtremalPoints];
+  // Look for the extremal points (myLExtremalPoints)
+  for (Standard_Integer anAxeInd = 0, aPrmInd = -1; anAxeInd < myNbInitAxes; ++anAxeInd)
   {
-    aParams[i] = RealLast();
-    aParams[i + 1] = RealFirst();
+    Standard_Integer aMinInd = ++aPrmInd, aMaxInd = ++aPrmInd;
+    aParams[aMinInd] = RealLast();
+    aParams[aMaxInd] = -RealLast();
+    Project (anInitialAxesArray[anAxeInd],
+             aParams[aMinInd], aParams[aMaxInd],
+             &myLExtremalPoints[aMinInd], &myLExtremalPoints[aMaxInd]);
   }
 
-  // Look for the extremal points (myLExtremalPoints)
-  for(Standard_Integer i = myPntsList.Lower() ; i <= myPntsList.Upper(); i++)
+  // For not optimal box it is necessary to compute the max axis
+  // created by the maximally distant extreme points
+  if (!myOptimal)
   {
-    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++;
+    for(Standard_Integer i = 0; i < 5; i++)
+      myTriIdx[i] = INT_MAX;
 
-      if(aParam > aParams[aPrmInd])
+    // 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)
       {
-        myLExtremalPoints[aPrmInd] = aCurrPoint;
-        aParams[aPrmInd] = aParam;
+        aMaxSqDist = aSqDist;
+        myTriIdx[0] = aPrmInd;
+        myTriIdx[1] = aPrmInd + 1;
       }
     }
-  }
 
-  // Compute myTriIdx[0] and myTriIdx[1].
-
-  Standard_Real aMaxSqDist = -1.0;
-  for(Standard_Integer aPrmInd = 0; aPrmInd < aNbPoints; 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();
   }
-
-  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);
 }
index e561d2b..ffed615 100644 (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)
diff --git a/tests/bugs/modalg_7/bug30595_1 b/tests/bugs/modalg_7/bug30595_1
new file mode 100644 (file)
index 0000000..61d7803
--- /dev/null
@@ -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
diff --git a/tests/bugs/modalg_7/bug30595_2 b/tests/bugs/modalg_7/bug30595_2
new file mode 100644 (file)
index 0000000..a3ae8a3
--- /dev/null
@@ -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
diff --git a/tests/bugs/modalg_7/bug30595_3 b/tests/bugs/modalg_7/bug30595_3
new file mode 100644 (file)
index 0000000..44095d1
--- /dev/null
@@ -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