TTK
Loading...
Searching...
No Matches
ttkMergeTreeBase.cpp
Go to the documentation of this file.
1#include <ttkMergeTreeBase.h>
2#include <ttkUtils.h>
3
7
9 const int cc,
10 vtkPoints *points,
11 vtkUnstructuredGrid *skeletonArcs,
12 ttk::ftm::ArcData &arcData) {
13 auto tree = ftmTree_[cc].tree.getTree(GetTreeType());
14 vtkDataArray *idMapper = connected_components_[cc]->GetPointData()->GetArray(
16 auto arc = tree->getSuperArc(arcId);
17 float point[3];
18 vtkIdType pointIds[2];
19
20 using ttk::SimplexId;
21 const SimplexId downNodeId = tree->getLowerNodeId(arc);
22 const SimplexId l_downVertexId = tree->getNode(downNodeId)->getVertexId();
23 const SimplexId g_downVertexId = idMapper->GetTuple1(l_downVertexId);
24 triangulation_[cc]->getVertexPoint(
25 l_downVertexId, point[0], point[1], point[2]);
26 const double scalarMin = inputScalars_[cc]->GetTuple1(l_downVertexId);
27
28 // Get or create first point of the arc
29 SimplexId nextPointId;
30 if(!arcData.hasPoint(g_downVertexId)) {
31 nextPointId = points->InsertNextPoint(point);
32 arcData.addPoint(g_downVertexId, nextPointId, scalarMin, false);
33 } else {
34 nextPointId = arcData.point_ids[g_downVertexId];
35 }
36
37 pointIds[0] = nextPointId;
38
39 for(const SimplexId vertexId : *arc) {
40 triangulation_[cc]->getVertexPoint(vertexId, point[0], point[1], point[2]);
41 pointIds[1] = points->InsertNextPoint(point);
42 const double scalar = inputScalars_[cc]->GetTuple1(vertexId);
43 arcData.setPoint(pointIds[1], scalar, true);
44
45 const SimplexId nextCell
46 = skeletonArcs->InsertNextCell(VTK_LINE, 2, pointIds);
47 arcData.fillArrayCell(
48 nextCell, arcId, ftmTree_[cc], triangulation_[cc], params_);
49
50 pointIds[0] = pointIds[1];
51 }
52
53 const SimplexId upNodeId = tree->getUpperNodeId(arc);
54 const SimplexId l_upVertexId = tree->getNode(upNodeId)->getVertexId();
55 const SimplexId g_upVertexId = idMapper->GetTuple1(l_upVertexId);
56 triangulation_[cc]->getVertexPoint(
57 l_upVertexId, point[0], point[1], point[2]);
58 const double scalarMax = inputScalars_[cc]->GetTuple1(l_upVertexId);
59
60 // Get or create last point of the arc
61 if(!arcData.hasPoint(g_upVertexId)) {
62 nextPointId = points->InsertNextPoint(point);
63 arcData.addPoint(g_upVertexId, nextPointId, scalarMax, false);
64 } else {
65 nextPointId = arcData.point_ids[g_upVertexId];
66 }
67
68 pointIds[1] = nextPointId;
69
70 const SimplexId nextCell
71 = skeletonArcs->InsertNextCell(VTK_LINE, 2, pointIds);
72 arcData.fillArrayCell(
73 nextCell, arcId, ftmTree_[cc], triangulation_[cc], params_);
74
75 return 1;
76}
77
79 const int cc,
80 vtkPoints *points,
81 vtkUnstructuredGrid *skeletonArcs,
82 ttk::ftm::ArcData &arcData) {
83 auto tree = ftmTree_[cc].tree.getTree(GetTreeType());
84 vtkDataArray *idMapper = connected_components_[cc]->GetPointData()->GetArray(
86 auto arc = tree->getSuperArc(arcId);
87 float point[3];
88 vtkIdType pointIds[2];
89
90 using ttk::SimplexId;
91 const SimplexId downNodeId = tree->getLowerNodeId(arc);
92 const SimplexId l_downVertexId = tree->getNode(downNodeId)->getVertexId();
93 const SimplexId g_downVertexId = idMapper->GetTuple1(l_downVertexId);
94 triangulation_[cc]->getVertexPoint(
95 l_downVertexId, point[0], point[1], point[2]);
96 const double scalarMin = inputScalars_[cc]->GetTuple1(l_downVertexId);
97 // Get or create first point of the arc
98 if(!arcData.hasPoint(g_downVertexId)) {
99 const SimplexId nextPointId = points->InsertNextPoint(point);
100 pointIds[0] = nextPointId;
101 arcData.addPoint(g_downVertexId, nextPointId, scalarMin, false);
102 } else {
103 pointIds[0] = arcData.point_ids[g_downVertexId];
104 }
105
106 const SimplexId upNodeId = tree->getUpperNodeId(arc);
107 const SimplexId l_upVertexId = tree->getNode(upNodeId)->getVertexId();
108 const SimplexId g_upVertexId = idMapper->GetTuple1(l_upVertexId);
109 triangulation_[cc]->getVertexPoint(
110 l_upVertexId, point[0], point[1], point[2]);
111 const double scalarMax = inputScalars_[cc]->GetTuple1(l_upVertexId);
112 // Get or create last point of the arc
113 if(!arcData.hasPoint(g_upVertexId)) {
114 const SimplexId nextPointId = points->InsertNextPoint(point);
115 pointIds[1] = nextPointId;
116 arcData.addPoint(g_upVertexId, nextPointId, scalarMax, false);
117 } else {
118 pointIds[1] = arcData.point_ids[g_upVertexId];
119 }
120
121 const SimplexId nextCell
122 = skeletonArcs->InsertNextCell(VTK_LINE, 2, pointIds);
123 arcData.fillArrayCell(
124 nextCell, arcId, ftmTree_[cc], triangulation_[cc], params_);
125
126 return 1;
127}
128
130 const int cc,
131 vtkPoints *points,
132 vtkUnstructuredGrid *skeletonArcs,
133 ttk::ftm::ArcData &arcData) {
134 auto tree = ftmTree_[cc].tree.getTree(GetTreeType());
135 vtkDataArray *idMapper = connected_components_[cc]->GetPointData()->GetArray(
137 auto arc = tree->getSuperArc(arcId);
138 float point[3];
139 vtkIdType pointIds[2];
140
141 using ttk::SimplexId;
142 const SimplexId downNodeId = tree->getLowerNodeId(arc);
143 const SimplexId l_downVertexId = tree->getNode(downNodeId)->getVertexId();
144 const SimplexId g_downVertexId = idMapper->GetTuple1(l_downVertexId);
145 triangulation_[cc]->getVertexPoint(
146 l_downVertexId, point[0], point[1], point[2]);
147 const double scalarMin = inputScalars_[cc]->GetTuple1(l_downVertexId);
148
149 // Get or create first point of the arc
150 SimplexId nextPointId;
151 if(!arcData.hasPoint(g_downVertexId)) {
152 nextPointId = points->InsertNextPoint(point);
153 arcData.addPoint(g_downVertexId, nextPointId, scalarMin, false);
154 } else {
155 nextPointId = arcData.point_ids[g_downVertexId];
156 }
157
158 pointIds[0] = nextPointId;
159
160 const SimplexId upNodeId = tree->getUpperNodeId(arc);
161 const SimplexId l_upVertexId = tree->getNode(upNodeId)->getVertexId();
162 const SimplexId g_upVertexId = idMapper->GetTuple1(l_upVertexId);
163 triangulation_[cc]->getVertexPoint(
164 l_upVertexId, point[0], point[1], point[2]);
165 const double scalarMax = inputScalars_[cc]->GetTuple1(l_upVertexId);
166
167 const double delta = (scalarMax - scalarMin) / (params_.samplingLvl + 1);
168 double scalarLimit = scalarMin + delta;
169 double scalarAvg = 0;
170
171 // Get or create last point of the arc
172 if(!arcData.hasPoint(g_upVertexId)) {
173 nextPointId = points->InsertNextPoint(point);
174 arcData.addPoint(g_upVertexId, nextPointId, scalarMax, false);
175 } else {
176 nextPointId = arcData.point_ids[g_upVertexId];
177 }
178
179 SimplexId c = 0;
180 float sum[3]{0, 0, 0};
181 for(const SimplexId vertexId : *arc) {
182 triangulation_[cc]->getVertexPoint(vertexId, point[0], point[1], point[2]);
183 const double scalarVertex = inputScalars_[cc]->GetTuple1(vertexId);
184
185 if(scalarVertex < scalarLimit) {
186 sum[0] += point[0];
187 sum[1] += point[1];
188 sum[2] += point[2];
189 scalarAvg += scalarVertex;
190 ++c;
191 } else {
192 if(c) {
193 sum[0] /= c;
194 sum[1] /= c;
195 sum[2] /= c;
196 scalarAvg /= c;
197
198 pointIds[1] = points->InsertNextPoint(sum);
199 arcData.setPoint(pointIds[1], scalarAvg, true);
200 const SimplexId nextCell
201 = skeletonArcs->InsertNextCell(VTK_LINE, 2, pointIds);
202 arcData.fillArrayCell(
203 nextCell, arcId, ftmTree_[cc], triangulation_[cc], params_);
204
205 pointIds[0] = pointIds[1];
206 }
207
208 scalarLimit += delta;
209 sum[0] = 0;
210 sum[1] = 0;
211 sum[2] = 0;
212 scalarAvg = 0;
213 c = 0;
214 }
215 }
216
217 // The up id
218 pointIds[1] = nextPointId;
219
220 const SimplexId nextCell
221 = skeletonArcs->InsertNextCell(VTK_LINE, 2, pointIds);
222 arcData.fillArrayCell(
223 nextCell, arcId, ftmTree_[cc], triangulation_[cc], params_);
224
225 return 1;
226}
227
228int ttkMergeTreeBase::getSegmentation(vtkDataSet *outputSegmentation) {
229 ttk::ftm::VertData vertData;
230 vertData.init(ftmTree_, params_);
231
232 for(int cc = 0; cc < nbCC_; cc++) {
233 auto tree = ftmTree_[cc].tree.getTree(GetTreeType());
234 vtkDataArray *idMapper
235 = connected_components_[cc]->GetPointData()->GetArray(
237 const auto numberOfSuperArcs = tree->getNumberOfSuperArcs();
238 // #pragma omp for
239 for(ttk::ftm::idSuperArc arcId = 0; arcId < numberOfSuperArcs; ++arcId) {
240 vertData.fillArrayPoint(
241 arcId, ftmTree_[cc], triangulation_[cc], idMapper, params_);
242 }
243 }
244
245 vtkPointData *pointData = outputSegmentation->GetPointData();
246 vertData.addArray(pointData, params_);
247 // vertex identifier field should not be propagated (it has caused
248 // downstream bugs in distanceField in oceanVortices)
249 pointData->RemoveArray(ttk::VertexScalarFieldName);
250
251 return 1;
252}
253
254int ttkMergeTreeBase::getSkeletonArcs(vtkUnstructuredGrid *outputSkeletonArcs) {
255 vtkNew<vtkUnstructuredGrid> skeletonArcs{};
256 vtkNew<vtkPoints> const points{};
257
258 ttk::ftm::ArcData arcData;
259 arcData.init(ftmTree_, params_);
260
261 const int samplingLevel = params_.samplingLvl;
262 for(int cc = 0; cc < nbCC_; cc++) {
263 auto tree = ftmTree_[cc].tree.getTree(GetTreeType());
264
265 using ttk::SimplexId;
266 const SimplexId numberOfSuperArcs = tree->getNumberOfSuperArcs();
267#ifndef TTK_ENABLE_KAMIKAZE
268 if(!numberOfSuperArcs) {
269 this->printErr("Error : tree has no super arcs.");
270 return 0;
271 }
272#endif
273
274 for(SimplexId arcId = 0; arcId < numberOfSuperArcs; ++arcId) {
275 const SimplexId numberOfRegularNodes = tree->getArcSize(arcId);
276 if(numberOfRegularNodes > 0 and samplingLevel > 0) {
277 addSampledSkeletonArc(arcId, cc, points, skeletonArcs, arcData);
278 } else if(samplingLevel == -1) {
279 addCompleteSkeletonArc(arcId, cc, points, skeletonArcs, arcData);
280 } else {
281 addDirectSkeletonArc(arcId, cc, points, skeletonArcs, arcData);
282 }
283 }
284 }
285
286 skeletonArcs->SetPoints(points);
287 arcData.addArray(skeletonArcs, params_);
288 outputSkeletonArcs->ShallowCopy(skeletonArcs);
289
290 // const SimplexId p_size = points->GetNumberOfPoints();
291 // const SimplexId s_size = tree->getNumberOfVertices();
292 // cout << "arcs points " << p_size << endl;
293 // cout << "scal points " << s_size << endl;
294 // cout << "nb arcs " << tree->getNumberOfSuperArcs()<< endl;
295 // if(p_size != s_size){
296 // exit(3);
297 // }
298
299 return 1;
300}
301
303 vtkUnstructuredGrid *outputSkeletonNodes) {
304 vtkNew<vtkUnstructuredGrid> skeletonNodes{};
305 vtkNew<vtkPoints> points{};
306
307 ttk::ftm::NodeData nodeData;
308 nodeData.init(ftmTree_, params_);
309 nodeData.setScalarType(inputScalars_[0]->GetDataType());
310
311 for(int cc = 0; cc < nbCC_; cc++) {
312 auto tree = ftmTree_[cc].tree.getTree(GetTreeType());
313 vtkDataArray *idMapper
314 = connected_components_[cc]->GetPointData()->GetArray(
316
317 const auto numberOfNodes = tree->getNumberOfNodes();
318#ifndef TTK_ENABLE_KAMIKAZE
319 if(!numberOfNodes) {
320 this->printErr("Error : tree has no nodes.");
321 return 0;
322 }
323#endif
324
325 for(ttk::ftm::idNode nodeId = 0; nodeId < numberOfNodes; ++nodeId) {
326 const auto node = tree->getNode(nodeId);
327#ifndef TTK_ENABLE_KAMIKAZE
328 if(!node) {
329 this->printMsg({"Error : node ", std::to_string(nodeId), " is null."},
331 return 0;
332 }
333#endif
334 const auto local_vertId = node->getVertexId();
335 float point[3];
336 triangulation_[cc]->getVertexPoint(
337 local_vertId, point[0], point[1], point[2]);
338 const auto nextPoint = points->InsertNextPoint(point);
339 nodeData.fillArrayPoint(
340 nextPoint, nodeId, ftmTree_[cc], idMapper, triangulation_[cc], params_);
341 }
342 }
343
344 ttkUtils::CellVertexFromPoints(skeletonNodes, points);
345
346 vtkPointData *pointData = skeletonNodes->GetPointData();
347 nodeData.addArray(pointData, params_);
348
349 outputSkeletonNodes->ShallowCopy(skeletonNodes);
350
351 return 1;
352}
353
354#ifdef TTK_ENABLE_FTM_TREE_STATS_TIME
355void ttkMergeTreeBase::printCSVStats() {
356 using namespace ttk;
357
358 for(auto &t : ftmTree_) {
359 switch(GetTreeType()) {
360 case ftm::TreeType::Join:
361 cout << "JT" << endl;
362 printCSVTree(t.tree.getTree(ftm::TreeType::Join));
363 break;
364 case ftm::TreeType::Split:
365 cout << "ST" << endl;
366 printCSVTree(t.tree.getTree(ftm::TreeType::Split));
367 break;
368 default:
369 cout << "JT" << endl;
370 printCSVTree(t.tree.getTree(ftm::TreeType::Join));
371 cout << "ST" << endl;
372 printCSVTree(t.tree.getTree(ftm::TreeType::Split));
373 break;
374 }
375 }
376}
377
378void ttkMergeTreeBase::printCSVTree(
379 const ttk::ftm::FTMTree_MT *const tree) const {
380 using namespace ttk::ftm;
381
382 const idSuperArc nbArc = tree->getNumberOfLeaves();
383 cout << "begin; ";
384 for(idSuperArc a = 0; a < nbArc; a++) {
385 cout << tree->getActiveTasks(a).begin << "; ";
386 }
387 cout << endl << "end; ";
388 for(idSuperArc a = 0; a < nbArc; a++) {
389 cout << tree->getActiveTasks(a).end << "; ";
390 }
391 cout << endl << "origing; ";
392 for(idSuperArc a = 0; a < nbArc; a++) {
393 cout << tree->getActiveTasks(a).origin << "; ";
394 }
395 cout << endl;
396}
397#endif
398
399// protected
400
401void ttkMergeTreeBase::identify(vtkDataSet *ds) const {
402 vtkNew<ttkSimplexIdTypeArray> identifiers{};
403 const auto nbPoints = ds->GetNumberOfPoints();
404 identifiers->SetName(ttk::VertexScalarFieldName);
405 identifiers->SetNumberOfComponents(1);
406 identifiers->SetNumberOfTuples(nbPoints);
407
408 for(int i = 0; i < nbPoints; i++) {
409 identifiers->SetTuple1(i, i);
410 }
411
412 ds->GetPointData()->AddArray(identifiers);
413}
ttk::ftm::TreeType GetTreeType() const
std::vector< ttk::Triangulation * > triangulation_
int addCompleteSkeletonArc(const ttk::ftm::idSuperArc arcId, const int cc, vtkPoints *points, vtkUnstructuredGrid *skeletonArcs, ttk::ftm::ArcData &arcData)
int addSampledSkeletonArc(const ttk::ftm::idSuperArc arcId, const int cc, vtkPoints *points, vtkUnstructuredGrid *skeletonArcs, ttk::ftm::ArcData &arcData)
std::vector< vtkDataArray * > inputScalars_
int getSegmentation(vtkDataSet *outputSegmentation)
ttk::ftm::Params params_
int getSkeletonNodes(vtkUnstructuredGrid *outputSkeletonNodes)
int addDirectSkeletonArc(const ttk::ftm::idSuperArc arcId, const int cc, vtkPoints *points, vtkUnstructuredGrid *skeletonArcs, ttk::ftm::ArcData &arcData)
int getSkeletonArcs(vtkUnstructuredGrid *outputSkeletonArcs)
void identify(vtkDataSet *ds) const
std::vector< vtkSmartPointer< vtkDataSet > > connected_components_
std::vector< ttk::ftm::LocalFTM > ftmTree_
static int CellVertexFromPoints(vtkDataSet *const dataSet, vtkPoints *const points)
Definition ttkUtils.cpp:328
void setDebugMsgPrefix(const std::string &prefix)
Definition Debug.h:364
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition Debug.h:149
idNode getNumberOfLeaves() const
Definition FTMTree_MT.h:403
long unsigned int idSuperArc
SuperArc index in vect_superArcs_.
unsigned int idNode
Node index in vect_nodes_.
The Topology ToolKit.
const char VertexScalarFieldName[]
default name for vertex scalar field
Definition DataTypes.h:35
int SimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:22
void fillArrayCell(const SimplexId pos, const idSuperArc arcId, LocalFTM &ftmTree, Triangulation *triangulation, Params params)
void addPoint(const SimplexId globalId, const SimplexId id, const float scalar, const bool reg)
int init(std::vector< LocalFTM > &ftmTree, Params params)
bool hasPoint(const SimplexId vertId)
std::vector< SimplexId > point_ids
void addArray(vtkUnstructuredGrid *skeletonArcs, Params params)
void setPoint(const SimplexId id, const float scalar, const bool reg)
int init(std::vector< LocalFTM > &ftmTree, Params params)
void fillArrayPoint(SimplexId arrIdx, const idNode nodeId, LocalFTM &ftmTree, vtkDataArray *idMapper, Triangulation *triangulation, Params params)
void addArray(vtkPointData *pointData, Params params)
void setScalarType(const int s)
void addArray(vtkPointData *pointData, Params params)
void fillArrayPoint(const idSuperArc arcId, LocalFTM &l_tree, Triangulation *triangulation, vtkDataArray *idMapper, Params params)
int init(std::vector< LocalFTM > &ftmTrees, Params params)
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)