Merge pull request #987 from atalon-lip6/topomap-clean

TopoMap

.github/workflows/check.yml

@@ -113,6 +113,7 @@ jobs:
           graphviz \
           python3-sklearn \
           zlib1g-dev \
+          libqhull-dev \
           doxygen \
           clang-tidy
 
@@ -197,6 +198,7 @@ jobs:
           graphviz \
           python3-sklearn \
           zlib1g-dev \
+          libqhull-dev \
           clang-tools \
           libomp-dev
 

.github/workflows/package.yml

@@ -43,6 +43,7 @@ jobs:
           graphviz \
           python3-sklearn \
           zlib1g-dev \
+          libqhull-dev \
           dpkg-dev
 
     - name: Install optional dependencies

.github/workflows/test.yml

@@ -64,6 +64,7 @@ jobs:
           graphviz \
           ninja-build \
           zlib1g-dev \
+          libqhull-dev \
           dpkg-dev
         sudo python3 -m pip install scikit-learn
 
@@ -281,7 +282,7 @@ jobs:
         brew install --cask xquartz
         brew install wget llvm libomp ninja open-mpi
         # TTK dependencies
-        brew install boost eigen graphviz numpy websocketpp ccache
+        brew install boost eigen graphviz numpy qhull websocketpp ccache
         # prepend ccache to system path
         echo "$(brew --prefix)/opt/ccache/libexec" >> $GITHUB_PATH
 
@@ -428,7 +429,7 @@ jobs:
       shell: bash
       run: |
         conda install -c conda-forge boost glew eigen spectralib zfp \
-          scikit-learn openmp graphviz ninja websocketpp sccache=0.4.2 python=3.10 zlib
+          scikit-learn openmp graphviz ninja websocketpp sccache=0.4.2 python=3.10 zlib qhull
         # add sccache to PATH
         echo "$CONDA_ROOT/bin" >> $GITHUB_PATH
         # add TTK & ParaView install folders to PATH
@@ -465,7 +466,7 @@ jobs:
     - name: Create & configure TTK build directory
       shell: cmd
       run: |
-        set CMAKE_PREFIX_PATH=%CONDA_ROOT%\Library\lib\cmake;%CONDA_ROOT%\Library\share\eigen3\cmake;%CONDA_ROOT%\Library\cmake;%ProgramFiles%\TTK-ParaView\lib\cmake
+        set CMAKE_PREFIX_PATH=%CONDA_ROOT%\Library\lib\cmake;%CONDA_ROOT%\Library\share\eigen3\cmake;%CONDA_ROOT%\Library\share\Qull\cmake;%CONDA_ROOT%\Library\cmake;%ProgramFiles%\TTK-ParaView\lib\cmake
         set CC=clang-cl.exe
         set CXX=clang-cl.exe
         call "%ProgramFiles%\Microsoft Visual Studio\2022\Enterprise\VC\Auxiliary\Build\vcvars64.bat"

CMake/Print.cmake

@@ -10,6 +10,7 @@ function(ttk_print_summary)
     message(STATUS "TTK_ENABLE_KAMIKAZE: ${TTK_ENABLE_KAMIKAZE}")
     message(STATUS "TTK_ENABLE_MPI: ${TTK_ENABLE_MPI}")
     message(STATUS "TTK_ENABLE_OPENMP: ${TTK_ENABLE_OPENMP}")
+    message(STATUS "TTK_ENABLE_QHULL: ${TTK_ENABLE_QHULL}")
     message(STATUS "TTK_ENABLE_SCIKIT_LEARN: ${TTK_ENABLE_SCIKIT_LEARN}")
     message(STATUS "TTK_ENABLE_SHARED_BASE_LIBRARIES: ${TTK_ENABLE_SHARED_BASE_LIBRARIES}")
     message(STATUS "TTK_ENABLE_SPECTRA: ${TTK_ENABLE_SPECTRA}")

CMake/config.cmake

@@ -295,6 +295,17 @@ else()
   message(STATUS "WebSocketPP not found, disabling WebSocketIO module in TTK.")
 endif()
 
