TTK
Loading...
Searching...
No Matches
QuadrangulationSubdivision.h
Go to the documentation of this file.
1
18
19#pragma once
20
21// base code includes
22#include <Dijkstra.h>
23#include <Geometry.h>
24#include <Quadrangulation.h>
25#include <Triangulation.h>
26
27#include <limits>
28#include <set>
29
30namespace ttk {
31
32 class QuadrangulationSubdivision : virtual public Debug {
33
34 public:
36 this->setDebugMsgPrefix("QuadrangulationSubdivision");
37 }
38
39 inline void setSubdivisionLevel(const unsigned int value) {
40 SubdivisionLevel = value;
41 }
42 inline void setRelaxationIterations(const unsigned int value) {
44 }
45 inline void setLockInputExtrema(const bool value) {
46 LockInputExtrema = value;
47 }
48 inline void setLockAllInputVertices(const bool value) {
50 }
51 inline void setShowResError(const bool value) {
52 ShowResError = value;
53 }
54 inline void setHausdorffLevel(const float value) {
55 HausdorffLevel = value;
56 }
57 inline void setInputQuads(void *const address, unsigned int size) {
58 inputQuads_ = static_cast<Quad *>(address);
59 inputQuadNumber_ = size;
60 }
61 inline void setInputVertices(void *const address, unsigned int size) {
62 inputVertices_ = static_cast<Point *>(address);
63 inputVertexNumber_ = size;
64 }
65 inline void setInputVertexIdentifiers(void *const address,
66 unsigned int size) {
67 nearestVertexIdentifier_.resize(size);
68 auto inputVertexIdentifiers = static_cast<SimplexId *>(address);
69 for(size_t i = 0; i < size; i++) {
70 nearestVertexIdentifier_[i] = inputVertexIdentifiers[i];
71 }
72 }
73 inline void
75 if(triangl != nullptr) {
78 }
79 }
80
81 template <typename triangulationType = AbstractTriangulation>
82 int execute(const triangulationType &triangulation);
83
84 private:
86 using Quad = Quadrangulation::Quad;
87
103 template <typename triangulationType>
104 int subdivise(Quadrangulation &qd, const triangulationType &triangulation);
105
117 template <typename triangulationType>
118 SimplexId findEdgeMiddle(const std::array<SimplexId, 2> &e,
119 const triangulationType &triangulation) const;
120
131 SimplexId findQuadBary(std::vector<float> &sum, const Quad &quad) const;
132
136 void clearData();
137
138 template <typename triangulationType>
139 float getBoundingBoxDiagonal(const triangulationType &triangulation) const;
140
141 template <typename triangulationType>
142 void computeHausdorff(std::vector<float> &hausdorff,
143 const Quadrangulation &qd,
144 const triangulationType &triangulation) const;
145
146 protected:
147 // number of vertices in the mesh
149
150 // number of subdivisions of the input quadrangles
151 unsigned int SubdivisionLevel{1};
152 // number of relaxation iterations
153 unsigned int RelaxationIterations{10};
154 // lock input extrema
155 bool LockInputExtrema{false};
156 // lock all input vertices
158 // display result despite error
159 bool ShowResError{false};
160 // Hausdorff warning level
161 float HausdorffLevel{200.F};
162
163 // number of input quadrangles
164 unsigned int inputQuadNumber_{};
165 // input quadrangles
166 Quad *inputQuads_{};
167
168 // number of input points (quad vertices)
169 unsigned int inputVertexNumber_{};
170 // input quadrangle vertices (3D coordinates)
172
173 // array of output quadrangles
174 std::vector<Quad> outputQuads_{};
175 // array of output quadrangle vertices
176 std::vector<Point> outputPoints_{};
177 // array mapping quadrangle neighbors
179 // array of nearest input vertex TTK identifier
180 std::vector<SimplexId> nearestVertexIdentifier_{};
181 // holds geodesic distance to every other quad vertex sharing a quad
182 std::vector<std::vector<float>> vertexDistance_{};
183
184 // array of output quadrangle vertex valences
185 std::vector<SimplexId> outputValences_{};
186 // density around vertices (exp minus euclidean distance between
187 // vertex and its closest neighbor)
188 std::vector<float> outputDensity_{};
189 // quad mesh difformity around vertices (exp minus ratio between
190 // smallest and largest euclidean distance to neighbors)
191 std::vector<float> outputDifformity_{};
192 // array of output quadrangle vertex type
193 // 0 - input (critical) point
194 // 1 - edge middle
195 // 2 - quadrangle barycenter
196 std::vector<SimplexId> outputVertType_{};
197 // array of output vertex subdivision level
198 std::vector<SimplexId> outputSubdivision_{};
199
200 // quadrangles statistics
201 std::vector<float> quadArea_{};
202 std::vector<float> quadDiagsRatio_{};
203 std::vector<float> quadEdgesRatio_{};
204 std::vector<float> quadAnglesRatio_{};
205 std::vector<float> hausdorff_{};
206 };
207} // namespace ttk
208
209// if the package is a pure template typename, uncomment the following line
210// #include <QuadrangulationSubdivision.cpp>
211
212template <typename triangulationType>
213ttk::SimplexId ttk::QuadrangulationSubdivision::findEdgeMiddle(
214 const std::array<SimplexId, 2> &e,
215 const triangulationType &triangulation) const {
216
217 SimplexId midId{};
218 float minValue{std::numeric_limits<float>::infinity()};
219
220 // euclidean barycenter of a and b
221 Point edgeEuclBary = (outputPoints_[e[0]] + outputPoints_[e[1]]) * 0.5F;
222
223 for(size_t i = 0; i < vertexDistance_[e[0]].size(); ++i) {
224 float const m = vertexDistance_[e[0]][i];
225 float const n = vertexDistance_[e[1]][i];
226 // stay on the shortest path between a and b
227 float sum = m + n;
228
229 // skip further computation
230 if(sum > minValue) {
231 continue;
232 }
233
234 if(m != std::numeric_limits<float>::infinity()
235 && n != std::numeric_limits<float>::infinity()) {
236 // try to get the middle of the shortest path
237 sum += std::abs(m - n);
238 }
239
240 // get the euclidean distance to AB
241 Point curr{};
242 triangulation.getVertexPoint(i, curr[0], curr[1], curr[2]);
243 // try to minimize the euclidean distance to AB too
244 sum += Geometry::distance(curr.data(), edgeEuclBary.data());
245
246 // search for the minimizing index
247 if(sum < minValue) {
248 minValue = sum;
249 midId = i;
250 }
251 }
252
253 return midId;
254}
255
256template <typename triangulationType>
257int ttk::QuadrangulationSubdivision::subdivise(
258 Quadrangulation &qd, const triangulationType &triangulation) {
259
260 // temp storage for quad subdivision
261 std::vector<Quad> tmp{};
262
263 Timer tm;
264
265 // avoid reallocation in loop, causing invalid pointers
266 outputPoints_.reserve(outputPoints_.size() * 5);
267
268 vertexDistance_.resize(outputPoints_.size());
269
270 // set & precondition quadrangulation object
271 qd.setInputPoints(this->outputPoints_.size(), this->outputPoints_.data());
272 qd.setInputCells(this->outputQuads_.size(), this->outputQuads_.data());
275
276 // compute shortest distance from every vertex to all other that share a quad
277#ifdef TTK_ENABLE_OPENMP
278#pragma omp parallel for num_threads(threadNumber_)
279#endif // TTK_ENABLE_OPENMP
280 for(size_t i = 0; i < outputPoints_.size(); ++i) {
281
282 // skip if already computed on a coarser subdivision
283 if(vertexDistance_[i].empty()) {
284
285 // do not propagate on the whole mesh
286 std::set<SimplexId> bounds{};
287 const auto ns{qd.getVertexStarNumber(i)};
288 for(SimplexId j = 0; j < ns; ++j) {
289 const auto cid{qd.getVertexStar(i, j)};
290 for(const auto v : this->outputQuads_[cid]) {
291 if(v == static_cast<LongSimplexId>(i)) {
292 continue;
293 }
294 bounds.emplace(nearestVertexIdentifier_[v]);
295 }
296 }
297
298 Dijkstra::shortestPath(nearestVertexIdentifier_[i], triangulation,
299 vertexDistance_[i],
300 {bounds.begin(), bounds.end()});
301 }
302 }
303
304 std::vector<SimplexId> quadBaryId(this->outputQuads_.size());
305 std::vector<float> sum{};
306
307#ifdef TTK_ENABLE_OPENMP
308#pragma omp parallel for num_threads(threadNumber_) firstprivate(sum)
309#endif // TTK_ENABLE_OPENMP
310 for(size_t i = 0; i < this->outputQuads_.size(); ++i) {
311 quadBaryId[i] = this->findQuadBary(sum, this->outputQuads_[i]);
312 }
313
314 std::vector<SimplexId> edgeMidId(qd.getNumberOfEdges());
315
316#ifdef TTK_ENABLE_OPENMP
317#pragma omp parallel for num_threads(threadNumber_)
318#endif // TTK_ENABLE_OPENMP
319 for(SimplexId i = 0; i < qd.getNumberOfEdges(); ++i) {
320 edgeMidId[i] = this->findEdgeMiddle(qd.getEdge(i), triangulation);
321 }
322
323 std::vector<SimplexId> processedEdges(qd.getNumberOfEdges(), -1);
324
325 for(size_t a = 0; a < this->outputQuads_.size(); ++a) {
326 const auto &q{this->outputQuads_[a]};
327
328 const auto processEdge = [&](const SimplexId e) -> SimplexId {
329 if(processedEdges[e] == -1) {
330 const auto midab{edgeMidId[e]};
331 Point pt{};
332 triangulation.getVertexPoint(midab, pt[0], pt[1], pt[2]);
333 /* add new point 3d coordinates to vector of output points */
334 this->outputPoints_.emplace_back(pt);
335 /* new point is an edge middle */
336 this->outputVertType_.emplace_back(1);
337 /* store also TTK identifier of triangular mesh vertex */
338 this->nearestVertexIdentifier_.emplace_back(midab);
339 // store in map
340 processedEdges[e] = this->outputPoints_.size() - 1;
341 }
342 return processedEdges[e];
343 };
344
345 const auto ij{processEdge(qd.getCellEdge(a, 0))};
346 const auto jk{processEdge(qd.getCellEdge(a, 1))};
347 const auto kl{processEdge(qd.getCellEdge(a, 2))};
348 const auto li{processEdge(qd.getCellEdge(a, 3))};
349
350 // barycenter TTK identifier
351 const auto baryid = quadBaryId[a];
352 // barycenter 3D coordinates
353 Point bary{};
354 triangulation.getVertexPoint(baryid, bary[0], bary[1], bary[2]);
355
356 // barycenter index in outputPoints_
357 const LongSimplexId baryIdx = outputPoints_.size();
358 outputPoints_.emplace_back(bary);
359 outputVertType_.emplace_back(2);
360 nearestVertexIdentifier_.emplace_back(baryid);
361
362 // add the four new quads
363 tmp.emplace_back(Quad{q[0], ij, baryIdx, li});
364 tmp.emplace_back(Quad{q[1], jk, baryIdx, ij});
365 tmp.emplace_back(Quad{q[2], kl, baryIdx, jk});
366 tmp.emplace_back(Quad{q[3], li, baryIdx, kl});
367 }
368
369 // output subdivision level
370 auto currSubd = outputSubdivision_.back() + 1;
371 auto subdBeg = outputSubdivision_.size();
372 outputSubdivision_.resize(outputPoints_.size());
373 std::fill(
374 outputSubdivision_.begin() + subdBeg, outputSubdivision_.end(), currSubd);
375
376 this->printMsg("Subdivised " + std::to_string(outputQuads_.size())
377 + " quads into " + std::to_string(tmp.size())
378 + " new quads (" + std::to_string(outputPoints_.size())
379 + " points)",
380 1.0, tm.getElapsedTime(), this->threadNumber_);
381
382 outputQuads_ = std::move(tmp);
383
384 return 0;
385}
386
387template <typename triangulationType>
388float ttk::QuadrangulationSubdivision::getBoundingBoxDiagonal(
389 const triangulationType &triangulation) const {
390
391 std::array<float, 3> pmin{std::numeric_limits<float>::max(),
392 std::numeric_limits<float>::max(),
393 std::numeric_limits<float>::max()};
394 std::array<float, 3> pmax{std::numeric_limits<float>::min(),
395 std::numeric_limits<float>::min(),
396 std::numeric_limits<float>::min()};
397
398 for(SimplexId i = 0; i < triangulation.getNumberOfVertices(); ++i) {
399 std::array<float, 3> p{};
400 triangulation.getVertexPoint(i, p[0], p[1], p[2]);
401 pmax[0] = std::max(pmax[0], p[0]);
402 pmax[1] = std::max(pmax[1], p[1]);
403 pmax[2] = std::max(pmax[2], p[2]);
404 pmin[0] = std::min(pmin[0], p[0]);
405 pmin[1] = std::min(pmin[1], p[1]);
406 pmin[2] = std::min(pmin[2], p[2]);
407 }
408
409 return Geometry::distance(pmin.data(), pmax.data());
410}
411
412template <typename triangulationType>
413void ttk::QuadrangulationSubdivision::computeHausdorff(
414 std::vector<float> &hausdorff,
415 const Quadrangulation &qd,
416 const triangulationType &triangulation) const {
417
418 Timer tm{};
419
420 hausdorff.resize(qd.getNumberOfVertices());
421
422 // compute the minimal distance from every triangulation point to
423 // every quadrangulation point
424
425 // compute triangulation bounding box diagonal
426 const auto bboxDiag = getBoundingBoxDiagonal(triangulation);
427
428 // closest quadrangulation vertex for every triangulation vertex
429 std::vector<SimplexId> nearestQuadVert(triangulation.getNumberOfVertices());
430
431#ifdef TTK_ENABLE_OPENMP
432#pragma omp parallel for num_threads(this->getThreadNumber())
433#endif // TTK_ENABLE_OPENMP
434 for(size_t i = 0; i < nearestQuadVert.size(); ++i) {
435 float minDist{std::numeric_limits<float>::infinity()};
436 std::array<float, 3> p{};
437 triangulation.getVertexPoint(i, p[0], p[1], p[2]);
438
439 for(SimplexId j = 0; j < qd.getNumberOfVertices(); ++j) {
440 std::array<float, 3> q{};
441 qd.getVertexPoint(j, q[0], q[1], q[2]);
442 auto dist = Geometry::distance(p.data(), q.data());
443 if(dist < minDist) {
444 minDist = dist;
445 nearestQuadVert[i] = j;
446 }
447 }
448 }
449
450 // nearest triangulation vertices for each quadrangulation vertex
451 std::vector<std::vector<SimplexId>> nearestTriVerts(qd.getNumberOfVertices());
452 for(SimplexId i = 0; i < triangulation.getNumberOfVertices(); ++i) {
453 nearestTriVerts[nearestQuadVert[i]].emplace_back(i);
454 }
455
456#ifdef TTK_ENABLE_OPENMP
457#pragma omp parallel for num_threads(this->getThreadNumber())
458#endif // TTK_ENABLE_OPENMP
459 for(size_t i = 0; i < nearestTriVerts.size(); ++i) {
460 std::array<float, 3> q{};
461 qd.getVertexPoint(i, q[0], q[1], q[2]);
462 float maxDist{};
463 for(const auto v : nearestTriVerts[i]) {
464 std::array<float, 3> p{};
465 triangulation.getVertexPoint(v, p[0], p[1], p[2]);
466 const auto dist = Geometry::distance(p.data(), q.data());
467 if(dist > maxDist) {
468 maxDist = dist;
469 }
470 }
471 hausdorff[i] = maxDist / bboxDiag / nearestQuadVert.size() * 1e8;
472 }
473
474 this->printMsg("Computed Hausdorff distance", 1.0, tm.getElapsedTime(),
475 this->threadNumber_, debug::LineMode::NEW,
477}
478
479// main routine
480template <typename triangulationType>
482 const triangulationType &triangulation) {
483
484 this->printMsg("Beginning computation...");
485 Timer tm;
486
487 // clear output variables
488 clearData();
489
490 // ensure consistency of dependent options
491 if(LockAllInputVertices) {
492 LockInputExtrema = true;
493 }
494
495 // store input points (MSC critical points)
496 for(size_t i = 0; i < inputVertexNumber_; i++) {
497 outputPoints_.emplace_back(inputVertices_[i]);
498 }
499
500 // copy input quads into vector
501 for(size_t i = 0; i < inputQuadNumber_; i++) {
502 outputQuads_.emplace_back(inputQuads_[i]);
503 }
504
505 // fill outputInfos_ with input data (critical points)
506 outputVertType_.resize(outputPoints_.size());
507 std::fill(outputVertType_.begin(), outputVertType_.end(), 0);
508
509 // fill outputSubdivision with input data
510 outputSubdivision_.resize(outputPoints_.size());
511 std::fill(outputSubdivision_.begin(), outputSubdivision_.end(), 0);
512
513 Quadrangulation qd{};
514 qd.setThreadNumber(this->threadNumber_);
515 qd.setDebugLevel(this->debugLevel_);
516
517 // main loop
518 for(size_t i = 0; i < SubdivisionLevel; i++) {
519 // subdivise each quadrangle by creating five new points, at the
520 // center of each edge (4) and at the barycenter of the four
521 // vertices (1).
522 subdivise(qd, triangulation);
523 }
524
525 qd.setInputPoints(this->outputPoints_.size(), this->outputPoints_.data());
526 qd.setInputCells(this->outputQuads_.size(), this->outputQuads_.data());
527
528 // also needed by computeStatistics
529 qd.preconditionVertexNeighbors();
530 qd.preconditionVertexStars();
531
532 if(this->RelaxationIterations > 0) {
533
534 // smoother mask
535 std::vector<char> mask(this->outputPoints_.size(), 1);
536 if(this->LockAllInputVertices) {
537 // all input vertices (before subdivision)
538 for(size_t i = 0; i < this->inputVertexNumber_; ++i) {
539 mask[i] = 0;
540 }
541 } else if(this->LockInputExtrema) {
542 // extraordinary vertices only (valence != 4)
543 for(SimplexId i = 0; i < qd.getNumberOfVertices(); ++i) {
544 if(qd.isVertexExtraordinary(i)) {
545 mask[i] = 0;
546 }
547 }
548 }
549
551 worker.setDebugLevel(this->debugLevel_);
552 worker.setThreadNumber(this->threadNumber_);
553 worker.execute(reinterpret_cast<float *>(this->outputPoints_.data()),
554 reinterpret_cast<float *>(this->outputPoints_.data()),
555 mask.data(), this->nearestVertexIdentifier_.data(),
556 this->RelaxationIterations, qd, triangulation);
557 }
558
559 qd.computeStatistics(this->outputValences_, this->outputDensity_,
560 this->outputDifformity_, this->quadArea_,
561 this->quadDiagsRatio_, this->quadEdgesRatio_,
562 this->quadAnglesRatio_);
563 this->computeHausdorff(this->hausdorff_, qd, triangulation);
564
565 bool criterion = false;
566 for(size_t i = 0; i < outputPoints_.size(); ++i) {
567 if(outputValences_[i] > 4) {
568 continue;
569 }
570 if(hausdorff_[i] > HausdorffLevel) {
571 criterion = true;
572 break;
573 }
574 }
575
576 if(criterion) {
577 // log, clean & early return
578 this->printErr("The output quadrangulation exceeds the provided Haussdorff "
579 "distance tolerance");
580 if(!ShowResError) {
581 clearData();
582 return 1;
583 }
584 }
585
586 this->printMsg("Produced " + std::to_string(outputQuads_.size())
587 + " quadrangles with " + std::to_string(outputPoints_.size())
588 + " points",
589 1.0, tm.getElapsedTime(), this->threadNumber_);
590
591 return 0;
592}
boost::tuple< double, double > Point
Definition TopoMap.cpp:56
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
virtual SimplexId getNumberOfVertices() const
virtual int setThreadNumber(const int threadNumber)
Definition BaseClass.h:80
Minimalist debugging class.
Definition Debug.h:88
void setDebugMsgPrefix(const std::string &prefix)
Definition Debug.h:364
virtual int setDebugLevel(const int &debugLevel)
Definition Debug.cpp:147
Replacement for std::vector<std::vector<SimplexId>>
TTK processing package for the topological simplification of scalar data.
void setInputVertices(void *const address, unsigned int size)
void setSubdivisionLevel(const unsigned int value)
std::vector< SimplexId > nearestVertexIdentifier_
int execute(const triangulationType &triangulation)
void preconditionTriangulation(AbstractTriangulation *const triangl)
void setInputVertexIdentifiers(void *const address, unsigned int size)
std::vector< std::vector< float > > vertexDistance_
void setRelaxationIterations(const unsigned int value)
void setInputQuads(void *const address, unsigned int size)
SimplexId getVertexStar(const SimplexId v, const int l) const
const std::array< SimplexId, 2 > & getEdge(const SimplexId e) const
ttk::SurfaceGeometrySmoother::Point Point
Ad-hoc vertex coordinates data structure (3 floats)
std::array< ttk::LongSimplexId, 4 > Quad
Ad-hoc quad data structure (4 vertex ids)
SimplexId getNumberOfEdges() const
void setInputCells(const SimplexId cellNumber, void *const quadCells)
SimplexId getVertexStarNumber(const SimplexId v) const
void setInputPoints(const SimplexId pointNumber, void *const pointCoords)
SimplexId getCellEdge(const SimplexId c, const int l)
TTK VTK-filter for smoothing meshes on surfaces.
void preconditionTriangulationSurface(AbstractTriangulation *const triangulation)
double getElapsedTime()
Definition Timer.h:15
int shortestPath(const SimplexId source, const triangulationType &triangulation, std::vector< T > &outputDists, const std::vector< SimplexId > &bounds=std::vector< SimplexId >(), const std::vector< bool > &mask=std::vector< bool >())
Compute the Dijkstra shortest path from source.
Definition Dijkstra.h:26
T distance(const T *p0, const T *p1, const int &dimension=3)
Definition Geometry.cpp:362
The Topology ToolKit.
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
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)