Some comments addressed.

topology-tool-kit · Nov 28, 2023 · 1d6d2cc · 1d6d2cc
1 parent 12851ce
commit 1d6d2cc
Show file tree

Hide file tree

Showing 5 changed files with 78 additions and 55 deletions.
diff --git a/core/base/dimensionReduction/DimensionReduction.cpp b/core/base/dimensionReduction/DimensionReduction.cpp
@@ -78,7 +78,7 @@ int DimensionReduction::execute(
     topomap.setDebugLevel(this->debugLevel_);
     topomap.setThreadNumber(this->threadNumber_);
 
-    std::vector<double> coordsTopomap(2*nRows);
+    std::vector<double> coordsTopomap(2 * nRows);
     topomap.execute<double>(coordsTopomap.data(), insertionTimeForTopomap,
                             inputMatrix, IsInputADistanceMatrix, nRows);
     outputEmbedding.resize(2);
@@ -89,8 +89,8 @@ int DimensionReduction::execute(
       outputEmbedding[1][i] = coordsTopomap[2 * i + 1];
     }
 
-    this->printMsg("Computed Topological Mapper", 1.0, t.getElapsedTime(),
-                   this->threadNumber_);
+    this->printMsg(
+      "Computed TopoMap", 1.0, t.getElapsedTime(), this->threadNumber_);
     return 0;
   }
 

diff --git a/core/base/dimensionReduction/DimensionReduction.h b/core/base/dimensionReduction/DimensionReduction.h
@@ -278,7 +278,7 @@ namespace ttk {
     float pca_Tolerance{0};
     std::string pca_MaxIteration{"auto"};
 
-    // Topological Mapper
+    // TopoMap
     size_t topomap_AngularSampleNb;
     bool topomap_CheckMST;
     TopoMap::STRATEGY topomap_Strategy{TopoMap::STRATEGY::KRUSKAL};

diff --git a/core/base/geometry/Geometry.cpp b/core/base/geometry/Geometry.cpp
@@ -16,6 +16,13 @@ T Geometry::angle(const T *vA0, const T *vA1, const T *vB0, const T *vB1) {
                 / (magnitude(vA0, vA1) * magnitude(vB0, vB1)));
 }
 
+template <typename T>
+T Geometry::angle2D(const T *vA0, const T *vA1, const T *vB0, const T *vB1) {
+  double vecA[2] = {vA1[0] - vA0[0], vA1[1] - vA0[1]};
+  double vecB[2] = {vB1[0] - vB0[0], vB1[1] - vB0[1]};
+  return atan2(vecB[1], vecB[0]) - atan2(vecA[1], vecA[0]);
+}
+
 template <typename T>
 bool Geometry::areVectorsColinear(const T *vA0,
                                   const T *vA1,
@@ -114,12 +121,23 @@ bool Geometry::areVectorsColinear(const T *vA0,
       return true;
     }
   }
-
   k[0] = k[1] = k[2] = 0;
 
   return false;
 }
 
