31 template <
typename POINT>
51 template <
typename InputIterator>
52 Grid(InputIterator first, InputIterator last,
53 std::size_t max_num_per_grid_cell = 512);
94 std::vector<std::vector<POINT*>
const*>
96 double half_len)
const;
98 std::vector<std::vector<POINT*>
const*>
100 Eigen::Vector3d
const& min_pnt, Eigen::Vector3d
const& max_pnt)
const;
126 std::array<double, 3>
const& extensions);
135 template <
typename T>
168 template <
typename P>
170 P
const& pnt, std::array<std::size_t, 3>
const& coords)
const;
172 template <
typename P>
174 std::array<std::size_t, 3>
const& coords,
175 double& sqr_min_dist,
176 POINT*& nearest_pnt)
const;
190 template <
typename POINT>
191 template <
typename InputIterator>
193 std::size_t max_num_per_grid_cell)
196 auto const n_pnts(std::distance(first, last));
200 std::array<double, 3> delta = {
201 {max[0] - min[0], max[1] - min[1], max[2] - min[2]}};
204 for (
auto& d : delta)
206 d = std::nextafter(d, std::numeric_limits<double>::max());
217 for (std::size_t k(0); k < 3; k++)
219 if (std::abs(delta[k]) < std::numeric_limits<double>::epsilon())
221 delta[k] = std::numeric_limits<double>::epsilon();
227 InputIterator it(first);
233 std::size_t
const pos(coords[0] + coords[1] *
_n_steps[0] +
234 coords[2] * n_plane);
240 ERR(
"Grid constructor: error computing indices [{:d}, {:d}, {:d}], "
241 "max indices [{:d}, {:d}, {:d}].",
249 template <
typename POINT>
250 template <
typename P>
251 std::vector<std::vector<POINT*>
const*>
253 double half_len)
const
255 Eigen::Vector3d
const c{center[0], center[1], center[2]};
256 std::array<std::size_t, 3> min_coords(getGridCoords(
c.array() - half_len));
257 std::array<std::size_t, 3> max_coords(getGridCoords(
c.array() + half_len));
259 std::vector<std::vector<POINT*>
const*> pnts;
261 (Eigen::Map<Eigen::Matrix<std::size_t, 3, 1>>(max_coords.data()) -
262 Eigen::Map<Eigen::Matrix<std::size_t, 3, 1>>(min_coords.data()))
265 std::size_t
const steps0_x_steps1(_n_steps[0] * _n_steps[1]);
266 for (std::size_t c0 = min_coords[0]; c0 < max_coords[0] + 1; c0++)
268 for (std::size_t c1 = min_coords[1]; c1 < max_coords[1] + 1; c1++)
270 const std::size_t coords0_p_coords1_x_steps0(c0 + c1 * _n_steps[0]);
271 for (std::size_t c2 = min_coords[2]; c2 < max_coords[2] + 1; c2++)
274 &(_grid_cell_nodes_map[coords0_p_coords1_x_steps0 +
275 c2 * steps0_x_steps1]));
282 template <
typename POINT>
283 std::vector<std::vector<POINT*>
const*>
285 Eigen::Vector3d
const& min_pnt, Eigen::Vector3d
const& max_pnt)
const
287 std::array<std::size_t, 3> min_coords(getGridCoords(min_pnt));
288 std::array<std::size_t, 3> max_coords(getGridCoords(max_pnt));
290 std::vector<std::vector<POINT*>
const*> pnts;
292 (Eigen::Map<Eigen::Matrix<std::size_t, 3, 1>>(max_coords.data()) -
293 Eigen::Map<Eigen::Matrix<std::size_t, 3, 1>>(min_coords.data()))
296 std::size_t
const steps0_x_steps1(_n_steps[0] * _n_steps[1]);
297 for (std::size_t c0 = min_coords[0]; c0 < max_coords[0] + 1; c0++)
299 for (std::size_t c1 = min_coords[1]; c1 < max_coords[1] + 1; c1++)
301 const std::size_t coords0_p_coords1_x_steps0(c0 + c1 * _n_steps[0]);
302 for (std::size_t c2 = min_coords[2]; c2 < max_coords[2] + 1; c2++)
305 &(_grid_cell_nodes_map[coords0_p_coords1_x_steps0 +
306 c2 * steps0_x_steps1]));
314 template <
typename POINT>
317 std::vector<std::string> grid_names;
322 const double dx((urb[0] - llf[0]) / _n_steps[0]);
323 const double dy((urb[1] - llf[1]) / _n_steps[1]);
324 const double dz((urb[2] - llf[2]) / _n_steps[2]);
327 for (std::size_t i(0); i < _n_steps[0]; i++)
329 for (std::size_t j(0); j < _n_steps[1]; j++)
331 for (std::size_t k(0); k < _n_steps[2]; k++)
333 std::string
name(
"Grid-");
334 name += std::to_string(i);
336 name += std::to_string(j);
338 name += std::to_string(k);
339 grid_names.push_back(
name);
343 std::make_unique<std::vector<GeoLib::Point*>>();
345 llf[0] + i * dx, llf[1] + j * dy, llf[2] + k * dz));
347 llf[1] + (j + 1) * dy,
350 llf[1] + (j + 1) * dy,
357 llf[2] + (k + 1) * dz));
359 llf[1] + (j + 1) * dy,
360 llf[2] + (k + 1) * dz));
362 llf[1] + (j + 1) * dy,
363 llf[2] + (k + 1) * dz));
366 llf[2] + (k + 1) * dz));
367 geo_obj->
addPointVec(std::move(points), grid_names.back(),
371 auto plys = std::make_unique<std::vector<GeoLib::Polyline*>>();
372 auto const& points = *geo_obj->
getPointVec(grid_names.back());
375 for (std::size_t l(0); l < 4; l++)
378 plys->push_back(ply0);
381 for (std::size_t l(4); l < 8; l++)
384 plys->push_back(ply1);
389 plys->push_back(ply2);
394 plys->push_back(ply3);
399 plys->push_back(ply4);
404 plys->push_back(ply5);
411 std::string merged_geo_name(
"Grid");
417 template <
typename POINT>
418 template <
typename T>
421 auto const& min_point{getMinPoint()};
422 auto const& max_point{getMinPoint()};
423 std::array<std::size_t, 3> coords{0, 0, 0};
424 for (std::size_t k(0); k < 3; k++)
426 if (pnt[k] < min_point[k])
430 if (pnt[k] >= max_point[k])
432 coords[k] = _n_steps[k] - 1;
435 coords[k] =
static_cast<std::size_t
>(
436 std::floor((pnt[k] - min_point[k])) /
437 std::nextafter(_step_sizes[k], std::numeric_limits<double>::max()));
442 template <
typename POINT>
443 template <
typename P>
445 P
const& pnt, std::array<std::size_t, 3>
const& coords)
const
447 auto const& min_point{getMinPoint()};
448 std::array<double, 6> dists{};
450 std::abs(pnt[2] - min_point[2] + coords[2] * _step_sizes[2]);
451 dists[5] = std::abs(pnt[2] - min_point[2] +
452 (coords[2] + 1) * _step_sizes[2]);
455 std::abs(pnt[1] - min_point[1] + coords[1] * _step_sizes[1]);
456 dists[3] = std::abs(pnt[1] - min_point[1] +
457 (coords[1] + 1) * _step_sizes[1]);
460 std::abs(pnt[0] - min_point[0] + coords[0] * _step_sizes[0]);
461 dists[2] = std::abs(pnt[0] - min_point[0] +
462 (coords[0] + 1) * _step_sizes[0]);
466 template <
typename POINT>
467 template <
typename P>
470 std::array<std::size_t, 3> coords(getGridCoords(pnt));
471 auto const& min_point{getMinPoint()};
472 auto const& max_point{getMaxPoint()};
474 double sqr_min_dist = (max_point - min_point).squaredNorm();
475 POINT* nearest_pnt(
nullptr);
477 std::array<double, 6> dists(getPointCellBorderDistances(pnt, coords));
479 if (calcNearestPointInGridCell(pnt, coords, sqr_min_dist, nearest_pnt))
481 double min_dist(std::sqrt(sqr_min_dist));
482 if (dists[0] >= min_dist && dists[1] >= min_dist &&
483 dists[2] >= min_dist && dists[3] >= min_dist &&
484 dists[4] >= min_dist && dists[5] >= min_dist)
492 double sqr_min_dist_tmp;
493 POINT* nearest_pnt_tmp(
nullptr);
494 std::size_t offset(1);
496 while (nearest_pnt ==
nullptr)
498 std::array<std::size_t, 3> ijk{
499 {coords[0] < offset ? 0 : coords[0] - offset,
500 coords[1] < offset ? 0 : coords[1] - offset,
501 coords[2] < offset ? 0 : coords[2] - offset}};
502 for (; ijk[0] < coords[0] + offset; ijk[0]++)
504 for (; ijk[1] < coords[1] + offset; ijk[1]++)
506 for (; ijk[2] < coords[2] + offset; ijk[2]++)
509 if (ijk[0] == coords[0] && ijk[1] == coords[1] &&
515 if (ijk[0] >= _n_steps[0] || ijk[1] >= _n_steps[1] ||
516 ijk[2] >= _n_steps[2])
521 if (calcNearestPointInGridCell(
522 pnt, ijk, sqr_min_dist_tmp, nearest_pnt_tmp))
524 if (sqr_min_dist_tmp < sqr_min_dist)
526 sqr_min_dist = sqr_min_dist_tmp;
527 nearest_pnt = nearest_pnt_tmp;
537 auto to_eigen = [](
auto const& point)
538 {
return Eigen::Map<Eigen::Vector3d const>(point.getCoords()); };
540 double len((to_eigen(pnt) - to_eigen(*nearest_pnt)).
norm());
542 std::vector<std::vector<POINT*>
const*> vecs_of_pnts(
543 getPntVecsOfGridCellsIntersectingCube(pnt, len));
545 const std::size_t n_vecs(vecs_of_pnts.size());
546 for (std::size_t j(0); j < n_vecs; j++)
548 std::vector<POINT*>
const& pnts(*(vecs_of_pnts[j]));
549 const std::size_t n_pnts(pnts.size());
550 for (std::size_t k(0); k < n_pnts; k++)
552 const double sqr_dist(
553 (to_eigen(pnt) - to_eigen(*pnts[k])).squaredNorm());
554 if (sqr_dist < sqr_min_dist)
556 sqr_min_dist = sqr_dist;
557 nearest_pnt = pnts[k];
565 template <
typename POINT>
567 std::array<double, 3>
const& extensions)
569 double const max_extension(
570 *std::max_element(extensions.cbegin(), extensions.cend()));
573 for (std::size_t k(0); k < 3; ++k)
576 if (extensions[k] >= 1e-4 * max_extension)
591 auto sc_ceil = [](
double v) {
return static_cast<std::size_t
>(ceil(v)); };
596 _n_steps[0] = sc_ceil(
597 std::cbrt(n_pnts * (extensions[0] / extensions[1]) *
598 (extensions[0] / extensions[2]) / n_per_cell));
599 _n_steps[1] = sc_ceil(
600 _n_steps[0] * std::min(extensions[1] / extensions[0], 100.0));
601 _n_steps[2] = sc_ceil(
602 _n_steps[0] * std::min(extensions[2] / extensions[0], 100.0));
605 if (dim[0] && dim[1])
607 _n_steps[0] = sc_ceil(std::sqrt(n_pnts * extensions[0] /
608 (n_per_cell * extensions[1])));
610 sc_ceil(_n_steps[0] *
611 std::min(extensions[1] / extensions[0], 100.0));
613 else if (dim[0] && dim[2])
615 _n_steps[0] = sc_ceil(std::sqrt(n_pnts * extensions[0] /
616 (n_per_cell * extensions[2])));
618 sc_ceil(_n_steps[0] *
619 std::min(extensions[2] / extensions[0], 100.0));
621 else if (dim[1] && dim[2])
623 _n_steps[1] = sc_ceil(std::sqrt(n_pnts * extensions[1] /
624 (n_per_cell * extensions[2])));
626 sc_ceil(std::min(extensions[2] / extensions[1], 100.0));
630 for (std::size_t k(0); k < 3; ++k)
635 sc_ceil(
static_cast<double>(n_pnts) / n_per_cell);
641 template <
typename POINT>
642 template <
typename P>
645 std::array<std::size_t, 3>
const& coords,
646 double& sqr_min_dist,
647 POINT*& nearest_pnt)
const
649 const std::size_t grid_idx(coords[0] + coords[1] * _n_steps[0] +
650 coords[2] * _n_steps[0] * _n_steps[1]);
651 std::vector<
typename std::add_pointer<
652 typename std::remove_pointer<POINT>::type>::type>
const&
653 pnts(_grid_cell_nodes_map[grid_idx]);
657 auto to_eigen = [](
auto const& point)
658 {
return Eigen::Map<Eigen::Vector3d const>(point.getCoords()); };
660 const std::size_t n_pnts(pnts.size());
661 sqr_min_dist = (to_eigen(*pnts[0]) - to_eigen(pnt)).squaredNorm();
662 nearest_pnt = pnts[0];
663 for (std::size_t i(1); i < n_pnts; i++)
665 const double sqr_dist(
666 (to_eigen(*pnts[i]) - to_eigen(pnt)).squaredNorm());
667 if (sqr_dist < sqr_min_dist)
669 sqr_min_dist = sqr_dist;
670 nearest_pnt = pnts[i];
676 template <
typename POINT>
677 template <
typename P>
679 P
const& pnt,
double eps)
const
681 std::vector<std::vector<POINT*>
const*> vec_pnts(
682 getPntVecsOfGridCellsIntersectingCube(pnt, eps));
684 double const sqr_eps(eps * eps);
686 auto to_eigen = [](
auto const& point)
687 {
return Eigen::Map<Eigen::Vector3d const>(point.getCoords()); };
689 std::vector<std::size_t> pnts;
690 for (
auto vec : vec_pnts)
692 for (
auto const p : *vec)
694 if ((to_eigen(*p) - to_eigen(pnt)).squaredNorm() <= sqr_eps)
696 pnts.push_back(p->getID());
Definition of the AABB class.
Definition of the GEOObjects class.
void ERR(char const *fmt, Args const &... args)
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Eigen::Vector3d const & getMinPoint() const
Eigen::Vector3d const & getMaxPoint() const
Container class for geometric objects.
const std::vector< Point * > * getPointVec(const std::string &name) const
void addPointVec(std::unique_ptr< std::vector< Point * >> points, std::string &name, std::unique_ptr< std::map< std::string, std::size_t >> pnt_id_name_map=nullptr, double eps=std::sqrt(std::numeric_limits< double >::epsilon()))
int mergeGeometries(std::vector< std::string > const &geo_names, std::string &merged_geo_name)
void addPolylineVec(std::unique_ptr< std::vector< Polyline * >> lines, const std::string &name, std::unique_ptr< std::map< std::string, std::size_t >> ply_names=nullptr)
POINT * getNearestPoint(P const &pnt) const
std::array< std::size_t, 3 > getGridCoords(T const &pnt) const
std::vector< std::vector< POINT * > const * > getPntVecsOfGridCellsIntersectingCuboid(Eigen::Vector3d const &min_pnt, Eigen::Vector3d const &max_pnt) const
static POINT * copyOrAddress(POINT &p)
static POINT * copyOrAddress(POINT *p)
std::array< double, 6 > getPointCellBorderDistances(P const &pnt, std::array< std::size_t, 3 > const &coords) const
void createGridGeometry(GeoLib::GEOObjects *geo_obj) const
bool calcNearestPointInGridCell(P const &pnt, std::array< std::size_t, 3 > const &coords, double &sqr_min_dist, POINT *&nearest_pnt) const
Grid(InputIterator first, InputIterator last, std::size_t max_num_per_grid_cell=512)
The constructor of the grid object takes a vector of points or nodes. Furthermore the user can specif...
std::vector< std::vector< POINT * > const * > getPntVecsOfGridCellsIntersectingCube(P const ¢er, double half_len) const
static POINT const * copyOrAddress(POINT const &p)
std::array< std::size_t, 3 > _n_steps
std::array< double, 3 > _step_sizes
std::vector< POINT * > * _grid_cell_nodes_map
std::vector< std::size_t > getPointsInEpsilonEnvironment(P const &pnt, double eps) const
void initNumberOfSteps(std::size_t n_per_cell, std::size_t n_pnts, std::array< double, 3 > const &extensions)
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
double norm(MatrixOrVector const &x, MathLib::VecNormType type)