Box.hxx

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2026, Davide Stocco and Enrico Bertolazzi.                  *
 *                                                                           *
 * The AABBtree project is distributed under the BSD 2-Clause License.       *
 *                                                                           *
 * Davide Stocco Enrico Bertolazzi                                           *
 * University of Trento University of Trento                                 *
 * davide.stocco@unitn.it                         enrico.bertolazzi@unitn.it *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
// DISCLAIMER: The code in this file is a modified version of the Eigen library.
 
#pragma once
 
#ifndef AABBTREE_BOX_HXX
#define AABBTREE_BOX_HXX
 
#include "AABBtree/Ray.hxx"
 
namespace AABBtree {

template <typename Real, Integer N> class Ray;

template <typename Real, Integer N> class Box {
  static_assert(std::is_floating_point<Real>::value,
                "Box Real type must be a floating-point type.");
  static_assert(N > 0, "Box dimension must be positive.");
 
  using Point = AABBtree::Point<Real, N>;
  using Vector = AABBtree::Vector<Real, N>;
  using BoxUniquePtr = AABBtree::BoxUniquePtr<Real, N>;
  using BoxUniquePtrList = AABBtree::BoxUniquePtrList<Real, N>;
 
  Point
      m_min; 
  Point
      m_max; 
 
public:
  ~Box() = default;

  Box() { this->set_empty(); }

  Box(Box const &b) : m_min(b.m_min), m_max(b.m_max) {}

  Box(Point const &t_min, Point const &t_max) : m_min(t_min), m_max(t_max) {}

  Box(Point const &p) : m_min(p), m_max(p) {}

  template <typename = std::enable_if<N == 1>>
  Box(Real const x_min, Real const x_max) : m_min(x_min), m_max(x_max) {}

  template <typename = std::enable_if<N == 2>>
  Box(Real const x_min, Real const y_min, Real const x_max, Real const y_max)
      : m_min(x_min, y_min), m_max(x_max, y_max) {}

  template <typename = std::enable_if<N == 3>>
  Box(Real const x_min, Real const y_min, Real const z_min, Real const x_max,
      Real const y_max, Real const z_max)
      : m_min(x_min, y_min, z_min), m_max(x_max, y_max, z_max) {}

  template <typename OtherReal>
  explicit Box(Box<OtherReal, N> const &b)
      : m_min(b.m_min.template cast<Real>()),
        m_max(b.m_max.template cast<Real>()) {}

  Box &operator=(Box const &b) = default;

  template <typename NewReal> Box<NewReal, N> cast() const {
    if constexpr (std::is_same<Real, NewReal>::value) {
      return *this;
    }
    return Box<NewReal, N>(m_min.template cast<NewReal>(),
                           m_max.template cast<NewReal>());
  }

  Point const &min() const { return m_min; }

  Point &min() { return m_min; }

  Point const &max() const { return m_max; }

  Point &max() { return m_max; }

  void reorder() {
    for (Integer i{0}; i < N; ++i) {
      if (m_min[i] > m_max[i]) {
        std::swap(m_min[i], m_max[i]);
      }
    }
  }

  bool is_empty() const { return (m_max.array() < m_min.array()).any(); }

  bool
  is_approx(Box const &b,
            Real const tol = Eigen::NumTraits<Real>::dummy_precision()) const {
    return m_min.isApprox(b.m_min, tol) && m_max.isApprox(b.m_max, tol);
  }

  bool is_degenerate(
      Real const tol = Eigen::NumTraits<Real>::dummy_precision()) const {
    return m_min.isApprox(m_max, tol);
  }

  Integer longest_axis() const {
    Vector sizes{m_max - m_min};
    return std::distance(sizes.data(),
                         std::max_element(sizes.data(), sizes.data() + N));
  }

  Integer longest_axis(Real &max_length, Real &mid_point) const {
    Integer axis{-1};
    max_length = -1;
    for (Integer i{0}; i < N; ++i) {
      if (Real const tmp_length{m_max(i) - m_min(i)}; max_length < tmp_length) {
        max_length = tmp_length;
        mid_point = 0.5 * (m_max(i) + m_min(i));
        axis = i;
      }
    }
    return axis;
  }

  void sort_axes_length(Vector &sizes, Eigen::Vector<Integer, N> &ipos) const {
    sizes = m_max - m_min;
    ipos = Eigen::Vector<Integer, N>::LinSpaced(N, 0, N - 1);
    std::sort(ipos.data(), ipos.data() + N,
              [&sizes](Integer i, Integer j) { return sizes[i] > sizes[j]; });
  }

  Point baricenter() const { return 0.5 * (m_min + m_max); }

  Real baricenter(Integer i) const { return 0.5 * (m_min(i) + m_max(i)); }

  Real volume() const { return (m_max - m_min).prod(); }

  Real surface() const {
    Vector sizes{m_max - m_min};
    Real area{0};
    for (Integer i{0}; i < N; ++i) {
      Real prod{1.0};
      for (Integer j{0}; j < N; ++j)
        if (j != i)
          prod *= sizes[j];
      area += 2.0 * prod;
    }
    return area;
  }

  Vector diagonal() const { return m_max - m_min; }

  void set_degenerate(Point const &p) {
    m_min = p;
    m_max = p;
  }

  void set_empty() {
    m_min.setConstant(+std::numeric_limits<Real>::infinity());
    m_max.setConstant(-std::numeric_limits<Real>::infinity());
  }

  bool intersects(Box const &b) const {
    // return (m_min.array() <= b.m_max.array()).all() && (b.m_min.array() <=
    // m_max.array()).all();
    for (int d = 0; d < m_min.size(); ++d)
      if (m_min[d] > b.m_max[d] || b.m_min[d] > m_max[d])
        return false;
    return true;
  }

  template <Integer d = 2> bool quasi_intersects(Box const &b) const {
    return (m_min.template head<d>().array() <=
            b.m_max.template head<d>().array())
               .all() &&
           (b.m_min.template head<d>().array() <=
            m_max.template head<d>().array())
               .all();
  }

  bool intersect(Box const &b_in, Box &b_out) const {
    b_out.m_min = m_min.cwiseMax(b_in.m_min);
    b_out.m_max = m_max.cwiseMin(b_in.m_max);
    return (b_out.m_min.array() <= b_out.m_max.array()).all();
  }

  Box merged(Box const &b) const {
    return Box(m_min.cwiseMin(b.m_min), m_max.cwiseMax(b.m_max));
  }

  Box &translate(Vector const &t) {
    m_min += t;
    m_max += t;
    return *this;
  }

  Box translated(Vector const &t) const { return Box(m_min + t, m_max + t); }

  template <typename Transform> Box transformed(Transform const &t) const {
    Box result(m_min.transform(t), m_max.transform(t));
    result.reorder();
    return result;
  }

  template <typename Transform> Box &transform(Transform const &t) {
    m_min.transform(t);
    m_max.transform(t);
    this->reorder();
    return *this;
  }

  enum class Side : Integer { LEFT = -1, INSIDE = 0, RIGHT = +1 };

  Side which_side(Real const x, Real const tol, Integer const dim) const {
    Real const x_min{m_min(dim)};
    Real const x_max{m_max(dim)};
    bool const on_left{x_max < x + tol};
    bool const on_right{x - tol < x_min};
    if (on_left && on_right) {
      return 0.5 * (x_min + x_max) < x ? Side::LEFT : Side::RIGHT;
    } else if (on_left) {
      return Side::LEFT;
    } else if (on_right) {
      return Side::RIGHT;
    } else {
      return Side::INSIDE;
    }
  }

  bool contains(Point const &p) const {
    return (m_min.array() <= p.array()).all() &&
           (p.array() <= m_max.array()).all();
  }

  bool intersects(Point const &p) const {
    return (m_min.array() <= p.array()).all() &&
           (p.array() <= m_max.array()).all();
  }

  template <Integer d = 2> bool quasi_intersects(Point const &p) const {
    return (m_min.template head<d>().array() <= p.template head<d>().array()) &&
           (p.template head<d>().array() <= m_max.template head<d>().array());
  }

  bool contains(Box const &b) const {
    return (m_min.array() <= b.m_min.array()).all() &&
           (b.m_max.array() <= m_max.array()).all();
  }

  template <typename Derived> Box &extend(Point const &p) {
    m_min = m_min.cwiseMin(p);
    m_max = m_max.cwiseMax(p);
    return *this;
  }

  Box &extend(Box const &b) {
    m_min = m_min.cwiseMin(b.m_min);
    m_max = m_max.cwiseMax(b.m_max);
    return *this;
  }

  Real squared_interior_distance(Point const &p) const {
    Real dist2{0.0};
    for (Integer i{0}; i < N; ++i) {
      if (m_min[i] > p[i]) {
        Real aux{m_min[i] - p[i]};
        dist2 += aux * aux;
      } else if (p[i] > m_max[i]) {
        Real aux{p[i] - m_max[i]};
        dist2 += aux * aux;
      }
    }
    return dist2;
  }

  Real squared_interior_distance(Point const &p, Point &c) const {
    if (this->contains(p))
      return 0;
    c = p.cwiseMax(m_min).cwiseMin(m_max);
    return (p - c).squaredNorm();
  }

  Real interior_distance(Point const &p) const {
    return std::sqrt(this->squared_interior_distance(p));
  }

  Real interior_distance(Point const &p, Point &c) const {
    return std::sqrt(this->squared_interior_distance(p, c));
  }

  Real squared_exterior_distance(Point const &p) const {
    Real dist2{0.0};
    for (Integer i{0}; i < N; ++i) {
      Real const aux1{std::abs(m_min[i] - p[i])};
      Real const aux2{std::abs(m_max[i] - p[i])};
      Real const aux3{std::max(aux1, aux2)};
      dist2 += aux3 * aux3;
    }
    return dist2;
  }

  Real squared_exterior_distance(Point const &p, Point &f) const {
    for (Integer i{0}; i < N; ++i) {
      Real const aux1{std::abs(m_min[i] - p[i])};
      Real const aux2{std::abs(m_max[i] - p[i])};
      f[i] = aux1 > aux2 ? m_min[i] : m_max[i];
    }
    return (f - p).squaredNorm();
  }

  Real exterior_distance(Point const &p) const {
    return std::sqrt(this->squared_exterior_distance(p));
  }

  Real exterior_distance(Point const &p, Point &f) const {
    return std::sqrt(this->squared_exterior_distance(p, f));
  }

  Real squared_interior_distance(Box const &b) const {
    if (this->intersects(b))
      return 0;
    Real dist2{0.0};
    for (Integer i{0}; i < N; ++i) {
      Real aux{m_min[i] - b.m_max[i]};
      if (aux < 0.0) {
        aux = b.m_min[i] - m_max[i];
        if (aux < 0.0) {
          continue;
        }
      }
      dist2 += aux * aux;
    }
    return dist2;
  }

  Real squared_interior_distance(Box const &b, Point &p1, Point &p2) const {
    if (this->intersects(b)) {
      return 0.0;
    }
    for (Integer i{0}; i < N; ++i) {
      if (m_min[i] > b.m_max[i]) {
        p1[i] = b.m_max[i];
        p2[i] = m_min[i];
      } else if (b.m_min[i] > m_max[i]) {
        p1[i] = m_max[i];
        p2[i] = b.m_min[i];
      } else {
        p1[i] = p2[i] = 0.5 * (std::min(m_max[i], b.m_max[i]) +
                               std::max(m_min[i], b.m_min[i]));
      }
    }
    return (p2 - p1).squaredNorm();
  }

  Real interior_distance(Box const &b) const {
    return std::sqrt(this->squared_interior_distance(b));
  }

  Real interior_distance(Box const &b, Point &p1, Point &p2) const {
    return std::sqrt(this->squared_interior_distance(b, p1, p2));
  }

  Real squared_exterior_distance(Box const &b) const {
    Real dist2{Real(0)};
    for (Integer i{0}; i < N; ++i) {
      Real const aux{std::max(m_max[i], b.m_max[i]) -
                     std::min(m_min[i], b.m_min[i])};
      dist2 += aux * aux;
    }
    return dist2;
  }

  Real squared_exterior_distance(Box const &b, Point &p1, Point &p2) const {
    p1 = b.m_min.cwiseMin(m_min);
    p2 = b.m_max.cwiseMax(m_max);
    return (p2 - p1).squaredNorm();
  }

  Real exterior_distance(Box const &b) const {
    return std::sqrt(this->squared_exterior_distance(b));
  }

  Real exterior_distance(Box const &b, Point &p1, Point &p2) const {
    return std::sqrt(this->squared_exterior_distance(b, p1, p2));
  }

  bool intersects(Ray<Real, N> const &r,
                  Real tol = Eigen::NumTraits<Real>::dummy_precision()) const {
    return r.intersects(*this, tol);
  }

  bool intersect(Ray<Real, N> const &r, Point &c, Point &f,
                 Real tol = Eigen::NumTraits<Real>::dummy_precision()) const {
    return r.intersect(*this, c, f, tol);
  }

  Real squared_interior_distance(
      Ray<Real, N> const &r,
      Real tol = Eigen::NumTraits<Real>::dummy_precision()) const {
    return r.squared_interior_distance(*this, tol);
  }

  Real squared_interior_distance(
      Ray<Real, N> const &r, Point &p1, Point &p2,
      Real tol = Eigen::NumTraits<Real>::dummy_precision()) const {
    return r.squared_interior_distance(*this, p2, p1, tol);
  }

  Real interior_distance(
      Ray<Real, N> const &r,
      Real tol = Eigen::NumTraits<Real>::dummy_precision()) const {
    return r.interior_distance(*this, tol);
  }

  Real interior_distance(
      Ray<Real, N> const &r, Point &p1, Point &p2,
      Real tol = Eigen::NumTraits<Real>::dummy_precision()) const {
    return r.interior_distance(*this, p2, p1, tol);
  }

  Real squared_exterior_distance(
      Ray<Real, N> const &r,
      Real tol = Eigen::NumTraits<Real>::dummy_precision()) const {
    return r.squared_exterior_distance(*this, tol);
  }

  Real squared_exterior_distance(
      Ray<Real, N> const &r, Point &p1, Point &p2,
      Real tol = Eigen::NumTraits<Real>::dummy_precision()) const {
    return r.squared_exterior_distance(*this, p2, p1, tol);
  }

  Real exterior_distance(
      Ray<Real, N> const &r,
      Real tol = Eigen::NumTraits<Real>::dummy_precision()) const {
    return r.exterior_distance(*this, tol);
  }

  Real exterior_distance(
      Ray<Real, N> const &r, Point &p1, Point &p2,
      Real tol = Eigen::NumTraits<Real>::dummy_precision()) const {
    return r.exterior_distance(*this, p2, p1, tol);
  }

  void print(std::ostream &os) const {
    os << "BOX INFO" << std::endl
       << "\tmin = " << m_min.transpose() << std::endl
       << "\tmax = " << m_max.transpose() << std::endl;
  }
 
}; // class Box

template <typename Real, Integer N>
std::ostream &operator<<(std::ostream &os, Box<Real, N> const &b) {
  b.print(os);
  return os;
}
 
} // namespace AABBtree
 
#endif // AABBTREE_BOX_HXX