Flow123d  DF_patch_fe_data_tables-b828b90
bih_tree.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file bih_tree.cc
15  * @brief
16  */
17 
18 #include "mesh/bih_tree.hh"
19 #include "mesh/bih_node.hh"
20 #include "mesh/mesh.h"
21 #include "system/global_defs.h"
22 #include <ctime>
23 #include <stack>
24 
25 /**
26  * Minimum reduction of box size to allow
27  * splitting of a node during tree creation.
28  */
29 const double BIHTree::size_reduce_factor = 0.8;
30 
31 const unsigned int BIHTree::default_leaf_size_limit = 20;
32 
33 
34 BIHTree::BIHTree(unsigned int soft_leaf_size_limit)
35 : leaf_size_limit(soft_leaf_size_limit) //, r_gen(123)
36 {}
37 
38 
40 }
41 
42 
44  if (elements_.size()==0) {
45  // For first call of method set vertices of main_box_ to valid value (default values set in constructor are NaNs)
46  main_box_ = BoundingBox( boxes[0].min() );
47  }
48  for(BoundingBox box : boxes) {
49  this->elements_.push_back(box);
50  main_box_.expand(box);
51  }
52 }
53 
54 
56  ASSERT_GT(elements_.size(), 0);
57 
58  max_n_levels = 2*log2(elements_.size());
59  nodes_.reserve(2*elements_.size() / leaf_size_limit);
60  in_leaves_.resize(elements_.size());
61  for(unsigned int i=0; i<in_leaves_.size(); i++) in_leaves_[i] = i;
62 
63  // make root node
64  nodes_.push_back(BIHNode());
65  nodes_.back().set_leaf(0, in_leaves_.size(), 0, 0);
66  uint height = make_node(main_box_, 0);
67 
68  node_stack_.reserve(2*height);
69 }
70 
71 
72 const BoundingBox& BIHTree::ele_bounding_box(unsigned int ele_idx) const
73 {
74  ASSERT(ele_idx < elements_.size());
75  return elements_[ele_idx];
76 }
77 
78 
79 void BIHTree::split_node(const BoundingBox &node_box, unsigned int node_idx) {
80  BIHNode &node = nodes_[node_idx];
81  ASSERT( node.is_leaf() ).error("Not leaf node.");
82  unsigned int axis = node_box.longest_axis();
83  double median = estimate_median(axis, node);
84 
85  // split elements in node according to the median
86  auto left = in_leaves_.begin() + node.leaf_begin(); // first of unresolved elements in @p in_leaves_
87  auto right = in_leaves_.begin() + node.leaf_end()-1; // last of unresolved elements in @p in_leaves_
88 
89  double left_bound=node_box.min(axis); // max bound of the left group
90  double right_bound=node_box.max(axis); // min bound of the right group
91 
92  while (left != right) {
93  if ( elements_[ *left ].projection_center(axis) < median) {
94  left_bound = std::max( left_bound, elements_[ *left ].max(axis) );
95  ++left;
96  }
97  else {
98  while ( left != right
99  && elements_[ *right ].projection_center(axis) >= median ) {
100  right_bound = std::min( right_bound, elements_[ *right ].min(axis) );
101  --right;
102  }
103  std::swap( *left, *right);
104  }
105  }
106  // in any case left==right is now the first element of the right group
107 
108  if ( elements_[ *left ].projection_center(axis) < median) {
109  left_bound = std::max( left_bound, elements_[ *left ].max(axis) );
110  ++left;
111  ++right;
112  } else {
113  right_bound = std::min( right_bound, elements_[ *right ].min(axis) );
114  }
115 
116  unsigned int left_begin = node.leaf_begin();
117  unsigned int left_end = left - in_leaves_.begin();
118  unsigned int right_end = node.leaf_end();
119  unsigned int depth = node.depth()+1;
120  // create new leaf nodes and possibly call split_node on them
121  // can not use node reference anymore
122  nodes_.push_back(BIHNode());
123  nodes_.back().set_leaf(left_begin, left_end, left_bound, depth);
124  nodes_.push_back(BIHNode());
125  nodes_.back().set_leaf(left_end, right_end, right_bound, depth);
126 
127  nodes_[node_idx].set_non_leaf(nodes_.size()-2, nodes_.size()-1, axis);
128 
129 // DebugOut().fmt("{} {} {} {} {} {} {}\n", node_idx, node_box.min(axis), left_bound, right_bound, node_box.max(axis),
130 // left_end - left_begin, right_end - left_end );
131 }
132 
133 
134 uint BIHTree::make_node(const BoundingBox &box, unsigned int node_idx) {
135  // we must refer to the node by index to prevent seg. fault due to nodes_ reallocation
136 
137  uint height = 0;
138  split_node(box,node_idx);
139 
140  {
141  BIHNode &node = nodes_[node_idx];
142  BIHNode &child = nodes_[ node.child(0) ];
143  BoundingBox node_box(box);
144  node_box.set_max(node.axis(), child.bound() );
145  if ( child.leaf_size() > leaf_size_limit
146  && child.depth() < max_n_levels)
147 // && ( node.axis() != node_box.longest_axis()
148 // || node_box.size(node_box.longest_axis()) < box.size(node.axis()) * size_reduce_factor )
149 // )
150  {
151  uint ht = make_node(node_box, node.child(0) );
152  height = max(height, ht);
153  }
154 // else{
155 // DebugOut().fmt("{} {} {} {}\n", node_idx, child.leaf_size(),
156 // node_box.size(node_box.longest_axis()),
157 // box.size(node.axis()));
158 // }
159  }
160 
161  {
162  BIHNode &node = nodes_[node_idx];
163  BIHNode &child = nodes_[ node.child(1) ];
164  BoundingBox node_box(box);
165  node_box.set_min(node.axis(), child.bound() );
166  if ( child.leaf_size() > leaf_size_limit
167  && child.depth() < max_n_levels)
168 // && ( node.axis() != node_box.longest_axis()
169 // || node_box.size(node_box.longest_axis()) < box.size(node.axis()) * size_reduce_factor )
170 // )
171  {
172  uint ht = make_node(node_box, node.child(1) );
173  height = max(height, ht);
174  }
175 // else{
176 // DebugOut().fmt("{} {} {} {}\n", node_idx, child.leaf_size(),
177 // node_box.size(node_box.longest_axis()),
178 // box.size(node.axis()));
179 // }
180  }
181  return height+1;
182 }
183 
184 
185 double BIHTree::estimate_median(unsigned char axis, const BIHNode &node)
186 {
187  unsigned int median_idx;
188  unsigned int n_elements = node.leaf_size();
189 
190  // TODO: possible optimizations:
191  // - try to apply nth_element directly to in_leaves_ array
192  // - if current approach is better (due to cache memory), check randomization of median for large meshes
193  // - good balancing of tree is crutial both for creation and find method
194 
195 // unsigned int sample_size = 50+n_elements/5;
196 // if (n_elements > sample_size) {
197 // // random sample
198 // std::uniform_int_distribution<unsigned int> distribution(node.leaf_begin(), node.leaf_end()-1);
199 // coors_.resize(sample_size);
200 // for (unsigned int i=0; i<coors_.size(); i++) {
201 // median_idx = distribution(this->r_gen);
202 //
203 // coors_[i] = elements_[ in_leaves_[ median_idx ] ].projection_center(axis);
204 // }
205 //
206 // } else
207  {
208  // all elements
209  coors_.resize(n_elements);
210  for (unsigned int i=0; i<coors_.size(); i++) {
211  median_idx = node.leaf_begin() + i;
212  coors_[i] = elements_[ in_leaves_[ median_idx ] ].projection_center(axis);
213  }
214 
215  }
216 
217  unsigned int median_position = (unsigned int)(coors_.size() / 2);
218  std::nth_element(coors_.begin(), coors_.begin()+median_position, coors_.end());
219 
220  return coors_[median_position];
221 }
222 
223 
224 unsigned int BIHTree::get_element_count() const {
225  return elements_.size();
226 }
227 
228 
230  return main_box_;
231 }
232 
233 
234 void BIHTree::find_bounding_box(const BoundingBox &box, std::vector<unsigned int> &result_list, bool full_list) const
235 {
236 
237  ASSERT_EQ(result_list.size() , 0);
238 
239  unsigned int counter = 0;
240  node_stack_.clear();
241  node_stack_.push_back(0);
242  while (! node_stack_.empty()) {
243  const BIHNode &node = nodes_[node_stack_.back()];
244  //DebugOut().fmt("node: {}\n", node_stack.top() );
245  node_stack_.pop_back();
246 
247 
248  if (node.is_leaf()) {
249 
250  counter ++;
251  //START_TIMER("leaf");
252  for (unsigned int i=node.leaf_begin(); i<node.leaf_end(); i++) {
253  if (full_list || elements_[ in_leaves_[i] ].intersect(box)) {
254 
255  result_list.push_back(in_leaves_[i]);
256  }
257  }
258  //END_TIMER("leaf");
259  } else {
260  //START_TIMER("recursion");
261  if ( ! box.projection_gt( node.axis(), nodes_[node.child(0)].bound() ) ) {
262  // box intersects left group
263  node_stack_.push_back( node.child(0) );
264  }
265  if ( ! box.projection_lt( node.axis(), nodes_[node.child(1)].bound() ) ) {
266  // box intersects right group
267  node_stack_.push_back( node.child(1) );
268  }
269  //END_TIMER("recursion");
270  }
271  }
272  //node_stack_.pop_back();
273  //cout << "stack size: " << node_stack_.size();
274 
275 // DebugOut().fmt("leaves: {}\n", counter);
276 
277 //#ifdef FLOW123D_DEBUG_ASSERTS
278 // // check uniqueness of element indexes
279 // std::vector<unsigned int> cpy(result_list);
280 // sort(cpy.begin(), cpy.end());
281 // std::vector<unsigned int>::iterator it = unique(cpy.begin(), cpy.end());
282 // ASSERT_PERMANENT_EQ(cpy.size() , it - cpy.begin());
283 //#endif
284 }
285 
286 
287 void BIHTree::find_point(const Space<3>::Point &point, std::vector<unsigned int> &result_list, bool full_list) const
288 {
289  find_bounding_box(BoundingBox(point), result_list, full_list);
290 }
291 
292 
293 
#define ASSERT(expr)
Definition: asserts.hh:351
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
#define ASSERT_GT(a, b)
Definition of comparative assert macro (Greater Than) only for debug mode.
Definition: asserts.hh:317
bool is_leaf() const
return true if node is leaf
Definition: bih_node.hh:56
unsigned int axis() const
return axes (coordination of splitting) of inner node
Definition: bih_node.hh:91
unsigned int leaf_begin() const
Definition: bih_node.hh:67
unsigned int leaf_size() const
Definition: bih_node.hh:84
unsigned char depth() const
return depth of leaf node
Definition: bih_node.hh:61
double bound() const
Definition: bih_node.hh:97
unsigned int child(unsigned int i_child) const
Return index of child node.
Definition: bih_node.hh:103
unsigned int leaf_end() const
Definition: bih_node.hh:73
void construct()
Definition: bih_tree.cc:55
double estimate_median(unsigned char axis, const BIHNode &node)
Definition: bih_tree.cc:185
unsigned int max_n_levels
Maximal count of BIH tree levels.
Definition: bih_tree.hh:143
std::vector< double > coors_
temporary vector stored values of coordinations for calculating median
Definition: bih_tree.hh:148
unsigned int get_element_count() const
Definition: bih_tree.cc:224
static const unsigned int default_leaf_size_limit
Default leaf size limit.
Definition: bih_tree.hh:45
std::vector< BoundingBox > elements_
mesh
Definition: bih_tree.hh:132
void add_boxes(const std::vector< BoundingBox > &boxes)
Definition: bih_tree.cc:43
unsigned int leaf_size_limit
Maximal number of elements stored in a leaf node of BIH tree.
Definition: bih_tree.hh:141
void split_node(const BoundingBox &node_box, unsigned int node_idx)
create bounding boxes of element
Definition: bih_tree.cc:79
uint make_node(const BoundingBox &box, unsigned int node_idx)
Definition: bih_tree.cc:134
void find_point(const Space< 3 >::Point &point, std::vector< unsigned int > &result_list, bool full_list=false) const
Definition: bih_tree.cc:287
void find_bounding_box(const BoundingBox &boundingBox, std::vector< unsigned int > &result_list, bool full_list=false) const
Definition: bih_tree.cc:234
BoundingBox main_box_
Main bounding box. (from mesh)
Definition: bih_tree.hh:134
std::vector< BIHNode > nodes_
vector of tree nodes
Definition: bih_tree.hh:139
BIHTree(unsigned int soft_leaf_size_limit=BIHTree::default_leaf_size_limit)
Definition: bih_tree.cc:34
const BoundingBox & ele_bounding_box(unsigned int ele_idx) const
Gets bounding box of element of given index ele_index.
Definition: bih_tree.cc:72
std::vector< unsigned int > in_leaves_
vector stored element indexes in leaf nodes
Definition: bih_tree.hh:146
std::vector< unsigned int > node_stack_
Stack for search algorithms.
Definition: bih_tree.hh:136
const BoundingBox & tree_box() const
Definition: bih_tree.cc:229
static const double size_reduce_factor
required reduction in size of box to allow further splitting
Definition: bih_tree.hh:106
~BIHTree()
Definition: bih_tree.cc:39
Bounding box in 3d ambient space.
Definition: bounding_box.hh:53
const Point & min() const
bool projection_lt(unsigned int axis, double value) const
unsigned char longest_axis() const
bool projection_gt(unsigned int axis, double value) const
void set_min(unsigned int axis, double min)
const Point & max() const
void set_max(unsigned int axis, double max)
void expand(const Point &point)
Armor::ArmaVec< double, spacedim > Point
Definition: point.hh:42
Global macros to enhance readability and debugging, general constants.
unsigned int uint
void swap(nlohmann::json &j1, nlohmann::json &j2) noexcept(is_nothrow_move_constructible< nlohmann::json >::value and is_nothrow_move_assignable< nlohmann::json >::value)
exchanges the values of two JSON objects
Definition: json.hpp:8688