TTK
Loading...
Searching...
No Matches
PersistentSimplexPairs.h
Go to the documentation of this file.
1
11
12#pragma once
13
15#include <Debug.h>
16#include <VisitedMask.h>
17
18#include <algorithm>
19#include <string>
20#include <vector>
21
22namespace ttk {
23
24 class PersistentSimplexPairs : virtual public Debug {
25 public:
27
34 int type;
35
37 : birth{b}, death{d}, type{t} {
38 }
39 };
40
46 if(data != nullptr) {
47 const auto dim = data->getDimensionality();
48 data->preconditionEdges();
49 if(dim == 2) {
51 } else if(dim == 3) {
55 }
56 this->nVerts_ = data->getNumberOfVertices();
57 this->nEdges_ = data->getNumberOfEdges();
58 this->nTri_ = dim > 1 ? data->getNumberOfTriangles() : 0;
59 this->nTetra_ = dim > 2 ? data->getNumberOfCells() : 0;
60 }
61 }
62
67 template <typename triangulationType>
68 int computePersistencePairs(std::vector<PersistencePair> &pairs,
69 const SimplexId *const orderField,
70 const triangulationType &triangulation) const;
71
72 private:
73 struct Simplex {
74 SimplexId dim_{-1}; // dimension
75 SimplexId id_{-1}; // id in triangulation (overlap between dimensions)
76 SimplexId cellId_{-1}; // cell id (unique)
77 // face (triangulation) indices
78 std::array<SimplexId, 4> faceIds_{-1, -1, -1, -1};
79 // order on vertices, sorted in descending order
80 std::array<SimplexId, 4> vertsOrder_{-1, -1, -1, -1};
81
82 friend bool operator<(const Simplex &lhs, const Simplex &rhs) {
83 return lhs.vertsOrder_ < rhs.vertsOrder_;
84 }
85
86 void fillVert(const SimplexId v, const SimplexId *const offset) {
87 this->dim_ = 0;
88 this->id_ = v;
89 this->cellId_ = v;
90 this->vertsOrder_[0] = offset[v];
91 }
92
93 template <typename triangulationType>
94 void fillEdge(const SimplexId e,
95 const SimplexId c,
96 const SimplexId *const offset,
97 const triangulationType &triangulation) {
98 this->dim_ = 1;
99 this->id_ = e;
100 this->cellId_ = c;
101 triangulation.getEdgeVertex(e, 0, this->faceIds_[0]);
102 triangulation.getEdgeVertex(e, 1, this->faceIds_[1]);
103 this->vertsOrder_[0] = offset[this->faceIds_[0]];
104 this->vertsOrder_[1] = offset[this->faceIds_[1]];
105 std::sort(this->vertsOrder_.rbegin(), this->vertsOrder_.rend());
106 }
107
108 template <typename triangulationType>
109 void fillTriangle(const SimplexId t,
110 const SimplexId c,
111 const SimplexId *const offset,
112 const triangulationType &triangulation) {
113 this->dim_ = 2;
114 this->id_ = t;
115 this->cellId_ = c;
116 triangulation.getTriangleEdge(t, 0, this->faceIds_[0]);
117 triangulation.getTriangleEdge(t, 1, this->faceIds_[1]);
118 triangulation.getTriangleEdge(t, 2, this->faceIds_[2]);
119 triangulation.getTriangleVertex(t, 0, this->vertsOrder_[0]);
120 triangulation.getTriangleVertex(t, 1, this->vertsOrder_[1]);
121 triangulation.getTriangleVertex(t, 2, this->vertsOrder_[2]);
122 this->vertsOrder_[0] = offset[this->vertsOrder_[0]];
123 this->vertsOrder_[1] = offset[this->vertsOrder_[1]];
124 this->vertsOrder_[2] = offset[this->vertsOrder_[2]];
125 std::sort(this->vertsOrder_.rbegin(), this->vertsOrder_.rend());
126 }
127
128 template <typename triangulationType>
129 void fillTetra(const SimplexId T,
130 const SimplexId c,
131 const SimplexId *const offset,
132 const triangulationType &triangulation) {
133 this->dim_ = 3;
134 this->id_ = T;
135 this->cellId_ = c;
136 triangulation.getCellTriangle(T, 0, this->faceIds_[0]);
137 triangulation.getCellTriangle(T, 1, this->faceIds_[1]);
138 triangulation.getCellTriangle(T, 2, this->faceIds_[2]);
139 triangulation.getCellTriangle(T, 3, this->faceIds_[3]);
140 triangulation.getCellVertex(T, 0, this->vertsOrder_[0]);
141 triangulation.getCellVertex(T, 1, this->vertsOrder_[1]);
142 triangulation.getCellVertex(T, 2, this->vertsOrder_[2]);
143 triangulation.getCellVertex(T, 3, this->vertsOrder_[3]);
144 this->vertsOrder_[0] = offset[this->vertsOrder_[0]];
145 this->vertsOrder_[1] = offset[this->vertsOrder_[1]];
146 this->vertsOrder_[2] = offset[this->vertsOrder_[2]];
147 this->vertsOrder_[3] = offset[this->vertsOrder_[3]];
148 std::sort(this->vertsOrder_.rbegin(), this->vertsOrder_.rend());
149 }
150 };
151
152 template <typename triangulationType>
153 std::vector<Simplex>
154 computeFiltrationOrder(const SimplexId *const offset,
155 const triangulationType &triangulation) const;
156
157 inline void addCellBoundary(const Simplex &c, VisitedMask &boundary) const {
158 for(SimplexId i = 0; i < c.dim_ + 1; ++i) {
159 const auto f{c.faceIds_[i]};
160 if(!boundary.isVisited_[f]) {
161 boundary.insert(f);
162 } else {
163 boundary.remove(f);
164 }
165 }
166 }
167
168 inline SimplexId getCellId(const SimplexId cdim,
169 const SimplexId cid) const {
170 if(cdim == 0) {
171 return cid;
172 } else if(cdim == 1) {
173 return cid + this->nVerts_;
174 } else if(cdim == 2) {
175 return cid + this->nVerts_ + this->nEdges_;
176 } else if(cdim == 3) {
177 return cid + this->nVerts_ + this->nEdges_ + this->nTri_;
178 }
179 return -1;
180 }
181
182 SimplexId eliminateBoundaries(const Simplex &c,
183 VisitedMask &boundary,
184 const std::vector<SimplexId> &filtOrder,
185 const std::vector<Simplex> &partners) const;
186
187 int pairCells(std::vector<PersistencePair> &pairs,
188 std::array<std::vector<bool>, 3> &boundaries,
189 const std::vector<Simplex> &filtration,
190 const std::vector<SimplexId> &filtOrder) const;
191
192 SimplexId nVerts_{0};
193 SimplexId nEdges_{0};
194 SimplexId nTri_{0};
195 SimplexId nTetra_{0};
196 };
197
198} // namespace ttk
199
200template <typename triangulationType>
202 std::vector<ttk::PersistentSimplexPairs::PersistencePair> &pairs,
203 const SimplexId *const orderField,
204 const triangulationType &triangulation) const {
205
206 Timer tm{};
207
208 // every simplex in the triangulation, sorted by filtration
209 const auto filtration
210 = this->computeFiltrationOrder(orderField, triangulation);
211
212 std::array<std::vector<bool>, 3> boundaries{};
213 boundaries[0].resize(this->nVerts_, false);
214 boundaries[1].resize(this->nEdges_, false);
215 boundaries[2].resize(this->nTri_, false);
216
217 // simplex id -> filtration order
218 std::vector<SimplexId> filtOrder(filtration.size());
219
220#ifdef TTK_ENABLE_OPENMP
221#pragma omp parallel for num_threads(threadNumber_)
222#endif // TTK_ENABLE_OPENMP
223 for(size_t i = 0; i < filtration.size(); ++i) {
224 filtOrder[filtration[i].cellId_] = i;
225 }
226
227 this->pairCells(pairs, boundaries, filtration, filtOrder);
228
229 this->printMsg("Computed " + std::to_string(pairs.size())
230 + " persistence pair" + (pairs.size() > 1 ? "s" : ""),
231 1.0, tm.getElapsedTime(), 1);
232
233 return 0;
234}
235
236template <typename triangulationType>
237std::vector<ttk::PersistentSimplexPairs::Simplex>
238 ttk::PersistentSimplexPairs::computeFiltrationOrder(
239 const SimplexId *const offset,
240 const triangulationType &triangulation) const {
241
242 Timer tm{};
243
244 std::vector<Simplex> res(this->nVerts_ + this->nEdges_ + this->nTri_
245 + this->nTetra_);
246
247#ifdef TTK_ENABLE_OPENMP
248#pragma omp parallel num_threads(threadNumber_)
249#endif // TTK_ENABLE_OPENMP
250 {
251#ifdef TTK_ENABLE_OPENMP
252#pragma omp for nowait
253#endif // TTK_ENABLE_OPENMP
254 for(SimplexId i = 0; i < this->nVerts_; ++i) {
255 res[i].fillVert(i, offset);
256 }
257
258#ifdef TTK_ENABLE_OPENMP
259#pragma omp for nowait
260#endif // TTK_ENABLE_OPENMP
261 for(SimplexId i = 0; i < this->nEdges_; ++i) {
262 const auto o = this->nVerts_ + i;
263 res[o].fillEdge(i, o, offset, triangulation);
264 }
265
266#ifdef TTK_ENABLE_OPENMP
267#pragma omp for nowait
268#endif // TTK_ENABLE_OPENMP
269 for(SimplexId i = 0; i < this->nTri_; ++i) {
270 const auto o = this->nVerts_ + this->nEdges_ + i;
271 res[o].fillTriangle(i, o, offset, triangulation);
272 }
273
274#ifdef TTK_ENABLE_OPENMP
275#pragma omp for
276#endif // TTK_ENABLE_OPENMP
277 for(SimplexId i = 0; i < this->nTetra_; ++i) {
278 const auto o = this->nVerts_ + this->nEdges_ + this->nTri_ + i;
279 res[o].fillTetra(i, o, offset, triangulation);
280 }
281 }
282
283 TTK_PSORT(this->threadNumber_, res.begin(), res.end());
284
285 this->printMsg(
286 "Computed filtration order", 1.0, tm.getElapsedTime(), this->threadNumber_);
287
288 return res;
289}
#define TTK_PSORT(NTHREADS,...)
Parallel sort macro.
Definition: OpenMP.h:46
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
virtual SimplexId getNumberOfCells() const
virtual SimplexId getNumberOfVertices() const
virtual int getDimensionality() const
virtual SimplexId getNumberOfTriangles() const
virtual SimplexId getNumberOfEdges() const
Minimalist debugging class.
Definition: Debug.h:88
int printMsg(const std::string &msg, const debug::Priority &priority=debug::Priority::INFO, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cout) const
Definition: Debug.h:118
Textbook algorithm to find persistence pairs.
int computePersistencePairs(std::vector< PersistencePair > &pairs, const SimplexId *const orderField, const triangulationType &triangulation) const
Compute the persistence pairs from the triangulation simplicial complex.
void preconditionTriangulation(AbstractTriangulation *const data)
Preprocess all the required connectivity requests on the triangulation.
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
PersistencePair(SimplexId b, SimplexId d, int t)
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)