TTK
Loading...
Searching...
No Matches
ContourAroundPoint.hpp
Go to the documentation of this file.
1
14#pragma once
15
17
18#include <cassert>
19#include <cmath>
20#include <limits> // for nan
21#include <map>
22#include <set>
23#include <vector>
24
25namespace ttk {
26 class ContourAroundPoint : virtual public Debug {
27
28 public:
30 this->setDebugMsgPrefix("ContourAroundPoint");
31 }
32
33 ~ContourAroundPoint() override = default;
34
59 template <class triangulationType = AbstractTriangulation>
60 int setInputField(triangulationType *triangulation,
61 void *scalars,
62 double sizeFilter,
63 double radius = 0.);
64
74 int setInputPoints(float *coords,
75 float *scalars,
76 float *isovals,
77 int *flags,
78 std::size_t np);
79
85 template <class scalarT>
86 int execute() const;
87
88 protected:
92 template <class triangulationType = AbstractTriangulation>
94
97 template <class triangulationType = AbstractTriangulation>
98 float compDist2(SimplexId v, SimplexId p) const;
99
100 template <typename scalarT, class triangulationType = AbstractTriangulation>
101 void handleOneInpPt(SimplexId vBeg,
102 float isoval,
103 int flag,
104 float scalar) const;
105
108 template <typename scalarT, class triangulationType = AbstractTriangulation>
109 void extendOutFld(const std::set<SimplexId> &inpEdges,
110 float isoval,
111 int flag) const;
112
115 template <typename scalarT, class triangulationType = AbstractTriangulation>
116 void extendOutPts(const std::vector<SimplexId> &vertices,
117 float isoval,
118 int flag,
119 float scalar) const;
120
122 template <class triangulationType = AbstractTriangulation>
123 double compRadius() {
124 if(_inpFldTriang) {
125 float x, y, z;
126 reinterpret_cast<triangulationType *>(_inpFldTriang)
127 ->getVertexPoint(0, x, y, z);
128 _radius = std::sqrt(x * x + y * y + z * z);
129 } else
130 _radius = std::numeric_limits<double>::signaling_NaN();
131 return _radius;
132 }
133
136 double _radius = 0.;
137
138 /* Input data (field and points) */
139
140 void *_inpFldTriang = nullptr; // scalar field domain
141 void *_inpFldScalars = nullptr; // scalar field image
142 // up to _inpDimMax-dimensional cells appear in _inpFldTriang;
143 // lower-dimensional cells may appear too.
145 // minimum required output region size,
146 // in number of contained input field vertices
147 std::size_t _sizeMin = 0;
148
149 float *_inpPtsCoords = nullptr;
150 float *_inpPtsScalars = nullptr;
151 float *_inpPtsIsovals = nullptr;
152 // 0 if the isovalue for a point is above the point's value;
153 // "point is a minimum is default"
154 int *_inpPtsFlags = nullptr;
155 std::size_t _inpPtsNum = 0;
156
157 /* Output data (contours and centroids) */
158
159 mutable std::vector<LongSimplexId> _outContoursCinfos;
161 mutable std::vector<float> _outContoursCoords;
162 mutable std::vector<float> _outContoursScalars;
163 mutable std::vector<int> _outContoursFlags;
164
165 mutable std::vector<float> _outCentroidsCoords;
166 mutable std::vector<float> _outCentroidsScalars;
167 mutable std::vector<int> _outCentroidsFlags;
168
171 template <class PointIterable, class WeightIterable>
172 static std::
173 array<double, std::tuple_size<typename PointIterable::value_type>::value>
174 average(const PointIterable &pts, const WeightIterable &ws);
175 };
176} // namespace ttk
177
178//----------------------------------------------------------------------------//
179//----------------------------------------------------------------------------//
180
181template <class Triang>
183 void *scalars,
184 double sizeFilter,
185 double radius) {
186
187 if(!triangulation)
188 return -1;
189 if(!scalars)
190 return -2;
191 _inpFldTriang = triangulation;
192 _inpFldScalars = scalars;
193
194 _inpDimMax = triangulation->getDimensionality();
195 if(_inpDimMax < 2 || _inpDimMax > 3)
196 return -3;
197
198 _sizeMin = triangulation->getNumberOfVertices() * sizeFilter / 10000. + 1;
199
200 if(radius == -1.)
201 compRadius();
202 else
203 _radius = radius;
204
205 // Call all the required precondition functions here!
206
207 // for getVertexEdgeNumber, getVertexEdge
208 triangulation->preconditionVertexEdges();
209 // for getVertexNeighbor
210 triangulation->preconditionVertexNeighbors();
211 // for getEdgeVertex
212 triangulation->preconditionEdges();
213 // for getEdgeTriangleNumber, getEdgeTriangle
214 triangulation->preconditionEdgeTriangles();
215 // for getTriangleEdge
216 triangulation->preconditionTriangleEdges();
217
218 return 0;
219}
220
221//----------------------------------------------------------------------------//
222
223template <class scalarT>
225 _outContoursCinfos.resize(0);
226 _outContoursNc = 0;
227 _outContoursCoords.resize(0);
228 _outContoursScalars.resize(0);
229 _outContoursFlags.resize(0);
230
231 _outCentroidsCoords.resize(0);
232 _outCentroidsScalars.resize(0);
233 _outCentroidsFlags.resize(0);
234
235#ifndef TTK_ENABLE_KAMIKAZE
236 if(!_inpFldTriang)
237 return -1;
238 if(!_inpFldScalars)
239 return -2;
240 if(!_inpPtsCoords)
241 return -3;
242 if(!_inpPtsScalars)
243 return -4;
244 if(!_inpPtsIsovals)
245 return -5;
246 if(!_inpPtsFlags)
247 return -6;
248#endif
249
250 Timer timUseObj;
251
252 // The following open-mp processing is only relevant for
253 // embarrassingly parallel algorithms.
254 //#ifdef TTK_ENABLE_OPENMP
255 //#pragma omp parallel for num_threads(threadNumber_)
256 //#endif
257 for(size_t p = 0; p < _inpPtsNum; ++p) {
258 handleOneInpPt<scalarT>(
259 findInpFldVert(p), _inpPtsIsovals[p], _inpPtsFlags[p], _inpPtsScalars[p]);
260 }
261
262 printMsg(ttk::debug::Separator::L2); // horizontal '-' separator
263 printMsg("Complete", 1, timUseObj.getElapsedTime());
264 printMsg(ttk::debug::Separator::L1); // horizontal '=' separator
265 return 1;
266}
267
268//----------------------------------------------------------------------------//
269template <class Triang>
271 float vx, vy, vz;
272 reinterpret_cast<Triang *>(_inpFldTriang)->getVertexPoint(v, vx, vy, vz);
273 const auto pCoords = &(_inpPtsCoords[p * 3]);
274 const float dx = pCoords[0] - vx;
275 const float dy = pCoords[1] - vy;
276 const float dz = pCoords[2] - vz;
277 return dx * dx + dy * dy + dz * dz;
278}
279
280//----------------------------------------------------------------------------//
281template <class Triang>
283 // This implementation is based on a naive nearest neighbor search
284 SimplexId minv = 0;
285 float mind = compDist2(minv, p);
286 const auto nv
287 = reinterpret_cast<Triang *>(_inpFldTriang)->getNumberOfVertices();
288 for(SimplexId v = 1; v < nv; ++v) {
289 const auto d = compDist2(v, p);
290 if(d < mind) {
291 minv = v;
292 mind = d;
293 }
294 }
295 return minv;
296}
297
298//----------------------------------------------------------------------------//
299template <typename scalarT, class Triang>
301 float isoval,
302 int flag,
303 float scalar) const {
304 auto triang = reinterpret_cast<Triang *>(_inpFldTriang);
305 auto inpScalars = reinterpret_cast<const scalarT *>(_inpFldScalars);
306 const bool vBegIsMin = flag == 0;
307
308 // The general idea of the algorithm is to start a (breadth-first) search
309 // on the triangulation starting at the input vertex.
310 // Until the isocontour is met, all vertices are collected to approximate
311 // the (weighted) centroid of the enclosed area/volume.
312 // The edges leading outside of the region are also collected, to compute
313 // the actual isocontour (that intersects these edges).
314
315 auto innerVerts = std::vector<SimplexId>{vBeg};
316 std::set<SimplexId> xEdges; // set data-structure is needed in addOutput
317
318 struct VertEdge {
319 SimplexId v; // where am i
320 SimplexId e; // how did i get here
321 };
322 std::vector<VertEdge> q;
323 std::set<SimplexId> handledEdges;
324
325 auto enqueueNeighbors = [triang, &q, &handledEdges](SimplexId vSrc) {
326 const auto n = triang->getVertexEdgeNumber(vSrc);
327 for(SimplexId i = 0; i < n; ++i) {
328 SimplexId e;
329 triang->getVertexEdge(vSrc, i, e);
330 if(handledEdges.count(e))
331 continue;
332 handledEdges.insert(e);
333 // We assume the local indexing of vertex neighbors is the same
334 // as for the edges (edge 0 leads to neighbor 0, etc.).
335 SimplexId vTgt;
336 triang->getVertexNeighbor(vSrc, i, vTgt);
337 q.push_back({vTgt, e});
338 }
339 };
340
341 enqueueNeighbors(vBeg);
342 while(!q.empty()) {
343 const auto ve = q.back();
344 q.pop_back();
345 const auto v = ve.v;
346 const bool vIsAboveIso = inpScalars[v] > isoval;
347 if(vBegIsMin == vIsAboveIso) { // crossed the contour
348 xEdges.insert(ve.e);
349 } else {
350 innerVerts.push_back(v);
351 enqueueNeighbors(v);
352 }
353 }
354
355 if(innerVerts.size() < _sizeMin)
356 return;
357 extendOutFld<scalarT>(xEdges, isoval, flag);
358 extendOutPts<scalarT>(innerVerts, isoval, flag, scalar);
359}
360
361//----------------------------------------------------------------------------//
362template <typename scalarT, class Triang>
363void ttk::ContourAroundPoint::extendOutFld(const std::set<SimplexId> &inpEdges,
364 float isoval,
365 int flag) const {
366 auto triang = reinterpret_cast<Triang *>(_inpFldTriang);
367 auto inpScalars = reinterpret_cast<const scalarT *>(_inpFldScalars);
368
369 // Every edge e of the given edges maps to one output vertex v;
370 // v is connected to a number of other such vertices.
371 // Zero if e is part of an "antenna", one if e is part of the boundary
372 // of the grid, two in the default 2D case, possibly more in 3D.
373 // These neighbors of v are reached by looking at the triangles that e is a
374 // part of.
375
376 std::map<SimplexId, SimplexId> inpe2outv;
377
378 auto getOrMakeOutVert = [&](SimplexId e) {
379 if(inpe2outv.count(e))
380 return inpe2outv[e];
381
382 // geometry of the vertex (linear interpolation on e)
383 const SimplexId v = _outContoursScalars.size();
384 inpe2outv[e] = v;
385
386 SimplexId p;
387 triang->getEdgeVertex(e, 0, p);
388 SimplexId q;
389 triang->getEdgeVertex(e, 1, q);
390
391 const float pVal = static_cast<float>(inpScalars[p]);
392 const float qVal = static_cast<float>(inpScalars[q]);
393 const double pFac = (qVal - isoval) / (qVal - pVal);
394 const double qFac = 1 - pFac;
395
396 float px, py, pz;
397 triang->getVertexPoint(p, px, py, pz);
398 float qx, qy, qz;
399 triang->getVertexPoint(q, qx, qy, qz);
400 const float x = px * pFac + qx * qFac;
401 const float y = py * pFac + qy * qFac;
402 const float z = pz * pFac + qz * qFac;
403 _outContoursCoords.push_back(x);
404 _outContoursCoords.push_back(y);
405 _outContoursCoords.push_back(z);
406
407 _outContoursScalars.push_back(isoval);
408 _outContoursFlags.push_back(flag);
409 return v;
410 };
411
412 // Each input triangle can contribute to at most one output edge
413 std::set<SimplexId> handledTris;
414
415 if(_inpDimMax == 3)
416 printWrn("Currently only the *edges* of contour surfaces are computed");
417
418 for(const auto e : inpEdges) {
419 const auto nTriLoc = triang->getEdgeTriangleNumber(e);
420 if(nTriLoc == 0)
421 continue; // we do not output "contour points" (0D cells)
422
423 const auto v = getOrMakeOutVert(e);
424
425 for(SimplexId tLoc = 0; tLoc < nTriLoc; ++tLoc) {
426 SimplexId tGlo;
427 triang->getEdgeTriangle(e, tLoc, tGlo);
428 if(handledTris.count(tGlo))
429 continue;
430
431 handledTris.insert(tGlo);
432 // Find the other edge that contains the isovalue
433 for(SimplexId eLoc = 0; eLoc < 3; ++eLoc) {
434 SimplexId eGlo;
435 triang->getTriangleEdge(tGlo, eLoc, eGlo);
436 if(eGlo == e || !inpEdges.count(eGlo))
437 continue;
438
439 // new output edge from v(e) to v(eGlo)
440 _outContoursCinfos.push_back(2);
441 _outContoursCinfos.push_back(v);
442 _outContoursCinfos.push_back(getOrMakeOutVert(eGlo));
443 ++_outContoursNc;
444 } // END target edge
445 } // END triangle
446 } // END source edge
447}
448
449//----------------------------------------------------------------------------//
450template <typename scalarT, class Triang>
452 const std::vector<SimplexId> &vertices,
453 float isoval,
454 int flag,
455 float extremeVal) const {
456
457 auto triang = reinterpret_cast<Triang *>(_inpFldTriang);
458 auto inpScalars = reinterpret_cast<const scalarT *>(_inpFldScalars);
459
460 const double radius = _radius;
461 const bool spherical = radius > 0.;
462
463 // We weight the vertices based on the difference to the extreme value.
464 // Vertices that are close to the extreme value get much weight,
465 // vertices that are close to the isovalue get little weight.
466 const double scalarDiffMax = std::abs(isoval - extremeVal);
467
468 // For now, we use a "full cosine" kernel
469 const double kernelInputScaler = M_PI / scalarDiffMax;
470 // d=0 ... cos(0)=1 ... w=1, d=scalarDiffMax ... cos(pi)=-1 ... w = 0
471 auto scalar_w = [kernelInputScaler](double d) {
472 return (std::cos(d * kernelInputScaler) + 1.) / 2.;
473 };
474
475 // Additionally, each vertex is weighted considering the size of the domain it
476 // "represents" (inversely proportional to the sampling density).
477 // NOTE Currently we assume that whenever we have a spherical domain,
478 // it is given on a regular lon-lat grid.
479 const bool regularLonLat = spherical;
480
481 // We collect the 4D samples and weights to compute the weighted average
482 // in a separate function
483 const auto numVerts = vertices.size();
484 auto pts4d = std::vector<std::array<float, 4>>(numVerts);
485 auto ws = std::vector<double>(numVerts);
486 for(std::size_t i = 0; i < numVerts; ++i) {
487
488 const auto v = vertices[i];
489 const scalarT vSca = inpScalars[v];
490 // Because we only use symmetric kernels for the weighting, the sign of
491 // the difference does not matter.
492 const double scalarDiff = vSca - extremeVal;
493 const double scalarW = scalar_w(scalarDiff);
494
495 float x, y, z;
496 triang->getVertexPoint(v, x, y, z);
497 pts4d[i] = {x, y, z, static_cast<float>(vSca)};
498
499 if(!regularLonLat) {
500 if(_inpDimMax > 2) {
501 // If someone has a lot of time, they can implement the case for a
502 // tetrahedralization.
503 printWrn("Vertex weighting does not consider vertex density of this "
504 "volumetric data");
505 ws[i] = scalarW;
506 } else {
507 // This is just *some* code to fill this gap until it is replaced by sth
508 // better; essentially this module was written for regular spherical
509 // grids.
510 // printMsg("Irregular triangulation");
511
512 // Compute the distance between points a and b.
513 static const auto compDist
514 = [](float ax, float ay, float az, float bx, float by, float bz) {
515 const auto dx = bx - ax;
516 const auto dy = by - ay;
517 const auto dz = bz - az;
518 return std::sqrt(dx * dx + dy * dy + dz * dz);
519 };
520
521 // We use Heron's formula for the area of a triangle, only using the
522 // side lengths a, b, c.
523 static const auto compArea = [](double a, double b, double c) {
524 const double s = (a + b + c) / 2.; // half the perimeter
525 // return std::sqrt(s*(s-a)*(s-b)*(s-c));
526 return std::sqrt(
527 3 * s * s - s * a - s * b
528 - s * c); // numerically more stable (with tested datasets)
529 };
530
531 // We assume a closed triangulation (i.e. no holes) and the vertex
532 // neighbors being ordered in a fan-like fashion.
533 double area = 0.;
534 float x2, y2, z2, x3, y3, z3;
535 SimplexId u;
536 triang->getVertexNeighbor(v, 0, u);
537 triang->getVertexPoint(u, x2, y2, z2);
538 double a = compDist(x, y, z, x2, y2, z2);
539
540 // Backup the first edge for the last triangle.
541 const float x2Bak = x2, y2Bak = y2, z2Bak = z2;
542 const double aBak = a;
543
544 const SimplexId n = triang->getVertexNeighborNumber(v);
545 for(SimplexId j = 1; j < n; ++j) {
546 triang->getVertexNeighbor(v, j, u);
547 triang->getVertexPoint(u, x3, y3, z3);
548 const double b = compDist(x, y, z, x3, y3, z3);
549 const double c = compDist(x2, y2, z2, x3, y3, z3);
550 area += compArea(a, b, c);
551 x2 = x3;
552 y2 = y3;
553 z2 = z3;
554 a = b;
555 }
556
557 const double b = aBak;
558 const double c = compDist(x2, y2, z2, x2Bak, y2Bak, z2Bak);
559 area += compArea(a, b, c);
560 // std::cout<<"area: "<<area<<std::endl;
561 ws[i] = scalarW * area;
562 }
563 } else {
564 // spatial weight 1 at equator, 0 at a pole
565 const double latInRad = std::asin(z / radius);
566 // std::asin is in [-pi/2,+pi/2] => equator at 0 => fits
567 const double spatialW = std::cos(latInRad);
568 ws[i] = scalarW * spatialW;
569 // Why product and not mean of the weights?
570 // -> The resulting weight should be 0 if the vertex is either outside
571 // of the sub/super-level set (scalar weight 0) or its cell size is 0
572 // (spatial weight 0)
573 }
574 }
575
576 const auto res = average(pts4d, ws);
577 auto x = res[0];
578 auto y = res[1];
579 auto z = res[2];
580
581 if(spherical) {
582 const auto radiusCur = std::sqrt(x * x + y * y + z * z);
583 const auto radiusScaler = radius / radiusCur;
584 x *= radiusScaler;
585 y *= radiusScaler;
586 z *= radiusScaler;
587 }
588
589 // std::cout<<x<<","<<y<<","<<z<<std::endl;
590 _outCentroidsCoords.push_back(static_cast<float>(x));
591 _outCentroidsCoords.push_back(static_cast<float>(y));
592 _outCentroidsCoords.push_back(static_cast<float>(z));
593 _outCentroidsScalars.push_back(static_cast<float>(res[3]));
594 _outCentroidsFlags.push_back(flag);
595}
596
597//----------------------------------------------------------------------------//
598template <class PointIterable, class WeightIterable>
599std::array<double, std::tuple_size<typename PointIterable::value_type>::value>
600 ttk::ContourAroundPoint::average(const PointIterable &pts,
601 const WeightIterable &ws) {
602
603 constexpr auto numComponents
604 = std::tuple_size<typename PointIterable::value_type>::value;
605 std::array<double, numComponents> res{};
606 // as of C++14 there is no (simple) compile time "constexpr for loop"
607#ifndef NDEBUG
608 for(std::size_t j = 0; j < numComponents; ++j)
609 // assert(std::get<j>(res) == 0.);
610 assert(res[j] == 0.);
611#endif
612 double wSum = 0.;
613
614 const auto numPts = pts.size();
615 assert(ws.size() == numPts);
616
617 for(std::size_t i = 0; i < numPts; ++i) {
618 auto &pt = pts[i];
619 const auto w = ws[i];
620 wSum += w;
621 for(std::size_t j = 0; j < numComponents; ++j)
622 // std::get<j>(res) += std::get<j>(pt) * w;
623 res[j] += pt[j] * w;
624 }
625
626 assert(wSum != 0.);
627 for(std::size_t j = 0; j < numComponents; ++j)
628 // std::get<j>(res) /= wSum;
629 res[j] /= wSum;
630
631 return res;
632}
int ImplicitTriangulationCRTP< Derived >::TTK_TRIANGULATION_INTERNAL() getVertexPoint(const SimplexId &vertexId, float &x, float &y, float &z) const
#define M_PI
Definition: Os.h:50
TTK contourAroundPoint processing package.
float compDist2(SimplexId v, SimplexId p) const
int setInputField(triangulationType *triangulation, void *scalars, double sizeFilter, double radius=0.)
std::vector< float > _outCentroidsCoords
~ContourAroundPoint() override=default
void extendOutFld(const std::set< SimplexId > &inpEdges, float isoval, int flag) const
std::vector< float > _outContoursScalars
static std::array< double, std::tuple_size< typename PointIterable::value_type >::value > average(const PointIterable &pts, const WeightIterable &ws)
void handleOneInpPt(SimplexId vBeg, float isoval, int flag, float scalar) const
std::vector< float > _outContoursCoords
std::vector< int > _outCentroidsFlags
std::vector< float > _outCentroidsScalars
double compRadius()
Set and return _radius.
std::vector< LongSimplexId > _outContoursCinfos
int setInputPoints(float *coords, float *scalars, float *isovals, int *flags, std::size_t np)
void extendOutPts(const std::vector< SimplexId > &vertices, float isoval, int flag, float scalar) const
SimplexId findInpFldVert(SimplexId p) const
std::vector< int > _outContoursFlags
Minimalist debugging class.
Definition: Debug.h:88
void setDebugMsgPrefix(const std::string &prefix)
Definition: Debug.h:364
double getElapsedTime()
Definition: Timer.h:15
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
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)