26 :
AABB(sfc->getAABB()), _n_steps({{1, 1, 1}})
28 auto min_point{getMinPoint()};
29 auto max_point{getMaxPoint()};
32 for (std::size_t k(0); k < 3; ++k)
34 max_point[k] += std::abs(max_point[k]) * 1e-6;
35 if (std::abs(max_point[k]) < std::numeric_limits<double>::epsilon())
37 max_point[k] = (max_point[k] - min_point[k]) * (1.0 + 1e-6);
41 Eigen::Vector3d delta = max_point - min_point;
43 if (delta[0] < std::numeric_limits<double>::epsilon())
45 const double max_delta(std::max(delta[1], delta[2]));
46 min_point[0] -= max_delta * 0.5e-3;
47 max_point[0] += max_delta * 0.5e-3;
48 delta[0] = max_point[0] - min_point[0];
51 if (delta[1] < std::numeric_limits<double>::epsilon())
53 const double max_delta(std::max(delta[0], delta[2]));
54 min_point[1] -= max_delta * 0.5e-3;
55 max_point[1] += max_delta * 0.5e-3;
56 delta[1] = max_point[1] - min_point[1];
59 if (delta[2] < std::numeric_limits<double>::epsilon())
61 const double max_delta(std::max(delta[0], delta[1]));
62 min_point[2] -= max_delta * 0.5e-3;
63 max_point[2] += max_delta * 0.5e-3;
64 delta[2] = max_point[2] - min_point[2];
70 const std::size_t n_tris(sfc->getNumberOfTriangles());
71 const std::size_t n_tris_per_cell(5);
73 Eigen::Matrix<bool, 3, 1> dim =
74 delta.array() >= std::numeric_limits<double>::epsilon();
83 auto sc_ceil = [](
double v)
84 {
return static_cast<std::size_t
>(std::ceil(v)); };
89 sc_ceil(std::cbrt(n_tris * delta[0] * delta[0] /
90 (n_tris_per_cell * delta[1] * delta[2])));
92 sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0));
94 sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0));
99 _n_steps[0] = sc_ceil(std::sqrt(n_tris * delta[0] /
100 (n_tris_per_cell * delta[2])));
101 _n_steps[2] = sc_ceil(_n_steps[0] * delta[2] / delta[0]);
103 else if (dim[0] && dim[1])
105 _n_steps[0] = sc_ceil(std::sqrt(n_tris * delta[0] /
106 (n_tris_per_cell * delta[1])));
107 _n_steps[1] = sc_ceil(_n_steps[0] * delta[1] / delta[0]);
109 else if (dim[1] && dim[2])
111 _n_steps[1] = sc_ceil(std::sqrt(n_tris * delta[1] /
112 (n_tris_per_cell * delta[2])));
114 sc_ceil(n_tris * delta[2] / (n_tris_per_cell * delta[1]));
118 for (std::size_t k(0); k < 3; ++k)
123 sc_ceil(
static_cast<double>(n_tris) / n_tris_per_cell);
129 for (std::size_t k(0); k < 3; k++)
131 _step_sizes[k] = delta[k] / _n_steps[k];
132 if (delta[k] > std::numeric_limits<double>::epsilon())
134 _inverse_step_sizes[k] = 1.0 / _step_sizes[k];
138 _inverse_step_sizes[k] = 0;
142 _triangles_in_grid_box.resize(_n_steps[0] * _n_steps[1] * _n_steps[2]);
143 sortTrianglesInGridCells(sfc);
152 Point const& p0(*((*sfc)[l]->getPoint(0)));
153 Point const& p1(*((*sfc)[l]->getPoint(1)));
154 Point const& p2(*((*sfc)[l]->getPoint(2)));
158 "Sorting triangle {:d} [({:f},{:f},{:f}), ({:f},{:f},{:f}), "
159 "({:f},{:f},{:f}) into grid. Bounding box is [{:f}, {:f}] x "
160 "[{:f}, {:f}] x [{:f}, {:f}].",
161 l, p0[0], p0[1], p0[2], p1[0], p1[1], p1[2], p2[0], p2[1],
162 p2[2], min[0], max[0], min[1], max[1], min[2], max[2]);
170 std::optional<std::array<std::size_t, 3>
const> c_p0(
176 std::optional<std::array<std::size_t, 3>
const> c_p1(
182 std::optional<std::array<std::size_t, 3>
const> c_p2(
190 std::size_t
const i_min(
191 std::min(std::min((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0]));
192 std::size_t
const i_max(
193 std::max(std::max((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0]));
194 std::size_t
const j_min(
195 std::min(std::min((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1]));
196 std::size_t
const j_max(
197 std::max(std::max((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1]));
198 std::size_t
const k_min(
199 std::min(std::min((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2]));
200 std::size_t
const k_max(
201 std::max(std::max((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2]));
206 for (std::size_t i(i_min); i <= i_max; i++)
208 for (std::size_t j(j_min); j <= j_max; j++)
210 for (std::size_t k(k_min); k <= k_max; k++)
213 .push_back(triangle);
225 std::array<std::size_t, 3> coords{
226 {
static_cast<std::size_t
>((p[0] - min_point[0]) *
228 static_cast<std::size_t
>((p[1] - min_point[1]) *
230 static_cast<std::size_t
>((p[2] - min_point[2]) *
237 "Computed indices ({:d},{:d},{:d}), max grid cell indices "
241 return std::optional<std::array<std::size_t, 3>>();
243 return std::optional<std::array<std::size_t, 3>>(coords);
249 std::optional<std::array<std::size_t, 3>> optional_c(
255 std::array<std::size_t, 3>
c(optional_c.value());
257 std::size_t
const grid_cell_idx(
c[0] +
c[1] *
_n_steps[0] +
259 std::vector<Triangle const*>
const& triangles(
261 auto const it = std::find_if(triangles.begin(), triangles.end(),
262 [eps, pnt](
auto const* triangle)
263 { return triangle->containsPoint(pnt, eps); });
264 return it != triangles.end();
void DBUG(char const *fmt, Args const &... args)
Definition of the Point3d class.
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
void sortTrianglesInGridCells(GeoLib::Surface const *const sfc)
std::vector< std::vector< GeoLib::Triangle const * > > _triangles_in_grid_box
bool isPointInSurface(MathLib::Point3d const &pnt, double eps=std::numeric_limits< double >::epsilon()) const
std::array< std::size_t, 3 > _n_steps
std::array< double, 3 > _inverse_step_sizes
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...
std::size_t getNumberOfTriangles() const
Class Triangle consists of a reference to a point vector and a vector that stores the indices in the ...
const Point * getPoint(std::size_t i) const
const access operator to access the i-th triangle Point