Revision 652ad60247b0023bd9837db6a29f8d47413bfee1 authored by Christoph Lehmann on 24 August 2023, 07:23:25 UTC, committed by Christoph Lehmann on 24 August 2023, 07:23:25 UTC
Draft: RTP acceleration: assemble only once See merge request ogs/ogs!4713
AABB.h
/**
* \file
* \author Thomas Fischer
* \date 2010-04-22
* \brief Definition of the AABB class.
*
* \copyright
* Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#pragma once
#include <fmt/ostream.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;
};
/**
* \brief Class AABB is an axis aligned bounding box around a given
* set of geometric points of (template) type PNT_TYPE.
*
* Let \f$P = \{p_k \in \mathbb{R}^3, \ k=1, \dotsc, n\}\f$ a set of 3d points.
* The bounding volume is described by its lower, left, front point \f$\ell\f$
* and the upper, right, back point \f$u\f$, i.e. the coordinates of \f$\ell\f$
* and \f$u\f$ are computed as follows \f$\ell_i = \min \limits_{p \in P}
* p_i\f$, \f$u_i = \max \limits_{p \in P} p_i\f$, respectively. The bounding
* box consists of half-open intervals \f$[\ell_x, u_x) \times [\ell_y, u_y)
* \times [\ell_z, u_z).\f$ The bounding box is enlarged up to the next
* available floating point number such that all input points are contained in
* the bounding box.
*
*/
class AABB
{
public:
/**
* construction of object, initialization the axis aligned bounding box
* @tparam PNT_TYPE a point type supporting accessing the coordinates via
* operator[]
* */
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();
}
/**
* Construction of object using input iterators. In contrast to give a
* vector this approach is more generic. You can use every (stl) container
* and C arrays as input for constructing the object.
* @attention{The constructor requires that std::distance(first, last) > 0.}
* @param first the input iterator to the initial position in the sequence
* @param last the input iterator to the final position in a container, i.e.
* [first, last).
* @attention{The iterator last must be reachable from first.}
*/
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();
}
/// Checks if the bounding box has to be updated.
/// @return true if AABB is updated.
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;
}
/**
* check if point is in the axis aligned bounding box
*/
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}; }
/**
* returns a point that coordinates are minimal for each dimension
* for the given point set
* @return a point
*/
Eigen::Vector3d const& getMinPoint() const { return _min_pnt; }
/**
* returns a point that coordinates are maximal for each dimension
* within the given point set
* @return a point
*/
Eigen::Vector3d const& getMaxPoint() const { return _max_pnt; }
/**
* Method checks if the given AABB object is contained within the
* AABB represented by this object.
* @param other_aabb the AABB to test with
* @return true if the other AABB is contained in the AABB
* represented by this object
*/
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:
/// Enlarge the bounding box the smallest possible amount (modifying the
/// unit in the last place). Only the coordinates of the maximum point are
/// changed such that the half-open property will be preserved.
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);
}
/// Private method that is used internally to update the min and max point
/// of the bounding box using point \f$p\f$ without enlarging the bounding
/// box. Using this method the bounding box of the initial point set is
/// enlarged only once.
/// @param p point that will possibly change the bounding box points
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
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...