+template <typename T>
+bool Geometry::isTriangleColinear2D(const T *pptA,
+                                    const T *pptB,
+                                    const T *pptC,
+                                    const T tolerance) {
+  double ptA[2] = {pptA[0], pptA[1]}, ptB[2] = {pptB[0], pptB[1]},
+         ptC[2] = {pptC[0], pptC[1]};
+  return fabs(ptA[0] * (ptB[1] - ptC[1]) + ptB[0] * (ptC[1] - ptA[1])
+              + ptC[0] * (ptA[1] - ptB[1]))
+         <= tolerance;
+}
+
 template <typename T>
 int Geometry::computeBarycentricCoordinates(const T *p0,
                                             const T *p1,
@@ -754,9 +772,13 @@ void Geometry::transposeMatrix(const std::vector<std::vector<T>> &a,
 #define GEOMETRY_SPECIALIZE(TYPE)                                              \
   template TYPE Geometry::angle<TYPE>(                                         \
     TYPE const *, TYPE const *, TYPE const *, TYPE const *);                   \
+  template TYPE Geometry::angle2D<TYPE>(                                       \
+    TYPE const *, TYPE const *, TYPE const *, TYPE const *);                   \
   template bool Geometry::areVectorsColinear<TYPE>(                            \
     TYPE const *, TYPE const *, TYPE const *, TYPE const *,                    \
     std::array<TYPE, 3> *, TYPE const *);                                      \
+  template bool Geometry::isTriangleColinear2D<TYPE>(                          \
+    TYPE const *, TYPE const *, TYPE const *, TYPE const);                     \
   template int Geometry::computeBarycentricCoordinates<TYPE>(                  \
     TYPE const *, TYPE const *, TYPE const *, std::array<TYPE, 2> &,           \
     int const &);                                                              \

diff --git a/core/base/geometry/Geometry.h b/core/base/geometry/Geometry.h
@@ -9,6 +9,7 @@
 
 #include <Debug.h>
 #include <array>
+#include <cmath>
 
 namespace ttk {
 
@@ -22,6 +23,20 @@ namespace ttk {
     template <typename T>
     T angle(const T *vA0, const T *vA1, const T *vB0, const T *vB1);
 
+    template <typename T>
+    T angle2D(const T *vA0, const T *vA1, const T *vB0, const T *vB1);
+
+    // Computes the BAC angle, in the interval [0, 2pi]
+    template <typename T>
+    inline double angle2DUndirected(const T *vA, const T *vB, const T *vC)
+    {
+      double angle = angle2D<T>(vA, vB, vA, vC);
+      if (angle < 0) {
+        angle += 2 * M_PI;
+      }
+      return angle;
+    }
+
     /// Check if two 3D std::vectors vA and vB are colinear.
     /// \param vA0 xyz coordinates of vA's origin
     /// \param vA1 xyz coordinates of vA's destination
@@ -39,6 +54,11 @@ namespace ttk {
                             const T *vB1,
                             std::array<T, 3> *coefficients = nullptr,
                             const T *tolerance = NULL);
+    template <typename T>
+    bool isTriangleColinear2D(const T *pptA,
+                              const T *pptB,
+                              const T *pptC,
+                              const T tolerance);
 
     /// Compute the barycentric coordinates of point \p p with regard to the
     /// edge defined by the 3D points \p p0 and \p p1.

diff --git a/core/base/topoMap/TopoMap.h b/core/base/topoMap/TopoMap.h
@@ -40,24 +40,11 @@
 static double Epsilon{ttk::Geometry::pow(10.0, -FLT_DIG + 2)};
 static double EpsilonDBL{ttk::Geometry::pow(10.0, -DBL_DIG + 2)};
 
-template <typename T>
-inline bool are_colinear(const T *pptA, const T *pptB, const T *pptC) {
-  double ptA[2] = {pptA[0], pptA[1]}, ptB[2] = {pptB[0], pptB[1]},
-         ptC[2] = {pptC[0], pptC[1]};
-  return fabs(ptA[0] * (ptB[1] - ptC[1]) + ptB[0] * (ptC[1] - ptA[1])
-              + ptC[0] * (ptA[1] - ptB[1]))
-         <= EpsilonDBL;
-}
-
 // Normalizes a given vector.
 template <typename T>
 static void
   computeUnitVector(const T *coordOrig, const T *coordDest, T *coordVect);
 
-// Computes the oriented angle BAC
-template <typename T>
-static double computeAngle(const T *ptA, const T *ptB, const T *ptC);
-
 // Rotates the first argument by angle around center.
 template <typename T>
 static void rotate(T *ptToRotate, const T *center, double angle);
@@ -99,10 +86,10 @@ namespace ttk {
     ~TopoMap() override;
 
     /**
-     * @brief Computes the topological mapper
+     * @brief Computes the TopoMap projection
      *
      * @param[out] outputCoords the final coordinates of the points, computed by
-     the topological mapper
+     TopoMap
      *
      * @param[out] insertionTime an optional array indicating for each point at
      what time it was inserted in the projection.
@@ -231,29 +218,31 @@ namespace ttk {
       printErr("Error, template type must be either double or float.\n");
       return 1;
     }
-#ifndef TTK_KAMIKAZE
+
     if(isDistMat) {
       this->printMsg("Input data: " + std::to_string(n) + " points.");
     } else {
       this->printMsg("Input data: " + std::to_string(n) + " points ("
                      + std::to_string(inputMatrix.size() / n)
                      + " dimensions).");
     }
-#endif
+
     std::vector<T> computedDistMatrix(n * n);
     if(!isDistMat) {
       size_t dim = inputMatrix.size() / n;
       if(n * dim != inputMatrix.size()) {
         this->printErr("Error, the coordinates input matrix has invalid size.");
         return 1;
       }
+
+#ifdef TTK_ENABLE_OPENMP
+#pragma omp parallel for num_threads(this->threadNumber_) \
+  shared(computedDistMatrix)
+#endif
       for(size_t i1 = 0; i1 < n; i1++) {
         for(size_t i2 = i1; i2 < n; i2++) {
-          T dist2 = 0;
-          for(size_t k = 0; k < dim; k++)
-            dist2 += (inputMatrix[dim * i1 + k] - inputMatrix[dim * i2 + k])
-                     * (inputMatrix[dim * i1 + k] - inputMatrix[dim * i2 + k]);
-          T dist = sqrt(dist2);
+          T dist = ttk::Geometry::distance<T>(
+            &inputMatrix[dim * i1], &inputMatrix[dim * i2], dim);
           computedDistMatrix[i1 * n + i2] = dist;
           computedDistMatrix[i2 * n + i1] = dist;
         }
@@ -420,8 +409,8 @@ namespace ttk {
 
       // Computing the angles.
       double angleSmall
-        = computeAngle(coordPtSmall, coordPrevSmall, coordPostSmall);
-      double angleBig = computeAngle(coordPtBig, coordPrevBig, coordPostBig);
+        = ttk::Geometry::angle2DUndirected(coordPtSmall, coordPrevSmall, coordPostSmall);
+      double angleBig = ttk::Geometry::angle2DUndirected(coordPtBig, coordPrevBig, coordPostBig);
       if(angleSmall - M_PI > EpsilonDBL || angleBig - M_PI > EpsilonDBL) {
         this->printErr("Error, angle out of bound (greater than pi).");
       }
@@ -450,8 +439,8 @@ namespace ttk {
         = {goalCoordChosenSmall[0] - outputCoords[idChosenSmall * 2],
            goalCoordChosenSmall[1] - outputCoords[idChosenSmall * 2 + 1]};
 
-      T distBaryPointSmall
-        = ttk::Geometry::distance(&outputCoords[idChosenSmall * 2], coordscenterSmall, 2);
+      T distBaryPointSmall = ttk::Geometry::distance(
+        &outputCoords[idChosenSmall * 2], coordscenterSmall, 2);
       T preFinalPosBarySmall[2] = {coordscenterSmall[0] + smallCompMoveVect[0],
                                    coordscenterSmall[1] + smallCompMoveVect[1]};
       T finalPosBarySmall[2]
@@ -466,7 +455,7 @@ namespace ttk {
         outputCoords[curIdSmall * 2 + 1] += smallCompMoveVect[1];
 
         // 2.e Performing the rotation.
-        double rotationAngle = computeAngle(
+        double rotationAngle = ttk::Geometry::angle2DUndirected(
           goalCoordChosenSmall, preFinalPosBarySmall, finalPosBarySmall);
         if(nSmall > 1 && std::isfinite(rotationAngle)) {
           rotate(
@@ -520,7 +509,7 @@ namespace ttk {
     this->printMsg("Non normalized distance matrix distortion: "
                    + ssDistortion.str());
     this->printMsg(
-      "(for normalized distortion values, use DistanceMatrixDistorsion).");
+      "(for normalized distortion values, use distanceMatrixDistorsion).");
 
     if(this->errorConvexHull) {
       this->printWrn(
@@ -543,9 +532,10 @@ namespace ttk {
       if(this->Strategy == STRATEGY::KRUSKAL) {
         for(size_t u1 = 0; u1 < n; u1++) {
           for(size_t u2 = u1 + 1; u2 < n; u2++) {
-            edgeHeapAfter.push(std::make_pair(
-              ttk::Geometry::distance(&outputCoords[2 * u1], &outputCoords[2 * u2], 2),
-              std::make_pair(u1, u2)));
+            edgeHeapAfter.push(
+              std::make_pair(ttk::Geometry::distance(
+                               &outputCoords[2 * u1], &outputCoords[2 * u2], 2),
+                             std::make_pair(u1, u2)));
           }
         }
       } else if(this->Strategy == STRATEGY::PRIM) {
@@ -604,7 +594,8 @@ namespace ttk {
           stillToDo.erase(v);
           for(size_t uToDo : stillToDo) {
             edgeHeapAfter.push(std::make_pair(
-              ttk::Geometry::distance(&outputCoords[2 * v], &outputCoords[2 * uToDo], 2),
+              ttk::Geometry::distance(
+                &outputCoords[2 * v], &outputCoords[2 * uToDo], 2),
               std::make_pair(v, uToDo)));
           }
         }
@@ -690,7 +681,7 @@ namespace ttk {
                                     T *allCoords,
                                     size_t n,
                                     size_t angularSampleNb,
-                                    size_t nThread) {
+                                    size_t nThread) { // TODO remove nThread
     TTK_FORCE_USE(nThread);
     // The distance between the two components.
     T shortestDistPossible
@@ -704,8 +695,8 @@ namespace ttk {
     getPrevNextEdges(hull1, iPt1, allCoords, coordPrev1, coordPost1);
     getPrevNextEdges(hull2, iPt2, allCoords, coordPrev2, coordPost2);
 
-    double angle1 = computeAngle(coordPt1, coordPrev1, coordPost1);
-    double angle2 = computeAngle(coordPt2, coordPrev2, coordPost2);
+    double angle1 = ttk::Geometry::angle2DUndirected(coordPt1, coordPrev1, coordPost1);
+    double angle2 = ttk::Geometry::angle2DUndirected(coordPt2, coordPrev2, coordPost2);
     if(angle1 - M_PI > EpsilonDBL || angle2 - M_PI > EpsilonDBL) {
       this->printErr("One angle of the convex hull is greater than pi. "
                      "Convexity error, aborting.");
@@ -860,7 +851,8 @@ namespace ttk {
     if(nbPoint <= 2) {
       idsInHull.push_back(compPtsIds[0]);
       if(nbPoint == 2) {
-        double dist = ttk::Geometry::distance(&compCoords[0], &compCoords[2], 2);
+        double dist
+          = ttk::Geometry::distance(&compCoords[0], &compCoords[2], 2);
 
         if(dist > Epsilon) {
           idsInHull.push_back(compPtsIds[1]);
@@ -905,7 +897,8 @@ namespace ttk {
         continue;
       }
       const double *ptCur = &compCoords[2 * iPt];
-      if(!are_colinear(pt0, ptDistinct, ptCur)) {
+      if(!ttk::Geometry::isTriangleColinear2D(
+           pt0, ptDistinct, ptCur, EpsilonDBL)) {
         areColinear = false;
         break;
       }
@@ -960,7 +953,7 @@ namespace ttk {
     ptsToSort.reserve(idsInHull.size());
     for(size_t u : idsInHull) {
       const double curPt[2] = {compCoords[2 * u], compCoords[2 * u + 1]};
-      double curAngle = computeAngle(bary, baryRight, curPt);
+      double curAngle = ttk::Geometry::angle2DUndirected(bary, baryRight, curPt);
       ptsToSort.emplace_back(std::make_pair(-curAngle, u));
     }
 
@@ -993,18 +986,6 @@ static void
   }
 }
 
-template <typename T>
-static double computeAngle(const T *ptA, const T *ptB, const T *ptC) {
-  double angle;
-  double vect1[2] = {ptB[0] - ptA[0], ptB[1] - ptA[1]};
-  double vect2[2] = {ptC[0] - ptA[0], ptC[1] - ptA[1]};
-
-  angle = atan2(vect2[1], vect2[0]) - atan2(vect1[1], vect1[0]);
-  if(angle < 0) {
-    angle += 2 * M_PI;
-  }
-  return angle;
-}
 
 template <typename T>
 static void rotate(T *ptToRotate, const T *center, double angle) {