TTK
Loading...
Searching...
No Matches
Octree.cpp
Go to the documentation of this file.
1#include <Octree.h>
2
3using namespace std;
4using namespace ttk;
5
7 initialize(t, 1000);
8}
9
10Octree::Octree(const AbstractTriangulation *t, const int k) {
11 initialize(t, k);
12}
13
14Octree::~Octree() = default;
15
16// Initialize the octree with the given triangulation and bucket threshold.
17void Octree::initialize(const AbstractTriangulation *t, const int k) {
18 this->setDebugMsgPrefix("PR Octree");
19 OctreeNode root(1);
20 allNodes_[1] = root;
21 capacity_ = k;
22 triangulation_ = t;
23
24 // find the minimum and maximum coordinate values
25 float mins[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
26 float maxs[3] = {FLT_MIN, FLT_MIN, FLT_MIN};
27 SimplexId vertexNum = t->getNumberOfVertices();
28 for(int i = 0; i < vertexNum; i++) {
29 float coord[3];
30 t->getVertexPoint(i, coord[0], coord[1], coord[2]);
31 for(int j = 0; j < 3; j++) {
32 if(coord[j] < mins[j])
33 mins[j] = coord[j];
34 if(coord[j] > maxs[j])
35 maxs[j] = coord[j];
36 }
37 }
38
39 for(int j = 0; j < 3; j++) {
40 center_[j] = (mins[j] + maxs[j]) / 2.0;
41 size_[j] = (maxs[j] - mins[j]) / 2.0;
42 }
43}
44
49 OctreeNode *root = lookupNode(1);
50 if(root != nullptr && root->vertexIds_.empty() && !root->childExists_) {
51 return true;
52 }
53
54 return false;
55}
56
61 // assert(node->locCode);
62 size_t depth = 0;
63 for(uint32_t lc = node->locCode_; lc != 1; lc >>= 3)
64 depth++;
65 return depth;
66}
67
72 assert(node->locCode_);
73 const uint32_t locCodeParent = node->locCode_ >> 3;
74 return lookupNode(locCodeParent);
75}
76
80void Octree::visitAll(const OctreeNode *node) {
81 if(node == nullptr)
82 return;
83 this->printMsg(to_string(node->locCode_));
84 for(int i = 0; i < 8; i++) {
85 if(node->childExists_ & (1 << i)) {
86 const uint32_t locCodeChild = (node->locCode_ << 3) | i;
87 const OctreeNode *child = lookupNode(locCodeChild);
88 visitAll(child);
89 }
90 }
91}
92
97 int vertexCount = 0;
98 unordered_map<uint32_t, OctreeNode>::iterator it;
99 for(it = allNodes_.begin(); it != allNodes_.end(); it++) {
100 if(it->second.childExists_ && it->second.vertexIds_.size() > 0) {
101 this->printErr("[Octree] WRONG! The internal node "
102 + to_string(it->second.locCode_)
103 + " should not contain any vertices!");
104 return -1;
105 }
106 if(it->second.vertexIds_.size() > 0) {
107 vertexCount += it->second.vertexIds_.size();
108 }
109 }
110 if(vertexCount != vertexNum) {
111 this->printErr("The number of vertices in the tree is "
112 + to_string(vertexCount) + ", which is not equal to "
113 + to_string(vertexNum));
114 return -1;
115 }
116
117 return 0;
118}
119
124 if(vertexId < 0 || vertexId >= triangulation_->getNumberOfVertices()) {
125 return -1;
126 }
127
128 OctreeNode *current = lookupNode(1);
129 if(current == nullptr) {
130 return -2;
131 }
132 uint32_t location = current->locCode_;
133 std::array<float, 3> ncenter{}, nsize{};
134
135 while(current != nullptr && current->childExists_ != 0) {
136 computeCenterSize(location, ncenter, nsize);
137 location = getChildLocation(location, vertexId, ncenter);
138 current = lookupNode(location);
139 }
140
141 if(current == nullptr) {
142 OctreeNode newnode(location);
143 newnode.vertexIds_ = vector<SimplexId>{vertexId};
144 OctreeNode *parent = lookupNode(location >> 3);
145 parent->childExists_ |= (1 << (location & 7));
146 allNodes_[location] = newnode;
147 } else {
148 current->vertexIds_.push_back(vertexId);
149 subdivide(current);
150 }
151 return 0;
152}
153
159 if(cellId < 0 || cellId >= triangulation_->getNumberOfCells()) {
160 return -1;
161 }
162
163 int dim = triangulation_->getCellVertexNumber(cellId);
164 std::array<float, 3> ncenter{}, nsize{};
165 for(int i = 0; i < dim; i++) {
166 SimplexId vertexId{};
167 OctreeNode *current = lookupNode(1);
168
169 triangulation_->getCellVertex(cellId, i, vertexId);
170
171 while(current != nullptr && current->childExists_ != 0) {
172 auto location = current->locCode_;
173 computeCenterSize(location, ncenter, nsize);
174 location = getChildLocation(location, vertexId, ncenter);
175 current = lookupNode(location);
176 }
177
178 if(current == nullptr) {
179 this->printErr("[Octree] insertCell(): Cannot find the vertex id ("
180 + to_string(vertexId) + ") in the tree!");
181 return -1;
182 } else {
183 if(current->cellIds_.empty()) {
184 current->cellIds_ = vector<SimplexId>{cellId};
185 } else if(current->cellIds_.back() != cellId) {
186 current->cellIds_.push_back(cellId);
187 }
188 }
189 }
190
191 return 0;
192}
193
197void Octree::reindex(vector<SimplexId> &vertices,
198 vector<SimplexId> &nodes,
199 vector<SimplexId> &cells) {
200 int totalCells = triangulation_->getNumberOfCells();
201 vector<int> cellMap(totalCells, -1);
202
203 OctreeNode *root = lookupNode(1);
204 int leafCount = 0;
205 int cellCount = 0;
206
207 // use the depth-first search to complete reindexing
208 stack<const OctreeNode *> nodeStack;
209 nodeStack.push(root);
210
211 while(!nodeStack.empty()) {
212 const OctreeNode *topNode = nodeStack.top();
213 nodeStack.pop();
214
215 if(topNode == nullptr) {
216 this->printErr("[Octree] reindex(): shouldn't get here!");
217 break;
218 }
219 if(topNode->childExists_) {
220 for(int i = 0; i < 8; i++) {
221 if(topNode->childExists_ & (1 << i)) {
222 const uint32_t locCodeChild = (topNode->locCode_ << 3) | i;
223 const OctreeNode *child = lookupNode(locCodeChild);
224 nodeStack.push(child);
225 }
226 }
227 } else {
228 vertices.insert(
229 vertices.end(), topNode->vertexIds_.begin(), topNode->vertexIds_.end());
230 vector<SimplexId> tmp(topNode->vertexIds_.size(), leafCount);
231 nodes.insert(nodes.end(), tmp.begin(), tmp.end());
232 leafCount++;
233
234 if(topNode->cellIds_.size() > 0) {
235 for(auto it = topNode->cellIds_.begin(); it != topNode->cellIds_.end();
236 it++) {
237 if(cellMap[*it] == -1) {
238 cellMap[*it] = cellCount++;
239 }
240 }
241 }
242 }
243 }
244
245 int count = 0;
246 cells.resize(totalCells);
247 for(int i = 0; i < totalCells; i++) {
248 if(cellMap[i] == -1) {
249 count++;
250 }
251 cells.at(cellMap[i]) = i;
252 }
253 this->printMsg("reindex(): There are " + to_string(count)
254 + " wrong entries!");
255}
256
257void Octree::computeCenterSize(uint32_t location,
258 std::array<float, 3> &centerArr,
259 std::array<float, 3> &sizeArr) {
260 int leftmost = 0;
261 uint32_t tmp = location;
262 while(tmp >>= 1) {
263 leftmost++;
264 }
265 if(leftmost % 3 != 0) {
266 this->printMsg("Location: " + std::to_string(location)
267 + ", leftmost: " + std::to_string(leftmost));
268 this->printErr("computeCenterSize(): the location seems not correct!");
269 this->printErr("Please try a larger bucket capacity!");
270 return;
271 }
272
273 // initialize the center and size arrays
274 centerArr = center_;
275 sizeArr = size_;
276
277 uint32_t mask;
278 for(int i = leftmost - 1; i >= 0; i -= 3) {
279 sizeArr[0] /= 2.0;
280 sizeArr[1] /= 2.0;
281 sizeArr[2] /= 2.0;
282 for(int j = 0; j < 3; j++) {
283 mask = 1 << (i - j);
284 if(location & mask) {
285 centerArr[j] += sizeArr[j];
286 } else {
287 centerArr[j] -= sizeArr[j];
288 }
289 }
290 }
291}
292
293uint32_t Octree::getChildLocation(uint32_t parLoc,
294 ttk::SimplexId vertexId,
295 const std::array<float, 3> &centerArr) {
296 float xval{}, yval{}, zval{};
297 if(triangulation_->getVertexPoint(vertexId, xval, yval, zval)) {
298 this->printErr("getChildLocation(): FAILED to get the coordinate values "
299 "of the vertex id "
300 + std::to_string(vertexId));
301 return -1;
302 }
303
304 if(xval < centerArr[0]) {
305 if(yval < centerArr[1]) {
306 if(zval < centerArr[2]) {
307 return (parLoc << 3);
308 } else {
309 return (parLoc << 3) + 1;
310 }
311 } else {
312 if(zval < centerArr[2]) {
313 return (parLoc << 3) + 2;
314 } else {
315 return (parLoc << 3) + 3;
316 }
317 }
318 } else {
319 if(yval < centerArr[1]) {
320 if(zval < centerArr[2]) {
321 return (parLoc << 3) + 4;
322 } else {
323 return (parLoc << 3) + 5;
324 }
325 } else {
326 if(zval < centerArr[2]) {
327 return (parLoc << 3) + 6;
328 } else {
329 return (parLoc << 3) + 7;
330 }
331 }
332 }
333
334 this->printErr("getChildLocation(): Shouldn't get here!");
335 return -1;
336}
337
338void Octree::subdivide(OctreeNode *node) {
339 if(node == nullptr)
340 return;
341
342 if((int)node->vertexIds_.size() > capacity_) {
343 uint32_t childCode = 0;
344 std::array<float, 3> ncenter{}, nsize{};
345
346 computeCenterSize(node->locCode_, ncenter, nsize);
347
348 for(int v : node->vertexIds_) {
349 childCode = getChildLocation(node->locCode_, v, ncenter);
350 OctreeNode *childNode = lookupNode(childCode);
351
352 if(childNode == nullptr) {
353 OctreeNode newnode(childCode);
354 newnode.vertexIds_.push_back(v);
355 allNodes_[childCode] = newnode;
356 node->childExists_ |= (1 << (childCode & 7));
357 } else {
358 childNode->vertexIds_.push_back(v);
359 }
360 }
361
362 node->vertexIds_.clear();
363
364 for(uint32_t i = 0; i < 8; i++) {
365 if(node->childExists_ & (1 << i)) {
366 subdivide(lookupNode((node->locCode_ << 3) | i));
367 }
368 }
369 }
370}
uint32_t locCode_
Definition: Octree.h:44
std::vector< ttk::SimplexId > cellIds_
Definition: Octree.h:47
std::vector< ttk::SimplexId > vertexIds_
Definition: Octree.h:46
uint8_t childExists_
Definition: Octree.h:45
Octree(const ttk::AbstractTriangulation *t)
Definition: Octree.cpp:6
bool empty()
Definition: Octree.cpp:48
size_t getNodeTreeDepth(const OctreeNode *node)
Definition: Octree.cpp:60
~Octree() override
void visitAll(const OctreeNode *node)
Definition: Octree.cpp:80
int insertVertex(ttk::SimplexId &vertexId)
Definition: Octree.cpp:123
void initialize(const ttk::AbstractTriangulation *t, const int k)
Definition: Octree.cpp:17
int insertCell(ttk::SimplexId &cellId)
Definition: Octree.cpp:158
void reindex(std::vector< ttk::SimplexId > &vertices, std::vector< ttk::SimplexId > &nodes, std::vector< ttk::SimplexId > &cells)
Definition: Octree.cpp:197
OctreeNode * getParentNode(OctreeNode *node)
Definition: Octree.cpp:71
int verifyTree(ttk::SimplexId &vertexNum)
Definition: Octree.cpp:96
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
virtual int getVertexPoint(const SimplexId &vertexId, float &x, float &y, float &z) const
virtual SimplexId getNumberOfCells() const
virtual SimplexId getNumberOfVertices() const
virtual int getCellVertex(const SimplexId &cellId, const int &localVertexId, SimplexId &vertexId) const
virtual SimplexId getCellVertexNumber(const SimplexId &cellId) const
void setDebugMsgPrefix(const std::string &prefix)
Definition: Debug.h:364
int printMsg(const std::string &msg, const debug::Priority &priority=debug::Priority::INFO, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cout) const
Definition: Debug.h:118
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition: Debug.h:149
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22