From a2f719b89c2aad2a8f262cceae46dd9b4b6aa0e6 Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Mon, 14 Jun 2021 10:10:04 +0200
Subject: [PATCH] [GL/OctTree] Use Eigen::Vector3d instead of MathLib::Point3d.
---
GeoLib/OctTree-impl.h | 35 +++++++++++++++++--------------
GeoLib/OctTree.h | 18 ++++++++--------
GeoLib/PointVec.cpp | 17 ++++++++++++---
Tests/GeoLib/TestOctTree.cpp | 40 +++++++++++++++++++++++-------------
4 files changed, 68 insertions(+), 42 deletions(-)
diff --git a/GeoLib/OctTree-impl.h b/GeoLib/OctTree-impl.h
index e044e7609e6..45ea89d9feb 100644
--- a/GeoLib/OctTree-impl.h
+++ b/GeoLib/OctTree-impl.h
@@ -13,9 +13,8 @@
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)
+ Eigen::Vector3d ll, Eigen::Vector3d ur, double eps)
{
// compute an axis aligned cube around the points ll and ur
const double dx(ur[0] - ll[0]);
@@ -60,7 +59,9 @@ OctTree<POINT, MAX_POINTS>* OctTree<POINT, MAX_POINTS>::createOctTree(
ur[k] += eps;
}
}
- return new OctTree<POINT, MAX_POINTS>(ll, ur, eps);
+ Eigen::Vector3d lower_left{ll[0], ll[1], ll[2]};
+ Eigen::Vector3d upper_right{ur[0], ur[1], ur[2]};
+ return new OctTree<POINT, MAX_POINTS>(lower_left, upper_right, eps);
}
template <typename POINT, std::size_t MAX_POINTS>
@@ -77,15 +78,19 @@ 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}});
+ Eigen::Vector3d const min =
+ Eigen::Map<Eigen::Vector3d>(pnt->getCoords()).array() - _eps;
+ Eigen::Vector3d const max =
+ Eigen::Map<Eigen::Vector3d>(pnt->getCoords()).array() + _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 (Eigen::Map<Eigen::Vector3d const>(p->getCoords()) -
+ Eigen::Map<Eigen::Vector3d const>(pnt->getCoords()))
+ .squaredNorm() < _eps * _eps;
+ });
if (it != query_pnts.end())
{
ret_pnt = *it;
@@ -164,8 +169,8 @@ void OctTree<POINT, MAX_POINTS>::getPointsInRange(
}
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(Eigen::Vector3d const& ll,
+ Eigen::Vector3d const& ur, double eps)
: _ll(ll), _ur(ur), _is_leaf(true), _eps(eps)
{
_children.fill(nullptr);
@@ -235,8 +240,8 @@ 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}});
+ Eigen::Vector3d p0{x_mid, y_mid, _ll[2]};
+ Eigen::Vector3d p1{_ur[0], _ur[1], z_mid};
// create child NEL
_children[static_cast<std::int8_t>(Quadrant::NEL)] =
diff --git a/GeoLib/OctTree.h b/GeoLib/OctTree.h
index 3b1d17ff42a..b780945234e 100644
--- a/GeoLib/OctTree.h
+++ b/GeoLib/OctTree.h
@@ -14,13 +14,11 @@
#pragma once
+#include <Eigen/Eigen>
#include <cstdint>
#include <limits>
#include <vector>
-#include "MathLib/MathTools.h"
-#include "MathLib/Point3d.h"
-
namespace GeoLib
{
/// @tparam POINT point data type the OctTree will use
@@ -45,9 +43,9 @@ public:
/// the memory requirements for the OctTree may be lower but the search
/// inside a OctTree leaf may be more expensive. The value should be
/// chosen application dependent. [default 8]
- template <typename T>
static OctTree<POINT, MAX_POINTS>* createOctTree(
- T ll, T ur, double eps = std::numeric_limits<double>::epsilon());
+ Eigen::Vector3d ll, Eigen::Vector3d ur,
+ double eps = std::numeric_limits<double>::epsilon());
/// Destroys the children of this node. @attention Does not destroy the
/// pointers to the managed objects.
@@ -76,8 +74,8 @@ public:
std::vector<POINT*>& pnts) const;
#ifndef NDEBUG
- MathLib::Point3d const& getLowerLeftCornerPoint() const { return _ll; }
- MathLib::Point3d const& getUpperRightCornerPoint() const { return _ur; }
+ Eigen::Vector3d const& getLowerLeftCornerPoint() const { return _ll; }
+ Eigen::Vector3d const& getUpperRightCornerPoint() const { return _ur; }
OctTree<POINT, MAX_POINTS> const* getChild(std::size_t i) const
{
return _children[i];
@@ -90,7 +88,7 @@ private:
/// @param ll lower left point
/// @param ur upper right point
/// @param eps the euclidean distance as a threshold to make objects unique
- OctTree(MathLib::Point3d const& ll, MathLib::Point3d const& ur, double eps);
+ OctTree(Eigen::Vector3d const& ll, Eigen::Vector3d const& ur, double eps);
enum class Quadrant : std::int8_t
{
@@ -140,9 +138,9 @@ private:
/// _children[7] is south east upper child
std::array<OctTree<POINT, MAX_POINTS>*, 8> _children;
/// lower left front face point of the cube
- MathLib::Point3d const _ll;
+ Eigen::Vector3d const _ll;
/// upper right back face point of the cube
- MathLib::Point3d const _ur;
+ Eigen::Vector3d const _ur;
/// vector of pointers to POINT objects
std::vector<POINT*> _pnts;
/// flag if this OctTree is a leaf
diff --git a/GeoLib/PointVec.cpp b/GeoLib/PointVec.cpp
index 6822171f434..5fb55a23751 100644
--- a/GeoLib/PointVec.cpp
+++ b/GeoLib/PointVec.cpp
@@ -21,6 +21,14 @@
namespace GeoLib
{
+namespace details
+{
+Eigen::Vector3d convertToEigen(MathLib::Point3d point3d)
+{
+ return Eigen::Vector3d{point3d[0], point3d[1], point3d[2]};
+}
+} // namespace details
+
PointVec::PointVec(
const std::string& name, std::unique_ptr<std::vector<Point*>> points,
std::unique_ptr<std::map<std::string, std::size_t>> name_id_map,
@@ -31,7 +39,8 @@ PointVec::PointVec(
_rel_eps(rel_eps * std::sqrt(MathLib::sqrDist(_aabb.getMinPoint(),
_aabb.getMaxPoint()))),
_oct_tree(OctTree<GeoLib::Point, 16>::createOctTree(
- _aabb.getMinPoint(), _aabb.getMaxPoint(), _rel_eps))
+ details::convertToEigen(_aabb.getMinPoint()),
+ details::convertToEigen(_aabb.getMaxPoint()), _rel_eps))
{
assert(_data_vec);
std::size_t const number_of_all_input_pnts(_data_vec->size());
@@ -180,7 +189,8 @@ std::size_t PointVec::uniqueInsert(Point* pnt)
_aabb.update(*pnt);
// recreate the (enlarged) OctTree
_oct_tree.reset(GeoLib::OctTree<GeoLib::Point, 16>::createOctTree(
- _aabb.getMinPoint(), _aabb.getMaxPoint(), _rel_eps));
+ details::convertToEigen(_aabb.getMinPoint()),
+ details::convertToEigen(_aabb.getMaxPoint()), _rel_eps));
// add all points that are already in the _data_vec
for (std::size_t k(0); k < _data_vec->size(); ++k)
{
@@ -267,7 +277,8 @@ void PointVec::resetInternalDataStructures()
_aabb.getMaxPoint()));
_oct_tree.reset(OctTree<GeoLib::Point, 16>::createOctTree(
- _aabb.getMinPoint(), _aabb.getMaxPoint(), _rel_eps));
+ details::convertToEigen(_aabb.getMinPoint()),
+ details::convertToEigen(_aabb.getMaxPoint()), _rel_eps));
GeoLib::Point* ret_pnt(nullptr);
for (auto const p : *_data_vec)
diff --git a/Tests/GeoLib/TestOctTree.cpp b/Tests/GeoLib/TestOctTree.cpp
index 62cb5930118..acd99a6a749 100644
--- a/Tests/GeoLib/TestOctTree.cpp
+++ b/Tests/GeoLib/TestOctTree.cpp
@@ -17,6 +17,14 @@
#include "GeoLib/OctTree.h"
#include "GeoLib/Point.h"
+namespace
+{
+Eigen::Vector3d convertToEigen(MathLib::Point3d point3d)
+{
+ return Eigen::Vector3d{point3d[0], point3d[1], point3d[2]};
+}
+} // namespace
+
class GeoLibOctTree : public testing::Test
{
public:
@@ -89,12 +97,13 @@ TEST_F(GeoLibOctTree, TestWithEquidistantPoints3d)
generateEquidistantPoints3d();
double const eps(10 * std::numeric_limits<double>::epsilon());
std::unique_ptr<GeoLib::OctTree<GeoLib::Point, 2>> oct_tree(
- GeoLib::OctTree<GeoLib::Point, 2>::createOctTree(*ps_ptr.front(),
- *ps_ptr.back(), eps));
+ GeoLib::OctTree<GeoLib::Point, 2>::createOctTree(
+ convertToEigen(*ps_ptr.front()), convertToEigen(*ps_ptr.back()),
+ eps));
#ifndef NDEBUG
- MathLib::Point3d const& ll(oct_tree->getLowerLeftCornerPoint());
- MathLib::Point3d const& ur(oct_tree->getUpperRightCornerPoint());
+ Eigen::Vector3d const& ll(oct_tree->getLowerLeftCornerPoint());
+ Eigen::Vector3d const& ur(oct_tree->getUpperRightCornerPoint());
EXPECT_EQ((*ps_ptr.front())[0], ll[0]);
EXPECT_EQ((*ps_ptr.front())[1], ll[1]);
@@ -277,8 +286,8 @@ TEST_F(GeoLibOctTree, TestWithAlternatingPoints3d)
ps_ptr.push_back(new GeoLib::Point(5 * small_displacement, 1, 0, 6));
GeoLib::AABB const aabb(ps_ptr.cbegin(), ps_ptr.cend());
- const MathLib::Point3d& min(aabb.getMinPoint());
- const MathLib::Point3d& max(aabb.getMaxPoint());
+ auto const min(convertToEigen(aabb.getMinPoint()));
+ auto const max(convertToEigen(aabb.getMaxPoint()));
std::unique_ptr<GeoLib::OctTree<GeoLib::Point, 8>> oct_tree(
GeoLib::OctTree<GeoLib::Point, 8>::createOctTree(min, max, eps));
@@ -332,8 +341,8 @@ TEST_F(GeoLibOctTree, TestSmallDistanceDifferentLeaves)
// create OctTree
GeoLib::AABB const aabb(ps_ptr.cbegin(), ps_ptr.cend());
- const MathLib::Point3d& min(aabb.getMinPoint());
- const MathLib::Point3d& max(aabb.getMaxPoint());
+ auto const min(convertToEigen(aabb.getMinPoint()));
+ auto const max(convertToEigen(aabb.getMaxPoint()));
std::unique_ptr<GeoLib::OctTree<GeoLib::Point, 2>> oct_tree(
GeoLib::OctTree<GeoLib::Point, 2>::createOctTree(min, max, eps));
@@ -364,9 +373,10 @@ TEST_F(GeoLibOctTree, TestOctTreeWithTwoEqualPoints)
double const eps(0.0);
GeoLib::AABB aabb(ps_ptr.begin(), ps_ptr.end());
+ auto const min(convertToEigen(aabb.getMinPoint()));
+ auto const max(convertToEigen(aabb.getMaxPoint()));
std::unique_ptr<GeoLib::OctTree<GeoLib::Point, 2>> oct_tree(
- GeoLib::OctTree<GeoLib::Point, 2>::createOctTree(
- aabb.getMinPoint(), aabb.getMaxPoint(), eps));
+ GeoLib::OctTree<GeoLib::Point, 2>::createOctTree(min, max, eps));
GeoLib::Point* ret_pnt(nullptr);
ASSERT_TRUE(oct_tree->addPoint(ps_ptr[0], ret_pnt));
@@ -382,9 +392,10 @@ TEST_F(GeoLibOctTree, TestOctTreeWithTwoEqualPointsOne)
double const eps(0.0);
GeoLib::AABB aabb(ps_ptr.begin(), ps_ptr.end());
+ auto const min(convertToEigen(aabb.getMinPoint()));
+ auto const max(convertToEigen(aabb.getMaxPoint()));
std::unique_ptr<GeoLib::OctTree<GeoLib::Point, 2>> oct_tree(
- GeoLib::OctTree<GeoLib::Point, 2>::createOctTree(
- aabb.getMinPoint(), aabb.getMaxPoint(), eps));
+ GeoLib::OctTree<GeoLib::Point, 2>::createOctTree(min, max, eps));
GeoLib::Point* ret_pnt(nullptr);
ASSERT_TRUE(oct_tree->addPoint(ps_ptr[0], ret_pnt));
@@ -400,9 +411,10 @@ TEST_F(GeoLibOctTree, TestOctTreeOnCubicDomain)
double const eps(0.0);
GeoLib::AABB aabb(ps_ptr.begin(), ps_ptr.end());
+ auto const min(convertToEigen(aabb.getMinPoint()));
+ auto const max(convertToEigen(aabb.getMaxPoint()));
std::unique_ptr<GeoLib::OctTree<GeoLib::Point, 2>> oct_tree(
- GeoLib::OctTree<GeoLib::Point, 2>::createOctTree(
- aabb.getMinPoint(), aabb.getMaxPoint(), eps));
+ GeoLib::OctTree<GeoLib::Point, 2>::createOctTree(min, max, eps));
GeoLib::Point* ret_pnt(nullptr);
ASSERT_TRUE(oct_tree->addPoint(ps_ptr[0], ret_pnt));
--
GitLab