33 :
virtual public Debug,
60 template <
class dataType>
63 bool const useMinMax =
true;
64 bool const cleanTree =
false;
66 std::vector<int> nodeCorr;
67 preprocessingPipeline<dataType>(barycenter, 0.0, 100.0, 100.0,
69 cleanTree, pt, nodeCorr,
false);
73 template <
class dataType>
77#ifdef TTK_ENABLE_OPENMP
78#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
80 for(
unsigned int i = 0; i < inputTrees.size(); ++i)
81 postprocessingPipeline<dataType>(&(inputTrees[i].tree));
84 template <
class dataType>
86 std::vector<std::vector<double *>> &vS,
87 std::vector<std::vector<double *>> &v2s,
89 std::array<double, 2> &middle) {
91 std::vector<double> alpha(2, 0.0);
92 for(
unsigned int i = 0; i < 2; ++i) {
94 for(
unsigned int j = 0; j < vSize; ++j) {
97 for(
unsigned int k = 0; k < 2; ++k) {
100 alpha[i] += (vS[i][k][j] / v2s[i][k][j]);
104 alpha[i] /= cptDivide;
107 for(
unsigned int i = 0; i < 2; ++i)
108 middle[i] = alpha[i] / (1 + alpha[i]);
111 template <
class dataType>
115 preprocessBarycenter<dataType>(barycenters[0]);
116 if(barycenters.size() > 1)
117 preprocessBarycenter<dataType>(barycenters[1]);
121#ifdef TTK_ENABLE_OPENMP
122#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
124 for(
unsigned int i = 0; i <
pVS_.size(); ++i) {
126 getInterpolation<dataType>(
128 getInterpolation<dataType>(
134 if(barycenters.size() > 1) {
136 getInterpolation<dataType>(barycenters[1],
pTrees2Vs_[i],
139 getInterpolation<dataType>(barycenters[1],
pTrees2Vs_[i],
152 postprocessingPipeline<dataType>(&(barycenters[0].tree));
153 if(barycenters.size() > 1)
154 postprocessingPipeline<dataType>(&(barycenters[1].tree));
160 template <
class dataType>
165 std::vector<double> &reconstructionErrors,
166 std::vector<std::vector<std::tuple<ftm::idNode, ftm::idNode, double>>>
168 std::vector<std::vector<std::tuple<ftm::idNode, ftm::idNode, double>>>
170 bool isSecondInput =
false) {
176 preprocessBarycenter<dataType>(barycenter);
177 if(inputTrees.size() != 0)
178 preprocessingTrees<dataType>(inputTrees);
183 recBaryMatchings.resize(reconstructedTrees.size());
184#ifdef TTK_ENABLE_OPENMP
185#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
187 for(
unsigned int i = 0; i < reconstructedTrees.size(); ++i) {
188 getMultiInterpolation<dataType>(barycenter, vSToUse, v2sToUse,
190 reconstructedTrees[i]);
193 computeOneDistance<dataType>(reconstructedTrees[i], barycenter,
194 recBaryMatchings[i], distance,
true);
199 if(inputTrees.size() != 0
202 barycenter, inputTrees, vSToUse, v2sToUse, vSizeToUse,
allTreesTs_,
203 reconstructionErrors, recInputMatchings);
205 std::stringstream ss;
206 ss <<
"Reconstruction Error = " << reconstructionError;
212 postprocessingPipeline<dataType>(&(barycenter.
tree));
213#ifdef TTK_ENABLE_OPENMP
214#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
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_)
221 for(
unsigned int i = 0; i < inputTrees.size(); ++i)
222 postprocessingPipeline<dataType>(&(inputTrees[i].tree));
225 for(
unsigned int i = 0; i < inputTrees.size(); ++i)
226 convertBranchDecompositionMatching<dataType>(
227 &(reconstructedTrees[i].tree), &(inputTrees[i].tree),
228 recInputMatchings[i]);
231 for(
unsigned int i = 0; i < reconstructedTrees.size(); ++i)
232 convertBranchDecompositionMatching<dataType>(
233 &(reconstructedTrees[i].tree), &(barycenter.
tree),
234 recBaryMatchings[i]);
237 template <
class dataType>
241 bool isSecondInput =
false) {
247 preprocessBarycenter<dataType>(barycenter);
249 std::array<double, 2> middle;
250 getGeodesicsMiddle<dataType>(
251 barycenter, vSToUse, v2sToUse, vSizeToUse, middle);
254 geodesicsTrees.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_)
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]);
269 int const i2 = (i + 1) % 2;
274 postprocessingPipeline<dataType>(&(barycenter.
tree));
275#ifdef TTK_ENABLE_OPENMP
276#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
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));
283 template <
class dataType>
287 bool isSecondInput =
false) {
293 preprocessBarycenter<dataType>(barycenter);
296 unsigned int noSample =
k_ * 2;
297 geodesicsEllipses.resize(noSample);
300 std::vector<std::vector<double *>> vS(2), v2s(2);
303 v2s[0] = v2sToUse[0];
304 v2s[1] = v2sToUse[1];
307 std::array<double, 2> middle;
308 getGeodesicsMiddle<dataType>(barycenter, vS, v2s, vSizeToUse, middle);
311#ifdef TTK_ENABLE_OPENMP
312#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
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);
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];
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]);
342 postprocessingPipeline<dataType>(&(barycenter.
tree));
343#ifdef TTK_ENABLE_OPENMP
344#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
346 for(
unsigned int i = 0; i < geodesicsEllipses.size(); ++i)
347 postprocessingPipeline<dataType>(&(geodesicsEllipses[i].tree));
351 return k_ * 4 * rectangleMultiplier;
354 template <
class dataType>
358 unsigned int rectangleMultiplier = 1,
359 bool isSecondInput =
false) {
365 preprocessBarycenter<dataType>(barycenter);
369 geodesicsRectangle.resize(noSample);
372 std::vector<std::vector<double *>> vS(2), v2s(2);
375 v2s[0] = v2sToUse[0];
376 v2s[1] = v2sToUse[1];
379#ifdef TTK_ENABLE_OPENMP
380#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
382 for(
unsigned int i = 0; i < noSample; ++i) {
384 int const quadrant = i / (noSample / 4);
386 double x = 0.0, y = 0.0;
387 double const offset = (i % (noSample / 4)) / ((noSample / 4.0) - 1.0);
408 std::vector<double> ts{x, y};
409 getMultiInterpolation<dataType>(
410 barycenter, vS, v2s, vSizeToUse, ts, geodesicsRectangle[i]);
415 postprocessingPipeline<dataType>(&(barycenter.
tree));
416#ifdef TTK_ENABLE_OPENMP
417#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
419 for(
unsigned int i = 0; i < geodesicsRectangle.size(); ++i)
420 postprocessingPipeline<dataType>(&(geodesicsRectangle[i].tree));
423 template <
class dataType>
427 bool isSecondInput =
false) {
433 preprocessBarycenter<dataType>(barycenter);
436 unsigned int noSample =
k_ *
k_;
437 geodesicsSurface.resize(noSample);
438 tSurface_.resize(geodesicsSurface.size());
442 std::vector<std::vector<double *>> vS(2), v2s(2);
445 v2s[0] = v2sToUse[0];
446 v2s[1] = v2sToUse[1];
449#ifdef TTK_ENABLE_OPENMP
450#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
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;
458 if(i == 0 or j == 0 or i ==
k_ - 1 or j ==
k_ - 1) {
466 boundaryID =
k_ * 2 + (
k_ - j);
468 boundaryID =
k_ * 3 + (
k_ - i);
472 std::vector<double> ts{x, y};
473 getMultiInterpolation<dataType>(
474 barycenter, vS, v2s, vSizeToUse, ts, geodesicsSurface[index]);
480 postprocessingPipeline<dataType>(&(barycenter.
tree));
481#ifdef TTK_ENABLE_OPENMP
482#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
484 for(
unsigned int i = 0; i < geodesicsSurface.size(); ++i)
485 postprocessingPipeline<dataType>(&(geodesicsSurface[i].tree));
Minimalist debugging class.
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
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_
bool branchDecomposition_
bool isPersistenceDiagram_
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)
std::vector< double * > pPersCorrelationMatrix_
std::vector< std::vector< double * > > pTrees2Vs_
std::vector< bool > surfaceIsBoundary_
std::vector< std::vector< double * > > pTrees2V2s_
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)
std::vector< double > geodesicsDistances_
unsigned int getNumberOfRectangles(unsigned int rectangleMultiplier=1)
std::vector< int > surfaceBoundaryID_
std::vector< std::vector< double > > tRectangle_
MergeTreePrincipalGeodesicsDecoding()
std::vector< double * > pBranchesCorrelationMatrix_
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)
std::vector< std::vector< double > > tSurface_
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)
std::vector< std::vector< double * > > pVS_
void processInputTrees(std::vector< ttk::ftm::MergeTree< dataType > > &inputTrees)
bool transferBarycenterInformation_
bool computeReconstructionError_
std::vector< std::vector< double * > > pV2s_
std::vector< std::vector< double > > tEllipses_
bool transferInputTreesInformation_
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)
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)