swh:1:snp:f521c49ab17ef7db6ec70b2430e1ed203f50383f
Tip revision: e655f492e7ec9806d2284bbb2343807b772bfec4 authored by rinkk on 26 May 2021, 12:26:20 UTC
[utils] adding arguments to create (incomplete) second order elements
[utils] adding arguments to create (incomplete) second order elements
Tip revision: e655f49
OctTree-impl.h
/**
* \date 2015-06-12
* \brief Implementation of the OctTree class.
*
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/
namespace GeoLib {
template <typename POINT, std::size_t MAX_POINTS>
template <typename T>
OctTree<POINT, MAX_POINTS>* OctTree<POINT, MAX_POINTS>::createOctTree(T ll, T ur,
double eps)
{
// compute an axis aligned cube around the points ll and ur
const double dx(ur[0] - ll[0]);
const double dy(ur[1] - ll[1]);
const double dz(ur[2] - ll[2]);
if (dx >= dy && dx >= dz) {
ll[1] -= (dx-dy)/2.0;
ur[1] += (dx-dy)/2.0;
ll[2] -= (dx-dz)/2.0;
ur[2] += (dx-dz)/2.0;
} else {
if (dy >= dx && dy >= dz) {
ll[0] -= (dy-dx)/2.0;
ur[0] += (dy-dx)/2.0;
ll[2] -= (dy-dz)/2.0;
ur[2] += (dy-dz)/2.0;
} else {
ll[0] -= (dz-dx)/2.0;
ur[0] += (dz-dx)/2.0;
ll[1] -= (dz-dy)/2.0;
ur[1] += (dz-dy)/2.0;
}
}
if (eps == 0.0)
{
eps = std::numeric_limits<double>::epsilon();
}
for (std::size_t k(0); k<3; ++k) {
if (ur[k] - ll[k] > 0.0) {
ur[k] += (ur[k] - ll[k]) * 1e-6;
} else {
ur[k] += eps;
}
}
return new OctTree<POINT, MAX_POINTS>(ll, ur, eps);
}
template <typename POINT, std::size_t MAX_POINTS>
OctTree<POINT, MAX_POINTS>::~OctTree()
{
for (auto c : _children)
{
delete c;
}
}
template <typename POINT, std::size_t MAX_POINTS>
bool OctTree<POINT, MAX_POINTS>::addPoint(POINT * pnt, POINT *& ret_pnt)
{
// first do a range query using a epsilon box around the point pnt
std::vector<POINT*> query_pnts;
MathLib::Point3d min(
std::array<double,3>{{(*pnt)[0]-_eps, (*pnt)[1]-_eps, (*pnt)[2]-_eps}});
MathLib::Point3d max(
std::array<double,3>{{(*pnt)[0]+_eps, (*pnt)[1]+_eps, (*pnt)[2]+_eps}});
getPointsInRange(min, max, query_pnts);
auto const it = std::find_if(
query_pnts.begin(), query_pnts.end(), [pnt, this](auto const* p) {
return MathLib::sqrDist(*p, *pnt) < _eps * _eps;
});
if (it != query_pnts.end())
{
ret_pnt = *it;
return false;
}
// the point pnt is not yet in the OctTree
if (isOutside(pnt)) {
ret_pnt = nullptr;
return false;
}
// at this place it holds true that the point is within [_ll, _ur]
if (!_is_leaf) {
for (auto c : _children) {
if (c->addPoint_(pnt, ret_pnt)) {
return true;
}
if (ret_pnt != nullptr)
{
return false;
}
}
}
ret_pnt = pnt;
if (_pnts.size () < MAX_POINTS) {
_pnts.push_back(pnt);
} else { // i.e. _pnts.size () == MAX_POINTS
splitNode(pnt);
_pnts.clear();
}
return true;
}
template <typename POINT, std::size_t MAX_POINTS>
template <typename T>
void OctTree<POINT, MAX_POINTS>::getPointsInRange(T const& min, T const& max,
std::vector<POINT*> &pnts) const
{
if (_ur[0] < min[0] || _ur[1] < min[1] || _ur[2] < min[2])
{
return;
}
if (max[0] < _ll[0] || max[1] < _ll[1] || max[2] < _ll[2])
{
return;
}
if (_is_leaf) {
std::copy_if(_pnts.begin(), _pnts.end(), std::back_inserter(pnts),
[&min, &max](auto const* p) {
return (min[0] <= (*p)[0] && (*p)[0] < max[0] &&
min[1] <= (*p)[1] && (*p)[1] < max[1] &&
min[2] <= (*p)[2] && (*p)[2] < max[2]);
});
} else {
for (std::size_t k(0); k<8; k++) {
_children[k]->getPointsInRange(min, max, pnts);
}
}
}
template <typename POINT, std::size_t MAX_POINTS>
OctTree<POINT, MAX_POINTS>::OctTree(
MathLib::Point3d const& ll, MathLib::Point3d const& ur, double eps)
: _ll(ll), _ur(ur), _is_leaf(true), _eps(eps)
{
_children.fill(nullptr);
}
template <typename POINT, std::size_t MAX_POINTS>
bool OctTree<POINT, MAX_POINTS>::addPoint_(POINT * pnt, POINT *& ret_pnt)
{
if (isOutside(pnt)) {
ret_pnt = nullptr;
return false;
}
// at this place it holds true that the point is within [_ll, _ur]
if (!_is_leaf) {
for (auto c : _children) {
if (c->addPoint_(pnt, ret_pnt)) {
return true;
}
if (ret_pnt != nullptr)
{
return false;
}
}
}
ret_pnt = pnt;
if (_pnts.size() < MAX_POINTS) {
_pnts.push_back(pnt);
} else { // i.e. _pnts.size () == MAX_POINTS
splitNode(pnt);
_pnts.clear();
}
return true;
}
template <typename POINT, std::size_t MAX_POINTS>
bool OctTree<POINT, MAX_POINTS>::addPointToChild(POINT * pnt)
{
if (isOutside(pnt))
{
return false;
}
if (_pnts.size() < MAX_POINTS) {
_pnts.push_back(pnt);
} else { // i.e. _pnts.size () == MAX_POINTS
splitNode(pnt);
_pnts.clear();
}
return true;
}
template <typename POINT, std::size_t MAX_POINTS>
void OctTree<POINT, MAX_POINTS>::splitNode(POINT * pnt)
{
const double x_mid((_ur[0] + _ll[0]) / 2.0);
const double y_mid((_ur[1] + _ll[1]) / 2.0);
const double z_mid((_ur[2] + _ll[2]) / 2.0);
MathLib::Point3d p0(std::array<double,3>{{x_mid, y_mid, _ll[2]}});
MathLib::Point3d p1(std::array<double,3>{{_ur[0], _ur[1], z_mid}});
// create child NEL
_children[static_cast<std::int8_t>(Quadrant::NEL)]
= new OctTree<POINT, MAX_POINTS> (p0, p1, _eps);
// create child NWL
p0[0] = _ll[0];
p1[0] = x_mid;
_children[static_cast<std::int8_t>(Quadrant::NWL)]
= new OctTree<POINT, MAX_POINTS> (p0, p1, _eps);
// create child SWL
p0[1] = _ll[1];
p1[1] = y_mid;
_children[static_cast<std::int8_t>(Quadrant::SWL)]
= new OctTree<POINT, MAX_POINTS> (_ll, p1, _eps);
// create child NEU
_children[static_cast<std::int8_t>(Quadrant::NEU)]
= new OctTree<POINT, MAX_POINTS> (p1, _ur, _eps);
// create child SEL
p0[0] = x_mid;
p1[0] = _ur[0];
_children[static_cast<std::int8_t>(Quadrant::SEL)]
= new OctTree<POINT, MAX_POINTS> (p0, p1, _eps);
// create child NWU
p0[0] = _ll[0];
p0[1] = y_mid;
p0[2] = z_mid;
p1[0] = x_mid;
p1[1] = _ur[1];
p1[2] = _ur[2];
_children[static_cast<std::int8_t>(Quadrant::NWU)]
= new OctTree<POINT, MAX_POINTS> (p0, p1, _eps);
// create child SWU
p0[1] = _ll[1];
p1[1] = y_mid;
_children[static_cast<std::int8_t>(Quadrant::SWU)]
= new OctTree<POINT, MAX_POINTS> (p0, p1, _eps);
// create child SEU
p0[0] = x_mid;
p1[0] = _ur[0];
p1[1] = y_mid;
p1[2] = _ur[2];
_children[static_cast<std::int8_t>(Quadrant::SEU)]
= new OctTree<POINT, MAX_POINTS> (p0, p1, _eps);
// add the passed point pnt to the children at first
for (std::size_t k(0); k < 8; k++) {
if (_children[k]->addPointToChild(pnt))
{
break;
}
}
// distribute points to sub quadtrees
const std::size_t n_pnts(_pnts.size());
for (std::size_t j(0); j < n_pnts; j++) {
for (auto c : _children) {
if (c->addPointToChild(_pnts[j])) {
break;
}
}
}
_is_leaf = false;
}
template <typename POINT, std::size_t MAX_POINTS>
bool OctTree<POINT, MAX_POINTS>::isOutside(POINT * pnt) const
{
if ((*pnt)[0] < _ll[0] || (*pnt)[1] < _ll[1] || (*pnt)[2] < _ll[2])
{
return true;
}
if ((*pnt)[0] >= _ur[0] || (*pnt)[1] >= _ur[1] || (*pnt)[2] >= _ur[2])
{
return true;
}
return false;
}
} // end namespace GeoLib