From bff0290de15f9cd9cfa5b3da28edf79cc24e1be4 Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Tue, 28 Aug 2012 15:18:27 +0200
Subject: [PATCH] changed implementation of class template GeoLib::Grid and
added a specialization for pointer types
---
GeoLib/Grid.h | 376 ++++++++++++++++++++++++++++++++++++++++++++++----
1 file changed, 351 insertions(+), 25 deletions(-)
diff --git a/GeoLib/Grid.h b/GeoLib/Grid.h
index b5c5f0f7cbe..3936e982621 100644
--- a/GeoLib/Grid.h
+++ b/GeoLib/Grid.h
@@ -32,18 +32,19 @@ class Grid : public GeoLib::AABB
public:
/**
* @brief The constructor of the grid object takes a vector of points or nodes. Furthermore the
- * use can specify the *average* maximum number of points per grid cell.
+ * user can specify the *average* maximum number of points per grid cell.
*
* The number of grid cells are computed with the following formula
* \f$\frac{n_{points}}{n_{cells}} \le n_{max\_per\_cell}\f$
*
* In order to limit the memory wasting the maximum number of points per grid cell
- * (in the average) should be a power of two.
+ * (in the average) should be a power of two (since std::vector objects resize itself
+ * with this step size).
*
* @param pnts (input) the points that are managed with the Grid
* @param max_num_per_grid_cell (input) max number per grid cell in the average (default 512)
*/
- Grid(std::vector<POINT*> const& pnts, size_t max_num_per_grid_cell = 512);
+ Grid(std::vector<POINT> const& pnts, size_t max_num_per_grid_cell = 512);
/**
* This is the destructor of the class. It deletes the internal data structures *not*
@@ -77,7 +78,7 @@ public:
* @param pnts (output) vector of vectors of points within grid cells that intersects
* the axis aligned cube
*/
- void getVecsOfGridCellsIntersectingCube(double const*const pnt, double half_len, std::vector<std::vector<POINT*> const*>& pnts) const;
+ void getVecsOfGridCellsIntersectingCube(double const*const pnt, double half_len, std::vector<std::vector<POINT> const*>& pnts) const;
#ifndef NDEBUG
/**
@@ -136,15 +137,14 @@ private:
std::vector<POINT*>* _grid_quad_to_node_map;
};
-
template <typename POINT>
-Grid<POINT>::Grid(std::vector<POINT*> const& pnts, size_t max_num_per_grid_cell) :
+Grid<POINT>::Grid(std::vector<POINT> const& pnts, size_t max_num_per_grid_cell) :
GeoLib::AABB(), _grid_quad_to_node_map(NULL)
{
// compute axis aligned bounding box
const size_t n_pnts(pnts.size());
for (size_t k(0); k < n_pnts; k++) {
- this->update(pnts[k]->getCoords());
+ this->update(pnts[k].getCoords());
}
double delta[3] = { 0.0, 0.0, 0.0 };
@@ -215,7 +215,7 @@ Grid<POINT>::Grid(std::vector<POINT*> const& pnts, size_t max_num_per_grid_cell)
// fill the grid vectors
for (size_t l(0); l < n_pnts; l++) {
- double const* const pnt(pnts[l]->getCoords());
+ double const* const pnt(pnts[l].getCoords());
const size_t i(static_cast<size_t> ((pnt[0] - _min_pnt[0]) * _inverse_step_sizes[0]));
const size_t j(static_cast<size_t> ((pnt[1] - _min_pnt[1]) * _inverse_step_sizes[1]));
const size_t k(static_cast<size_t> ((pnt[2] - _min_pnt[2]) * _inverse_step_sizes[2]));
@@ -224,7 +224,7 @@ Grid<POINT>::Grid(std::vector<POINT*> const& pnts, size_t max_num_per_grid_cell)
std::cout << "error computing indices " << std::endl;
}
- _grid_quad_to_node_map[i + j * _n_steps[0] + k * n_plane].push_back(pnts[l]);
+ _grid_quad_to_node_map[i + j * _n_steps[0] + k * n_plane].push_back(&pnts[l]);
}
#ifndef NDEBUG
@@ -309,30 +309,27 @@ POINT const* Grid<POINT>::getNearestPoint(double const*const pnt) const
return nearest_pnt;
}
-template <typename POINT>
-void Grid<POINT>::getVecsOfGridCellsIntersectingCube(double const*const pnt, double half_len, std::vector<std::vector<POINT*> const*>& pnts) const
+template<typename POINT>
+void Grid<POINT>::getVecsOfGridCellsIntersectingCube(double const* const pnt, double half_len,
+ std::vector<std::vector<POINT> const*>& pnts) const
{
- double tmp_pnt[3] = {pnt[0]-half_len, pnt[1]-half_len, pnt[2]-half_len}; // min
+ double tmp_pnt[3] = { pnt[0] - half_len, pnt[1] - half_len, pnt[2] - half_len }; // min
size_t min_coords[3];
getGridCoords(tmp_pnt, min_coords);
- tmp_pnt[0] = pnt[0]+half_len;
- tmp_pnt[1] = pnt[1]+half_len;
- tmp_pnt[2] = pnt[2]+half_len;
+ tmp_pnt[0] = pnt[0] + half_len;
+ tmp_pnt[1] = pnt[1] + half_len;
+ tmp_pnt[2] = pnt[2] + half_len;
size_t max_coords[3];
getGridCoords(tmp_pnt, max_coords);
size_t coords[3], steps0_x_steps1(_n_steps[0] * _n_steps[1]);
- for (coords[0] = min_coords[0]; coords[0] < max_coords[0]+1; coords[0]++) {
- for (coords[1] = min_coords[1]; coords[1] < max_coords[1]+1; coords[1]++) {
- const size_t coords0_p_coords1_x_steps0 (coords[0] + coords[1] * _n_steps[0]);
- for (coords[2] = min_coords[2]; coords[2] < max_coords[2]+1; coords[2]++) {
- pnts.push_back (&(_grid_quad_to_node_map[coords0_p_coords1_x_steps0 + coords[2] * steps0_x_steps1]));
-// // copy pnts
-// const size_t n_pnts(_grid_quad_to_node_map[coords0_p_coords1_x_steps0 + coords[2] * steps0_x_steps1].size());
-// for (size_t k(0); k<n_pnts; k++) {
-// pnts.push_back(_grid_quad_to_node_map[coords0_p_coords1_x_steps0 + coords[2] * steps0_x_steps1][k]);
-// }
+ for (coords[0] = min_coords[0]; coords[0] < max_coords[0] + 1; coords[0]++) {
+ for (coords[1] = min_coords[1]; coords[1] < max_coords[1] + 1; coords[1]++) {
+ const size_t coords0_p_coords1_x_steps0(coords[0] + coords[1] * _n_steps[0]);
+ for (coords[2] = min_coords[2]; coords[2] < max_coords[2] + 1; coords[2]++) {
+ pnts.push_back(&(_grid_quad_to_node_map[coords0_p_coords1_x_steps0 + coords[2]
+ * steps0_x_steps1]));
}
}
}
@@ -471,6 +468,335 @@ void Grid<POINT>::getPointCellBorderDistances(double const*const pnt, double dis
dists[2] = (_step_sizes[0] - dists[4]); // right
}
+/**
+ * this is a specialization of class template Grid for pointer types
+ */
+template <typename POINT>
+class Grid<POINT*> : public GeoLib::AABB
+{
+public:
+ Grid(std::vector<POINT> const& pnts, size_t max_num_per_grid_cell = 512) :
+ GeoLib::AABB(), _grid_quad_to_node_map(NULL)
+ {
+ // compute axis aligned bounding box
+ const size_t n_pnts(pnts.size());
+ for (size_t k(0); k < n_pnts; k++) {
+ this->update(pnts[k]->getCoords());
+ }
+
+ double delta[3] = { 0.0, 0.0, 0.0 };
+ for (size_t k(0); k < 3; k++) {
+ // make the bounding box a little bit bigger,
+ // such that the node with maximal coordinates fits into the grid
+ _max_pnt[k] *= (1.0 + 1e-6);
+ if (fabs(_max_pnt[k]) < std::numeric_limits<double>::epsilon()) {
+ _max_pnt[k] = (_max_pnt[k] - _min_pnt[k]) * (1.0 + 1e-6);
+ }
+ delta[k] = _max_pnt[k] - _min_pnt[k];
+ }
+
+ // *** condition: n_pnts / (_n_steps[0] * _n_steps[1] * _n_steps[2]) < max_num_per_grid_cell
+ // *** with _n_steps[1] = _n_steps[0] * delta[1]/delta[0], _n_steps[2] = _n_steps[0] * delta[2]/delta[0]
+ if (fabs(delta[1]) < std::numeric_limits<double>::epsilon() ||
+ fabs(delta[2]) < std::numeric_limits<double>::epsilon()) {
+ // 1d case y = z = 0
+ if (fabs(delta[1]) < std::numeric_limits<double>::epsilon() &&
+ fabs(delta[2]) < std::numeric_limits<double>::epsilon()) {
+ _n_steps[0] = static_cast<size_t> (ceil(n_pnts / (double) max_num_per_grid_cell));
+ _n_steps[1] = 1;
+ _n_steps[2] = 1;
+ } else {
+ // 1d case x = z = 0
+ if (fabs(delta[0]) < std::numeric_limits<double>::epsilon() &&
+ fabs(delta[2]) < std::numeric_limits<double>::epsilon()) {
+ _n_steps[0] = 1;
+ _n_steps[1] = static_cast<size_t> (ceil(n_pnts / (double) max_num_per_grid_cell));
+ _n_steps[2] = 1;
+ } else {
+ // 1d case x = y = 0
+ if (fabs(delta[0]) < std::numeric_limits<double>::epsilon() && fabs(delta[1])
+ < std::numeric_limits<double>::epsilon()) {
+ _n_steps[0] = 1;
+ _n_steps[1] = 1;
+ _n_steps[2] = static_cast<size_t> (ceil(n_pnts / (double) max_num_per_grid_cell));
+ } else {
+ // 2d case
+ if (fabs(delta[1]) < std::numeric_limits<double>::epsilon()) {
+ _n_steps[0] = static_cast<size_t> (ceil(sqrt(n_pnts * delta[0] / (max_num_per_grid_cell * delta[2]))));
+ _n_steps[1] = 1;
+ _n_steps[2] = static_cast<size_t> (ceil(_n_steps[0] * delta[2] / delta[0]));
+ } else {
+ _n_steps[0] = static_cast<size_t> (ceil(sqrt(n_pnts * delta[0] / (max_num_per_grid_cell * delta[1]))));
+ _n_steps[1] = static_cast<size_t> (ceil(_n_steps[0] * delta[1] / delta[0]));
+ _n_steps[2] = 1;
+ }
+ }
+ }
+ }
+ } else {
+ // 3d case
+ _n_steps[0] = static_cast<size_t> (ceil(pow(n_pnts * delta[0] * delta[0]
+ / (max_num_per_grid_cell * delta[1] * delta[2]), 1. / 3.)));
+ _n_steps[1] = static_cast<size_t> (ceil(_n_steps[0] * delta[1] / delta[0]));
+ _n_steps[2] = static_cast<size_t> (ceil(_n_steps[0] * delta[2] / delta[0]));
+ }
+
+ const size_t n_plane(_n_steps[0] * _n_steps[1]);
+ _grid_quad_to_node_map = new std::vector<POINT>[n_plane * _n_steps[2]];
+
+ // some frequently used expressions to fill the grid vectors
+ for (size_t k(0); k < 3; k++) {
+ _step_sizes[k] = delta[k] / _n_steps[k];
+ _inverse_step_sizes[k] = 1.0 / _step_sizes[k];
+ }
+
+ // fill the grid vectors
+ for (size_t l(0); l < n_pnts; l++) {
+ double const* const pnt(pnts[l]->getCoords());
+ const size_t i(static_cast<size_t> ((pnt[0] - _min_pnt[0]) * _inverse_step_sizes[0]));
+ const size_t j(static_cast<size_t> ((pnt[1] - _min_pnt[1]) * _inverse_step_sizes[1]));
+ const size_t k(static_cast<size_t> ((pnt[2] - _min_pnt[2]) * _inverse_step_sizes[2]));
+
+ if (i >= _n_steps[0] || j >= _n_steps[1] || k >= _n_steps[2]) {
+ std::cout << "error computing indices " << std::endl;
+ }
+
+ _grid_quad_to_node_map[i + j * _n_steps[0] + k * n_plane].push_back(pnts[l]);
+ }
+
+#ifndef NDEBUG
+ size_t pnts_cnt(0);
+ for (size_t k(0); k < n_plane * _n_steps[2]; k++)
+ pnts_cnt += _grid_quad_to_node_map[k].size();
+
+ assert(n_pnts==pnts_cnt);
+#endif
+ }
+
+ /**
+ * This is the destructor of the class. It deletes the internal data structures *not*
+ * including the pointers to the points.
+ */
+ virtual ~Grid()
+ {
+ delete [] _grid_quad_to_node_map;
+ }
+
+ /**
+ * The method calculates the grid cell the given point is belonging to, i.e.,
+ * the (internal) coordinates of the grid cell are computed. The method searches the actual
+ * grid cell and all its neighbors for the POINT object which has the smallest
+ * distance. A pointer to this object is returned.
+ *
+ * If there is not such a point, i.e., all the searched grid cells do not contain any
+ * POINT object a NULL pointer is returned.
+ *
+ * @param pnt a field that holds the coordinates of the point
+ * @return a point with the smallest distance within the grid cells that are outlined above or NULL
+ */
+ POINT const* getNearestPoint(double const*const pnt) const
+ {
+ size_t coords[3];
+ getGridCoords(pnt, coords);
+
+ double sqr_min_dist (MathLib::sqrDist(&_min_pnt, &_max_pnt));
+ POINT* nearest_pnt(NULL);
+
+ double dists[6];
+ getPointCellBorderDistances(pnt, dists, coords);
+
+ if (calcNearestPointInGridCell(pnt, coords, sqr_min_dist, nearest_pnt)) {
+ double min_dist(sqrt(sqr_min_dist));
+ if (dists[0] >= min_dist
+ && dists[1] >= min_dist
+ && dists[2] >= min_dist
+ && dists[3] >= min_dist
+ && dists[4] >= min_dist
+ && dists[5] >= min_dist) {
+ return nearest_pnt;
+ }
+ } else {
+ // search in all border cells for at least one neighbor
+ double sqr_min_dist_tmp;
+ POINT* nearest_pnt_tmp(NULL);
+ size_t offset(1);
+
+ while (nearest_pnt == NULL) {
+ size_t tmp_coords[3];
+ for (tmp_coords[0] = coords[0]-offset; tmp_coords[0]<coords[0]+offset; tmp_coords[0]++) {
+ for (tmp_coords[1] = coords[1]-offset; tmp_coords[1]<coords[1]+offset; tmp_coords[1]++) {
+ for (tmp_coords[2] = coords[2]-offset; tmp_coords[2]<coords[2]+offset; tmp_coords[2]++) {
+ // do not check the origin grid cell twice
+ if (!(tmp_coords[0] == coords[0] && tmp_coords[1] == coords[1] && tmp_coords[2] == coords[2])) {
+ // check if temporary grid cell coordinates are valid
+ if (tmp_coords[0] < _n_steps[0] && tmp_coords[1] < _n_steps[1] && tmp_coords[2] < _n_steps[2]) {
+ if (calcNearestPointInGridCell(pnt, tmp_coords, sqr_min_dist_tmp, nearest_pnt_tmp)) {
+ if (sqr_min_dist_tmp < sqr_min_dist) {
+ sqr_min_dist = sqr_min_dist_tmp;
+ nearest_pnt = nearest_pnt_tmp;
+ }
+ }
+ } // valid grid cell coordinates
+ } // same element
+ } // end k
+ } // end j
+ } // end i
+ offset++;
+ } // end while
+ } // end else
+
+ double len (sqrt(MathLib::sqrDist(pnt, nearest_pnt->getCoords())));
+ // search all other grid cells within the cube with the edge nodes
+ std::vector<std::vector<POINT*> const*> vecs_of_pnts;
+ getVecsOfGridCellsIntersectingCube(pnt, len, vecs_of_pnts);
+
+ const size_t n_vecs(vecs_of_pnts.size());
+ for (size_t j(0); j<n_vecs; j++) {
+ std::vector<POINT*> const& pnts(vecs_of_pnts[j]);
+ const size_t n_pnts(pnts.size());
+ for (size_t k(0); k<n_pnts; k++) {
+ const double sqr_dist (MathLib::sqrDist(pnt, pnts[k]->getCoords()));
+ if (sqr_dist < sqr_min_dist) {
+ sqr_min_dist = sqr_dist;
+ nearest_pnt = pnts[k];
+ }
+ }
+ }
+
+ return nearest_pnt;
+ }
+
+ /**
+ * Method fetches the vectors of all grid cells intersecting the axis aligned cube
+ * defined by its center and half edge length.
+ *
+ * @param pnt (input) the center point of the axis aligned cube
+ * @param half_len (input) half of the edge length of the axis aligned cube
+ * @param pnts (output) vector of vectors of points within grid cells that intersects
+ * the axis aligned cube
+ */
+ void getVecsOfGridCellsIntersectingCube(double const*const pnt, double half_len, std::vector<std::vector<POINT*> const*>& pnts) const
+ {
+ double tmp_pnt[3] = {pnt[0]-half_len, pnt[1]-half_len, pnt[2]-half_len}; // min
+ size_t min_coords[3];
+ getGridCoords(tmp_pnt, min_coords);
+
+ tmp_pnt[0] = pnt[0]+half_len;
+ tmp_pnt[1] = pnt[1]+half_len;
+ tmp_pnt[2] = pnt[2]+half_len;
+ size_t max_coords[3];
+ getGridCoords(tmp_pnt, max_coords);
+
+ size_t coords[3], steps0_x_steps1(_n_steps[0] * _n_steps[1]);
+ for (coords[0] = min_coords[0]; coords[0] < max_coords[0]+1; coords[0]++) {
+ for (coords[1] = min_coords[1]; coords[1] < max_coords[1]+1; coords[1]++) {
+ const size_t coords0_p_coords1_x_steps0 (coords[0] + coords[1] * _n_steps[0]);
+ for (coords[2] = min_coords[2]; coords[2] < max_coords[2]+1; coords[2]++) {
+ pnts.push_back (&(_grid_quad_to_node_map[coords0_p_coords1_x_steps0 + coords[2] * steps0_x_steps1]));
+ }
+ }
+ }
+ }
+
+#ifndef NDEBUG
+ /**
+ * Method creates a geometry for every mesh grid box. Additionally it
+ * creates one geometry containing all the box geometries.
+ * @param geo_obj
+ */
+ void createGridGeometry(GeoLib::GEOObjects* geo_obj) const;
+#endif
+
+private:
+ /**
+ * Method calculates the grid cell coordinates for the given point pnt. If
+ * the point is located outside of the bounding box the coordinates of the
+ * grid cell on the border that is nearest to given point will be returned.
+ * @param pnt (input) the coordinates of the point
+ * @param coords (output) the coordinates of the grid cell
+ */
+ inline void getGridCoords(double const*const pnt, size_t* coords) const
+ {
+ for (size_t k(0); k<3; k++) {
+ if (pnt[k] < _min_pnt[k]) {
+ coords[k] = 0;
+ } else {
+ if (pnt[k] > _max_pnt[k]) {
+ coords[k] = _n_steps[k]-1;
+ } else {
+ coords[k] = static_cast<size_t>((pnt[k]-_min_pnt[k]) * _inverse_step_sizes[k]);
+ }
+ }
+ }
+ }
+
+
+ /**
+ *
+ * point numbering of the grid cell is as follow
+ * 7 -------- 6
+ * /: /|
+ * / : / |
+ * / : / |
+ * / : / |
+ * 4 -------- 5 |
+ * | 3 ....|... 2
+ * | . | /
+ * | . | /
+ * | . | /
+ * |. |/
+ * 0 -------- 1
+ * the face numbering is as follow:
+ * face nodes
+ * 0 0,3,2,1 bottom
+ * 1 0,1,5,4 front
+ * 2 1,2,6,5 right
+ * 3 2,3,7,6 back
+ * 4 3,0,4,7 left
+ * 5 4,5,6,7 top
+ * @param pnt (input) coordinates of the point
+ * @param dists (output) squared distances of the point to the faces
+ * ordered in the same sequence as above described
+ * @param coords coordinates of the grid cell
+ */
+ void getPointCellBorderDistances(double const*const pnt, double dists[6], size_t const* const coords) const
+ {
+ dists[0] = (pnt[2] - _min_pnt[2] + coords[2]*_step_sizes[2]); // bottom
+ dists[5] = (_step_sizes[2] - dists[0]); // top
+
+ dists[1] = (pnt[1] - _min_pnt[1] + coords[1]*_step_sizes[1]); // front
+ dists[3] = (_step_sizes[1] - dists[1]); // back
+
+ dists[4] = (pnt[0] - _min_pnt[0] + coords[0]*_step_sizes[0]); // left
+ dists[2] = (_step_sizes[0] - dists[4]); // right
+ }
+
+ bool calcNearestPointInGridCell(double const* const pnt, size_t const* const coords,
+ double &sqr_min_dist, POINT* &nearest_pnt) const
+ {
+ const size_t grid_idx (coords[0] + coords[1] * _n_steps[0] + coords[2] * _n_steps[0] * _n_steps[1]);
+ std::vector<POINT*> const& pnts(_grid_quad_to_node_map[grid_idx]);
+ if (pnts.empty()) return false;
+
+ const size_t n_pnts(pnts.size());
+ sqr_min_dist = MathLib::sqrDist(pnts[0]->getCoords(), pnt);
+ nearest_pnt = pnts[0];
+ for (size_t i(1); i < n_pnts; i++) {
+ const double sqr_dist(MathLib::sqrDist(pnts[i]->getCoords(), pnt));
+ if (sqr_dist < sqr_min_dist) {
+ sqr_min_dist = sqr_dist;
+ nearest_pnt = pnts[i];
+ }
+ }
+ return true;
+ };
+
+ double _step_sizes[3];
+ double _inverse_step_sizes[3];
+ size_t _n_steps[3];
+ std::vector<POINT>* _grid_quad_to_node_map;
+};
+
} // end namespace GeoLib
#endif /* MESHGRID_H_ */
--
GitLab