MVE - Multi-View Environment mve-devel
Loading...
Searching...
No Matches
octree.cc
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015, Simon Fuhrmann
3 * TU Darmstadt - Graphics, Capture and Massively Parallel Computing
4 * All rights reserved.
5 *
6 * This software may be modified and distributed under the terms
7 * of the BSD 3-Clause license. See the LICENSE.txt file for details.
8 */
9
10#include <stdexcept>
11#include <list>
12#include <iostream>
13#include <algorithm>
14
15#include "util/timer.h"
16#include "mve/mesh_io.h"
17#include "fssr/octree.h"
18
20
21Octree::Node*
22Octree::Iterator::first_node (void)
23{
24 this->current = this->root;
25 this->level = 0;
26 this->path = 0;
27 return this->current;
28}
29
31Octree::Iterator::first_leaf (void)
32{
33 this->first_node();
34 while (this->current->children != nullptr)
35 {
36 this->current = this->current->children;
37 this->level = this->level + 1;
38 this->path = this->path << 3;
39 }
40 return this->current;
41}
42
44Octree::Iterator::next_node (void)
45{
46 if (this->current->children == nullptr)
47 return this->next_branch();
48
49 this->current = this->current->children;
50 this->level = this->level + 1;
51 this->path = this->path << 3;
52 return this->current;
53}
54
56Octree::Iterator::next_branch (void)
57{
58 if (this->current->parent == nullptr)
59 {
60 this->current = nullptr;
61 return nullptr;
62 }
63
64 if (this->current - this->current->parent->children == 7)
65 {
66 this->current = this->current->parent;
67 this->level = this->level - 1;
68 this->path = this->path >> 3;
69 return this->next_branch();
70 }
71
72 this->current += 1;
73 this->path += 1;
74 return this->current;
75}
76
78Octree::Iterator::next_leaf (void)
79{
80 if (this->current->children != nullptr)
81 {
82 while (this->current->children != nullptr)
83 {
84 this->current = this->current->children;
85 this->level = this->level + 1;
86 this->path = this->path << 3;
87 }
88 return this->current;
89 }
90
91 this->next_branch();
92 if (this->current == nullptr)
93 return nullptr;
94 while (this->current->children != nullptr)
95 {
96 this->current = this->current->children;
97 this->level = this->level + 1;
98 this->path = this->path << 3;
99 }
100 return this->current;
101}
102
104Octree::Iterator::descend (int octant) const
105{
106 Iterator iter(*this);
107 iter.current = iter.current->children + octant;
108 iter.level = iter.level + 1;
109 iter.path = (iter.path << 3) | octant;
110 return iter;
111}
112
114Octree::Iterator::descend (uint8_t level, uint64_t path) const
115{
116 Iterator iter;
117 iter.root = this->root;
118 iter.current = this->root;
119 iter.path = 0;
120 iter.level = 0;
121 for (int i = 0; i < level; ++i)
122 {
123 if (iter.current->children == nullptr)
124 {
125 iter.current = nullptr;
126 return iter;
127 }
128
129 int const octant = (path >> ((level - i - 1) * 3)) & 7;
130 iter = iter.descend(octant);
131 }
132
133 if (iter.path != path || iter.level != level)
134 throw std::runtime_error("descend(): failed");
135
136 return iter;
137}
138
140Octree::Iterator::ascend (void) const
141{
142 Iterator iter;
143 iter.root = this->root;
144 iter.current = this->current->parent;
145 iter.path = this->path >> 3;
146 iter.level = this->level - 1;
147 return iter;
148}
149
150/* -------------------------------------------------------------------- */
151
152void
153Octree::insert_samples (SampleList const& samples)
154{
155 for (std::size_t i = 0; i < samples.size(); i++)
156 this->insert_sample(samples[i]);
157}
158
159void
160Octree::insert_sample (Sample const& sample)
161{
162 if (this->root == nullptr)
163 {
164 this->root = new Node();
165 this->root_center = sample.pos;
166 this->root_size = sample.scale;
167 this->num_nodes = 1;
168 }
169
170 /* Expand octree root if sample is outside the octree. */
171 while (!this->is_inside_octree(sample.pos))
172 this->expand_root_for_point(sample.pos);
173
174 /* Find node by expanding the root or descending the tree. */
175 Node* node = nullptr;
176 if (sample.scale >= this->root_size * 2.0)
177 node = this->find_node_expand(sample);
178 else
179 node = this->find_node_descend(sample, this->get_iterator_for_root());
180
181 node->samples.push_back(sample);
182 this->num_samples += 1;
183}
184
185void
186Octree::create_children (Node* node)
187{
188 if (node->children != nullptr)
189 throw std::runtime_error("create_children(): Children exist!");
190 node->children = new Node[8];
191 this->num_nodes += 8;
192 for (int i = 0; i < 8; ++i)
193 node->children[i].parent = node;
194}
195
196bool
197Octree::is_inside_octree (math::Vec3d const& pos)
198{
199 double const len2 = this->root_size / 2.0;
200 for (int i = 0; i < 3; ++i)
201 if (pos[i] < this->root_center[i] - len2
202 || pos[i] > this->root_center[i] + len2)
203 return false;
204 return true;
205}
206
207void
208Octree::expand_root_for_point (math::Vec3d const& pos)
209{
210 /* Compute old root octant and new root center and size. */
211 int octant = 0;
212 for (int i = 0; i < 3; ++i)
213 if (pos[i] > this->root_center[i])
214 {
215 this->root_center[i] += this->root_size / 2.0;
216 }
217 else
218 {
219 octant |= (1 << i);
220 this->root_center[i] -= this->root_size / 2.0;
221 }
222 this->root_size *= 2.0;
223
224 /* Create new root. */
225 Node* new_root = new Node();
226 this->create_children(new_root);
227 std::swap(new_root->children[octant].children, this->root->children);
228 std::swap(new_root->children[octant].samples, this->root->samples);
229 delete this->root;
230 this->root = new_root;
231
232 /* Fix parent pointers of old child nodes. */
233 if (this->root->children[octant].children != nullptr)
234 {
235 Node* children = this->root->children[octant].children;
236 for (int i = 0; i < 8; ++i)
237 children[i].parent = this->root->children + octant;
238 }
239}
240
241Octree::Node*
242Octree::find_node_descend (Sample const& sample, Iterator const& iter)
243{
244 math::Vec3d node_center;
245 double node_size;
246 this->node_center_and_size(iter, &node_center, &node_size);
247
248 if (sample.scale > node_size * 2.0)
249 throw std::runtime_error("find_node_descend(): Sanity check failed!");
250
251 /*
252 * The current level l is appropriate if sample scale s is
253 * scale(l) <= s < scale(l) * 2. As a sanity check, this function
254 * must not be called if s >= scale(l) * 2. If the current level is
255 * the maximum allowed level, return this node also. Descend otherwise.
256 */
257 if (node_size <= sample.scale || iter.level >= this->max_level)
258 return iter.current;
259
260 /* Descend octree. Find octant and create children if required. */
261 int octant = 0;
262 for (int i = 0; i < 3; ++i)
263 if (sample.pos[i] > node_center[i])
264 octant |= (1 << i);
265 if (iter.current->children == nullptr)
266 this->create_children(iter.current);
267 return this->find_node_descend(sample, iter.descend(octant));
268}
269
270Octree::Node*
271Octree::find_node_expand (Sample const& sample)
272{
273 if (this->root_size > sample.scale)
274 throw std::runtime_error("find_node_expand(): Sanity check failed!");
275
276 /*
277 * The current level l is appropriate if sample scale s is
278 * scale(l) <= scale < scale(l) * 2. As a sanity check, this function
279 * must not be called if scale(l) > s. Otherwise expand.
280 */
281 if (sample.scale < this->root_size * 2.0)
282 return this->root;
283
284 this->expand_root_for_point(sample.pos);
285 return this->find_node_expand(sample);
286}
287
288int
289Octree::get_num_levels (Node const* node) const
290{
291 if (node == nullptr)
292 return 0;
293 if (node->children == nullptr)
294 return 1;
295 int depth = 0;
296 for (int i = 0; i < 8; ++i)
297 depth = std::max(depth, this->get_num_levels(node->children + i));
298 return depth + 1;
299}
300
301void
302Octree::get_samples_per_level (std::vector<std::size_t>* stats,
303 Node const* node, std::size_t level) const
304{
305 if (node == nullptr)
306 return;
307 if (stats->size() <= level)
308 stats->resize(level + 1, 0);
309 stats->at(level) += node->samples.size();
310
311 /* Descend into octree. */
312 if (node->children == nullptr)
313 return;
314 for (int i = 0; i < 8; ++i)
315 this->get_samples_per_level(stats, node->children + i, level + 1);
316}
317
318void
319Octree::node_center_and_size (Iterator const& iter,
320 math::Vec3d* center, double *size) const
321{
322 *center = this->root_center;
323 *size = this->root_size;
324 for (int i = 0; i < iter.level; ++i)
325 {
326 int const octant = iter.path >> ((iter.level - i - 1) * 3);
327 double const offset = *size / 4.0;
328 for (int j = 0; j < 3; ++j)
329 (*center)[j] += ((octant & (1 << j)) ? offset : -offset);
330 *size /= 2.0;
331 }
332}
333
335Octree::get_iterator_for_root (void) const
336{
337 if (this->root == nullptr)
338 throw std::logic_error("Iterator request on empty octree");
339
340 Iterator iter;
341 iter.root = this->root;
342 iter.first_node();
343 return iter;
344}
345
346void
347Octree::influence_query (math::Vec3d const& pos, double factor,
348 std::vector<Sample const*>* result, Iterator const& iter,
349 math::Vec3d const& parent_node_center) const
350{
351 if (iter.current == nullptr)
352 return;
353
354 /*
355 * Strategy is the following: Try to rule out this octree node. Assume
356 * the largest scale sample (node_size * 2) in this node and compute
357 * an estimate for the closest possible distance of any sample in the
358 * node to the query. If 'factor' times the largest scale is less than
359 * the closest distance, the node can be skipped and traversal stops.
360 * Otherwise all samples in the node have to be tested.
361 *
362 * - Note: The 'factor' depends on the basis/weighting function. In this
363 * implementation, factor is always 3.0.
364 * - Note: Nodes can contain samples with scale values much smaller than
365 * node_size. This is because the octree depth is limited.
366 */
367 /* Compute current node center based on parent's. */
368 uint32_t x = (iter.path & 1) >> 0;
369 uint32_t y = (iter.path & 2) >> 1;
370 uint32_t z = (iter.path & 4) >> 2;
371 double node_size = this->root_size / (1 << iter.level);
372 double offset = (iter.level > 0) * node_size / 2.0;
373 math::Vec3d node_center(
374 parent_node_center[0] - offset + x * node_size,
375 parent_node_center[1] - offset + y * node_size,
376 parent_node_center[2] - offset + z * node_size);
377
378 /* Estimate for the minimum distance. No sample is closer to pos. */
379 double const min_distance = (pos - node_center).norm()
380 - MATH_SQRT3 * node_size / 2.0;
381 double const max_scale = node_size * 2.0;
382 if (min_distance > max_scale * factor)
383 return;
384
385 /* Node could not be ruled out. Test all samples. */
386 for (std::size_t i = 0; i < iter.current->samples.size(); ++i)
387 {
388 Sample const& s = iter.current->samples[i];
389 if ((pos - s.pos).square_norm() > MATH_POW2(factor * s.scale))
390 continue;
391 result->push_back(&s);
392 }
393
394 /* Descend into octree. */
395 if (iter.current->children == nullptr)
396 return;
397 for (int i = 0; i < 8; ++i)
398 this->influence_query(pos, factor, result, iter.descend(i),
399 node_center);
400}
401
402void
403Octree::refine_octree (void)
404{
405 if (this->root == nullptr)
406 return;
407
408 std::list<Node*> queue;
409 queue.push_back(this->root);
410 while (!queue.empty())
411 {
412 Node* node = queue.front();
413 queue.pop_front();
414
415 if (node->children == nullptr)
416 this->create_children(node);
417 else
418 for (int i = 0; i < 8; ++i)
419 queue.push_back(node->children + i);
420 }
421}
422
423void
424Octree::limit_octree_level (void)
425{
426 std::cout << "Limiting octree to "
427 << this->max_level << " levels..." << std::endl;
428
429 if (this->root == nullptr)
430 return;
431 this->limit_octree_level(this->root, nullptr, 0);
432}
433
434void
435Octree::limit_octree_level (Node* node, Node* parent, int level)
436{
437 if (level == this->max_level)
438 parent = node;
439
440 if (level > this->max_level)
441 {
442 parent->samples.insert(parent->samples.end(),
443 node->samples.begin(), node->samples.end());
444 node->samples.clear();
445 }
446
447 if (node->children != nullptr)
448 for (int i = 0; i < 8; ++i)
449 this->limit_octree_level(node->children + i, parent, level + 1);
450
451 if (level > this->max_level)
452 this->num_nodes -= 1;
453
454 if (level == this->max_level)
455 {
456 delete [] node->children;
457 node->children = nullptr;
458 }
459}
460
461void
462Octree::print_stats (std::ostream& out)
463{
464 out << "Octree contains " << this->get_num_samples()
465 << " samples in " << this->get_num_nodes() << " nodes on "
466 << this->get_num_levels() << " levels." << std::endl;
467
468 std::vector<std::size_t> octree_stats;
469 this->get_samples_per_level(&octree_stats);
470
471 bool printed = false;
472 for (std::size_t i = 0; i < octree_stats.size(); ++i)
473 {
474 if (!printed && octree_stats[i] == 0)
475 continue;
476 else
477 {
478 out << " Level " << i << ": "
479 << octree_stats[i] << " samples" << std::endl;
480 printed = true;
481 }
482 }
483}
484
#define FSSR_NAMESPACE_END
Definition defines.h:14
#define FSSR_NAMESPACE_BEGIN
Definition defines.h:13
#define MATH_SQRT3
Definition defines.h:50
#define MATH_POW2(x)
Definition defines.h:68
std::vector< Sample > SampleList
Representation of a list of samples.
Definition sample.h:31
void swap(mve::Image< T > &a, mve::Image< T > &b)
Specialization of std::swap for efficient image swapping.
Definition image.h:478
Octree iterator that keeps track of level and path through the octree.
Definition octree.h:58
Iterator descend(int octant) const
Definition octree.cc:104
Node * first_node(void)
Definition octree.cc:22
Simple recursive octree node that stores samples in a vector.
Definition octree.h:39
std::vector< Sample > samples
Definition octree.h:48
Node * children
Definition octree.h:45
Representation of a point sample.
Definition sample.h:22
math::Vec3f pos
Definition sample.h:23
float scale
Definition sample.h:26