diff --git a/core/base/dimensionReduction/DimensionReduction.cpp b/core/base/dimensionReduction/DimensionReduction.cpp index aaa9ed6d31..08f530292d 100644 --- 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 coordsTopomap(2*nRows); + std::vector coordsTopomap(2 * nRows); topomap.execute(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 index cf0bfe3ea1..dca11f65f2 100644 --- 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 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..8f14f596b0 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,20 @@ 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 +54,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. diff --git a/core/base/topoMap/TopoMap.h b/core/base/topoMap/TopoMap.h index 39860cbf84..9a97110fa7 100644 --- 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 -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 static void computeUnitVector(const T *coordOrig, const T *coordDest, T *coordVect); -// Computes the oriented angle BAC -template -static double computeAngle(const T *ptA, const T *ptB, const T *ptC); - // Rotates the first argument by angle around center. template 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,7 +218,7 @@ 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 { @@ -239,7 +226,7 @@ namespace ttk { + std::to_string(inputMatrix.size() / n) + " dimensions)."); } -#endif + std::vector computedDistMatrix(n * n); if(!isDistMat) { size_t dim = inputMatrix.size() / n; @@ -247,13 +234,15 @@ namespace ttk { 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( + &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 -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 static void rotate(T *ptToRotate, const T *center, double angle) {