Super Fibonacci and Oriented Bounding Boxes (#2472) [skip ci]

* add codesign for mac execs

* better ifdef guard

* simple brute force obb

* wrapper on cgals obb

* tutorial for OBB

* minimal test

* minimal test

* note

* use igl::PI

* expose quantity to optimize over

cmake/igl/igl_add_tutorial.cmake

@@ -17,4 +17,13 @@ function(igl_add_tutorial name)
     )
 
     set_target_properties(${name} PROPERTIES FOLDER Libigl_Tutorials)
+
+    # Do this codesign only on macOS
+    # add_custom_command(TARGET your_target POST_BUILD COMMAND codesign -s - $<TARGET_FILE:your_target>
+    if(APPLE)
+      add_custom_command(TARGET ${name} POST_BUILD
+        COMMAND codesign -s - $<TARGET_FILE:${name}>
+        COMMENT "Codesigning ${name}"
+      )
+    endif()
 endfunction()

include/igl/copyleft/cgal/oriented_bounding_box.cpp

@@ -0,0 +1,44 @@
+#include "oriented_bounding_box.h"
+#include <CGAL/Aff_transformation_3.h>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Optimal_bounding_box/oriented_bounding_box.h>
+
+#include <array>
+#include <iostream>
+
+template <typename DerivedP, typename DerivedR>
+  IGL_INLINE void igl::copyleft::cgal::oriented_bounding_box(
+    const Eigen::MatrixBase<DerivedP>& P,
+    Eigen::PlainObjectBase<DerivedR> & R)
+{
+  typedef typename DerivedP::Scalar Scalar;
+  typedef CGAL::Simple_cartesian<Scalar> K;
+  typedef typename K::Point_3 Point;
+  std::vector<Point> points(P.rows());
+  for (int i = 0; i < P.rows(); ++i)
+  {
+    points[i] = Point(P(i, 0), P(i, 1), P(i, 2));
+  }
+
+  std::array<Point, 8> obb_points;
+  CGAL::Aff_transformation_3<K> affine;
+  CGAL::oriented_bounding_box(
+    points, 
+    affine,
+    CGAL::parameters::use_convex_hull(true));
+  // Convert to Eigen affine transformation
+  R.resize(3,3);
+  for(int i = 0; i < 3; ++i)
+  {
+    for (int j = 0; j < 3; ++j)
+    {
+      R(i, j) = affine.m(i, j);
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit instantiation of template function
+template void igl::copyleft::cgal::oriented_bounding_box<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 3, 3, 0, 3, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 3, 0, 3, 3>>&);
+#endif
+

include/igl/copyleft/cgal/oriented_bounding_box.h

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2025 Alec Jacobson <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/
+
+#ifndef IGL_COPYLEFT_CGAL_ORIENTED_BOUNDING_BOX_H
+#define IGL_COPYLEFT_CGAL_ORIENTED_BOUNDING_BOX_H
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+#include <vector>
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      /// Given a set of points compute the rotation transformation of them such
+      /// that their axis-aligned bounding box is as small as possible.
+      ///
+      /// igl::oriented_bounding_box is often faster and better
+      ///
+      /// @param[in] P  #P by 3 list of point locations
+      /// @param[out] R  rotation matrix
+      template <typename DerivedP, typename DerivedR>
+      IGL_INLINE void oriented_bounding_box(
+        const Eigen::MatrixBase<DerivedP>& P,
+        Eigen::PlainObjectBase<DerivedR> & R);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "oriented_bounding_box.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/point_areas.h

@@ -6,8 +6,8 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/
 
-#ifndef IGL_POINT_AREAS_H
-#define IGL_POINT_AREAS_H
+#ifndef IGL_COPYLEFT_CGAL_POINT_AREAS_H
+#define IGL_COPYLEFT_CGAL_POINT_AREAS_H
 #include "../../igl_inline.h"
 #include <Eigen/Core>
 namespace igl

include/igl/oriented_bounding_box.cpp

@@ -0,0 +1,50 @@
+#include "oriented_bounding_box.h"
+
+#include "super_fibonacci.h"
+#include "parallel_for.h"
+
+template <typename DerivedP, typename DerivedR>
+IGL_INLINE void igl::oriented_bounding_box(
+  const Eigen::MatrixBase<DerivedP>& P,
+  const int n,
+  const OrientedBoundingBoxMinimizeType minimize_type,
+  Eigen::PlainObjectBase<DerivedR> & R)
+{
+  typedef typename DerivedP::Scalar Scalar;
+  Eigen::Matrix<Scalar,Eigen::Dynamic,4,Eigen::RowMajor> Q;
+  igl::super_fibonacci(n-1, Q);
+  Q.conservativeResize(n, 4);
+  Q.row(Q.rows()-1) << 0, 0, 0, 1;
+  Eigen::Matrix<Scalar,Eigen::Dynamic,1> losses(Q.rows());
+  igl::parallel_for(Q.rows(),[&](const int i)
+  {
+    Eigen::Quaternion<Scalar> q(Q(i,3), Q(i,0), Q(i,1), Q(i,2));
+    const auto R = q.toRotationMatrix();
+    const auto PR = (P*R).eval();
+    Eigen::Matrix<Scalar,1,3> min_corner = PR.colwise().minCoeff();
+    Eigen::Matrix<Scalar,1,3> max_corner = PR.colwise().maxCoeff();
+    Eigen::Matrix<Scalar,1,3> diagonal = max_corner - min_corner;
+    switch(minimize_type)
+    {
+      case ORIENTED_BOUNDING_BOX_MINIMIZE_VOLUME:
+        losses(i) = diagonal.prod();
+        break;
+     case ORIENTED_BOUNDING_BOX_MINIMIZE_SURFACE_AREA:
+        losses(i) = 2 * (diagonal(0) * diagonal(1) + diagonal(1) * diagonal(2) + diagonal(0) * diagonal(2));
+        break;
+      case ORIENTED_BOUNDING_BOX_MINIMIZE_DIAGONAL_LENGTH:
+        losses(i) = diagonal.squaredNorm();
+        break;
+      default:
+        assert(false && "Unknown minimize type");
+    }
+  });
+  int bi;
+  const Scalar best_loss = losses.minCoeff(&bi);
+  R = Eigen::Quaternion<Scalar>(Q(bi,3), Q(bi,0), Q(bi,1), Q(bi,2)).toRotationMatrix().eval();
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit instantiation of template function
+template void igl::oriented_bounding_box<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 3, 3, 0, 3, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, int, igl::OrientedBoundingBoxMinimizeType, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 3, 0, 3, 3>>&);
+#endif 

include/igl/oriented_bounding_box.h

@@ -0,0 +1,44 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2025 Alec Jacobson <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/
+
+#ifndef IGL_ORIENTED_BOUNDING_BOX_H
+#define IGL_ORIENTED_BOUNDING_BOX_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+#include <vector>
+namespace igl
+{
+  enum OrientedBoundingBoxMinimizeType
+  {
+    ORIENTED_BOUNDING_BOX_MINIMIZE_VOLUME = 0,
+    ORIENTED_BOUNDING_BOX_MINIMIZE_SURFACE_AREA = 1,
+    ORIENTED_BOUNDING_BOX_MINIMIZE_DIAGONAL_LENGTH = 2,
+    NUM_ORIENTED_BOUNDING_BOX_MINIMIZE_TYPES = 3,
+  };
+  /// Given a set of points compute the rotation transformation of them such
+  /// that their axis-aligned bounding box is as small as possible.
+  ///
+  /// Consider passing the points on the convex hull of original list of points.
+  ///
+  /// @param[in] P  #P by 3 list of point locations
+  /// @param[in] n  number of rotations to try
+  /// @param[in] minimize_type  which quantity to minimize
+  /// @param[out] R  rotation matrix
+  template <typename DerivedP, typename DerivedR>
+  IGL_INLINE void oriented_bounding_box(
+    const Eigen::MatrixBase<DerivedP>& P,
+    const int n,
+    const OrientedBoundingBoxMinimizeType minimize_type,
+    Eigen::PlainObjectBase<DerivedR> & R);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "oriented_bounding_box.cpp"
+#endif
+#endif

include/igl/super_fibonacci.cpp

@@ -0,0 +1,37 @@
+#include "super_fibonacci.h"
+#include "PI.h"
+#include <cmath>
+
+template <typename DerivedQ>
+IGL_INLINE void igl::super_fibonacci(
+  const int n,
+  Eigen::PlainObjectBase<DerivedQ> & Q)
+{
+  typedef typename DerivedQ::Scalar Scalar;
+  // https://marcalexa.github.io/superfibonacci/
+  Scalar dn = 1.0 / (Scalar)n;
+  Scalar mc0 = 1.0 / std::sqrt(2.0);
+  Scalar mc1 = 1.0 / 1.533751168755204288118041;
+
+  Q.resize(n,4);
+  for (int i = 0; i < n; i++)
+  {
+    Scalar s = (Scalar)i+0.5;
+    Scalar ab = 2.0 * igl::PI * s;
+    Scalar alpha = ab * mc0;
+    Scalar beta = ab * mc1;
+    s *= dn;
+    Scalar r = std::sqrt(s);
+    Scalar R = std::sqrt(1.0-s);
+    Q(i,0) = r*std::sin(alpha);
+    Q(i,1) = r*std::cos(alpha);
+    Q(i,2) = R*std::sin(beta);
+    Q(i,3) = R*std::cos(beta);
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::super_fibonacci<Eigen::Matrix<double, -1, 4, 1, -1, 4>>(int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 1, -1, 4>>&);
+template void igl::super_fibonacci<Eigen::Matrix<double, -1, -1, 0, -1, -1>>(int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&);
+#endif

include/igl/super_fibonacci.h

@@ -0,0 +1,33 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2025 Alec Jacobson <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/
+
+#ifndef IGL_SUPER_FIBONACCI_H
+#define IGL_SUPER_FIBONACCI_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+#include <vector>
+
+namespace igl
+{
+  ///   super_fibonacci Generate n quaternions according to "Super-Fibonacci
+  /// Spirals: Fast, Low-Discrepancy Sampling of SO(3)" [Alexa 2021]
+  /// 
+  /// 
+  /// @param[in]  n  number of rotations to generate
+  /// @param[out] q  n by 4list of unit quaternions
+  template <typename DerivedQ>
+  IGL_INLINE void super_fibonacci(
+    const int n,
+    Eigen::PlainObjectBase<DerivedQ> & Q);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "super_fibonacci.cpp"
+#endif
+#endif

tests/include/igl/super_fibonacci.cpp

@@ -0,0 +1,11 @@
+#include <test_common.h>
+#include <igl/super_fibonacci.h>
+
+TEST_CASE("super_fibonacci: simple", "[igl]")
+{
+  Eigen::MatrixXd Q;
+  igl::super_fibonacci(10,Q);
+  REQUIRE(Q.rows() == 10);
+  REQUIRE((Q.rowwise().norm()).eval().maxCoeff() == Approx(1.0));
+  REQUIRE((Q.rowwise().norm()).eval().minCoeff() == Approx(1.0));
+}

tutorial/910_OrientedBoundingBox/main.cpp

@@ -0,0 +1,175 @@
+#include <igl/read_triangle_mesh.h>
+#include <igl/copyleft/cgal/oriented_bounding_box.h>
+#include <igl/copyleft/cgal/convex_hull.h>
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/oriented_bounding_box.h>
+#include <igl/moments.h>
+#include <igl/super_fibonacci.h>
+#include <igl/matlab_format.h>
+#include <igl/write_triangle_mesh.h>
+#include <igl/get_seconds.h>
+#include <igl/parallel_for.h>
+#include <Eigen/Core>
+#include <iostream>
+#include <limits>
+
+
+int main(int argc, char * argv[])
+{
+  IGL_TICTOC_LAMBDA;
+
+  Eigen::MatrixXi F;
+  Eigen::MatrixXd V;
+  // Read in inputs as double precision floating point meshes
+  igl::read_triangle_mesh(
+      argc>1?argv[1]: TUTORIAL_SHARED_PATH "/hand.mesh",V,F);
+
+  V.array() += 1;
+
+  tictoc();
+  {
+    Eigen::RowVector3d min_corner = V.colwise().minCoeff();
+    Eigen::RowVector3d max_corner = V.colwise().maxCoeff();
+  };
+  double t_aabb = tictoc();
+
+  // PCA
+  tictoc();
+  Eigen::RowVector3d mean = V.colwise().mean();
+  Eigen::MatrixXd V_centered = V.rowwise() - mean;
+  Eigen::Matrix3d cov = (V_centered.transpose() * V_centered) / (V.rows() - 1);
+  Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eig(cov);
+  Eigen::Matrix3d PR = eig.eigenvectors();
+  double t_pca = tictoc();
+
+
+  Eigen::Matrix3d R;
+  tictoc();
+  igl::copyleft::cgal::oriented_bounding_box(V, R);
+  double t_cgal = tictoc();
+
+  tictoc();
+  Eigen::MatrixXd W;
+  Eigen::MatrixXi G;
+  igl::copyleft::cgal::convex_hull(V,W,G);
+  Eigen::Matrix3d BR;
+  igl::oriented_bounding_box(W, 50000, igl::ORIENTED_BOUNDING_BOX_MINIMIZE_VOLUME, BR);
+  double t_best = tictoc();
+
+  // Cube mesh
+  Eigen::MatrixXd CV(8,3);
+  CV<< 
+    0,0,0,
+    0,1,0,
+    1,0,0,
+    1,1,0,
+    1,0,1,
+    1,1,1,
+    0,0,1,
+    0,1,1;
+  Eigen::MatrixXi CF(12,3);
+  CF<< 
+    1,2,0,
+    1,3,2,
+    3,4,2,
+    3,5,4,
+    0,4,6,
+    0,2,4,
+    7,3,1,
+    7,5,3,
+    7,0,6,
+    7,1,0,
+    5,6,4,
+    5,7,6;
+  Eigen::MatrixXi CE(12,2);
+  CE<< 
+    0,1,
+    0,2,
+    0,6,
+    1,3,
+    1,7,
+    2,3,
+    2,4,
+    3,5,
+    4,5,
+    4,6,
+    5,7,
+    6,7;
+
+
+  const auto transform_cube = [](const Eigen::MatrixXd & V, const Eigen::MatrixXd & CV)
+  {
+    Eigen::RowVector3d min_corner = V.colwise().minCoeff();
+    Eigen::RowVector3d max_corner = V.colwise().maxCoeff();
+    Eigen::RowVector3d diff = max_corner - min_corner;
+    Eigen::MatrixXd AV = CV;
+    AV.array().rowwise() *= diff.array();
+    AV.rowwise() += min_corner;
+    return AV;
+  };
+  Eigen::MatrixXd AV = transform_cube(V, CV);
+
+  Eigen::MatrixXd OV = transform_cube((V*R).eval(), CV);
+  OV = (OV*R.transpose()).eval();
+
+  Eigen::MatrixXd PV = transform_cube((V*PR).eval(), CV);
+  PV = (PV*PR.transpose()).eval();
+
+  Eigen::MatrixXd BV = transform_cube((V*BR).eval(), CV);
+  BV = (BV*BR.transpose()).eval();
+
+  const auto volume = [](const Eigen::MatrixXd & V, const Eigen::MatrixXi & F)
+  {
+    Eigen::Vector3d m1;
+    Eigen::Matrix3d m2;
+    double vol = 0;
+    igl::moments(V,F,vol,m1,m2);
+    return vol;
+  };
+  double Avol = volume(AV,CF);
+  double Pvol = volume(PV,CF);
+  double Ovol = volume(OV,CF);
+  double Bvol = volume(BV,CF);
+  printf("method: %20s %20s  %20s %20s\n", "AABB", "PCA", "CGAL", "super_fibonacci");
+  printf("volume: %20g %20g  %20g %20g\n", Avol, Pvol, Ovol, Bvol);
+  printf("time:   %20g %20g  %20g %20g\n", t_aabb, t_pca, t_cgal, t_best);
+
+
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(V,F);
+  viewer.data().set_face_based(true);
+
+  viewer.append_mesh();
+  viewer.data().set_mesh(AV,CF);
+  viewer.data().set_face_based(true);
+  viewer.data().set_colors(Eigen::RowVector4d(0.43064,0.77072,0.51013,0.5));
+  viewer.data().set_edges(AV,CE,Eigen::RowVector3d(0,0,0));
+  viewer.data().show_lines = false;
+
+  viewer.append_mesh();
+  viewer.data().set_mesh(PV,CF);
+  viewer.data().set_face_based(true);
+  viewer.data().set_colors(Eigen::RowVector4d(0.9374,0.54164,0.66868,0.5));
+  viewer.data().set_edges(PV,CE,Eigen::RowVector3d(0,0,0));
+  viewer.data().show_lines = false;
+
+  viewer.append_mesh();
+  viewer.data().set_mesh(OV,CF);
+  viewer.data().set_face_based(true);
+  viewer.data().set_colors(Eigen::RowVector4d(0.91191,0.60119,0.32995,0.5));
+  viewer.data().show_lines = false;
+
+  viewer.append_mesh();
+  viewer.data().set_mesh(BV,CF);
+  viewer.data().set_face_based(true);
+  viewer.data().set_colors(Eigen::RowVector4d(0.34481,0.72126,0.96535,0.5));
+  viewer.data().show_lines = false;
+  viewer.data().set_edges(OV,CE,Eigen::RowVector3d(0,0,0));
+
+  // set background to white
+  viewer.core().background_color = Eigen::Vector4f(1,1,1,1);
+
+  viewer.data().show_lines = false;
+  viewer.launch();
+}
+

tutorial/CMakeLists.txt

@@ -129,4 +129,5 @@ if(LIBIGL_TUTORIALS_CHAPTER9)
     igl_add_tutorial(907_DynamicAABB igl::glfw)
     igl_add_tutorial(908_IntersectionBlockingDecimation igl::glfw igl::predicates)
     igl_add_tutorial(909_BatchMarchingCubes igl::glfw)
+    igl_add_tutorial(910_OrientedBoundingBox igl::glfw igl_copyleft::cgal)
 endif()