Commit 8576730c authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov
Browse files

Merge branch 'SmallFixes_2021-06' into 'master'

Small unrelated fixes almost all in GeoLib

See merge request ogs/ogs!3661
parents 6f17c15a eb59a31f
......@@ -25,8 +25,8 @@
#include <vector>
#include "BaseLib/Error.h"
#include "MathLib/Point3d.h"
#include "MathLib/MathTools.h"
#include "MathLib/Point3d.h"
namespace GeoLib
{
......@@ -54,11 +54,13 @@ public:
* operator[]
* */
template <typename PNT_TYPE>
AABB(std::vector<PNT_TYPE*> const& pnts, std::vector<std::size_t> const& ids)
AABB(std::vector<PNT_TYPE*> const& pnts,
std::vector<std::size_t> const& ids)
{
assert(! ids.empty());
assert(!ids.empty());
init(pnts[ids[0]]);
for (std::size_t i=1; i<ids.size(); ++i) {
for (std::size_t i = 1; i < ids.size(); ++i)
{
updateWithoutEnlarge(*(pnts[ids[i]]));
}
enlarge();
......@@ -67,24 +69,28 @@ public:
AABB(AABB const&) = default;
/**
* 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.
* 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).
* @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)
if (std::distance(first, last) <= 0)
{
OGS_FATAL("AABB::AABB(InputIterator first, InputIterator last): first > last");
OGS_FATAL(
"AABB::AABB(InputIterator first, InputIterator last): first > "
"last");
}
init(*first);
InputIterator it(first);
while (it != last) {
while (it != last)
{
updateWithoutEnlarge(*it);
it++;
}
......@@ -94,27 +100,31 @@ public:
/// 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)
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++) {
for (std::size_t k(0); k < 3; k++)
{
// if the minimum point is updated pair.first==true
if (p[k] < _min_pnt[k]) {
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]) {
if (p[k] >= _max_pnt[k])
{
_max_pnt[k] = p[k];
updated.second[k] = true;
}
}
if (updated.second.any()) {
if (updated.second.any())
{
enlarge(updated.second);
return true;
}
......@@ -125,7 +135,7 @@ public:
* check if point is in the axis aligned bounding box
*/
template <typename T>
bool containsPoint(T const & pnt, double eps) const
bool containsPoint(T const& pnt, double eps) const
{
if (pnt[0] < _min_pnt[0] - eps || _max_pnt[0] + eps <= pnt[0])
{
......@@ -143,7 +153,7 @@ public:
}
template <typename T>
bool containsPointXY(T const & pnt) const
bool containsPointXY(T const& pnt) const
{
if (pnt[0] < _min_pnt[0] || _max_pnt[0] <= pnt[0])
{
......@@ -184,30 +194,32 @@ public:
}
protected:
MathLib::Point3d _min_pnt = MathLib::Point3d{std::array<double,3>{{
std::numeric_limits<double>::max(),
std::numeric_limits<double>::max(),
std::numeric_limits<double>::max()}}};
MathLib::Point3d _max_pnt = MathLib::Point3d{std::array<double,3>{{
std::numeric_limits<double>::lowest(),
std::numeric_limits<double>::lowest(),
std::numeric_limits<double>::lowest()}}};
MathLib::Point3d _min_pnt = MathLib::Point3d{std::array<double, 3>{
{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(),
std::numeric_limits<double>::max()}}};
MathLib::Point3d _max_pnt = MathLib::Point3d{
std::array<double, 3>{{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());
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)
void init(PNT_TYPE const& pnt)
{
_min_pnt[0] = _max_pnt[0] = pnt[0];
_min_pnt[1] = _max_pnt[1] = pnt[1];
......@@ -215,7 +227,7 @@ private:
}
template <typename PNT_TYPE>
void init(PNT_TYPE * const & pnt)
void init(PNT_TYPE* const& pnt)
{
init(*pnt);
}
......@@ -226,26 +238,29 @@ private:
/// enlarged only once.
/// @param p point that will possibly change the bounding box points
template <typename PNT_TYPE>
void updateWithoutEnlarge(PNT_TYPE const & p)
void updateWithoutEnlarge(PNT_TYPE const& p)
{
for (std::size_t k(0); k<3; k++) {
if (p[k] < _min_pnt[k]) {
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]) {
if (p[k] >= _max_pnt[k])
{
_max_pnt[k] = p[k];
}
}
}
template <typename PNT_TYPE>
void updateWithoutEnlarge(PNT_TYPE * const & pnt)
void updateWithoutEnlarge(PNT_TYPE* const& pnt)
{
updateWithoutEnlarge(*pnt);
}
template <typename PNT_TYPE>
void update(PNT_TYPE const * pnt)
void update(PNT_TYPE const* pnt)
{
update(*pnt);
}
......
This diff is collapsed.
......@@ -194,6 +194,9 @@ void BoostXmlGmlInterface::readPolylines(
{
// polyline has no name, ignore it.
pl.ignoreConfigParameterAll("pnt");
WARN(
"Polyline name is required! Polylines without a name are "
"ignored.");
}
}
}
......
......@@ -10,43 +10,53 @@
* http://www.opengeosys.org/project/license
*/
namespace GeoLib {
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)
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 (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) {
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 {
}
else
{
ur[k] += eps;
}
}
......@@ -63,19 +73,19 @@ OctTree<POINT, MAX_POINTS>::~OctTree()
}
template <typename POINT, std::size_t MAX_POINTS>
bool OctTree<POINT, MAX_POINTS>::addPoint(POINT * pnt, POINT *& ret_pnt)
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}});
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;
});
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;
......@@ -83,15 +93,19 @@ bool OctTree<POINT, MAX_POINTS>::addPoint(POINT * pnt, POINT *& ret_pnt)
}
// the point pnt is not yet in the OctTree
if (isOutside(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)) {
if (!_is_leaf)
{
for (auto c : _children)
{
if (c->addPoint_(pnt, ret_pnt))
{
return true;
}
if (ret_pnt != nullptr)
......@@ -103,9 +117,12 @@ bool OctTree<POINT, MAX_POINTS>::addPoint(POINT * pnt, POINT *& ret_pnt)
ret_pnt = pnt;
if (_pnts.size () < MAX_POINTS) {
if (_pnts.size() < MAX_POINTS)
{
_pnts.push_back(pnt);
} else { // i.e. _pnts.size () == MAX_POINTS
}
else
{ // i.e. _pnts.size () == MAX_POINTS
splitNode(pnt);
_pnts.clear();
}
......@@ -114,8 +131,8 @@ bool OctTree<POINT, MAX_POINTS>::addPoint(POINT * pnt, POINT *& ret_pnt)
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
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])
{
......@@ -127,40 +144,49 @@ void OctTree<POINT, MAX_POINTS>::getPointsInRange(T const& min, T const& max,
return;
}
if (_is_leaf) {
if (_is_leaf)
{
std::copy_if(_pnts.begin(), _pnts.end(), std::back_inserter(pnts),
[&min, &max](auto const* p) {
[&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++) {
}
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)
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)
bool OctTree<POINT, MAX_POINTS>::addPoint_(POINT* pnt, POINT*& ret_pnt)
{
if (isOutside(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)) {
if (!_is_leaf)
{
for (auto c : _children)
{
if (c->addPoint_(pnt, ret_pnt))
{
return true;
}
if (ret_pnt != nullptr)
......@@ -171,9 +197,12 @@ bool OctTree<POINT, MAX_POINTS>::addPoint_(POINT * pnt, POINT *& ret_pnt)
}
ret_pnt = pnt;
if (_pnts.size() < MAX_POINTS) {
if (_pnts.size() < MAX_POINTS)
{
_pnts.push_back(pnt);
} else { // i.e. _pnts.size () == MAX_POINTS
}
else
{ // i.e. _pnts.size () == MAX_POINTS
splitNode(pnt);
_pnts.clear();
}
......@@ -181,16 +210,19 @@ bool OctTree<POINT, MAX_POINTS>::addPoint_(POINT * pnt, POINT *& ret_pnt)
}
template <typename POINT, std::size_t MAX_POINTS>
bool OctTree<POINT, MAX_POINTS>::addPointToChild(POINT * pnt)
bool OctTree<POINT, MAX_POINTS>::addPointToChild(POINT* pnt)
{
if (isOutside(pnt))
{
return false;
}
if (_pnts.size() < MAX_POINTS) {
if (_pnts.size() < MAX_POINTS)
{
_pnts.push_back(pnt);
} else { // i.e. _pnts.size () == MAX_POINTS
}
else
{ // i.e. _pnts.size () == MAX_POINTS
splitNode(pnt);
_pnts.clear();
}
......@@ -198,39 +230,39 @@ bool OctTree<POINT, MAX_POINTS>::addPointToChild(POINT * pnt)
}
template <typename POINT, std::size_t MAX_POINTS>
void OctTree<POINT, MAX_POINTS>::splitNode(POINT * pnt)
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}});
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);
_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);
_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);
_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);
_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);
_children[static_cast<std::int8_t>(Quadrant::SEL)] =
new OctTree<POINT, MAX_POINTS>(p0, p1, _eps);
// create child NWU
p0[0] = _ll[0];
......@@ -239,25 +271,26 @@ void OctTree<POINT, MAX_POINTS>::splitNode(POINT * pnt)
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);
_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);
_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);
_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++) {
for (std::size_t k(0); k < 8; k++)
{
if (_children[k]->addPointToChild(pnt))
{
break;
......@@ -266,9 +299,12 @@ void OctTree<POINT, MAX_POINTS>::splitNode(POINT * pnt)
// 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])) {
for (std::size_t j(0); j < n_pnts; j++)
{
for (auto c : _children)
{
if (c->addPointToChild(_pnts[j]))
{
break;
}
}
......@@ -277,7 +313,7 @@ void OctTree<POINT, MAX_POINTS>::splitNode(POINT * pnt)
}
template <typename POINT, std::size_t MAX_POINTS>
bool OctTree<POINT, MAX_POINTS>::isOutside(POINT * pnt) const
bool OctTree<POINT, MAX_POINTS>::isOutside(POINT* pnt) const
{
if ((*pnt)[0] < _ll[0] || (*pnt)[1] < _ll[1] || (*pnt)[2] < _ll[2])
{
......@@ -289,5 +325,4 @@ bool OctTree<POINT, MAX_POINTS>::isOutside(POINT * pnt) const
}
return false;
}
} // end namespace GeoLib
} // end namespace GeoLib
......@@ -15,18 +15,19 @@
#pragma once
#include <cstdint>
#include <vector>
#include <limits>
#include <vector>
#include "MathLib/MathTools.h"
#include "MathLib/Point3d.h"
namespace GeoLib {
namespace GeoLib