TTK
Loading...
Searching...
No Matches
RangeDrivenOctree.h
Go to the documentation of this file.
1
18
19#pragma once
20
21// base code includes
22#include <Debug.h>
23#include <Triangulation.h>
24
25namespace ttk {
26
27 class RangeDrivenOctree : virtual public Debug {
28 public:
30
31 template <class dataTypeU, class dataTypeV, typename triangulationType>
32 inline int build(const triangulationType *const triangulation);
33
34 inline bool empty() const {
35 return nodeList_.empty();
36 }
37
38 void flush();
39
40 int getTet2NodeMap(std::vector<SimplexId> &map,
41 const bool &forSegmentation = false) const;
42
43 int rangeSegmentQuery(const std::pair<double, double> &p0,
44 const std::pair<double, double> &p1,
45 std::vector<SimplexId> &cellList) const;
46
47 inline void setCellList(const SimplexId *cellList) {
48 cellList_ = cellList;
49 }
50
51 inline void setCellNumber(const SimplexId &cellNumber) {
52 cellNumber_ = cellNumber;
53 }
54
55 inline void setLeafMinimumCellNumber(const SimplexId &number) {
57 }
58
59 inline void setLeafMinimumDomainVolumeRatio(const float &ratio) {
61 }
62
63 inline void setLeafMinimumRangeAreaRatio(const float &ratio) {
65 }
66
67 inline void setPointList(const float *pointList) {
68 pointList_ = pointList;
69 }
70
71 inline void setRange(const void *u, const void *v) {
72 u_ = u;
73 v_ = v;
74 }
75
76 inline void setVertexNumber(const SimplexId &vertexNumber) {
77 vertexNumber_ = vertexNumber;
78 }
79
80 int stats(std::ostream &stream);
81
82 int statNode(const SimplexId &nodeId, std::ostream &stream);
83
84 protected:
85 struct OctreeNode {
86 std::pair<std::pair<double, double>, std::pair<double, double>>
88 std::vector<SimplexId> cellList_{};
89 std::vector<SimplexId> childList_{};
90 std::array<std::pair<float, float>, 3> domainBox_{};
91 };
92
93 template <class dataTypeU, class dataTypeV>
94 int buildNode(const std::vector<SimplexId> &cellList,
95 const std::array<std::pair<float, float>, 3> &domainBox,
96 const std::pair<std::pair<double, double>,
97 std::pair<double, double>> &rangeBox,
98 SimplexId &nodeId);
99
100 int rangeSegmentQuery(const std::pair<double, double> &p0,
101 const std::pair<double, double> &p1,
102 const SimplexId &nodeId,
103 std::vector<SimplexId> &cellList) const;
104
105 bool segmentIntersection(const std::pair<double, double> &p0,
106 const std::pair<double, double> &p1,
107 const std::pair<double, double> &q0,
108 const std::pair<double, double> &q1) const;
109
110 const void *u_{};
111 const void *v_{};
112 const float *pointList_{};
119 std::vector<OctreeNode> nodeList_{};
120 std::vector<std::array<std::pair<float, float>, 3>> cellDomainBox_{};
121 std::vector<std::pair<std::pair<double, double>, std::pair<double, double>>>
123 };
124} // namespace ttk
125
126// if the package is not a template, comment the following line
127// #include <RangeDrivenOctree.cpp>
128
129template <class dataTypeU, class dataTypeV, typename triangulationType>
131 const triangulationType *const triangulation) {
132
133 Timer t;
134
135 const dataTypeU *u = (const dataTypeU *)u_;
136 const dataTypeV *v = (const dataTypeV *)v_;
137
138 if(triangulation) {
139 cellNumber_ = triangulation->getNumberOfCells();
140 }
141
142 if(triangulation) {
143 vertexNumber_ = triangulation->getNumberOfVertices();
144 }
145
147
149
150 // WARNING: assuming tets only here
151#ifdef TTK_ENABLE_OPENMP
152#pragma omp parallel for num_threads(threadNumber_)
153#endif
154 for(SimplexId i = 0; i < cellNumber_; i++) {
155
156 for(int j = 0; j < 3; j++) {
157 cellDomainBox_[i][j].first = FLT_MAX;
158 cellDomainBox_[i][j].second = -FLT_MAX;
159 }
160
161 const SimplexId *cell = nullptr;
162 if(!triangulation) {
163 cell = &(cellList_[5 * i + 1]);
164 }
165
166 for(int j = 0; j < 4; j++) {
167
168 // update the domain bounding box for that tet
169 SimplexId vertexId = 0;
170
171 if(triangulation) {
172 triangulation->getCellVertex(i, j, vertexId);
173 } else if(cell != nullptr) {
174 vertexId = cell[j];
175 }
176
177 std::array<float, 3> p{};
178 if(triangulation) {
179 triangulation->getVertexPoint(vertexId, p[0], p[1], p[2]);
180 } else {
181 p[0] = pointList_[3 * vertexId];
182 p[1] = pointList_[3 * vertexId + 1];
183 p[2] = pointList_[3 * vertexId + 2];
184 }
185
186 // update x-range
187 if(p[0] < cellDomainBox_[i][0].first)
188 cellDomainBox_[i][0].first = p[0];
189 if(p[0] > cellDomainBox_[i][0].second)
190 cellDomainBox_[i][0].second = p[0];
191
192 // update y-range
193 if(p[1] < cellDomainBox_[i][1].first)
194 cellDomainBox_[i][1].first = p[1];
195 if(p[1] > cellDomainBox_[i][1].second)
196 cellDomainBox_[i][1].second = p[1];
197
198 // update z-range
199 if(p[2] < cellDomainBox_[i][2].first)
200 cellDomainBox_[i][2].first = p[2];
201 if(p[2] > cellDomainBox_[i][2].second)
202 cellDomainBox_[i][2].second = p[2];
203
204 // update the range bounding box
205 if(!j) {
206 cellRangeBox_[i].first.first = u[vertexId];
207 cellRangeBox_[i].first.second = u[vertexId];
208
209 cellRangeBox_[i].second.first = v[vertexId];
210 cellRangeBox_[i].second.second = v[vertexId];
211 } else {
212
213 // update u
214 if(u[vertexId] < cellRangeBox_[i].first.first)
215 cellRangeBox_[i].first.first = u[vertexId];
216 if(u[vertexId] > cellRangeBox_[i].first.second)
217 cellRangeBox_[i].first.second = u[vertexId];
218
219 // update v
220 if(v[vertexId] < cellRangeBox_[i].second.first)
221 cellRangeBox_[i].second.first = v[vertexId];
222 if(v[vertexId] > cellRangeBox_[i].second.second)
223 cellRangeBox_[i].second.second = v[vertexId];
224 }
225 }
226 }
227
228 std::vector<SimplexId> domain(cellNumber_);
229 for(SimplexId i = 0; i < cellNumber_; i++)
230 domain[i] = i;
231
232 // get global bBoxes
233 std::array<std::pair<float, float>, 3> domainBox{};
234 std::pair<std::pair<double, double>, std::pair<double, double>> rangeBox;
235
236 for(SimplexId i = 0; i < vertexNumber_; i++) {
237
238 // domain one
239 std::array<float, 3> p{};
240 if(triangulation) {
241 triangulation->getVertexPoint(i, p[0], p[1], p[2]);
242 } else {
243 p[0] = pointList_[3 * i];
244 p[1] = pointList_[3 * i + 1];
245 p[2] = pointList_[3 * i + 2];
246 }
247
248 for(int j = 0; j < 3; j++) {
249 if(!i) {
250 domainBox[j].first = domainBox[j].second = p[j];
251 } else {
252 if(p[j] < domainBox[j].first)
253 domainBox[j].first = p[j];
254 if(p[j] > domainBox[j].second)
255 domainBox[j].second = p[j];
256 }
257 }
258
259 if(!i) {
260 rangeBox.first.first = rangeBox.first.second = u[i];
261 rangeBox.second.first = rangeBox.second.second = v[i];
262 } else {
263 if(u[i] < rangeBox.first.first)
264 rangeBox.first.first = u[i];
265 if(u[i] > rangeBox.first.second)
266 rangeBox.first.second = u[i];
267
268 if(v[i] < rangeBox.second.first)
269 rangeBox.second.first = v[i];
270 if(v[i] > rangeBox.second.second)
271 rangeBox.second.second = v[i];
272 }
273 }
274
275 rangeArea_ = (rangeBox.first.second - rangeBox.first.first)
276 * (rangeBox.second.second - rangeBox.second.first);
277 domainVolume_ = (domainBox[0].second - domainBox[0].first)
278 * (domainBox[1].second - domainBox[1].first)
279 * (domainBox[2].second - domainBox[2].first);
280
281 // special case for tets obtained from regular grid subdivision (assuming 6)
284
285 leafMinimumDomainVolumeRatio_ = (1.0 / ((float)cellNumber_)) / 2.0;
286
287 this->printMsg(
288 "Range area ratio: " + std::to_string(leafMinimumRangeAreaRatio_),
290
291 buildNode<dataTypeU, dataTypeV>(domain, domainBox, rangeBox, rootId_);
292
293 this->printMsg("Octree built", 1.0, t.getElapsedTime(), this->threadNumber_);
294
295 // debug
296 // stats(std::cout);
297 // end of debug
298
299 return 0;
300}
301
302template <class dataTypeU, class dataTypeV>
304 const std::vector<SimplexId> &cellList,
305 const std::array<std::pair<float, float>, 3> &domainBox,
306 const std::pair<std::pair<double, double>, std::pair<double, double>>
307 &rangeBox,
308 SimplexId &nodeId) {
309
310 nodeId = nodeList_.size();
311
312 nodeList_.emplace_back();
313
314 nodeList_.back().rangeBox_ = rangeBox;
315 nodeList_.back().domainBox_ = domainBox;
316
317 float const rangeArea = (rangeBox.first.second - rangeBox.first.first)
318 * (rangeBox.second.second - rangeBox.second.first);
319
320 float const domainVolume = (domainBox[0].second - domainBox[0].first)
321 * (domainBox[1].second - domainBox[1].first)
322 * (domainBox[2].second - domainBox[2].first);
323
324 if(((SimplexId)cellList.size() > leafMinimumCellNumber_)
325 && (rangeArea > leafMinimumRangeAreaRatio_ * rangeArea_)
326 && (domainVolume > leafMinimumDomainVolumeRatio_ * domainVolume_)) {
327
328 nodeList_.back().childList_.resize(8);
329
330 std::array<std::array<std::pair<float, float>, 3>, 8> childDomainBox{};
331 std::array<std::pair<std::pair<dataTypeU, dataTypeU>,
332 std::pair<dataTypeV, dataTypeV>>,
333 8>
334 childRangeBox{};
335 std::array<std::vector<SimplexId>, 8> childCellList{};
336
337 float const midX
338 = domainBox[0].first + (domainBox[0].second - domainBox[0].first) / 2.0;
339 float const midY
340 = domainBox[1].first + (domainBox[1].second - domainBox[1].first) / 2.0;
341 float const midZ
342 = domainBox[2].first + (domainBox[2].second - domainBox[2].first) / 2.0;
343
344 // 0 - - -
345 childDomainBox[0][0].first = domainBox[0].first;
346 childDomainBox[0][0].second = midX;
347 childDomainBox[0][1].first = domainBox[1].first;
348 childDomainBox[0][1].second = midY;
349 childDomainBox[0][2].first = domainBox[2].first;
350 childDomainBox[0][2].second = midZ;
351
352 // 1 - - +
353 childDomainBox[1][0].first = domainBox[0].first;
354 childDomainBox[1][0].second = midX;
355 childDomainBox[1][1].first = domainBox[1].first;
356 childDomainBox[1][1].second = midY;
357 childDomainBox[1][2].first = midZ;
358 childDomainBox[1][2].second = domainBox[2].second;
359
360 // 2 - + -
361 childDomainBox[2][0].first = domainBox[0].first;
362 childDomainBox[2][0].second = midX;
363 childDomainBox[2][1].first = midY;
364 childDomainBox[2][1].second = domainBox[1].second;
365 childDomainBox[2][2].first = domainBox[2].first;
366 childDomainBox[2][2].second = midZ;
367
368 // 3 - + +
369 childDomainBox[3][0].first = domainBox[0].first;
370 childDomainBox[3][0].second = midX;
371 childDomainBox[3][1].first = midY;
372 childDomainBox[3][1].second = domainBox[1].second;
373 childDomainBox[3][2].first = midZ;
374 childDomainBox[3][2].second = domainBox[2].second;
375
376 // 4 + - -
377 childDomainBox[4][0].first = midX;
378 childDomainBox[4][0].second = domainBox[0].second;
379 childDomainBox[4][1].first = domainBox[1].first;
380 childDomainBox[4][1].second = midY;
381 childDomainBox[4][2].first = domainBox[2].first;
382 childDomainBox[4][2].second = midZ;
383
384 // 5 + - +
385 childDomainBox[5][0].first = midX;
386 childDomainBox[5][0].second = domainBox[0].second;
387 childDomainBox[5][1].first = domainBox[1].first;
388 childDomainBox[5][1].second = midY;
389 childDomainBox[5][2].first = midZ;
390 childDomainBox[5][2].second = domainBox[2].second;
391
392 // 6 + + -
393 childDomainBox[6][0].first = midX;
394 childDomainBox[6][0].second = domainBox[0].second;
395 childDomainBox[6][1].first = midY;
396 childDomainBox[6][1].second = domainBox[1].second;
397 childDomainBox[6][2].first = domainBox[2].first;
398 childDomainBox[6][2].second = midZ;
399
400 // 7 + + +
401 childDomainBox[7][0].first = midX;
402 childDomainBox[7][0].second = domainBox[0].second;
403 childDomainBox[7][1].first = midY;
404 childDomainBox[7][1].second = domainBox[1].second;
405 childDomainBox[7][2].first = midZ;
406 childDomainBox[7][2].second = domainBox[2].second;
407
408 for(SimplexId i = 0; i < (SimplexId)cellList.size(); i++) {
409
410 SimplexId childId = 0;
411
412 for(int j = 0; j < 8; j++) {
413 if((cellDomainBox_[cellList[i]][0].first >= childDomainBox[j][0].first)
414 && (cellDomainBox_[cellList[i]][0].first
415 < childDomainBox[j][0].second)
416 && (cellDomainBox_[cellList[i]][1].first
417 >= childDomainBox[j][1].first)
418 && (cellDomainBox_[cellList[i]][1].first
419 < childDomainBox[j][1].second)
420 && (cellDomainBox_[cellList[i]][2].first
421 >= childDomainBox[j][2].first)
422 && (cellDomainBox_[cellList[i]][2].first
423 < childDomainBox[j][2].second)) {
424
425 childId = j;
426 break;
427 }
428 }
429
430 // update child's range box
431 if(childCellList[childId].empty()) {
432 childRangeBox[childId].first.first
433 = cellRangeBox_[cellList[i]].first.first;
434 childRangeBox[childId].first.second
435 = cellRangeBox_[cellList[i]].first.second;
436
437 childRangeBox[childId].second.first
438 = cellRangeBox_[cellList[i]].second.first;
439 childRangeBox[childId].second.second
440 = cellRangeBox_[cellList[i]].second.second;
441 } else {
442 if(cellRangeBox_[cellList[i]].first.first
443 < childRangeBox[childId].first.first) {
444 childRangeBox[childId].first.first
445 = cellRangeBox_[cellList[i]].first.first;
446 }
447 if(cellRangeBox_[cellList[i]].first.second
448 > childRangeBox[childId].first.second) {
449 childRangeBox[childId].first.second
450 = cellRangeBox_[cellList[i]].first.second;
451 }
452
453 if(cellRangeBox_[cellList[i]].second.first
454 < childRangeBox[childId].second.first) {
455 childRangeBox[childId].second.first
456 = cellRangeBox_[cellList[i]].second.first;
457 }
458 if(cellRangeBox_[cellList[i]].second.second
459 > childRangeBox[childId].second.second) {
460 childRangeBox[childId].second.second
461 = cellRangeBox_[cellList[i]].second.second;
462 }
463 }
464
465 childCellList[childId].push_back(cellList[i]);
466 }
467
468 for(int i = 0; i < 8; i++) {
469 buildNode<dataTypeU, dataTypeV>(childCellList[i], childDomainBox[i],
470 childRangeBox[i],
471 nodeList_[nodeId].childList_[i]);
472 }
473 } else {
474 // leaf
475 nodeList_[nodeId].cellList_ = cellList;
476 }
477
478 return 0;
479}
Minimalist debugging class.
Definition Debug.h:88
TTK optional package for octree based range queries in bivariate volumetric data.
int rangeSegmentQuery(const std::pair< double, double > &p0, const std::pair< double, double > &p1, std::vector< SimplexId > &cellList) const
void setLeafMinimumRangeAreaRatio(const float &ratio)
const SimplexId * cellList_
std::vector< std::array< std::pair< float, float >, 3 > > cellDomainBox_
void setCellList(const SimplexId *cellList)
int build(const triangulationType *const triangulation)
void setLeafMinimumCellNumber(const SimplexId &number)
bool segmentIntersection(const std::pair< double, double > &p0, const std::pair< double, double > &p1, const std::pair< double, double > &q0, const std::pair< double, double > &q1) const
std::vector< OctreeNode > nodeList_
void setVertexNumber(const SimplexId &vertexNumber)
void setRange(const void *u, const void *v)
void setPointList(const float *pointList)
int stats(std::ostream &stream)
void setCellNumber(const SimplexId &cellNumber)
int getTet2NodeMap(std::vector< SimplexId > &map, const bool &forSegmentation=false) const
int buildNode(const std::vector< SimplexId > &cellList, const std::array< std::pair< float, float >, 3 > &domainBox, const std::pair< std::pair< double, double >, std::pair< double, double > > &rangeBox, SimplexId &nodeId)
int statNode(const SimplexId &nodeId, std::ostream &stream)
void setLeafMinimumDomainVolumeRatio(const float &ratio)
std::vector< std::pair< std::pair< double, double >, std::pair< double, double > > > cellRangeBox_
double getElapsedTime()
Definition Timer.h:15
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:22
std::array< std::pair< float, float >, 3 > domainBox_
std::vector< SimplexId > cellList_
std::pair< std::pair< double, double >, std::pair< double, double > > rangeBox_
std::vector< SimplexId > childList_
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)