diff --git a/README.md b/README.md index b27ed92..b8c4edd 100644 --- a/README.md +++ b/README.md @@ -9,8 +9,10 @@ Granulometry from 3D Point clouds **G3Point** is a tool which aims at automatically measuring the size, shape, and orientation of a large number of individual grains as detected from any type of 3D point clouds describing the topography of surfaces covered by sediments. + The tool has been developped initially in *Matlab* https://github.com/philippesteer/G3Point This repository aims at making a plugin for CloudCompare from the original Matlab tool, also including later developments of the tool, either in *Matlab* (https://github.com/philippesteer/G3Point_dev) or in *Python* (https://github.com/p-leroy/g3point_python). + This algorithm relies on 3 main phases: 1. Grain **segmentation** using a waterhsed algorithm 2. Grain **merging and cleaning** diff --git a/include/ActionA.h b/include/ActionA.h index 3f56124..2349d1f 100644 --- a/include/ActionA.h +++ b/include/ActionA.h @@ -2,14 +2,17 @@ #include "ccPointCloud.h" +#include "Eigen/Dense" + #pragma once class ccMainAppInterface; namespace G3Point { -int getBestOctreeLevel(); -void get_neighbors(unsigned index); +void add_to_stack(int index, const Eigen::ArrayXi& n_donors, const Eigen::ArrayXXi& donors, std::vector& stack); +int segment_labels(bool useParallelStrategy=true); +void get_neighbors_distances_slopes(unsigned index); bool query_neighbors(ccPointCloud* cloud, ccMainAppInterface* appInterface, bool useParallelStrategy=true); void performActionA( ccMainAppInterface *appInterface ); } diff --git a/src/ActionA.cpp b/src/ActionA.cpp index ca951d5..12063ab 100644 --- a/src/ActionA.cpp +++ b/src/ActionA.cpp @@ -6,37 +6,138 @@ #include "ccOctree.h" #include "ccProgressDialog.h" #include "ccQtHelpers.h" +#include "CCGeom.h" #include "QMainWindow" #include "QCoreApplication" #include "QThreadPool" #include "QtConcurrent" -#include "Eigen/Dense" - #include "algorithm" +#include "iostream" namespace G3Point { -Eigen::MatrixXi neighbors_indexes; +Eigen::ArrayXXi neighbors_indexes; +Eigen::ArrayXXf neighbors_distances; +Eigen::ArrayXXf neighbors_slopes; int kNN = 20; ccOctree::Shared octree; ccPointCloud* cloud; unsigned char bestOctreeLevel = 0; CCCoreLib::DgmOctree::NearestNeighboursSearchStruct nNSS; -void get_neighbors(unsigned index) +void add_to_stack(int index, const Eigen::ArrayXi& n_donors, const Eigen::ArrayXXi& donors, std::vector& stack) +{ + stack.push_back(index); + + for (int k = 0; k < n_donors(index); k++) + { + add_to_stack(donors(index, k), n_donors, donors, stack); + } +} + +int segment_labels(bool useParallelStrategy) +{ + // for each point, find in the neighborhood the point with the minimum slope (the receiver) + Eigen::ArrayXf min_slopes(neighbors_slopes.rowwise().minCoeff()); + Eigen::ArrayXi index_of_min_slope = Eigen::ArrayXi::Zero(cloud->size()); + Eigen::ArrayXi receivers(cloud->size()); + + for (unsigned index = 0; index < cloud->size(); index++) + { + float min_slope = min_slopes(index); + for (int k = 0; k < kNN; k++) + { + if (neighbors_slopes(index, k) == min_slope) + { + index_of_min_slope(index) = k; + } + } + receivers(index) = neighbors_indexes(index, index_of_min_slope(index)); + } + + // 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); + int l = 0; + for (unsigned int k = 0; k < cloud->size(); k++) + { + if (min_slopes(k) > 0) + { + localMaximumIndexes(l) = k; + receivers(k) = k; + l++; + } + } + + // identify the donors for each receiver + Eigen::ArrayXi nDonors = Eigen::ArrayXi::Zero(cloud->size()); + Eigen::ArrayXXi donors = Eigen::ArrayXXi::Zero(cloud->size(), kNN); + + for (unsigned int k = 0; k < cloud->size(); k++) + { + int receiver = receivers(k); + if (receiver != k) // this receiver is not a local maximum + { + nDonors(receiver) = nDonors(receiver) + 1; + donors(receiver, nDonors(receiver) - 1) = k; + } + } + + // build the stacks + Eigen::ArrayXi labels = Eigen::ArrayXi::Zero(cloud->size()); + Eigen::ArrayXi labelsk = Eigen::ArrayXi::Zero(cloud->size()); + Eigen::ArrayXi labelsnpoint = Eigen::ArrayXi::Zero(cloud->size()); + std::vector> stacks; + + for (int k = 0; k < localMaximumIndexes.size(); k++) + { + int localMaximumIndex = localMaximumIndexes(k); + std::vector stack; + add_to_stack(localMaximumIndex, nDonors, donors, stack); + stacks.push_back(stack); + // labels + for (auto i : stack) + { + labels(i) = k; + labelsnpoint(i) = stack.size(); + } + } + + int nLabels = localMaximumIndexes.size(); + + return nLabels; +} + +void get_neighbors_distances_slopes(unsigned index) { const CCVector3* P = cloud->getPoint(index); - nNSS.level = bestOctreeLevel; + double maxSquareDist = 0; + int neighborhoodSize = 0; + CCCoreLib::ReferenceCloud Yk(cloud); - octree->findNearestNeighborsStartingFromCell(); + // get the nearest neighbors + if (octree->findPointNeighbourhood(P, &Yk, static_cast(kNN + 1), bestOctreeLevel, maxSquareDist, 0, &neighborhoodSize) >= static_cast(kNN)) + { + for (int k = 0; k < kNN; k++) + { + // store the index of the neighbor + 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()); + neighbors_distances(index, k) = distance; + // compute the slope to the neighbor + neighbors_slopes(index, k) = (P->z - neighbor->z) / distance; + } + } } -bool query_neighbors(ccMainAppInterface* appInterface, bool useParallelStrategy) +bool query_neighbors(ccPointCloud* cloud, ccMainAppInterface* appInterface, bool useParallelStrategy) { QString errorStr; @@ -67,6 +168,7 @@ bool query_neighbors(ccMainAppInterface* appInterface, bool useParallelStrategy) int maxThreadCount = 0; CCCoreLib::DgmOctree::NearestNeighboursSearchStruct nNSS; std::vector pointsIndexes; + pointsIndexes.resize(nPoints); if (useParallelStrategy) { @@ -79,14 +181,14 @@ bool query_neighbors(ccMainAppInterface* appInterface, bool useParallelStrategy) maxThreadCount = ccQtHelpers::GetMaxThreadCount(); } QThreadPool::globalInstance()->setMaxThreadCount(maxThreadCount); - QtConcurrent::blockingMap(pointsIndexes, get_neighbors); + QtConcurrent::blockingMap(pointsIndexes, get_neighbors_distances_slopes); } else { //manually call the static per-point method! for (unsigned i = 0; i < nPoints; ++i) { - get_neighbors(i); + get_neighbors_distances_slopes(i); } } @@ -103,8 +205,6 @@ void performActionA( ccMainAppInterface *appInterface ) return; } - /*** HERE STARTS THE ACTION ***/ - //we need one point cloud if (!appInterface->haveOneSelection()) { @@ -124,21 +224,21 @@ void performActionA( ccMainAppInterface *appInterface ) cloud = ccHObjectCaster::ToPointCloud(ent); + // initialize the matrix which will contain the results + neighbors_indexes.resize(cloud->size(), kNN); + neighbors_distances.resize(cloud->size(), kNN); + neighbors_slopes.resize(cloud->size(), kNN); + // Find neighbors of each point of the cloud - query_neighbors(cloud, appInterface); + query_neighbors(cloud, appInterface, true); // Perform initial segmentation + int nLabels = segment_labels(); - // This is how you can output messages - // Display a standard message in the console - appInterface->dispToConsole( "[ExamplePlugin] Hello world!", ccMainAppInterface::STD_CONSOLE_MESSAGE ); + appInterface->dispToConsole( "[G3Point] initial segmentation: " + QString::number(nLabels) + " labels", ccMainAppInterface::STD_CONSOLE_MESSAGE ); - // Display a warning message in the console - appInterface->dispToConsole( "[ExamplePlugin] Warning: example plugin shouldn't be used as is", ccMainAppInterface::WRN_CONSOLE_MESSAGE ); - - // Display an error message in the console AND pop-up an error box - appInterface->dispToConsole( "Example plugin shouldn't be used - it doesn't do anything!", ccMainAppInterface::ERR_CONSOLE_MESSAGE ); - - /*** HERE ENDS THE ACTION ***/ + neighbors_indexes.resize(0, 0); + neighbors_distances.resize(0, 0); + neighbors_slopes.resize(0, 0); } }