TTK
Searching...
No Matches
TopoMap.h
Go to the documentation of this file.
1
20
22
23#pragma once
24
25// ttk common includes
26#include <Debug.h>
27#include <Geometry.h>
28#include <UnionFind.h>
29
30// STL includes
31#include <map>
32#include <queue>
33#include <set>
34#include <sstream>
35#include <utility> // For std::pair
36#include <vector>
37
38// cos and sin computations, so we are cautious on testing equalities between
39// double.
40static double Epsilon{ttk::Geometry::powInt(10.0, -FLT_DIG + 2)};
41static double EpsilonDBL{ttk::Geometry::powInt(10.0, -DBL_DIG + 2)};
42
43// Normalizes a given vector.
44template <typename T>
45static void
46 computeUnitVector(const T *coordOrig, const T *coordDest, T *coordVect);
47
48// Rotates the first argument by angle around center.
49template <typename T>
50static void rotate(T *ptToRotate, const T *center, double angle);
51
52// Rotate the set of points by the given angle, considering the given center as
53// center of the rotation.
54template <typename T>
55static void
56 rotatePolygon(std::vector<T> &coords, T *centerCoords, const double angle);
57
58// Real call to the convex hull engine: either Qhull or Boost.
59bool computeConvexHull_aux(const std::vector<double> &coords,
60 std::vector<size_t> &res,
61 std::string &errMsg);
62
63namespace ttk {
64
69 class TopoMap : virtual public Debug {
70
71 public:
75 enum class STRATEGY { KRUSKAL = 0, PRIM = 1 };
76
77 TopoMap();
78 TopoMap(size_t angularSampleNb, bool checkMST, STRATEGY strategy)
79 : AngularSampleNb(angularSampleNb), CheckMST(checkMST),
80 Strategy(strategy) {
81
82 // inherited from Debug: prefix will be printed at the beginning of every
83 // msg
84 this->setDebugMsgPrefix("TopoMap");
85 }
86 ~TopoMap() override;
87
107 template <typename T>
108 int execute(T *outputCoords,
109 int *insertionTime,
110 const std::vector<T> &inputMatrix,
111 bool isDistMat,
112 size_t n);
113
114 protected:
116 bool CheckMST{false};
117 bool errorConvexHull{false};
119
120 private:
121 // Fills its last two arguments with the position of the end of the previous
122 // (resp. next) edge so that idCenter is the point at the angle at the
123 // intersection of the previous and next edges.
124 template <typename T>
125 bool getPrevNextEdges(const std::vector<size_t> &idsPtsPolygon,
126 size_t idCenter,
127 const T *allCoords,
128 T *coordPrev,
129 T *coordPost) const;
130
131 // Tries to find the best angle of rotation for the two components. Updates
132 // the coordiates of their vertices accordingly.
133 template <typename T>
134 T rotateMergingCompsBest(const std::vector<size_t> &hull1,
135 const std::vector<size_t> &hull2,
136 const std::vector<size_t> &comp1,
137 const std::vector<size_t> &comp2,
138 size_t iPt1,
139 size_t iPt2,
140 const std::vector<T> &distMatrix,
141 T *allCoords,
142 size_t n,
143 size_t angularSampleNb);
144
145 template <typename T>
146 bool computeConvexHull(T *allCoords,
147 const std::vector<size_t> &compPtsIds,
148 std::vector<size_t> &idsInHull) const;
149 };
150
151 // Sketch of the algorithm:
152 // 1. Sort the edges
153 // 2. For each edge uv with cost c_uv
154 // a. Get connected components comp(u), comp(v)
155 //
156 // b. Compute convex hull of comp(u) and of (comp(v))
157 //
158 // c. Find in hull(u) (resp. hull(v)) the point pU (resp. pV) in the
159 // convex hull which is closest to u (resp. v) in high dimension, u if
160 // possible.
161 //
162 // d. Translate comp(u) such that pU lies on the line of the bissector of
163 // angleV, at distance c_uv of pV, and on the opposite direction of
164 // hull(v). angleU (resp. angleV) is the angle of hull(u) (resp. hull(v))
165 // at vertex pU (resp. pV).
166 //
167 // e. Rotate comp(u) such that the bissector of angleU prolongates the
168 // one of angleU, and such that pU and pV are closest than any pair of
169 // points in comp(u),comp(v): the two components are facing each other at
170 // pU and pV.
171 //
172 //
173 // f. Try several rotations of comp(u) and comp(v) such that pU and pV
174 // still are the closest pair of vertices from the two components. We
175 // keep the angle of rotation which minimizes the difference between the
176 // high dimension distance matrix and the new distance matrix between the
177 // two components.
178 //
179 // g. comp(u) and comp(v) were merged, we now consider them as a single
180 // component.
181
182 template <typename T>
183 int ttk::TopoMap::execute(T *outputCoords,
184 int *insertionTime,
185 const std::vector<T> &inputMatrix,
186 bool isDistMat,
187 size_t n) {
188 ttk::Timer timer;
189
190#ifdef TTK_ENABLE_QHULL
191#ifndef Qhull_FOUND
192 this->printWrn(
193 "Qhull was enabled but it is not installed or was not found.");
194 this->printWrn("Defaulting to Boost support instead.");
195#endif
196#endif
197
198#ifndef Qhull_FOUND
199 this->printMsg("Using Boost for convex hulls.");
200 this->printWrn("Bugs have been reported in Boost's implementation.");
202#endif
203
204 if(AngularSampleNb < 2) {
205 this->printWrn("The number of angular samples set is less than 2.");
206 this->printWrn("Setting it to 2 (minimum possible value).");
207 this->AngularSampleNb = 2;
208 }
209 // We first set the minimum value we consider as null. We multiply by 100
210 // the smallest significative value because we use sqrt/cos/sin so we take
211 // precautionary measures.
212 if(std::is_same<T, double>::value) {
213 Epsilon = ttk::Geometry::powInt(10.0, -DBL_DIG + 2);
214 } else if(std::is_same<T, float>::value) {
215 Epsilon = ttk::Geometry::powInt(10.0, -FLT_DIG + 2);
216 } else {
217 printErr("Error, template type must be either double or float.\n");
218 return 1;
219 }
220
221 if(isDistMat) {
222 this->printMsg("Input data: " + std::to_string(n) + " points.");
223 } else {
224 this->printMsg("Input data: " + std::to_string(n) + " points ("
225 + std::to_string(inputMatrix.size() / n)
226 + " dimensions).");
227 }
228
229 std::vector<T> computedDistMatrix(n * n);
230 if(!isDistMat) {
231 size_t dim = inputMatrix.size() / n;
232 if(n * dim != inputMatrix.size()) {
233 this->printErr("Error, the coordinates input matrix has invalid size.");
234 return 1;
235 }
236
237#ifdef TTK_ENABLE_OPENMP
239 shared(computedDistMatrix)
240#endif
241 for(size_t i1 = 0; i1 < n; i1++) {
242 for(size_t i2 = i1; i2 < n; i2++) {
243 T dist = ttk::Geometry::distance<T>(
244 &inputMatrix[dim * i1], &inputMatrix[dim * i2], dim);
245 computedDistMatrix[i1 * n + i2] = dist;
246 computedDistMatrix[i2 * n + i1] = dist;
247 }
248 }
249 } else {
250 if(inputMatrix.size() != n * n) {
251 this->printErr("Invalid size for the distance matrix.");
252 return 1;
253 }
254 }
255
256 std::vector<double> edgesMSTBefore,
257 edgesMSTAfter; // Only used for minimum spanning tree checking;
258 edgesMSTBefore.reserve(n);
259 const std::vector<T> &distMatrix
260 = isDistMat ? inputMatrix : computedDistMatrix;
261
262 if(insertionTime != nullptr)
263 for(size_t i = 0; i < n; i++)
264 insertionTime[i] = -1;
265
266 // 1. Sorting the edges
267 using heapType = std::pair<T, std::pair<size_t, size_t>>;
268 std::priority_queue<heapType, std::vector<heapType>, std::greater<heapType>>
269 edgeHeap;
270 if(this->Strategy == STRATEGY::KRUSKAL) {
271 for(size_t u1 = 0; u1 < n; u1++) {
272 for(size_t u2 = u1 + 1; u2 < n; u2++) {
273 edgeHeap.push(
274 std::make_pair(distMatrix[u1 * n + u2], std::make_pair(u1, u2)));
275 }
276 }
277 } else if(this->Strategy == STRATEGY::PRIM) {
278 for(size_t u1 = 0; u1 < n; u1++) {
279 edgeHeap.push(std::make_pair(distMatrix[u1], std::make_pair(0, u1)));
280 }
281 } else {
282 this->printErr("Invalid stategy for the MST: only Kruskal (0) or Prim "
283 "(1) algorithm can be used.");
284 return 1;
285 }
286
287 // ufVectors is the union find vector
288 std::map<UnionFind *, std::vector<size_t>> ufToComp;
289 std::vector<UnionFind> ufVector(n);
290 for(size_t i = 0; i < n; i++) {
291 ufToComp[&ufVector[i]].emplace_back(i);
292 }
293
294 for(size_t i = 0; i < 2 * n; i++) {
295 outputCoords[i] = 0;
296 }
297
298 T finalDistortion = 0;
299 // 2. Processing the edges
300 size_t nbEdgesMerged = 0;
301 double MSTWeight = 0;
302
303 std::set<size_t>
304 stillToDo; // List of vertices not yet embedded, used for Prim strategy.
305 if(this->Strategy == STRATEGY::PRIM) {
306 for(size_t i = 1; i < n; i++) {
307 stillToDo.insert(i);
308 }
309 }
310 while(nbEdgesMerged != n - 1) {
311 if(edgeHeap.empty()) {
312 this->printErr("Error building the MST. Aborting.");
313 return 1;
314 }
315 const auto &elt = edgeHeap.top();
316 T edgeCost = elt.first;
317 size_t u = elt.second.first;
318 size_t v = elt.second.second;
319 edgeHeap.pop();
320
321 // 2.a Getting the components
322 UnionFind *reprU = ufVector[u].find();
323 UnionFind *reprV = ufVector[v].find();
324 if(reprU == reprV) { // Already in the same component
325 continue;
326 }
327 nbEdgesMerged++;
328 MSTWeight += edgeCost;
329
330 if(insertionTime != nullptr) {
331 if(insertionTime[u] == -1) {
332 insertionTime[u] = nbEdgesMerged;
333 }
334 if(insertionTime[v] == -1) {
335 insertionTime[v] = nbEdgesMerged;
336 }
337 }
338 edgesMSTBefore.push_back(edgeCost);
339 std::vector<size_t> &compU = ufToComp[reprU], &compV = ufToComp[reprV];
340 size_t idSmall = compU.size() < compV.size() ? u : v;
341 size_t idBig = idSmall == u ? v : u;
342 std::vector<size_t> &compSmall = idSmall == u ? compU : compV;
343 std::vector<size_t> &compBig = idSmall == u ? compV : compU;
344 size_t nBig = compBig.size();
345 size_t nSmall = compSmall.size();
346 std::vector<size_t> idsInHullSmall, idsInHullBig;
347
348 // 2.b Computing the convex hull.
349 // We retrieve the current coordinates of the big component..
350 bool statusHull = true;
351 statusHull
352 = statusHull | computeConvexHull(outputCoords, compBig, idsInHullBig);
353 statusHull = statusHull
354 | computeConvexHull(outputCoords, compSmall, idsInHullSmall);
355 if(!statusHull) {
356 this->printErr("Aborting.");
357 return 1;
358 }
359
360 size_t idChosenBig = idsInHullBig[0], idChosenSmall = idsInHullSmall[0];
361
362 // 2.c We want to select, among all vertices in the convex hull, the one
363 // which is closest to the vertex of the edge we work on.
364 for(size_t vert : idsInHullBig) {
365 T dist = distMatrix[vert * n + idBig];
366 if(dist < distMatrix[idChosenBig * n + idBig]) {
367 idChosenBig = vert;
368 }
369 }
370
371 for(size_t vert : idsInHullSmall) {
372 T dist = distMatrix[vert * n + idSmall];
373 if(dist < distMatrix[idChosenSmall * n + idSmall]) {
374 idChosenSmall = vert;
375 }
376 }
377
378 size_t sizeBigHull = idsInHullBig.size(),
379 sizeSmallHull = idsInHullSmall.size();
380 std::vector<T> coordsBigHull(sizeBigHull * 2),
381 coordsSmallHull(sizeSmallHull * 2);
382 for(size_t iHull = 0; iHull < sizeBigHull; iHull++) {
383 size_t vert = idsInHullBig[iHull];
384 coordsBigHull[iHull * 2] = outputCoords[vert * 2];
385 coordsBigHull[iHull * 2 + 1] = outputCoords[vert * 2 + 1];
386 }
387 for(size_t iHull = 0; iHull < sizeSmallHull; iHull++) {
388 size_t vert = idsInHullSmall[iHull];
389 coordsSmallHull[iHull * 2] = outputCoords[vert * 2];
390 coordsSmallHull[iHull * 2 + 1] = outputCoords[vert * 2 + 1];
391 }
392
393 // Identifying the angles we are working on from the convex hulls.
394 T coordPrevBig[2], coordPostBig[2];
395 T coordPrevSmall[2], coordPostSmall[2];
396 if(!getPrevNextEdges(idsInHullSmall, idChosenSmall, outputCoords,
397 coordPrevSmall, coordPostSmall)) {
398 return 1;
399 }
400 if(!getPrevNextEdges(idsInHullBig, idChosenBig, outputCoords,
401 coordPrevBig, coordPostBig)) {
402 return 1;
403 }
404 T coordPtSmall[2] = {
405 outputCoords[2 * idChosenSmall], outputCoords[2 * idChosenSmall + 1]};
406 T coordPtBig[2]
407 = {outputCoords[2 * idChosenBig], outputCoords[2 * idChosenBig + 1]};
408
409 // Computing the angles.
410 double angleSmall = ttk::Geometry::angle2DUndirected(
411 coordPtSmall, coordPrevSmall, coordPostSmall);
412 double angleBig = ttk::Geometry::angle2DUndirected(
413 coordPtBig, coordPrevBig, coordPostBig);
414 if(angleSmall - M_PI > EpsilonDBL || angleBig - M_PI > EpsilonDBL) {
415 this->printErr("Error, angle out of bound (greater than pi).");
416 }
417 T coordscenterBig[2] = {coordPrevBig[0], coordPrevBig[1]};
418 T coordscenterSmall[2] = {coordPrevSmall[0], coordPrevSmall[1]};
419 // Computing the coordinates of the bisectors.
420 rotate(coordscenterSmall, coordPtSmall, angleSmall / 2);
421 rotate(coordscenterBig, coordPtBig, angleBig / 2);
422
423 T unitcenterBigVect[2], unitcenterSmallVect[2];
424 computeUnitVector(
425 &outputCoords[idChosenBig * 2], coordscenterBig, unitcenterBigVect);
426 computeUnitVector(&outputCoords[idChosenSmall * 2], coordscenterSmall,
427 unitcenterSmallVect);
428
429 // Computing the new coordinates for the chosen point on the small
430 // components, and then translating the whole component.
431 T goalCoordChosenSmall[2]
432 = {outputCoords[idChosenBig * 2] - edgeCost * unitcenterBigVect[0],
433 outputCoords[idChosenBig * 2 + 1] - edgeCost * unitcenterBigVect[1]};
434 if(sizeBigHull == 1) {
435 goalCoordChosenSmall[0] = outputCoords[idChosenBig * 2] + edgeCost;
436 goalCoordChosenSmall[1] = outputCoords[idChosenBig * 2 + 1];
437 }
438 T smallCompMoveVect[2]
439 = {goalCoordChosenSmall[0] - outputCoords[idChosenSmall * 2],
440 goalCoordChosenSmall[1] - outputCoords[idChosenSmall * 2 + 1]};
441
442 T distBaryPointSmall = ttk::Geometry::distance2D(
443 &outputCoords[idChosenSmall * 2], coordscenterSmall);
444 T preFinalPosBarySmall[2] = {coordscenterSmall[0] + smallCompMoveVect[0],
445 coordscenterSmall[1] + smallCompMoveVect[1]};
446 T finalPosBarySmall[2]
447 = {goalCoordChosenSmall[0] - unitcenterBigVect[0] * distBaryPointSmall,
448 goalCoordChosenSmall[1] - unitcenterBigVect[1] * distBaryPointSmall};
449
450 // Performing the translation + rotation such that the bisectors in
451 // compSmall and compBig are aligned, and the two components are facing
452 // each other, not crossing.
453 for(size_t curIdSmall : compSmall) {
454 outputCoords[curIdSmall * 2] += smallCompMoveVect[0];
455 outputCoords[curIdSmall * 2 + 1] += smallCompMoveVect[1];
456
457 // 2.e Performing the rotation.
458 double rotationAngle = ttk::Geometry::angle2DUndirected(
459 goalCoordChosenSmall, preFinalPosBarySmall, finalPosBarySmall);
460 if(nSmall > 1 && std::isfinite(rotationAngle)) {
461 rotate(
462 &outputCoords[curIdSmall * 2], goalCoordChosenSmall, rotationAngle);
463 }
464 }
465
466 // 2.f Trying several rotations of the two components and keeping the one
467 // which makes the new distance matrix closest to the one provided in the
468 // input.
469 if(nBig > 1) {
470 finalDistortion = rotateMergingCompsBest(
471 idsInHullSmall, idsInHullBig, compSmall, compBig, idChosenSmall,
472 idChosenBig, distMatrix, outputCoords, n, this->AngularSampleNb);
473 if(finalDistortion < -1) {
474 return 1;
475 }
476 }
477
478 T finalDist = ttk::Geometry::distance2D(
479 &outputCoords[2 * idChosenSmall], &outputCoords[2 * idChosenBig]);
480 if(fabs(finalDist - edgeCost) > Epsilon) {
481 this->printErr(
482 "The distance we set is too far from the goal distance.");
483 }
484
485 UnionFind *unionRepr = UnionFind::makeUnion(reprU, reprV);
486 UnionFind *otherRepr = (unionRepr == reprU) ? reprV : reprU;
487 if(unionRepr != reprU && unionRepr != reprV) {
488 printErr("Error with the union find module results. Aborting.");
489 return 1;
490 }
491 std::vector<size_t> &unionSet = ufToComp[unionRepr];
492 std::vector<size_t> &otherSet = ufToComp[otherRepr];
493
494 unionSet.insert(unionSet.end(), otherSet.begin(), otherSet.end());
495
496 ufToComp.erase(otherRepr);
497
498 if(this->Strategy == STRATEGY::PRIM) {
499 stillToDo.erase(v);
500 for(size_t uToDo : stillToDo) {
501 edgeHeap.push(std::make_pair(
502 distMatrix[v * n + uToDo], std::make_pair(v, uToDo)));
503 }
504 }
505 }
506 std::stringstream ssDistortion;
507 ssDistortion << std::scientific << 2 * finalDistortion;
508 this->printMsg("Non normalized distance matrix distortion: "
509 + ssDistortion.str());
510 this->printMsg(
511 "(for normalized distortion values, use distanceMatrixDistorsion).");
512
513 if(this->errorConvexHull) {
514 this->printWrn(
515 "Warning, somme convex hull computations were wrong, so the projection "
516 "may be less optimal than it could have been. If you are using Boost, "
517 "consider switching to Qhull to avoid this.");
518 }
519 this->printMsg("The minimum spanning tree has weight "
520 + std::to_string(MSTWeight) + ".");
521
522 // We check that the lengths of the edges selected to build a minimum
523 // spanning tree are preserved by our algorithm, as should be.
524 if(CheckMST) {
525
526 edgesMSTAfter.reserve(n);
527 this->printMsg("Checking the new minimum spanning tree.");
528 std::priority_queue<heapType, std::vector<heapType>,
529 std::greater<heapType>>
530 edgeHeapAfter;
531 if(this->Strategy == STRATEGY::KRUSKAL) {
532 for(size_t u1 = 0; u1 < n; u1++) {
533 for(size_t u2 = u1 + 1; u2 < n; u2++) {
534 edgeHeapAfter.push(
535 std::make_pair(ttk::Geometry::distance2D(
536 &outputCoords[2 * u1], &outputCoords[2 * u2]),
537 std::make_pair(u1, u2)));
538 }
539 }
540 } else if(this->Strategy == STRATEGY::PRIM) {
541 for(size_t u1 = 0; u1 < n; u1++) {
542 edgeHeapAfter.push(std::make_pair(
543 ttk::Geometry::distance2D(&outputCoords[0], &outputCoords[2 * u1]),
544 std::make_pair(0, u1)));
545 }
546 }
547
548 std::map<UnionFind *, std::vector<size_t>> ufToCompAfter;
549 std::vector<UnionFind> ufVectorAfter(n);
550 std::vector<UnionFind *> ufPtrVectorAfter(n);
551 for(size_t i = 0; i < n; i++) {
552 ufPtrVectorAfter[i] = &ufVectorAfter[i];
553 ufToCompAfter[ufPtrVectorAfter[i]].push_back(i);
554 }
555
556 nbEdgesMerged = 0;
557 stillToDo.clear();
558 for(size_t i = 1; i < n; i++)
559 stillToDo.insert(i);
560 while(nbEdgesMerged != n - 1) {
561 if(edgeHeapAfter.empty()) {
562 this->printErr("Error building the MST for the check. Aborting.");
563 return 1;
564 }
565 const auto &elt = edgeHeapAfter.top();
566 T edgeCost = elt.first;
567 size_t u = elt.second.first;
568 size_t v = elt.second.second;
569 edgeHeapAfter.pop();
570
571 UnionFind *reprU = ufPtrVectorAfter[u]->find();
572 UnionFind *reprV = ufPtrVectorAfter[v]->find();
573 if(reprU == reprV) { // Already in the same component
574 continue;
575 }
576 UnionFind *unionRepr = UnionFind::makeUnion(reprU, reprV);
577 UnionFind *otherRepr = (unionRepr == reprU) ? reprV : reprU;
578 if(unionRepr != reprU && unionRepr != reprV) {
579 printErr("Error with the union find module results. Aborting.");
580 return 1;
581 }
582 nbEdgesMerged++;
583 std::vector<size_t> &unionSet = ufToCompAfter[unionRepr];
584 std::vector<size_t> &otherSet = ufToCompAfter[otherRepr];
585
586 unionSet.insert(unionSet.end(), otherSet.begin(), otherSet.end());
587
588 ufPtrVectorAfter[u] = unionRepr;
589 ufPtrVectorAfter[v] = unionRepr;
590
591 ufToCompAfter.erase(otherRepr);
592 if(this->Strategy == STRATEGY::PRIM) {
593 stillToDo.erase(v);
594 for(size_t uToDo : stillToDo) {
595 edgeHeapAfter.push(
596 std::make_pair(ttk::Geometry::distance2D(
597 &outputCoords[2 * v], &outputCoords[2 * uToDo]),
598 std::make_pair(v, uToDo)));
599 }
600 }
601
602 edgesMSTAfter.push_back(edgeCost);
603 }
604
605 const T epsilonCheck = Epsilon * 5;
606 double sumDiff = 0, maxDiff = 0;
607 size_t nbProblematic = 0;
608 if(this->Strategy == STRATEGY::PRIM) {
609 sort(edgesMSTBefore.begin(), edgesMSTBefore.end());
610 sort(edgesMSTAfter.begin(), edgesMSTAfter.end());
611 }
612 for(size_t i = 0; i < edgesMSTBefore.size(); i++) {
613 if(fabs(edgesMSTBefore[i] - edgesMSTAfter[i]) >= epsilonCheck) {
614 nbProblematic++;
615 double diff = fabs(edgesMSTBefore[i] - edgesMSTAfter[i]);
616 sumDiff += diff;
617 maxDiff = std::max(diff, maxDiff);
618 }
619 }
620 if(nbProblematic != 0) {
621 this->printWrn("Error on " + std::to_string(nbProblematic)
622 + " edge(s) of the MST.");
623
624 std::stringstream ssMax, ssAvg;
625 ssMax << std::scientific << maxDiff;
626 ssAvg << std::scientific << sumDiff / nbProblematic;
627 this->printWrn("The maximum error is " + ssMax.str()
628 + " and the average (on the number of errors) error is "
629 + ssAvg.str() + ".");
630 }
631 }
632
633 this->printMsg(ttk::debug::Separator::L2); // horizontal '-' separator
634 this->printMsg("Complete", 1, timer.getElapsedTime());
635 this->printMsg(ttk::debug::Separator::L1); // horizontal '=' separator
636
637 return 0;
638 }
639
640 template <typename T>
641 bool TopoMap::getPrevNextEdges(const std::vector<size_t> &idsPtsPolygon,
642 size_t idCenter,
643 const T *allCoords,
644 T *coordPrev,
645 T *coordPost) const {
646 size_t n = idsPtsPolygon.size();
647 size_t iPtPrev = 0, iPtPost = 0;
648 bool found = false;
649
650 for(size_t i = 0; i < n; i++) {
651 if(idsPtsPolygon[i] == idCenter) {
652 found = true;
653 iPtPost = idsPtsPolygon[(i + 1) % n];
654 iPtPrev = idsPtsPolygon[(i + n - 1) % n];
655 break;
656 }
657 }
658 if(!found) {
659 printErr("Error, we could not find the edges incident to the point we "
660 "chose for the rotation of the component. Aborting.");
661 return false;
662 }
663
664 coordPrev[0] = allCoords[2 * iPtPrev];
665 coordPrev[1] = allCoords[2 * iPtPrev + 1];
666 coordPost[0] = allCoords[2 * iPtPost];
667 coordPost[1] = allCoords[2 * iPtPost + 1];
668
669 return true;
670 }
671
672 template <typename T>
673 T TopoMap::rotateMergingCompsBest(const std::vector<size_t> &hull1,
674 const std::vector<size_t> &hull2,
675 const std::vector<size_t> &comp1,
676 const std::vector<size_t> &comp2,
677 size_t iPt1,
678 size_t iPt2,
679 const std::vector<T> &distMatrix,
680 T *allCoords,
681 size_t n,
682 size_t angularSampleNb) {
683 // The distance between the two components.
684 T shortestDistPossible
685 = ttk::Geometry::distance2D(&allCoords[2 * iPt1], &allCoords[2 * iPt2]);
686 T coordPt1[2] = {allCoords[2 * iPt1], allCoords[2 * iPt1 + 1]};
687 T coordPt2[2] = {allCoords[2 * iPt2], allCoords[2 * iPt2 + 1]};
688 size_t comp1Size = comp1.size(), comp2Size = comp2.size();
689
690 T coordPrev1[2], coordPost1[2];
691 T coordPrev2[2], coordPost2[2];
692 getPrevNextEdges(hull1, iPt1, allCoords, coordPrev1, coordPost1);
693 getPrevNextEdges(hull2, iPt2, allCoords, coordPrev2, coordPost2);
694
695 double angle1
696 = ttk::Geometry::angle2DUndirected(coordPt1, coordPrev1, coordPost1);
697 double angle2
698 = ttk::Geometry::angle2DUndirected(coordPt2, coordPrev2, coordPost2);
699 if(angle1 - M_PI > EpsilonDBL || angle2 - M_PI > EpsilonDBL) {
700 this->printErr("One angle of the convex hull is greater than pi. "
701 "Convexity error, aborting.");
702 return -2;
703 }
704 T coordBissect1[2] = {coordPrev1[0], coordPrev1[1]};
705 T coordBissect2[2] = {coordPrev2[0], coordPrev2[1]};
706 rotate(coordBissect1, coordPt1, angle1 / 2);
707 rotate(coordBissect2, coordPt2, angle2 / 2);
708
709 // If a component has only one vertex, no need to try several angles.
710 double semiAngle1 = comp1Size == 1 ? M_PI / 2 : angle1 / 2;
711 double semiAngle2 = comp2Size == 1 ? M_PI / 2 : angle2 / 2;
712
713 double angleMax1 = M_PI / 2 - semiAngle1, angleMin1 = -angleMax1;
714 double angleMax2 = M_PI / 2 - semiAngle2, angleMin2 = -angleMax2;
715 double step1 = (angleMax1 - angleMin1) / angularSampleNb,
716 step2 = (angleMax2 - angleMin2) / angularSampleNb;
717 double bestAnglePair[2] = {0, 0};
718 T bestScore = 3.e38; // Near the maximum value for float
719
720 std::vector<std::vector<T>> origDistMatrix(comp1Size);
721 for(size_t i = 0; i < comp1Size; i++) {
722 origDistMatrix[i].resize(comp2Size);
723 for(size_t j = 0; j < comp2Size; j++) {
724 origDistMatrix[i][j] = distMatrix[comp1[i] * n + comp2[j]];
725 }
726 }
727 std::vector<T> initialCoords1(2 * comp1Size), initialCoords2(2 * comp2Size);
728 for(size_t i = 0; i < comp1.size(); i++) {
729 size_t iPt = comp1[i];
730 initialCoords1[2 * i] = allCoords[2 * iPt];
731 initialCoords1[2 * i + 1] = allCoords[2 * iPt + 1];
732 }
733 for(size_t i = 0; i < comp2.size(); i++) {
734 size_t iPt = comp2[i];
735 initialCoords2[2 * i] = allCoords[2 * iPt];
736 initialCoords2[2 * i + 1] = allCoords[2 * iPt + 1];
737 }
738
739 size_t nbIter1 = std::isfinite(step1) ? angularSampleNb + 1 : 1;
740 size_t nbIter2 = std::isfinite(step2) ? angularSampleNb + 1 : 1;
741
742 if(step1 * nbIter1 < 0.001) { // No need to split such a small angle
743 nbIter1 = 1;
744 }
745 if(step2 * nbIter2 < 0.001) { // No need to split such a small angle
746 nbIter2 = 1;
747 }
748
749 bool errorDuringLoop = false;
750 bool foundValidAngle = false;
751#ifdef TTK_ENABLE_OPENMP
753#endif
754 // We enumerate the rotation angles for the first component.
755 for(size_t i1 = 0; i1 < nbIter1; i1++) {
756 std::vector<T> coords1Test(2 * comp1Size), coords2Test(2 * comp2Size);
757 coords1Test = initialCoords1;
758 double testAngle1 = angleMin1 + step1 * i1;
759 rotatePolygon(coords1Test, coordPt1, testAngle1);
760
761 // For each first rotation angle, we enumerate the rotations angles for
762 // the second componenent.
763 for(size_t i2 = 0; i2 < nbIter2; i2++) {
764 bool curAngleError = false;
765 coords2Test = initialCoords2;
766 double testAngle2 = angleMin2 + i2 * step2;
767 rotatePolygon(coords2Test, coordPt2, testAngle2);
768
769 // We compute a distortion score if we rotate the two components by the
770 // current angles.
771 T curScore = 0;
772 for(size_t i = 0; i < comp1Size; i++) {
773 T coordARotate[2] = {coords1Test[2 * i], coords1Test[2 * i + 1]};
774 for(size_t j = 0; j < comp2Size; j++) {
775 T coordBRotate[2] = {coords2Test[2 * j], coords2Test[2 * j + 1]};
776 T newDist = ttk::Geometry::distance2D(coordARotate, coordBRotate);
777 curScore += (newDist - origDistMatrix[i][j])
778 * (newDist - origDistMatrix[i][j]);
779 if(newDist + Epsilon < shortestDistPossible) {
780 curAngleError = true;
781 errorDuringLoop = true;
782 break;
783 }
784 }
785 }
786
787 if(curAngleError) {
788 continue;
789 }
790
791 // Now if the score is the best so far, we keep this pair angle as the
792 // best one we tried so far.
793#ifdef TTK_ENABLE_OPENMP
794#pragma omp critical
795#endif
796 {
797 // In case of errors in convex hull computation, we still have a valid
798 // angle (to preserve the MST).
799 foundValidAngle = true;
800 // This pair of angle minimises the distortion and the choice is
802 if(curScore < bestScore
803 || (curScore == bestScore
804 && (testAngle1 < bestAnglePair[0]
805 || (testAngle1 == bestAnglePair[0]
806 && testAngle2 < bestAnglePair[1])))) {
807 bestScore = curScore;
808 bestAnglePair[0] = std::isfinite(testAngle1) ? testAngle1 : 0;
809 bestAnglePair[1] = std::isfinite(testAngle2) ? testAngle2 : 0;
810 }
811 }
812 }
813 }
814
815 // No valid angle were found, because the convex hull is invalid;
816 if(!foundValidAngle) {
817 this->printErr(
818 "No valid angle was found, due to errors in the convex hull "
820 return -10;
821 }
822
823 // Applying the chosen rotations to the two components.
824 for(size_t i1 : comp1) {
825 rotate(&allCoords[2 * i1], coordPt1, bestAnglePair[0]);
826 }
827 for(size_t i2 : comp2) {
828 rotate(&allCoords[2 * i2], coordPt2, bestAnglePair[1]);
829 }
830
831 if(errorDuringLoop) {
832 this->errorConvexHull = true;
833 }
834
835 return bestScore;
836 }
837
838 template <typename T>
839 bool TopoMap::computeConvexHull(T *allCoords,
840 const std::vector<size_t> &compPtsIds,
841 std::vector<size_t> &idsInHull) const {
842 size_t nbPoint = compPtsIds.size();
843 std::vector<double> compCoords(2 * nbPoint);
844 for(size_t i = 0; i < nbPoint; i++) {
845 size_t vertId = compPtsIds[i];
846 compCoords[2 * i] = static_cast<double>(allCoords[2 * vertId]);
847 compCoords[2 * i + 1] = static_cast<double>(allCoords[2 * vertId + 1]);
848 }
849
850 if(nbPoint <= 2) {
851 idsInHull.push_back(compPtsIds[0]);
852 if(nbPoint == 2) {
853 double dist = ttk::Geometry::distance2D(&compCoords[0], &compCoords[2]);
854
855 if(dist > Epsilon) {
856 idsInHull.push_back(compPtsIds[1]);
857 }
858 }
859 return true;
860 }
861
862 // Testing if all points are colinear. If so, we take extremal ones. We do
863 // this to avoid that the convex hull algorithm does not detect it and gives
864 // a wrong result.
865 double dirVect[2] = {0, 0};
866 bool areColinear = true;
867
868 size_t idFirstDistinct = 0;
869 while(idFirstDistinct < nbPoint && fabs(dirVect[0]) < Epsilon
870 && fabs(dirVect[1]) < Epsilon) {
871 idFirstDistinct++;
872 dirVect[0] = compCoords[2 * idFirstDistinct] - compCoords[0];
873 dirVect[1] = compCoords[2 * idFirstDistinct + 1] - compCoords[1];
874 if(fabs(dirVect[0]) < Epsilon) {
875 dirVect[0] = 0;
876 }
877 if(fabs(dirVect[1]) < Epsilon) {
878 dirVect[1] = 0;
879 }
880 }
881
882 // Minimum/maximum values for x and y coordinates.
883 int idMins[2] = {compCoords[0] < compCoords[2] ? 0 : 1,
884 compCoords[1] < compCoords[3] ? 0 : 1};
885 int idMaxs[2] = {compCoords[0] > compCoords[2] ? 0 : 1,
886 compCoords[1] > compCoords[3] ? 0 : 1};
887
888 const double *pt0 = &compCoords[0];
889 const double *ptDistinct = &compCoords[2 * idFirstDistinct];
890
891 for(size_t iPt = idFirstDistinct + 1; iPt < nbPoint; iPt++) {
892 double curVect[2] = {compCoords[2 * iPt] - compCoords[0],
893 compCoords[2 * iPt + 1] - compCoords[1]};
894 if(fabs(curVect[0]) < Epsilon && fabs(curVect[1]) < Epsilon) {
895 continue;
896 }
897 const double *ptCur = &compCoords[2 * iPt];
899 pt0, ptDistinct, ptCur, EpsilonDBL)) {
900 areColinear = false;
901 break;
902 }
903
904 // Updating the min and max values for x and y coordinates.
905 if(compCoords[2 * iPt] < compCoords[2 * idMins[0]]) {
906 idMins[0] = iPt;
907 }
908 if(compCoords[2 * iPt + 1] < compCoords[2 * idMins[1]]) {
909 idMins[1] = iPt;
910 }
911 if(compCoords[2 * iPt] > compCoords[2 * idMaxs[0]]) {
912 idMaxs[0] = iPt;
913 }
914 if(compCoords[2 * iPt + 1] > compCoords[2 * idMaxs[1]]) {
915 idMaxs[1] = iPt;
916 }
917 }
918
919 if(areColinear) {
920 if(fabs(dirVect[0]) > Epsilon) {
921 idsInHull.push_back(compPtsIds[idMins[0]]);
922 idsInHull.push_back(compPtsIds[idMaxs[0]]);
923 } else {
924 idsInHull.push_back(compPtsIds[idMins[1]]);
925 idsInHull.push_back(compPtsIds[idMaxs[1]]);
926 }
927 return true;
928 }
929
930 // Computing the convex hull using Boost or Qhull according to the build.
931 std::string errMsg;
932 bool status = computeConvexHull_aux(compCoords, idsInHull, errMsg);
933 if(!status) {
934 this->printErr(errMsg);
935 return false;
936 }
937
938 // In its C++ API, Qhull does not return the points of the convex hull in
939 // (anti-)clockwise order. Boost does not seem consistent with the
940 // orientation of the hull. So we reorder the points by sorting them
941 // according to the angle they form with the center and a fictitious point
942 // at its right.
943 double sumX = 0, sumY = 0;
944 for(size_t u : idsInHull) {
945 sumX += compCoords[2 * u];
946 sumY += compCoords[2 * u + 1];
947 }
948 double bary[2] = {sumX / idsInHull.size(), sumY / idsInHull.size()};
949 double baryRight[2] = {bary[0] + 2, bary[1]};
950 std::vector<std::pair<double, size_t>> ptsToSort;
951 ptsToSort.reserve(idsInHull.size());
952 for(size_t u : idsInHull) {
953 const double curPt[2] = {compCoords[2 * u], compCoords[2 * u + 1]};
954 double curAngle
955 = ttk::Geometry::angle2DUndirected(bary, baryRight, curPt);
956 ptsToSort.emplace_back(std::make_pair(-curAngle, u));
957 }
958
959 // We sort the points according to the angle, keeping the indices
960 // associated.
961 sort(ptsToSort.begin(), ptsToSort.end());
962 for(size_t i = 0; i < ptsToSort.size(); i++)
963 idsInHull[i] = ptsToSort[i].second;
964
965 // We obtained the indices of the points belonging to the hull. We replace
966 // them by the global ids of the chosen points.
967 for(size_t &x : idsInHull) {
968 x = compPtsIds[x];
969 }
970 return true;
971 }
972} // namespace ttk
973
974template <typename T>
975static void
976 computeUnitVector(const T *coordOrig, const T *coordDest, T *coordVect) {
977 T tmp[2] = {coordDest[0] - coordOrig[0], coordDest[1] - coordOrig[1]};
978 T dist = sqrt(tmp[0] * tmp[0] + tmp[1] * tmp[1]);
979 if(dist < Epsilon) {
980 coordVect[0] = 1;
981 coordVect[1] = 0;
982 } else {
983 coordVect[0] = tmp[0] / dist;
984 coordVect[1] = tmp[1] / dist;
985 }
986}
987
988template <typename T>
989static void rotate(T *ptToRotate, const T *center, double angle) {
990 const double &xCtr = center[0], &yCtr = center[1];
991 T &xPt = ptToRotate[0], &yPt = ptToRotate[1];
992 const double dx = xPt - xCtr, dy = yPt - yCtr;
993 xPt = dx * cos(angle) - dy * sin(angle) + xCtr;
994 yPt = dx * sin(angle) + dy * cos(angle) + yCtr;
995}
996
997template <typename T>
998static void
999 rotatePolygon(std::vector<T> &coords, T *centerCoords, const double angle) {
1000 double xCenter = centerCoords[0], yCenter = centerCoords[1];
1001 size_t nbPoint = coords.size() / 2;
1002 for(size_t iPt = 0; iPt < nbPoint; iPt++) {
1003 T &x = coords[iPt * 2], &y = coords[iPt * 2 + 1];
1004 double xNew, yNew;
1005 xNew = (x - xCenter) * cos(angle) - (y - yCenter) * sin(angle) + xCenter;
1006 yNew = (y - yCenter) * cos(angle) + (x - xCenter) * sin(angle) + yCenter;
1007 x = xNew;
1008 y = yNew;
1009 }
1010}
#define M_PI
Definition Os.h:50
bool computeConvexHull_aux(const std::vector< double > &coords, std::vector< size_t > &res, std::string &errMsg)
Definition TopoMap.cpp:60
bool computeConvexHull_aux(const std::vector< double > &coords, std::vector< size_t > &res, std::string &errMsg)
Definition TopoMap.cpp:60
Minimalist debugging class.
Definition Debug.h:88
void setDebugMsgPrefix(const std::string &prefix)
Definition Debug.h:364
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition Debug.h:149
double getElapsedTime()
Definition Timer.h:15
STRATEGY Strategy
Definition TopoMap.h:118
int execute(T *outputCoords, int *insertionTime, const std::vector< T > &inputMatrix, bool isDistMat, size_t n)
Computes the TopoMap projection.
Definition TopoMap.h:183
bool errorConvexHull
Definition TopoMap.h:117
bool CheckMST
Definition TopoMap.h:116
size_t AngularSampleNb
Definition TopoMap.h:115
TopoMap(size_t angularSampleNb, bool checkMST, STRATEGY strategy)
Definition TopoMap.h:78
~TopoMap() override
Union Find implementation for connectivity tracking.
Definition UnionFind.h:17
static UnionFind * makeUnion(UnionFind *uf0, UnionFind *uf1)
Definition UnionFind.h:40
UnionFind * find()
Definition UnionFind.h:99
std::string to_string(__int128)
Definition ripserpy.cpp:99
double angle2DUndirected(const T *vA, const T *vB, const T *vC)
Definition Geometry.h:31
bool isTriangleColinear2D(const T *pptA, const T *pptB, const T *pptC, const T tolerance)
Definition Geometry.cpp:130
T distance2D(const T *p0, const T *p1)
Definition Geometry.h:198
T powInt(const T val, const int n)
Definition Geometry.h:381
The Topology ToolKit.
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)