+
+option(TTK_ENABLE_QHULL "Use Qhull instead of Boost for convex hulls" ON)
+if (TTK_ENABLE_QHULL)
+  find_package(Qhull QUIET)
+  if(Qhull_FOUND)
+    message(STATUS "Found Qhull ${Qhull_VERSION} (${Qhull_DIR})")
+  else()
+    message(STATUS "Qhull not found, disabling Qhull support in TTK.")
+  endif()
+endif()
+
 # --- Install path
 
 include(GNUInstallDirs)

core/base/TTKBaseConfig.cmake.in

@@ -13,6 +13,10 @@ if (@TTK_ENABLE_OPENMP@)
   find_dependency(OpenMP REQUIRED)
 endif()
 
+if (@TTK_ENABLE_QHULL@ AND @Qhull_FOUND@)
+  find_dependency(Qhull REQUIRED)
+endif()
+
 if (@TTK_ENABLE_ZFP@ AND NOT @TTK_ENABLE_SHARED_BASE_LIBRARIES@)
   set(ZFP_DIR "@ZFP_DIR@" CACHE PATH "Use TTK ZFP dir" FORCE)
   find_dependency(ZFP REQUIRED @ZFP_DIR@)

core/base/dimensionReduction/CMakeLists.txt

@@ -5,6 +5,7 @@ ttk_add_base_library(dimensionReduction
     DimensionReduction.h
   DEPENDS
     triangulation
+    topoMap
   )
 
 install(

core/base/dimensionReduction/DimensionReduction.cpp

@@ -1,4 +1,5 @@
 #include <DimensionReduction.h>
+#include <TopoMap.h>
 
 #include <map>
 
@@ -51,6 +52,7 @@ bool DimensionReduction::isPythonFound() const {
 
 int DimensionReduction::execute(
   std::vector<std::vector<double>> &outputEmbedding,
+  int *insertionTimeForTopomap,
   const std::vector<double> &inputMatrix,
   const int nRows,
   const int nColumns) const {
@@ -70,6 +72,28 @@ int DimensionReduction::execute(
   Timer t;
 
   const int numberOfComponents = std::max(2, this->NumberOfComponents);
+  if(this->Method == METHOD::TOPOMAP) {
+    TopoMap topomap(
+      this->topomap_AngularSampleNb, topomap_CheckMST, topomap_Strategy);
+    topomap.setDebugLevel(this->debugLevel_);
+    topomap.setThreadNumber(this->threadNumber_);
+
+    std::vector<double> coordsTopomap(2 * nRows);
+    topomap.execute<double>(coordsTopomap.data(), insertionTimeForTopomap,
+                            inputMatrix, IsInputADistanceMatrix, nRows);
+    outputEmbedding.resize(2);
+    outputEmbedding[0].resize(nRows);
+    outputEmbedding[1].resize(nRows);
+    for(int i = 0; i < nRows; i++) {
+      outputEmbedding[0][i] = coordsTopomap[2 * i];
+      outputEmbedding[1][i] = coordsTopomap[2 * i + 1];
+    }
+
+    this->printMsg(
+      "Computed TopoMap", 1.0, t.getElapsedTime(), this->threadNumber_);
+    return 0;
+  }
+
   const int numberOfNeighbors = std::max(1, this->NumberOfNeighbors);
 
   // declared here to avoid crossing initialization with goto

core/base/dimensionReduction/DimensionReduction.h

@@ -38,9 +38,16 @@
 ///   Generators Periodic Picture example</a> \n
 ///
 
+/// \b Related \b publication: \n
+/// "Topomap: A 0-dimensional homology preserving projection of high-dimensional
+/// data"\n Harish Doraiswamy, Julien Tierny, Paulo J. S. Silva, Luis Gustavo
+/// Nonato, and Claudio Silva\n Proc. of IEEE VIS 2020.\n IEEE Transactions on
+/// Visualization and Computer Graphics 27(2): 561-571, 2020.
+
 #pragma once
 
 #include <Debug.h>
+#include <TopoMap.h>
 
 namespace ttk {
 
@@ -63,6 +70,8 @@ namespace ttk {
       ISOMAP = 4,
       /** Principal Component Analysis */
       PCA = 5,
+      /** TopoMap */
+      TOPOMAP = 6,
     };
 
     inline void setSEParameters(const std::string &Affinity,
@@ -161,6 +170,10 @@ namespace ttk {
       pca_Tolerance = Tolerance;
       pca_MaxIteration = MaxIteration;
     }
+    inline void setTopoParameters(const size_t AngularSampleNb, bool CheckMST) {
+      topomap_AngularSampleNb = AngularSampleNb;
+      topomap_CheckMST = CheckMST;
+    }
 
     inline void setInputModulePath(const std::string &modulePath) {
       ModulePath = modulePath;
@@ -191,6 +204,7 @@ namespace ttk {
     }
 
     inline void setIsInputDistanceMatrix(const bool data) {
+      this->IsInputADistanceMatrix = data;
       if(data) {
         this->se_Affinity = "precomputed";
         this->mds_Dissimilarity = "precomputed";
@@ -207,6 +221,7 @@ namespace ttk {
     bool isPythonFound() const;
 
     int execute(std::vector<std::vector<double>> &outputEmbedding,
+                int *insertionTimeForTopoMap,
                 const std::vector<double> &inputMatrix,
                 const int nRows,
                 const int nColumns) const;
@@ -263,6 +278,11 @@ namespace ttk {
     float pca_Tolerance{0};
     std::string pca_MaxIteration{"auto"};
 
+    // TopoMap
+    size_t topomap_AngularSampleNb;
+    bool topomap_CheckMST;
+    TopoMap::STRATEGY topomap_Strategy{TopoMap::STRATEGY::KRUSKAL};
+
     // testing
     std::string ModulePath{"default"};
     std::string ModuleName{"dimensionReduction"};
@@ -273,5 +293,6 @@ namespace ttk {
     int NumberOfNeighbors{5};
     int IsDeterministic{true};
     char majorVersion_{'0'};
+    bool IsInputADistanceMatrix{false};
   };
 } // namespace ttk

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 &);                                                              \

core/base/geometry/Geometry.h

@@ -9,6 +9,7 @@
 
 #include <Debug.h>
 #include <array>
+#include <cmath>
 
 namespace ttk {
 
@@ -22,6 +23,19 @@ 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 +53,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.
@@ -175,6 +194,30 @@ namespace ttk {
     template <typename T>
     T distance(const std::vector<T> &p0, const std::vector<T> &p1);
 
+    template <typename T>
+    inline T distance2D(const T *p0, const T *p1) {
+      T distance = 0;
+      T di0 = (p0[0] - p1[0]);
+      T di1 = (p0[1] - p1[1]);
+      if(std::is_same<T, float>::value) { // TODO constexpr when c++17
+        distance = fmaf(di0, di0, distance);
+        distance = fmaf(di1, di1, distance);
+      } else {
+        distance = fma(di0, di0, distance);
+        distance = fma(di1, di1, distance);
+      }
+
+      return sqrt(distance);
+    }
+
+    /// Compute the Euclidean distance between two vectors by first flattening
+    /// them
+    /// \param p0 xyz coordinates of the first input point.
+    /// \param p1 xyz coordinates of the second input point.
+    template <typename T>
+    T distanceFlatten(const std::vector<std::vector<T>> &p0,
+                      const std::vector<std::vector<T>> &p1);
+
     /// Compute the Euclidean distance between two vectors by first flattening
     /// them
     /// \param p0 xyz coordinates of the first input point.

core/base/topoMap/CMakeLists.txt

@@ -0,0 +1,20 @@
+ttk_add_base_library(topoMap
+  SOURCES
+    TopoMap.cpp
+  HEADERS
+    TopoMap.h
+  DEPENDS
+    geometry
+    unionFind
+)
+
+if(TTK_ENABLE_QHULL)
+  if (Qhull_FOUND)
+    target_link_libraries(topoMap PUBLIC Qhull::qhullcpp)
+    target_link_libraries(topoMap PUBLIC Qhull::qhull_r)
+    target_compile_definitions(topoMap PUBLIC Qhull_FOUND)
+  endif()
+  target_compile_definitions(topoMap PUBLIC TTK_ENABLE_QHULL)
+endif()
+
+