00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef LOOSEOCTREE_H
00019 #define LOOSEOCTREE_H
00020
00021 #include <vector>
00022 #include "octant.h"
00023
00024 using namespace std;
00025
00026 namespace mesh {
00027
00037 template <typename NodeType>
00038 class LooseOctree {
00039 public:
00042 typedef bool (NodeType::*filter)(const Vector3&) const;
00043 public:
00047 LooseOctree(const Vector3& center, double length, int max_depth);
00050 ~LooseOctree();
00055 Octant<NodeType>* add_node(NodeType* node, Octant<NodeType>* start=0, int depth=0);
00064 void find_nodes_by_point(const Vector3& coord,
00065 std::vector<NodeType *>& node_list,
00066 int& nof_visited,
00067 LooseOctree::filter filter,
00068 Octant<NodeType>* octant=0);
00077 NodeType* find_one_node_by_point(const Vector3& coord,
00078 int& nof_visited,
00079 LooseOctree::filter filter,
00080 Octant<NodeType>* octant=0);
00081
00085 void get_statistics(int& max_depth,
00086 double& avg_node_depth,
00087 int& nof_octants,
00088 int& nof_nodes);
00089
00095 void export_vtk();
00098 void export_vtk_level(int level, std::string filename);
00099
00100 protected:
00101 void get_statistics_recursive(Octant<NodeType>* octant,
00102 int depth,
00103 int& max_depth,
00104 double& sum_depth,
00105 int& nof_octants,
00106 int& nof_nodes);
00107
00108 void export_vtk_recursive(Octant<NodeType>* octant,int depth,
00109 int x,
00110 int y,
00111 int z,
00112 int target_depth,
00113 std::vector<int>& cells);
00114
00115 private:
00117 Octant<NodeType>* _root;
00119 int _max_depth;
00120 };
00121
00122
00123 template <typename NodeType>
00124 LooseOctree<NodeType>::LooseOctree(const Vector3& center, double length, int max_depth)
00125 : _max_depth(max_depth)
00126 {
00127 _root = new Octant<NodeType>(center, length);
00128 }
00129
00130 template <typename NodeType>
00131 LooseOctree<NodeType>::~LooseOctree() {
00132 delete _root;
00133 _root = 0;
00134 }
00135
00136 template <typename NodeType>
00137 Octant<NodeType>* LooseOctree<NodeType>::add_node(NodeType* node, Octant<NodeType>* octant, int depth) {
00138
00139 if (octant == 0)
00140 octant = _root;
00141
00142 Vector3 node_extent, node_center;
00143 node->get_bounding_box(node_center, node_extent);
00144
00145
00146
00147 if ((depth < _max_depth ) &&
00148 octant->_length > 2*node_extent.x &&
00149 octant->_length > 2*node_extent.y &&
00150 octant->_length > 2*node_extent.z
00151 ) {
00152
00153 int child_idx = octant->get_child_index_for_point(node_center);
00154
00155
00156 if (octant->_children[child_idx] == 0)
00157 octant->_children[child_idx] = new Octant<NodeType>(octant, child_idx);
00158
00159
00160 return add_node(node, octant->_children[child_idx], depth + 1);
00161
00162 } else {
00163
00164 octant->add_node(node);
00165 return octant;
00166 }
00167 }
00168
00169 template <typename NodeType>
00170 void LooseOctree<NodeType>::find_nodes_by_point(const Vector3& coord,
00171 std::vector<NodeType *>& nodes,
00172 int& nof_visited,
00173 LooseOctree<NodeType>::filter filter,
00174 Octant<NodeType>* octant) {
00175
00176 if (octant == 0)
00177 octant = _root;
00178
00179 Vector3 dist = coord - octant->_center;
00180 double ext_length = 2.0 * octant->_length;
00181 if (fabs(dist.x) > ext_length ||
00182 fabs(dist.y) > ext_length ||
00183 fabs(dist.z) > ext_length)
00184 return;
00185
00186 typename std::vector<NodeType *>::iterator it = octant->_nodes.begin();
00187 while (it != octant->_nodes.end()) {
00188 NodeType* node = *it;
00189 if ((node->*filter)(coord))
00190 nodes.push_back(node);
00191 nof_visited ++;
00192 it ++;
00193 }
00194
00195 for (int i = 0; i < 8; i ++)
00196 if (octant->_children[i])
00197 find_nodes_by_point(coord, nodes, nof_visited, filter, octant->_children[i]);
00198 }
00199
00200 template <typename NodeType>
00201 NodeType* LooseOctree<NodeType>::find_one_node_by_point(const Vector3& coord,
00202 int& nof_visited,
00203 LooseOctree<NodeType>::filter filter,
00204 Octant<NodeType>* octant) {
00205 NodeType* node;
00206
00207
00208 if (octant == 0)
00209 octant = _root;
00210
00211 Vector3 dist = coord - octant->_center;
00212 double ext_length = 2.0 * octant->_length;
00213 if (fabs(dist.x) > ext_length ||
00214 fabs(dist.y) > ext_length ||
00215 fabs(dist.z) > ext_length)
00216 return 0;
00217
00218 typename std::vector<NodeType *>::iterator it = octant->_nodes.begin();
00219 while (it != octant->_nodes.end()) {
00220 node = *it;
00221 nof_visited ++;
00222 it ++;
00223 if ((node->*filter)(coord))
00224 return node;
00225 }
00226
00227 for (int i = 0; i < 8; i ++) {
00228 if (octant->_children[i]) {
00229 node = find_one_node_by_point(coord, nof_visited, filter, octant->_children[i]);
00230
00231 if (node)
00232 return node;
00233 }
00234 }
00235
00236 return 0;
00237 }
00238
00239 template <typename NodeType>
00240 void LooseOctree<NodeType>::get_statistics(int& max_depth,
00241 double& avg_node_depth,
00242 int& nof_octants,
00243 int& nof_nodes) {
00244 Octant<NodeType>* octant = _root;
00245 int depth = 0;
00246 double sum_depth = 0.0;
00247
00248 max_depth = 0;
00249 nof_octants = 0;
00250 nof_nodes = 0;
00251 get_statistics_recursive(octant, depth, max_depth, sum_depth,
00252 nof_octants, nof_nodes);
00253 avg_node_depth = sum_depth / nof_nodes;
00254 }
00255
00256 template <typename NodeType>
00257 void LooseOctree<NodeType>::get_statistics_recursive(Octant<NodeType>* octant,
00258 int depth,
00259 int& max_depth,
00260 double& sum_depth,
00261 int& nof_octants,
00262 int& nof_nodes) {
00263
00264 if (depth > max_depth)
00265 max_depth = depth;
00266 sum_depth += octant->_nodes.size() * depth;
00267 nof_octants ++;
00268 nof_nodes += octant->_nodes.size();
00269
00270
00271 for (int i = 0; i < 8; i ++)
00272 if (octant->_children[i])
00273 get_statistics_recursive(octant->_children[i], depth+1, max_depth, sum_depth,
00274 nof_octants, nof_nodes);
00275 }
00276
00277 template <typename NodeType>
00278 void LooseOctree<NodeType>::export_vtk() {
00279
00280 int max_depth, nof_octants, nof_nodes;
00281 double avg_node_depth;
00282 get_statistics(max_depth, avg_node_depth, nof_octants, nof_nodes);
00283
00284
00285 for (int level = 0; level <= max_depth; ++ level) {
00286 std::ostringstream buf;
00287 buf << "octree_l" << level << ".vtk";
00288 export_vtk_level(level, buf.str());
00289 }
00290 }
00291
00292 template <typename NodeType>
00293 void LooseOctree<NodeType>::export_vtk_level(int level, std::string filename) {
00294 std::vector<int> cells;
00295
00296 export_vtk_recursive(_root, 0, 0, 0, 0, level, cells);
00297 int nof_octants = cells.size() / 9;
00298
00299 std::ofstream of;
00300 of.open(filename.c_str());
00301 rAssert(of.is_open());
00302
00303
00304 of.precision(6);
00305 of << scientific;
00306 of << "# vtk DataFile Version 2.0" << endl;
00307 of << "generated using VtkExport::export_eigenfields" << endl;
00308 of << "ASCII" << endl << endl;
00309 of << "DATASET UNSTRUCTURED_GRID" << endl;
00310
00311
00312 int n = (1 << level) + 1;
00313 Vector3 origin = _root->_center;
00314 origin -= _root->_length;
00315 double side = 2.0*_root->_length / (1 << level);
00316 of << "POINTS " << n*n*n << " float" << endl;
00317
00318 for (int z = 0; z < n; ++ z)
00319 for (int y = 0; y < n; ++ y)
00320 for (int x = 0; x < n; ++ x)
00321 of << origin.x + x*side << " "
00322 << origin.y + y*side << " "
00323 << origin.z + z*side << endl;
00324 of << endl;
00325
00326
00327 of << "CELLS " << nof_octants << " " << cells.size() << endl;
00328 std::vector<int>::iterator it = cells.begin();
00329 for (int i = 0; it != cells.end(); ++ it, ++ i) {
00330 of << *it;
00331 if (i % 9 == 8)
00332 of << endl;
00333 else
00334 of << " ";
00335 }
00336 of << endl;
00337
00338
00339 of << "CELL_TYPES " << nof_octants << endl;
00340 for (int i = 0; i < nof_octants; ++ i)
00341 of << "11" << endl;
00342 of << endl;
00343
00344
00345 of << "CELL_DATA " << nof_octants << endl;
00346 of << "SCALARS cell_scalars int 1" << endl;
00347 of << "LOOKUP_TABLE default" << endl;
00348 for (int i = 0; i < nof_octants; ++ i)
00349 of << "0" << endl;
00350 of << endl;
00351 }
00352
00353 template <typename NodeType>
00354 void LooseOctree<NodeType>::export_vtk_recursive(Octant<NodeType>* octant,
00355 int depth,
00356 int x,
00357 int y,
00358 int z,
00359 int target_depth,
00360 std::vector<int>& cells) {
00361 if (depth == target_depth) {
00362
00363 int n = (1 << depth) + 1;
00364 cells.push_back(8);
00365 cells.push_back((x+0) + (y+0)*n + (z+0)*n*n);
00366 cells.push_back((x+1) + (y+0)*n + (z+0)*n*n);
00367 cells.push_back((x+0) + (y+1)*n + (z+0)*n*n);
00368 cells.push_back((x+1) + (y+1)*n + (z+0)*n*n);
00369 cells.push_back((x+0) + (y+0)*n + (z+1)*n*n);
00370 cells.push_back((x+1) + (y+0)*n + (z+1)*n*n);
00371 cells.push_back((x+0) + (y+1)*n + (z+1)*n*n);
00372 cells.push_back((x+1) + (y+1)*n + (z+1)*n*n);
00373 } else {
00374
00375 for (int i = 0; i < 8; i ++) {
00376 if (octant->_children[i])
00377 export_vtk_recursive(octant->_children[i],
00378 depth+1,
00379 x + ((i & 0x1) == 0x1) * (1 << (target_depth - depth - 1)),
00380 y + ((i & 0x2) == 0x2) * (1 << (target_depth - depth - 1)),
00381 z + ((i & 0x4) == 0x4) * (1 << (target_depth - depth - 1)),
00382 target_depth,
00383 cells);
00384 }
00385 }
00386 }
00387
00388 }
00389
00390 #endif