OGS
Grid.h
Go to the documentation of this file.
1
15#pragma once
16
17#include <array>
18#include <bitset>
19#include <vector>
20
21#include "BaseLib/Logging.h"
22
23// GeoLib
24#include "AABB.h"
25#ifndef NDEBUG
26#include "GEOObjects.h"
27#endif
28
29namespace GeoLib
30{
31template <typename POINT>
32class Grid final : public GeoLib::AABB
33{
34public:
51 template <typename InputIterator>
52 Grid(InputIterator first, InputIterator last,
53 std::size_t max_num_per_grid_cell = 1);
54
55 Grid(Grid const&) = delete;
56 Grid& operator=(Grid const&) = delete;
57
62 virtual ~Grid() { delete[] _grid_cell_nodes_map; }
63
78 template <typename P>
79 POINT* getNearestPoint(P const& pnt) const;
80
81 template <typename P>
82 std::vector<std::size_t> getPointsInEpsilonEnvironment(P const& pnt,
83 double eps) const;
84
96 template <typename P>
97 std::vector<std::vector<POINT*> const*>
99 double half_len) const;
100
101 std::vector<std::vector<POINT*> const*>
103 Eigen::Vector3d const& min_pnt, Eigen::Vector3d const& max_pnt) const;
104
105#ifndef NDEBUG
112#endif
113
114private:
128 void initNumberOfSteps(std::size_t n_per_cell, std::size_t n_pnts,
129 std::array<double, 3> const& extensions);
130
138 template <typename T>
139 std::array<std::size_t, 3> getGridCoords(T const& pnt) const;
140
171 template <typename P>
172 std::array<double, 6> getPointCellBorderDistances(
173 P const& pnt, std::array<std::size_t, 3> const& coords) const;
174
175 template <typename P>
176 bool calcNearestPointInGridCell(P const& pnt,
177 std::array<std::size_t, 3> const& coords,
178 double& sqr_min_dist,
179 POINT*& nearest_pnt) const;
180
181 static POINT* copyOrAddress(POINT& p) { return &p; }
182 static POINT const* copyOrAddress(POINT const& p) { return &p; }
183 static POINT* copyOrAddress(POINT* p) { return p; }
184
185 std::array<std::size_t, 3> _n_steps = {{1, 1, 1}};
186 std::array<double, 3> _step_sizes = {{0.0, 0.0, 0.0}};
190 std::vector<POINT*>* _grid_cell_nodes_map = nullptr;
191};
192
193template <typename POINT>
194template <typename InputIterator>
195Grid<POINT>::Grid(InputIterator first, InputIterator last,
196 std::size_t max_num_per_grid_cell)
197 : GeoLib::AABB(first, last)
198{
199 auto const n_pnts(std::distance(first, last));
200
201 auto const& max{getMaxPoint()};
202 auto const& min{getMinPoint()};
203 std::array<double, 3> delta = {
204 {max[0] - min[0], max[1] - min[1], max[2] - min[2]}};
205
206 // enlarge delta
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); });
210
211 assert(n_pnts > 0);
212 initNumberOfSteps(max_num_per_grid_cell, static_cast<std::size_t>(n_pnts),
213 delta);
214
215 const std::size_t n_plane(_n_steps[0] * _n_steps[1]);
216 _grid_cell_nodes_map = new std::vector<POINT*>[n_plane * _n_steps[2]];
217
218 // some frequently used expressions to fill the grid vectors
219 for (std::size_t k(0); k < 3; k++)
220 {
221 if (std::abs(delta[k]) < std::numeric_limits<double>::epsilon())
222 {
223 delta[k] = std::numeric_limits<double>::epsilon();
224 }
225 _step_sizes[k] = delta[k] / _n_steps[k];
226 }
227
228 // fill the grid vectors
229 InputIterator it(first);
230 while (it != last)
231 {
232 std::array<std::size_t, 3> coords(getGridCoords(*copyOrAddress(*it)));
233 if (coords < _n_steps)
234 {
235 std::size_t const pos(coords[0] + coords[1] * _n_steps[0] +
236 coords[2] * n_plane);
237 _grid_cell_nodes_map[pos].push_back(
238 const_cast<POINT*>(copyOrAddress(*it)));
239 }
240 else
241 {
242 ERR("Grid constructor: error computing indices [{:d}, {:d}, {:d}], "
243 "max indices [{:d}, {:d}, {:d}].",
244 coords[0], coords[1], coords[2], _n_steps[0], _n_steps[1],
245 _n_steps[2]);
246 }
247 it++;
248 }
249}
250
251template <typename POINT>
252template <typename P>
253std::vector<std::vector<POINT*> const*>
255 double half_len) const
256{
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));
260
261 std::vector<std::vector<POINT*> const*> pnts;
262 pnts.reserve(
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()))
265 .prod());
266
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++)
269 {
270 for (std::size_t c1 = min_coords[1]; c1 < max_coords[1] + 1; c1++)
271 {
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++)
274 {
275 pnts.push_back(
276 &(_grid_cell_nodes_map[coords0_p_coords1_x_steps0 +
277 c2 * steps0_x_steps1]));
278 }
279 }
280 }
281 return pnts;
282}
283
284template <typename POINT>
285std::vector<std::vector<POINT*> const*>
287 Eigen::Vector3d const& min_pnt, Eigen::Vector3d const& max_pnt) const
288{
289 std::array<std::size_t, 3> min_coords(getGridCoords(min_pnt));
290 std::array<std::size_t, 3> max_coords(getGridCoords(max_pnt));
291
292 std::vector<std::vector<POINT*> const*> pnts;
293 pnts.reserve(
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()))
296 .prod());
297
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++)
300 {
301 for (std::size_t c1 = min_coords[1]; c1 < max_coords[1] + 1; c1++)
302 {
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++)
305 {
306 pnts.push_back(
307 &(_grid_cell_nodes_map[coords0_p_coords1_x_steps0 +
308 c2 * steps0_x_steps1]));
309 }
310 }
311 }
312 return pnts;
313}
314
315#ifndef NDEBUG
316template <typename POINT>
318{
319 std::vector<std::string> grid_names;
320
321 GeoLib::Point const& llf(getMinPoint());
322 GeoLib::Point const& urb(getMaxPoint());
323
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]);
327
328 // create grid names and grid boxes as geometry
329 for (std::size_t i(0); i < _n_steps[0]; i++)
330 {
331 for (std::size_t j(0); j < _n_steps[1]; j++)
332 {
333 for (std::size_t k(0); k < _n_steps[2]; k++)
334 {
335 std::string name("Grid-");
336 name += std::to_string(i);
337 name += "-";
338 name += std::to_string(j);
339 name += "-";
340 name += std::to_string(k);
341 grid_names.push_back(name);
342
343 {
344 std::vector<GeoLib::Point*> points;
345 points.push_back(new GeoLib::Point(
346 llf[0] + i * dx, llf[1] + j * dy, llf[2] + k * dz));
347 points.push_back(new GeoLib::Point(llf[0] + i * dx,
348 llf[1] + (j + 1) * dy,
349 llf[2] + k * dz));
350 points.push_back(new GeoLib::Point(llf[0] + (i + 1) * dx,
351 llf[1] + (j + 1) * dy,
352 llf[2] + k * dz));
353 points.push_back(new GeoLib::Point(llf[0] + (i + 1) * dx,
354 llf[1] + j * dy,
355 llf[2] + k * dz));
356 points.push_back(new GeoLib::Point(llf[0] + i * dx,
357 llf[1] + j * dy,
358 llf[2] + (k + 1) * dz));
359 points.push_back(new GeoLib::Point(llf[0] + i * dx,
360 llf[1] + (j + 1) * dy,
361 llf[2] + (k + 1) * dz));
362 points.push_back(new GeoLib::Point(llf[0] + (i + 1) * dx,
363 llf[1] + (j + 1) * dy,
364 llf[2] + (k + 1) * dz));
365 points.push_back(new GeoLib::Point(llf[0] + (i + 1) * dx,
366 llf[1] + j * dy,
367 llf[2] + (k + 1) * dz));
368 geo_obj->addPointVec(std::move(points), grid_names.back());
369 }
370
371 std::vector<GeoLib::Polyline*> plys;
372 auto const& points = *geo_obj->getPointVec(grid_names.back());
373
374 auto* ply0(new GeoLib::Polyline(points));
375 for (std::size_t l(0); l < 4; l++)
376 ply0->addPoint(l);
377 ply0->addPoint(0);
378 plys.push_back(ply0);
379
380 auto* ply1(new GeoLib::Polyline(points));
381 for (std::size_t l(4); l < 8; l++)
382 ply1->addPoint(l);
383 ply1->addPoint(4);
384 plys.push_back(ply1);
385
386 auto* ply2(new GeoLib::Polyline(points));
387 ply2->addPoint(0);
388 ply2->addPoint(4);
389 plys.push_back(ply2);
390
391 auto* ply3(new GeoLib::Polyline(points));
392 ply3->addPoint(1);
393 ply3->addPoint(5);
394 plys.push_back(ply3);
395
396 auto* ply4(new GeoLib::Polyline(points));
397 ply4->addPoint(2);
398 ply4->addPoint(6);
399 plys.push_back(ply4);
400
401 auto* ply5(new GeoLib::Polyline(points));
402 ply5->addPoint(3);
403 ply5->addPoint(7);
404 plys.push_back(ply5);
405
406 geo_obj->addPolylineVec(std::move(plys), grid_names.back(),
408 }
409 }
410 }
411 std::string merged_geo_name("Grid");
412
413 geo_obj->mergeGeometries(grid_names, merged_geo_name);
414}
415#endif
416
417template <typename POINT>
418template <typename T>
419std::array<std::size_t, 3> Grid<POINT>::getGridCoords(T const& pnt) const
420{
421 auto const [min_point, max_point] = getMinMaxPoints();
422
423 std::array<std::size_t, 3> coords{0, 0, 0};
424 for (std::size_t k(0); k < 3; k++)
425 {
426 if (pnt[k] < min_point[k])
427 {
428 continue;
429 }
430 if (pnt[k] >= max_point[k])
431 {
432 coords[k] = _n_steps[k] - 1;
433 continue;
434 }
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())));
439 }
440 return coords;
441}
442
443template <typename POINT>
444template <typename P>
446 P const& pnt, std::array<std::size_t, 3> const& coords) const
447{
448 auto const& min_point{getMinPoint()};
449 std::array<double, 6> dists{};
450 dists[0] =
451 std::abs(pnt[2] - min_point[2] + coords[2] * _step_sizes[2]); // bottom
452 dists[5] = std::abs(pnt[2] - min_point[2] +
453 (coords[2] + 1) * _step_sizes[2]); // top
454
455 dists[1] =
456 std::abs(pnt[1] - min_point[1] + coords[1] * _step_sizes[1]); // front
457 dists[3] = std::abs(pnt[1] - min_point[1] +
458 (coords[1] + 1) * _step_sizes[1]); // back
459
460 dists[4] =
461 std::abs(pnt[0] - min_point[0] + coords[0] * _step_sizes[0]); // left
462 dists[2] = std::abs(pnt[0] - min_point[0] +
463 (coords[0] + 1) * _step_sizes[0]); // right
464 return dists;
465}
466
467template <typename POINT>
468template <typename P>
470{
471 std::array<std::size_t, 3> coords(getGridCoords(pnt));
472 auto const& min_point{getMinPoint()};
473 auto const& max_point{getMaxPoint()};
474
475 double sqr_min_dist = (max_point - min_point).squaredNorm();
476 POINT* nearest_pnt(nullptr);
477
478 std::array<double, 6> const dists(getPointCellBorderDistances(pnt, coords));
479
480 if (calcNearestPointInGridCell(pnt, coords, sqr_min_dist, nearest_pnt))
481 {
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)
486 {
487 return nearest_pnt;
488 }
489 }
490 else
491 {
492 // search in all border cells for at least one neighbor
493 double sqr_min_dist_tmp;
494 POINT* nearest_pnt_tmp(nullptr);
495 std::size_t offset(1);
496
497 while (nearest_pnt == nullptr)
498 {
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]++)
504 {
505 for (; ijk[1] < coords[1] + offset; ijk[1]++)
506 {
507 for (; ijk[2] < coords[2] + offset; ijk[2]++)
508 {
509 // do not check the origin grid cell twice
510 if (ijk[0] == coords[0] && ijk[1] == coords[1] &&
511 ijk[2] == coords[2])
512 {
513 continue;
514 }
515 // check if temporary grid cell coordinates are valid
516 if (ijk[0] >= _n_steps[0] || ijk[1] >= _n_steps[1] ||
517 ijk[2] >= _n_steps[2])
518 {
519 continue;
520 }
521
522 if (calcNearestPointInGridCell(
523 pnt, ijk, sqr_min_dist_tmp, nearest_pnt_tmp))
524 {
525 if (sqr_min_dist_tmp < sqr_min_dist)
526 {
527 sqr_min_dist = sqr_min_dist_tmp;
528 nearest_pnt = nearest_pnt_tmp;
529 }
530 }
531 } // end k
532 } // end j
533 } // end i
534 offset++;
535 } // end while
536 } // end else
537
538 double const len(
539 (pnt.asEigenVector3d() - nearest_pnt->asEigenVector3d()).norm());
540 // search all other grid cells within the cube with the edge nodes
541 std::vector<std::vector<POINT*> const*> vecs_of_pnts(
542 getPntVecsOfGridCellsIntersectingCube(pnt, len));
543
544 const std::size_t n_vecs(vecs_of_pnts.size());
545 for (std::size_t j(0); j < n_vecs; j++)
546 {
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++)
550 {
551 const double sqr_dist(
552 (pnt.asEigenVector3d() - pnts[k]->asEigenVector3d())
553 .squaredNorm());
554 if (sqr_dist < sqr_min_dist)
555 {
556 sqr_min_dist = sqr_dist;
557 nearest_pnt = pnts[k];
558 }
559 }
560 }
561
562 return nearest_pnt;
563}
564
565template <typename POINT>
566void Grid<POINT>::initNumberOfSteps(std::size_t n_per_cell, std::size_t n_pnts,
567 std::array<double, 3> const& extensions)
568{
569 double const max_extension(
570 *std::max_element(extensions.cbegin(), extensions.cend()));
571
572 std::bitset<3> dim; // all bits set to zero
573 for (std::size_t k(0); k < 3; ++k)
574 {
575 // set dimension if the ratio kth-extension/max_extension >= 1e-4
576 if (extensions[k] >= 1e-4 * max_extension)
577 {
578 dim[k] = true;
579 }
580 }
581
582 // structured grid: n_cells = _n_steps[0] * _n_steps[1] * _n_steps[2]
583 // *** condition: n_pnts / n_cells < n_per_cell
584 // => n_pnts / n_per_cell < n_cells
585 // _n_steps[1] = _n_steps[0] * extensions[1]/extensions[0],
586 // _n_steps[2] = _n_steps[0] * extensions[2]/extensions[0],
587 // => n_cells = _n_steps[0]^3 * extensions[1]/extensions[0] *
588 // extensions[2]/extensions[0],
589 // => _n_steps[0] = cbrt(n_cells * extensions[0]^2 /
590 // (extensions[1]*extensions[2]))
591 auto sc_ceil = [](double v) { return static_cast<std::size_t>(ceil(v)); };
592
593 switch (dim.count())
594 {
595 case 3: // 3d case
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));
603 break;
604 case 2: // 2d cases
605 if (dim[0] && dim[1])
606 { // xy
607 _n_steps[0] = sc_ceil(std::sqrt(n_pnts * extensions[0] /
608 (n_per_cell * extensions[1])));
609 _n_steps[1] =
610 sc_ceil(_n_steps[0] *
611 std::min(extensions[1] / extensions[0], 100.0));
612 }
613 else if (dim[0] && dim[2])
614 { // xz
615 _n_steps[0] = sc_ceil(std::sqrt(n_pnts * extensions[0] /
616 (n_per_cell * extensions[2])));
617 _n_steps[2] =
618 sc_ceil(_n_steps[0] *
619 std::min(extensions[2] / extensions[0], 100.0));
620 }
621 else if (dim[1] && dim[2])
622 { // yz
623 _n_steps[1] = sc_ceil(std::sqrt(n_pnts * extensions[1] /
624 (n_per_cell * extensions[2])));
625 _n_steps[2] =
626 sc_ceil(std::min(extensions[2] / extensions[1], 100.0));
627 }
628 break;
629 case 1: // 1d cases
630 for (std::size_t k(0); k < 3; ++k)
631 {
632 if (dim[k])
633 {
634 _n_steps[k] =
635 sc_ceil(static_cast<double>(n_pnts) / n_per_cell);
636 }
637 }
638 }
639}
640
641template <typename POINT>
642template <typename P>
644 P const& pnt,
645 std::array<std::size_t, 3> const& coords,
646 double& sqr_min_dist,
647 POINT*& nearest_pnt) const
648{
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]);
654 if (pnts.empty())
655 return false;
656
657 const std::size_t n_pnts(pnts.size());
658 sqr_min_dist =
659 (pnts[0]->asEigenVector3d() - pnt.asEigenVector3d()).squaredNorm();
660 nearest_pnt = pnts[0];
661 for (std::size_t i(1); i < n_pnts; i++)
662 {
663 const double sqr_dist(
664 (pnts[i]->asEigenVector3d() - pnt.asEigenVector3d()).squaredNorm());
665 if (sqr_dist < sqr_min_dist)
666 {
667 sqr_min_dist = sqr_dist;
668 nearest_pnt = pnts[i];
669 }
670 }
671 return true;
672}
673
674template <typename POINT>
675template <typename P>
677 P const& pnt, double eps) const
678{
679 std::vector<std::vector<POINT*> const*> vec_pnts(
680 getPntVecsOfGridCellsIntersectingCube(pnt, eps));
681
682 double const sqr_eps(eps * eps);
683
684 std::vector<std::size_t> pnts;
685 for (auto vec : vec_pnts)
686 {
687 for (auto const p : *vec)
688 {
689 if ((p->asEigenVector3d() - pnt.asEigenVector3d()).squaredNorm() <=
690 sqr_eps)
691 {
692 pnts.push_back(p->getID());
693 }
694 }
695 }
696
697 return pnts;
698}
699
700} // end namespace GeoLib
Definition of the AABB class.
Definition of the GEOObjects class.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Definition AABB.h:56
Eigen::Vector3d const & getMaxPoint() const
Definition AABB.h:187
Eigen::Vector3d const & getMinPoint() const
Definition AABB.h:180
Container class for geometric objects.
Definition GEOObjects.h:57
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
Definition Grid.h:286
POINT * getNearestPoint(P const &pnt) const
Definition Grid.h:469
std::array< std::size_t, 3 > getGridCoords(T const &pnt) const
Definition Grid.h:419
static POINT * copyOrAddress(POINT &p)
Definition Grid.h:181
Grid(Grid const &)=delete
static POINT const * copyOrAddress(POINT const &p)
Definition Grid.h:182
std::array< double, 6 > getPointCellBorderDistances(P const &pnt, std::array< std::size_t, 3 > const &coords) const
Definition Grid.h:445
std::vector< std::vector< POINT * > const * > getPntVecsOfGridCellsIntersectingCube(P const &center, double half_len) const
Definition Grid.h:254
void createGridGeometry(GeoLib::GEOObjects *geo_obj) const
Definition Grid.h:317
bool calcNearestPointInGridCell(P const &pnt, std::array< std::size_t, 3 > const &coords, double &sqr_min_dist, POINT *&nearest_pnt) const
Definition Grid.h:643
std::array< std::size_t, 3 > _n_steps
Definition Grid.h:185
std::array< double, 3 > _step_sizes
Definition Grid.h:186
std::vector< POINT * > * _grid_cell_nodes_map
Definition Grid.h:190
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...
Definition Grid.h:195
virtual ~Grid()
Definition Grid.h:62
static POINT * copyOrAddress(POINT *p)
Definition Grid.h:183
std::vector< std::size_t > getPointsInEpsilonEnvironment(P const &pnt, double eps) const
Definition Grid.h:676
void initNumberOfSteps(std::size_t n_per_cell, std::size_t n_pnts, std::array< double, 3 > const &extensions)
Definition Grid.h:566
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
Definition Polyline.h:40
std::map< std::string, std::size_t > NameIdMap
Definition TemplateVec.h:41