diff --git a/CMakeLists.txt b/CMakeLists.txt index 409531e..77f51cc 100644 --- a/CMakeLists.txt +++ b/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 ) diff --git a/src/G3PointAction.cpp b/src/G3PointAction.cpp index 19751dc..c9f6b83 100644 --- a/src/G3PointAction.cpp +++ b/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); - indexInWolmanSelection++; + if (condition(ix, iy)) + { + iWolmanSelection(indexInWolmanSelection) = iWolman(ix, iy); + indexInWolmanSelection++; + } } } Eigen::ArrayXi wolmanSelection; wolmanSelection = pointsLabels(iWolmanSelection); + // std::cout << "[" << k << "] wolmanSelection\n" << wolmanSelection << std::endl; std::unordered_set setOfA(wolmanSelection.begin(), wolmanSelection.end()); std::vector y_ind; - // // A is wolmanSelection - // // B is ellipsoidsLabels - for (int i = 0; i < nEllipsoids; ++i) + for (auto label : wolmanSelection) { - int ellipsoidLabel = ellipsoidLabels[i]; - if (setOfA.find(ellipsoidLabel) != setOfA.end()) + // some grains does not have an associated ellipdsoid, they shall not be considered in the Wolman statistics + 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; diff --git a/src/G3PointPlots.cpp b/src/G3PointPlots.cpp index bd627d4..fc7f807 100644 --- a/src/G3PointPlots.cpp +++ b/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 diff --git a/src/GrainsAsEllipsoids.cpp b/src/GrainsAsEllipsoids.cpp index 91ff19f..8ce1dc1 100644 --- a/src/GrainsAsEllipsoids.cpp +++ b/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)) - { - continue; - } + // if (m_fitNotOK.count(idx)) + // { + // 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 - { - continue; - } + // if (m_fitNotOK.count(index)) // when the fit was not successful, the point is not exported + // { + // 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)) - { - continue; - } + // if (m_fitNotOK.count(index)) + // { + // 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)) - { - continue; - } + // if (m_fitNotOK.count(index)) + // { + // 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; diff --git a/src/WolmanCustomPlot.cpp b/src/WolmanCustomPlot.cpp index 49b18d4..efebd63 100644 --- a/src/WolmanCustomPlot.cpp +++ b/src/WolmanCustomPlot.cpp @@ -1,6 +1,8 @@ #include "WolmanCustomPlot.h" #include "ui_WolmanCustomPlot.h" +#include + 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 x_data(d_sample.size()); - QVector y_data(d_sample.size()); - for (int k = 0; k < d_sample.size(); k++) + + // build x_data + int nSamples = d_sample.size(); + QVector x_data(nSamples); + for (int k = 0; k < nSamples; k++) { x_data[k] = d_sample(k); - y_data[k] = (static_cast(k)) / static_cast(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 y_data(nValidSamples); + for (int k = 0; k < nValidSamples; k++) + { + y_data[k] = (static_cast(k)) / static_cast(nValidSamples); + } + m_graph->setData(x_data, y_data); m_graph->rescaleAxes(); // give the axes some labels: