diff --git a/.github/workflows/check.yml b/.github/workflows/check.yml index b961b1189f..696f876531 100644 --- a/.github/workflows/check.yml +++ b/.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 diff --git a/.github/workflows/package.yml b/.github/workflows/package.yml index aa538b98cb..d8f74cda5c 100644 --- a/.github/workflows/package.yml +++ b/.github/workflows/package.yml @@ -43,6 +43,7 @@ jobs: graphviz \ python3-sklearn \ zlib1g-dev \ + libqhull-dev \ dpkg-dev - name: Install optional dependencies diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 36bbfea636..a13a578124 100644 --- a/.github/workflows/test.yml +++ b/.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 @@ -276,7 +277,7 @@ jobs: brew install --cask xquartz brew install wget llvm libomp ninja python 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 @@ -423,7 +424,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 @@ -460,7 +461,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" diff --git a/CMake/Print.cmake b/CMake/Print.cmake index 95f4c2bb92..4416894bc9 100644 --- a/CMake/Print.cmake +++ b/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}") diff --git a/CMake/config.cmake b/CMake/config.cmake index 46e9f4e8d9..19d0ad9687 100644 --- a/CMake/config.cmake +++ b/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) diff --git a/core/base/TTKBaseConfig.cmake.in b/core/base/TTKBaseConfig.cmake.in index de8f634273..e8719cfa18 100644 --- a/core/base/TTKBaseConfig.cmake.in +++ b/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@) diff --git a/core/base/dimensionReduction/CMakeLists.txt b/core/base/dimensionReduction/CMakeLists.txt index e30372419c..d6aee3db22 100644 --- a/core/base/dimensionReduction/CMakeLists.txt +++ b/core/base/dimensionReduction/CMakeLists.txt @@ -5,6 +5,7 @@ ttk_add_base_library(dimensionReduction DimensionReduction.h DEPENDS triangulation + topoMap ) install( diff --git a/core/base/dimensionReduction/DimensionReduction.cpp b/core/base/dimensionReduction/DimensionReduction.cpp index 6528f1a99e..08f530292d 100644 --- a/core/base/dimensionReduction/DimensionReduction.cpp +++ b/core/base/dimensionReduction/DimensionReduction.cpp @@ -1,4 +1,5 @@ #include +#include #include @@ -51,6 +52,7 @@ bool DimensionReduction::isPythonFound() const { int DimensionReduction::execute( std::vector> &outputEmbedding, + int *insertionTimeForTopomap, const std::vector &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 coordsTopomap(2 * nRows); + topomap.execute(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 diff --git a/core/base/dimensionReduction/DimensionReduction.h b/core/base/dimensionReduction/DimensionReduction.h index fb1c27e986..dca11f65f2 100644 --- a/core/base/dimensionReduction/DimensionReduction.h +++ b/core/base/dimensionReduction/DimensionReduction.h @@ -38,9 +38,16 @@ /// Generators Periodic Picture example \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 +#include 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> &outputEmbedding, + int *insertionTimeForTopoMap, const std::vector &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 diff --git a/core/base/geometry/Geometry.cpp b/core/base/geometry/Geometry.cpp index 7c0b388be3..2e5eab851c 100644 --- 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 +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 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 +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 int Geometry::computeBarycentricCoordinates(const T *p0, const T *p1, @@ -754,9 +772,13 @@ void Geometry::transposeMatrix(const std::vector> &a, #define GEOMETRY_SPECIALIZE(TYPE) \ template TYPE Geometry::angle( \ TYPE const *, TYPE const *, TYPE const *, TYPE const *); \ + template TYPE Geometry::angle2D( \ + TYPE const *, TYPE const *, TYPE const *, TYPE const *); \ template bool Geometry::areVectorsColinear( \ TYPE const *, TYPE const *, TYPE const *, TYPE const *, \ std::array *, TYPE const *); \ + template bool Geometry::isTriangleColinear2D( \ + TYPE const *, TYPE const *, TYPE const *, TYPE const); \ template int Geometry::computeBarycentricCoordinates( \ TYPE const *, TYPE const *, TYPE const *, std::array &, \ int const &); \ diff --git a/core/base/geometry/Geometry.h b/core/base/geometry/Geometry.h index 2b997c6301..21c47add1c 100644 --- a/core/base/geometry/Geometry.h +++ b/core/base/geometry/Geometry.h @@ -9,6 +9,7 @@ #include #include +#include namespace ttk { @@ -22,6 +23,19 @@ namespace ttk { template T angle(const T *vA0, const T *vA1, const T *vB0, const T *vB1); + template + T angle2D(const T *vA0, const T *vA1, const T *vB0, const T *vB1); + + // Computes the BAC angle, in the interval [0, 2pi] + template + inline double angle2DUndirected(const T *vA, const T *vB, const T *vC) { + double angle = angle2D(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 *coefficients = nullptr, const T *tolerance = NULL); + template + 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 T distance(const std::vector &p0, const std::vector &p1); + template + 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::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 + T distanceFlatten(const std::vector> &p0, + const std::vector> &p1); + /// Compute the Euclidean distance between two vectors by first flattening /// them /// \param p0 xyz coordinates of the first input point. diff --git a/core/base/topoMap/CMakeLists.txt b/core/base/topoMap/CMakeLists.txt new file mode 100644 index 0000000000..3f85c6cd76 --- /dev/null +++ b/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() + + diff --git a/core/base/topoMap/TopoMap.cpp b/core/base/topoMap/TopoMap.cpp new file mode 100644 index 0000000000..cc34e3c533 --- /dev/null +++ b/core/base/topoMap/TopoMap.cpp @@ -0,0 +1,96 @@ +#include + +ttk::TopoMap::TopoMap() { + // inherited from Debug: prefix will be printed at the beginning of every msg + this->setDebugMsgPrefix("TopoMap"); +} +ttk::TopoMap::~TopoMap() = default; + +#if defined(TTK_ENABLE_QHULL) && defined(Qhull_FOUND) + +#include +bool computeConvexHull_aux(const std::vector &coords, + std::vector &res, + std::string &errMsg) { + size_t nbPoint = coords.size() / 2; + char qHullFooStr[1] = ""; // We use no Qhull options. + orgQhull::Qhull qhull; + try { + qhull.runQhull(qHullFooStr, 2, nbPoint, coords.data(), qHullFooStr); + } catch(orgQhull::QhullError &e) { + errMsg = "Error with qHull module: " + std::string(e.what()); + return false; + } + + // Qhull gives us the coordinates of the points in the convex hull. Here we + // retrive the indices of this points in the list we provided. We will also + // compute the barycenter of the points in the convex hull. + for(const auto &u : qhull.vertexList()) { + const orgQhull::QhullPoint &qhullPt = u.point(); + auto coordsCur = qhullPt.coordinates(); + for(size_t j = 0; j < coords.size() / 2; j++) { + if(fabs(coords[2 * j] - coordsCur[0]) + + fabs(coords[2 * j + 1] - coordsCur[1]) + < EpsilonDBL) { + res.push_back(j); + break; + } + } + } + + if(res.size() != qhull.vertexList().size()) { + errMsg = "Error : could not retrieve all vertices in the convex hull."; + return false; + } + return true; +} + +#else // Using boost. + +#include +#include +#include +namespace bg = boost::geometry; +BOOST_GEOMETRY_REGISTER_BOOST_TUPLE_CS(cs::cartesian) + +using Point = boost::tuple; +using Polygon = boost::geometry::model::polygon; +using Mpoints = boost::geometry::model::multi_point; + +bool computeConvexHull_aux(const std::vector &coords, + std::vector &res, + std::string &errMsg) { + Mpoints multiPoints; + Polygon hull; + size_t nbPoint = coords.size() / 2; + for(size_t i = 0; i < nbPoint; i++) { + boost::geometry::append( + multiPoints, Point(coords[2 * i], coords[2 * i + 1])); + } + + boost::geometry::convex_hull(multiPoints, hull); + + // From the coordinates of the points in the convex hull, we find their + // indices. + for(const auto &boostPt : hull.outer()) { + double coordsCur[2] = {boostPt.get<0>(), boostPt.get<1>()}; + for(size_t j = 0; j < coords.size() / 2; j++) { + if(fabs(coords[2 * j] - coordsCur[0]) + + fabs(coords[2 * j + 1] - coordsCur[1]) + < Epsilon) { + res.push_back(j); + break; + } + } + } + if(res.size() != hull.outer().size()) { + errMsg = "Error : could not retrieve all vertices in the convex hull."; + return false; + } + + // Boost closes the polygon, hence the first and the last vertices are the + // identical. We do not want a duplicate vertex. + res.pop_back(); + return true; +} +#endif diff --git a/core/base/topoMap/TopoMap.h b/core/base/topoMap/TopoMap.h new file mode 100644 index 0000000000..6569d0fbc9 --- /dev/null +++ b/core/base/topoMap/TopoMap.h @@ -0,0 +1,1010 @@ +/// \ingroup base +/// \class ttk::TopoMap +/// \author Alexandre Talon +/// \date March 2023 +/// +/// +/// This module defines the %TopoMap class that embeds points into 2D. The input +/// is a distance matrix representing the distances between any pair of points +/// in high dimension. It computes the 2D coordinates of a projection of the +/// points which preserves the distances of a minimum spanning tree computed by +/// Kruskal's algorithm. It also preserves the components Kruskal's algorithm +/// builds iteratively. +/// +/// +/// \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. + +/// \sa TopoMap + +#pragma once + +// ttk common includes +#include +#include +#include + +// STL includes +#include +#include +#include +#include +#include // For std::pair +#include + +// cos and sin computations, so we are cautious on testing equalities between +// double. +static double Epsilon{ttk::Geometry::pow(10.0, -FLT_DIG + 2)}; +static double EpsilonDBL{ttk::Geometry::pow(10.0, -DBL_DIG + 2)}; + +// Normalizes a given vector. +template +static void + computeUnitVector(const T *coordOrig, const T *coordDest, T *coordVect); + +// Rotates the first argument by angle around center. +template +static void rotate(T *ptToRotate, const T *center, double angle); + +// Rotate the set of points by the given angle, considering the given center as +// center of the rotation. +template +static void + rotatePolygon(std::vector &coords, T *centerCoords, const double angle); + +// Real call to the convex hull engine: either Qhull or Boost. +bool computeConvexHull_aux(const std::vector &coords, + std::vector &res, + std::string &errMsg); + +namespace ttk { + + /** + * The TopoMap class computes a 2D projection of the points given as input, + * either by their distance matrix or by their coordinates. + */ + class TopoMap : virtual public Debug { + + public: + /** We rely on constructing a minimal spanning tree. We allow for the two + * main algorithms to build one: Kruskal and Prim. + */ + enum class STRATEGY { KRUSKAL = 0, PRIM = 1 }; + + TopoMap(); + TopoMap(size_t angularSampleNb, bool checkMST, STRATEGY strategy) + : AngularSampleNb(angularSampleNb), CheckMST(checkMST), + Strategy(strategy) { + + // inherited from Debug: prefix will be printed at the beginning of every + // msg + this->setDebugMsgPrefix("TopoMap"); + } + ~TopoMap() override; + + /** + * @brief Computes the TopoMap projection + * + * @param[out] outputCoords the final coordinates of the points, computed by + TopoMap + * + * @param[out] insertionTime an optional array indicating for each point at + what time it was inserted in the projection. + * + * @param[in] inputMatrix the high-dimension input distance matrix or + coordinates of the points + * + * @param[in] isDistMat equal to true if inputMatrix is a distance matrix, + false if it stores input point coordinates + * + * @param[in] n the number of input points + + * @return 0 in case of success. + */ + template + int execute(T *outputCoords, + int *insertionTime, + const std::vector &inputMatrix, + bool isDistMat, + size_t n); + + protected: + size_t AngularSampleNb{2}; + bool CheckMST{false}; + bool errorConvexHull{false}; + STRATEGY Strategy{STRATEGY::KRUSKAL}; + + private: + // Fills its last two arguments with the position of the end of the previous + // (resp. next) edge so that idCenter is the point at the angle at the + // intersection of the previous and next edges. + template + bool getPrevNextEdges(const std::vector &idsPtsPolygon, + size_t idCenter, + const T *allCoords, + T *coordPrev, + T *coordPost) const; + + // Tries to find the best angle of rotation for the two components. Updates + // the coordiates of their vertices accordingly. + template + T rotateMergingCompsBest(const std::vector &hull1, + const std::vector &hull2, + const std::vector &comp1, + const std::vector &comp2, + size_t iPt1, + size_t iPt2, + const std::vector &distMatrix, + T *allCoords, + size_t n, + size_t angularSampleNb); + + template + bool computeConvexHull(T *allCoords, + const std::vector &compPtsIds, + std::vector &idsInHull) const; + }; + + // Sketch of the algorithm: + // 1. Sort the edges + // 2. For each edge uv with cost c_uv + // a. Get connected components comp(u), comp(v) + // + // b. Compute convex hull of comp(u) and of (comp(v)) + // + // c. Find in hull(u) (resp. hull(v)) the point pU (resp. pV) in the + // convex hull which is closest to u (resp. v) in high dimension, u if + // possible. + // + // d. Translate comp(u) such that pU lies on the line of the bissector of + // angleV, at distance c_uv of pV, and on the opposite direction of + // hull(v). angleU (resp. angleV) is the angle of hull(u) (resp. hull(v)) + // at vertex pU (resp. pV). + // + // e. Rotate comp(u) such that the bissector of angleU prolongates the + // one of angleU, and such that pU and pV are closest than any pair of + // points in comp(u),comp(v): the two components are facing each other at + // pU and pV. + // + // + // f. Try several rotations of comp(u) and comp(v) such that pU and pV + // still are the closest pair of vertices from the two components. We + // keep the angle of rotation which minimizes the difference between the + // high dimension distance matrix and the new distance matrix between the + // two components. + // + // g. comp(u) and comp(v) were merged, we now consider them as a single + // component. + + template + int ttk::TopoMap::execute(T *outputCoords, + int *insertionTime, + const std::vector &inputMatrix, + bool isDistMat, + size_t n) { + ttk::Timer timer; + +#ifdef TTK_ENABLE_QHULL +#ifndef Qhull_FOUND + this->printWrn( + "Qhull was enabled but it is not installed or was not found."); + this->printWrn("Defaulting to Boost support instead."); +#endif +#endif + +#ifndef Qhull_FOUND + this->printMsg("Using Boost for convex hulls."); + this->printWrn("Bugs have been reported in Boost's implementation."); + this->printWrn("Consider enabling Qhull instead."); +#endif + + if(AngularSampleNb < 2) { + this->printWrn("The number of angular samples set is less than 2."); + this->printWrn("Setting it to 2 (minimum possible value)."); + this->AngularSampleNb = 2; + } + // We first set the minimum value we consider as null. We multiply by 100 + // the smallest significative value because we use sqrt/cos/sin so we take + // precautionary measures. + if(std::is_same::value) { + Epsilon = ttk::Geometry::pow(10.0, -DBL_DIG + 2); + } else if(std::is_same::value) { + Epsilon = ttk::Geometry::pow(10.0, -FLT_DIG + 2); + } else { + printErr("Error, template type must be either double or float.\n"); + return 1; + } + + 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)."); + } + + std::vector 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 dist = ttk::Geometry::distance( + &inputMatrix[dim * i1], &inputMatrix[dim * i2], dim); + computedDistMatrix[i1 * n + i2] = dist; + computedDistMatrix[i2 * n + i1] = dist; + } + } + } else { + if(inputMatrix.size() != n * n) { + this->printErr("Invalid size for the distance matrix."); + return 1; + } + } + + std::vector edgesMSTBefore, + edgesMSTAfter; // Only used for minimum spanning tree checking; + edgesMSTBefore.reserve(n); + const std::vector &distMatrix + = isDistMat ? inputMatrix : computedDistMatrix; + + if(insertionTime != nullptr) + for(size_t i = 0; i < n; i++) + insertionTime[i] = -1; + + // 1. Sorting the edges + using heapType = std::pair>; + std::priority_queue, std::greater> + edgeHeap; + if(this->Strategy == STRATEGY::KRUSKAL) { + for(size_t u1 = 0; u1 < n; u1++) { + for(size_t u2 = u1 + 1; u2 < n; u2++) { + edgeHeap.push( + std::make_pair(distMatrix[u1 * n + u2], std::make_pair(u1, u2))); + } + } + } else if(this->Strategy == STRATEGY::PRIM) { + for(size_t u1 = 0; u1 < n; u1++) { + edgeHeap.push(std::make_pair(distMatrix[u1], std::make_pair(0, u1))); + } + } else { + this->printErr("Invalid stategy for the MST: only Kruskal (0) or Prim " + "(1) algorithm can be used."); + return 1; + } + + // ufVectors is the union find vector + std::map> ufToComp; + std::vector ufVector(n); + for(size_t i = 0; i < n; i++) { + ufToComp[&ufVector[i]].emplace_back(i); + } + + for(size_t i = 0; i < 2 * n; i++) { + outputCoords[i] = 0; + } + + T finalDistortion = 0; + // 2. Processing the edges + size_t nbEdgesMerged = 0; + double MSTWeight = 0; + + std::set + stillToDo; // List of vertices not yet embedded, used for Prim strategy. + if(this->Strategy == STRATEGY::PRIM) { + for(size_t i = 1; i < n; i++) { + stillToDo.insert(i); + } + } + while(nbEdgesMerged != n - 1) { + if(edgeHeap.empty()) { + this->printErr("Error building the MST. Aborting."); + return 1; + } + const auto &elt = edgeHeap.top(); + T edgeCost = elt.first; + size_t u = elt.second.first; + size_t v = elt.second.second; + edgeHeap.pop(); + + // 2.a Getting the components + UnionFind *reprU = ufVector[u].find(); + UnionFind *reprV = ufVector[v].find(); + if(reprU == reprV) { // Already in the same component + continue; + } + nbEdgesMerged++; + MSTWeight += edgeCost; + + if(insertionTime != nullptr) { + if(insertionTime[u] == -1) { + insertionTime[u] = nbEdgesMerged; + } + if(insertionTime[v] == -1) { + insertionTime[v] = nbEdgesMerged; + } + } + edgesMSTBefore.push_back(edgeCost); + std::vector &compU = ufToComp[reprU], &compV = ufToComp[reprV]; + size_t idSmall = compU.size() < compV.size() ? u : v; + size_t idBig = idSmall == u ? v : u; + std::vector &compSmall = idSmall == u ? compU : compV; + std::vector &compBig = idSmall == u ? compV : compU; + size_t nBig = compBig.size(); + size_t nSmall = compSmall.size(); + std::vector idsInHullSmall, idsInHullBig; + + // 2.b Computing the convex hull. + // We retrieve the current coordinates of the big component.. + bool statusHull = true; + statusHull + = statusHull | computeConvexHull(outputCoords, compBig, idsInHullBig); + statusHull = statusHull + | computeConvexHull(outputCoords, compSmall, idsInHullSmall); + if(!statusHull) { + this->printErr("Aborting."); + return 1; + } + + size_t idChosenBig = idsInHullBig[0], idChosenSmall = idsInHullSmall[0]; + + // 2.c We want to select, among all vertices in the convex hull, the one + // which is closest to the vertex of the edge we work on. + for(size_t vert : idsInHullBig) { + T dist = distMatrix[vert * n + idBig]; + if(dist < distMatrix[idChosenBig * n + idBig]) { + idChosenBig = vert; + } + } + + for(size_t vert : idsInHullSmall) { + T dist = distMatrix[vert * n + idSmall]; + if(dist < distMatrix[idChosenSmall * n + idSmall]) { + idChosenSmall = vert; + } + } + + size_t sizeBigHull = idsInHullBig.size(), + sizeSmallHull = idsInHullSmall.size(); + std::vector coordsBigHull(sizeBigHull * 2), + coordsSmallHull(sizeSmallHull * 2); + for(size_t iHull = 0; iHull < sizeBigHull; iHull++) { + size_t vert = idsInHullBig[iHull]; + coordsBigHull[iHull * 2] = outputCoords[vert * 2]; + coordsBigHull[iHull * 2 + 1] = outputCoords[vert * 2 + 1]; + } + for(size_t iHull = 0; iHull < sizeSmallHull; iHull++) { + size_t vert = idsInHullSmall[iHull]; + coordsSmallHull[iHull * 2] = outputCoords[vert * 2]; + coordsSmallHull[iHull * 2 + 1] = outputCoords[vert * 2 + 1]; + } + + // Identifying the angles we are working on from the convex hulls. + T coordPrevBig[2], coordPostBig[2]; + T coordPrevSmall[2], coordPostSmall[2]; + if(!getPrevNextEdges(idsInHullSmall, idChosenSmall, outputCoords, + coordPrevSmall, coordPostSmall)) { + return 1; + } + if(!getPrevNextEdges(idsInHullBig, idChosenBig, outputCoords, + coordPrevBig, coordPostBig)) { + return 1; + } + T coordPtSmall[2] = { + outputCoords[2 * idChosenSmall], outputCoords[2 * idChosenSmall + 1]}; + T coordPtBig[2] + = {outputCoords[2 * idChosenBig], outputCoords[2 * idChosenBig + 1]}; + + // Computing the angles. + double angleSmall = 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)."); + } + T coordscenterBig[2] = {coordPrevBig[0], coordPrevBig[1]}; + T coordscenterSmall[2] = {coordPrevSmall[0], coordPrevSmall[1]}; + // Computing the coordinates of the bisectors. + rotate(coordscenterSmall, coordPtSmall, angleSmall / 2); + rotate(coordscenterBig, coordPtBig, angleBig / 2); + + T unitcenterBigVect[2], unitcenterSmallVect[2]; + computeUnitVector( + &outputCoords[idChosenBig * 2], coordscenterBig, unitcenterBigVect); + computeUnitVector(&outputCoords[idChosenSmall * 2], coordscenterSmall, + unitcenterSmallVect); + + // Computing the new coordinates for the chosen point on the small + // components, and then translating the whole component. + T goalCoordChosenSmall[2] + = {outputCoords[idChosenBig * 2] - edgeCost * unitcenterBigVect[0], + outputCoords[idChosenBig * 2 + 1] - edgeCost * unitcenterBigVect[1]}; + if(sizeBigHull == 1) { + goalCoordChosenSmall[0] = outputCoords[idChosenBig * 2] + edgeCost; + goalCoordChosenSmall[1] = outputCoords[idChosenBig * 2 + 1]; + } + T smallCompMoveVect[2] + = {goalCoordChosenSmall[0] - outputCoords[idChosenSmall * 2], + goalCoordChosenSmall[1] - outputCoords[idChosenSmall * 2 + 1]}; + + T distBaryPointSmall = ttk::Geometry::distance2D( + &outputCoords[idChosenSmall * 2], coordscenterSmall); + T preFinalPosBarySmall[2] = {coordscenterSmall[0] + smallCompMoveVect[0], + coordscenterSmall[1] + smallCompMoveVect[1]}; + T finalPosBarySmall[2] + = {goalCoordChosenSmall[0] - unitcenterBigVect[0] * distBaryPointSmall, + goalCoordChosenSmall[1] - unitcenterBigVect[1] * distBaryPointSmall}; + + // Performing the translation + rotation such that the bisectors in + // compSmall and compBig are aligned, and the two components are facing + // each other, not crossing. + for(size_t curIdSmall : compSmall) { + outputCoords[curIdSmall * 2] += smallCompMoveVect[0]; + outputCoords[curIdSmall * 2 + 1] += smallCompMoveVect[1]; + + // 2.e Performing the rotation. + double rotationAngle = ttk::Geometry::angle2DUndirected( + goalCoordChosenSmall, preFinalPosBarySmall, finalPosBarySmall); + if(nSmall > 1 && std::isfinite(rotationAngle)) { + rotate( + &outputCoords[curIdSmall * 2], goalCoordChosenSmall, rotationAngle); + } + } + + // 2.f Trying several rotations of the two components and keeping the one + // which makes the new distance matrix closest to the one provided in the + // input. + if(nBig > 1) { + finalDistortion = rotateMergingCompsBest( + idsInHullSmall, idsInHullBig, compSmall, compBig, idChosenSmall, + idChosenBig, distMatrix, outputCoords, n, this->AngularSampleNb); + if(finalDistortion < -1) { + return 1; + } + } + + T finalDist = ttk::Geometry::distance2D( + &outputCoords[2 * idChosenSmall], &outputCoords[2 * idChosenBig]); + if(fabs(finalDist - edgeCost) > Epsilon) { + this->printErr( + "The distance we set is too far from the goal distance."); + } + + UnionFind *unionRepr = UnionFind::makeUnion(reprU, reprV); + UnionFind *otherRepr = (unionRepr == reprU) ? reprV : reprU; + if(unionRepr != reprU && unionRepr != reprV) { + printErr("Error with the union find module results. Aborting."); + return 1; + } + std::vector &unionSet = ufToComp[unionRepr]; + std::vector &otherSet = ufToComp[otherRepr]; + + unionSet.insert(unionSet.end(), otherSet.begin(), otherSet.end()); + + ufToComp.erase(otherRepr); + + if(this->Strategy == STRATEGY::PRIM) { + stillToDo.erase(v); + for(size_t uToDo : stillToDo) { + edgeHeap.push(std::make_pair( + distMatrix[v * n + uToDo], std::make_pair(v, uToDo))); + } + } + } + std::stringstream ssDistortion; + ssDistortion << std::scientific << 2 * finalDistortion; + this->printMsg("Non normalized distance matrix distortion: " + + ssDistortion.str()); + this->printMsg( + "(for normalized distortion values, use distanceMatrixDistorsion)."); + + if(this->errorConvexHull) { + this->printWrn( + "Warning, somme convex hull computations were wrong, so the projection " + "may be less optimal than it could have been. If you are using Boost, " + "consider switching to Qhull to avoid this."); + } + this->printMsg("The minimum spanning tree has weight " + + std::to_string(MSTWeight) + "."); + + // We check that the lengths of the edges selected to build a minimum + // spanning tree are preserved by our algorithm, as should be. + if(CheckMST) { + + edgesMSTAfter.reserve(n); + this->printMsg("Checking the new minimum spanning tree."); + std::priority_queue, + std::greater> + edgeHeapAfter; + 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::distance2D( + &outputCoords[2 * u1], &outputCoords[2 * u2]), + std::make_pair(u1, u2))); + } + } + } else if(this->Strategy == STRATEGY::PRIM) { + for(size_t u1 = 0; u1 < n; u1++) { + edgeHeapAfter.push(std::make_pair( + ttk::Geometry::distance2D(&outputCoords[0], &outputCoords[2 * u1]), + std::make_pair(0, u1))); + } + } + + std::map> ufToCompAfter; + std::vector ufVectorAfter(n); + std::vector ufPtrVectorAfter(n); + for(size_t i = 0; i < n; i++) { + ufPtrVectorAfter[i] = &ufVectorAfter[i]; + ufToCompAfter[ufPtrVectorAfter[i]].push_back(i); + } + + nbEdgesMerged = 0; + stillToDo.clear(); + for(size_t i = 1; i < n; i++) + stillToDo.insert(i); + while(nbEdgesMerged != n - 1) { + if(edgeHeapAfter.empty()) { + this->printErr("Error building the MST for the check. Aborting."); + return 1; + } + const auto &elt = edgeHeapAfter.top(); + T edgeCost = elt.first; + size_t u = elt.second.first; + size_t v = elt.second.second; + edgeHeapAfter.pop(); + + UnionFind *reprU = ufPtrVectorAfter[u]->find(); + UnionFind *reprV = ufPtrVectorAfter[v]->find(); + if(reprU == reprV) { // Already in the same component + continue; + } + UnionFind *unionRepr = UnionFind::makeUnion(reprU, reprV); + UnionFind *otherRepr = (unionRepr == reprU) ? reprV : reprU; + if(unionRepr != reprU && unionRepr != reprV) { + printErr("Error with the union find module results. Aborting."); + return 1; + } + nbEdgesMerged++; + std::vector &unionSet = ufToCompAfter[unionRepr]; + std::vector &otherSet = ufToCompAfter[otherRepr]; + + unionSet.insert(unionSet.end(), otherSet.begin(), otherSet.end()); + + ufPtrVectorAfter[u] = unionRepr; + ufPtrVectorAfter[v] = unionRepr; + + ufToCompAfter.erase(otherRepr); + if(this->Strategy == STRATEGY::PRIM) { + stillToDo.erase(v); + for(size_t uToDo : stillToDo) { + edgeHeapAfter.push( + std::make_pair(ttk::Geometry::distance2D( + &outputCoords[2 * v], &outputCoords[2 * uToDo]), + std::make_pair(v, uToDo))); + } + } + + edgesMSTAfter.push_back(edgeCost); + } + + const T epsilonCheck = Epsilon * 5; + double sumDiff = 0, maxDiff = 0; + size_t nbProblematic = 0; + if(this->Strategy == STRATEGY::PRIM) { + sort(edgesMSTBefore.begin(), edgesMSTBefore.end()); + sort(edgesMSTAfter.begin(), edgesMSTAfter.end()); + } + for(size_t i = 0; i < edgesMSTBefore.size(); i++) { + if(fabs(edgesMSTBefore[i] - edgesMSTAfter[i]) >= epsilonCheck) { + nbProblematic++; + double diff = fabs(edgesMSTBefore[i] - edgesMSTAfter[i]); + sumDiff += diff; + maxDiff = std::max(diff, maxDiff); + } + } + if(nbProblematic != 0) { + this->printWrn("Error on " + std::to_string(nbProblematic) + + " edge(s) of the MST."); + + std::stringstream ssMax, ssAvg; + ssMax << std::scientific << maxDiff; + ssAvg << std::scientific << sumDiff / nbProblematic; + this->printWrn("The maximum error is " + ssMax.str() + + " and the average (on the number of errors) error is " + + ssAvg.str() + "."); + } + } + + this->printMsg(ttk::debug::Separator::L2); // horizontal '-' separator + this->printMsg("Complete", 1, timer.getElapsedTime()); + this->printMsg(ttk::debug::Separator::L1); // horizontal '=' separator + + return 0; + } + + template + bool TopoMap::getPrevNextEdges(const std::vector &idsPtsPolygon, + size_t idCenter, + const T *allCoords, + T *coordPrev, + T *coordPost) const { + size_t n = idsPtsPolygon.size(); + size_t iPtPrev = 0, iPtPost = 0; + bool found = false; + + for(size_t i = 0; i < n; i++) { + if(idsPtsPolygon[i] == idCenter) { + found = true; + iPtPost = idsPtsPolygon[(i + 1) % n]; + iPtPrev = idsPtsPolygon[(i + n - 1) % n]; + break; + } + } + if(!found) { + printErr("Error, we could not find the edges incident to the point we " + "chose for the rotation of the component. Aborting."); + return false; + } + + coordPrev[0] = allCoords[2 * iPtPrev]; + coordPrev[1] = allCoords[2 * iPtPrev + 1]; + coordPost[0] = allCoords[2 * iPtPost]; + coordPost[1] = allCoords[2 * iPtPost + 1]; + + return true; + } + + template + T TopoMap::rotateMergingCompsBest(const std::vector &hull1, + const std::vector &hull2, + const std::vector &comp1, + const std::vector &comp2, + size_t iPt1, + size_t iPt2, + const std::vector &distMatrix, + T *allCoords, + size_t n, + size_t angularSampleNb) { + // The distance between the two components. + T shortestDistPossible + = ttk::Geometry::distance2D(&allCoords[2 * iPt1], &allCoords[2 * iPt2]); + T coordPt1[2] = {allCoords[2 * iPt1], allCoords[2 * iPt1 + 1]}; + T coordPt2[2] = {allCoords[2 * iPt2], allCoords[2 * iPt2 + 1]}; + size_t comp1Size = comp1.size(), comp2Size = comp2.size(); + + T coordPrev1[2], coordPost1[2]; + T coordPrev2[2], coordPost2[2]; + getPrevNextEdges(hull1, iPt1, allCoords, coordPrev1, coordPost1); + getPrevNextEdges(hull2, iPt2, allCoords, 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."); + return -2; + } + T coordBissect1[2] = {coordPrev1[0], coordPrev1[1]}; + T coordBissect2[2] = {coordPrev2[0], coordPrev2[1]}; + rotate(coordBissect1, coordPt1, angle1 / 2); + rotate(coordBissect2, coordPt2, angle2 / 2); + + // If a component has only one vertex, no need to try several angles. + double semiAngle1 = comp1Size == 1 ? M_PI / 2 : angle1 / 2; + double semiAngle2 = comp2Size == 1 ? M_PI / 2 : angle2 / 2; + + double angleMax1 = M_PI / 2 - semiAngle1, angleMin1 = -angleMax1; + double angleMax2 = M_PI / 2 - semiAngle2, angleMin2 = -angleMax2; + double step1 = (angleMax1 - angleMin1) / angularSampleNb, + step2 = (angleMax2 - angleMin2) / angularSampleNb; + double bestAnglePair[2] = {0, 0}; + T bestScore = 3.e38; // Near the maximum value for float + + std::vector> origDistMatrix(comp1Size); + for(size_t i = 0; i < comp1Size; i++) { + origDistMatrix[i].resize(comp2Size); + for(size_t j = 0; j < comp2Size; j++) { + origDistMatrix[i][j] = distMatrix[comp1[i] * n + comp2[j]]; + } + } + std::vector initialCoords1(2 * comp1Size), initialCoords2(2 * comp2Size); + for(size_t i = 0; i < comp1.size(); i++) { + size_t iPt = comp1[i]; + initialCoords1[2 * i] = allCoords[2 * iPt]; + initialCoords1[2 * i + 1] = allCoords[2 * iPt + 1]; + } + for(size_t i = 0; i < comp2.size(); i++) { + size_t iPt = comp2[i]; + initialCoords2[2 * i] = allCoords[2 * iPt]; + initialCoords2[2 * i + 1] = allCoords[2 * iPt + 1]; + } + + size_t nbIter1 = std::isfinite(step1) ? angularSampleNb + 1 : 1; + size_t nbIter2 = std::isfinite(step2) ? angularSampleNb + 1 : 1; + + if(step1 * nbIter1 < 0.001) { // No need to split such a small angle + nbIter1 = 1; + } + if(step2 * nbIter2 < 0.001) { // No need to split such a small angle + nbIter2 = 1; + } + + bool errorDuringLoop = false; + bool foundValidAngle = false; +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(this->threadNumber_) shared(allCoords) +#endif + // We enumerate the rotation angles for the first component. + for(size_t i1 = 0; i1 < nbIter1; i1++) { + std::vector coords1Test(2 * comp1Size), coords2Test(2 * comp2Size); + coords1Test = initialCoords1; + double testAngle1 = angleMin1 + step1 * i1; + rotatePolygon(coords1Test, coordPt1, testAngle1); + + // For each first rotation angle, we enumerate the rotations angles for + // the second componenent. + for(size_t i2 = 0; i2 < nbIter2; i2++) { + bool curAngleError = false; + coords2Test = initialCoords2; + double testAngle2 = angleMin2 + i2 * step2; + rotatePolygon(coords2Test, coordPt2, testAngle2); + + // We compute a distortion score if we rotate the two components by the + // current angles. + T curScore = 0; + for(size_t i = 0; i < comp1Size; i++) { + T coordARotate[2] = {coords1Test[2 * i], coords1Test[2 * i + 1]}; + for(size_t j = 0; j < comp2Size; j++) { + T coordBRotate[2] = {coords2Test[2 * j], coords2Test[2 * j + 1]}; + T newDist = ttk::Geometry::distance2D(coordARotate, coordBRotate); + curScore += (newDist - origDistMatrix[i][j]) + * (newDist - origDistMatrix[i][j]); + if(newDist + Epsilon < shortestDistPossible) { + curAngleError = true; + errorDuringLoop = true; + break; + } + } + } + + if(curAngleError) { + continue; + } + + // Now if the score is the best so far, we keep this pair angle as the + // best one we tried so far. +#ifdef TTK_ENABLE_OPENMP +#pragma omp critical +#endif + { + // In case of errors in convex hull computation, we still have a valid + // angle (to preserve the MST). + foundValidAngle = true; + // This pair of angle minimises the distortion and the choice is + // deterministic with threads. + if(curScore < bestScore + || (curScore == bestScore + && (testAngle1 < bestAnglePair[0] + || (testAngle1 == bestAnglePair[0] + && testAngle2 < bestAnglePair[1])))) { + bestScore = curScore; + bestAnglePair[0] = std::isfinite(testAngle1) ? testAngle1 : 0; + bestAnglePair[1] = std::isfinite(testAngle2) ? testAngle2 : 0; + } + } + } + } + + // No valid angle were found, because the convex hull is invalid; + if(!foundValidAngle) { + this->printErr( + "No valid angle was found, due to errors in the convex hull " + "computation. Please consider using Qhull instead of Boost. Aborting."); + return -10; + } + + // Applying the chosen rotations to the two components. + for(size_t i1 : comp1) { + rotate(&allCoords[2 * i1], coordPt1, bestAnglePair[0]); + } + for(size_t i2 : comp2) { + rotate(&allCoords[2 * i2], coordPt2, bestAnglePair[1]); + } + + if(errorDuringLoop) { + this->errorConvexHull = true; + } + + return bestScore; + } + + template + bool TopoMap::computeConvexHull(T *allCoords, + const std::vector &compPtsIds, + std::vector &idsInHull) const { + size_t nbPoint = compPtsIds.size(); + std::vector compCoords(2 * nbPoint); + for(size_t i = 0; i < nbPoint; i++) { + size_t vertId = compPtsIds[i]; + compCoords[2 * i] = static_cast(allCoords[2 * vertId]); + compCoords[2 * i + 1] = static_cast(allCoords[2 * vertId + 1]); + } + + if(nbPoint <= 2) { + idsInHull.push_back(compPtsIds[0]); + if(nbPoint == 2) { + double dist = ttk::Geometry::distance2D(&compCoords[0], &compCoords[2]); + + if(dist > Epsilon) { + idsInHull.push_back(compPtsIds[1]); + } + } + return true; + } + + // Testing if all points are colinear. If so, we take extremal ones. We do + // this to avoid that the convex hull algorithm does not detect it and gives + // a wrong result. + double dirVect[2] = {0, 0}; + bool areColinear = true; + + size_t idFirstDistinct = 0; + while(idFirstDistinct < nbPoint && fabs(dirVect[0]) < Epsilon + && fabs(dirVect[1]) < Epsilon) { + idFirstDistinct++; + dirVect[0] = compCoords[2 * idFirstDistinct] - compCoords[0]; + dirVect[1] = compCoords[2 * idFirstDistinct + 1] - compCoords[1]; + if(fabs(dirVect[0]) < Epsilon) { + dirVect[0] = 0; + } + if(fabs(dirVect[1]) < Epsilon) { + dirVect[1] = 0; + } + } + + // Minimum/maximum values for x and y coordinates. + int idMins[2] = {compCoords[0] < compCoords[2] ? 0 : 1, + compCoords[1] < compCoords[3] ? 0 : 1}; + int idMaxs[2] = {compCoords[0] > compCoords[2] ? 0 : 1, + compCoords[1] > compCoords[3] ? 0 : 1}; + + const double *pt0 = &compCoords[0]; + const double *ptDistinct = &compCoords[2 * idFirstDistinct]; + + for(size_t iPt = idFirstDistinct + 1; iPt < nbPoint; iPt++) { + double curVect[2] = {compCoords[2 * iPt] - compCoords[0], + compCoords[2 * iPt + 1] - compCoords[1]}; + if(fabs(curVect[0]) < Epsilon && fabs(curVect[1]) < Epsilon) { + continue; + } + const double *ptCur = &compCoords[2 * iPt]; + if(!ttk::Geometry::isTriangleColinear2D( + pt0, ptDistinct, ptCur, EpsilonDBL)) { + areColinear = false; + break; + } + + // Updating the min and max values for x and y coordinates. + if(compCoords[2 * iPt] < compCoords[2 * idMins[0]]) { + idMins[0] = iPt; + } + if(compCoords[2 * iPt + 1] < compCoords[2 * idMins[1]]) { + idMins[1] = iPt; + } + if(compCoords[2 * iPt] > compCoords[2 * idMaxs[0]]) { + idMaxs[0] = iPt; + } + if(compCoords[2 * iPt + 1] > compCoords[2 * idMaxs[1]]) { + idMaxs[1] = iPt; + } + } + + if(areColinear) { + if(fabs(dirVect[0]) > Epsilon) { + idsInHull.push_back(compPtsIds[idMins[0]]); + idsInHull.push_back(compPtsIds[idMaxs[0]]); + } else { + idsInHull.push_back(compPtsIds[idMins[1]]); + idsInHull.push_back(compPtsIds[idMaxs[1]]); + } + return true; + } + + // Computing the convex hull using Boost or Qhull according to the build. + std::string errMsg; + bool status = computeConvexHull_aux(compCoords, idsInHull, errMsg); + if(!status) { + this->printErr(errMsg); + return false; + } + + // In its C++ API, Qhull does not return the points of the convex hull in + // (anti-)clockwise order. Boost does not seem consistent with the + // orientation of the hull. So we reorder the points by sorting them + // according to the angle they form with the center and a fictitious point + // at its right. + double sumX = 0, sumY = 0; + for(size_t u : idsInHull) { + sumX += compCoords[2 * u]; + sumY += compCoords[2 * u + 1]; + } + double bary[2] = {sumX / idsInHull.size(), sumY / idsInHull.size()}; + double baryRight[2] = {bary[0] + 2, bary[1]}; + std::vector> ptsToSort; + ptsToSort.reserve(idsInHull.size()); + for(size_t u : idsInHull) { + const double curPt[2] = {compCoords[2 * u], compCoords[2 * u + 1]}; + double curAngle + = ttk::Geometry::angle2DUndirected(bary, baryRight, curPt); + ptsToSort.emplace_back(std::make_pair(-curAngle, u)); + } + + // We sort the points according to the angle, keeping the indices + // associated. + sort(ptsToSort.begin(), ptsToSort.end()); + for(size_t i = 0; i < ptsToSort.size(); i++) + idsInHull[i] = ptsToSort[i].second; + + // We obtained the indices of the points belonging to the hull. We replace + // them by the global ids of the chosen points. + for(size_t &x : idsInHull) { + x = compPtsIds[x]; + } + return true; + } +} // namespace ttk + +template +static void + computeUnitVector(const T *coordOrig, const T *coordDest, T *coordVect) { + T tmp[2] = {coordDest[0] - coordOrig[0], coordDest[1] - coordOrig[1]}; + T dist = sqrt(tmp[0] * tmp[0] + tmp[1] * tmp[1]); + if(dist < Epsilon) { + coordVect[0] = 1; + coordVect[1] = 0; + } else { + coordVect[0] = tmp[0] / dist; + coordVect[1] = tmp[1] / dist; + } +} + +template +static void rotate(T *ptToRotate, const T *center, double angle) { + const double &xCtr = center[0], &yCtr = center[1]; + T &xPt = ptToRotate[0], &yPt = ptToRotate[1]; + const double dx = xPt - xCtr, dy = yPt - yCtr; + xPt = dx * cos(angle) - dy * sin(angle) + xCtr; + yPt = dx * sin(angle) + dy * cos(angle) + yCtr; +} + +template +static void + rotatePolygon(std::vector &coords, T *centerCoords, const double angle) { + double xCenter = centerCoords[0], yCenter = centerCoords[1]; + size_t nbPoint = coords.size() / 2; + for(size_t iPt = 0; iPt < nbPoint; iPt++) { + T &x = coords[iPt * 2], &y = coords[iPt * 2 + 1]; + double xNew, yNew; + xNew = (x - xCenter) * cos(angle) - (y - yCenter) * sin(angle) + xCenter; + yNew = (y - yCenter) * cos(angle) + (x - xCenter) * sin(angle) + yCenter; + x = xNew; + y = yNew; + } +} diff --git a/core/vtk/ttkDimensionReduction/ttkDimensionReduction.cpp b/core/vtk/ttkDimensionReduction/ttkDimensionReduction.cpp index 9571e7b7f4..212ebd9f15 100644 --- a/core/vtk/ttkDimensionReduction/ttkDimensionReduction.cpp +++ b/core/vtk/ttkDimensionReduction/ttkDimensionReduction.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -78,8 +79,17 @@ int ttkDimensionReduction::RequestData(vtkInformation *ttkNotUsed(request), outputData_.clear(); - const int errorCode = this->execute( - this->outputData_, inputData, numberOfRows, numberOfColumns); + vtkNew insertionTimeCol{}; + int *insertionPtr = nullptr; + if(this->Method == DimensionReduction::METHOD::TOPOMAP) { + insertionTimeCol->SetNumberOfTuples(numberOfRows); + insertionTimeCol->SetName("InsertionTime"); + insertionPtr = ttkUtils::GetPointer(insertionTimeCol); + } + + const int errorCode + = this->execute(this->outputData_, insertionPtr, inputData, numberOfRows, + numberOfColumns); if(!errorCode) { if(KeepAllDataArrays) @@ -92,6 +102,9 @@ int ttkDimensionReduction::RequestData(vtkInformation *ttkNotUsed(request), arr->SetName(s.data()); output->AddColumn(arr); } + if(this->Method == DimensionReduction::METHOD::TOPOMAP) { + output->AddColumn(insertionTimeCol); + } } } else { output->ShallowCopy(input); diff --git a/core/vtk/ttkDimensionReduction/ttkDimensionReduction.h b/core/vtk/ttkDimensionReduction/ttkDimensionReduction.h index 4e7666d360..0c4d247957 100644 --- a/core/vtk/ttkDimensionReduction/ttkDimensionReduction.h +++ b/core/vtk/ttkDimensionReduction/ttkDimensionReduction.h @@ -43,6 +43,12 @@ /// Generators Periodic Picture example \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 // VTK Module @@ -50,6 +56,7 @@ // TTK includes #include +#include #include #include @@ -221,6 +228,16 @@ class TTKDIMENSIONREDUCTION_EXPORT ttkDimensionReduction vtkSetMacro(pca_MaxIteration, const std::string &); vtkGetMacro(pca_MaxIteration, std::string); + // TopoMap + vtkSetMacro(topomap_AngularSampleNb, unsigned long int); + vtkGetMacro(topomap_AngularSampleNb, unsigned long int); + + vtkSetMacro(topomap_CheckMST, bool); + vtkGetMacro(topomap_CheckMST, bool); + + ttkSetEnumMacro(topomap_Strategy, ttk::TopoMap::STRATEGY); + vtkGetEnumMacro(topomap_Strategy, ttk::TopoMap::STRATEGY); + // testing vtkSetMacro(ModulePath, const std::string &); vtkGetMacro(ModulePath, std::string); diff --git a/paraview/xmls/DimensionReduction.xml b/paraview/xmls/DimensionReduction.xml index 80b0768a40..c9f55cedf0 100644 --- a/paraview/xmls/DimensionReduction.xml +++ b/paraview/xmls/DimensionReduction.xml @@ -12,7 +12,8 @@ - TTK dimensionReduction plugin documentation. + TTK filter for generic dimension reduction methods. + This filter supports various methods via the scikit-learn third party dependency (spectral embedding, local linear embedding, multi-dimensional scaling, t-SNE, isomap, PCA) as well as TopoMap (IEEE VIS 2020) Online examples: @@ -30,6 +31,11 @@ - https://topology-tool-kit.github.io/examples/persistentGenerators_periodicPicture/ + Related publication: + + "Topomap: A 0-dimensional homology preserving projection of high-dimensional data". + H.Doraiswamy, J. Tierny, P. J. S. Silva, L. G. Nonato and C. Silva. +IEEE Transactions on Visualization and Computer Graphics 27(2): 561-571, 2020. + default_values=".*" > + @@ -130,7 +136,38 @@ number_of_elements="1" default_values="2" panel_visibility="default"> - + + + + + + + + + + + + + + Set the number of output components (i.e. dimensions). @@ -142,7 +179,38 @@ number_of_elements="1" default_values="5" panel_visibility="default"> - + + + + + + + + + + + + + + Set the number of neighbors. @@ -185,7 +253,11 @@ mode="visibility" property="Method" value="4" /> - + + @@ -663,29 +735,77 @@ - - + + + + + - Set the path of the Python module. - - - - - + + + - - + + + When projecting in 2D, we regularly need to merge two already projected subsets of the points. We put them apart the right distance, but we may rotate each by some angle between two bounds. We test several angles and keep the best one. This parameter defines how many tests we do for each rotation. Beware, this parameter appears squared in the algorithm complexity. + + + + + + + + + + + This way to embedd points in 2D is based on preserving the lengths of the edges selected by Kruskal Algortihm to build a spanning tree. Checking this box enables testing that these lengths are indeed preserved. This test is very quick. + + + + + + + + + Set the path of the Python module. + + + + + + + + + Set the name of the Python module. @@ -716,6 +836,7 @@ ${DEBUG_WIDGETS} + @@ -726,7 +847,6 @@ - @@ -819,15 +939,56 @@ property="Method" value="5" /> - + + + + + + + + + + + + + + + + + + + + + + -