initial segmentation

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**

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<int>& 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 );
 }

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<int>& 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<std::vector<int>> stacks;
+
+	for (int k = 0; k < localMaximumIndexes.size(); k++)
+	{
+		int localMaximumIndex = localMaximumIndexes(k);
+		std::vector<int> 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<unsigned>(kNN + 1), bestOctreeLevel, maxSquareDist, 0, &neighborhoodSize) >= static_cast<unsigned>(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<unsigned> 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);
 }
 }