OGS
SurfaceGrid.cpp
Go to the documentation of this file.
1/*
2 * \file
3 * \date 2012-09-22
4 * \brief Definition of the SurfaceGrid class.
5 *
7 * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
9 * See accompanying file LICENSE.txt or
10 */
11
12#include "SurfaceGrid.h"
13
14#include <algorithm>
15#include <bitset>
16
17#include "BaseLib/Error.h"
18#include "BaseLib/Logging.h"
19#include "MathLib/Point3d.h"
20#include "Surface.h"
21#include "Triangle.h"
22
23namespace GeoLib
24{
26 : AABB(sfc->getAABB()), _n_steps({{1, 1, 1}})
27{
28 auto [min_point, max_point] = getMinMaxPoints();
29 // enlarge the bounding, such that the points with maximal coordinates
30 // fits into the grid
31 for (std::size_t k(0); k < 3; ++k)
32 {
33 max_point[k] += std::abs(max_point[k]) * 1e-6;
34 if (std::abs(max_point[k]) < std::numeric_limits<double>::epsilon())
35 {
36 max_point[k] = (max_point[k] - min_point[k]) * (1.0 + 1e-6);
37 }
38 }
39
40 Eigen::Vector3d delta = max_point - min_point;
41
42 if (delta[0] < std::numeric_limits<double>::epsilon())
43 {
44 const double max_delta(std::max(delta[1], delta[2]));
45 min_point[0] -= max_delta * 0.5e-3;
46 max_point[0] += max_delta * 0.5e-3;
47 delta[0] = max_point[0] - min_point[0];
48 }
49
50 if (delta[1] < std::numeric_limits<double>::epsilon())
51 {
52 const double max_delta(std::max(delta[0], delta[2]));
53 min_point[1] -= max_delta * 0.5e-3;
54 max_point[1] += max_delta * 0.5e-3;
55 delta[1] = max_point[1] - min_point[1];
56 }
57
58 if (delta[2] < std::numeric_limits<double>::epsilon())
59 {
60 const double max_delta(std::max(delta[0], delta[1]));
61 min_point[2] -= max_delta * 0.5e-3;
62 max_point[2] += max_delta * 0.5e-3;
63 delta[2] = max_point[2] - min_point[2];
64 }
65
66 update(min_point);
67 update(max_point);
68
69 const std::size_t n_tris(sfc->getNumberOfTriangles());
70 const std::size_t n_tris_per_cell(5);
71
72 Eigen::Matrix<bool, 3, 1> dim =
73 delta.array() >= std::numeric_limits<double>::epsilon();
74
75 // *** condition: n_tris / n_cells < n_tris_per_cell
76 // where n_cells = _n_steps[0] * _n_steps[1] * _n_steps[2]
77 // *** with _n_steps[0] =
78 // ceil(pow(n_tris*delta[0]*delta[0]/(n_tris_per_cell*delta[1]*delta[2]),
79 // 1/3.)));
80 // _n_steps[1] = _n_steps[0] * delta[1]/delta[0],
81 // _n_steps[2] = _n_steps[0] * delta[2]/delta[0]
82 auto sc_ceil = [](double v)
83 { return static_cast<std::size_t>(std::ceil(v)); };
84 switch (dim.count())
85 {
86 case 3: // 3d case
87 _n_steps[0] =
88 sc_ceil(std::cbrt(n_tris * delta[0] * delta[0] /
89 (n_tris_per_cell * delta[1] * delta[2])));
90 _n_steps[1] =
91 sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0));
92 _n_steps[2] =
93 sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0));
94 break;
95 case 2: // 2d cases
96 if (dim[0] && dim[2])
97 { // 2d case: xz plane, y = const
98 _n_steps[0] = sc_ceil(std::sqrt(n_tris * delta[0] /
99 (n_tris_per_cell * delta[2])));
100 _n_steps[2] = sc_ceil(_n_steps[0] * delta[2] / delta[0]);
101 }
102 else if (dim[0] && dim[1])
103 { // 2d case: xy plane, z = const
104 _n_steps[0] = sc_ceil(std::sqrt(n_tris * delta[0] /
105 (n_tris_per_cell * delta[1])));
106 _n_steps[1] = sc_ceil(_n_steps[0] * delta[1] / delta[0]);
107 }
108 else if (dim[1] && dim[2])
109 { // 2d case: yz plane, x = const
110 _n_steps[1] = sc_ceil(std::sqrt(n_tris * delta[1] /
111 (n_tris_per_cell * delta[2])));
112 _n_steps[2] =
113 sc_ceil(n_tris * delta[2] / (n_tris_per_cell * delta[1]));
114 }
115 break;
116 case 1: // 1d cases
117 for (std::size_t k(0); k < 3; ++k)
118 {
119 if (dim[k])
120 {
121 _n_steps[k] =
122 sc_ceil(static_cast<double>(n_tris) / n_tris_per_cell);
123 }
124 }
125 }
126
127 // some frequently used expressions to fill the grid vectors
128 for (std::size_t k(0); k < 3; k++)
129 {
130 _step_sizes[k] = delta[k] / _n_steps[k];
131 if (delta[k] > std::numeric_limits<double>::epsilon())
132 {
133 _inverse_step_sizes[k] = 1.0 / _step_sizes[k];
134 }
135 else
136 {
137 _inverse_step_sizes[k] = 0;
138 }
139 }
140
141 _triangles_in_grid_box.resize(_n_steps[0] * _n_steps[1] * _n_steps[2]);
142 sortTrianglesInGridCells(sfc);
143}
144
146{
147 for (std::size_t l(0); l < sfc->getNumberOfTriangles(); l++)
148 {
149 if (!sortTriangleInGridCells((*sfc)[l]))
150 {
151 Point const& p0(*((*sfc)[l]->getPoint(0)));
152 Point const& p1(*((*sfc)[l]->getPoint(1)));
153 Point const& p2(*((*sfc)[l]->getPoint(2)));
154 auto const [min, max] = getMinMaxPoints();
155 OGS_FATAL(
156 "Sorting triangle {:d} [({:f},{:f},{:f}), ({:f},{:f},{:f}), "
157 "({:f},{:f},{:f}) into grid. Bounding box is [{:f}, {:f}] x "
158 "[{:f}, {:f}] x [{:f}, {:f}].",
159 l, p0[0], p0[1], p0[2], p1[0], p1[1], p1[2], p2[0], p2[1],
160 p2[2], min[0], max[0], min[1], max[1], min[2], max[2]);
161 }
162 }
163}
164
166{
167 // compute grid coordinates for each triangle point
168 std::optional<std::array<std::size_t, 3> const> c_p0(
169 getGridCellCoordinates(*(triangle->getPoint(0))));
170 if (!c_p0)
171 {
172 return false;
173 }
174 std::optional<std::array<std::size_t, 3> const> c_p1(
175 getGridCellCoordinates(*(triangle->getPoint(1))));
176 if (!c_p1)
177 {
178 return false;
179 }
180 std::optional<std::array<std::size_t, 3> const> c_p2(
181 getGridCellCoordinates(*(triangle->getPoint(2))));
182 if (!c_p2)
183 {
184 return false;
185 }
186
187 // determine interval in grid (grid cells the triangle will be inserted)
188 std::size_t const i_min(
189 std::min(std::min((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0]));
190 std::size_t const i_max(
191 std::max(std::max((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0]));
192 std::size_t const j_min(
193 std::min(std::min((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1]));
194 std::size_t const j_max(
195 std::max(std::max((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1]));
196 std::size_t const k_min(
197 std::min(std::min((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2]));
198 std::size_t const k_max(
199 std::max(std::max((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2]));
200
201 const std::size_t n_plane(_n_steps[0] * _n_steps[1]);
202
203 // insert the triangle into the grid cells
204 for (std::size_t i(i_min); i <= i_max; i++)
205 {
206 for (std::size_t j(j_min); j <= j_max; j++)
207 {
208 for (std::size_t k(k_min); k <= k_max; k++)
209 {
210 _triangles_in_grid_box[i + j * _n_steps[0] + k * n_plane]
211 .push_back(triangle);
212 }
213 }
214 }
215
216 return true;
217}
218
219std::optional<std::array<std::size_t, 3>> SurfaceGrid::getGridCellCoordinates(
220 MathLib::Point3d const& p) const
221{
222 auto const& min_point{getMinPoint()};
223 std::array<std::size_t, 3> coords{
224 {static_cast<std::size_t>((p[0] - min_point[0]) *
226 static_cast<std::size_t>((p[1] - min_point[1]) *
228 static_cast<std::size_t>((p[2] - min_point[2]) *
230
231 if (coords[0] >= _n_steps[0] || coords[1] >= _n_steps[1] ||
232 coords[2] >= _n_steps[2])
233 {
234 DBUG(
235 "Computed indices ({:d},{:d},{:d}), max grid cell indices "
236 "({:d},{:d},{:d})",
237 coords[0], coords[1], coords[2], _n_steps[0], _n_steps[1],
238 _n_steps[2]);
239 return {};
240 }
241 return coords;
242}
243
245 double eps) const
246{
247 auto const grid_cell_coordinates = getGridCellCoordinates(pnt);
248 if (!grid_cell_coordinates)
249 {
250 return false;
251 }
252
253 std::size_t const grid_cell_index =
254 (*grid_cell_coordinates)[0] +
255 (*grid_cell_coordinates)[1] * _n_steps[0] +
256 (*grid_cell_coordinates)[2] * _n_steps[0] * _n_steps[1];
257
258 std::vector<Triangle const*> const& triangles(
259 _triangles_in_grid_box[grid_cell_index]);
260 return std::any_of(triangles.begin(), triangles.end(),
261 [eps, pnt](auto const* triangle)
262 { return triangle->containsPoint(pnt, eps); });
263}
264
265} // end namespace GeoLib
#define OGS_FATAL(...)
Definition Error.h:26
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Definition of the Point3d class.
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 & getMinPoint() const
Definition AABB.h:180
MinMaxPoints getMinMaxPoints() const
Definition AABB.h:174
void sortTrianglesInGridCells(GeoLib::Surface const *const sfc)
std::vector< std::vector< GeoLib::Triangle const * > > _triangles_in_grid_box
Definition SurfaceGrid.h:46
bool isPointInSurface(MathLib::Point3d const &pnt, double eps=std::numeric_limits< double >::epsilon()) const
std::array< std::size_t, 3 > _n_steps
Definition SurfaceGrid.h:45
std::array< double, 3 > _inverse_step_sizes
Definition SurfaceGrid.h:44
bool sortTriangleInGridCells(GeoLib::Triangle const *const triangle)
SurfaceGrid(GeoLib::Surface const *const sfc)
std::optional< std::array< std::size_t, 3 > > getGridCellCoordinates(MathLib::Point3d const &p) const
A Surface is represented by Triangles. It consists of a reference to a vector of (pointers to) points...
Definition Surface.h:33
std::size_t getNumberOfTriangles() const
Definition Surface.cpp:82
Class Triangle consists of a reference to a point vector and a vector that stores the indices in the ...
Definition Triangle.h:27
const Point * getPoint(std::size_t i) const
const access operator to access the i-th triangle Point
Definition Triangle.h:50