Revision a57c08dd8891d9efae6fbbad9a343f7149da3aa1 authored by Wenqing Wang on 04 June 2021, 13:34:00 UTC, committed by Wenqing Wang on 04 June 2021, 13:34:00 UTC
coordinate system
1 parent ddd427b
QuadTree.h
/**
* \file
* \author Thomas Fischer
* \date 2010-11-09
* \brief Definition of the QuadTree class.
*
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#pragma once
#include <cassert>
#include <limits>
#include <list>
#include <utility>
#include "BaseLib/Logging.h"
namespace GeoLib
{
/**
* A quadtree is a rooted tree in which every internal
* node has four children. Every node corresponds to a square.
* (see Computational Geometry - Algorithms and Applications [Mark de Berg,
* Otfried Cheong, Marc van Kreveld, Mark Overmars] - Chapter 14)
*
* One can instantiate the class template with a point type and a value for
* the maximal number of points per node. The point-type have to provide
* the access to its coordinates via operator[] and for debugging
* purposes operator<<)
*/
template <typename POINT> class QuadTree
{
public:
enum class Quadrant {
NE = 0, //!< north east
NW, //!< north west
SW, //!< south west
SE //!< south east
};
/**
* This is the constructor for class QuadTree. It takes two points
* (lower left and the upper right points).
* @param ll lower left point of the square
* @param ur upper right point of the square
* @param max_points_per_leaf maximum number of points per leaf
*/
QuadTree(POINT ll, POINT ur, std::size_t max_points_per_leaf)
: _father(nullptr),
_ll(std::move(ll)),
_ur(std::move(ur)),
_max_points_per_leaf(max_points_per_leaf)
{
assert (_max_points_per_leaf > 0);
if ((_ur[0] - _ll[0]) > (_ur[1] - _ll[1]))
{
_ur[1] = _ll[1] + _ur[0] - _ll[0];
}
else
{
_ur[0] = _ll[0] + _ur[1] - _ll[1];
}
DBUG(
"QuadTree(): lower left: ({:f},{:f},{:f}), upper right: "
"({:f},{:f},{:f}), depth: {:d}",
_ll[0],
_ll[1],
_ll[2],
_ur[0],
_ur[1],
_ur[2],
_depth);
}
/**
* destructor
*/
~QuadTree()
{
for (auto& child : _children)
{
delete child;
}
}
/**
* This method adds the given point to the quadtree. If necessary,
* the quadtree will be extended.
* @param pnt the point
* @return If the point can be inserted the method returns true, else false.
*/
bool addPoint (POINT const* pnt)
{
if ((*pnt)[0] < _ll[0])
{
return false;
}
if ((*pnt)[0] >= _ur[0])
{
return false;
}
if ((*pnt)[1] < _ll[1])
{
return false;
}
if ((*pnt)[1] >= _ur[1])
{
return false;
}
if (!_is_leaf) {
for (auto& child : _children)
{
if (child->addPoint(pnt))
{
return true;
}
}
return false;
}
// check if point is already in quadtree
bool pnt_in_quadtree (false);
for (std::size_t k(0); k < _pnts.size() && !pnt_in_quadtree; k++) {
const double v0((*(_pnts[k]))[0] - (*pnt)[0]);
const double v1((*(_pnts[k]))[1] - (*pnt)[1]);
const double sqr_dist (v0*v0 + v1*v1);
if (sqr_dist < std::numeric_limits<double>::epsilon())
{
pnt_in_quadtree = true;
}
}
if (!pnt_in_quadtree)
{
_pnts.push_back (pnt);
}
else
{
return false;
}
if (_pnts.size() > _max_points_per_leaf)
{
splitNode();
}
return true;
}
/**
* This method balances the quadtree, i.e., it will be inserted nodes
* such that the depth between neighbored leafs is at most one. If you want
* to create a mesh (for instance with GMSH) you can use this method to
* improve the mesh quality. The balance method should be used after
* inserting all points.
*/
void balance ()
{
std::list<QuadTree<POINT>*> leaf_list;
getLeafs (leaf_list);
while (!leaf_list.empty())
{
QuadTree<POINT>* node (leaf_list.front());
leaf_list.pop_front ();
if (node->isLeaf())
{
if (needToRefine (node))
{
node->splitNode ();
leaf_list.push_back (node->getChild(Quadrant::NE));
leaf_list.push_back (node->getChild(Quadrant::NW));
leaf_list.push_back (node->getChild(Quadrant::SW));
leaf_list.push_back (node->getChild(Quadrant::SE));
// check if north neighbor has to be refined
QuadTree<POINT>* north_neighbor (node->getNorthNeighbor());
if (north_neighbor != nullptr)
{
if (north_neighbor->getDepth() < node->getDepth())
{
if (north_neighbor->isLeaf())
{
leaf_list.push_back(north_neighbor);
}
}
}
// check if west neighbor has to be refined
QuadTree<POINT>* west_neighbor (node->getWestNeighbor());
if (west_neighbor != nullptr)
{
if (west_neighbor->getDepth() < node->getDepth())
{
if (west_neighbor->isLeaf())
{
leaf_list.push_back(west_neighbor);
}
}
}
// check if south neighbor has to be refined
QuadTree<POINT>* south_neighbor (node->getSouthNeighbor());
if (south_neighbor != nullptr)
{
if (south_neighbor->getDepth() < node->getDepth())
{
if (south_neighbor->isLeaf())
{
leaf_list.push_back(south_neighbor);
}
}
}
// check if east neighbor has to be refined
QuadTree<POINT>* east_neighbor (node->getEastNeighbor());
if (east_neighbor != nullptr)
{
if (east_neighbor->getDepth() < node->getDepth())
{
if (east_neighbor->isLeaf())
{
leaf_list.push_back(east_neighbor);
}
}
}
}
}
}
}
/**
* add all leafs of the quadtree to the list
* @param leaf_list list of leafs
*/
void getLeafs (std::list<QuadTree<POINT>*>& leaf_list)
{
if (_is_leaf)
{
leaf_list.push_back (this);
}
else
{
for (auto& child : _children)
{
child->getLeafs(leaf_list);
}
}
}
const std::vector<POINT const*>& getPoints () const { return _pnts; }
void getSquarePoints (POINT& ll, POINT& ur) const
{
ll = _ll;
ur = _ur;
}
void getLeaf (const POINT& pnt, POINT& ll, POINT& ur)
{
if (this->isLeaf())
{
ll = _ll;
ur = _ur;
}
else
{
if (pnt[0] <= 0.5 * (_ur[0] + _ll[0])) // WEST
{
if (pnt[1] <= 0.5 * (_ur[1] + _ll[1]))
{ // SOUTH
_children[static_cast<int>(Quadrant::SW)]->getLeaf (pnt, ll, ur);
}
else
{ // NORTH
_children[static_cast<int>(Quadrant::NW)]->getLeaf(
pnt, ll, ur);
}
}
else // EAST
{
if (pnt[1] <= 0.5 * (_ur[1] + _ll[1]))
{ // SOUTH
_children[static_cast<int>(Quadrant::SE)]->getLeaf (pnt, ll, ur);
}
else
{ // NORTH
_children[static_cast<int>(Quadrant::NE)]->getLeaf(
pnt, ll, ur);
}
}
}
}
QuadTree<POINT> const* getFather() const { return _father; }
QuadTree<POINT> const* getChild (Quadrant quadrant) const
{
return _children[quadrant];
}
/**
* Method calculates the maximum depth of the QuadTree instance. It is used within
* the method GMSHAdaptiveMeshDensity::getSteinerPoints().
* @param max_depth (input/output) at the entry max_depth contains the maximum depth up to now
*/
void getMaxDepth (std::size_t &max_depth) const
{
if (max_depth < _depth)
{
max_depth = _depth;
}
for (auto& child : _children)
{
if (child)
{
child->getMaxDepth(max_depth);
}
}
}
/**
* Method returns the depth of the current QuadTree node.
* @return the depth of the current QuadTree node
*/
std::size_t getDepth () const { return _depth; }
private:
QuadTree<POINT>* getChild (Quadrant quadrant)
{
return _children[static_cast<int>(quadrant)];
}
bool isLeaf () const { return _is_leaf; }
bool isChild (QuadTree<POINT> const* const tree, Quadrant quadrant) const
{
return _children[static_cast<int>(quadrant)] == tree;
}
QuadTree<POINT>* getNorthNeighbor () const
{
if (this->_father == nullptr)
{ // root of QuadTree
return nullptr;
}
if (this->_father->isChild(this, Quadrant::SW))
{
return this->_father->getChild(Quadrant::NW);
}
if (this->_father->isChild(this, Quadrant::SE))
{
return this->_father->getChild(Quadrant::NE);
}
QuadTree<POINT>* north_neighbor (this->_father->getNorthNeighbor ());
if (north_neighbor == nullptr)
{
return nullptr;
}
if (north_neighbor->isLeaf())
{
return north_neighbor;
}
if (this->_father->isChild(this, Quadrant::NW))
{
return north_neighbor->getChild(Quadrant::SW);
}
return north_neighbor->getChild(Quadrant::SE);
}
QuadTree<POINT>* getSouthNeighbor () const
{
if (this->_father == nullptr)
{ // root of QuadTree
return nullptr;
}
if (this->_father->isChild(this, Quadrant::NW))
{
return this->_father->getChild(Quadrant::SW);
}
if (this->_father->isChild(this, Quadrant::NE))
{
return this->_father->getChild(Quadrant::SE);
}
QuadTree<POINT>* south_neighbor (this->_father->getSouthNeighbor ());
if (south_neighbor == nullptr)
{
return nullptr;
}
if (south_neighbor->isLeaf())
{
return south_neighbor;
}
if (this->_father->isChild(this, Quadrant::SW))
{
return south_neighbor->getChild(Quadrant::NW);
}
return south_neighbor->getChild(Quadrant::NE);
}
QuadTree<POINT>* getEastNeighbor () const
{
if (this->_father == nullptr)
{ // root of QuadTree
return nullptr;
}
if (this->_father->isChild(this, Quadrant::NW))
{
return this->_father->getChild(Quadrant::NE);
}
if (this->_father->isChild(this, Quadrant::SW))
{
return this->_father->getChild(Quadrant::SE);
}
QuadTree<POINT>* east_neighbor (this->_father->getEastNeighbor ());
if (east_neighbor == nullptr)
{
return nullptr;
}
if (east_neighbor->isLeaf())
{
return east_neighbor;
}
if (this->_father->isChild(this, Quadrant::SE))
{
return east_neighbor->getChild(Quadrant::SW);
}
return east_neighbor->getChild(Quadrant::NW);
}
QuadTree<POINT>* getWestNeighbor () const
{
if (this->_father == nullptr)
{ // root of QuadTree
return nullptr;
}
if (this->_father->isChild(this, Quadrant::NE))
{
return this->_father->getChild(Quadrant::NW);
}
if (this->_father->isChild(this, Quadrant::SE))
{
return this->_father->getChild(Quadrant::SW);
}
QuadTree<POINT>* west_neighbor (this->_father->getWestNeighbor ());
if (west_neighbor == nullptr)
{
return nullptr;
}
if (west_neighbor->isLeaf())
{
return west_neighbor;
}
if (this->_father->isChild(this, Quadrant::SW))
{
return west_neighbor->getChild(Quadrant::SE);
}
return west_neighbor->getChild(Quadrant::NE);
}
/**
* private constructor
* @param ll lower left point
* @param ur upper right point
* @param father father in the tree
* @param depth depth of the node
* @param max_points_per_leaf maximum number of points per leaf
* @return
*/
QuadTree(POINT ll,
POINT ur,
QuadTree* father,
std::size_t depth,
std::size_t max_points_per_leaf)
: _father(father),
_ll(std::move(ll)),
_ur(std::move(ur)),
_depth(depth),
_max_points_per_leaf(max_points_per_leaf)
{
}
void splitNode ()
{
// create children
POINT mid_point(_ll);
mid_point[0] += (_ur[0] - _ll[0]) / 2.0;
mid_point[1] += (_ur[1] - _ll[1]) / 2.0;
assert(_children[0] == nullptr);
_children[0] = new QuadTree<POINT> (mid_point, _ur, this, _depth + 1, _max_points_per_leaf); // north east
POINT h_ll(mid_point), h_ur(mid_point);
h_ll[0] = _ll[0];
h_ur[1] = _ur[1];
assert(_children[1] == nullptr);
_children[1] = new QuadTree<POINT> (h_ll, h_ur, this, _depth + 1, _max_points_per_leaf); // north west
assert(_children[2] == nullptr);
_children[2] = new QuadTree<POINT> (_ll, mid_point, this, _depth + 1, _max_points_per_leaf); // south west
h_ll = _ll;
h_ll[0] = mid_point[0];
h_ur = _ur;
h_ur[1] = mid_point[1];
assert(_children[3] == nullptr);
_children[3] = new QuadTree<POINT> (h_ll, h_ur, this, _depth + 1, _max_points_per_leaf); // south east
// distribute points to sub quadtrees
for (std::size_t j(0); j < _pnts.size(); j++) {
bool nfound(true);
for (std::size_t k(0); k < 4 && nfound; k++)
{
if (_children[k]->addPoint(_pnts[j]))
{
nfound = false;
}
}
}
_pnts.clear();
_is_leaf = false;
}
static bool needToRefine (QuadTree<POINT>* node)
{
QuadTree<POINT>* north_neighbor (node->getNorthNeighbor ());
if (north_neighbor != nullptr)
{
if (north_neighbor->getDepth() == node->getDepth())
{
if (!north_neighbor->isLeaf ())
{
if (!(north_neighbor->getChild(Quadrant::SW))->isLeaf())
{
return true;
}
if (!(north_neighbor->getChild(Quadrant::SE))->isLeaf())
{
return true;
}
}
}
}
QuadTree<POINT>* west_neighbor (node->getWestNeighbor ());
if (west_neighbor != nullptr)
{
if (west_neighbor->getDepth() == node->getDepth())
{
if (!west_neighbor->isLeaf ())
{
if (!(west_neighbor->getChild(Quadrant::SE))->isLeaf())
{
return true;
}
if (!(west_neighbor->getChild(Quadrant::NE))->isLeaf())
{
return true;
}
}
}
}
QuadTree<POINT>* south_neighbor (node->getSouthNeighbor ());
if (south_neighbor != nullptr)
{
if (south_neighbor->getDepth() == node->getDepth())
{
if (!south_neighbor->isLeaf())
{
if (!(south_neighbor->getChild(Quadrant::NE))->isLeaf())
{
return true;
}
if (!(south_neighbor->getChild(Quadrant::NW))->isLeaf())
{
return true;
}
}
}
}
QuadTree<POINT>* east_neighbor (node->getEastNeighbor ());
if (east_neighbor != nullptr)
{
if (east_neighbor->getDepth() == node->getDepth())
{
if (!east_neighbor->isLeaf ())
{
if (!(east_neighbor->getChild(Quadrant::NW))->isLeaf())
{
return true;
}
if (!(east_neighbor->getChild(Quadrant::SW))->isLeaf())
{
return true;
}
}
}
}
return false;
}
QuadTree<POINT>* _father;
/**
* children are sorted:
* _children[0] is north east child
* _children[1] is north west child
* _children[2] is south west child
* _children[3] is south east child
*/
std::array<QuadTree<POINT>*, 4> _children = {nullptr, nullptr, nullptr,
nullptr};
/**
* lower left point of the square
*/
POINT _ll;
/**
* upper right point of the square
*/
POINT _ur;
std::size_t _depth = 0;
std::vector<POINT const*> _pnts;
bool _is_leaf = true;
/**
* maximum number of points per leaf
*/
const std::size_t _max_points_per_leaf;
};
} // namespace GeoLib
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...