AABB.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <spdlog/fmt/ostr.h>
 
#include <Eigen/Core>
#include <bitset>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdlib>
#include <iterator>
#include <limits>
#include <tuple>
#include <vector>
 
#include "BaseLib/Error.h"
 
namespace GeoLib
{
struct MinMaxPoints final
{
    Eigen::Vector3d min;
    Eigen::Vector3d max;
};

class AABB
{
public:
    template <typename PNT_TYPE>
    AABB(std::vector<PNT_TYPE*> const& pnts,
         std::vector<std::size_t> const& ids)
    {
        assert(!ids.empty());
        init(pnts[ids[0]]);
        for (std::size_t i = 1; i < ids.size(); ++i)
        {
            updateWithoutEnlarge(*(pnts[ids[i]]));
        }
        enlarge();
    }

    template <typename InputIterator>
    AABB(InputIterator first, InputIterator last)
    {
        if (std::distance(first, last) <= 0)
        {
            OGS_FATAL(
                "AABB::AABB(InputIterator first, InputIterator last): first > "
                "last");
        }
        init(*first);
        InputIterator it(first);
        while (it != last)
        {
            updateWithoutEnlarge(*it);
            it++;
        }
        enlarge();
    }

    template <typename PNT_TYPE>
    bool update(PNT_TYPE const& p)
    {
        // First component of the pair signals if the minimum point is changed
        // Second component signals not only if the max point is changed.
        // Furthermore it is signaled what coordinate (0,1,2) is changed.
        std::pair<bool, std::bitset<3>> updated(false, 0);
        for (std::size_t k(0); k < 3; k++)
        {
            // if the minimum point is updated pair.first==true
            if (p[k] < _min_pnt[k])
            {
                _min_pnt[k] = p[k];
                updated.first = true;
            }
            // if the kth coordinate of the maximum point is updated
            // pair.second[k]==true
            if (p[k] >= _max_pnt[k])
            {
                _max_pnt[k] = p[k];
                updated.second[k] = true;
            }
        }
 
        if (updated.second.any())
        {
            enlarge(updated.second);
            return true;
        }
        return updated.first;
    }

    template <typename T>
    bool containsPoint(T const& pnt, double eps) const
    {
        if (pnt[0] < _min_pnt[0] - eps || _max_pnt[0] + eps <= pnt[0])
        {
            return false;
        }
        if (pnt[1] < _min_pnt[1] - eps || _max_pnt[1] + eps <= pnt[1])
        {
            return false;
        }
        if (pnt[2] < _min_pnt[2] - eps || _max_pnt[2] + eps <= pnt[2])
        {
            return false;
        }
        return true;
    }
 
    template <typename T>
    bool containsPointXY(T const& pnt) const
    {
        if (pnt[0] < _min_pnt[0] || _max_pnt[0] <= pnt[0])
        {
            return false;
        }
        if (pnt[1] < _min_pnt[1] || _max_pnt[1] <= pnt[1])
        {
            return false;
        }
        return true;
    }
 
    MinMaxPoints getMinMaxPoints() const { return {_min_pnt, _max_pnt}; }
    Eigen::Vector3d const& getMinPoint() const { return _min_pnt; }

    Eigen::Vector3d const& getMaxPoint() const { return _max_pnt; }

    bool containsAABB(AABB const& other_aabb) const
    {
        return containsPoint(other_aabb.getMinPoint(), 0) &&
               containsPoint(other_aabb.getMaxPoint(), 0);
    }
 
    friend std::ostream& operator<<(std::ostream& os, AABB const& aabb)
    {
        auto const [min, max] = aabb.getMinMaxPoints();
        os << "\tx [" << min[0] << ", " << max[0] << ") (extent "
           << max[0] - min[0] << ")\n";
        os << "\ty [" << min[1] << ", " << max[1] << ") (extent "
           << max[1] - min[1] << ")\n";
        os << "\tz [" << min[2] << ", " << max[2] << ") (extent "
           << max[2] - min[2] << ")";
        return os;
    }
 
private:
    Eigen::Vector3d _min_pnt{std::numeric_limits<double>::max(),
                             std::numeric_limits<double>::max(),
                             std::numeric_limits<double>::max()};
    Eigen::Vector3d _max_pnt{std::numeric_limits<double>::lowest(),
                             std::numeric_limits<double>::lowest(),
                             std::numeric_limits<double>::lowest()};
 
private:
    void enlarge(std::bitset<3> to_update = 7)
    {
        for (std::size_t k = 0; k < 3; ++k)
        {
            if (to_update[k])
            {
                _max_pnt[k] = std::nextafter(
                    _max_pnt[k], std::numeric_limits<double>::max());
            }
        }
    }
 
    template <typename PNT_TYPE>
    void init(PNT_TYPE const& pnt)
    {
        _min_pnt[0] = _max_pnt[0] = pnt[0];
        _min_pnt[1] = _max_pnt[1] = pnt[1];
        _min_pnt[2] = _max_pnt[2] = pnt[2];
    }
 
    template <typename PNT_TYPE>
    void init(PNT_TYPE* const& pnt)
    {
        init(*pnt);
    }

    template <typename PNT_TYPE>
    void updateWithoutEnlarge(PNT_TYPE const& p)
    {
        for (std::size_t k(0); k < 3; k++)
        {
            if (p[k] < _min_pnt[k])
            {
                _min_pnt[k] = p[k];
            }
            if (p[k] >= _max_pnt[k])
            {
                _max_pnt[k] = p[k];
            }
        }
    }
 
    template <typename PNT_TYPE>
    void updateWithoutEnlarge(PNT_TYPE* const& pnt)
    {
        updateWithoutEnlarge(*pnt);
    }
 
    template <typename PNT_TYPE>
    void update(PNT_TYPE const* pnt)
    {
        update(*pnt);
    }
};
}  // namespace GeoLib
 
namespace fmt
{
template <>
struct formatter<::GeoLib::AABB> : ostream_formatter
{
};
}  // namespace fmt