TTK
Loading...
Searching...
No Matches
TrackingFromOverlap.h
Go to the documentation of this file.
1
25
26#pragma once
27
28#include <Debug.h>
29#include <algorithm>
30#include <boost/variant.hpp>
31#include <map>
32#include <unordered_map>
33
34using topologyType = unsigned char;
35using idType = long long int;
36
37using labelTypeVariant = boost::variant<double,
38 float,
39 long long,
40 unsigned long long,
41 long,
42 unsigned long,
43 int,
44 unsigned int,
45 short,
46 unsigned short,
47 char,
48 signed char,
49 unsigned char>;
50using sizeType = float;
51
52namespace ttk {
53 class TrackingFromOverlap : virtual public Debug {
54 public:
56 this->setDebugMsgPrefix("TrackingFromOverlap");
57 }
58 ~TrackingFromOverlap() override = default;
59
60 struct Node {
63 float x{};
64 float y{};
65 float z{};
66
70
71 Node() = default;
72 };
73
74 using Edges = std::vector<idType>; // [index0, index1, overlap, branch,...]
75 using Nodes = std::vector<Node>;
76
78 const float *coordinates;
79
80 CoordinateComparator(const float *coords) : coordinates(coords) {
81 }
82
83 inline bool operator()(const size_t &i, const size_t &j) {
84 size_t const ic = i * 3;
85 size_t const jc = j * 3;
86 return coordinates[ic] == coordinates[jc]
87 ? coordinates[ic + 1] == coordinates[jc + 1]
88 ? coordinates[ic + 2] < coordinates[jc + 2]
89 : coordinates[ic + 1] < coordinates[jc + 1]
90 : coordinates[ic] < coordinates[jc];
91 }
92 };
93
94 // This function sorts points based on their x, y, and then z coordinate
95 int sortCoordinates(const float *pointCoordinates,
96 const size_t nPoints,
97 std::vector<size_t> &sortedIndices) const {
98 printMsg("Sorting coordinates ... ", debug::Priority::PERFORMANCE);
99 Timer t;
100
101 sortedIndices.resize(nPoints);
102 for(size_t i = 0; i < nPoints; i++)
103 sortedIndices[i] = i;
104 CoordinateComparator const c = CoordinateComparator(pointCoordinates);
105 sort(sortedIndices.begin(), sortedIndices.end(), c);
106
107 std::stringstream msg;
108 msg << "done (" << t.getElapsedTime() << " s).";
110
111 return 1;
112 }
113
114 int computeBranches(std::vector<Edges> &timeEdgesMap,
115 std::vector<Nodes> &timeNodesMap) const {
116 printMsg("Computing branches ... ", debug::Priority::PERFORMANCE);
117 Timer tm;
118
119 size_t const nT = timeNodesMap.size();
120
121 // Compute max pred and succ
122 for(size_t t = 1; t < nT; t++) {
123 auto &nodes0 = timeNodesMap[t - 1];
124 auto &nodes1 = timeNodesMap[t];
125 auto &edges = timeEdgesMap[t - 1];
126
127 size_t const nE = edges.size();
128
129 for(size_t i = 0; i < nE; i += 4) {
130 auto n0Index = edges[i];
131 auto n1Index = edges[i + 1];
132 auto &n0 = nodes0[n0Index];
133 auto &n1 = nodes1[n1Index];
134
135 sizeType const n0MaxSuccSize
136 = n0.maxSuccID != -1 ? nodes1[n0.maxSuccID].size : 0;
137 sizeType const n1MaxPredSize
138 = n1.maxPredID != -1 ? nodes0[n1.maxPredID].size : 0;
139 if(n0MaxSuccSize < n1.size)
140 n0.maxSuccID = n1Index;
141 if(n1MaxPredSize < n0.size)
142 n1.maxPredID = n0Index;
143 }
144 }
145
146 // Label first nodes of branches
147 idType branchCounter = 0;
148
149 for(size_t t = 0; t < nT; t++)
150 for(auto &n : timeNodesMap[t])
151 n.branchID = n.maxPredID == -1 ? branchCounter++ : -1;
152
153 for(size_t t = 1; t < nT; t++) {
154 auto &nodes0 = timeNodesMap[t - 1];
155 auto &nodes1 = timeNodesMap[t];
156
157 for(size_t i = 0; i < nodes1.size(); i++) {
158 auto &n1 = nodes1[i];
159 if(n1.maxPredID != -1
160 && ((idType)i) != nodes0[n1.maxPredID].maxSuccID)
161 n1.branchID = branchCounter++;
162 }
163 }
164
165 // Propagate branch labels
166 for(size_t t = 1; t < nT; t++) {
167 auto &nodes0 = timeNodesMap[t - 1];
168 auto &nodes1 = timeNodesMap[t];
169 auto &edges = timeEdgesMap[t - 1];
170
171 size_t const nE = edges.size();
172
173 for(size_t i = 0; i < nE; i += 4) {
174 auto n0Index = edges[i];
175 auto n1Index = edges[i + 1];
176 auto &n0 = nodes0[n0Index];
177 auto &n1 = nodes1[n1Index];
178
179 if(n1.branchID == -1 && n0Index == n1.maxPredID)
180 n1.branchID = n0.branchID;
181 }
182 }
183
184 // Label edges
185 for(size_t t = 1; t < nT; t++) {
186 auto &nodes0 = timeNodesMap[t - 1];
187 auto &nodes1 = timeNodesMap[t];
188 auto &edges = timeEdgesMap[t - 1];
189
190 size_t const nE = edges.size();
191
192 for(size_t i = 0; i < nE; i += 4) {
193 auto n0Index = edges[i];
194 auto n1Index = edges[i + 1];
195 auto &n0 = nodes0[n0Index];
196 auto &n1 = nodes1[n1Index];
197
198 edges[i + 3] = n0.branchID == n1.branchID ? n0.branchID
199 : n0.maxSuccID == n1Index ? n0.branchID
200 : n1.branchID;
201 }
202 }
203
204 std::stringstream msg;
205 msg << "done (" << tm.getElapsedTime() << " s).";
207
208 return 1;
209 }
210
211 // This function sorts all unique labels of a point set and then maps these
212 // labels to their respective index in the sorted list
213 template <typename labelType>
214 int computeLabelIndexMap(const labelType *pointLabels,
215 const size_t nPoints,
216 std::map<labelType, size_t> &labelIndexMap) const;
217
218 // This function computes all nodes and their properties based on a labeled
219 // point set
220 template <typename labelType>
221 int computeNodes(const float *pointCoordinates,
222 const labelType *pointLabels,
223 const size_t nPoints,
224 Nodes &nodes) const;
225
226 // This function computes the overlap between two labeled point sets
227 template <typename labelType>
228 int computeOverlap(const float *pointCoordinates0,
229 const float *pointCoordinates1,
230 const labelType *pointLabels0,
231 const labelType *pointLabels1,
232 const size_t nPoints0,
233 const size_t nPoints1,
234
235 Edges &edges) const;
236
237 private:
238 };
239} // namespace ttk
240
241// =============================================================================
242// Compute LabelIndexMap
243// =============================================================================
244template <typename labelType>
246 const labelType *pointLabels,
247 const size_t nPoints,
248 std::map<labelType, size_t> &labelIndexMap) const {
249 for(size_t i = 0; i < nPoints; i++)
250 labelIndexMap[pointLabels[i]] = 0;
251 size_t i = 0;
252 for(auto &it : labelIndexMap)
253 it.second = i++;
254 return 1;
255}
256
257// =============================================================================
258// Identify Nodes
259// =============================================================================
260template <typename labelType>
261int ttk::TrackingFromOverlap::computeNodes(const float *pointCoordinates,
262 const labelType *pointLabels,
263 const size_t nPoints,
264 Nodes &nodes) const {
265 printMsg("Identifying nodes ..... ", debug::Priority::PERFORMANCE);
266
267 Timer t;
268
269 std::map<labelType, size_t> labelIndexMap;
270 this->computeLabelIndexMap(pointLabels, nPoints, labelIndexMap);
271
272 size_t const nNodes = labelIndexMap.size();
273
274 nodes.resize(nNodes);
275 for(size_t i = 0, q = 0; i < nPoints; i++) {
276 labelType label = pointLabels[i];
277 Node &n = nodes[labelIndexMap[label]];
278 n.label = label;
279 n.size++;
280 n.x += pointCoordinates[q++];
281 n.y += pointCoordinates[q++];
282 n.z += pointCoordinates[q++];
283 }
284
285 for(size_t i = 0; i < nNodes; i++) {
286 Node &n = nodes[i];
287 float const size = (float)n.size;
288 n.x /= size;
289 n.y /= size;
290 n.z /= size;
291 }
292
293 // Print Status
294 {
295 std::stringstream msg;
296 msg << "done (#" << nNodes << " in " << t.getElapsedTime() << " s).";
298 }
299
300 return 1;
301}
302
303// =============================================================================
304// Track Nodes
305// =============================================================================
306template <typename labelType>
307int ttk::TrackingFromOverlap::computeOverlap(const float *pointCoordinates0,
308 const float *pointCoordinates1,
309 const labelType *pointLabels0,
310 const labelType *pointLabels1,
311 const size_t nPoints0,
312 const size_t nPoints1,
313
314 Edges &edges) const {
315 // -------------------------------------------------------------------------
316 // Compute labelIndexMaps
317 // -------------------------------------------------------------------------
318 std::map<labelType, size_t> labelIndexMap0;
319 std::map<labelType, size_t> labelIndexMap1;
320 this->computeLabelIndexMap<labelType>(pointLabels0, nPoints0, labelIndexMap0);
321 this->computeLabelIndexMap<labelType>(pointLabels1, nPoints1, labelIndexMap1);
322
323 // -------------------------------------------------------------------------
324 // Sort coordinates
325 // -------------------------------------------------------------------------
326 std::vector<size_t> sortedIndices0;
327 std::vector<size_t> sortedIndices1;
328 this->sortCoordinates(pointCoordinates0, nPoints0, sortedIndices0);
329 this->sortCoordinates(pointCoordinates1, nPoints1, sortedIndices1);
330
331 // -------------------------------------------------------------------------
332 // Track Nodes
333 // -------------------------------------------------------------------------
334 printMsg("Tracking .............. ", debug::Priority::PERFORMANCE);
335 Timer t;
336
337 /* Function that determines configuration of point p0 and p1:
338 0: p0Coords = p1Coords
339 >0: p0Coords < p1Coords
340 <0: p0Coords > p1Coords
341 */
342 auto compare = [&](size_t p0, size_t p1) {
343 size_t p0CoordIndex = p0 * 3;
344 size_t p1CoordIndex = p1 * 3;
345
346 float const p0_X = pointCoordinates0[p0CoordIndex++];
347 float const p0_Y = pointCoordinates0[p0CoordIndex++];
348 float const p0_Z = pointCoordinates0[p0CoordIndex];
349
350 float const p1_X = pointCoordinates1[p1CoordIndex++];
351 float const p1_Y = pointCoordinates1[p1CoordIndex++];
352 float const p1_Z = pointCoordinates1[p1CoordIndex];
353
354 return p0_X == p1_X ? p0_Y == p1_Y ? p0_Z == p1_Z ? 0
355 : p0_Z < p1_Z ? -1
356 : 1
357 : p0_Y < p1_Y ? -1
358 : 1
359 : p0_X < p1_X ? -1
360 : 1;
361 };
362
363 size_t i = 0; // iterator for 0
364 size_t j = 0; // iterator for 1
365
366 size_t nEdges = 0;
367 std::unordered_map<size_t, std::unordered_map<size_t, size_t>> edgesMap;
368 // Iterate over both point sets synchronously using comparison function
369 while(i < nPoints0 && j < nPoints1) {
370 size_t const pointIndex0 = sortedIndices0[i];
371 size_t const pointIndex1 = sortedIndices1[j];
372
373 // Determine point configuration
374 int const c = compare(pointIndex0, pointIndex1);
375
376 if(c == 0) { // Points have same coordinates -> track
377 labelType label0 = pointLabels0[pointIndex0];
378 labelType label1 = pointLabels1[pointIndex1];
379
380 size_t const &nodeIndex0 = labelIndexMap0[label0];
381 size_t const &nodeIndex1 = labelIndexMap1[label1];
382
383 // Find edge and increase overlap counter
384 auto edges0 = edgesMap.find(nodeIndex0); // Edges from label0 to nodes1
385
386 // If map does not exist then create it
387 if(edges0 == edgesMap.end()) {
388 edgesMap[nodeIndex0] = std::unordered_map<size_t, size_t>();
389 edges0 = edgesMap.find(nodeIndex0);
390 }
391
392 // Find edge label0 -> label1
393 auto edge = edges0->second.find(nodeIndex1);
394
395 // If edge does not exist then create it
396 if(edge == edges0->second.end()) {
397 edges0->second[nodeIndex1] = 0;
398 edge = edges0->second.find(nodeIndex1);
399 nEdges++;
400 }
401
402 // Increase overlap
403 edge->second++;
404
405 i++;
406 j++;
407 } else if(c > 0) { // p0 in front of p1 -> let p1 catch up
408 j++;
409 } else { // p1 in front of p0 -> let p0 catch up
410 i++;
411 }
412 }
413
414 // -------------------------------------------------------------------------
415 // Pack Output
416 // -------------------------------------------------------------------------
417 {
418 edges.resize(nEdges * 4);
419 size_t q = 0;
420 for(auto &it0 : edgesMap) {
421 for(auto &it1 : it0.second) {
422 edges[q++] = it0.first;
423 edges[q++] = it1.first;
424 edges[q++] = it1.second;
425 edges[q++] = -1;
426 }
427 }
428 }
429
430 // Print Status
431 {
432 std::stringstream msg;
433 msg << "done (#" << nEdges << " in " << t.getElapsedTime() << " s).";
435 }
436
437 return 0;
438}
boost::variant< double, float, long long, unsigned long long, long, unsigned long, int, unsigned int, short, unsigned short, char, signed char, unsigned char > labelTypeVariant
float sizeType
unsigned char topologyType
long long int idType
Minimalist debugging class.
Definition Debug.h:88
void setDebugMsgPrefix(const std::string &prefix)
Definition Debug.h:364
double getElapsedTime()
Definition Timer.h:15
TTK trackingFromOverlap processing package that tracks labeled point sets.
std::vector< idType > Edges
int computeNodes(const float *pointCoordinates, const labelType *pointLabels, const size_t nPoints, Nodes &nodes) const
int sortCoordinates(const float *pointCoordinates, const size_t nPoints, std::vector< size_t > &sortedIndices) const
int computeBranches(std::vector< Edges > &timeEdgesMap, std::vector< Nodes > &timeNodesMap) const
~TrackingFromOverlap() override=default
int computeOverlap(const float *pointCoordinates0, const float *pointCoordinates1, const labelType *pointLabels0, const labelType *pointLabels1, const size_t nPoints0, const size_t nPoints1, Edges &edges) const
int computeLabelIndexMap(const labelType *pointLabels, const size_t nPoints, std::map< labelType, size_t > &labelIndexMap) const
The Topology ToolKit.
bool operator()(const size_t &i, const size_t &j)
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)