All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Octree.h
Go to the documentation of this file.
1 /*
2  * Octree.h
3  *
4  * Created on: Apr 21, 2016
5  * Author: Michael Wegner
6  */
7 
8 #ifndef NETWORKIT_CPP_VIZ_OCTREE_H_
9 #define NETWORKIT_CPP_VIZ_OCTREE_H_
10 
11 #include <vector>
12 
13 #include "Point.h"
14 #include "../algebraic/Vector.h"
15 
16 #include "../auxiliary/Log.h"
17 
18 namespace NetworKit {
19 
23 template<typename T>
24 struct BoundingBox {
25 public:
29  BoundingBox() : center(Point<T>()), sideLength(0), halfSideLength(0), sqSideLength(0), dimension(0) {}
30 
36  BoundingBox(const Point<T>& center, const T sideLength) : center(center), sideLength(sideLength), halfSideLength(sideLength/2.0), sqSideLength(sideLength*sideLength), dimension(center.getDimensions()) {}
37 
41  BoundingBox(const BoundingBox<T>& other) : center(other.center), sideLength(other.sideLength), halfSideLength(other.halfSideLength), sqSideLength(other.sqSideLength), dimension(other.dimension) {}
42 
47  void setCenter(const Point<T>& center) {
48  this->center = center;
49  dimension = center.getDimensions();
50  }
51 
55  inline Point<T>& getCenter() {
56  return center;
57  }
58 
63  void setSideLength(T sideLength) {
64  this->sideLength = sideLength;
65  this->halfSideLength = sideLength/2.0;
66  this->sqSideLength = sideLength * sideLength;
67  }
68 
72  inline T getSideLength() const {
73  return sideLength;
74  }
75 
79  inline T getHalfSideLength() const {
80  return halfSideLength;
81  }
82 
86  inline T getSqSideLength() const {
87  return sqSideLength;
88  }
89 
93  bool contains(const Point<T>& point) const {
94  for (index d = 0; d < dimension; ++d) {
95  if (center[d] - halfSideLength > point[d] || point[d] > center[d] + halfSideLength) {
96  return false;
97  }
98  }
99 
100  return true;
101  }
102 
103 private:
104  Point<T> center;
105  T sideLength;
106  T halfSideLength;
107  T sqSideLength;
108  count dimension;
109 };
110 
114 template<typename T>
115 struct OctreeNode {
118  std::vector<OctreeNode> children;
120 
121  OctreeNode() : weight(0), centerOfMass({0,0}), children({}), bBox(BoundingBox<T>()) {}
122  OctreeNode(BoundingBox<T>& bBox) : weight(0), centerOfMass(Point<T>(bBox.getCenter().getDimensions())), children({}), bBox(bBox) {}
123 
127  inline bool isLeaf() const {
128  return children.size() == 0;
129  }
130 
134  inline bool isEmpty() const {
135  return weight == 0;
136  }
137 
141  inline bool contains(const Point<T>& point) const {
142  return bBox.contains(point);
143  }
144 
149  if (!isLeaf()) {
150  centerOfMass.scale(1.0/(double) weight);
151 
152  // remove empty childs
153  children.erase(std::remove_if(children.begin(), children.end(), [&](OctreeNode<T>& child){return child.isEmpty();}), children.end());
154 
155  for (auto &child : children) {
156  child.computeCenterOfMass();
157  }
158  }
159  }
160 
161  void compress() {
162  if (!isLeaf()) {
163  count numNonEmptyChilren = 0;
164  while (numNonEmptyChilren <= 1) {
165  numNonEmptyChilren = 0;
166  OctreeNode<T> lastChild;
167  for (auto &child : children) {
168  if (!child.isEmpty()) {
169  numNonEmptyChilren++;
170  lastChild = child;
171  }
172  }
173 
174  if (numNonEmptyChilren == 1) { // compress
175  children = lastChild.children;
176  bBox = lastChild.bBox;
177  }
178  }
179 
180  for (auto &child : children) {
181  child.compress();
182  }
183  }
184  }
185 
189  std::string toString() {
190  std::string str;
191  str += bBox.getCenter().toString() + " sL=" + std::to_string(bBox.getSideLength()) + "w=" + std::to_string(weight);
192  str += "(";
193  for (auto &child : children) {
194  str += "[ ";
195  str += child.toString();
196  str += " ]";
197  }
198  str += ")";
199  return str;
200  }
201 
205  void split(count dimensions, count numChildren) {
206  children = std::vector<OctreeNode<T>>(numChildren, OctreeNode<T>(bBox));
207  for (index i = 0; i < numChildren; ++i) { // 0-bit => center - halfSideLength, 1-bit => center + halfSideLength, least-significant bit is lowest dimension
208  children[i].bBox.setSideLength(bBox.getHalfSideLength());
209  for (index d = 0; d < dimensions; ++d) {
210  if (((i & ~(~static_cast<index>(0) << (d+1))) >> d) == 0) { // 0-bit
211  children[i].bBox.getCenter()[d] -= children[i].bBox.getHalfSideLength();
212  } else {
213  children[i].bBox.getCenter()[d] += children[i].bBox.getHalfSideLength();
214  }
215  }
216  }
217  }
218 
225  void addPoint(const Point<T>& point, count dimensions, count numChildrenPerNode) {
226  if (weight == 0) { // empty leaf
227  weight++; // we add a point
228  centerOfMass = point; // center of mass of a single point is the point
229  } else {
230  if (isLeaf()) { // split the leaf!
231  if (point.distance(centerOfMass) < 1e-3) {
232  centerOfMass += point;
233  weight++;
234  return;
235  }
236  split(dimensions, numChildrenPerNode);
237  for (auto &child : children) { // add leaf point to one of this node's new children
238  if (child.contains(centerOfMass)) {
239  child.addPoint(centerOfMass, dimensions, numChildrenPerNode);
240  break;
241  }
242  }
243  }
244 
245  for (auto &child : children) {
246  if (child.contains(point)) {
247  child.addPoint(point, dimensions, numChildrenPerNode);
248  break;
249  }
250  }
251 
252  weight++;
253  centerOfMass += point;
254  }
255  }
256 };
257 
258 
264 template<typename T>
265 class Octree {
266 public:
270  Octree() = default;
271 
276  Octree(const std::vector<Vector>& points);
277 
282  void recomputeTree(const std::vector<Vector> &points);
283 
284  inline std::vector<std::pair<count, Point<T>>> approximateDistance(const Point<T>& p, const double theta) const {
285  return approximateDistance(root, p, theta);
286  }
287 
288  inline void approximateDistance(const Point<T>& p, const double theta, std::vector<std::pair<count, Point<T>>>& result) const {
289  approximateDistance(root, p, theta, result);
290  }
291 
292  template<typename L>
293  inline void approximateDistance(const Point<T>& p, const double theta, L& handle) const {
294  approximateDistance(root, p, theta*theta, handle);
295  }
296 
300  std::string toString() {
301  return root.toString();
302  }
303 
304 private:
305  OctreeNode<T> root;
306  count dimensions;
307  count numChildrenPerNode;
308 
313  void batchInsert(const std::vector<Vector>& points);
314 
315 
316  std::vector<std::pair<count, Point<T>>> approximateDistance(const OctreeNode<T>& node, const Point<T>& p, const double theta) const;
317  void approximateDistance(const OctreeNode<T>& node, const Point<T>& p, const double theta, std::vector<std::pair<count, Point<T>>>& result) const;
318 
319  template<typename L>
320  void approximateDistance(const OctreeNode<T>& node, const Point<T>& p, const double sqTheta, L& handle) const;
321 };
322 
323 template<typename T>
324 Octree<T>::Octree(const std::vector<Vector>& points) {
325  dimensions = points.size();
326  numChildrenPerNode = pow(2, dimensions);
327  batchInsert(points);
328 // root.compress();
329 }
330 
331 template<typename T>
332 void Octree<T>::recomputeTree(const std::vector<Vector>& points) {
333  root.children.clear();
334  batchInsert(points);
335 }
336 
337 template<typename T>
338 void Octree<T>::batchInsert(const std::vector<Vector>& points) {
339  Point<T> center(dimensions);
340  T sideLength = 0;
341  for (count d = 0; d < dimensions; ++d) {
342  T minVal = points[d][0];
343  T maxVal = points[d][0];
344 
345  for (index i = 1; i < points[d].getDimension(); ++i) {
346  minVal = std::min(minVal, points[d][i]);
347  maxVal = std::max(maxVal, points[d][i]);
348  }
349 
350  sideLength = std::max(sideLength, fabs(maxVal - minVal) * 1.005); // add 0.5% to bounding box
351  center[d] = (minVal + maxVal)/2.0;
352  }
353 
354  root.bBox = {center, sideLength};
355 
356  for (index i = 0; i < points[0].getDimension(); ++i) {
357  Point<T> p(points.size());
358  for (count d = 0; d < dimensions; ++d) {
359  p[d] = points[d][i];
360  }
361 
362  root.addPoint(p, dimensions, numChildrenPerNode);
363  }
364 
365  root.computeCenterOfMass();
366 }
367 
368 template<typename T>
369 std::vector<std::pair<count, Point<T>>> Octree<T>::approximateDistance(const OctreeNode<T>& node, const Point<T>& p, const double theta) const {
370  if (node.isEmpty()) return {};
371  if (node.isLeaf()) {
372  if (node.centerOfMass == p) return {};
373  return {std::make_pair(node.weight, node.centerOfMass)};
374  } else {
375  double dist = p.distance(node.centerOfMass);
376  if (dist == 0 || node.bBox.getSideLength() <= theta * dist) {
377  return {std::make_pair(node.weight, node.centerOfMass)};
378  } else { // split further
379  std::vector<std::pair<count, Point<T>>> points;
380  points.reserve(node.weight);
381  for (auto &child : node.children) {
382  std::vector<std::pair<count, Point<T>>> childPoints = approximateDistance(child, p, theta);
383  points.insert(points.end(), childPoints.begin(), childPoints.end());
384  }
385  return points;
386  }
387  }
388 
389  return {};
390 }
391 
392 template<typename T>
393 void Octree<T>::approximateDistance(const OctreeNode<T>& node, const Point<T>& p, const double theta, std::vector<std::pair<count, Point<T>>>& result) const {
394  if (node.isEmpty()) return;
395  if (node.isLeaf()) {
396  if (node.centerOfMass != p) {
397  result.push_back(std::make_pair(node.weight, node.centerOfMass));
398  }
399  } else {
400  double dist = p.distance(node.centerOfMass);
401  if (dist == 0 || node.bBox.getSideLength() <= theta * dist) {
402  result.push_back(std::make_pair(node.weight, node.centerOfMass));
403  } else { // split further and recurse
404  for (auto &child : node.children) {
405  approximateDistance(child, p, theta, result);
406  }
407  }
408  }
409 }
410 
411 template<typename T> template<typename L>
412 void Octree<T>::approximateDistance(const OctreeNode<T>& node, const Point<T>& p, const double sqTheta, L& handle) const {
413  //if (node.isEmpty()) return;
414  if (!node.isLeaf()) {
415  double sqDist = p.squaredDistance(node.centerOfMass);
416  if (sqDist == 0 || node.bBox.getSqSideLength() <= sqTheta * sqDist) {
417  handle(node.weight, node.centerOfMass, sqDist);
418  } else { // split further and recurse
419  for (auto &child : node.children) {
420  approximateDistance(child, p, sqTheta, handle);
421  }
422  }
423  } else if (node.centerOfMass != p) { // node is leaf and non-empty since octree only stores non-empty nodes
424  handle(node.weight, node.centerOfMass, p.squaredDistance(node.centerOfMass));
425  }
426 }
427 
428 
429 } /* namespace NetworKit */
430 
431 #endif /* NETWORKIT_CPP_VIZ_OCTREE_H_ */
Implementation of a k-dimensional octree for the purpose of Barnes-Hut approximation.
Definition: Octree.h:265
void compress()
Definition: Octree.h:161
void approximateDistance(const Point< T > &p, const double theta, std::vector< std::pair< count, Point< T >>> &result) const
Definition: Octree.h:288
void setCenter(const Point< T > &center)
Sets the center of the bounding box.
Definition: Octree.h:47
Bounding box used by the Octree class.
Definition: Octree.h:24
OctreeNode(BoundingBox< T > &bBox)
Definition: Octree.h:122
BoundingBox()
Constructor for creating an empty bounding box of size 0.
Definition: Octree.h:29
bool isLeaf() const
Definition: Octree.h:127
BoundingBox< T > bBox
Definition: Octree.h:119
uint64_t index
Typedefs.
Definition: Globals.h:20
T getSideLength() const
Definition: Octree.h:72
BoundingBox(const Point< T > &center, const T sideLength)
Constructor for creating a bounding box.
Definition: Octree.h:36
void addPoint(const Point< T > &point, count dimensions, count numChildrenPerNode)
Adds point to octree node.
Definition: Octree.h:225
T getHalfSideLength() const
Definition: Octree.h:79
Point< T > centerOfMass
Definition: Octree.h:117
Formerly marked as deprecated: To take advantage of automatic mapping between C++ and Python data str...
Definition: Point.h:39
OctreeNode()
Definition: Octree.h:121
void setSideLength(T sideLength)
Sets the side length of the bounding box.
Definition: Octree.h:63
count weight
Definition: Octree.h:116
void recomputeTree(const std::vector< Vector > &points)
Clears current content and inserts points in points into the octree.
Definition: Octree.h:332
T distance(const Point< T > &p) const
Definition: Point.h:130
bool contains(const Point< T > &point) const
Definition: Octree.h:93
T getSqSideLength() const
Definition: Octree.h:86
uint64_t count
Definition: Globals.h:21
index node
Definition: Globals.h:23
bool isEmpty() const
Definition: Octree.h:134
std::vector< OctreeNode > children
Definition: Octree.h:118
void split(count dimensions, count numChildren)
Split area corresponding octree node so as to obtain numChildren octree node children.
Definition: Octree.h:205
bool contains(const Point< T > &point) const
Definition: Octree.h:141
BoundingBox(const BoundingBox< T > &other)
Definition: Octree.h:41
Point< T > & getCenter()
Definition: Octree.h:55
void computeCenterOfMass()
Computes octree node's (possibly weighted) center of mass.
Definition: Octree.h:148
std::string toString()
Definition: Octree.h:300
std::string toString()
Definition: Octree.h:189
count getDimensions() const
Definition: Point.h:50
std::vector< std::pair< count, Point< T > > > approximateDistance(const Point< T > &p, const double theta) const
Definition: Octree.h:284
Octree()=default
Default constructor.
void approximateDistance(const Point< T > &p, const double theta, L &handle) const
Definition: Octree.h:293
Node in the octree data structure.
Definition: Octree.h:115