debug of Wolman computation + not fitted grains are exported in results (used in Python and Matlab)

CMakeLists.txt

@@ -16,7 +16,7 @@ target_sources(G3PointPlugin
     ui/WolmanCustomPlot.ui
 )
 
-	set(QG3POINT_PLUGIN_VERSION "0.5")
+	set(QG3POINT_PLUGIN_VERSION "0.6")
 
 	add_subdirectory( include )
 	add_subdirectory( src )

src/G3PointAction.cpp

@@ -1135,7 +1135,7 @@ bool G3PointAction::wolman()
 		ellipsoidLabels(i) = i;
 	}
 
-	// rebuild a matrix with the coordinates the labels of the original cloud
+	// rebuild vectors with the coordinates and the labels of the original points
 	int n_points = m_grainsAsEllipsoids->m_cloud->size();
 	Eigen::ArrayXf x(n_points);
 	Eigen::ArrayXf y(n_points);
@@ -1178,13 +1178,17 @@ bool G3PointAction::wolman()
 		Eigen::ArrayXf distances;
 		Eigen::Index minLoc;
 
-		if (k % 20 == 0)
+		if (k % 20 == 0) // display progress
 		{
 			std::cout << k << std::endl;
 		}
+
 		float r0 = dist(urbg);
 		float r1 = dist(urbg);
-		// std::cout << r0 << ", " << r1 << "," << std::endl;
+// #pragma omp critical
+// 		{
+// 		std::cout << r0 << ", " << r1 << "," << std::endl;
+// 		}
 		x_grid = arange(x.minCoeff() - r0 * dx, x.maxCoeff(), dx);
 		y_grid = arange(y.minCoeff() - r1 * dx, y.maxCoeff(), dx);
 		int nx = x_grid.size();
@@ -1200,36 +1204,57 @@ bool G3PointAction::wolman()
 				iWolman(ix, iy) = minLoc;
 			}
 		}
+		// std::cout << dist.block(0, 0, 5, 5) << std::endl;
+		// std::cout << iWolman.block(0, 0, 5, 5) << std::endl;
 		XXb condition = (dist < dx / 10);
 		Eigen::ArrayXf iWolmanSelection(condition.count());
 		int indexInWolmanSelection = 0;
-		for (int k = 0; k < condition.size(); k++)
+		// for (int k = 0; k < condition.size(); k++)
+		// {
+		// 	if (condition(k))
+		// 	{
+		// 		iWolmanSelection(indexInWolmanSelection) = iWolman(k);
+		// 		indexInWolmanSelection++;
+		// 	}
+		// }
+		for (int ix = 0; ix < nx; ix++)
 		{
-			if (condition(k))
+			for (int iy = 0; iy < ny; iy++)
 			{
-				iWolmanSelection(indexInWolmanSelection) = iWolman(k);
+				if (condition(ix, iy))
-				indexInWolmanSelection++;
+				{
+					iWolmanSelection(indexInWolmanSelection) = iWolman(ix, iy);
+					indexInWolmanSelection++;
+				}
 			}
 		}
 		Eigen::ArrayXi wolmanSelection;
 		wolmanSelection = pointsLabels(iWolmanSelection);
+		// std::cout << "[" << k << "] wolmanSelection\n" << wolmanSelection << std::endl;
 		std::unordered_set<int> setOfA(wolmanSelection.begin(), wolmanSelection.end());
 		std::vector<int> y_ind;
-	// 	// A is wolmanSelection
+		for (auto label : wolmanSelection)
-	// 	// B is ellipsoidsLabels
-		for (int i = 0; i < nEllipsoids; ++i)
 		{
-			int ellipsoidLabel = ellipsoidLabels[i];
+			// some grains does not have an associated ellipdsoid, they shall not be considered in the Wolman statistics
-			if (setOfA.find(ellipsoidLabel) != setOfA.end())
+			if (m_grainsAsEllipsoids->m_fitNotOK.count(label))
 			{
-				y_ind.push_back(ellipsoidLabel);
+				// std::cout << "discard ellipsoid with label: "<< ellipsoidLabel << " (fitNotOK)" << std::endl;
+				continue;
+			}
+			else
+			{
+				y_ind.push_back(label);
 			}
 		}
 		Eigen::ArrayXf d_item(y_ind.size());
+		// std::cout << "y_ind\n";
 		for (int i = 0; i < y_ind.size(); i++)
 		{
+			// std::cout << y_ind[i] << " " ;
 			d_item(i) = b_axis(y_ind[i]) * 1000; // conversion to mm
 		}
