TTK
Loading...
Searching...
No Matches
PathCompression.h
Go to the documentation of this file.
1
27
28#pragma once
29
30// base code includes
31#include <Triangulation.h>
32
33using ttk::SimplexId;
34
35namespace ttk {
36 namespace pcp {
37#ifdef TTK_ENABLE_64BIT_IDS
46 constexpr unsigned long long int getHash(const unsigned long long int a,
47 const unsigned long long int b) {
48 return (a * b + (a * a) + (b * b) + (a * a * a) * (b * b * b))
49 % ULONG_LONG_MAX;
50 }
51#else
60 constexpr unsigned int getHash(const unsigned int a, const unsigned int b) {
61 return (a * b + (a * a) + (b * b) + (a * a * a) * (b * b * b)) % UINT_MAX;
62 }
63#endif
64 } // namespace pcp
65 class PathCompression : public virtual Debug {
66 public:
68
74 };
75
81 }
82
96 template <typename triangulationType>
97 inline int execute(OutputSegmentation &outSegmentation,
98 const SimplexId *const orderArray,
99 const triangulationType &triangulation);
100
117 template <typename triangulationType>
118 int computePathCompression(SimplexId *const ascSegmentation,
119 SimplexId *const dscSegmentation,
120 const SimplexId *const orderArray,
121 const triangulationType &triangulation) const;
122
140 template <typename triangulationType>
142 SimplexId *const segmentation,
143 const bool computeAscending,
144 const SimplexId *const orderArray,
145 const triangulationType &triangulation) const;
146
160 template <typename triangulationType>
161 int computeMSHash(SimplexId *const morseSmaleSegmentation,
162 const SimplexId *const ascSegmentation,
163 const SimplexId *const dscSegmentation,
164 const triangulationType &triangulation) const;
165
166 protected:
167 // Compute ascending segmentation?
169
170 // Compute descending segmentation?
172
173 // Compute Morse-Smale segmentation hash?
175 };
176} // namespace ttk
177
178template <typename triangulationType>
180 const SimplexId *const orderArray,
181 const triangulationType &triangulation) {
182 if(orderArray == nullptr)
183 return this->printErr("Input offset field pointer is null.");
184
185 Timer t;
186
187 this->printMsg("Start computing segmentations", 0.0, t.getElapsedTime(),
188 this->threadNumber_);
189
192 computePathCompression(outSegmentation.ascending_,
193 outSegmentation.descending_, orderArray,
194 triangulation);
197 outSegmentation.ascending_, true, orderArray, triangulation);
200 outSegmentation.descending_, false, orderArray, triangulation);
201 }
202
204 computeMSHash(outSegmentation.morseSmale_, outSegmentation.ascending_,
205 outSegmentation.descending_, triangulation);
206 }
207
208 this->printMsg("Data-set ("
209 + std::to_string(triangulation.getNumberOfVertices())
210 + " points) processed",
211 1.0, t.getElapsedTime(), this->threadNumber_);
212
213 return 0;
214}
215
216template <typename triangulationType>
218 SimplexId *const ascSegmentation,
219 SimplexId *const dscSegmentation,
220 const SimplexId *const orderArray,
221 const triangulationType &triangulation) const {
222
223 ttk::Timer localTimer;
224
225 const SimplexId nVertices = triangulation.getNumberOfVertices();
226 std::vector<SimplexId> lActiveVertices;
227
228#ifdef TTK_ENABLE_OPENMP
229#pragma omp parallel num_threads(threadNumber_) private(lActiveVertices)
230 {
231 lActiveVertices.reserve(std::ceil(nVertices / threadNumber_));
232
233#pragma omp for schedule(static)
234#else // TTK_ENABLE_OPENMP
235 lActiveVertices.reserve(nVertices);
236#endif // TTK_ENABLE_OPENMP
237 // find the largest and smallest neighbor for each vertex
238 for(SimplexId i = 0; i < nVertices; i++) {
239 SimplexId neighborId{0};
240 SimplexId numNeighbors = triangulation.getVertexNeighborNumber(i);
241
242 bool hasLargerNeighbor = false;
243 SimplexId &dmi = dscSegmentation[i];
244 dmi = i;
245
246 bool hasSmallerNeighbor = false;
247 SimplexId &ami = ascSegmentation[i];
248 ami = i;
249
250 // check all neighbors
251 for(SimplexId n = 0; n < numNeighbors; n++) {
252 triangulation.getVertexNeighbor(i, n, neighborId);
253
254 if(orderArray[neighborId] < orderArray[ami]) {
255 ami = neighborId;
256 hasSmallerNeighbor = true;
257 } else if(orderArray[neighborId] > orderArray[dmi]) {
258 dmi = neighborId;
259 hasLargerNeighbor = true;
260 }
261 }
262
263 if(hasLargerNeighbor || hasSmallerNeighbor) {
264 lActiveVertices.push_back(i);
265 }
266 }
267
268 size_t lnActiveVertices = lActiveVertices.size();
269 size_t currentIndex = 0;
270
271 // compress paths until no changes occur
272 while(lnActiveVertices > 0) {
273 for(size_t i = 0; i < lnActiveVertices; i++) {
274 SimplexId &v = lActiveVertices[i];
275 SimplexId &vDsc = dscSegmentation[v];
276 SimplexId &vAsc = ascSegmentation[v];
277
278// compress paths
279#ifdef TTK_ENABLE_OPENMP
280#pragma omp atomic read
281#endif // TTK_ENABLE_OPENMP
282 vDsc = dscSegmentation[vDsc];
283
284#ifdef TTK_ENABLE_OPENMP
285#pragma omp atomic read
286#endif // TTK_ENABLE_OPENMP
287 vAsc = ascSegmentation[vAsc];
288
289 // check if fully compressed
290 if(vDsc != dscSegmentation[vDsc] || vAsc != ascSegmentation[vAsc]) {
291 lActiveVertices[currentIndex++] = v;
292 }
293 }
294
295 lnActiveVertices = currentIndex;
296 currentIndex = 0;
297 }
298#ifdef TTK_ENABLE_OPENMP
299 }
300#endif // TTK_ENABLE_OPENMP
301
302 this->printMsg("Asc. and Desc. segmentation computed", 1.0,
303 localTimer.getElapsedTime(), this->threadNumber_,
305
306 return 0; // return success
307}
308
309template <typename triangulationType>
311 SimplexId *const segmentation,
312 const bool computeAscending,
313 const SimplexId *const orderArray,
314 const triangulationType &triangulation) const {
315
316 ttk::Timer localTimer;
317
318 const SimplexId nVertices = triangulation.getNumberOfVertices();
319 std::vector<SimplexId> lActiveVertices;
320
321#ifdef TTK_ENABLE_OPENMP
322#pragma omp parallel num_threads(threadNumber_)
323 {
324 lActiveVertices.reserve(std::ceil(nVertices / threadNumber_));
325
326#pragma omp for schedule(static)
327#else // TTK_ENABLE_OPENMP
328 lActiveVertices.reserve(nVertices);
329#endif // TTK_ENABLE_OPENMP
330 // find the largest neighbor for each vertex
331 for(SimplexId i = 0; i < nVertices; i++) {
332 SimplexId neighborId{0};
333 SimplexId numNeighbors = triangulation.getVertexNeighborNumber(i);
334
335 bool hasLargerNeighbor = false;
336 SimplexId &mi = segmentation[i];
337 mi = i;
338
339 // check all neighbors
340 for(SimplexId n = 0; n < numNeighbors; n++) {
341 triangulation.getVertexNeighbor(i, n, neighborId);
342
343 if(computeAscending) {
344 if(orderArray[neighborId] < orderArray[mi]) {
345 mi = neighborId;
346 hasLargerNeighbor = true;
347 }
348 } else {
349 if(orderArray[neighborId] > orderArray[mi]) {
350 mi = neighborId;
351 hasLargerNeighbor = true;
352 }
353 }
354 }
355
356 if(hasLargerNeighbor) {
357 lActiveVertices.push_back(i);
358 }
359 }
360
361 size_t lnActiveVertices = lActiveVertices.size();
362 size_t currentIndex = 0;
363
364 // compress paths until no changes occur
365 while(lnActiveVertices > 0) {
366 for(size_t i = 0; i < lnActiveVertices; i++) {
367 SimplexId &v = lActiveVertices[i];
368 SimplexId &vMan = segmentation[v];
369
370// compress paths
371#ifdef TTK_ENABLE_OPENMP
372#pragma omp atomic read
373#endif // TTK_ENABLE_OPENMP
374 vMan = segmentation[vMan];
375
376 // check if fully compressed
377 if(vMan != segmentation[vMan]) {
378 lActiveVertices[currentIndex++] = v;
379 }
380 }
381
382 lnActiveVertices = currentIndex;
383 currentIndex = 0;
384 }
385#ifdef TTK_ENABLE_OPENMP
386 }
387#endif // TTK_ENABLE_OPENMP
388
389 if(computeAscending) {
390 this->printMsg("Ascending segmentation computed", 1.0,
391 localTimer.getElapsedTime(), this->threadNumber_,
393 } else {
394 this->printMsg("Descending segmentation computed", 1.0,
395 localTimer.getElapsedTime(), this->threadNumber_,
397 }
398
399 return 0; // return success
400}
401
402template <typename triangulationType>
404 SimplexId *const morseSmaleSegmentation,
405 const SimplexId *const ascSegmentation,
406 const SimplexId *const dscSegmentation,
407 const triangulationType &triangulation) const {
408
409 ttk::Timer localTimer;
410
411 const size_t nVerts = triangulation.getNumberOfVertices();
412
413#ifdef TTK_ENABLE_OPENMP
414#pragma omp parallel for schedule(static) num_threads(threadNumber_)
415#endif // TTK_ENABLE_OPENMP
416 for(size_t i = 0; i < nVerts; ++i) {
417 morseSmaleSegmentation[i]
418 = pcp::getHash(ascSegmentation[i], dscSegmentation[i]);
419 }
420
421 this->printMsg("Morse-Smale segmentation hash computed", 1.0,
422 localTimer.getElapsedTime(), this->threadNumber_,
424
425 return 0;
426}
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
Minimalist debugging class.
Definition: Debug.h:88
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
TTK processing package for the computation of Morse-Smale segmentations using Path Compression.
void preconditionTriangulation(AbstractTriangulation *const data)
int execute(OutputSegmentation &outSegmentation, const SimplexId *const orderArray, const triangulationType &triangulation)
Main function for computing the Morse-Smale complex.
int computeMSHash(SimplexId *const morseSmaleSegmentation, const SimplexId *const ascSegmentation, const SimplexId *const dscSegmentation, const triangulationType &triangulation) const
Computes a MS segmentation hash.
int computePathCompressionSingle(SimplexId *const segmentation, const bool computeAscending, const SimplexId *const orderArray, const triangulationType &triangulation) const
Compute the ascending or descending segmentation.
int computePathCompression(SimplexId *const ascSegmentation, SimplexId *const dscSegmentation, const SimplexId *const orderArray, const triangulationType &triangulation) const
Compute the ascending and descending segmentation in one run.
double getElapsedTime()
Definition: Timer.h:15
constexpr unsigned int getHash(const unsigned int a, const unsigned int b)
Get a hash value from two keys.
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
Pointers to pre-allocated segmentation point data arrays.
printMsg(debug::output::GREEN+"                           "+debug::output::ENDCOLOR+debug::output::GREEN+"▒"+debug::output::ENDCOLOR+debug::output::GREEN+"▒▒▒▒▒▒▒▒▒▒▒▒▒░"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)