bounding_box.hpp

Go to the documentation of this file.

 
#ifndef NIHU_BOUNDING_BOX_HPP_INCLUDED
#define NIHU_BOUNDING_BOX_HPP_INCLUDED
 
#include <Eigen/Dense>
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <iostream>
 
namespace NiHu
{
namespace fmm
{
template <size_t Dim, class Scalar = double>
class bounding_box
{
public:
    static size_t const dimension = Dim;
    using scalar_t = Scalar;
    using location_t = Eigen::Matrix<scalar_t, dimension, 1>;
 
    bounding_box(
        location_t const &center = location_t::Zero(),
        scalar_t diameter = 2.0
    )
        : m_center(center)
        , m_diameter(diameter)
    {
    }
 
    template <class NodesDerived>
    explicit bounding_box(Eigen::DenseBase<NodesDerived> const &nodes)
        : m_diameter(0.)
    {
        for (size_t d = 0; d < dimension; ++d)
        {
            scalar_t a = nodes.row(d).minCoeff();
            scalar_t b = nodes.row(d).maxCoeff();
            m_center(d) = (a + b) / 2.;
            m_diameter = std::max(m_diameter, b-a);
        }
    }
 
    location_t const &get_center(void) const
    {
        return m_center;
    }
 
    void set_center(location_t const &c)
    {
        m_center = c;
    }
 
    scalar_t get_diameter(void) const
    {
        return m_diameter;
    }
 
    void set_diameter(scalar_t const &d)
    {
        m_diameter = d;
    }
 
    bounding_box get_child(size_t idx) const
    {
        location_t offset = idx2dist(idx) * (m_diameter / 4.0);
        return bounding_box(m_center + offset, m_diameter / 2.0);
    }
 
    bool is_adjacent(bounding_box const &other, double tol = 1e-3) const
    {
        location_t dist = other.get_center() - this->get_center();
        auto D = (this->get_diameter() + other.get_diameter()) / 2.;
        for (auto d = 0; d < dimension; ++d)
            if (std::abs(dist(d)) > D * (1.0 + tol))
                return false;
        return true;
    }
 
    static unsigned dist2idx(location_t const &child, location_t const &parent)
    {
        auto idx = 0;
        for (auto d = 0; d < dimension; ++d)
            if (child(d) > parent(d))
                idx |= (1 << d);
        return idx;
    }
 
    static location_t idx2dist(size_t idx)
    {
        if (idx >= (1 << dimension))
            throw std::invalid_argument("bounding_box child index too large");
 
        location_t dst = location_t::Zero();
        for (auto d = 0; d < dimension; ++d)
            dst(d) = ((idx >> d) & 1) == 1 ? 1.0 : -1.0;
        return dst;
    }
 
    void print_debug(std::ostream &os = std::cout) const
    {
        os << "(" << this->get_center().transpose() << ") " << this->get_diameter();
    }
 
private:
    location_t m_center;
    scalar_t m_diameter;
};
 
template <size_t Dim, class T>
std::ostream &operator<<(std::ostream &os, bounding_box<Dim, T> const &bb)
{
    bb.print_debug(os);
    return os;
}
} // end of namespace fmm
} // end of namespace NiHu
 
#endif /* NIHU_BOUNDING_BOX_HPP_INCLUDED */