TTK
Loading...
Searching...
No Matches
MergeTreePrincipalGeodesicsDecoding.h
Go to the documentation of this file.
1
18
19#pragma once
20
21// ttk common includes
22#include <Debug.h>
24
25namespace ttk {
26
33 : virtual public Debug,
35
36 protected:
40
41 std::vector<std::vector<double *>> pVS_, pV2s_, pTrees2Vs_, pTrees2V2s_;
42 size_t vSize_, vSize2_;
44
45 // Filled
46 std::vector<std::vector<double>> tEllipses_, tRectangle_, tSurface_;
47 std::vector<std::vector<std::vector<double>>> tGeodesics_;
48
49 std::vector<double> geodesicsDistances_; // distance between extremities
50
51 std::vector<bool> surfaceIsBoundary_;
52 std::vector<int> surfaceBoundaryID_;
53
54 public:
56
57 //----------------------------------------------------------------------------
58 // Utils
59 //----------------------------------------------------------------------------
60 template <class dataType>
62 if(not isPersistenceDiagram_) {
63 bool const useMinMax = true;
64 bool const cleanTree = false;
65 bool const pt = 0.0;
66 std::vector<int> nodeCorr;
67 preprocessingPipeline<dataType>(barycenter, 0.0, 100.0, 100.0,
68 branchDecomposition_, useMinMax,
69 cleanTree, pt, nodeCorr, false);
70 }
71 }
72
73 template <class dataType>
75 std::vector<ttk::ftm::MergeTree<dataType>> &inputTrees) {
76 preprocessingTrees<dataType>(inputTrees, treesNodeCorr_);
77#ifdef TTK_ENABLE_OPENMP
78#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
79#endif
80 for(unsigned int i = 0; i < inputTrees.size(); ++i)
81 postprocessingPipeline<dataType>(&(inputTrees[i].tree));
82 }
83
84 template <class dataType>
86 std::vector<std::vector<double *>> &vS,
87 std::vector<std::vector<double *>> &v2s,
88 size_t vSize,
89 std::array<double, 2> &middle) {
90 int cptDivide = 0;
91 std::vector<double> alpha(2, 0.0);
92 for(unsigned int i = 0; i < 2; ++i) {
93 cptDivide = 0;
94 for(unsigned int j = 0; j < vSize; ++j) {
95 if(barycenter.tree.isNodeAlone(j))
96 continue;
97 for(unsigned int k = 0; k < 2; ++k) {
98 if(std::abs(v2s[i][k][j]) < Geometry::pow(10.0, -DBL_DIG))
99 continue;
100 alpha[i] += (vS[i][k][j] / v2s[i][k][j]);
101 ++cptDivide;
102 }
103 }
104 alpha[i] /= cptDivide;
105 }
106
107 for(unsigned int i = 0; i < 2; ++i)
108 middle[i] = alpha[i] / (1 + alpha[i]);
109 }
110
111 template <class dataType>
113 std::vector<ttk::ftm::MergeTree<dataType>> &barycenters) {
114 // Preprocessing
115 preprocessBarycenter<dataType>(barycenters[0]);
116 if(barycenters.size() > 1)
117 preprocessBarycenter<dataType>(barycenters[1]);
118
119 // Compute
120 geodesicsDistances_.resize(pVS_.size());
121#ifdef TTK_ENABLE_OPENMP
122#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
123#endif
124 for(unsigned int i = 0; i < pVS_.size(); ++i) {
125 ftm::MergeTree<dataType> extremityV1, extremityV2;
126 getInterpolation<dataType>(
127 barycenters[0], pVS_[i], pV2s_[i], vSize_, 0.0, extremityV1);
128 getInterpolation<dataType>(
129 barycenters[0], pVS_[i], pV2s_[i], vSize_, 1.0, extremityV2);
130 // Get distance
131 dataType distance;
133 extremityV1, extremityV2, distance, true, useDoubleInput_);
134 if(barycenters.size() > 1) {
135 ftm::MergeTree<dataType> extremityV1_2, extremityV2_2;
136 getInterpolation<dataType>(barycenters[1], pTrees2Vs_[i],
137 pTrees2V2s_[i], vSize2_, 0.0,
138 extremityV1_2);
139 getInterpolation<dataType>(barycenters[1], pTrees2Vs_[i],
140 pTrees2V2s_[i], vSize2_, 1.0,
141 extremityV2_2);
142 // Get distance
143 dataType distance2;
144 computeOneDistance(extremityV1_2, extremityV2_2, distance2, true,
145 useDoubleInput_, false);
146 distance = mixDistances(distance, distance2);
147 }
148 geodesicsDistances_[i] = distance;
149 }
150
151 // Postprocessing
152 postprocessingPipeline<dataType>(&(barycenters[0].tree));
153 if(barycenters.size() > 1)
154 postprocessingPipeline<dataType>(&(barycenters[1].tree));
155 }
156
157 //----------------------------------------------------------------------------
158 // Construct functions
159 //----------------------------------------------------------------------------
160 template <class dataType>
162 ftm::MergeTree<dataType> &barycenter,
163 std::vector<ftm::MergeTree<dataType>> &inputTrees,
164 std::vector<ftm::MergeTree<dataType>> &reconstructedTrees,
165 std::vector<double> &reconstructionErrors,
166 std::vector<std::vector<std::tuple<ftm::idNode, ftm::idNode, double>>>
167 &recInputMatchings,
168 std::vector<std::vector<std::tuple<ftm::idNode, ftm::idNode, double>>>
169 &recBaryMatchings,
170 bool isSecondInput = false) {
171 auto &vSToUse = (isSecondInput ? pTrees2Vs_ : pVS_);
172 auto &v2sToUse = (isSecondInput ? pTrees2V2s_ : pV2s_);
173 auto vSizeToUse = (isSecondInput ? vSize2_ : vSize_);
174
175 // Preprocessing
176 preprocessBarycenter<dataType>(barycenter);
177 if(inputTrees.size() != 0)
178 preprocessingTrees<dataType>(inputTrees);
179
180 // Reconstruction
181 reconstructedTrees.resize(allTreesTs_.size());
183 recBaryMatchings.resize(reconstructedTrees.size());
184#ifdef TTK_ENABLE_OPENMP
185#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
186#endif
187 for(unsigned int i = 0; i < reconstructedTrees.size(); ++i) {
188 getMultiInterpolation<dataType>(barycenter, vSToUse, v2sToUse,
189 vSizeToUse, allTreesTs_[i],
190 reconstructedTrees[i]);
192 dataType distance;
193 computeOneDistance<dataType>(reconstructedTrees[i], barycenter,
194 recBaryMatchings[i], distance, true);
195 }
196 }
197
198 // Compute reconstruction error (if input trees are provided)
199 if(inputTrees.size() != 0
201 auto reconstructionError = computeReconstructionError(
202 barycenter, inputTrees, vSToUse, v2sToUse, vSizeToUse, allTreesTs_,
203 reconstructionErrors, recInputMatchings);
205 std::stringstream ss;
206 ss << "Reconstruction Error = " << reconstructionError;
207 printMsg(ss.str());
208 }
209 }
210
211 // Postprocessing
212 postprocessingPipeline<dataType>(&(barycenter.tree));
213#ifdef TTK_ENABLE_OPENMP
214#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
215#endif
216 for(unsigned int i = 0; i < reconstructedTrees.size(); ++i)
217 postprocessingPipeline<dataType>(&(reconstructedTrees[i].tree));
218#ifdef TTK_ENABLE_OPENMP
219#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
220#endif
221 for(unsigned int i = 0; i < inputTrees.size(); ++i)
222 postprocessingPipeline<dataType>(&(inputTrees[i].tree));
223
224 if(inputTrees.size() != 0 and transferInputTreesInformation_) {
225 for(unsigned int i = 0; i < inputTrees.size(); ++i)
226 convertBranchDecompositionMatching<dataType>(
227 &(reconstructedTrees[i].tree), &(inputTrees[i].tree),
228 recInputMatchings[i]);
229 }
231 for(unsigned int i = 0; i < reconstructedTrees.size(); ++i)
232 convertBranchDecompositionMatching<dataType>(
233 &(reconstructedTrees[i].tree), &(barycenter.tree),
234 recBaryMatchings[i]);
235 }
236
237 template <class dataType>
239 ftm::MergeTree<dataType> &barycenter,
240 std::vector<std::vector<ftm::MergeTree<dataType>>> &geodesicsTrees,
241 bool isSecondInput = false) {
242 auto &vSToUse = (isSecondInput ? pTrees2Vs_ : pVS_);
243 auto &v2sToUse = (isSecondInput ? pTrees2V2s_ : pV2s_);
244 auto vSizeToUse = (isSecondInput ? vSize2_ : vSize_);
245
246 // Preprocessing
247 preprocessBarycenter<dataType>(barycenter);
248
249 std::array<double, 2> middle;
250 getGeodesicsMiddle<dataType>(
251 barycenter, vSToUse, v2sToUse, vSizeToUse, middle);
252
253 // Construct geodesics trees
254 geodesicsTrees.resize(
255 vSToUse.size(), std::vector<ftm::MergeTree<dataType>>(k_));
256 tGeodesics_.resize(
257 geodesicsTrees.size(),
258 std::vector<std::vector<double>>(
259 geodesicsTrees[0].size(), std::vector<double>(vSToUse.size(), 0.0)));
260#ifdef TTK_ENABLE_OPENMP
261#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
262#endif
263 for(unsigned int i = 0; i < geodesicsTrees.size(); ++i)
264 for(unsigned int j = 0; j < geodesicsTrees[i].size(); ++j) {
265 double t = 1.0 / (k_ - 1) * j;
266 getInterpolation<dataType>(barycenter, vSToUse[i], v2sToUse[i],
267 vSizeToUse, t, geodesicsTrees[i][j]);
268 tGeodesics_[i][j][i] = t;
269 int const i2 = (i + 1) % 2;
270 tGeodesics_[i][j][i2] = middle[i2];
271 }
272
273 // Postprocessing
274 postprocessingPipeline<dataType>(&(barycenter.tree));
275#ifdef TTK_ENABLE_OPENMP
276#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
277#endif
278 for(unsigned int i = 0; i < geodesicsTrees.size(); ++i)
279 for(unsigned int j = 0; j < geodesicsTrees[i].size(); ++j)
280 postprocessingPipeline<dataType>(&(geodesicsTrees[i][j].tree));
281 }
282
283 template <class dataType>
285 ftm::MergeTree<dataType> &barycenter,
286 std::vector<ftm::MergeTree<dataType>> &geodesicsEllipses,
287 bool isSecondInput = false) {
288 auto &vSToUse = (isSecondInput ? pTrees2Vs_ : pVS_);
289 auto &v2sToUse = (isSecondInput ? pTrees2V2s_ : pV2s_);
290 auto vSizeToUse = (isSecondInput ? vSize2_ : vSize_);
291
292 // Preprocessing
293 preprocessBarycenter<dataType>(barycenter);
294
295 // Init
296 unsigned int noSample = k_ * 2;
297 geodesicsEllipses.resize(noSample);
298 tEllipses_.resize(geodesicsEllipses.size());
299
300 std::vector<std::vector<double *>> vS(2), v2s(2);
301 vS[0] = vSToUse[0];
302 vS[1] = vSToUse[1];
303 v2s[0] = v2sToUse[0];
304 v2s[1] = v2sToUse[1];
305
306 // Get middle of geodesic
307 std::array<double, 2> middle;
308 getGeodesicsMiddle<dataType>(barycenter, vS, v2s, vSizeToUse, middle);
309
310 // Get ellipses
311#ifdef TTK_ENABLE_OPENMP
312#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
313#endif
314 for(unsigned int i = 0; i < noSample; ++i) {
315 double const angle = 360.0 / noSample * i;
316 double const pi = M_PI;
317 double const radius = 1.0;
318 double x = -1 * radius * std::cos(-1 * angle * pi / 180);
319 double y = -1 * radius * std::sin(-1 * angle * pi / 180);
320
321 // 0: upper-left ; 1: upper-right ; 2: bottom-right ; 3: bottom-left
322 int const quadrant = (x < 0.0 ? (y > 0.0 ? 0 : 3) : (y > 0.0 ? 1 : 2));
323 if(quadrant == 0 or quadrant == 3)
324 x = (x + 1.0) * middle[0];
325 if(quadrant == 1 or quadrant == 2)
326 x = x * (1.0 - middle[0]) + middle[0];
327 if(quadrant == 0 or quadrant == 1)
328 y = y * (1.0 - middle[1]) + middle[1];
329 if(quadrant == 2 or quadrant == 3)
330 y = (y + 1.0) * middle[1];
331
332 // Get interpolation
333 if(x > 1.0 or y > 1.0 or x < 0.0 or y < 0.0)
334 printErr("[constructGeodesicsEllipses] extrapolation.");
335 std::vector<double> ts{x, y};
336 getMultiInterpolation<dataType>(
337 barycenter, vS, v2s, vSizeToUse, ts, geodesicsEllipses[i]);
338 tEllipses_[i] = ts;
339 }
340
341 // Postprocessing
342 postprocessingPipeline<dataType>(&(barycenter.tree));
343#ifdef TTK_ENABLE_OPENMP
344#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
345#endif
346 for(unsigned int i = 0; i < geodesicsEllipses.size(); ++i)
347 postprocessingPipeline<dataType>(&(geodesicsEllipses[i].tree));
348 }
349
350 unsigned int getNumberOfRectangles(unsigned int rectangleMultiplier = 1) {
351 return k_ * 4 * rectangleMultiplier;
352 }
353
354 template <class dataType>
356 ftm::MergeTree<dataType> &barycenter,
357 std::vector<ftm::MergeTree<dataType>> &geodesicsRectangle,
358 unsigned int rectangleMultiplier = 1,
359 bool isSecondInput = false) {
360 auto &vSToUse = (isSecondInput ? pTrees2Vs_ : pVS_);
361 auto &v2sToUse = (isSecondInput ? pTrees2V2s_ : pV2s_);
362 auto vSizeToUse = (isSecondInput ? vSize2_ : vSize_);
363
364 // Preprocessing
365 preprocessBarycenter<dataType>(barycenter);
366
367 // Init
368 unsigned int noSample = getNumberOfRectangles(rectangleMultiplier);
369 geodesicsRectangle.resize(noSample);
370 tRectangle_.resize(geodesicsRectangle.size());
371
372 std::vector<std::vector<double *>> vS(2), v2s(2);
373 vS[0] = vSToUse[0];
374 vS[1] = vSToUse[1];
375 v2s[0] = v2sToUse[0];
376 v2s[1] = v2sToUse[1];
377
378 // Get ellipses
379#ifdef TTK_ENABLE_OPENMP
380#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
381#endif
382 for(unsigned int i = 0; i < noSample; ++i) {
383 // 0: left ; 1: up; 2: right ; 3: bottom
384 int const quadrant = i / (noSample / 4);
385
386 double x = 0.0, y = 0.0;
387 double const offset = (i % (noSample / 4)) / ((noSample / 4.0) - 1.0);
388 switch(quadrant) {
389 case 0: // Left
390 x = 0.0;
391 y = 0.0 + offset;
392 break;
393 case 1: // Up
394 x = 0.0 + offset;
395 y = 1.0;
396 break;
397 case 2: // Right
398 x = 1.0;
399 y = 1.0 - offset;
400 break;
401 case 3: // Bottom
402 default:
403 x = 1.0 - offset;
404 y = 0.0;
405 }
406
407 // Get interpolation
408 std::vector<double> ts{x, y};
409 getMultiInterpolation<dataType>(
410 barycenter, vS, v2s, vSizeToUse, ts, geodesicsRectangle[i]);
411 tRectangle_[i] = ts;
412 }
413
414 // Postprocessing
415 postprocessingPipeline<dataType>(&(barycenter.tree));
416#ifdef TTK_ENABLE_OPENMP
417#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
418#endif
419 for(unsigned int i = 0; i < geodesicsRectangle.size(); ++i)
420 postprocessingPipeline<dataType>(&(geodesicsRectangle[i].tree));
421 }
422
423 template <class dataType>
425 ftm::MergeTree<dataType> &barycenter,
426 std::vector<ftm::MergeTree<dataType>> &geodesicsSurface,
427 bool isSecondInput = false) {
428 auto &vSToUse = (isSecondInput ? pTrees2Vs_ : pVS_);
429 auto &v2sToUse = (isSecondInput ? pTrees2V2s_ : pV2s_);
430 auto vSizeToUse = (isSecondInput ? vSize2_ : vSize_);
431
432 // Preprocessing
433 preprocessBarycenter<dataType>(barycenter);
434
435 // Init
436 unsigned int noSample = k_ * k_;
437 geodesicsSurface.resize(noSample);
438 tSurface_.resize(geodesicsSurface.size());
439 surfaceIsBoundary_.resize(geodesicsSurface.size(), false);
440 surfaceBoundaryID_.resize(geodesicsSurface.size(), -1);
441
442 std::vector<std::vector<double *>> vS(2), v2s(2);
443 vS[0] = vSToUse[0];
444 vS[1] = vSToUse[1];
445 v2s[0] = v2sToUse[0];
446 v2s[1] = v2sToUse[1];
447
448 // Get surface
449#ifdef TTK_ENABLE_OPENMP
450#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
451#endif
452 for(unsigned int i = 0; i < k_; ++i) {
453 for(unsigned int j = 0; j < k_; ++j) {
454 double const x = 1.0 / (k_ - 1) * i;
455 double const y = 1.0 / (k_ - 1) * j;
456 int const index = i * k_ + j;
457 // Get Boundary ID
458 if(i == 0 or j == 0 or i == k_ - 1 or j == k_ - 1) {
459 surfaceIsBoundary_[index] = true;
460 int boundaryID = 0;
461 if(i == 0)
462 boundaryID = j;
463 else if(j == k_ - 1)
464 boundaryID = k_ + i;
465 else if(i == k_ - 1)
466 boundaryID = k_ * 2 + (k_ - j);
467 else // if(j == 0)
468 boundaryID = k_ * 3 + (k_ - i);
469 surfaceBoundaryID_[index] = boundaryID;
470 }
471 // Get interpolation
472 std::vector<double> ts{x, y};
473 getMultiInterpolation<dataType>(
474 barycenter, vS, v2s, vSizeToUse, ts, geodesicsSurface[index]);
475 tSurface_[index] = ts;
476 }
477 }
478
479 // Postprocessing
480 postprocessingPipeline<dataType>(&(barycenter.tree));
481#ifdef TTK_ENABLE_OPENMP
482#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
483#endif
484 for(unsigned int i = 0; i < geodesicsSurface.size(); ++i)
485 postprocessingPipeline<dataType>(&(geodesicsSurface[i].tree));
486 }
487 }; // MergeTreePrincipalGeodesicsDecoding class
488
489} // namespace ttk
#define M_PI
Definition Os.h:50
Minimalist debugging class.
Definition Debug.h:88
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition Debug.h:149
void computeOneDistance(ftm::MergeTree< dataType > &tree1, ftm::MergeTree< dataType > &tree2, std::vector< std::tuple< ftm::idNode, ftm::idNode, double > > &matching, dataType &distance, bool isCalled=false, bool useDoubleInput=false, bool isFirstInput=true)
double mixDistances(dataType distance1, dataType distance2)
std::vector< std::vector< int > > treesNodeCorr_
std::vector< std::vector< double > > allTreesTs_
dataType computeReconstructionError(ftm::MergeTree< dataType > &barycenter, std::vector< ftm::MergeTree< dataType > > &inputTrees, std::vector< std::vector< double * > > &vS, std::vector< std::vector< double * > > &v2s, size_t vSize, std::vector< std::vector< double > > &allTreesTs, std::vector< double > &reconstructionErrors, std::vector< std::vector< std::tuple< ftm::idNode, ftm::idNode, double > > > &matchings, bool transposeVector=true)
void computeGeodesicsDistance(std::vector< ttk::ftm::MergeTree< dataType > > &barycenters)
std::vector< std::vector< std::vector< double > > > tGeodesics_
void constructGeodesicsEllipses(ftm::MergeTree< dataType > &barycenter, std::vector< ftm::MergeTree< dataType > > &geodesicsEllipses, bool isSecondInput=false)
void constructGeodesicsSurface(ftm::MergeTree< dataType > &barycenter, std::vector< ftm::MergeTree< dataType > > &geodesicsSurface, bool isSecondInput=false)
unsigned int getNumberOfRectangles(unsigned int rectangleMultiplier=1)
void reconstruction(ftm::MergeTree< dataType > &barycenter, std::vector< ftm::MergeTree< dataType > > &inputTrees, std::vector< ftm::MergeTree< dataType > > &reconstructedTrees, std::vector< double > &reconstructionErrors, std::vector< std::vector< std::tuple< ftm::idNode, ftm::idNode, double > > > &recInputMatchings, std::vector< std::vector< std::tuple< ftm::idNode, ftm::idNode, double > > > &recBaryMatchings, bool isSecondInput=false)
void preprocessBarycenter(ftm::MergeTree< dataType > &barycenter)
void constructGeodesicsRectangle(ftm::MergeTree< dataType > &barycenter, std::vector< ftm::MergeTree< dataType > > &geodesicsRectangle, unsigned int rectangleMultiplier=1, bool isSecondInput=false)
void constructGeodesicsTrees(ftm::MergeTree< dataType > &barycenter, std::vector< std::vector< ftm::MergeTree< dataType > > > &geodesicsTrees, bool isSecondInput=false)
void processInputTrees(std::vector< ttk::ftm::MergeTree< dataType > > &inputTrees)
void getGeodesicsMiddle(ftm::MergeTree< dataType > &barycenter, std::vector< std::vector< double * > > &vS, std::vector< std::vector< double * > > &v2s, size_t vSize, std::array< double, 2 > &middle)
bool isNodeAlone(idNode nodeId)
T1 pow(const T1 val, const T2 n)
Definition Geometry.h:454
The Topology ToolKit.
ftm::FTMTree_MT tree
Definition FTMTree_MT.h:901
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)