diff --git a/CMakeLists.txt b/CMakeLists.txt index db1a142..6a7822e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -24,6 +24,7 @@ if ( PLUGIN_G3POINT ) C:/Users/PaulLeroy/miniconda3/envs/env_4_CloudCompare/include C:/opt/open3d-devel-windows-amd64-0.18.0/include C:/opt/GeometricTools/GTE + C:/opt/boost_1_77_0 ) # may be needed for debug diff --git a/include/ActionA.h b/include/ActionA.h index 82ac36b..c336834 100644 --- a/include/ActionA.h +++ b/include/ActionA.h @@ -37,6 +37,7 @@ public: void fit(); void exportResults(); bool wolman(); + bool angles(); bool processNewStacks(std::vector>& newStacks, int pointCount); bool buildStacksFromG3PointLabelSF(CCCoreLib::ScalarField *g3PointLabel); bool merge(XXb& condition); diff --git a/include/G3PointDialog.h b/include/G3PointDialog.h index dceb809..d17f40b 100644 --- a/include/G3PointDialog.h +++ b/include/G3PointDialog.h @@ -29,6 +29,7 @@ public: void emitFit(){emit fit();} void emitExportResults(){emit exportResults();} void emitWolman(){emit wolman();} + void emitAngles(){emit angles();} void emitTransparencyChanged(double transparency){emit transparencyChanged(transparency);} void emitDrawSurfaces(bool state){emit drawSurfaces(state);} void emitDrawLines(bool state){emit drawLines(state);} @@ -60,6 +61,7 @@ signals: void fit(); void exportResults(); void wolman(); + void angles(); void allClicked(bool state); void onlyOneClicked(bool state); diff --git a/src/ActionA.cpp b/src/ActionA.cpp index 11526ae..174d6ff 100644 --- a/src/ActionA.cpp +++ b/src/ActionA.cpp @@ -1084,6 +1084,28 @@ double std_dev(const T &vec) return std::sqrt((vec - vec.mean()).square().sum() / (vec.size() - 1)); } +void showWolman(const Eigen::ArrayXf& d_sample) +{ + // QCustomPlot + WolmanCustomPlot* wolmanCustomPlot = new WolmanCustomPlot(); + QCPGraph* graph = wolmanCustomPlot->addGraph(); + QVector x_data(d_sample.size()); + QVector y_data(d_sample.size()); + for (int k = 0; k < d_sample.size(); 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()); + graph->setData(x_data, y_data); + graph->rescaleAxes(); + // give the axes some labels: + wolmanCustomPlot->xAxis->setScaleType(QCPAxis::stLogarithmic); + wolmanCustomPlot->xAxis->setLabel("Diameter [mm]"); + wolmanCustomPlot->yAxis->setLabel("CDF"); + wolmanCustomPlot->show(); +} + bool G3PointAction::wolman() { int n_iter = 10; @@ -1219,24 +1241,147 @@ bool G3PointAction::wolman() quant(d_sample, 0.5), quant(d_sample, 0.9)}; - // QCustomPlot - WolmanCustomPlot* wolmanCustomPlot = new WolmanCustomPlot(); - QCPGraph* graph = wolmanCustomPlot->addGraph(); - QVector x_data(d_sample.size()); - QVector y_data(d_sample.size()); - for (int k = 0; k < d_sample.size(); k++) + showWolman(d_sample); + + return true; +} + +//! Default number of classes for associated histogram +static const unsigned MAX_HISTOGRAM_SIZE = 512; + +bool computeHistogram(const QVector& data, QVector& axis, QVector& histogram) +{ + double minData = *std::min_element(data.begin(), data.end()); + double maxData = *std::max_element(data.begin(), data.end()); + double range = maxData - minData; + + if (range == 0 || data.size() == 0) { - x_data[k] = d_sample(k); - y_data[k] = (static_cast(k)) / static_cast(d_sample.size()); + //can't build histogram of a flat field + return false; } - std::sort(x_data.begin(), x_data.end()); - graph->setData(x_data, y_data); - graph->rescaleAxes(); - // give the axes some labels: - wolmanCustomPlot->xAxis->setScaleType(QCPAxis::stLogarithmic); - wolmanCustomPlot->xAxis->setLabel("Diameter [mm]"); - wolmanCustomPlot->yAxis->setLabel("CDF"); - wolmanCustomPlot->show(); + else + { + unsigned count = data.size(); + unsigned numberOfBins = static_cast(ceil(sqrt(static_cast(count)))); + // numberOfBins = std::max(std::min(numberOfBins, MAX_HISTOGRAM_SIZE), 4); + numberOfBins = 10; + + axis.resize(numberOfBins); + for (int i = 0; i < numberOfBins; i++) + { + axis[i] = minData + i * range / numberOfBins; + } + + //reserve memory + try + { + histogram.resize(numberOfBins); + } + catch (const std::bad_alloc&) + { + ccLog::Warning("[computeHistogram] Failed to allocate histogram!"); + } + + std::fill(histogram.begin(), histogram.end(), 0); + + //compute histogram + ScalarType step = static_cast(numberOfBins) / range; + for (unsigned i = 0; i < count; ++i) + { + const ScalarType& val = data[i]; + + if (CCCoreLib::ScalarField::ValidValue(val)) + { + unsigned bin = static_cast((val - minData) * step); + ++histogram[std::min(bin, numberOfBins - 1)]; + } + } + } + + return true; +} + +void showHistogram(const QVector& data) +{ + // QCustomPlot + WolmanCustomPlot* anglesCustomPlot = new WolmanCustomPlot(); + QVector axis; + QVector histogram; + computeHistogram(data, axis, histogram); + QCPBars *regen = new QCPBars(anglesCustomPlot->xAxis, anglesCustomPlot->yAxis); + regen->setData(axis, histogram); + regen->rescaleAxes(); + regen->setPen(QPen(QColor(0, 168, 140).lighter(130))); + regen->setBrush(QColor(0, 168, 140)); + anglesCustomPlot->xAxis->setLabel("Azimut [°]"); + anglesCustomPlot->yAxis->setLabel("Counts"); + anglesCustomPlot->show(); +} + +bool G3PointAction::angles() +{ + if (m_grainsAsEllipsoids.isNull()) + { + ccLog::Error("[G3PointAction::angles] no ellipsoid, not possible to do angles analysis"); + return false; + } + + float delta = 1e32; + int n_ellipsoids = m_grainsAsEllipsoids->m_rotationMatrix.size(); + QVector granuloAngleMView(n_ellipsoids); + QVector granuloAngleXView(n_ellipsoids); + + for (int i = 0; i < n_ellipsoids; i++) + { + float u, v, w; + + Eigen::Vector3f p2 {m_grainsAsEllipsoids->m_rotationMatrix[i](0, 0), + m_grainsAsEllipsoids->m_rotationMatrix[i](1, 0), + m_grainsAsEllipsoids->m_rotationMatrix[i](2, 0)}; + + // x-y plot - mapview (angle with y axis) + Eigen::Vector3f p1 {m_grainsAsEllipsoids->m_center[i].x(), + m_grainsAsEllipsoids->m_center[i].y() + delta, + m_grainsAsEllipsoids->m_center[i].z()}; + float angle = atan2(p1.cross(p2).norm(), p1.dot(p2)); + u = p2(0); + v = p2(1); + if ((angle > M_PI / 2) || (angle < - M_PI / 2)) + { + u = -u; + v = -v; + } + granuloAngleMView[i] = (atan(v / u) + M_PI / 2) * 180 / M_PI; + + // x-z plot + p1 << m_grainsAsEllipsoids->m_center[i].x(), + m_grainsAsEllipsoids->m_center[i].y(), + m_grainsAsEllipsoids->m_center[i].z() + delta; + angle = atan2(p1.cross(p2).norm(), p1.dot(p2)); + v = p2(0); + w = p2(1); + if ((angle > M_PI / 2) || (angle < - M_PI / 2)) + { + v = -v; + w = -w; + } + granuloAngleXView[i] = (atan(v / w) + M_PI / 2) * 180 / M_PI; + } + + // QCustomPlot + WolmanCustomPlot* anglesCustomPlot = new WolmanCustomPlot(); + QVector axis; + QVector histogram; + computeHistogram(granuloAngleMView, axis, histogram); + QCPBars *regen = new QCPBars(anglesCustomPlot->xAxis, anglesCustomPlot->yAxis); + regen->setData(axis, histogram); + regen->rescaleAxes(); + regen->setPen(QPen(QColor(0, 168, 140).lighter(130))); + regen->setBrush(QColor(0, 168, 140)); + anglesCustomPlot->xAxis->setLabel("Azimut [°]"); + anglesCustomPlot->yAxis->setLabel("Counts"); + anglesCustomPlot->show(); return true; } @@ -1872,6 +2017,7 @@ void G3PointAction::showDlg() connect(m_dlg, &G3PointDialog::fit, s_g3PointAction.get(), &G3Point::G3PointAction::fit); connect(m_dlg, &G3PointDialog::exportResults, s_g3PointAction.get(), &G3Point::G3PointAction::exportResults); connect(m_dlg, &G3PointDialog::wolman, s_g3PointAction.get(), &G3Point::G3PointAction::wolman); + connect(m_dlg, &G3PointDialog::angles, s_g3PointAction.get(), &G3Point::G3PointAction::angles); connect(m_dlg, &QDialog::finished, s_g3PointAction.get(), &G3Point::G3PointAction::clean); connect(m_dlg, &QDialog::finished, s_g3PointAction.get(), &G3Point::G3PointAction::resetDlg); // dialog is defined with Qt::WA_DeleteOnClose diff --git a/src/G3PointDialog.cpp b/src/G3PointDialog.cpp index 9dff6e4..ac42445 100644 --- a/src/G3PointDialog.cpp +++ b/src/G3PointDialog.cpp @@ -26,6 +26,7 @@ G3PointDialog::G3PointDialog(QString cloudName, QWidget *parent) connect(this->ui->pushButtonFit, &QPushButton::clicked, this, &G3PointDialog::emitFit); connect(this->ui->pushButtonExportResults, &QPushButton::clicked, this, &G3PointDialog::emitExportResults); connect(this->ui->pushButtonWolman, &QPushButton::clicked, this, &G3PointDialog::emitWolman); + connect(this->ui->pushButtonAngles, &QPushButton::clicked, this, &G3PointDialog::emitAngles); connect(this->ui->checkBoxSurfaces, &QCheckBox::clicked, this, &G3PointDialog::emitDrawSurfaces); connect(this->ui->checkBoxWireframes, &QCheckBox::clicked, this, &::G3PointDialog::emitDrawLines); connect(this->ui->checkBoxPoints, &QCheckBox::clicked, this, &G3PointDialog::emitDrawPoints); diff --git a/src/GrainsAsEllipsoids.cpp b/src/GrainsAsEllipsoids.cpp index c532bf5..3a2d432 100644 --- a/src/GrainsAsEllipsoids.cpp +++ b/src/GrainsAsEllipsoids.cpp @@ -432,6 +432,15 @@ bool GrainsAsEllipsoids::explicitToImplicit(const Eigen::Array3f& center, const Eigen::Matrix3f& rotationMatrix, Eigen::ArrayXd& parameters) { + // INSPIRED BY MATLAB CODE + + // Cast ellipsoid defined with explicit parameters to implicit vector form. + // + // Examples: + // p = ellipse_ex2im([xc,yc,zc],[xr,yr,zr],eye(3,3)); + + // Matlab code => Copyright 2011 Levente Hunyadi + float xrr = 1 / radii(0); float yrr = 1 / radii(1); float zrr = 1 / radii(2); @@ -506,6 +515,29 @@ bool GrainsAsEllipsoids::implicitToExplicit(const Eigen::ArrayXd& parameters, Eigen::Array3f& radii, Eigen::Matrix3f& rotationMatrix) { + // INSPIRED BY MATLAB CODE + + // Cast ellipsoid defined with implicit parameter vector to explicit form. + // The implicit equation of a general ellipse is + // F(x,y,z) = Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz - 1 = 0 + // + // Input arguments: + // v: + // the 10 parameters describing the ellipsoid algebraically + // Output arguments: + // center: + // ellispoid center coordinates [cx; cy; cz] + // ax: + // ellipsoid semi-axes (radii) [a; b; c] + // quat: NOT IN THIS CPP VERSION, ONLY MATLAB VERSION + // ellipsoid rotation in quaternion representation + // R: + // ellipsoid rotation (radii directions as rows of the 3x3 matrix) + // + // See also: ellipse_im2ex + + // Matlab code => Copyright 2011 Levente Hunyadi + Eigen::ArrayXd p = parameters; p(3) = 0.5 * p(3); @@ -554,6 +586,26 @@ bool GrainsAsEllipsoids::implicitToExplicit(const Eigen::ArrayXd& parameters, bool GrainsAsEllipsoids::directFit(const Eigen::ArrayX3d& xyz, Eigen::ArrayXd& parameters) { + // INSPIRED BY MATLAB CODE + + // Direct least squares fitting of ellipsoids under the constraint 4J - I^2 > 0. + // The constraint confines the class of ellipsoids to fit to those whose smallest radius + // is at least half of the largest radius. + // + // Input arguments: + // x,y,z; + // x, y and z coodinates of 3D points + // + // Output arguments: + // p: + // a 10-parameter vector of the algebraic ellipsoid fit + // + // References: + // Qingde Li and John G. Griffiths, "Least Squares Ellipsoid Specific Fitting", + // Proceedings of the Geometric Modeling and Processing, 2004. + + // Matlab code reference => Copyright 2011 Levente Hunyadi + Eigen::MatrixXd d(xyz.rows(), 10); d << xyz(Eigen::all, 0).pow(2).matrix() diff --git a/ui/qG3PointDialog.ui b/ui/qG3PointDialog.ui index 0d125cf..8a729fc 100644 --- a/ui/qG3PointDialog.ui +++ b/ui/qG3PointDialog.ui @@ -17,7 +17,7 @@ - 0 + 1 @@ -360,6 +360,13 @@ + + + + Angles + + +