Minimum volume oriented bounding ellipsoid (#7377)

This PR adds support for Oriented Bounding Ellipsoids (OBE) to Open3D, complementing the existing bounding volume primitives (AABB and OBB). Oriented bounding ellipsoids provide a tighter fit for certain geometries compared to boxes, which can be beneficial for:

    More accurate spatial representations of point clouds and meshes
    Improved collision detection and computational geometry applications
    Better geometric analysis and understanding of 3D shapes

The implementation uses Khachiyan's algorithm to compute the minimum volume enclosing ellipsoid, providing an efficient approximation of the optimal bounding ellipsoid.

MVOBE support for C++ and Python for point clouds and triangle meshes (both legacy and tensor). All computation is done on CPU. Also add ellipsoid primitive for rendering as either mesh or lineset.

Co-authored-by: github-actions[bot] <github-actions[bot]@users.noreply.github.com>
Co-authored-by: Copilot <copilot@github.com>
Co-authored-by: Sameer Sheorey <41028320+ssheorey@users.noreply.github.com>

Author: 4 people

cpp/open3d/geometry/BoundingVolume.cpp

@@ -15,11 +15,107 @@
 #include "open3d/geometry/Qhull.h"
 #include "open3d/geometry/TriangleMesh.h"
 #include "open3d/t/geometry/kernel/MinimumOBB.h"
+#include "open3d/t/geometry/kernel/MinimumOBE.h"
 #include "open3d/utility/Logging.h"
 
 namespace open3d {
 namespace geometry {
 
+OrientedBoundingEllipsoid& OrientedBoundingEllipsoid::Clear() {
+    center_.setZero();
+    radii_.setZero();
+    R_ = Eigen::Matrix3d::Identity();
+    color_.setOnes();
+    return *this;
+}
+
+bool OrientedBoundingEllipsoid::IsEmpty() const { return Volume() <= 0; }
+
+Eigen::Vector3d OrientedBoundingEllipsoid::GetMinBound() const {
+    auto points = GetEllipsoidPoints();
+    return ComputeMinBound(points);
+}
+
+Eigen::Vector3d OrientedBoundingEllipsoid::GetMaxBound() const {
+    auto points = GetEllipsoidPoints();
+    return ComputeMaxBound(points);
+}
+
+Eigen::Vector3d OrientedBoundingEllipsoid::GetCenter() const { return center_; }
+
+AxisAlignedBoundingBox OrientedBoundingEllipsoid::GetAxisAlignedBoundingBox()
+        const {
+    return AxisAlignedBoundingBox::CreateFromPoints(GetEllipsoidPoints());
+}
+
+OrientedBoundingBox OrientedBoundingEllipsoid::GetOrientedBoundingBox(
+        bool) const {
+    return OrientedBoundingBox::CreateFromPoints(GetEllipsoidPoints());
+}
+
+OrientedBoundingBox OrientedBoundingEllipsoid::GetMinimalOrientedBoundingBox(
+        bool robust) const {
+    return OrientedBoundingBox::CreateFromPoints(GetEllipsoidPoints());
+}
+
+OrientedBoundingEllipsoid& OrientedBoundingEllipsoid::Transform(
+        const Eigen::Matrix4d& transformation) {
+    utility::LogError(
+            "A general transform of an OrientedBoundingEllipsoid is not "
+            "implemented. "
+            "Call Translate, Scale, and Rotate.");
+    return *this;
+}
+
+OrientedBoundingEllipsoid& OrientedBoundingEllipsoid::Translate(
+        const Eigen::Vector3d& translation, bool relative) {
+    if (relative) {
+        center_ += translation;
+    } else {
+        center_ = translation;
+    }
+    return *this;
+}
+
+OrientedBoundingEllipsoid& OrientedBoundingEllipsoid::Scale(
+        const double scale, const Eigen::Vector3d& center) {
+    radii_ *= scale;
+    center_ = scale * (center_ - center) + center;
+    return *this;
+}
+
+OrientedBoundingEllipsoid& OrientedBoundingEllipsoid::Rotate(
+        const Eigen::Matrix3d& R, const Eigen::Vector3d& center) {
+    R_ = R * R_;
+    center_ = R * (center_ - center) + center;
+    return *this;
+}
+
+double OrientedBoundingEllipsoid::Volume() const {
+    return 4 * M_PI * radii_(0) * radii_(1) * radii_(2) / 3.0;
+}
+
+std::vector<Eigen::Vector3d> OrientedBoundingEllipsoid::GetEllipsoidPoints()
+        const {
+    Eigen::Vector3d x_axis = R_ * Eigen::Vector3d(radii_(0), 0, 0);
+    Eigen::Vector3d y_axis = R_ * Eigen::Vector3d(0, radii_(1), 0);
+    Eigen::Vector3d z_axis = R_ * Eigen::Vector3d(0, 0, radii_(2));
+    std::vector<Eigen::Vector3d> points(6);
+    points[0] = center_ + x_axis;
+    points[1] = center_ - x_axis;
+    points[2] = center_ + y_axis;
+    points[3] = center_ - y_axis;
+    points[4] = center_ + z_axis;
+    points[5] = center_ - z_axis;
+    return points;
+}
+
+OrientedBoundingEllipsoid OrientedBoundingEllipsoid::CreateFromPoints(
+        const std::vector<Eigen::Vector3d>& points, bool robust) {
+    return t::geometry::kernel::minimum_obe::ComputeMinimumOBEKhachiyan(points,
+                                                                        robust);
+}
+
 OrientedBoundingBox& OrientedBoundingBox::Clear() {
     center_.setZero();
     extent_.setZero();
@@ -363,4 +459,4 @@ std::vector<size_t> AxisAlignedBoundingBox::GetPointIndicesWithinBoundingBox(
 }
 
 }  // namespace geometry
-}  // namespace open3d
+}  // namespace open3d

cpp/open3d/geometry/BoundingVolume.h

@@ -15,6 +15,119 @@ namespace open3d {
 namespace geometry {
 
 class AxisAlignedBoundingBox;
+class OrientedBoundingBox;
+
+class OrientedBoundingEllipsoid : public Geometry3D {
+public:
+    /// \brief Default constructor.
+    ///
+    /// Creates an empty Oriented Bounding Ellipsoid.
+    OrientedBoundingEllipsoid()
+        : Geometry3D(Geometry::GeometryType::OrientedBoundingEllipsoid),
+          center_(0, 0, 0),
+          R_(Eigen::Matrix3d::Identity()),
+          radii_(0, 0, 0),
+          color_(1, 1, 1) {}
+
+    /// \brief Parameterized constructor.
+    ///
+    /// \param center Specifies the center position of the bounding ellipsoid.
+    /// \param R The rotation matrix specifying the orientation of the
+    /// bounding ellipsoid with the original frame of reference.
+    /// \param radii The radii of the bounding ellipsoid.
+    OrientedBoundingEllipsoid(const Eigen::Vector3d& center,
+                              const Eigen::Matrix3d& R,
+                              const Eigen::Vector3d& radii)
+        : Geometry3D(Geometry::GeometryType::OrientedBoundingEllipsoid),
+          center_(center),
+          R_(R),
+          radii_(radii),
+          color_(1, 1, 1) {}
+    ~OrientedBoundingEllipsoid() override {}
+
+public:
+    OrientedBoundingEllipsoid& Clear() override;
+    bool IsEmpty() const override;
+    virtual Eigen::Vector3d GetMinBound() const override;
+    virtual Eigen::Vector3d GetMaxBound() const override;
+    virtual Eigen::Vector3d GetCenter() const override;
+
+    /// Creates an axis-aligned bounding box around the object.
+    virtual AxisAlignedBoundingBox GetAxisAlignedBoundingBox() const override;
+
+    /// Returns an oriented bounding box around the ellipsoid.
+    virtual OrientedBoundingBox GetOrientedBoundingBox(
+            bool robust) const override;
+
+    /// Returns an oriented bounding box around the ellipsoid.
+    virtual OrientedBoundingBox GetMinimalOrientedBoundingBox(
+            bool robust) const override;
+
+    virtual OrientedBoundingEllipsoid& Transform(
+            const Eigen::Matrix4d& transformation) override;
+    virtual OrientedBoundingEllipsoid& Translate(
+            const Eigen::Vector3d& translation, bool relative = true) override;
+    virtual OrientedBoundingEllipsoid& Scale(
+            const double scale, const Eigen::Vector3d& center) override;
+    virtual OrientedBoundingEllipsoid& Rotate(
+            const Eigen::Matrix3d& R, const Eigen::Vector3d& center) override;
+
+    /// Returns the volume of the bounding box.
+    double Volume() const;
+
+    /** Returns the six critical points of the bounding ellipsoid.
+     * \verbatim
+     *      ------- x
+     *     /|
+     *    / |
+     *   /  | z
+     *  y
+     *                            2
+     *                         .--|---.
+     *                    .--'    |     '--.
+     *               .--'         |          '--.
+     *            .'              |   4           '.
+     *           /                |  /              \
+     *          /                 | /                \
+     *       0 |------------------|-------------------| 1
+     *          \               / |                  /
+     *           \             /  |                 /
+     *            '.         5    |               .'
+     *               '--.         |          .--'
+     *                    '--.    |     .--'
+     *                         '--|---'
+     *                            3
+     * \endverbatim
+     */
+    std::vector<Eigen::Vector3d> GetEllipsoidPoints() const;
+
+    /// Return indices to points that are within the bounding ellipsoid.
+    std::vector<size_t> GetPointIndicesWithinBoundingEllipsoid(
+            const std::vector<Eigen::Vector3d>& points) const;
+
+    /// Creates an oriented bounding ellipsoid using a PCA.
+    /// Note, that this is only an approximation to the minimum oriented
+    /// bounding box that could be computed for example with O'Rourke's
+    /// algorithm (cf. http://cs.smith.edu/~jorourke/Papers/MinVolBox.pdf,
+    /// https://www.geometrictools.com/Documentation/MinimumVolumeBox.pdf)
+    /// \param points The input points
+    /// \param robust If set to true uses a more robust method which works
+    ///               in degenerate cases but introduces noise to the points
+    ///               coordinates.
+    static OrientedBoundingEllipsoid CreateFromPoints(
+            const std::vector<Eigen::Vector3d>& points, bool robust = false);
+
+public:
+    /// The center point of the bounding ellipsoid.
+    Eigen::Vector3d center_;
+    /// The rotation matrix of the bounding ellipsoid to transform the original
+    /// frame of reference to the frame of this ellipsoid.
+    Eigen::Matrix3d R_;
+    /// The radii of the bounding ellipsoid in its frame of reference.
+    Eigen::Vector3d radii_;
+    /// The color of the bounding ellipsoid in RGB.
+    Eigen::Vector3d color_;
+};
 
 /// \class OrientedBoundingBox
 ///
@@ -206,6 +319,7 @@ class AxisAlignedBoundingBox : public Geometry3D {
     /// and orientation as the object
     virtual OrientedBoundingBox GetMinimalOrientedBoundingBox(
             bool robust = false) const override;
+
     virtual AxisAlignedBoundingBox& Transform(
             const Eigen::Matrix4d& transformation) override;
     virtual AxisAlignedBoundingBox& Translate(
@@ -288,4 +402,4 @@ class AxisAlignedBoundingBox : public Geometry3D {
 };
 
 }  // namespace geometry
-}  // namespace open3d
+}  // namespace open3d

cpp/open3d/geometry/Geometry.h

@@ -47,6 +47,8 @@ class Geometry {
         OrientedBoundingBox = 11,
         /// AxisAlignedBoundingBox
         AxisAlignedBoundingBox = 12,
+        /// OrientedBoundingEllipsoid
+        OrientedBoundingEllipsoid = 13,
     };
 
 public:

cpp/open3d/geometry/Geometry3D.h

@@ -19,6 +19,7 @@ namespace geometry {
 
 class AxisAlignedBoundingBox;
 class OrientedBoundingBox;
+class OrientedBoundingEllipsoid;
 
 /// \class Geometry3D
 ///

cpp/open3d/geometry/LineSet.h

@@ -17,6 +17,7 @@ namespace open3d {
 namespace geometry {
 
 class PointCloud;
+class OrientedBoundingEllipsoid;
 class OrientedBoundingBox;
 class AxisAlignedBoundingBox;
 class TriangleMesh;
@@ -133,6 +134,13 @@ class LineSet : public Geometry3D {
     static std::shared_ptr<LineSet> CreateFromAxisAlignedBoundingBox(
             const AxisAlignedBoundingBox &box);
 
+    /// \brief Factory function to create a LineSet from an
+    /// OrientedBoundingEllipsoid.
+    ///
+    /// \param ellipsoid The input bounding ellipsoid.
+    static std::shared_ptr<LineSet> CreateFromOrientedBoundingEllipsoid(
+            const OrientedBoundingEllipsoid &ellipsoid);
+
     /// Factory function to create a LineSet from edges of a triangle mesh.
     ///
     /// \param mesh The input triangle mesh.

cpp/open3d/geometry/LineSetFactory.cpp

@@ -17,9 +17,9 @@ namespace open3d {
 namespace geometry {
 
 std::shared_ptr<LineSet> LineSet::CreateFromPointCloudCorrespondences(
-        const PointCloud &cloud0,
-        const PointCloud &cloud1,
-        const std::vector<std::pair<int, int>> &correspondences) {
+        const PointCloud& cloud0,
+        const PointCloud& cloud1,
+        const std::vector<std::pair<int, int>>& correspondences) {
     auto lineset_ptr = std::make_shared<LineSet>();
     size_t point0_size = cloud0.points_.size();
     size_t point1_size = cloud1.points_.size();
@@ -39,7 +39,7 @@ std::shared_ptr<LineSet> LineSet::CreateFromPointCloudCorrespondences(
 }
 
 std::shared_ptr<LineSet> LineSet::CreateFromTriangleMesh(
-        const TriangleMesh &mesh) {
+        const TriangleMesh& mesh) {
     auto line_set = std::make_shared<LineSet>();
     line_set->points_ = mesh.vertices_;
 
@@ -52,7 +52,7 @@ std::shared_ptr<LineSet> LineSet::CreateFromTriangleMesh(
             line_set->lines_.push_back(Eigen::Vector2i(vidx0, vidx1));
         }
     };
-    for (const auto &triangle : mesh.triangles_) {
+    for (const auto& triangle : mesh.triangles_) {
         InsertEdge(triangle(0), triangle(1));
         InsertEdge(triangle(1), triangle(2));
         InsertEdge(triangle(2), triangle(0));
@@ -61,8 +61,21 @@ std::shared_ptr<LineSet> LineSet::CreateFromTriangleMesh(
     return line_set;
 }
 
+std::shared_ptr<LineSet> LineSet::CreateFromOrientedBoundingEllipsoid(
+        const OrientedBoundingEllipsoid& ellipsoid) {
+    std::shared_ptr<TriangleMesh> obel =
+            geometry::TriangleMesh::CreateEllipsoid(ellipsoid.radii_(0),
+                                                    ellipsoid.radii_(1),
+                                                    ellipsoid.radii_(2));
+    obel->Rotate(ellipsoid.R_, Eigen::Vector3d::Zero());
+    obel->Translate(ellipsoid.center_);
+    auto line_set = CreateFromTriangleMesh(*obel);
+    line_set->PaintUniformColor(ellipsoid.color_);
+    return line_set;
+}
+
 std::shared_ptr<LineSet> LineSet::CreateFromOrientedBoundingBox(
-        const OrientedBoundingBox &box) {
+        const OrientedBoundingBox& box) {
     auto line_set = std::make_shared<LineSet>();
     line_set->points_ = box.GetBoxPoints();
     line_set->lines_.push_back(Eigen::Vector2i(0, 1));
@@ -82,7 +95,7 @@ std::shared_ptr<LineSet> LineSet::CreateFromOrientedBoundingBox(
 }
 
 std::shared_ptr<LineSet> LineSet::CreateFromAxisAlignedBoundingBox(
-        const AxisAlignedBoundingBox &box) {
+        const AxisAlignedBoundingBox& box) {
     auto line_set = std::make_shared<LineSet>();
     line_set->points_ = box.GetBoxPoints();
     line_set->lines_.push_back(Eigen::Vector2i(0, 1));
@@ -101,7 +114,7 @@ std::shared_ptr<LineSet> LineSet::CreateFromAxisAlignedBoundingBox(
     return line_set;
 }
 
-std::shared_ptr<LineSet> LineSet::CreateFromTetraMesh(const TetraMesh &mesh) {
+std::shared_ptr<LineSet> LineSet::CreateFromTetraMesh(const TetraMesh& mesh) {
     auto line_set = std::make_shared<LineSet>();
     line_set->points_ = mesh.vertices_;
 
@@ -114,7 +127,7 @@ std::shared_ptr<LineSet> LineSet::CreateFromTetraMesh(const TetraMesh &mesh) {
             line_set->lines_.push_back(Eigen::Vector2i(vidx0, vidx1));
         }
     };
-    for (const auto &tetra : mesh.tetras_) {
+    for (const auto& tetra : mesh.tetras_) {
         InsertEdge(tetra(0), tetra(1));
         InsertEdge(tetra(1), tetra(2));
         InsertEdge(tetra(2), tetra(0));
@@ -129,8 +142,8 @@ std::shared_ptr<LineSet> LineSet::CreateFromTetraMesh(const TetraMesh &mesh) {
 std::shared_ptr<LineSet> LineSet::CreateCameraVisualization(
         int view_width_px,
         int view_height_px,
-        const Eigen::Matrix3d &intrinsic,
-        const Eigen::Matrix4d &extrinsic,
+        const Eigen::Matrix3d& intrinsic,
+        const Eigen::Matrix4d& extrinsic,
         double scale) {
     Eigen::Matrix4d intrinsic4;
     intrinsic4 << intrinsic(0, 0), intrinsic(0, 1), intrinsic(0, 2), 0.0,
@@ -140,8 +153,8 @@ std::shared_ptr<LineSet> LineSet::CreateCameraVisualization(
     Eigen::Matrix4d m = (intrinsic4 * extrinsic).inverse();
     auto lines = std::make_shared<geometry::LineSet>();
 
-    auto mult = [](const Eigen::Matrix4d &m,
-                   const Eigen::Vector3d &v) -> Eigen::Vector3d {
+    auto mult = [](const Eigen::Matrix4d& m,
+                   const Eigen::Vector3d& v) -> Eigen::Vector3d {
         Eigen::Vector4d v4(v.x(), v.y(), v.z(), 1.0);
         auto result = m * v4;
         return Eigen::Vector3d{result.x() / result.w(), result.y() / result.w(),