TTK
Loading...
Searching...
No Matches
TrackingFromPersistenceDiagrams.cpp
Go to the documentation of this file.
2
4 this->setDebugMsgPrefix("TrackingFromPersistenceDiagrams");
5}
6
9
10 // Check the consistency of the variables
11#ifndef TTK_ENABLE_KAMIKAZE
12 if(!numberOfInputs_)
13 return -1;
14 if(!inputData_)
15 return -3;
16
17 for(int i = 0; i < numberOfInputs_; i++) {
18 if(!inputData_[i])
19 return -4;
20 }
21#endif
22
23 this->printMsg("Complete", 1, t.getElapsedTime(), threadNumber_);
24
25 return 0;
26}
27
29 int i,
30 std::vector<ttk::DiagramType> &inputPersistenceDiagrams,
31 std::vector<std::vector<MatchingType>> &outputMatchings,
32 const std::string &algorithm,
33 const std::string &wasserstein,
34 double tolerance,
35 double px,
36 double py,
37 double pz,
38 double ps,
39 double pe) {
40 ttk::BottleneckDistance bottleneckDistance_;
41 // bottleneckDistance_.setWrapper(this);
42 bottleneckDistance_.setPersistencePercentThreshold(tolerance);
43 bottleneckDistance_.setPX(px);
44 bottleneckDistance_.setPY(py);
45 bottleneckDistance_.setPZ(pz);
46 bottleneckDistance_.setPS(ps);
47 bottleneckDistance_.setPE(pe);
48 bottleneckDistance_.setAlgorithm(algorithm);
49 bottleneckDistance_.setWasserstein(wasserstein);
50
51 bottleneckDistance_.execute(inputPersistenceDiagrams[i],
52 inputPersistenceDiagrams[i + 1],
53 outputMatchings[i]);
54
55 return 0;
56}
57
59 int numInputs,
60 std::vector<ttk::DiagramType> &inputPersistenceDiagrams,
61 std::vector<std::vector<MatchingType>> &outputMatchings,
62 const std::string &algorithm,
63 const std::string &wasserstein,
64 double tolerance,
65 double px,
66 double py,
67 double pz,
68 double ps,
69 double pe) {
70
71#ifdef TTK_ENABLE_OPENMP
72#pragma omp parallel for num_threads(threadNumber_)
73#endif // TTK_ENABLE_OPENMP
74 for(int i = 0; i < numInputs - 1; ++i) {
75 performSingleMatching(
76 i, inputPersistenceDiagrams, outputMatchings,
77 algorithm, // Not from paraview, from enclosing tracking plugin
78 wasserstein, tolerance, px, py, pz, ps, pe // Coefficients
79 );
80 }
81
82 // Never from PV,
83 // Alaways activate output matchings.
84 return 0;
85}
86
88 std::vector<ttk::DiagramType> &allDiagrams,
89 std::vector<std::vector<MatchingType>> &allMatchings,
90 std::vector<trackingTuple> &trackings) {
91 auto numPersistenceDiagramsInput = (int)allDiagrams.size();
92
93 for(int in = 1; in < numPersistenceDiagramsInput - 1; ++in) {
94 std::vector<MatchingType> matchings1 = allMatchings[in - 1];
95 std::vector<MatchingType> matchings2 = allMatchings[in];
96
97 auto matchingsSize1 = (int)matchings1.size();
98 auto matchingsSize2 = (int)matchings2.size();
99 int endIndex = numPersistenceDiagramsInput - 2;
100
101 for(int i = 0; i < matchingsSize1; ++i) {
102 const auto m1ai0 = std::get<0>(matchings1[i]);
103 const auto m1ai1 = std::get<1>(matchings1[i]);
104
105 for(int j = 0; j < matchingsSize2; ++j) {
106 const auto m2aj0 = std::get<0>(matchings2[j]);
107 const auto m2aj1 = std::get<1>(matchings2[j]);
108
109 if(m1ai1 != m2aj0)
110 continue;
111
112 // Detect in trackings and push.
113 bool found = false;
114 for(trackingTuple &tt : trackings) {
115 int chainStart = std::get<0>(tt);
116 int chainEnd = std::get<1>(tt);
117 std::vector<SimplexId> &chain = std::get<2>(tt);
118
119 if(chainEnd == -1) {
120 auto chainSize = (int)chain.size();
121 if(chainSize == 0) {
122 // Should not happen
123 this->printErr("Brain error.");
124 } else if(chainStart + chainSize == in
125 && chain.at((unsigned long)chainSize - 1) == m1ai0) {
126 found = true;
127 chain.push_back(m1ai1);
128 int numEnd = in == endIndex ? endIndex : -1;
129 if(in == endIndex) {
130 chain.push_back(m2aj1);
131 std::get<1>(tt) = numEnd;
132 }
133 std::get<2>(tt) = chain;
134 }
135 }
136 tt = std::make_tuple(chainStart, chainEnd, chain);
137 }
138
139 if(!found) {
140 std::vector<SimplexId> chain;
141 chain.push_back(m1ai0);
142 chain.push_back(m1ai1);
143 if(in == endIndex) {
144 chain.push_back(m2aj1);
145 }
146 int numEnd = in == endIndex ? endIndex : -1;
147 trackingTuple tt = std::make_tuple(in - 1, numEnd, chain);
148 trackings.push_back(tt);
149 }
150 // Create new.
151 }
152 }
153
154 // End non-matched chains.
155 for(trackingTuple &tt : trackings) {
156 int chainStart = std::get<0>(tt);
157 int chainEnd = std::get<1>(tt);
158 if(chainEnd == -1) {
159 std::vector<SimplexId> &chain = std::get<2>(tt);
160 auto chainSize = (int)chain.size();
161 if(chainStart + chainSize - 1 < in)
162 std::get<1>(tt) = in - 1;
163 }
164 }
165 }
166
167 // Post-processing
168 std::sort(trackings.begin(), trackings.end(),
169 [](const trackingTuple &a, const trackingTuple &b) -> bool {
170 return std::get<0>(a) < std::get<0>(b);
171 });
172
173 return 0;
174}
175
177 const std::vector<ttk::DiagramType> &allDiagrams,
178 const std::vector<trackingTuple> &trackings,
179 std::vector<std::set<int>> &trackingTupleToMerged,
180 const double postProcThresh) {
181
182 const int numPersistenceDiagramsInput = allDiagrams.size();
183
184 // Merge close connected components with threshold.
185 for(size_t k = 0; k < trackings.size(); ++k) {
186 const auto &tk = trackings[k];
187 int startK = std::get<0>(tk);
188 int endK = std::get<1>(tk);
189 if(endK < 0)
190 endK = numPersistenceDiagramsInput - 1;
191 const std::vector<SimplexId> &chainK = std::get<2>(tk);
192 const ttk::DiagramType &diagramStartK = allDiagrams[startK];
193 const ttk::DiagramType &diagramEndK = allDiagrams[endK];
194
195 const auto n1 = chainK.front();
196 const auto n2 = chainK.back();
197 const auto &tuple1 = diagramStartK[n1];
198 const auto &tuple2 = diagramEndK[n2];
199
200 double x1, y1, z1, x2, y2, z2;
201
202 const auto point1Type1 = tuple1.birth.type;
203 const auto point1Type2 = tuple1.death.type;
204 bool t11Min = point1Type1 == CriticalType::Local_minimum;
205 bool t11Max = point1Type1 == CriticalType::Local_maximum;
206 bool t12Min = point1Type2 == CriticalType::Local_minimum;
207 bool t12Max = point1Type2 == CriticalType::Local_maximum;
208 // bool bothEx1 = t11Ex && t12Ex;
209 bool t1Max = t11Max || t12Max;
210 bool t1Min = !t1Max && (t11Min || t12Min);
211
212 x1 = t1Max ? tuple1.death.coords[0] : t1Min ? tuple1.birth.coords[0] : 0;
213 y1 = t1Max ? tuple1.death.coords[1] : t1Min ? tuple1.birth.coords[1] : 0;
214 z1 = t1Max ? tuple1.death.coords[2] : t1Min ? tuple1.birth.coords[2] : 0;
215
216 const auto point2Type1 = tuple2.birth.type;
217 const auto point2Type2 = tuple2.death.type;
218 bool t21Min = point2Type1 == CriticalType::Local_minimum;
219 bool t21Max = point2Type1 == CriticalType::Local_maximum;
220 bool t22Min = point2Type2 == CriticalType::Local_minimum;
221 bool t22Max = point2Type2 == CriticalType::Local_maximum;
222 // bool bothEx2 = t21Ex && t22Ex;
223 bool t2Max = t21Max || t22Max;
224 bool t2Min = !t2Max && (t21Min || t22Min);
225
226 // if (bothEx2) {
227 x2 = t2Max ? tuple2.death.coords[0] : t2Min ? tuple2.birth.coords[0] : 0;
228 y2 = t2Max ? tuple2.death.coords[1] : t2Min ? tuple2.birth.coords[1] : 0;
229 z2 = t2Max ? tuple2.death.coords[2] : t2Min ? tuple2.birth.coords[2] : 0;
230 // }
231 // if (!bothEx1 && !bothEx2)
232 // continue;
233
234 // Saddle-saddle matching not supported.
235 if(!t1Min && !t2Min && !t1Max && !t2Max)
236 continue;
237
238 // Check every other tracking trajectory.
239 for(size_t m = k + 1; m < trackings.size(); ++m) {
240 const auto &tm = trackings[m];
241 int startM = std::get<0>(tm);
242 int endM = std::get<1>(tm);
243 const std::vector<SimplexId> &chainM = std::get<2>(tm);
244 if((endK > 0 && startM > endK) || (endM > 0 && startK > endM))
245 continue;
246
247 for(int c = 0; c < (int)chainM.size(); ++c) {
248 bool doMatch1 = startM + c == startK;
249 bool doMatch2 = startM + c == endK;
250
251 // if (startM + c != startK && startM + c != endK) continue;
252 if(!doMatch1 && !doMatch2)
253 continue;
254
256 const auto n3 = chainM[c];
257 const ttk::DiagramType &diagramM = allDiagrams[startM + c];
258 const auto &tuple3 = diagramM[n3];
259 double x3, y3, z3;
260 const auto point3Type1 = tuple3.birth.type;
261 const auto point3Type2 = tuple3.death.type;
262 bool t31Min = point3Type1 == CriticalType::Local_minimum;
263 bool t31Max = point3Type1 == CriticalType::Local_maximum;
264 bool t32Min = point3Type2 == CriticalType::Local_minimum;
265 bool t32Max = point3Type2 == CriticalType::Local_maximum;
266 // bool bothEx3 = t31Ex && t32Ex;
267 // if (!bothEx3)
268 // continue;
269 bool t3Max = t31Max || t32Max;
270 bool t3Min = !t3Max && (t31Min || t32Min);
271
272 x3 = t3Max ? tuple3.death.coords[0]
273 : t3Min ? tuple3.birth.coords[0]
274 : 0;
275 y3 = t3Max ? tuple3.death.coords[1]
276 : t3Min ? tuple3.birth.coords[1]
277 : 0;
278 z3 = t3Max ? tuple3.death.coords[2]
279 : t3Min ? tuple3.birth.coords[2]
280 : 0;
281
282 double dist = 0;
283 bool hasMatched = false;
284 const auto square = [](const double a) { return a * a; };
285 if(doMatch1 && ((t3Max && t1Max) || (t3Min && t1Min))) {
286 double dist13
287 = std::sqrt(square(x1 - x3) + square(y1 - y3) + square(z1 - z3));
288 dist = dist13;
289 if(dist13 >= postProcThresh)
290 continue;
291 hasMatched = true;
292 }
293
294 if(doMatch2 && ((t3Max && t2Max) || (t3Min && t2Min))) {
295 double dist23
296 = std::sqrt(square(x2 - x3) + square(y2 - y3) + square(z2 - z3));
297 dist = dist23;
298 if(dist23 >= postProcThresh)
299 continue;
300 hasMatched = true;
301 }
302
303 if(!hasMatched)
304 continue;
305
307 std::stringstream msg;
308 msg << "Merged " << m << " with " << k << ": d = " << dist << ".";
309 printMsg(msg.str());
310
311 // Get every other tracking trajectory.
312 std::set<int> &mergedM = trackingTupleToMerged[m];
313 // std::set<int> mergedK = trackingTupleToMerged[k];
314
315 // Push for others to merge.
316 // for (auto& i : mergedM) mergedK.insert(i);
317 // for (auto& i : mergedK) mergedM.insert(i);
318 // mergedK.insert(m);
319 mergedM.insert(k);
320 break;
321 }
322 }
323 }
324
325 return 0;
326}
void setWasserstein(const std::string &wasserstein)
void setPY(const double py)
int execute(const ttk::DiagramType &diag0, const ttk::DiagramType &diag1, std::vector< MatchingType > &matchings)
void setPX(const double px)
void setAlgorithm(const std::string &algorithm)
void setPZ(const double pz)
void setPersistencePercentThreshold(const double t)
void setPS(const double ps)
void setPE(const double pe)
void setDebugMsgPrefix(const std::string &prefix)
Definition: Debug.h:364
double getElapsedTime()
Definition: Timer.h:15
int performSingleMatching(int i, std::vector< ttk::DiagramType > &inputPersistenceDiagrams, std::vector< std::vector< MatchingType > > &outputMatchings, const std::string &algorithm, const std::string &wasserstein, double tolerance, double px, double py, double pz, double ps, double pe)
int performMatchings(int numInputs, std::vector< ttk::DiagramType > &inputPersistenceDiagrams, std::vector< std::vector< MatchingType > > &outputMatchings, const std::string &algorithm, const std::string &wasserstein, double tolerance, double px, double py, double pz, double ps, double pe)
int performTracking(std::vector< ttk::DiagramType > &allDiagrams, std::vector< std::vector< MatchingType > > &allMatchings, std::vector< trackingTuple > &trackings)
int performPostProcess(const std::vector< ttk::DiagramType > &allDiagrams, const std::vector< trackingTuple > &trackings, std::vector< std::set< int > > &trackingTupleToMerged, const double postProcThresh)
std::tuple< int, int, std::vector< SimplexId > > trackingTuple
std::vector< PersistencePair > DiagramType
Persistence Diagram type as a vector of Persistence pairs.
printMsg(debug::output::GREEN+"                           "+debug::output::ENDCOLOR+debug::output::GREEN+"▒"+debug::output::ENDCOLOR+debug::output::GREEN+"▒▒▒▒▒▒▒▒▒▒▒▒▒░"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)