multithreading is OK

CMakeLists.txt

@@ -24,14 +24,6 @@ target_sources(G3PointPlugin
 
 	# 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
-		C:/opt/open3d-devel-windows-amd64-0.19.0/include
-		C:/opt/GeometricTools/GTE
-		C:/opt/boost_1_77_0
-	)
-	
 	# may be needed for debug
 	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /bigobj")
 
@@ -42,19 +34,40 @@ target_sources(G3PointPlugin
 		Qt5::PrintSupport
 		)
 
-	# Find installed Open3D, which exports Open3D::Open3D
+	# EIGEN
-	find_package( Open3D REQUIRED )
+	set( EIGEN_ROOT_DIR "" CACHE PATH "Eigen root (contains the Eigen directory)" )
-	target_link_libraries( ${PROJECT_NAME} Open3D::Open3D )
+	if ( NOT EIGEN_ROOT_DIR )
+		message( SEND_ERROR "No Eigen root directory specified (EIGEN_ROOT_DIR)" )
+	else()
+		message( STATUS "EIGEN_ROOT_DIR " ${EIGEN_ROOT_DIR} )
+	endif()
 
-	message( "Open3D_DIR ${Open3D_DIR}" )
+	# OPEN3D
+	option(USE_OPEN3D_WITH_G3POINT "Enable the use of Open3D with G3Point" OFF)
+	if (USE_OPEN3D_WITH_G3POINT)
+		find_package( Open3D REQUIRED ) # Find installed Open3D, which exports Open3D::Open3D
+		target_link_libraries( ${PROJECT_NAME} Open3D::Open3D )
+		if(CMAKE_BUILD_TYPE STREQUAL "Debug")
 			set( OPENCV_DEP_DLL_FILES
 				${Open3D_DIR}/../bin/Open3D.dll
-		${Open3D_DIR}/../bin/tbb12.dll
+				${Open3D_DIR}/../bin/tbb12_debug.dll)
-	)
+		elseif(CMAKE_BUILD_TYPE STREQUAL "Release" OR CMAKE_BUILD_TYPE STREQUAL "RelWithDebInfo")
+			set( OPENCV_DEP_DLL_FILES
+				${Open3D_DIR}/../bin/Open3D.dll
+				${Open3D_DIR}/../bin/tbb12.dll)
+		endif()
+		message( "Open3D_DIR ${Open3D_DIR}" )
 		copy_files( "${OPENCV_DEP_DLL_FILES}" "${CLOUDCOMPARE_DEST_FOLDER}" 1) #mind the quotes!
+	endif()
+
 	copy_files( "${CMAKE_CURRENT_SOURCE_DIR}/shaders/DrawGrains.vs" "${CLOUDCOMPARE_DEST_FOLDER}/shaders/G3Point" 1) #mind the quotes!
 	copy_files( "${CMAKE_CURRENT_SOURCE_DIR}/shaders/DrawGrains.fs" "${CLOUDCOMPARE_DEST_FOLDER}/shaders/G3Point" 1) #mind the quotes!
 
+	target_include_directories( ${PROJECT_NAME} PRIVATE
+		../../../../libs/qCC_db/extern/CCCoreLib/extern/nanoflann/include
+		${EIGEN_ROOT_DIR}
+	)
+
 	#================
 	# git commit hash
 	# Get the current working branch

include/G3PointAction.h

@@ -1,24 +1,27 @@
 #include "Eigen/Dense"
+#include <nanoflann.hpp>
 
+// qCC_db
 #include <ccOctree.h>
 #include <ccScalarField.h>
+#include <ccPointCloud.h>
+
+// CCCoreLib
+#include <DgmOctree.h>
+#include <Neighbourhood.h>
 
 #include <vector>
 
 #include <QObject>
 
 #include <G3PointDialog.h>
-
 #include <GrainsAsEllipsoids.h>
-
 #include <AnglesCustomPlot.h>
-
 #include <G3PointPlots.h>
 
 #pragma once
 
 class ccMainAppInterface;
-class ccPointCloud;
 
 namespace G3Point
 {
@@ -53,6 +56,38 @@ public:
 
 	template<typename T> static bool EigenArrayToFile(QString name, T array);
 
+
+    // A small adaptor to let nanoflann access ccPointCloud data
+    struct CloudAdaptor
+    {
+        const ccPointCloud* cloud;
+
+        CloudAdaptor(const ccPointCloud* c) : cloud(c) {}
+
+        // Must return the number of data points
+        inline size_t kdtree_get_point_count() const { return cloud->size(); }
+
+        // Returns the dim'th component of the idx'th point
+        inline float kdtree_get_pt(const size_t idx, int dim) const
+        {
+            if (dim == 0)
+                return cloud->getPoint(static_cast<unsigned>(idx))->x;
+            else if (dim == 1)
+                return cloud->getPoint(static_cast<unsigned>(idx))->y;
+            else
+                return cloud->getPoint(static_cast<unsigned>(idx))->z;
+        }
+
+        // Optional bounding-box computation: return false to default to a standard bbox computation loop.
+        template<class BBOX>
+        bool kdtree_get_bbox(BBOX&) const { return false; }
+    };
+
+    // Typedef for a 3D KD-tree index
+    using KDTree = nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<float, CloudAdaptor>,
+                                                       CloudAdaptor,
+                                                       3 /* dim */>;
+
 private:
 	bool sfConvertToRandomRGB(const ccHObject::Container &selectedEntities, QWidget* parent);
 	void addToStack(int index, const Eigen::ArrayXi& n_donors, const Eigen::ArrayXXi& donors, std::vector<int>& stack);
@@ -70,6 +105,10 @@ private:
 	bool computeNormalsAndOrientThemWithCloudCompare();
 	void orientNormals(const Eigen::Vector3d &sensorCenter);
 	bool computeNormalsWithOpen3D();
+	bool findNearestNeighborsNanoFlann(const unsigned int globalIndex, CCCoreLib::ReferenceCloud *points, const KDTree *kdTree);
+	bool computeNormWithFlann(unsigned int index, NormsTableType* theNorms, const KDTree *kdTree);
+    bool computeNormalsWithCloudCompare();
+    bool computeNormals();
 	bool queryNeighbors(ccPointCloud* cloud, ccMainAppInterface* appInterface, bool useParallelStrategy=true);
 	void init();
 	void showDlg();

src/G3Point.cpp

@@ -15,22 +15,6 @@
 //#                                                                        #
 //##########################################################################
 
-// First:
-//	Replace all occurrences of 'G3PointPlugin' by your own plugin class name in this file.
-//	This includes the resource path to info.json in the constructor.
-
-// Second:
-//	Open G3PointPlugin.qrc, change the "prefix" and the icon filename for your plugin.
-//	Change the name of the file to <yourPluginName>.qrc
-
-// Third:
-//	Open the info.json file and fill in the information about the plugin.
-//	 "type" should be one of: "Standard", "GL", or "I/O" (required)
-//	 "name" is the name of the plugin (required)
-//	 "icon" is the Qt resource path to the plugin's icon (from the .qrc file)
-//	 "description" is used as a tootip if the plugin has actions and is displayed in the plugin dialog
-//	 "authors", "maintainers", and "references" show up in the plugin dialog as well
-
 #include <QtGui>
 
 #include "G3Point.h"

src/G3PointAction.cpp

@@ -1,4 +1,6 @@
 #include "G3PointAction.h"
+#include "DgmOctreeReferenceCloud.h"
+#include "Neighbourhood.h"
 
 // CCPluginAPI
 #include <ccMainAppInterface.h>
@@ -33,7 +35,9 @@
 #include <random>
 
 // Open3D
+#ifdef USE_OPEN3D_WITH_G3POINT
 #include <open3d/geometry/PointCloud.h>
+#endif
 
 // Eigen
 #include <Eigen/Geometry>
@@ -369,7 +373,7 @@ Eigen::ArrayXXd G3PointAction::computeMeanAngleBetweenNormalsAtBorders()
 	Eigen::ArrayXXi duplicated_labels(m_cloud->size(), m_kNN);
 	for (int n = 0; n < m_kNN; n++)
 	{
-		duplicated_labels(Eigen::all, n) = m_labels;
+		duplicated_labels(Eigen::placeholders::all, n) = m_labels;
 	}
 	Eigen::ArrayXXi labels_of_neighbors(m_cloud->size(), m_kNN);
 	for (int index = 0; index < static_cast<int>(m_cloud->size()); index++)
@@ -403,12 +407,12 @@ Eigen::ArrayXXd G3PointAction::computeMeanAngleBetweenNormalsAtBorders()
 
 	for (auto i : indborder)
 	{
-		auto neighbors = m_neighborsIndexes(i, Eigen::all);  // indexes of the neighbors of i
+		auto neighbors = m_neighborsIndexes(i, Eigen::placeholders::all);  // indexes of the neighbors of i
-		Eigen::Vector3d N1(m_normals(i, Eigen::all)); // normal at i
+		Eigen::Vector3d N1(m_normals(i, Eigen::placeholders::all)); // normal at i
 		for (auto j : neighbors)
 		{
 			// Take the normals vector for i and j
-			Eigen::Vector3d N2(m_normals(j, Eigen::all)); // normal at j
+			Eigen::Vector3d N2(m_normals(j, Eigen::placeholders::all)); // normal at j
 			double angle = angleRot2VecMat(N1, N2);
 			if ((m_labels(i) != -1) && (m_labels(j) != -1))  // points which belong to the discarded grains have the -1 label
 			{
@@ -1263,21 +1267,21 @@ bool G3PointAction::wolman()
 
 	Eigen::ArrayXXf dq(n_iter, 3);
 	Eigen::ArrayXf d_sample = d[0];
-	dq(0, Eigen::all) << quant(d[0], 0.1), quant(d[0], 0.5), quant(d[0], 0.9);
+	dq(0, Eigen::placeholders::all) << quant(d[0], 0.1), quant(d[0], 0.5), quant(d[0], 0.9);
 	for (int i = 1; i < n_iter; i++)
 	{
 		Eigen::ArrayXf tmp(d_sample.size() + d[i].size());
 		tmp << d_sample, d[i];
 		d_sample = tmp;
-		dq(i, Eigen::all) << quant(d[i], 0.1), quant(d[i], 0.5), quant(d[i], 0.9);
+		dq(i, Eigen::placeholders::all) << quant(d[i], 0.1), quant(d[i], 0.5), quant(d[i], 0.9);
 	}
 
 	// std::cout << "d_sample " << d_sample << std::endl;
 
 	// compute standard deviation
-	Eigen::Array3d edq {std_dev(dq(Eigen::all, 0)),
+	Eigen::Array3d edq {std_dev(dq(Eigen::placeholders::all, 0)),
-					   std_dev(dq(Eigen::all, 1)),
+					   std_dev(dq(Eigen::placeholders::all, 1)),
-					   std_dev(dq(Eigen::all, 2))};
+					   std_dev(dq(Eigen::placeholders::all, 2))};
 	Eigen::Array3d dq_final {quant(d_sample, 0.1),
 							quant(d_sample, 0.5),
 							quant(d_sample, 0.9)};
@@ -1430,11 +1434,11 @@ bool G3PointAction::cleanLabels()
 			Eigen::RowVector3d centroid = points.colwise().mean();
 			points.rowwise() -= centroid;
 			// SVD decomposition A = U S V∗
-			s(k, Eigen::all) = points.jacobiSvd().singularValues();
+			s(k, Eigen::placeholders::all) = points.jacobiSvd().singularValues();
 		}
 		// filtering condition: (l2 / l0 > min_flatness) or (l1 / l0 > 2 * min_flatness)
-		Xb condition = (s(Eigen::all, 2) / s(Eigen::all, 0) > m_minFlatness)
+		Xb condition = (s(Eigen::placeholders::all, 2) / s(Eigen::placeholders::all, 0) > m_minFlatness)
-					   || (s(Eigen::all, 1) / s(Eigen::all, 0) > 2. * m_minFlatness);
+					   || (s(Eigen::placeholders::all, 1) / s(Eigen::placeholders::all, 0) > 2. * m_minFlatness);
 		size_t numberOfGrainsToKeep = condition.count();
 		if (numberOfGrainsToKeep == m_stacks.size())
 		{
@@ -1742,7 +1746,7 @@ void G3PointAction::orientNormals(const Eigen::Vector3d& sensorCenter)
 	{
 		const CCVector3 *point = m_cloud->getPoint(i);
 		Eigen::Vector3d P1 = sensorCenter - Eigen::Vector3d(point->x, point->y, point->z);
-		Eigen::Vector3d P2 = m_normals(i, Eigen::all);
+		Eigen::Vector3d P2 = m_normals(i, Eigen::placeholders::all);
 		double angle = atan2(P1.cross(P2).norm(), P1.dot(P2));
 		if ((angle < - M_PI / 2) || (angle > M_PI / 2))
 		{
@@ -1755,6 +1759,7 @@ void G3PointAction::orientNormals(const Eigen::Vector3d& sensorCenter)
 
 bool G3PointAction::computeNormalsWithOpen3D()
 {
+#ifdef USE_OPEN3D_WITH_G3POINT
 	// create an open3D point cloud from the original point cloud
 	std::vector<Eigen::Vector3d> points(m_cloud->size());
 	for (int index =0; index < points.size(); index++) // copy all points
@@ -1807,6 +1812,170 @@ bool G3PointAction::computeNormalsWithOpen3D()
 	}
 
 	return true;
+#else
+	return false;
+#endif
+}
+
+bool G3PointAction::findNearestNeighborsNanoFlann(const unsigned globalIndex,
+	                                              CCCoreLib::ReferenceCloud* points,
+	                                              const KDTree* kdTree)
+{
+	// Prepare query
+	const CCVector3* Q = m_cloud->getPoint(globalIndex);
+	float query[3] = { Q->x, Q->y, Q->z };
+
+	std::vector<size_t> retIndexes(m_kNN);
+	std::vector<float> outDistsSqr(m_kNN);
+
+		   // Perform search
+	nanoflann::KNNResultSet<float> resultSet(m_kNN);
+	resultSet.init(&retIndexes[0], &outDistsSqr[0]);
+	if(kdTree->findNeighbors(resultSet, &query[0]))
+	{
+		points->resize(m_kNN);
+		for (int i = 0; i < m_kNN; ++i)
+		{
+			points->setPointIndex(i, retIndexes[i]);
+		}
+		return true;
+	}
+	else
+	{
+		return false;
+	}
+}
+
+bool G3PointAction::computeNormWithFlann(unsigned index,
+										 NormsTableType* theNorms,
+										 const G3PointAction::KDTree* kdTree)
+{
+	CCVector3 N;
+
+	QScopedPointer<CCCoreLib::ReferenceCloud> points(new CCCoreLib::ReferenceCloud(m_cloud));
+	if(findNearestNeighborsNanoFlann(index, points.data(), kdTree))
+	{
+		CCCoreLib::Neighbourhood neighbourhood(points.data());
+		N = *neighbourhood.getLSPlaneNormal();
+	}
+	else
+	{
+		return false;
+	}
+
+	theNorms->setValue(index, N);
+
+	return true;
+}
+
+bool G3PointAction::computeNormalsWithCloudCompare()
+{
+	unsigned pointCount = m_cloud->size();
+
+	if (!m_cloud || m_cloud->size() == 0)
+	{
+		ccLog::Error("Invalid cloud.");
+		return false;
+	}
+
+	CloudAdaptor adaptor(m_cloud);
+
+	// Build KD-tree (parameter: number of leaf nodes to inspect per query)
+
+	size_t leaf_max_size = 10;
+	nanoflann::KDTreeSingleIndexAdaptorFlags flags = nanoflann::KDTreeSingleIndexAdaptorFlags::None;
+	unsigned int n_thread_build = 0; // 0 => nanoflann automatically determines the number of threads to use
+
+	nanoflann::KDTreeSingleIndexAdaptorParams params(leaf_max_size, flags, n_thread_build);
+	QSharedPointer<KDTree> m_kdTree(new KDTree(3, adaptor, params));
+	m_kdTree->buildIndex();
+
+	// we instantiate 3D normal vectors
+	QSharedPointer<NormsTableType> theNorms(new NormsTableType);
+	QScopedPointer<NormsIndexesTableType> normsIndexes(new NormsIndexesTableType);
+	static const CCVector3 blankN(0, 0, 0);
+	if (!theNorms->resizeSafe(pointCount, true, &blankN))
+	{
+		normsIndexes->resize(0);
+		return false;
+	}
+
+	ccLog::Print("[computeNormalsWithCloudCompare]");
+#ifdef QT_DEBUG
+	//manually call the static per-point method!
+	for (unsigned index = 0; index < pointCount; ++index)
+	{
+		computeNormWithFlann(index, theNorms.data(), m_kNN, m_kdTree.data(), m_cloud);
+	}
+#else
+	std::vector<unsigned> pointsIndexes;
+	pointsIndexes.resize(pointCount);
+	for (unsigned i = 0; i < pointCount; ++i)
+	{
+		pointsIndexes[i] = i;
+	}
+	int threadCount = std::max(1, ccQtHelpers::GetMaxThreadCount() - 2);
+	ccLog::Print("[computeNormalsWithCloudCompare] parallel strategy, thread count " + QString::number(threadCount));
+	QThreadPool::globalInstance()->setMaxThreadCount(threadCount);
+	QtConcurrent::blockingMap(pointsIndexes, [=](int index){computeNormWithFlann(index, theNorms.data(), m_kdTree.data());});
+#endif
+
+	if (!m_cloud->hasNormals())
+	{
+		if (!m_cloud->resizeTheNormsTable())
+		{
+			ccLog::Error(QString("Not enough memory to compute normals on cloud '%1'").arg(m_cloud->getName()));
+			return false;
+		}
+	}
+
+	// we hide normals during process
+	m_cloud->showNormals(false);
+
+	// compress the normals
+	for (unsigned i = 0; i < theNorms->currentSize(); i++)
+	{
+		const CCVector3&   N     = theNorms->at(i);
+		const CompressedNormType nCode = ccNormalVectors::GetNormIndex(N);
+		m_cloud->setPointNormalIndex(i, nCode);
+	}
+
+	// preferred orientation
+	ccLog::Print("[computeNormalsWithCloudCompare] orient normals, PLUS_Z ");
+	ccNormalVectors::UpdateNormalOrientations(m_cloud, *m_cloud->normals(), ccNormalVectors::PLUS_Z);
+
+	return true;
+}
+
+bool G3PointAction::computeNormals()
+{
+	// if there are normals, already, propose to keep them
+	if (m_cloud->hasNormals())
+	{
+		QMessageBox msgBox;
+		msgBox.setInformativeText("Recompute normals?");
+		msgBox.setText("There are existing normals, keep them or recompute.");
+		QPushButton *keepButton = msgBox.addButton(tr("Keep"), QMessageBox::ActionRole);
+		msgBox.addButton(tr("Recompute"), QMessageBox::AcceptRole);
+		QPushButton *cancelButton = msgBox.addButton(tr("Cancel"), QMessageBox::AcceptRole);
+
+		msgBox.exec();
+
+		if (msgBox.clickedButton() == keepButton)
+		{
+			return true;
+		}
+		else if (msgBox.clickedButton() == cancelButton)
+		{
+			return false;
+		}
+	}
+
+#ifdef USE_OPEN3D_WITH_G3POINT
+	return computeNormalsWithOpen3D();
+#else
+	return computeNormalsWithCloudCompare();
+#endif
 }
 
 bool G3PointAction::queryNeighbors(ccPointCloud* cloud, ccMainAppInterface* appInterface, bool useParallelStrategy)
@@ -1879,7 +2048,7 @@ void G3PointAction::segment()
 
 	computeNodeSurfaces();
 
-	computeNormalsWithOpen3D();
+	computeNormals();
 
 	// compute the centroid
 	unsigned pointCount = m_cloud->size();
@@ -1950,7 +2119,7 @@ void G3PointAction::getBorders()
 	Eigen::ArrayXXi duplicatedLabelsInColumns(m_cloud->size(), m_kNN);
 	for (int n = 0; n < m_kNN; n++)
 	{
-		duplicatedLabelsInColumns(Eigen::all, n) = m_labels;
+		duplicatedLabelsInColumns(Eigen::placeholders::all, n) = m_labels;
 	}
 	Eigen::ArrayXXi labelsOfNeighbors(m_cloud->size(), m_kNN);
 	for (int index = 0; index < static_cast<float>(m_cloud->size()); index++)

src/GrainsAsEllipsoids.cpp

@@ -626,15 +626,15 @@ bool GrainsAsEllipsoids::directFit(const Eigen::ArrayX3d& xyz, Eigen::ArrayXd& p
 
 	Eigen::MatrixXd d(xyz.rows(), 10);
 
-	d << xyz(Eigen::all, 0).pow(2).matrix()
+	d << xyz(Eigen::placeholders::all, 0).pow(2).matrix()
-		, xyz(Eigen::all, 1).pow(2).matrix()
+		, xyz(Eigen::placeholders::all, 1).pow(2).matrix()
-		, xyz(Eigen::all, 2).pow(2).matrix()
+		, xyz(Eigen::placeholders::all, 2).pow(2).matrix()
-		, (2 * xyz(Eigen::all, 1) * xyz(Eigen::all, 2)).matrix()
+		, (2 * xyz(Eigen::placeholders::all, 1) * xyz(Eigen::placeholders::all, 2)).matrix()
-		, (2 * xyz(Eigen::all, 0) * xyz(Eigen::all, 2)).matrix()
+		, (2 * xyz(Eigen::placeholders::all, 0) * xyz(Eigen::placeholders::all, 2)).matrix()
-		, (2 * xyz(Eigen::all, 0) * xyz(Eigen::all, 1)).matrix()
+		, (2 * xyz(Eigen::placeholders::all, 0) * xyz(Eigen::placeholders::all, 1)).matrix()
-		, (2 * xyz(Eigen::all, 0)).matrix()
+		, (2 * xyz(Eigen::placeholders::all, 0)).matrix()
-		, (2 * xyz(Eigen::all, 1)).matrix()
+		, (2 * xyz(Eigen::placeholders::all, 1)).matrix()
-		, (2 * xyz(Eigen::all, 2)).matrix()
+		, (2 * xyz(Eigen::placeholders::all, 2)).matrix()
 		, Eigen::MatrixXd::Ones(xyz.rows(), 1);
 
 	Eigen::MatrixXd s = d.transpose() * d;
@@ -686,7 +686,7 @@ bool GrainsAsEllipsoids::directFit(const Eigen::ArrayX3d& xyz, Eigen::ArrayXd& p
 		{
 			if (eigenValues(k) == eigenValue)
 			{
-				v = eigensolver.eigenvectors()(Eigen::all, k).real();
+				v = eigensolver.eigenvectors()(Eigen::placeholders::all, k).real();
 				break;
 			}
 		}
@@ -698,7 +698,7 @@ bool GrainsAsEllipsoids::directFit(const Eigen::ArrayX3d& xyz, Eigen::ArrayXd& p
 		{
 			if (abs(eigenValues(k)) == eigenValue)
 			{
-				v = eigensolver.eigenvectors()(Eigen::all, k).real();
+				v = eigensolver.eigenvectors()(Eigen::placeholders::all, k).real();
 				break;
 			}
 		}
@@ -710,7 +710,7 @@ bool GrainsAsEllipsoids::directFit(const Eigen::ArrayX3d& xyz, Eigen::ArrayXd& p
 		{
 			if (eigenValues(k) == eigenValue)
 			{
-				v = eigensolver.eigenvectors()(Eigen::all, k).real();
+				v = eigensolver.eigenvectors()(Eigen::placeholders::all, k).real();
 				break;
 			}
 		}