31template <
typename POINT>
51 template <
typename InputIterator>
52 Grid(InputIterator first, InputIterator last,
53 std::size_t max_num_per_grid_cell = 1);
97 std::vector<std::vector<POINT*>
const*>
99 double half_len)
const;
101 std::vector<std::vector<POINT*>
const*>
103 Eigen::Vector3d
const& min_pnt, Eigen::Vector3d
const& max_pnt)
const;
129 std::array<double, 3>
const& extensions);
138 template <
typename T>
171 template <
typename P>
173 P
const& pnt, std::array<std::size_t, 3>
const& coords)
const;
175 template <
typename P>
177 std::array<std::size_t, 3>
const& coords,
178 double& sqr_min_dist,
179 POINT*& nearest_pnt)
const;
193template <
typename POINT>
194template <
typename InputIterator>
196 std::size_t max_num_per_grid_cell)
199 auto const n_pnts(std::distance(first, last));
203 std::array<double, 3> delta = {
204 {max[0] - min[0], max[1] - min[1], max[2] - min[2]}};
207 constexpr double direction = std::numeric_limits<double>::max();
208 std::transform(begin(delta), end(delta), begin(delta),
209 [](
double const d) {
return std::nextafter(d, direction); });
219 for (std::size_t k(0); k < 3; k++)
221 if (std::abs(delta[k]) < std::numeric_limits<double>::epsilon())
223 delta[k] = std::numeric_limits<double>::epsilon();
229 InputIterator it(first);
235 std::size_t
const pos(coords[0] + coords[1] *
_n_steps[0] +
236 coords[2] * n_plane);
242 ERR(
"Grid constructor: error computing indices [{:d}, {:d}, {:d}], "
243 "max indices [{:d}, {:d}, {:d}].",
251template <
typename POINT>
253std::vector<std::vector<POINT*>
const*>
255 double half_len)
const
257 Eigen::Vector3d
const c{center[0], center[1], center[2]};
258 std::array<std::size_t, 3> min_coords(getGridCoords(c.array() - half_len));
259 std::array<std::size_t, 3> max_coords(getGridCoords(c.array() + half_len));
261 std::vector<std::vector<POINT*>
const*> pnts;
263 (Eigen::Map<Eigen::Matrix<std::size_t, 3, 1>>(max_coords.data()) -
264 Eigen::Map<Eigen::Matrix<std::size_t, 3, 1>>(min_coords.data()))
267 std::size_t
const steps0_x_steps1(_n_steps[0] * _n_steps[1]);
268 for (std::size_t c0 = min_coords[0]; c0 < max_coords[0] + 1; c0++)
270 for (std::size_t c1 = min_coords[1]; c1 < max_coords[1] + 1; c1++)
272 const std::size_t coords0_p_coords1_x_steps0(c0 + c1 * _n_steps[0]);
273 for (std::size_t c2 = min_coords[2]; c2 < max_coords[2] + 1; c2++)
276 &(_grid_cell_nodes_map[coords0_p_coords1_x_steps0 +
277 c2 * steps0_x_steps1]));
284template <
typename POINT>
285std::vector<std::vector<POINT*>
const*>
287 Eigen::Vector3d
const& min_pnt, Eigen::Vector3d
const& max_pnt)
const
289 std::array<std::size_t, 3> min_coords(getGridCoords(min_pnt));
290 std::array<std::size_t, 3> max_coords(getGridCoords(max_pnt));
292 std::vector<std::vector<POINT*>
const*> pnts;
294 (Eigen::Map<Eigen::Matrix<std::size_t, 3, 1>>(max_coords.data()) -
295 Eigen::Map<Eigen::Matrix<std::size_t, 3, 1>>(min_coords.data()))
298 std::size_t
const steps0_x_steps1(_n_steps[0] * _n_steps[1]);
299 for (std::size_t c0 = min_coords[0]; c0 < max_coords[0] + 1; c0++)
301 for (std::size_t c1 = min_coords[1]; c1 < max_coords[1] + 1; c1++)
303 const std::size_t coords0_p_coords1_x_steps0(c0 + c1 * _n_steps[0]);
304 for (std::size_t c2 = min_coords[2]; c2 < max_coords[2] + 1; c2++)
307 &(_grid_cell_nodes_map[coords0_p_coords1_x_steps0 +
308 c2 * steps0_x_steps1]));
316template <
typename POINT>
319 std::vector<std::string> grid_names;
324 const double dx((urb[0] - llf[0]) / _n_steps[0]);
325 const double dy((urb[1] - llf[1]) / _n_steps[1]);
326 const double dz((urb[2] - llf[2]) / _n_steps[2]);
329 for (std::size_t i(0); i < _n_steps[0]; i++)
331 for (std::size_t j(0); j < _n_steps[1]; j++)
333 for (std::size_t k(0); k < _n_steps[2]; k++)
335 std::string name(
"Grid-");
336 name += std::to_string(i);
338 name += std::to_string(j);
340 name += std::to_string(k);
341 grid_names.push_back(name);
344 std::vector<GeoLib::Point*> points;
346 llf[0] + i * dx, llf[1] + j * dy, llf[2] + k * dz));
348 llf[1] + (j + 1) * dy,
351 llf[1] + (j + 1) * dy,
358 llf[2] + (k + 1) * dz));
360 llf[1] + (j + 1) * dy,
361 llf[2] + (k + 1) * dz));
363 llf[1] + (j + 1) * dy,
364 llf[2] + (k + 1) * dz));
367 llf[2] + (k + 1) * dz));
368 geo_obj->
addPointVec(std::move(points), grid_names.back());
371 std::vector<GeoLib::Polyline*> plys;
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");
417template <
typename POINT>
421 auto const [min_point, max_point] = getMinMaxPoints();
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],
438 std::numeric_limits<double>::max())));
443template <
typename POINT>
446 P
const& pnt, std::array<std::size_t, 3>
const& coords)
const
448 auto const& min_point{getMinPoint()};
449 std::array<double, 6> dists{};
451 std::abs(pnt[2] - min_point[2] + coords[2] * _step_sizes[2]);
452 dists[5] = std::abs(pnt[2] - min_point[2] +
453 (coords[2] + 1) * _step_sizes[2]);
456 std::abs(pnt[1] - min_point[1] + coords[1] * _step_sizes[1]);
457 dists[3] = std::abs(pnt[1] - min_point[1] +
458 (coords[1] + 1) * _step_sizes[1]);
461 std::abs(pnt[0] - min_point[0] + coords[0] * _step_sizes[0]);
462 dists[2] = std::abs(pnt[0] - min_point[0] +
463 (coords[0] + 1) * _step_sizes[0]);
467template <
typename POINT>
471 std::array<std::size_t, 3> coords(getGridCoords(pnt));
472 auto const& min_point{getMinPoint()};
473 auto const& max_point{getMaxPoint()};
475 double sqr_min_dist = (max_point - min_point).squaredNorm();
476 POINT* nearest_pnt(
nullptr);
478 std::array<double, 6>
const dists(getPointCellBorderDistances(pnt, coords));
480 if (calcNearestPointInGridCell(pnt, coords, sqr_min_dist, nearest_pnt))
482 double min_dist(std::sqrt(sqr_min_dist));
483 if (dists[0] >= min_dist && dists[1] >= min_dist &&
484 dists[2] >= min_dist && dists[3] >= min_dist &&
485 dists[4] >= min_dist && dists[5] >= min_dist)
493 double sqr_min_dist_tmp;
494 POINT* nearest_pnt_tmp(
nullptr);
495 std::size_t offset(1);
497 while (nearest_pnt ==
nullptr)
499 std::array<std::size_t, 3> ijk{
500 {coords[0] < offset ? 0 : coords[0] - offset,
501 coords[1] < offset ? 0 : coords[1] - offset,
502 coords[2] < offset ? 0 : coords[2] - offset}};
503 for (; ijk[0] < coords[0] + offset; ijk[0]++)
505 for (; ijk[1] < coords[1] + offset; ijk[1]++)
507 for (; ijk[2] < coords[2] + offset; ijk[2]++)
510 if (ijk[0] == coords[0] && ijk[1] == coords[1] &&
516 if (ijk[0] >= _n_steps[0] || ijk[1] >= _n_steps[1] ||
517 ijk[2] >= _n_steps[2])
522 if (calcNearestPointInGridCell(
523 pnt, ijk, sqr_min_dist_tmp, nearest_pnt_tmp))
525 if (sqr_min_dist_tmp < sqr_min_dist)
527 sqr_min_dist = sqr_min_dist_tmp;
528 nearest_pnt = nearest_pnt_tmp;
539 (pnt.asEigenVector3d() - nearest_pnt->asEigenVector3d()).norm());
541 std::vector<std::vector<POINT*>
const*> vecs_of_pnts(
542 getPntVecsOfGridCellsIntersectingCube(pnt, len));
544 const std::size_t n_vecs(vecs_of_pnts.size());
545 for (std::size_t j(0); j < n_vecs; j++)
547 std::vector<POINT*>
const& pnts(*(vecs_of_pnts[j]));
548 const std::size_t n_pnts(pnts.size());
549 for (std::size_t k(0); k < n_pnts; k++)
551 const double sqr_dist(
552 (pnt.asEigenVector3d() - pnts[k]->asEigenVector3d())
554 if (sqr_dist < sqr_min_dist)
556 sqr_min_dist = sqr_dist;
557 nearest_pnt = pnts[k];
565template <
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);
641template <
typename POINT>
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 const std::size_t n_pnts(pnts.size());
659 (pnts[0]->asEigenVector3d() - pnt.asEigenVector3d()).squaredNorm();
660 nearest_pnt = pnts[0];
661 for (std::size_t i(1); i < n_pnts; i++)
663 const double sqr_dist(
664 (pnts[i]->asEigenVector3d() - pnt.asEigenVector3d()).squaredNorm());
665 if (sqr_dist < sqr_min_dist)
667 sqr_min_dist = sqr_dist;
668 nearest_pnt = pnts[i];
674template <
typename POINT>
677 P
const& pnt,
double eps)
const
679 std::vector<std::vector<POINT*>
const*> vec_pnts(
680 getPntVecsOfGridCellsIntersectingCube(pnt, eps));
682 double const sqr_eps(eps * eps);
684 std::vector<std::size_t> pnts;
685 for (
auto vec : vec_pnts)
687 for (
auto const p : *vec)
689 if ((p->asEigenVector3d() - pnt.asEigenVector3d()).squaredNorm() <=
692 pnts.push_back(p->getID());
Definition of the AABB class.
Definition of the GEOObjects class.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Eigen::Vector3d const & getMaxPoint() const
Eigen::Vector3d const & getMinPoint() const
Container class for geometric objects.
void addPolylineVec(std::vector< Polyline * > &&lines, std::string const &name, PolylineVec::NameIdMap &&ply_names)
const std::vector< Point * > * getPointVec(const std::string &name) const
void addPointVec(std::vector< Point * > &&points, std::string &name, PointVec::NameIdMap &&pnt_id_name_map, double const eps=std::sqrt(std::numeric_limits< double >::epsilon()))
int mergeGeometries(std::vector< std::string > const &geo_names, std::string &merged_geo_name)
std::vector< std::vector< POINT * > const * > getPntVecsOfGridCellsIntersectingCuboid(Eigen::Vector3d const &min_pnt, Eigen::Vector3d const &max_pnt) const
POINT * getNearestPoint(P const &pnt) const
std::array< std::size_t, 3 > getGridCoords(T const &pnt) const
static POINT * copyOrAddress(POINT &p)
Grid(Grid const &)=delete
static POINT const * copyOrAddress(POINT const &p)
std::array< double, 6 > getPointCellBorderDistances(P const &pnt, std::array< std::size_t, 3 > const &coords) const
std::vector< std::vector< POINT * > const * > getPntVecsOfGridCellsIntersectingCube(P const ¢er, double half_len) 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
std::array< std::size_t, 3 > _n_steps
std::array< double, 3 > _step_sizes
std::vector< POINT * > * _grid_cell_nodes_map
Grid & operator=(Grid const &)=delete
Grid(InputIterator first, InputIterator last, std::size_t max_num_per_grid_cell=1)
The constructor of the grid object takes a vector of points or nodes. Furthermore the user can specif...
static POINT * copyOrAddress(POINT *p)
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 ...
std::map< std::string, std::size_t > NameIdMap