From 63005fc6957d8b03517e40a68e77657f7b569a7e Mon Sep 17 00:00:00 2001 From: Paul Leroy Date: Tue, 6 Feb 2024 10:36:20 +0100 Subject: [PATCH 1/4] begin cluster_labels function --- include/ActionA.h | 6 ++++- src/ActionA.cpp | 58 ++++++++++++++++++++++++++++------------------- 2 files changed, 40 insertions(+), 24 deletions(-) diff --git a/include/ActionA.h b/include/ActionA.h index e1ac8bf..a6eda93 100644 --- a/include/ActionA.h +++ b/include/ActionA.h @@ -26,6 +26,7 @@ private: bool sfConvertToRandomRGB(const ccHObject::Container &selectedEntities, QWidget* parent); void add_to_stack(int index, const Eigen::ArrayXi& n_donors, const Eigen::ArrayXXi& donors, std::vector& stack); int segment_labels(bool useParallelStrategy=true); + int cluster_labels(); int segment_labels_steepest_slope(bool useParallelStrategy=true); void add_to_stack_braun_willett(int index, const Eigen::ArrayXi& delta, const Eigen::ArrayXi &Di, std::vector& stack, int local_maximum); int segment_labels_braun_willett(bool useParallelStrategy=true); @@ -43,7 +44,10 @@ private: unsigned char m_bestOctreeLevel = 0; CCCoreLib::DgmOctree::NearestNeighboursSearchStruct m_nNSS; ccMainAppInterface *m_app; - Eigen::ArrayXi m_stack; + std::vector> m_stacks; + Eigen::ArrayXi m_labels; + Eigen::ArrayXi m_labelsnpoint; + Eigen::ArrayXi m_localMaximumIndexes; qG3PointDialog* m_dlg; static G3PointAction* s_g3PointAction; diff --git a/src/ActionA.cpp b/src/ActionA.cpp index 8b7a144..39401e9 100644 --- a/src/ActionA.cpp +++ b/src/ActionA.cpp @@ -214,7 +214,6 @@ int G3PointAction::segment_labels(bool useParallelStrategy) // build the stacks std::cout << "[segment_labels] build the stacks" << std::endl; Eigen::ArrayXi labels = Eigen::ArrayXi::Zero(m_cloud->size()); - Eigen::ArrayXi labelsk = Eigen::ArrayXi::Zero(m_cloud->size()); Eigen::ArrayXi labelsnpoint = Eigen::ArrayXi::Zero(m_cloud->size()); std::vector> stacks; @@ -246,7 +245,7 @@ int G3PointAction::segment_labels(bool useParallelStrategy) for (auto i : stack) { labels(i) = k; - labelsnpoint(i) = m_stack.size(); + labelsnpoint(i) = stack.size(); if (g3point_label) { g3point_label->setValue(i, k); @@ -287,6 +286,25 @@ int G3PointAction::segment_labels(bool useParallelStrategy) return nLabels; } +int G3PointAction::cluster_labels() +{ + ccLog::Print("[cluster_labels]"); + int nlabels = m_stacks.size(); + + // Compute the distances between sinks associated to each label + Eigen::ArrayXXd D1; + D1.resize(nlabels, nlabels); + for (int i = 0; i < nlabels; i++) + { + for (int j = 0; j < nlabels; j++) + { + D1(i, j) = (*m_cloud->getPoint(m_localMaximumIndexes(i)) - *m_cloud->getPoint(m_localMaximumIndexes(j))).norm(); + } + } + + return 0; +} + void G3PointAction::add_to_stack_braun_willett(int index, const Eigen::ArrayXi& delta, const Eigen::ArrayXi& Di, std::vector& stack, int local_maximum) { stack.push_back(index); @@ -324,13 +342,13 @@ int G3PointAction::segment_labels_braun_willett(bool useParallelStrategy) // if the minimum slope is positive, the receiver is a local maximum int nb_maxima = (min_slopes > 0).count(); - Eigen::ArrayXi localMaximumIndexes = Eigen::ArrayXi::Zero(nb_maxima); + Eigen::ArrayXi m_localMaximumIndexes = Eigen::ArrayXi::Zero(nb_maxima); int l = 0; for (unsigned int k = 0; k < m_cloud->size(); k++) { if (min_slopes(k) > 0) { - localMaximumIndexes(l) = k; + m_localMaximumIndexes(l) = k; receivers(k) = k; l++; } @@ -346,13 +364,11 @@ int G3PointAction::segment_labels_braun_willett(bool useParallelStrategy) Dij.push_back(list_of_donors); } std::cout << "[segment_labels_braun_willett] create di and Dij" << std::endl; - std::vector di_vec(m_cloud->size()); for (unsigned int k = 0; k < m_cloud->size(); k++) { int receiver = receivers(k); - di[receiver] = di[receiver] + 1; - di_vec[receiver] = di[receiver]; - Dij[receiver].push_back(k); + di[receiver] = di[receiver] + 1; // increment the number of donors of the receiver + Dij[receiver].push_back(k); // add the donor the the list of donors of the receiver } // build Di, the list of donors @@ -380,10 +396,6 @@ int G3PointAction::segment_labels_braun_willett(bool useParallelStrategy) // build the stacks std::cout << "[segment_labels_braun_willett] build the stacks" << std::endl; - Eigen::ArrayXi labels = Eigen::ArrayXi::Zero(m_cloud->size()); - Eigen::ArrayXi labelsk = Eigen::ArrayXi::Zero(m_cloud->size()); - Eigen::ArrayXi labelsnpoint = Eigen::ArrayXi::Zero(m_cloud->size()); - std::vector> stacks; int sfIdx = m_cloud->getScalarFieldIndexByName("g3point_label"); if (sfIdx == -1) @@ -394,9 +406,9 @@ int G3PointAction::segment_labels_braun_willett(bool useParallelStrategy) ccLog::Error("[G3Point::segment_labels] impossible to create scalar field g3point_label"); } } - CCCoreLib::ScalarField* g3point_label = m_cloud->getScalarField(sfIdx); - RGBAColorsTableType randomColors = getRandomColors(localMaximumIndexes.size()); + + RGBAColorsTableType randomColors = getRandomColors(m_localMaximumIndexes.size()); if (!m_cloud->resizeTheRGBTable(false)) { @@ -404,23 +416,23 @@ int G3PointAction::segment_labels_braun_willett(bool useParallelStrategy) return -1; } - for (int k = 0; k < localMaximumIndexes.size(); k++) + for (int k = 0; k < m_localMaximumIndexes.size(); k++) { - int localMaximumIndex = localMaximumIndexes(k); + int localMaximumIndex = m_localMaximumIndexes(k); std::vector stack; add_to_stack_braun_willett(localMaximumIndex, delta, Di, stack, localMaximumIndex); // labels for (auto i : stack) { - labels(i) = k; - labelsnpoint(i) = m_stack.size(); + m_labels(i) = k; + (i) = stack.size(); if (g3point_label) { g3point_label->setValue(i, k); m_cloud->setPointColor(i, randomColors.getValue(k)); } } - stacks.push_back(stack); + m_stacks.push_back(stack); } if (g3point_label) @@ -449,7 +461,7 @@ int G3PointAction::segment_labels_braun_willett(bool useParallelStrategy) m_app->updateUI(); } - int nLabels = localMaximumIndexes.size(); + int nLabels = m_localMaximumIndexes.size(); return nLabels; } @@ -519,7 +531,6 @@ int G3PointAction::segment_labels_steepest_slope(bool useParallelStrategy) // build the stacks std::cout << "[segment_labels] build the stacks" << std::endl; Eigen::ArrayXi labels = Eigen::ArrayXi::Zero(m_cloud->size()); - Eigen::ArrayXi labelsk = Eigen::ArrayXi::Zero(m_cloud->size()); Eigen::ArrayXi labelsnpoint = Eigen::ArrayXi::Zero(m_cloud->size()); std::vector> stacks; @@ -551,7 +562,7 @@ int G3PointAction::segment_labels_steepest_slope(bool useParallelStrategy) for (auto i : stack) { labels(i) = k; - labelsnpoint(i) = m_stack.size(); + labelsnpoint(i) = stack.size(); if (g3point_label) { g3point_label->setValue(i, k); @@ -682,7 +693,8 @@ void G3PointAction::run() m_neighbors_indexes.resize(m_cloud->size(), m_kNN); m_neighbors_distances.resize(m_cloud->size(), m_kNN); m_neighbors_slopes.resize(m_cloud->size(), m_kNN); - m_stack.resize(m_cloud->size()); + m_labels = Eigen::ArrayXi::Zero(m_cloud->size()); + m_labelsnpoint = Eigen::ArrayXi::Zero(m_cloud->size()); // Find neighbors of each point of the cloud query_neighbors(m_cloud, m_app, true); From c4ae428e0589a88386f23768e80dcb1f97316294 Mon Sep 17 00:00:00 2001 From: Paul Leroy Date: Tue, 6 Feb 2024 10:46:41 +0100 Subject: [PATCH 2/4] Update ActionA.cpp --- src/ActionA.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ActionA.cpp b/src/ActionA.cpp index 39401e9..9380698 100644 --- a/src/ActionA.cpp +++ b/src/ActionA.cpp @@ -342,7 +342,7 @@ int G3PointAction::segment_labels_braun_willett(bool useParallelStrategy) // if the minimum slope is positive, the receiver is a local maximum int nb_maxima = (min_slopes > 0).count(); - Eigen::ArrayXi m_localMaximumIndexes = Eigen::ArrayXi::Zero(nb_maxima); + m_localMaximumIndexes = Eigen::ArrayXi::Zero(nb_maxima); int l = 0; for (unsigned int k = 0; k < m_cloud->size(); k++) { @@ -425,7 +425,7 @@ int G3PointAction::segment_labels_braun_willett(bool useParallelStrategy) for (auto i : stack) { m_labels(i) = k; - (i) = stack.size(); + m_labelsnpoint(i) = stack.size(); if (g3point_label) { g3point_label->setValue(i, k); From 8c20ba9d4cc64bdcd1df211a56d792123561c432 Mon Sep 17 00:00:00 2001 From: Paul Leroy Date: Wed, 7 Feb 2024 08:32:31 +0100 Subject: [PATCH 3/4] compute_mean_angle --- include/ActionA.h | 8 +++- src/ActionA.cpp | 112 +++++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 114 insertions(+), 6 deletions(-) diff --git a/include/ActionA.h b/include/ActionA.h index a6eda93..6f6876b 100644 --- a/include/ActionA.h +++ b/include/ActionA.h @@ -26,19 +26,24 @@ private: bool sfConvertToRandomRGB(const ccHObject::Container &selectedEntities, QWidget* parent); void add_to_stack(int index, const Eigen::ArrayXi& n_donors, const Eigen::ArrayXXi& donors, std::vector& stack); int segment_labels(bool useParallelStrategy=true); + int compute_mean_angle(); int cluster_labels(); int segment_labels_steepest_slope(bool useParallelStrategy=true); void add_to_stack_braun_willett(int index, const Eigen::ArrayXi& delta, const Eigen::ArrayXi &Di, std::vector& stack, int local_maximum); int segment_labels_braun_willett(bool useParallelStrategy=true); void get_neighbors_distances_slopes(unsigned index); + void compute_node_surfaces(); bool query_neighbors(ccPointCloud* cloud, ccMainAppInterface* appInterface, bool useParallelStrategy=true); void run(); void setkNN(int kNN); + int m_kNN = 20; + double rad_factor = 0.6; + Eigen::ArrayXXi m_neighbors_indexes; Eigen::ArrayXXd m_neighbors_distances; Eigen::ArrayXXd m_neighbors_slopes; - int m_kNN = 20; + Eigen::ArrayXd m_area; ccOctree::Shared m_octree; ccPointCloud* m_cloud; unsigned char m_bestOctreeLevel = 0; @@ -48,6 +53,7 @@ private: Eigen::ArrayXi m_labels; Eigen::ArrayXi m_labelsnpoint; Eigen::ArrayXi m_localMaximumIndexes; + Eigen::ArrayXi m_ndon; qG3PointDialog* m_dlg; static G3PointAction* s_g3PointAction; diff --git a/src/ActionA.cpp b/src/ActionA.cpp index 9380698..619aca3 100644 --- a/src/ActionA.cpp +++ b/src/ActionA.cpp @@ -286,14 +286,54 @@ int G3PointAction::segment_labels(bool useParallelStrategy) return nLabels; } +int G3PointAction::compute_mean_angle() +{ + // Determine if the normals at the border of labels are similar + // Find the indexborder nodes (no donor and many other labels in the neighbourhood) + Eigen::ArrayXXi duplicated_labels(m_cloud->size(), m_kNN); + for (int n = 0; n < m_kNN; n++) + { + duplicated_labels(Eigen::all, n) = m_labels; + } + Eigen::ArrayXXi labels_of_neighbors(m_cloud->size(), m_kNN); + for (int index = 0; index < m_cloud->size(); index++) + { + for (int n = 0; n < m_kNN; n++) + { + labels_of_neighbors(index, n) = m_labels(m_neighbors_indexes(index, n)); + } + } + + Eigen::ArrayXi temp = m_kNN - (labels_of_neighbors == duplicated_labels).cast().rowwise().sum(); + auto condition = ((temp >= m_kNN / 4) && (m_ndon == 0)); + Eigen::ArrayXi indborder(condition.count()); + std::cout << "condition.count() " << condition.count() << std::endl; + int l = 0; + for (int c = 0; c < condition.size(); c++) + { + if (condition(c)) + { + indborder(l) = c; + l++; + } + } + + std::cout << "temp" << std::endl; + std::cout << condition.block(10, 0, 20, 1) << std::endl; + + std::cout << "indborder" << std::endl; + std::cout << indborder.block(0, 0, 10, 1) << std::endl; +} + int G3PointAction::cluster_labels() { ccLog::Print("[cluster_labels]"); int nlabels = m_stacks.size(); - // Compute the distances between sinks associated to each label - Eigen::ArrayXXd D1; - D1.resize(nlabels, nlabels); + std::cout << "COMPARE VALUES " << nlabels << " " << m_localMaximumIndexes.size() << std::endl; + + // Compute the distances between the sinks associated to each label + Eigen::ArrayXXd D1(nlabels, nlabels); for (int i = 0; i < nlabels; i++) { for (int j = 0; j < nlabels; j++) @@ -302,6 +342,56 @@ int G3PointAction::cluster_labels() } } + // Estimate the distances between labels using the areas + int k = 0; + Eigen::ArrayXXd D2 = Eigen::ArrayXXd::Zero(nlabels, nlabels); + Eigen::ArrayXd radius = Eigen::ArrayXd::Zero(nlabels); + for (auto &stack : m_stacks) // Radius of each label (assuming the surface corresponds to a disk) + { + radius(k) = sqrt(m_area(stack).sum() / M_PI); + k++; + } + for(int i = 0; i < nlabels; i++) // Compute inter-distances by summing radius + { + for(int j = 0; j < nlabels; j++) + { + D2(i, j) = radius(i) + radius(j); + } + } + + // If the radius of the sink is above the distance to the other sink (by a factor of rad_factor), set Dist to 1 + Eigen::ArrayXXi Dist = Eigen::ArrayXXi::Zero(nlabels, nlabels); + Dist = (rad_factor * D2 > D1).select(1, Dist); + std::cout << "Dist" << std::endl; + for (int i = 0; i < 10; i++) // set the values of the diagonal to 0 + { + Dist(i, i) = 0; + } + + // If labels are neighbours, set Nneigh to 1 + Eigen::ArrayXXi Nneigh = Eigen::ArrayXXi::Zero(nlabels, nlabels); + k = 0; + for (auto &stack : m_stacks) + { + Eigen::ArrayXXi labels(stack.size(), m_kNN); + for (int index = 0; index < stack.size(); index++) + { + for (int n = 0; n < m_kNN; n++) + { + labels(index, n) = m_labels(m_neighbors_indexes(stack[index], n)); + } + } + auto reshaped = labels.reshaped(); + std::set unique_elements(reshaped.begin(), reshaped.end()); + for (auto unique : unique_elements) + { + Nneigh(k, unique) = 1; + } + k++; + } + + compute_mean_angle(); + return 0; } @@ -368,9 +458,11 @@ int G3PointAction::segment_labels_braun_willett(bool useParallelStrategy) { int receiver = receivers(k); di[receiver] = di[receiver] + 1; // increment the number of donors of the receiver - Dij[receiver].push_back(k); // add the donor the the list of donors of the receiver + Dij[receiver].push_back(k); // add the donor to the list of donors of the receiver } + m_ndon = di; + // build Di, the list of donors Eigen::ArrayXi Di = Eigen::ArrayXi::Zero(m_cloud->size()); // list of donors int idx = 0; @@ -621,7 +713,7 @@ void G3PointAction::get_neighbors_distances_slopes(unsigned index) m_neighbors_indexes(index, k) = Yk.getPointGlobalIndex(k + 1); // compute the distance to the neighbor const CCVector3* neighbor = Yk.getPoint(k + 1); - float distance = sqrt((*P - *neighbor).norm2()); + float distance = (*P - *neighbor).norm(); m_neighbors_distances(index, k) = distance; // compute the slope to the neighbor m_neighbors_slopes(index, k) = (P->z - neighbor->z) / distance; @@ -629,6 +721,11 @@ void G3PointAction::get_neighbors_distances_slopes(unsigned index) } } +void G3PointAction::compute_node_surfaces() +{ + m_area = M_PI * m_neighbors_distances.rowwise().minCoeff().square(); +} + bool G3PointAction::query_neighbors(ccPointCloud* cloud, ccMainAppInterface* appInterface, bool useParallelStrategy) { std::cout << "[query_neighbor]" << std::endl; @@ -695,15 +792,20 @@ void G3PointAction::run() m_neighbors_slopes.resize(m_cloud->size(), m_kNN); m_labels = Eigen::ArrayXi::Zero(m_cloud->size()); m_labelsnpoint = Eigen::ArrayXi::Zero(m_cloud->size()); + m_stacks.clear(); // needed in case of several runs // Find neighbors of each point of the cloud query_neighbors(m_cloud, m_app, true); + compute_node_surfaces(); + // Perform initial segmentation // int nLabels = segment_labels(); int nLabels = segment_labels_braun_willett(); // int nLabels = segment_labels_steepest_slope(); + cluster_labels(); + m_app->dispToConsole( "[G3Point] initial segmentation: " + QString::number(nLabels) + " labels", ccMainAppInterface::STD_CONSOLE_MESSAGE ); m_neighbors_indexes.resize(0, 0); From 6830f2eda5453b0d050cb64b9ede8272be350124 Mon Sep 17 00:00:00 2001 From: Paul Leroy Date: Thu, 8 Feb 2024 00:22:02 +0100 Subject: [PATCH 4/4] try to use mlpack for dbscan --- CMakeLists.txt | 21 +++++++++++++-- include/ActionA.h | 2 ++ src/ActionA.cpp | 68 ++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 88 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 57a551b..64cfde4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -15,10 +15,27 @@ if ( PLUGIN_G3POINT ) add_subdirectory( src ) add_subdirectory( ui ) + target_compile_features(${PROJECT_NAME} PRIVATE cxx_std_17) # for mlpack + target_include_directories( ${PROJECT_NAME} PRIVATE C:/opt/eigen-3.4.0 + C:/Users/PaulLeroy/miniconda3/envs/env_4_CloudCompare/include ) - # set dependencies to necessary libraries - # target_link_libraries( ${PROJECT_NAME} LIB1 ) + # Find installed Open3D, which exports Open3D::Open3D +# find_package(Open3D REQUIRED) +# target_link_libraries( ${PROJECT_NAME} Open3D::Open3D) + + # On Windows if BUILD_SHARED_LIBS is enabled, copy .dll files to the executable directory +# if(WIN32) +# get_target_property(open3d_type Open3D::Open3D TYPE) +# if(open3d_type STREQUAL "SHARED_LIBRARY") +# message(STATUS "Copying Open3D.dll to ${CMAKE_CURRENT_BINARY_DIR}/$") +# add_custom_command(TARGET Draw POST_BUILD +# COMMAND ${CMAKE_COMMAND} -E copy +# ${CMAKE_INSTALL_PREFIX}/bin/Open3D.dll +# ${CMAKE_CURRENT_BINARY_DIR}/$) +# endif() +# endif() + endif() diff --git a/include/ActionA.h b/include/ActionA.h index 6f6876b..72e90a3 100644 --- a/include/ActionA.h +++ b/include/ActionA.h @@ -33,6 +33,8 @@ private: int segment_labels_braun_willett(bool useParallelStrategy=true); void get_neighbors_distances_slopes(unsigned index); void compute_node_surfaces(); + void orient_normals(); + void compute_normals_and_orient_them(); bool query_neighbors(ccPointCloud* cloud, ccMainAppInterface* appInterface, bool useParallelStrategy=true); void run(); void setkNN(int kNN); diff --git a/src/ActionA.cpp b/src/ActionA.cpp index 619aca3..f0815a3 100644 --- a/src/ActionA.cpp +++ b/src/ActionA.cpp @@ -22,6 +22,8 @@ #include #include +#include + namespace G3Point { G3PointAction* G3PointAction::s_g3PointAction; @@ -286,9 +288,27 @@ int G3PointAction::segment_labels(bool useParallelStrategy) return nLabels; } +class mySearch +{ +public: + mySearch() {} + + void Search(const arma::mat& queryPoints, + const mlpack::math::Range& range, + std::vector>& neighbors, + std::vector>& distances); +}; + +void mySearch::Search(const arma::mat& queryPoints, + const mlpack::math::Range& range, + std::vector>& neighbors, + std::vector>& distances) +{ + +} + int G3PointAction::compute_mean_angle() { - // Determine if the normals at the border of labels are similar // Find the indexborder nodes (no donor and many other labels in the neighbourhood) Eigen::ArrayXXi duplicated_labels(m_cloud->size(), m_kNN); for (int n = 0; n < m_kNN; n++) @@ -323,6 +343,26 @@ int G3PointAction::compute_mean_angle() std::cout << "indborder" << std::endl; std::cout << indborder.block(0, 0, 10, 1) << std::endl; + + // Compute the angle of the normal vector between the neighbours of each grain / label + int nlabels = m_stacks.size(); + Eigen::ArrayXXi A = Eigen::ArrayXXi::Zero(nlabels, nlabels); + Eigen::ArrayXXi N = Eigen::ArrayXXi::Zero(nlabels, nlabels); + + for (auto i : indborder) + { + auto j = m_neighbors_indexes(i, Eigen::all); // indexes of the neighbourhood of i + // Take the normals vector for i and j (duplicate the normal vector of i to have the same size as for j) +// P1 = numpy.tile(normals[i, :], (params.knn, 1)); +// P2 = m_normals(j, Eigen::all); + // Compute the angle between the normal of i and the normals of j + // Add this angle to the angle matrix between each label +// A[labels[i], labels[j]] = A[labels[i], labels[j]] + angle_rot_2_vec_mat(P1, P2) + // Number of occurrences +// N[labels[i], labels[j]] = N[labels[i], labels[j]] + 1 + } + + mlpack::DBSCAN(1, 1); } int G3PointAction::cluster_labels() @@ -726,6 +766,32 @@ void G3PointAction::compute_node_surfaces() m_area = M_PI * m_neighbors_distances.rowwise().minCoeff().square(); } +void G3PointAction::orient_normals() +{ + // Flip the normals so they are oriented towards the sensor center + // x,y,z: points + // u,v,w: normals + // ox,oy,oz: sensor center + +// p1 = sensor_center - points +// p2 = normals + +// Flip the normals if they are not pointing towards the sensor +//angle = np.arctan2(np.linalg.norm(np.cross(p1, p2), axis=1), np.sum(p1 * p2, axis=1)) +//index = (angle > np.pi / 2) | (angle < -np.pi / 2) +//normals[index] = -normals[index] # invert normal + +//return normals +} + +void G3PointAction::compute_normals_and_orient_them() +{ +// pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamKNN(params.knn)) +// centroid = np.mean(xyz, axis=0) +// sensor_center = np.array([centroid[0], centroid[1], 1000]) +// normals = orient_normals(xyz, np.asarray(pcd.normals), sensor_center) +} + bool G3PointAction::query_neighbors(ccPointCloud* cloud, ccMainAppInterface* appInterface, bool useParallelStrategy) { std::cout << "[query_neighbor]" << std::endl;