OGS 6.2.0-97-g4a610c866
SurfaceGrid.cpp
Go to the documentation of this file.
1 /*
2  * \date 2012-09-22
3  * \brief Definition of the SurfaceGrid class.
4  *
5  * \copyright
6  * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
7  * Distributed under a Modified BSD License.
8  * See accompanying file LICENSE.txt or
9  */
10 
11 #include "SurfaceGrid.h"
12 
13 #include <algorithm>
14 #include <bitset>
15 
16 #include <logog/include/logog.hpp>
17 
18 #include "BaseLib/Error.h"
19 
20 #include "MathLib/Point3d.h"
21 
22 #include "Surface.h"
23 #include "Triangle.h"
24 
25 namespace GeoLib {
26 
28  AABB(sfc->getAABB()), _n_steps({{1,1,1}})
29 {
30  // enlarge the bounding, such that the points with maximal coordinates
31  // fits into the grid
32  for (std::size_t k(0); k<3; ++k) {
33  _max_pnt[k] += std::abs(_max_pnt[k]) * 1e-6;
34  if (std::abs(_max_pnt[k]) < std::numeric_limits<double>::epsilon()) {
35  _max_pnt[k] = (_max_pnt[k] - _min_pnt[k]) * (1.0 + 1e-6);
36  }
37  }
38 
39  std::array<double, 3> delta{{ _max_pnt[0] - _min_pnt[0],
40  _max_pnt[1] - _min_pnt[1], _max_pnt[2] - _min_pnt[2] }};
41 
42  if (delta[0] < std::numeric_limits<double>::epsilon()) {
43  const double max_delta(std::max(delta[1], delta[2]));
44  _min_pnt[0] -= max_delta * 0.5e-3;
45  _max_pnt[0] += max_delta * 0.5e-3;
46  delta[0] = _max_pnt[0] - _min_pnt[0];
47  }
48 
49  if (delta[1] < std::numeric_limits<double>::epsilon()) {
50  const double max_delta(std::max(delta[0], delta[2]));
51  _min_pnt[1] -= max_delta * 0.5e-3;
52  _max_pnt[1] += max_delta * 0.5e-3;
53  delta[1] = _max_pnt[1] - _min_pnt[1];
54  }
55 
56  if (delta[2] < std::numeric_limits<double>::epsilon()) {
57  const double max_delta(std::max(delta[0], delta[1]));
58  _min_pnt[2] -= max_delta * 0.5e-3;
59  _max_pnt[2] += max_delta * 0.5e-3;
60  delta[2] = _max_pnt[2] - _min_pnt[2];
61  }
62 
63  const std::size_t n_tris(sfc->getNumberOfTriangles());
64  const std::size_t n_tris_per_cell(5);
65 
66  std::bitset<3> dim; // all bits are set to zero.
67  for (std::size_t k(0); k < 3; ++k)
68  {
69  if (std::abs(delta[k]) >= std::numeric_limits<double>::epsilon())
70  {
71  dim[k] = true;
72  }
73  }
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] = ceil(pow(n_tris*delta[0]*delta[0]/(n_tris_per_cell*delta[1]*delta[2]), 1/3.)));
78  // _n_steps[1] = _n_steps[0] * delta[1]/delta[0],
79  // _n_steps[2] = _n_steps[0] * delta[2]/delta[0]
80  auto sc_ceil = [](double v){
81  return static_cast<std::size_t>(std::ceil(v));
82  };
83  switch (dim.count()) {
84  case 3: // 3d case
85  _n_steps[0] = sc_ceil(std::cbrt(
86  n_tris*delta[0]*delta[0]/(n_tris_per_cell*delta[1]*delta[2])));
87  _n_steps[1] = sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0));
88  _n_steps[2] = sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0));
89  break;
90  case 2: // 2d cases
91  if (dim[0] && dim[2]) { // 2d case: xz plane, y = const
92  _n_steps[0] = sc_ceil(std::sqrt(n_tris*delta[0]/(n_tris_per_cell*delta[2])));
93  _n_steps[2] = sc_ceil(_n_steps[0]*delta[2]/delta[0]);
94  }
95  else if (dim[0] && dim[1]) { // 2d case: xy plane, z = const
96  _n_steps[0] = sc_ceil(std::sqrt(n_tris*delta[0]/(n_tris_per_cell*delta[1])));
97  _n_steps[1] = sc_ceil(_n_steps[0] * delta[1]/delta[0]);
98  }
99  else if (dim[1] && dim[2]) { // 2d case: yz plane, x = const
100  _n_steps[1] = sc_ceil(std::sqrt(n_tris*delta[1]/(n_tris_per_cell*delta[2])));
101  _n_steps[2] = sc_ceil(n_tris * delta[2] / (n_tris_per_cell*delta[1]));
102  }
103  break;
104  case 1: // 1d cases
105  for (std::size_t k(0); k<3; ++k) {
106  if (dim[k]) {
107  _n_steps[k] = sc_ceil(static_cast<double>(n_tris)/n_tris_per_cell);
108  }
109  }
110  }
111 
112  // some frequently used expressions to fill the grid vectors
113  for (std::size_t k(0); k<3; k++) {
114  _step_sizes[k] = delta[k] / _n_steps[k];
115  _inverse_step_sizes[k] = 1.0 / _step_sizes[k];
116  }
117 
120 }
121 
123 {
124  for (std::size_t l(0); l<sfc->getNumberOfTriangles(); l++) {
125  if (! sortTriangleInGridCells((*sfc)[l])) {
126  Point const& p0(*((*sfc)[l]->getPoint(0)));
127  Point const& p1(*((*sfc)[l]->getPoint(1)));
128  Point const& p2(*((*sfc)[l]->getPoint(2)));
129  ERR("Sorting triangle %d [(%f,%f,%f), (%f,%f,%f), (%f,%f,%f) into "
130  "grid.",
131  l, p0[0], p0[1], p0[2], p1[0], p1[1], p1[2], p2[0], p2[1], p2[2]
132  );
133  OGS_FATAL("");
134  }
135  }
136 }
137 
139 {
140  // compute grid coordinates for each triangle point
141  boost::optional<std::array<std::size_t, 3> const> c_p0(
142  getGridCellCoordinates(*(triangle->getPoint(0))));
143  if (!c_p0)
144  {
145  return false;
146  }
147  boost::optional<std::array<std::size_t, 3> const> c_p1(
148  getGridCellCoordinates(*(triangle->getPoint(1))));
149  if (!c_p1)
150  {
151  return false;
152  }
153  boost::optional<std::array<std::size_t, 3> const> c_p2(
154  getGridCellCoordinates(*(triangle->getPoint(2))));
155  if (!c_p2)
156  {
157  return false;
158  }
159 
160  // determine interval in grid (grid cells the triangle will be inserted)
161  std::size_t const i_min(std::min(std::min((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0]));
162  std::size_t const i_max(std::max(std::max((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0]));
163  std::size_t const j_min(std::min(std::min((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1]));
164  std::size_t const j_max(std::max(std::max((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1]));
165  std::size_t const k_min(std::min(std::min((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2]));
166  std::size_t const k_max(std::max(std::max((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2]));
167 
168  const std::size_t n_plane(_n_steps[0]*_n_steps[1]);
169 
170  // insert the triangle into the grid cells
171  for (std::size_t i(i_min); i<=i_max; i++) {
172  for (std::size_t j(j_min); j<=j_max; j++) {
173  for (std::size_t k(k_min); k<=k_max; k++) {
174  _triangles_in_grid_box[i+j*_n_steps[0]+k*n_plane].push_back(triangle);
175  }
176  }
177  }
178 
179  return true;
180 }
181 
182 boost::optional<std::array<std::size_t, 3>>
184 {
185  std::array<std::size_t, 3> coords{{
186  static_cast<std::size_t>((p[0]-_min_pnt[0]) * _inverse_step_sizes[0]),
187  static_cast<std::size_t>((p[1]-_min_pnt[1]) * _inverse_step_sizes[1]),
188  static_cast<std::size_t>((p[2]-_min_pnt[2]) * _inverse_step_sizes[2])
189  }};
190 
191  if (coords[0]>=_n_steps[0] || coords[1]>=_n_steps[1] || coords[2]>=_n_steps[2]) {
192  DBUG("Computed indices (%d,%d,%d), max grid cell indices (%d,%d,%d)",
193  coords[0], coords[1], coords[2], _n_steps[0], _n_steps[1], _n_steps[2]);
194  return boost::optional<std::array<std::size_t, 3>>();
195  }
196  return boost::optional<std::array<std::size_t, 3>>(coords);
197 }
198 
199 bool SurfaceGrid::isPointInSurface(MathLib::Point3d const& p, double eps) const
200 {
201  boost::optional<std::array<std::size_t,3>> optional_c(getGridCellCoordinates(p));
202  if (!optional_c) {
203  return false;
204  }
205  std::array<std::size_t,3> c(optional_c.get());
206 
207  std::size_t const grid_cell_idx(c[0]+c[1]*_n_steps[0]+c[2]*_n_steps[0]*_n_steps[1]);
208  std::vector<Triangle const*> const& triangles(_triangles_in_grid_box[grid_cell_idx]);
209  bool nfound(true);
210  const std::size_t n_triangles(triangles.size());
211  for (std::size_t k(0); k < n_triangles && nfound; k++) {
212  if (triangles[k]->containsPoint(p, eps)) {
213  nfound = false;
214  }
215  }
216 
217  return !nfound;
218 }
219 
220 } // end namespace GeoLib
bool containsPoint(T const &pnt, double eps) const
Definition: AABB.h:130
std::array< std::size_t, 3 > _n_steps
Definition: SurfaceGrid.h:42
std::vector< std::vector< GeoLib::Triangle const * > > _triangles_in_grid_box
Definition: SurfaceGrid.h:43
boost::optional< std::array< std::size_t, 3 > > getGridCellCoordinates(MathLib::Point3d const &p) const
std::array< double, 3 > _inverse_step_sizes
Definition: SurfaceGrid.h:41
MathLib::Point3d _min_pnt
Definition: AABB.h:189
Definition of the Point3d class.
std::size_t getNumberOfTriangles() const
Definition: Surface.cpp:90
bool sortTriangleInGridCells(GeoLib::Triangle const *const triangle)
Definition of the GEOObjects class.
Definition: BaseItem.h:20
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Definition: AABB.h:50
std::array< double, 3 > _step_sizes
Definition: SurfaceGrid.h:40
Class Triangle consists of a reference to a point vector and a vector that stores the indices in the ...
Definition: Triangle.h:26
A Surface is represented by Triangles. It consists of a reference to a vector of (pointers to) points...
Definition: Surface.h:34
SurfaceGrid(GeoLib::Surface const *const sfc)
Definition: SurfaceGrid.cpp:27
const Point * getPoint(std::size_t i) const
const access operator to access the i-th triangle Point
Definition: Triangle.h:50
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
void sortTrianglesInGridCells(GeoLib::Surface const *const surface)
bool isPointInSurface(MathLib::Point3d const &pnt, double eps=std::numeric_limits< double >::epsilon()) const
MathLib::Point3d _max_pnt
Definition: AABB.h:193