TTK
Loading...
Searching...
No Matches
RangeDrivenOctree.cpp
Go to the documentation of this file.
1// NOTE
2// this class has been implemented recursively. this is bad programming.
3// a non-recursive version will be implemented asap.
4
5#include <RangeDrivenOctree.h>
6
7using namespace ttk;
8
10 this->setDebugMsgPrefix("RangeDrivenOctree");
11}
12
14 nodeList_.clear();
15 cellDomainBox_.clear();
16 cellRangeBox_.clear();
17}
18
19int RangeDrivenOctree::getTet2NodeMap(std::vector<SimplexId> &map,
20 const bool &forSegmentation) const {
21
22 std::vector<SimplexId> randomMap;
23 if(forSegmentation) {
24 randomMap.resize(nodeList_.size());
25 for(size_t i = 0; i < randomMap.size(); i++) {
26 randomMap[i] = rand() % (randomMap.size());
27 }
28 }
29
30 map.resize(cellNumber_);
31 for(size_t i = 0; i < nodeList_.size(); i++) {
32 for(size_t j = 0; j < nodeList_[i].cellList_.size(); j++) {
33 if(forSegmentation) {
34 map[nodeList_[i].cellList_[j]] = randomMap[i];
35 } else {
36 map[nodeList_[i].cellList_[j]] = i;
37 }
38 }
39 }
40
41 return 0;
42}
43
45 const std::pair<double, double> &p0,
46 const std::pair<double, double> &p1,
47 std::vector<SimplexId> &cellList) const {
48
49 Timer t;
50
52 cellList.clear();
53
54 SimplexId const ret = rangeSegmentQuery(p0, p1, rootId_, cellList);
55
56 if(debugLevel_ >= static_cast<int>(debug::Priority::DETAIL)) {
57 this->printMsg("Query done", 1.0, t.getElapsedTime(), this->threadNumber_,
59 this->printMsg(
60 std::vector<std::vector<std::string>>{
61 {"#Non empty leaves", std::to_string(cellList.size())}},
63 }
64
65 return ret;
66}
67
69 const std::pair<double, double> &p0,
70 const std::pair<double, double> &p1,
71 const SimplexId &nodeId,
72 std::vector<SimplexId> &cellList) const {
73
74 // check for intersection for each segment of the range bounding box
75 std::pair<double, double> q0, q1;
76
77 // bottom range segment (min, min) (max, min)
78 q0.first = nodeList_[nodeId].rangeBox_.first.first;
79 q0.second = nodeList_[nodeId].rangeBox_.second.first;
80
81 q1.first = nodeList_[nodeId].rangeBox_.first.second;
82 q1.second = q0.second;
83
84 if(segmentIntersection(p0, p1, q0, q1)) {
85 if(nodeList_[nodeId].childList_.size()) {
86 for(size_t i = 0; i < nodeList_[nodeId].childList_.size(); i++) {
87 rangeSegmentQuery(p0, p1, nodeList_[nodeId].childList_[i], cellList);
88 }
89 return 0;
90 } else {
91 // terminal leaf
92 // return our cells
93 if(debugLevel_ >= static_cast<int>(debug::Priority::VERBOSE)) {
94 this->printMsg("Node #" + std::to_string(nodeId) + " returns its "
95 + std::to_string(nodeList_[nodeId].cellList_.size())
96 + " cells(s).",
98 }
99
100 cellList.insert(cellList.end(), nodeList_[nodeId].cellList_.begin(),
101 nodeList_[nodeId].cellList_.end());
103 return 0;
104 }
105 }
106
107 // right segment (max, min) (max, max)
108 q0.first = nodeList_[nodeId].rangeBox_.first.second;
109 q0.second = nodeList_[nodeId].rangeBox_.second.first;
110
111 q1.first = nodeList_[nodeId].rangeBox_.first.second;
112 q1.second = nodeList_[nodeId].rangeBox_.second.second;
113
114 if(segmentIntersection(p0, p1, q0, q1)) {
115 if(nodeList_[nodeId].childList_.size()) {
116 for(size_t i = 0; i < nodeList_[nodeId].childList_.size(); i++) {
117 rangeSegmentQuery(p0, p1, nodeList_[nodeId].childList_[i], cellList);
118 }
119 return 0;
120 } else {
121 // terminal leaf
122 // return our cells
123 if(debugLevel_ >= static_cast<int>(debug::Priority::VERBOSE)) {
124 this->printMsg("Node #" + std::to_string(nodeId) + " returns its "
125 + std::to_string(nodeList_[nodeId].cellList_.size())
126 + " cells(s).",
128 }
129
130 cellList.insert(cellList.end(), nodeList_[nodeId].cellList_.begin(),
131 nodeList_[nodeId].cellList_.end());
133 return 0;
134 }
135 }
136
137 // top segment (min, max) (max, max)
138 q0.first = nodeList_[nodeId].rangeBox_.first.first;
139 q0.second = nodeList_[nodeId].rangeBox_.second.second;
140
141 q1.first = nodeList_[nodeId].rangeBox_.first.second;
142 q1.second = nodeList_[nodeId].rangeBox_.second.second;
143
144 if(segmentIntersection(p0, p1, q0, q1)) {
145 if(nodeList_[nodeId].childList_.size()) {
146 for(size_t i = 0; i < nodeList_[nodeId].childList_.size(); i++) {
147 rangeSegmentQuery(p0, p1, nodeList_[nodeId].childList_[i], cellList);
148 }
149 return 0;
150 } else {
151 // terminal leaf
152 // return our cells
153 if(debugLevel_ >= static_cast<int>(debug::Priority::VERBOSE)) {
154 this->printMsg("Node #" + std::to_string(nodeId) + " returns its "
155 + std::to_string(nodeList_[nodeId].cellList_.size())
156 + " cells(s).",
158 }
159
160 cellList.insert(cellList.end(), nodeList_[nodeId].cellList_.begin(),
161 nodeList_[nodeId].cellList_.end());
163 return 0;
164 }
165 }
166
167 // left segment (min, min) (min, max)
168 q0.first = nodeList_[nodeId].rangeBox_.first.first;
169 q0.second = nodeList_[nodeId].rangeBox_.second.first;
170
171 q1.first = nodeList_[nodeId].rangeBox_.first.first;
172 q1.second = nodeList_[nodeId].rangeBox_.second.second;
173
174 if(segmentIntersection(p0, p1, q0, q1)) {
175 if(nodeList_[nodeId].childList_.size()) {
176 for(size_t i = 0; i < nodeList_[nodeId].childList_.size(); i++) {
177 rangeSegmentQuery(p0, p1, nodeList_[nodeId].childList_[i], cellList);
178 }
179 return 0;
180 } else {
181 // terminal leaf
182 // return our cells
183 if(debugLevel_ >= static_cast<int>(debug::Priority::VERBOSE)) {
184 this->printMsg("Node #" + std::to_string(nodeId) + " returns its "
185 + std::to_string(nodeList_[nodeId].cellList_.size())
186 + " cells(s).",
188 }
189
190 cellList.insert(cellList.end(), nodeList_[nodeId].cellList_.begin(),
191 nodeList_[nodeId].cellList_.end());
193 return 0;
194 }
195 }
196
197 // is the segment completely included in the range bounding box?
198 if((p0.first >= nodeList_[nodeId].rangeBox_.first.first)
199 && (p0.first < nodeList_[nodeId].rangeBox_.first.second)
200 && (p0.second >= nodeList_[nodeId].rangeBox_.second.first)
201 && (p0.second < nodeList_[nodeId].rangeBox_.second.second)) {
202
203 // p0 is in there
204 if(nodeList_[nodeId].childList_.size()) {
205 for(size_t i = 0; i < nodeList_[nodeId].childList_.size(); i++) {
206 rangeSegmentQuery(p0, p1, nodeList_[nodeId].childList_[i], cellList);
207 }
208 return 0;
209 } else {
210 // terminal leaf
211 // return our cells
212 if(debugLevel_ >= static_cast<int>(debug::Priority::VERBOSE)) {
213 this->printMsg("Node #" + std::to_string(nodeId) + " returns its "
214 + std::to_string(nodeList_[nodeId].cellList_.size())
215 + " cells(s).",
217 }
218
219 cellList.insert(cellList.end(), nodeList_[nodeId].cellList_.begin(),
220 nodeList_[nodeId].cellList_.end());
222 return 0;
223 }
224 }
225 if((p1.first >= nodeList_[nodeId].rangeBox_.first.first)
226 && (p1.first < nodeList_[nodeId].rangeBox_.first.second)
227 && (p1.second >= nodeList_[nodeId].rangeBox_.second.first)
228 && (p1.second < nodeList_[nodeId].rangeBox_.second.second)) {
229
230 // p1 is in there
231 if(nodeList_[nodeId].childList_.size()) {
232 for(size_t i = 0; i < nodeList_[nodeId].childList_.size(); i++) {
233 rangeSegmentQuery(p0, p1, nodeList_[nodeId].childList_[i], cellList);
234 }
235 return 0;
236 } else {
237 // terminal leaf
238 // return our cells
239 if(debugLevel_ >= static_cast<int>(debug::Priority::VERBOSE)) {
240 this->printMsg("Node #" + std::to_string(nodeId) + " returns its "
241 + std::to_string(nodeList_[nodeId].cellList_.size())
242 + " cells(s).",
244 }
245
246 cellList.insert(cellList.end(), nodeList_[nodeId].cellList_.begin(),
247 nodeList_[nodeId].cellList_.end());
249 return 0;
250 }
251 }
252
253 return 0;
254}
255
256int RangeDrivenOctree::statNode(const SimplexId &nodeId, std::ostream &stream) {
257
258 stream << "[RangeDrivenOctree]" << std::endl;
259 stream << "[RangeDrivenOctree] Node #" << nodeId << std::endl;
260 stream << "[RangeDrivenOctree] Domain box: ["
261 << nodeList_[nodeId].domainBox_[0].first << " "
262 << nodeList_[nodeId].domainBox_[0].second << "] ["
263 << nodeList_[nodeId].domainBox_[1].first << " "
264 << nodeList_[nodeId].domainBox_[1].second << "] ["
265 << nodeList_[nodeId].domainBox_[2].first << " "
266 << nodeList_[nodeId].domainBox_[2].second << "] "
267 << " volume="
268 << (nodeList_[nodeId].domainBox_[0].second
269 - nodeList_[nodeId].domainBox_[0].first)
270 * (nodeList_[nodeId].domainBox_[1].second
271 - nodeList_[nodeId].domainBox_[1].first)
272 * (nodeList_[nodeId].domainBox_[2].second
273 - nodeList_[nodeId].domainBox_[2].first)
274 << " threshold=" << leafMinimumDomainVolumeRatio_ * domainVolume_
275 << std::endl;
276 stream << "[RangeDrivenOctree] Range box: ["
277 << nodeList_[nodeId].rangeBox_.first.first << " "
278 << nodeList_[nodeId].rangeBox_.first.second << "] ["
279 << nodeList_[nodeId].rangeBox_.second.first << " "
280 << nodeList_[nodeId].rangeBox_.second.second << "] "
281 << " area="
282 << (nodeList_[nodeId].rangeBox_.first.second
283 - nodeList_[nodeId].rangeBox_.first.first)
284 * (nodeList_[nodeId].rangeBox_.second.second
285 - nodeList_[nodeId].rangeBox_.second.first)
286 << " threshold=" << leafMinimumRangeAreaRatio_ * rangeArea_
287 << std::endl;
288 stream << "[RangeDrivenOctree] Number of cells: "
289 << nodeList_[nodeId].cellList_.size() << std::endl;
290
291 return 0;
292}
293
294int RangeDrivenOctree::stats(std::ostream &stream) {
295
296 SimplexId leafNumber = 0, nonEmptyLeafNumber = 0, minCellNumber = -1,
297 maxCellNumber = -1, storedCellNumber = 0;
298 float averageCellNumber = 0;
299 SimplexId maxCellId = 0;
300
301 for(size_t i = 0; i < nodeList_.size(); i++) {
302 if(!nodeList_[i].childList_.size()) {
303 // leaf
304 leafNumber++;
305 if(nodeList_[i].cellList_.size()) {
306 nonEmptyLeafNumber++;
307 storedCellNumber += nodeList_[i].cellList_.size();
308
309 averageCellNumber += nodeList_[i].cellList_.size();
310 if((minCellNumber == -1)
311 || (nodeList_[i].cellList_.size()
312 < static_cast<size_t>(minCellNumber)))
313 minCellNumber = nodeList_[i].cellList_.size();
314 if((maxCellNumber == -1)
315 || (nodeList_[i].cellList_.size()
316 > static_cast<size_t>(maxCellNumber))) {
317 maxCellNumber = nodeList_[i].cellList_.size();
318 maxCellId = i;
319 }
320 }
321 }
322 }
323 averageCellNumber /= nonEmptyLeafNumber;
324
325 stream << "[RangeDrivenOctree] Domain volume: " << domainVolume_ << std::endl;
326 stream << "[RangeDrivenOctree] Range area: " << rangeArea_ << std::endl;
327 stream << "[RangeDrivenOctree] Leaf number: " << leafNumber << std::endl;
328 stream << "[RangeDrivenOctree] Non empty leaf number: " << nonEmptyLeafNumber
329 << std::endl;
330 stream << "[RangeDrivenOctree] Average cell number: " << averageCellNumber
331 << std::endl;
332 stream << "[RangeDrivenOctree] Min cell number: " << minCellNumber
333 << std::endl;
334 stream << "[RangeDrivenOctree] Max cell number: " << maxCellNumber
335 << std::endl;
336 stream << "[RangeDrivenOctree] Stored cell number: " << storedCellNumber
337 << " [input=" << cellNumber_ << "]" << std::endl;
338
339 stream << "[RangeDrivenOctree] Max-cell nodeId: " << maxCellId << std::endl;
340
341 if(debugLevel_ > 5) {
342 for(size_t i = 0; i < nodeList_.size(); i++) {
343 if(nodeList_[i].cellList_.size())
344 statNode(i, stream);
345 }
346 }
347
348 return 0;
349}
350
352 const std::pair<double, double> &p0,
353 const std::pair<double, double> &p1,
354 const std::pair<double, double> &q0,
355 const std::pair<double, double> &q1) const {
356
357 // NOTE: [q0, q1] is axis aligned
358 bool horizontal = true;
359 double x = 0, y = 0;
360
361 if(q0.first == q1.first)
362 horizontal = false;
363
364 double denP = p1.first - p0.first;
365 if(!denP)
366 denP = DBL_EPSILON;
367
368 double P = (p1.second - p0.second) / denP;
369 if(!P)
370 P = DBL_EPSILON;
371 double const bP = p1.second - P * p1.first;
372
373 if(horizontal) {
374 y = q0.second;
375 x = (y - bP) / P;
376 } else {
377 x = q0.first;
378 y = P * x + bP;
379 }
380
381 // does the intersection lands on the range segment?
382 if(horizontal) {
383 if((x < fmin(q0.first, q1.first)) || (x > fmax(q0.first, q1.first))) {
384 return false;
385 }
386 } else {
387 if((y < fmin(q0.second, q1.second)) || (y > fmax(q0.second, q1.second))) {
388 return false;
389 }
390 }
391
392 // does it lands on the polygon segment's projection?
393 if((x < fmin(p0.first, p1.first) || (x > fmax(p0.first, p1.first)))
394 && ((y < fmin(p0.second, p1.second)
395 || (y > fmax(p0.second, p1.second))))) {
396 return false;
397 }
398
399 return true;
400}
int debugLevel_
Definition Debug.h:379
void setDebugMsgPrefix(const std::string &prefix)
Definition Debug.h:364
int rangeSegmentQuery(const std::pair< double, double > &p0, const std::pair< double, double > &p1, std::vector< SimplexId > &cellList) const
const SimplexId * cellList_
std::vector< std::array< std::pair< float, float >, 3 > > cellDomainBox_
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_
int stats(std::ostream &stream)
int getTet2NodeMap(std::vector< SimplexId > &map, const bool &forSegmentation=false) const
int statNode(const SimplexId &nodeId, std::ostream &stream)
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
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)