TTK
Loading...
Searching...
No Matches
MorseSmaleQuadrangulation.cpp
Go to the documentation of this file.
2
3int ttk::MorseSmaleQuadrangulation::findSepsVertices(
4 const std::vector<size_t> &seps,
5 std::vector<LongSimplexId> &srcs,
6 std::vector<LongSimplexId> &dsts) const {
7
8 srcs.resize(seps.size());
9 dsts.resize(seps.size());
10
11 for(size_t i = 0; i < seps.size(); ++i) {
12 auto src = sepCellIds_[sepBegs_[seps[i]]];
13 auto dst = sepCellIds_[sepEnds_[seps[i]]];
14 auto src_dim = sepCellDims_[sepBegs_[seps[i]]];
15 auto dst_dim = sepCellDims_[sepEnds_[seps[i]]];
16 for(LongSimplexId j = 0; j < criticalPointsNumber_; ++j) {
17 if(criticalPointsCellIds_[j] == src
18 && criticalPointsType_[j] == src_dim) {
19 srcs[i] = j;
20 }
21 if(criticalPointsCellIds_[j] == dst
22 && criticalPointsType_[j] == dst_dim) {
23 dsts[i] = j;
24 }
25 }
26 }
27
28 return 0;
29}
30
31int ttk::MorseSmaleQuadrangulation::dualQuadrangulate() {
32
33 // iterate over separatrices middles to build quadrangles around
34 // them: the separatrix vertices and two barycenters
35
36 std::vector<Quad> dualQuads{};
37
38 auto quadNeighbors = [&](const LongSimplexId a) {
39 std::set<LongSimplexId> neighs{};
40 for(const auto &q : outputCells_) {
41 if(a == q[0] || a == q[2]) {
42 neighs.emplace(q[1]);
43 neighs.emplace(q[3]);
44 }
45 if(a == q[1] || a == q[3]) {
46 neighs.emplace(q[0]);
47 neighs.emplace(q[2]);
48 }
49 }
50 return neighs;
51 };
52
53 for(size_t i = 0; i < outputPointsIds_.size(); ++i) {
54 // only keep sep middles
55 if(outputPointsTypes_[i] != 1) {
56 continue;
57 }
58
59 // get a list of direct neighbors
60 auto neighs = quadNeighbors(i);
61
62 // skip sep middle between single extremum and saddle in
63 // degenerate cell (point not used)
64 if(neighs.empty()) {
65 continue;
66 }
67
68 // neighs should contain exactly 4 indices, two corresponding to
69 // critical points, and two for generated points, mostly cell
70 // barycenters
71 std::vector<LongSimplexId> crit{}, gen{};
72 for(const auto n : neighs) {
73 if(n < criticalPointsNumber_) {
74 crit.emplace_back(n);
75 } else {
76 gen.emplace_back(n);
77 }
78 }
79
80 dualQuads.emplace_back(Quad{crit[0], gen[0], crit[1], gen[1]});
81 }
82
83 // for degenerate quadrangles, iterate over the single extremum, and
84 // use the three generated points v0, v1 and v2
85 for(size_t i = 0; i < outputPointsIds_.size(); ++i) {
86
87 // only keep v0 points in degenerate cells
88 if(outputPointsTypes_[i] != 3) {
89 continue;
90 }
91
92 // v0 has 3 neighbors: v1, v2 and the double extremum
93
94 auto v0neighs = quadNeighbors(i);
95 std::vector<LongSimplexId> crit{}, gen{};
96 for(const auto n : v0neighs) {
97 if(n < criticalPointsNumber_) {
98 crit.emplace_back(n);
99 } else {
100 gen.emplace_back(n);
101 }
102 }
103
104 // v1 has 3 neighbors: v0, a sep middle and the single extremum
105
106 // follow v1 to get the single extremum
107 auto v1neighs = quadNeighbors(gen[0]);
108 LongSimplexId single{};
109 for(const auto n : v1neighs) {
110 if(n < criticalPointsNumber_) {
111 single = n;
112 }
113 }
114
115 // the single extremum has 3 neighbors: v1, v2 and the saddle
116
117 // follow single to find saddle
118 auto sneighs = quadNeighbors(single);
119 for(const auto n : sneighs) {
120 if(n < criticalPointsNumber_) {
121 crit.emplace_back(n);
122 }
123 }
124
125 dualQuads.emplace_back(Quad{crit[0], gen[0], crit[1], gen[1]});
126 }
127
128 // clean output points, remove points with no neighbors
129 std::vector<bool> hasNeighs(outputPointsIds_.size(), false);
130 for(const auto &q : dualQuads) {
131 for(const auto &qv : q) {
132 hasNeighs[qv] = true;
133 }
134 }
135
136 std::vector<SimplexId> newPos(outputPointsIds_.size(), -1);
137 SimplexId pos{0};
138 for(size_t i = 0; i < outputPointsIds_.size(); ++i) {
139 if(hasNeighs[i]) {
140 newPos[i] = pos++;
141 }
142 }
143
144 decltype(this->outputPoints_) newPoints(3 * pos);
145 decltype(this->outputPointsIds_) newPointsIds(pos);
146 decltype(this->outputPointsTypes_) newPointsTypes(pos);
147 decltype(this->outputPointsCells_) newPointsCells(pos);
148
149 for(size_t i = 0; i < outputPointsIds_.size(); ++i) {
150 if(hasNeighs[i]) {
151 newPoints[3 * newPos[i] + 0] = this->outputPoints_[3 * i + 0];
152 newPoints[3 * newPos[i] + 1] = this->outputPoints_[3 * i + 1];
153 newPoints[3 * newPos[i] + 2] = this->outputPoints_[3 * i + 2];
154 newPointsIds[newPos[i]] = this->outputPointsIds_[i];
155 newPointsTypes[newPos[i]] = this->outputPointsTypes_[i];
156 newPointsCells[newPos[i]] = this->outputPointsCells_[i];
157 }
158 }
159
160 this->outputPoints_ = std::move(newPoints);
161 this->outputPointsIds_ = std::move(newPointsIds);
162 this->outputPointsTypes_ = std::move(newPointsTypes);
163 this->outputPointsCells_ = std::move(newPointsCells);
164
165 // compact cells
166 for(auto &q : dualQuads) {
167 for(auto &qv : q) {
168 qv = newPos[qv];
169 }
170 }
171
172 outputCells_ = std::move(dualQuads);
173
174 return 0;
175}
176
177void ttk::MorseSmaleQuadrangulation::clearData() {
178 outputCells_.clear();
179 outputPoints_.clear();
180 outputPointsIds_.clear();
181 outputPointsCells_.clear();
182}
long long int LongSimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:15
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22