+		// std::cout << std::endl;
+		// std::cout << "[" << k << "] " << d_item << std::endl;
 #pragma omp critical
 		{
 		d.push_back(d_item);
@@ -1247,6 +1272,8 @@ bool G3PointAction::wolman()
 		dq(i, Eigen::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)),
 					   std_dev(dq(Eigen::all, 1)),
@@ -1255,6 +1282,18 @@ bool G3PointAction::wolman()
 							quant(d_sample, 0.5),
 							quant(d_sample, 0.9)};
 
+	// std::cout << "d_sample\n" << d_sample << std::endl;
+	std::cout << "quantl 0.1 0.5 0.9 \n"
+			  << dq_final(0) << " "
+			  << dq_final(1) << " "
+			  << dq_final(2) << std::endl;
+
+	std::cout << "std_dev 0.1 0.5 0.9 \n"
+			  << edq(0) << " "
+			  << edq(1) << " "
+			  << edq(2) << std::endl;
+
+
 	showWolman(d_sample, dq_final, edq);
 
 	return true;

src/G3PointPlots.cpp

@@ -89,9 +89,9 @@ bool G3PointPlots::exportToCSV(QString filename, SharedDataContainer container,
 	if (dq_final) // WolmanCustomPlot
 	{
 		stream << "# D10 [mm], D50 [mm], D90 [mm]" << endl;
-		stream << (*dq_final)(0) << (*dq_final)(0) << (*dq_final)(0) << endl;
+		stream << (*dq_final)(0) << " " << (*dq_final)(1) << " " << (*dq_final)(2) << endl;
 		stream << "# std(D10) [mm], std(D50) [mm], std(D90) [mm]" << endl;
-		stream << (*edq)(0) << (*edq)(0) << (*edq)(0) << endl;
+		stream << (*edq)(0) << " " <<  (*edq)(1) << " " <<  (*edq)(2 ) << endl;
 		stream << "diameter [m], pdf" << endl;
 	}
 	else // AnglesCustomPlot

src/GrainsAsEllipsoids.cpp

@@ -74,11 +74,25 @@ GrainsAsEllipsoids::GrainsAsEllipsoids(ccPointCloud *cloud, ccMainAppInterface *
 	}
 
 	// remove data corresponding to stacks were the fit was not successful
-	Eigen::Array3f nullArray;
-	nullArray.fill(NAN);
 	for (auto el : m_fitNotOK)
 	{
-		m_center[el] = nullArray;
+		// if the fit is not OK, we use the centroid as a center
+		int nPoints = m_stacks[el].size();
+		Eigen::MatrixX3d points(nPoints, 3);
+		for (int index = 0; index < nPoints; index++)
+		{
+			const CCVector3* point = m_cloud->getPoint(m_stacks[el][index]);
+			points(index, 0) = point->x;
+			points(index, 1) = point->y;
+			points(index, 2) = point->z;
+		}
+		// compute the centroid of the label
+		Eigen::RowVector3d centroid = points.colwise().mean();
+		m_center[el] << centroid.x(), centroid.y(), centroid.z();
+
+		m_radii[el].fill(0);
+
+		m_rotationMatrix[el].fill(NAN);
 	}
 
 	m_ccBBoxAll.setValidity(true);
@@ -113,10 +127,10 @@ bool GrainsAsEllipsoids::exportResultsAsCloud()
 
 	for (int idx = 0; idx < m_center.size(); idx++)
 	{
-		if (m_fitNotOK.count(idx))
+		// if (m_fitNotOK.count(idx))
-		{
+		// {
-			continue;
+		// 	continue;
-		}
+		// }
 		Eigen::Vector3f center {m_center[idx].x(), m_center[idx].y(), m_center[idx].z()};
 		Eigen::Vector3f point = center;
 		CCVector3 ccPoint(point(0), point(1), point(2));
@@ -149,10 +163,10 @@ bool GrainsAsEllipsoids::exportResultsAsCloud()
 	int indexInResults = 0;
 	for (int index = 0; index < m_center.size(); index++)
 	{
-		if (m_fitNotOK.count(index)) // when the fit was not successful, the point is not exported
+		// if (m_fitNotOK.count(index)) // when the fit was not successful, the point is not exported
-		{
+		// {
-			continue;
+		// 	continue;
-		}
+		// }
 		sf->setValue(indexInResults, index);
 		indexInResults++;
 	}
@@ -172,10 +186,10 @@ bool GrainsAsEllipsoids::exportResultsAsCloud()
 	CCCoreLib::ScalarField* sfRadiusZ = cloud->getScalarField(sfIdxRadiusZ);
 	for (unsigned int index = 0; index < cloud->size(); index++)
 	{
-		if (m_fitNotOK.count(index))
+		// if (m_fitNotOK.count(index))
-		{
+		// {
-			continue;
+		// 	continue;
-		}
+		// }
 		sfRadiusX->setValue(index, m_radii[index].x());
 		sfRadiusY->setValue(index, m_radii[index].y());
 		sfRadiusZ->setValue(index, m_radii[index].z());
@@ -213,10 +227,10 @@ bool GrainsAsEllipsoids::exportResultsAsCloud()
 	CCCoreLib::ScalarField* sfR22 = cloud->getScalarField(sfIdxR22);
 	for (unsigned int index = 0; index < cloud->size(); index++)
 	{
-		if (m_fitNotOK.count(index))
+		// if (m_fitNotOK.count(index))
-		{
+		// {
-			continue;
+		// 	continue;
-		}
+		// }
 		sfR00->setValue(index, m_rotationMatrix[index](0, 0));
 		sfR01->setValue(index, m_rotationMatrix[index](0, 1));
 		sfR02->setValue(index, m_rotationMatrix[index](0, 2));
@@ -241,8 +255,6 @@ bool GrainsAsEllipsoids::exportResultsAsCloud()
 	cloud->showColors(true);
 	cloud->setPointSize(9);
 
-	// m_cloud->getParent()->addChild(cloud, ccHObject::DP_PARENT_OF_OTHER, 0);
-	// m_app->addToDB(cloud);
 	m_cloud->addChild(cloud);
 	m_app->addToDB(cloud);
 
@@ -1250,7 +1262,7 @@ bool GrainsAsEllipsoids::fromFile_MeOnly(QFile& in, short dataVersion, int flags
 	{
 		float maxRadius = m_radii[idx].maxCoeff();
 		CCVector3 center(m_center[idx](0), m_center[idx](1), m_center[idx](2));
-		if (center.x != center.x)
+		if (m_radii[idx].x() != -1) // all radii are equal to zero when the fit was not successful
 		{
 			m_fitNotOK.insert(idx);
 			continue;

src/WolmanCustomPlot.cpp

@@ -1,6 +1,8 @@
 #include "WolmanCustomPlot.h"
 #include "ui_WolmanCustomPlot.h"
 
+#include <iostream>
+
 WolmanCustomPlot::WolmanCustomPlot(const Eigen::ArrayXf &d_sample, const Eigen::Array3d& dq_final, const Eigen::Array3d& edq):
 	m_dq_final(dq_final),
 	m_edq(edq),
@@ -15,14 +17,30 @@ WolmanCustomPlot::WolmanCustomPlot(const Eigen::ArrayXf &d_sample, const Eigen::
 	QPen pen;
 
 	m_graph = this->addGraph();
-	QVector<double> x_data(d_sample.size());
+
-	QVector<double> y_data(d_sample.size());
+	// build x_data
-	for (int k = 0; k < d_sample.size(); k++)
+	int nSamples = d_sample.size();
+	QVector<double> x_data(nSamples);
+	for (int k = 0; k < nSamples; k++)
 	{
 		x_data[k] = d_sample(k);
-		y_data[k] = (static_cast<double>(k)) / static_cast<double>(d_sample.size());
 	}
-	std::sort(x_data.begin(), x_data.end());
+	while(x_data.contains(0.))
+	{
+		int index = x_data.indexOf(0.);
+		std::cout << "[WolmanCustomPlot::WolmanCustomPlot] remove null diameter at index " << QString::number(index).toStdString() << std::endl;
+		x_data.remove(index);
+	}
+	std::sort(x_data.begin(), x_data.end()); // sort diameters
+
+	// build y_data
+	int nValidSamples = x_data.size();
+	QVector<double> y_data(nValidSamples);
+	for (int k = 0; k < nValidSamples; k++)
+	{
+		y_data[k] = (static_cast<double>(k)) / static_cast<double>(nValidSamples);
+	}
+
 	m_graph->setData(x_data, y_data);
 	m_graph->rescaleAxes();
 	// give the axes some labels: