TTK
Loading...
Searching...
No Matches
MergeTreePrincipalGeodesicsBase.h
Go to the documentation of this file.
1
10
11#pragma once
12
14#include <cmath>
15
16#define POINTS_BELOW_DIAG_TOLERANCE 1e-4
17
18namespace ttk {
19
20 class MergeTreePrincipalGeodesicsBase : virtual public Debug,
22 protected:
23 // Filled by the algorithm
24 std::vector<std::vector<std::vector<double>>> vS_, v2s_, trees2Vs_,
26 std::vector<std::vector<double>> allTs_, allScaledTs_, allTreesTs_;
27 std::vector<std::vector<double>> branchesCorrelationMatrix_,
29 std::vector<std::vector<std::tuple<ftm::idNode, ftm::idNode, double>>>
31
32 public:
35 "MergeTreePrincipalGeodesicsBase"); // inherited from Debug: prefix will
36 // be printed at the beginning of
37 // every msg
38 }
39
40 //----------------------------------------------------------------------------
41 // Matching / Distance
42 //----------------------------------------------------------------------------
43 template <class dataType>
45 ftm::MergeTree<dataType> &barycenter,
46 std::vector<ftm::MergeTree<dataType>> &inputTrees,
47 std::vector<std::vector<double *>> &vS,
48 std::vector<std::vector<double *>> &v2s,
49 size_t vSize,
50 std::vector<std::vector<double>> &allTreesTs,
51 std::vector<double> &reconstructionErrors,
52 std::vector<std::vector<std::tuple<ftm::idNode, ftm::idNode, double>>>
53 &matchings,
54 bool transposeVector = true) {
55 reconstructionErrors.resize(inputTrees.size());
56 matchings.resize(inputTrees.size());
57 dataType reconstructionError = 0.0;
58#ifdef TTK_ENABLE_OPENMP
59#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_) reduction(+:reconstructionError)
60#endif
61 for(unsigned int i = 0; i < inputTrees.size(); ++i) {
62 ftm::MergeTree<dataType> reconstructedTree;
63 getMultiInterpolation<dataType>(barycenter, vS, v2s, vSize,
64 allTreesTs[i], reconstructedTree,
65 transposeVector);
66 dataType error;
67 computeOneDistance<dataType>(
68 reconstructedTree, inputTrees[i], matchings[i], error, true);
69 reconstructionError += error / inputTrees.size();
70 reconstructionErrors[i] = error;
71 }
72 return reconstructionError;
73 }
74
75 template <class dataType>
77 ftm::MergeTree<dataType> &barycenter,
78 std::vector<ftm::MergeTree<dataType>> &inputTrees,
79 std::vector<std::vector<std::vector<double>>> &vS,
80 std::vector<std::vector<std::vector<double>>> &v2s,
81 std::vector<std::vector<double>> &allTreesTs,
82 std::vector<double> &reconstructionErrors) {
83 std::vector<std::vector<double *>> pVS, pV2s;
86 std::vector<std::vector<std::tuple<ftm::idNode, ftm::idNode, double>>>
87 matchings;
88 return computeReconstructionError(barycenter, inputTrees, pVS, pV2s,
89 vS[0].size(), allTreesTs,
90 reconstructionErrors, matchings, false);
91 }
92
93 template <class dataType>
95 ftm::MergeTree<dataType> &barycenter,
96 std::vector<ftm::MergeTree<dataType>> &inputTrees,
97 std::vector<std::vector<std::vector<double>>> &vS,
98 std::vector<std::vector<std::vector<double>>> &v2s,
99 std::vector<std::vector<double>> &allTreesTs) {
100 std::vector<double> reconstructionErrors;
102 barycenter, inputTrees, vS, v2s, allTreesTs, reconstructionErrors);
103 }
104
105 //----------------------------------------------------------------------------
106 // Interpolation
107 //----------------------------------------------------------------------------
108 double getGeodesicVectorMiddle(std::vector<double> &v,
109 std::vector<double> &v2) {
110 std::vector<double> vProj, v2Proj;
111 ttk::Geometry::addVectorsProjection(v, v2, vProj, v2Proj);
112
113 double alpha = 0.0;
114 int cptDivide = 0;
115 for(unsigned int i = 0; i < vProj.size(); ++i) {
116 if(std::abs(v2Proj[i]) < Geometry::pow(10.0, -DBL_DIG))
117 continue;
118 alpha += (vProj[i] / v2Proj[i]);
119 ++cptDivide;
120 }
121 alpha /= cptDivide;
122 double m = alpha / (1 + alpha);
123 return m;
124 }
125
126 double getAdjustedTMax(double tMin, double m) {
127 if(std::isnan(m))
128 return tMin;
129 return (m - tMin) * (1 - m) / m + m;
130 }
131
132 double getAdjustedTMin(double tMax, double m) {
133 if(std::isnan(m))
134 return tMax;
135 return (m - tMax) * m / (1 - m) + m;
136 }
137
139 double newT, double m, double &tMin, double &tMax, bool updateTMin) {
140 if(updateTMin) {
141 tMin = std::max(newT, tMin);
142 auto adjustedTMax = getAdjustedTMax(newT, m);
143 if(adjustedTMax > tMin)
144 tMax = std::min(adjustedTMax, tMax);
145 } else {
146 tMax = std::min(newT, tMax);
147 auto adjustedTMin = getAdjustedTMin(newT, m);
148 if(adjustedTMin < tMax)
149 tMin = std::max(adjustedTMin, tMin);
150 }
151 }
152
153 template <class dataType>
154 bool adjustDiagonalT(dataType baryBirth,
155 dataType baryDeath,
156 ftm::idNode node,
157 std::vector<std::vector<double>> &vNew,
158 std::vector<std::vector<double>> &v2New,
159 double &tMin,
160 double &tMax) {
161 bool shortener = false;
162
163 // Compute V1 extremity
164 auto newBirthV1 = baryBirth - vNew[node][0];
165 auto newDeathV1 = baryDeath - vNew[node][1];
166
167 // Compute V2 extremity
168 auto newBirthV2 = baryBirth + v2New[node][0];
169 auto newDeathV2 = baryDeath + v2New[node][1];
170
171 // Adjust t and v
172 if(newBirthV1 > newDeathV1 and newBirthV2 > newDeathV2) {
173 shortener = true;
174
175 double divisor
176 = (vNew[node][1] - vNew[node][0]) / (baryDeath - baryBirth);
177 vNew[node][0] /= divisor;
178 vNew[node][1] /= divisor;
179 double divisor2
180 = (v2New[node][0] - v2New[node][1]) / (baryDeath - baryBirth);
181 v2New[node][0] /= divisor2;
182 v2New[node][1] /= divisor2;
183 } else if(newBirthV1 > newDeathV1 or newBirthV2 > newDeathV2) {
184 double newT = (baryDeath - vNew[node][1] - baryBirth + vNew[node][0])
185 / (vNew[node][0] + v2New[node][0]
186 - (vNew[node][1] + v2New[node][1]));
187 double m = getGeodesicVectorMiddle(vNew[node], v2New[node]);
188 if(newBirthV1 > newDeathV1)
189 updateT(newT, m, tMin, tMax, true);
190 if(newBirthV2 > newDeathV2)
191 updateT(newT, m, tMin, tMax, false);
192 }
193
194 return shortener;
195 }
196
197 template <class dataType>
199 dataType baryBirth,
200 dataType baryDeath,
201 ftm::idNode node,
202 std::vector<std::vector<double>> &vNew,
203 std::vector<std::vector<double>> &v2New,
204 double &tMin,
205 double &tMax) {
206 bool shortener = false;
207
208 dataType birthParent = 0.0;
209 dataType deathParent = 1.0;
210
211 auto interpolant = [&](int i, double t) {
212 return -vNew[node][i] + t * (vNew[node][i] + v2New[node][i]);
213 };
214
215 // Compute V1 extremity
216 auto newBirthV1 = baryBirth + interpolant(0, tMin);
217 auto newDeathV1 = baryDeath + interpolant(1, tMin);
218 auto parentBirthV1 = 0.0;
219 auto parentDeathV1 = 1.0;
220
221 // Compute V2 extremity
222 auto newBirthV2 = baryBirth + interpolant(0, tMax);
223 auto newDeathV2 = baryDeath + interpolant(1, tMax);
224 auto parentBirthV2 = 0.0;
225 auto parentDeathV2 = 1.0;
226
227 //
228 bool extremOutPathIn
229 = ((newDeathV1 > parentDeathV1 and not(newBirthV1 < parentBirthV1)
230 and not(newDeathV2 > parentDeathV2) and newBirthV2 < parentBirthV2)
231 or (not(newDeathV1 > parentDeathV1) and newBirthV1 < parentBirthV1
232 and newDeathV2 > parentDeathV2
233 and not(newBirthV2 < parentBirthV2)));
234
235 // --- Adjust V
236 // Both extremities outside and path outside
237 if(((newDeathV1 > parentDeathV1 and newDeathV2 > parentDeathV2)
238 or (newBirthV1 < parentBirthV1 and newBirthV2 < parentBirthV2))
239 and not(extremOutPathIn)) {
240 shortener = true;
241
242 // Shorten V1
243 double coef1 = 1.0;
244 if(newDeathV1 > parentDeathV1 and newBirthV1 < parentBirthV1)
245 if(newBirthV1 < newDeathV1)
246 coef1 = (parentBirthV1 - baryBirth) / interpolant(0, tMin);
247 else
248 coef1 = (parentDeathV1 - baryDeath) / interpolant(1, tMin);
249 else if(newBirthV1 < parentBirthV1)
250 coef1 = (parentBirthV1 - baryBirth) / interpolant(0, tMin);
251 else if(newDeathV1 > parentDeathV1)
252 coef1 = (parentDeathV1 - baryDeath) / interpolant(1, tMin);
253 vNew[node][0] *= coef1;
254 vNew[node][1] *= coef1;
255
256 // Shorten V2
257 double coef2 = 1.0;
258 if(newDeathV2 > parentDeathV2 and newBirthV2 < parentBirthV2)
259 if(newBirthV2 < newDeathV2)
260 coef2 = (parentBirthV2 - baryBirth) / interpolant(0, tMax);
261 else
262 coef2 = (parentDeathV2 - baryDeath) / interpolant(1, tMax);
263 else if(newBirthV2 < parentBirthV2)
264 coef2 = (parentBirthV2 - baryBirth) / interpolant(0, tMax);
265 else if(newDeathV2 > parentDeathV2)
266 coef2 = (parentDeathV2 - baryDeath) / interpolant(1, tMax);
267 v2New[node][0] *= coef2;
268 v2New[node][1] *= coef2;
269 }
270
271 // --- Adjust T
272 // One extremity outside or both extremities outside but path inside
273 if(((newDeathV1 > parentDeathV1 or newBirthV1 < parentBirthV1)
274 != (newDeathV2 > parentDeathV2 or newBirthV2 < parentBirthV2))
275 or extremOutPathIn) {
276 double m = getGeodesicVectorMiddle(vNew[node], v2New[node]);
277 if(newDeathV1 > parentDeathV1 or newDeathV2 > parentDeathV2) {
278 double newT = (deathParent - baryDeath + vNew[node][1])
279 / (vNew[node][1] + v2New[node][1]);
280 if(newDeathV1 > parentDeathV1)
281 updateT(newT, m, tMin, tMax, true);
282 if(newDeathV2 > parentDeathV2)
283 updateT(newT, m, tMin, tMax, false);
284 }
285 if(newBirthV1 < parentBirthV1 or newBirthV2 < parentBirthV2) {
286 double newT = (birthParent - baryBirth + vNew[node][0])
287 / (vNew[node][0] + v2New[node][0]);
288 if(newBirthV1 < parentBirthV1)
289 updateT(newT, m, tMin, tMax, true);
290 if(newBirthV2 < parentBirthV2)
291 updateT(newT, m, tMin, tMax, false);
292 }
293 }
294
295 return shortener;
296 }
297
298 template <class dataType>
300 dataType death,
301 std::vector<std::vector<double>> &vNew,
302 std::vector<std::vector<double>> &v2New,
303 int i,
304 double t,
305 double tMin,
306 double tMax,
307 const std::string &msg,
308 std::stringstream &ssT) {
309 double tNew = (t * (tMax - tMin) + tMin);
310 std::streamsize sSize = std::cout.precision();
311 ssT << std::setprecision(12) << std::endl << msg << std::endl;
312 ssT << "interBirth : "
313 << birth - vNew[i][0] + tNew * (vNew[i][0] + v2New[i][0]) << " _ "
314 << "interDeath : "
315 << death - vNew[i][1] + tNew * (vNew[i][1] + v2New[i][1])
316 << std::endl;
317 ssT << "v : " << vNew[i][0] << " _ " << vNew[i][1] << std::endl;
318 ssT << "v2 : " << v2New[i][0] << " _ " << v2New[i][1] << std::endl;
319 ssT << "ts : " << tMin << " _ " << tMax << std::endl;
320 ssT << "t : " << t << " _ tNew : " << tNew << std::endl;
321 std::cout.precision(sSize);
322 }
323
324 template <class dataType>
326 std::vector<std::vector<double>> &v,
327 std::vector<std::vector<double>> &v2,
328 int i,
329 double t) {
330 ftm::FTMTree_MT *baryTree = &(barycenter.tree);
331 double tMin = 0.0, tMax = 1.0;
332
333 // - Get node information
334 auto birthDeath = baryTree->getBirthDeath<dataType>(i);
335 auto birth = std::get<0>(birthDeath);
336 auto death = std::get<1>(birthDeath);
337
338 // - Get parent information
339 auto parent = baryTree->getParentSafe(i);
340 auto parentBirthDeath = baryTree->getBirthDeath<dataType>(parent);
341 auto parentBirth = std::get<0>(parentBirthDeath);
342 auto parentDeath = std::get<1>(parentBirthDeath);
343
344 if(normalizedWasserstein_ and !baryTree->notNeedToNormalize(i)) {
345 birth = (birth - parentBirth) / (parentDeath - parentBirth);
346 death = (death - parentBirth) / (parentDeath - parentBirth);
347 }
348
349 // - Adjust Diagonal T
350 bool diagonalShortener
351 = adjustDiagonalT(birth, death, i, v, v2, tMin, tMax);
352
353 // - Adjust Nesting T
354 bool nestingShortener = false;
355 if(normalizedWasserstein_ and !baryTree->notNeedToNormalize(i)) {
356 nestingShortener
357 = adjustNestingT(barycenter, birth, death, i, v, v2, tMin, tMax);
358 }
359
360 // - Compute new t
361 double tNew = t * (tMax - tMin) + tMin;
362
363 if(diagonalShortener or nestingShortener)
364 printWrn("[getTNew] shortener");
365
366 return tNew;
367 }
368
369 template <class dataType>
371 std::vector<double *> &v,
372 std::vector<double *> &v2,
373 size_t vSize,
374 double t,
375 std::vector<dataType> &interpolationVectorT,
376 bool transposeVector) {
377 // --- Init
378 ftm::FTMTree_MT *baryTree = &(barycenter.tree);
379 std::vector<dataType> scalarsVector;
380 ttk::ftm::getTreeScalars<dataType>(barycenter, scalarsVector);
381 std::vector<double> interpolationVector(scalarsVector.size());
382 std::vector<std::vector<double>> vNew, v2New;
383 pointersToVectors(v, vSize, vNew);
384 pointersToVectors(v2, vSize, v2New);
385 if(transposeVector) {
386 std::vector<std::vector<double>> out;
388 vNew = out;
390 v2New = out;
391 }
392
393 // --- Tree traversal
394 int cptBelowDiagonal = 0, cptNotNesting = 0;
395 std::queue<ftm::idNode> queue;
396 queue.emplace(baryTree->getRoot());
397 int noNestingShortener = 0, noDiagonalShortener = 0;
398 while(!queue.empty()) {
399 ftm::idNode i = queue.front();
400 queue.pop();
401
402 // - Get node information
403 auto nodeBirthDeath = baryTree->getBirthDeathNode<dataType>(i);
404 auto nodeBirth = std::get<0>(nodeBirthDeath);
405 auto nodeDeath = std::get<1>(nodeBirthDeath);
406 auto birth = scalarsVector[nodeBirth];
407 auto death = scalarsVector[nodeDeath];
408
409 // - Get parent information
410 auto parent = baryTree->getParentSafe(i);
411 auto parentNodeBirthDeath
412 = baryTree->getBirthDeathNode<dataType>(parent);
413 auto parentNodeBirth = std::get<0>(parentNodeBirthDeath);
414 auto parentNodeDeath = std::get<1>(parentNodeBirthDeath);
415 auto parentBirth = scalarsVector[parentNodeBirth];
416 auto parentDeath = scalarsVector[parentNodeDeath];
417 auto interParentBirth = interpolationVector[parentNodeBirth];
418 auto interParentDeath = interpolationVector[parentNodeDeath];
419
420 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
421 std::stringstream ssT;
422 std::streamsize sSize2 = std::cout.precision();
423 ssT << std::setprecision(12);
424 ssT << "info" << std::endl;
425 ssT << "birth death : " << birth << " ___ " << death << std::endl;
426 ssT << "par. birth death : " << parentBirth << " ___ " << parentDeath
427 << std::endl;
428 ssT << "par. inter birth death : " << interParentBirth << " ___ "
429 << interParentDeath << std::endl;
430
431 if(normalizedWasserstein_ and !baryTree->notNeedToNormalize(i)) {
432 birth = (birth - parentBirth) / (parentDeath - parentBirth);
433 death = (death - parentBirth) / (parentDeath - parentBirth);
434 ssT << "norm. birth death : " << birth << " ___ " << death
435 << std::endl;
436 }
437 std::cout.precision(sSize2);
438 // --------------------------------------------------------------------
439
440 // - Adjust tMin and tMax
441 double tMin = 0.0, tMax = 1.0;
442
443 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
444 getInterpolationVectorDebugMsg<dataType>(birth, death, vNew, v2New, i,
445 t, tMin, tMax,
446 "before adjustDiagonalT", ssT);
447 // --------------------------------------------------------------------
448
449 // - Adjust Diagonal T
450 bool diagonalShortener
451 = adjustDiagonalT(birth, death, i, vNew, v2New, tMin, tMax);
452 noDiagonalShortener += diagonalShortener;
453
454 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
455 getInterpolationVectorDebugMsg<dataType>(
456 birth, death, vNew, v2New, i, t, tMin, tMax,
457 "after adjustDiagonalT/before adjustNestingT", ssT);
458 // --------------------------------------------------------------------
459
460 // - Adjust Nesting T
461 if(normalizedWasserstein_ and !baryTree->notNeedToNormalize(i)) {
462 bool nestingShortener = adjustNestingT(
463 barycenter, birth, death, i, vNew, v2New, tMin, tMax);
464 noNestingShortener += nestingShortener;
465
466 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
467 getInterpolationVectorDebugMsg<dataType>(birth, death, vNew, v2New, i,
468 t, tMin, tMax,
469 "after adjustNestingT", ssT);
470 // ------------------------------------------------------------------
471 }
472
473 // - Compute interpolation
474 double tNew = t * (tMax - tMin) + tMin;
475 double interBirth
476 = birth - vNew[i][0] + tNew * (vNew[i][0] + v2New[i][0]);
477 double interDeath
478 = death - vNew[i][1] + tNew * (vNew[i][1] + v2New[i][1]);
479 if(normalizedWasserstein_ and !baryTree->notNeedToNormalize(i)) {
480 double coef = (interParentDeath - interParentBirth);
481 interBirth *= coef;
482 interBirth += interParentBirth;
483 interDeath *= coef;
484 interDeath += interParentBirth;
485 }
486 interpolationVector[nodeBirth] = interBirth;
487 interpolationVector[nodeDeath] = interDeath;
488
489 // - Verify NaN
490 if(std::isnan(interpolationVector[nodeBirth])
491 or std::isnan(interpolationVector[nodeDeath])) {
492 printMsg(ssT.str());
493 printErr("NOT A NUMBER");
494 }
495 // - Verify points below diagonal
496 double eps = POINTS_BELOW_DIAG_TOLERANCE;
497 if(interpolationVector[nodeBirth] > interpolationVector[nodeDeath]
498 and interpolationVector[nodeBirth] - interpolationVector[nodeDeath]
499 < eps) {
500 interpolationVector[nodeBirth] = interpolationVector[nodeDeath];
501 }
502 if(interpolationVector[nodeBirth] > interpolationVector[nodeDeath]) {
503 ++cptBelowDiagonal;
504 printMsg(ssT.str());
505 std::streamsize sSize = std::cout.precision();
506 std::stringstream ss;
507 ss << std::setprecision(12) << interpolationVector[nodeBirth] << " > "
508 << interpolationVector[nodeDeath];
509 printErr(ss.str());
510 std::cout.precision(sSize);
511 }
512 // - Verify nesting condition
513 bool verifyNesting = !baryTree->notNeedToNormalize(i);
514 if(normalizedWasserstein_ and verifyNesting) {
515 if(interpolationVector[nodeBirth] < interParentBirth
516 and (interParentBirth - interpolationVector[nodeBirth]) < eps)
517 interpolationVector[nodeBirth] = interParentBirth;
518 if(interpolationVector[nodeBirth] > interParentDeath
519 and (interpolationVector[nodeBirth] - interParentDeath) < eps) {
520 interpolationVector[nodeBirth] = interParentDeath;
521 }
522 if(interpolationVector[nodeDeath] < interParentBirth
523 and (interParentBirth - interpolationVector[nodeDeath]) < eps) {
524 interpolationVector[nodeDeath] = interParentBirth;
525 }
526 if(interpolationVector[nodeDeath] > interParentDeath
527 and (interpolationVector[nodeDeath] - interParentDeath) < eps)
528 interpolationVector[nodeDeath] = interParentDeath;
529 }
530 if(normalizedWasserstein_ and verifyNesting
531 and (interpolationVector[nodeBirth] < interParentBirth
532 or interpolationVector[nodeBirth] > interParentDeath
533 or interpolationVector[nodeDeath] < interParentBirth
534 or interpolationVector[nodeDeath] > interParentDeath)) {
535 ++cptNotNesting;
536 printMsg(ssT.str());
537 std::streamsize sSize = std::cout.precision();
538 std::stringstream ss;
539 ss << std::setprecision(12) << interpolationVector[nodeBirth] << " _ "
540 << interpolationVector[nodeDeath] << " --- " << interParentBirth
541 << " _ " << interParentDeath;
542 printErr(ss.str());
543 std::cout.precision(sSize);
544 }
545
546 // - Push children to the queue
547 std::vector<ftm::idNode> children;
548 baryTree->getChildren(i, children);
549 for(auto child : children)
550 queue.emplace(child);
551 }
552 if(noNestingShortener != 0)
553 printWrn("[getInterpolationVector] "
554 + std::to_string(noNestingShortener)
555 + " nesting vector shortener.");
556 if(noDiagonalShortener != 0)
557 printMsg("[getInterpolationVector] "
558 + std::to_string(noDiagonalShortener)
559 + " diagonal vector shortener.",
561 if(cptBelowDiagonal != 0)
562 printErr("[getInterpolationVector] point below diagonal.");
563 if(cptNotNesting != 0)
564 printErr("[getInterpolationVector] nesting condition not ok.");
565
566 interpolationVectorT.resize(scalarsVector.size());
567 for(unsigned int i = 0; i < interpolationVectorT.size(); ++i)
568 interpolationVectorT[i] = interpolationVector[i];
569 }
570
571 template <class dataType>
573 std::vector<double *> &v,
574 std::vector<double *> &v2,
575 size_t vSize,
576 double t,
577 ftm::MergeTree<dataType> &interpolated,
578 bool transposeVector = true) {
579 ftm::FTMTree_MT *barycenterTree = &(barycenter.tree);
580
581 // Compute new scalars vector
582 std::vector<dataType> newScalarsVector;
583 getInterpolationVector<dataType>(
584 barycenter, v, v2, vSize, t, newScalarsVector, transposeVector);
585
586 // Create new tree
587 interpolated
588 = ttk::ftm::createEmptyMergeTree<dataType>(newScalarsVector.size());
589 ttk::ftm::setTreeScalars<dataType>(interpolated, newScalarsVector);
590 ftm::FTMTree_MT *treeNew = &(interpolated.tree);
591
592 // Copy the barycenter tree structure
593 treeNew->copyMergeTreeStructure(barycenterTree);
594
595 // Delete diagonal and almost-diagonal points
596 persistenceThresholding<dataType>(treeNew, 0.001); // 0.001 %
597 }
598
599 template <class dataType>
601 std::vector<std::vector<double>> &v,
602 std::vector<std::vector<double>> &v2,
603 double t,
604 ftm::MergeTree<dataType> &interpolated) {
605 std::vector<double *> pV, pV2;
606 vectorsToPointers(v, pV);
607 vectorsToPointers(v2, pV2);
608 return getInterpolation<dataType>(
609 barycenter, pV, pV2, v[0].size(), t, interpolated, false);
610 }
611
612 template <class dataType>
614 std::vector<std::vector<double *>> &vS,
615 std::vector<std::vector<double *>> &v2s,
616 size_t vSize,
617 std::vector<double> &ts,
618 ftm::MergeTree<dataType> &interpolated,
619 bool transposeVector = true) {
620 getInterpolation<dataType>(
621 barycenter, vS[0], v2s[0], vSize, ts[0], interpolated, transposeVector);
622 for(unsigned int i = 1; i < vS.size(); ++i) {
623 ftm::MergeTree<dataType> newInterpolated;
624 getInterpolation<dataType>(interpolated, vS[i], v2s[i], vSize, ts[i],
625 newInterpolated, transposeVector);
626 interpolated = newInterpolated;
627 }
628 }
629
630 template <class dataType>
631 void
633 std::vector<std::vector<std::vector<double>>> &vS,
634 std::vector<std::vector<std::vector<double>>> &v2s,
635 std::vector<double> &ts,
636 ftm::MergeTree<dataType> &interpolated) {
637 std::vector<std::vector<double *>> pVS, pV2s;
639 vectorOfVectorsToPointers(v2s, pV2s);
641 barycenter, pVS, pV2s, vS[0][0].size(), ts, interpolated, false);
642 }
643
644 template <class dataType>
645 void
647 std::vector<std::vector<std::vector<double>>> &vS,
648 std::vector<std::vector<std::vector<double>>> &v2s,
649 std::vector<std::vector<double>> &v,
650 std::vector<std::vector<double>> &v2,
651 std::vector<double> &ts,
652 double t,
653 ftm::MergeTree<dataType> &interpolated) {
654 std::vector<std::vector<std::vector<double>>> vSTemp = vS, v2sTemp = v2s;
655 vSTemp.push_back(v);
656 v2sTemp.push_back(v2);
657 std::vector<double> tsTemp = ts;
658 if(tsTemp.size() != 0)
659 tsTemp[vSTemp.size() - 1] = t;
660 else
661 tsTemp.push_back(t);
662 getMultiInterpolation(barycenter, vSTemp, v2sTemp, tsTemp, interpolated);
663 }
664
665 // ----------------------------------------------------------------------------
666 // Vector Utils
667 // ----------------------------------------------------------------------------
668 void callGramSchmidt(std::vector<std::vector<double>> &vS,
669 std::vector<double> &v,
670 std::vector<double> &newV);
671
672 void vectorToPointer(std::vector<double> &vec, double *&pVec);
673
674 void vectorsToPointers(std::vector<std::vector<double>> &vec,
675 std::vector<double *> &pVec);
676
678 std::vector<std::vector<std::vector<double>>> &vS,
679 std::vector<std::vector<double *>> &pVS);
680
681 void pointerToVector(double *pVec, size_t size, std::vector<double> &vec);
682
683 void pointersToVectors(std::vector<double *> &pVec,
684 std::vector<size_t> sizes,
685 std::vector<std::vector<double>> &vec);
686
687 void pointersToVectors(std::vector<double *> &pVec,
688 size_t size,
689 std::vector<std::vector<double>> &vec);
690 }; // MergeTreePrincipalGeodesicsBase class
691
692} // namespace ttk
#define ttkNotUsed(x)
Mark function/method parameters that are not used in the function body at all.
Definition: BaseClass.h:47
#define POINTS_BELOW_DIAG_TOLERANCE
Minimalist debugging class.
Definition: Debug.h:88
int printWrn(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition: Debug.h:159
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
double getGeodesicVectorMiddle(std::vector< double > &v, std::vector< double > &v2)
dataType computeReconstructionError(ftm::MergeTree< dataType > &barycenter, std::vector< ftm::MergeTree< dataType > > &inputTrees, std::vector< std::vector< std::vector< double > > > &vS, std::vector< std::vector< std::vector< double > > > &v2s, std::vector< std::vector< double > > &allTreesTs)
void getInterpolationVectorDebugMsg(dataType birth, dataType death, std::vector< std::vector< double > > &vNew, std::vector< std::vector< double > > &v2New, int i, double t, double tMin, double tMax, const std::string &msg, std::stringstream &ssT)
std::vector< std::vector< double > > allTreesTs_
std::vector< std::vector< std::vector< double > > > v2s_
void vectorToPointer(std::vector< double > &vec, double *&pVec)
void pointersToVectors(std::vector< double * > &pVec, std::vector< size_t > sizes, std::vector< std::vector< double > > &vec)
double getTNew(ftm::MergeTree< dataType > &barycenter, std::vector< std::vector< double > > &v, std::vector< std::vector< double > > &v2, int i, double t)
std::vector< std::vector< std::vector< double > > > trees2V2s_
void getMultiInterpolation(ftm::MergeTree< dataType > &barycenter, std::vector< std::vector< double * > > &vS, std::vector< std::vector< double * > > &v2s, size_t vSize, std::vector< double > &ts, ftm::MergeTree< dataType > &interpolated, bool transposeVector=true)
std::vector< std::vector< std::vector< double > > > vS_
void getInterpolation(ftm::MergeTree< dataType > &barycenter, std::vector< double * > &v, std::vector< double * > &v2, size_t vSize, double t, ftm::MergeTree< dataType > &interpolated, bool transposeVector=true)
void vectorsToPointers(std::vector< std::vector< double > > &vec, std::vector< double * > &pVec)
bool adjustNestingT(ftm::MergeTree< dataType > &ttkNotUsed(barycenter), dataType baryBirth, dataType baryDeath, ftm::idNode node, std::vector< std::vector< double > > &vNew, std::vector< std::vector< double > > &v2New, double &tMin, double &tMax)
std::vector< std::vector< double > > allTs_
std::vector< std::vector< double > > branchesCorrelationMatrix_
void getMultiInterpolation(ftm::MergeTree< dataType > &barycenter, std::vector< std::vector< std::vector< double > > > &vS, std::vector< std::vector< std::vector< double > > > &v2s, std::vector< double > &ts, ftm::MergeTree< dataType > &interpolated)
std::vector< std::vector< double > > persCorrelationMatrix_
std::vector< std::vector< std::tuple< ftm::idNode, ftm::idNode, double > > > baryMatchings_
std::vector< std::vector< std::tuple< ftm::idNode, ftm::idNode, double > > > baryMatchings2_
bool adjustDiagonalT(dataType baryBirth, dataType baryDeath, ftm::idNode node, std::vector< std::vector< double > > &vNew, std::vector< std::vector< double > > &v2New, double &tMin, double &tMax)
void getMultiInterpolation(ftm::MergeTree< dataType > &barycenter, std::vector< std::vector< std::vector< double > > > &vS, std::vector< std::vector< std::vector< double > > > &v2s, std::vector< std::vector< double > > &v, std::vector< std::vector< double > > &v2, std::vector< double > &ts, double t, ftm::MergeTree< dataType > &interpolated)
void pointerToVector(double *pVec, size_t size, std::vector< double > &vec)
void getInterpolation(ftm::MergeTree< dataType > &barycenter, std::vector< std::vector< double > > &v, std::vector< std::vector< double > > &v2, double t, ftm::MergeTree< dataType > &interpolated)
void getInterpolationVector(ftm::MergeTree< dataType > &barycenter, std::vector< double * > &v, std::vector< double * > &v2, size_t vSize, double t, std::vector< dataType > &interpolationVectorT, bool transposeVector)
std::vector< std::vector< std::vector< double > > > trees2Vs_
std::vector< std::vector< double > > allScaledTs_
void vectorOfVectorsToPointers(std::vector< std::vector< std::vector< double > > > &vS, std::vector< std::vector< double * > > &pVS)
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 updateT(double newT, double m, double &tMin, double &tMax, bool updateTMin)
dataType computeReconstructionError(ftm::MergeTree< dataType > &barycenter, std::vector< ftm::MergeTree< dataType > > &inputTrees, std::vector< std::vector< std::vector< double > > > &vS, std::vector< std::vector< std::vector< double > > > &v2s, std::vector< std::vector< double > > &allTreesTs, std::vector< double > &reconstructionErrors)
void callGramSchmidt(std::vector< std::vector< double > > &vS, std::vector< double > &v, std::vector< double > &newV)
bool notNeedToNormalize(idNode nodeId)
std::tuple< dataType, dataType > getBirthDeath(idNode nodeId)
std::tuple< ftm::idNode, ftm::idNode > getBirthDeathNode(idNode nodeId)
idNode getParentSafe(idNode nodeId)
void getChildren(idNode nodeId, std::vector< idNode > &res)
void copyMergeTreeStructure(FTMTree_MT *tree)
T1 pow(const T1 val, const T2 n)
Definition: Geometry.h:411
void addVectorsProjection(const std::vector< T > &a, const std::vector< T > &b, std::vector< T > &a_out, std::vector< T > &b_out)
Definition: Geometry.cpp:618
void transposeMatrix(const std::vector< std::vector< T > > &a, std::vector< std::vector< T > > &out)
Definition: Geometry.cpp:746
unsigned int idNode
Node index in vect_nodes_.
Definition: FTMDataTypes.h:39
The Topology ToolKit.
ftm::FTMTree_MT tree
Definition: FTMTree_MT.h:901