36 :
virtual public Debug,
63 template <
class dataType>
66 bool const useMinMax =
true;
67 bool const cleanTree =
false;
69 std::vector<int> nodeCorr;
70 preprocessingPipeline<dataType>(barycenter, 0.0, 100.0, 100.0,
72 cleanTree, pt, nodeCorr,
false);
76 template <
class dataType>
80#ifdef TTK_ENABLE_OPENMP
81#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
83 for(
unsigned int i = 0; i < inputTrees.size(); ++i)
84 postprocessingPipeline<dataType>(&(inputTrees[i].tree));
87 template <
class dataType>
89 std::vector<std::vector<double *>> &vS,
90 std::vector<std::vector<double *>> &v2s,
92 std::array<double, 2> &middle) {
94 std::vector<double> alpha(2, 0.0);
95 for(
unsigned int i = 0; i < 2; ++i) {
97 for(
unsigned int j = 0; j < vSize; ++j) {
100 for(
unsigned int k = 0; k < 2; ++k) {
103 alpha[i] += (vS[i][k][j] / v2s[i][k][j]);
107 alpha[i] /= cptDivide;
110 for(
unsigned int i = 0; i < 2; ++i)
111 middle[i] = alpha[i] / (1 + alpha[i]);
114 template <
class dataType>
118 preprocessBarycenter<dataType>(barycenters[0]);
119 if(barycenters.size() > 1)
120 preprocessBarycenter<dataType>(barycenters[1]);
124#ifdef TTK_ENABLE_OPENMP
125#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
127 for(
unsigned int i = 0; i <
pVS_.size(); ++i) {
129 getInterpolation<dataType>(
131 getInterpolation<dataType>(
137 if(barycenters.size() > 1) {
139 getInterpolation<dataType>(barycenters[1],
pTrees2Vs_[i],
142 getInterpolation<dataType>(barycenters[1],
pTrees2Vs_[i],
155 postprocessingPipeline<dataType>(&(barycenters[0].tree));
156 if(barycenters.size() > 1)
157 postprocessingPipeline<dataType>(&(barycenters[1].tree));
163 template <
class dataType>
168 std::vector<double> &reconstructionErrors,
169 std::vector<std::vector<std::tuple<ftm::idNode, ftm::idNode, double>>>
171 std::vector<std::vector<std::tuple<ftm::idNode, ftm::idNode, double>>>
173 bool isSecondInput =
false) {
179 preprocessBarycenter<dataType>(barycenter);
180 if(inputTrees.size() != 0)
181 preprocessingTrees<dataType>(inputTrees);
186 recBaryMatchings.resize(reconstructedTrees.size());
187#ifdef TTK_ENABLE_OPENMP
188#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
190 for(
unsigned int i = 0; i < reconstructedTrees.size(); ++i) {
191 getMultiInterpolation<dataType>(barycenter, vSToUse, v2sToUse,
193 reconstructedTrees[i]);
196 computeOneDistance<dataType>(reconstructedTrees[i], barycenter,
197 recBaryMatchings[i], distance,
true);
202 if(inputTrees.size() != 0
205 barycenter, inputTrees, vSToUse, v2sToUse, vSizeToUse,
allTreesTs_,
206 reconstructionErrors, recInputMatchings);
208 std::stringstream ss;
209 ss <<
"Reconstruction Error = " << reconstructionError;
215 postprocessingPipeline<dataType>(&(barycenter.
tree));
216#ifdef TTK_ENABLE_OPENMP
217#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
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_)
224 for(
unsigned int i = 0; i < inputTrees.size(); ++i)
225 postprocessingPipeline<dataType>(&(inputTrees[i].tree));
228 for(
unsigned int i = 0; i < inputTrees.size(); ++i)
229 convertBranchDecompositionMatching<dataType>(
230 &(reconstructedTrees[i].tree), &(inputTrees[i].tree),
231 recInputMatchings[i]);
234 for(
unsigned int i = 0; i < reconstructedTrees.size(); ++i)
235 convertBranchDecompositionMatching<dataType>(
236 &(reconstructedTrees[i].tree), &(barycenter.
tree),
237 recBaryMatchings[i]);
240 template <
class dataType>
244 bool isSecondInput =
false) {
250 preprocessBarycenter<dataType>(barycenter);
252 std::array<double, 2> middle;
253 getGeodesicsMiddle<dataType>(
254 barycenter, vSToUse, v2sToUse, vSizeToUse, middle);
257 geodesicsTrees.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_)
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]);
272 int const i2 = (i + 1) % 2;
277 postprocessingPipeline<dataType>(&(barycenter.
tree));
278#ifdef TTK_ENABLE_OPENMP
279#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
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));
286 template <
class dataType>
290 bool isSecondInput =
false) {
296 preprocessBarycenter<dataType>(barycenter);
299 unsigned int noSample =
k_ * 2;
300 geodesicsEllipses.resize(noSample);
303 std::vector<std::vector<double *>> vS(2), v2s(2);
306 v2s[0] = v2sToUse[0];
307 v2s[1] = v2sToUse[1];
310 std::array<double, 2> middle;
311 getGeodesicsMiddle<dataType>(barycenter, vS, v2s, vSizeToUse, middle);
314#ifdef TTK_ENABLE_OPENMP
315#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
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);
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];
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]);
345 postprocessingPipeline<dataType>(&(barycenter.
tree));
346#ifdef TTK_ENABLE_OPENMP
347#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
349 for(
unsigned int i = 0; i < geodesicsEllipses.size(); ++i)
350 postprocessingPipeline<dataType>(&(geodesicsEllipses[i].tree));
354 return k_ * 4 * rectangleMultiplier;
357 template <
class dataType>
361 unsigned int rectangleMultiplier = 1,
362 bool isSecondInput =
false) {
368 preprocessBarycenter<dataType>(barycenter);
372 geodesicsRectangle.resize(noSample);
375 std::vector<std::vector<double *>> vS(2), v2s(2);
378 v2s[0] = v2sToUse[0];
379 v2s[1] = v2sToUse[1];
382#ifdef TTK_ENABLE_OPENMP
383#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
385 for(
unsigned int i = 0; i < noSample; ++i) {
387 int const quadrant = i / (noSample / 4);
389 double x = 0.0, y = 0.0;
390 double const offset = (i % (noSample / 4)) / ((noSample / 4.0) - 1.0);
411 std::vector<double> ts{x, y};
412 getMultiInterpolation<dataType>(
413 barycenter, vS, v2s, vSizeToUse, ts, geodesicsRectangle[i]);
418 postprocessingPipeline<dataType>(&(barycenter.
tree));
419#ifdef TTK_ENABLE_OPENMP
420#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
422 for(
unsigned int i = 0; i < geodesicsRectangle.size(); ++i)
423 postprocessingPipeline<dataType>(&(geodesicsRectangle[i].tree));
426 template <
class dataType>
430 bool isSecondInput =
false) {
436 preprocessBarycenter<dataType>(barycenter);
439 unsigned int noSample =
k_ *
k_;
440 geodesicsSurface.resize(noSample);
441 tSurface_.resize(geodesicsSurface.size());
445 std::vector<std::vector<double *>> vS(2), v2s(2);
448 v2s[0] = v2sToUse[0];
449 v2s[1] = v2sToUse[1];
452#ifdef TTK_ENABLE_OPENMP
453#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
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;
461 if(i == 0 or j == 0 or i ==
k_ - 1 or j ==
k_ - 1) {
469 boundaryID =
k_ * 2 + (
k_ - j);
471 boundaryID =
k_ * 3 + (
k_ - i);
475 std::vector<double> ts{x, y};
476 getMultiInterpolation<dataType>(
477 barycenter, vS, v2s, vSizeToUse, ts, geodesicsSurface[index]);
483 postprocessingPipeline<dataType>(&(barycenter.
tree));
484#ifdef TTK_ENABLE_OPENMP
485#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
487 for(
unsigned int i = 0; i < geodesicsSurface.size(); ++i)
488 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)