saberbud/qG3Point
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

histograms + code commenting

Browse Source
This commit is contained in:
Paul Leroy
2025-02-24 17:31:45 +01:00
parent 9c8b34f4c9
commit 0f8bd0a329
7 changed files with 227 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CMakeLists.txt
+1
View File
@@ -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
include/ActionA.h
+1
View File
@@ -37,6 +37,7 @@ public:
void fit();
void exportResults();
bool wolman();
bool angles();
bool processNewStacks(std::vector<std::vector<int>>& newStacks, int pointCount);
bool buildStacksFromG3PointLabelSF(CCCoreLib::ScalarField *g3PointLabel);
bool merge(XXb& condition);
include/G3PointDialog.h
+2
View File
@@ -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);
src/ActionA.cpp
+162 -16
View File
@@ -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<double> x_data(d_sample.size());
QVector<double> 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<double>(k)) / static_cast<double>(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<double> x_data(d_sample.size());
QVector<double> 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<double>& data, QVector<double>& axis, QVector<double>& 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<double>(k)) / static_cast<double>(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<unsigned>(ceil(sqrt(static_cast<double>(count))));
// numberOfBins = std::max<unsigned>(std::min<unsigned>(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<ScalarType>(numberOfBins) / range;
for (unsigned i = 0; i < count; ++i)
{
const ScalarType& val = data[i];
if (CCCoreLib::ScalarField::ValidValue(val))
{
unsigned bin = static_cast<unsigned>((val - minData) * step);
++histogram[std::min(bin, numberOfBins - 1)];
}
}
}
return true;
}
void showHistogram(const QVector<double>& data)
{
// QCustomPlot
WolmanCustomPlot* anglesCustomPlot = new WolmanCustomPlot();
QVector<double> axis;
QVector<double> 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<double> granuloAngleMView(n_ellipsoids);
QVector<double> 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<double> axis;
QVector<double> 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
src/G3PointDialog.cpp
+1
View File
@@ -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);
src/GrainsAsEllipsoids.cpp
+52
View File
@@ -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()
ui/qG3PointDialog.ui
+8 -1
View File
@@ -17,7 +17,7 @@
<item>
<widget class="QTabWidget" name="tabWidget">
<property name="currentIndex">
<number>0</number>
<number>1</number>
</property>
<widget class="QWidget" name="tab">
<attribute name="title">
@@ -360,6 +360,13 @@
</property>
</widget>
</item>
<item>
<widget class="QPushButton" name="pushButtonAngles">
<property name="text">
<string>Angles</string>
</property>
</widget>
</item>
<item>
<spacer name="verticalSpacer">
<property name="orientation">