Box.hxx

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, 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 *
 * e-mail: davide.stocco@unitn.it                             e-mail: 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.");
 
    // Constants for numerical computations
    constexpr static Real EPS{std::numeric_limits<Real>::epsilon()};  
    constexpr static Real INF{std::numeric_limits<Real>::infinity()}; 
    constexpr static Real DUMMY_TOL{EPS*static_cast<Real>(100.0)};    
 
    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 /*= DUMMY_TOL*/) const
    {return m_min.isApprox(b.m_min, tol) && m_max.isApprox(b.m_max, tol);}

    bool is_degenerate(Real const tol /*= DUMMY_TOL*/) 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(+INF);
      m_max.setConstant(-INF);
    }

    bool intersects(Box const & b) const {
      return (m_min.array() <= b.m_max.array()).all() && (b.m_min.array() <= m_max.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 mi       { m_min(dim)   };
      Real const ma       { m_max(dim)   };
      bool const on_left  { ma < x + tol };
      bool const on_right { x - tol < mi };
      if (on_left && on_right) {
        return (mi+ma)/2 < x ? Side::LEFT : Side::RIGHT;
      }
      if (on_left) return Side::LEFT;
      if (on_right) return Side::RIGHT;
      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();
    }

    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 = DUMMY_TOL) const
    {return r.intersects(*this, tol);}

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

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

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

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

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

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

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

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

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

    void print(std::ostream & os) const {
      os <<
        "────────────────────────────────────────────────────────────────────────────────\n" <<
        "BOX INFO\n" <<
        "\tmin = " << m_min.transpose() << '\n' <<
        "\tmax = " << m_max.transpose() << '\n' <<
        "────────────────────────────────────────────────────────────────────────────────\n";
    }
 
  }; // 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