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