TTK
Loading...
Searching...
No Matches
MetricDistortion.h
Go to the documentation of this file.
1
11
12#pragma once
13
14// ttk common includes
15#include <Debug.h>
16#include <Geometry.h>
17#include <Triangulation.h>
18
19namespace ttk {
20
27 class MetricDistortion : virtual public Debug {
28
29 public:
31
33 // Pre-condition functions.
34 if(triangulation) {
35 triangulation->preconditionBoundaryVertices();
36 triangulation->preconditionVertexNeighbors();
37 }
38 return 0;
39 }
40
41 template <class triangulationType, class tableDataType>
42 void computeSurfaceCurvature(const triangulationType *triangulation,
43 std::vector<tableDataType *> &distanceMatrix,
44 std::vector<double> &surfaceCurvature,
45 std::vector<double> &metricCurvature,
46 std::vector<double> &diffCurvature) {
47 unsigned int dim = triangulation->getNumberOfVertices();
48 surfaceCurvature = std::vector<double>(dim, std::nan(""));
49 metricCurvature = std::vector<double>(dim, std::nan(""));
50 diffCurvature = std::vector<double>(dim, std::nan(""));
51
52 std::vector<std::vector<std::tuple<int, int>>> point2CellPoints(dim);
53 unsigned int noTriangle = 0, noQuad = 0;
54 for(int i = 0; i < triangulation->getNumberOfCells(); ++i) {
55 auto cellNoPoints = triangulation->getCellVertexNumber(i);
56 if(cellNoPoints < 3 or cellNoPoints > 4)
57 continue;
58 for(int j = 0; j < cellNoPoints; ++j) {
59 auto first = (j + 1) % cellNoPoints;
60 auto second = (j + (cellNoPoints - 1)) % cellNoPoints;
61
62 SimplexId i0, i1;
63 triangulation->getCellVertex(i, first, i0);
64 triangulation->getCellVertex(i, second, i1);
65
66 std::tuple<int, int> tup = std::make_tuple(i0, i1);
67
68 SimplexId ij;
69 triangulation->getCellVertex(i, j, ij);
70 point2CellPoints[ij].push_back(tup);
71 }
72 noTriangle += (cellNoPoints == 3);
73 noQuad += (cellNoPoints == 4);
74 }
75
76#ifdef TTK_ENABLE_OPENMP
77#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
78#endif
79 for(unsigned int i = 0; i < dim; ++i) {
80 double sumAngleSurface = 0.0, sumAngleMetric = 0.0;
81
82 for(unsigned int j = 0; j < point2CellPoints[i].size(); ++j) {
83 auto tup = point2CellPoints[i][j];
84 int i0 = std::get<0>(tup);
85 int i1 = std::get<1>(tup);
86
87 float p_i0[3], p_i1[3], p_i[3];
88 triangulation->getVertexPoint(i0, p_i0[0], p_i0[1], p_i0[2]);
89 triangulation->getVertexPoint(i1, p_i1[0], p_i1[1], p_i1[2]);
90 triangulation->getVertexPoint(i, p_i[0], p_i[1], p_i[2]);
91
92 double dist_i_i0 = Geometry::distance(&p_i[0], &p_i0[0]);
93 double dist_i_i1 = Geometry::distance(&p_i[0], &p_i1[0]);
94 double dist_i0_i1 = Geometry::distance(&p_i0[0], &p_i1[0]);
95 double angleSurface;
97 dist_i_i0, dist_i_i1, dist_i0_i1, angleSurface);
98 sumAngleSurface += angleSurface;
99 if(distanceMatrix.size() != 0) {
100 double angleMetric;
101 double distMat_i_i0 = distanceMatrix[i][i0];
102 double distMat_i_i1 = distanceMatrix[i][i1];
103 double distMat_i0_i1 = distanceMatrix[i0][i1];
105 distMat_i_i0, distMat_i_i1, distMat_i0_i1, angleMetric);
106 sumAngleMetric += angleMetric;
107 }
108 }
109
110 unsigned int cornerNoCell = (noTriangle > noQuad ? 2 : 1);
111 double coef = (point2CellPoints[i].size() <= cornerNoCell
112 ? 0.5
113 : (triangulation->isVertexOnBoundary(i) ? 1 : 2));
114 surfaceCurvature[i] = coef * M_PI - sumAngleSurface;
115 // surfaceCurvature[i] *= std::pow(coef, -1);
116
117 if(distanceMatrix.size() != 0) {
118 metricCurvature[i] = coef * M_PI - sumAngleMetric;
119 // metricCurvature[i] *= std::pow(coef, -1);
120 diffCurvature[i] = metricCurvature[i] - surfaceCurvature[i];
121 }
122 }
123 }
124
125 void initIndicators(double &min, double &max, double &avg) {
126 min = std::numeric_limits<double>::max();
127 max = std::numeric_limits<double>::lowest();
128 avg = 0.0;
129 }
130
132 double &min, double &max, double &avg, double value, int noValue) {
133 min = std::min(min, value);
134 max = std::max(max, value);
135 avg += value / noValue;
136 }
137
138 template <class triangulationType, class tableDataType>
140 const triangulationType *triangulation,
141 std::vector<tableDataType *> &distanceMatrix,
142 std::vector<double> &surfaceDistance,
143 std::vector<double> &metricDistance,
144 std::vector<double> &ratioDistance,
145 std::vector<std::array<double, 3>> &surfacePointDistance,
146 std::vector<std::array<double, 3>> &metricPointDistance,
147 std::vector<std::array<double, 3>> &ratioPointDistance) {
148 unsigned int dim = triangulation->getNumberOfCells();
149 surfaceDistance = std::vector<double>(dim, std::nan(""));
150 metricDistance = std::vector<double>(dim, std::nan(""));
151 ratioDistance = std::vector<double>(dim, std::nan(""));
152#ifdef TTK_ENABLE_OPENMP
153#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
154#endif
155 for(unsigned int i = 0; i < dim; ++i) {
156 if(triangulation->getCellVertexNumber(i) != 2)
157 continue;
158
159 SimplexId i0, i1;
160 triangulation->getCellVertex(i, 0, i0);
161 triangulation->getCellVertex(i, 1, i1);
162
163 float p0[3], p1[3];
164 triangulation->getVertexPoint(i0, p0[0], p0[1], p0[2]);
165 triangulation->getVertexPoint(i1, p1[0], p1[1], p1[2]);
166
167 surfaceDistance[i] = Geometry::distance(&p0[0], &p1[0]);
168
169 if(distanceMatrix.size() != 0) {
170 metricDistance[i] = distanceMatrix[i0][i1];
171 ratioDistance[i] = metricDistance[i] / surfaceDistance[i];
172 }
173 }
174
175 dim = triangulation->getNumberOfVertices();
176 surfacePointDistance = std::vector<std::array<double, 3>>(dim);
177 metricPointDistance = std::vector<std::array<double, 3>>(dim);
178 ratioPointDistance = std::vector<std::array<double, 3>>(dim);
179#ifdef TTK_ENABLE_OPENMP
180#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
181#endif
182 for(unsigned int i = 0; i < dim; ++i) {
183 double minDistanceS, maxDistanceS, avgDistanceS;
184 initIndicators(minDistanceS, maxDistanceS, avgDistanceS);
185 double minDistanceM, maxDistanceM, avgDistanceM;
186 initIndicators(minDistanceM, maxDistanceM, avgDistanceM);
187
188 auto neighborNum = triangulation->getVertexNeighborNumber(i);
189 for(int j = 0; j < neighborNum; ++j) {
190 SimplexId neighbor;
191 triangulation->getVertexNeighbor(i, j, neighbor);
192
193 float p0[3], p1[3];
194 triangulation->getVertexPoint(i, p0[0], p0[1], p0[2]);
195 triangulation->getVertexPoint(neighbor, p1[0], p1[1], p1[2]);
196
197 double distance = Geometry::distance(&p0[0], &p1[0]);
199 minDistanceS, maxDistanceS, avgDistanceS, distance, neighborNum);
200
201 if(distanceMatrix.size() != 0) {
202 double distanceM = distanceMatrix[i][neighbor];
204 minDistanceM, maxDistanceM, avgDistanceM, distanceM, neighborNum);
205 }
206 }
207 surfacePointDistance[i][0] = minDistanceS;
208 surfacePointDistance[i][1] = maxDistanceS;
209 surfacePointDistance[i][2] = avgDistanceS;
210 if(distanceMatrix.size() != 0) {
211 metricPointDistance[i][0] = minDistanceM;
212 metricPointDistance[i][1] = maxDistanceM;
213 metricPointDistance[i][2] = avgDistanceM;
214 for(unsigned int j = 0; j < 3; ++j)
215 ratioPointDistance[i][j]
216 = metricPointDistance[i][j] / surfacePointDistance[i][j];
217 }
218 }
219 }
220
221 template <class triangulationType, class tableDataType>
222 void computeSurfaceArea(const triangulationType *triangulation,
223 std::vector<tableDataType *> &distanceMatrix,
224 std::vector<double> &surfaceArea,
225 std::vector<double> &metricArea,
226 std::vector<double> &ratioArea) {
227 unsigned int dim = triangulation->getNumberOfCells();
228 surfaceArea = std::vector<double>(dim, std::nan(""));
229 metricArea = std::vector<double>(dim, std::nan(""));
230 ratioArea = std::vector<double>(dim, std::nan(""));
231#ifdef TTK_ENABLE_OPENMP
232#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
233#endif
234 for(unsigned int i = 0; i < dim; ++i) {
235 auto cellNoPoints = triangulation->getCellVertexNumber(i);
236 if(cellNoPoints < 3 or cellNoPoints > 4)
237 continue;
238
239 SimplexId i0, i1, i2;
240 triangulation->getCellVertex(i, 0, i0);
241 triangulation->getCellVertex(i, 1, i1);
242 triangulation->getCellVertex(i, 2, i2);
243
244 float p0[3], p1[3], p2[3];
245 triangulation->getVertexPoint(i0, p0[0], p0[1], p0[2]);
246 triangulation->getVertexPoint(i1, p1[0], p1[1], p1[2]);
247 triangulation->getVertexPoint(i2, p2[0], p2[1], p2[2]);
248
249 float area;
250 Geometry::computeTriangleArea(&p0[0], &p1[0], &p2[0], area);
251 surfaceArea[i] = area;
252
253 if(distanceMatrix.size() != 0) {
254 double areaMetric;
255 double distMat_i0_i1 = distanceMatrix[i0][i1];
256 double distMat_i1_i2 = distanceMatrix[i1][i2];
257 double distMat_i2_i0 = distanceMatrix[i2][i0];
259 distMat_i0_i1, distMat_i1_i2, distMat_i2_i0, areaMetric);
260 metricArea[i] = areaMetric;
261 }
262
263 if(cellNoPoints == 4) {
264 SimplexId i3;
265 triangulation->getCellVertex(i, 3, i3);
266 float p3[3];
267 triangulation->getVertexPoint(i3, p3[0], p3[1], p3[2]);
268
269 float areaSurface;
270 Geometry::computeTriangleArea(&p1[0], &p2[0], &p3[0], areaSurface);
271 surfaceArea[i] += areaSurface;
272
273 if(distanceMatrix.size() != 0) {
274 double areaMetric;
275 double distMat_i1_i2 = distanceMatrix[i1][i2];
276 double distMat_i2_i3 = distanceMatrix[i2][i3];
277 double distMat_i3_i1 = distanceMatrix[i3][i1];
279 distMat_i1_i2, distMat_i2_i3, distMat_i3_i1, areaMetric);
280 metricArea[i] += areaMetric;
281 }
282 }
283
284 if(distanceMatrix.size() != 0)
285 ratioArea[i] = metricArea[i] / surfaceArea[i];
286 }
287 }
288
289 }; // MetricDistortion class
290
291} // namespace ttk
#define M_PI
Definition: Os.h:50
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
Minimalist debugging class.
Definition: Debug.h:88
int preconditionTriangulation(AbstractTriangulation *triangulation)
void initIndicators(double &min, double &max, double &avg)
void updateIndicators(double &min, double &max, double &avg, double value, int noValue)
void computeSurfaceArea(const triangulationType *triangulation, std::vector< tableDataType * > &distanceMatrix, std::vector< double > &surfaceArea, std::vector< double > &metricArea, std::vector< double > &ratioArea)
void computeSurfaceDistance(const triangulationType *triangulation, std::vector< tableDataType * > &distanceMatrix, std::vector< double > &surfaceDistance, std::vector< double > &metricDistance, std::vector< double > &ratioDistance, std::vector< std::array< double, 3 > > &surfacePointDistance, std::vector< std::array< double, 3 > > &metricPointDistance, std::vector< std::array< double, 3 > > &ratioPointDistance)
void computeSurfaceCurvature(const triangulationType *triangulation, std::vector< tableDataType * > &distanceMatrix, std::vector< double > &surfaceCurvature, std::vector< double > &metricCurvature, std::vector< double > &diffCurvature)
int computeTriangleArea(const T *p0, const T *p1, const T *p2, T &area)
Definition: Geometry.cpp:263
int computeTriangleAreaFromSides(const T s0, const T s1, const T s2, T &area)
Definition: Geometry.cpp:278
int computeTriangleAngleFromSides(const T s0, const T s1, const T s2, T &angle)
Definition: Geometry.cpp:303
T distance(const T *p0, const T *p1, const int &dimension=3)
Definition: Geometry.cpp:344
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22