Mesh: [skip ci] implement initial guess of cylinder axis
This commit is contained in:
@@ -1027,6 +1027,7 @@ CylinderFit::CylinderFit()
|
||||
: _vBase(0,0,0)
|
||||
, _vAxis(0,0,1)
|
||||
, _fRadius(0)
|
||||
, _initialGuess(false)
|
||||
{
|
||||
}
|
||||
|
||||
@@ -1034,6 +1035,21 @@ CylinderFit::~CylinderFit()
|
||||
{
|
||||
}
|
||||
|
||||
Base::Vector3f CylinderFit::GetInitialAxisFromNormals(const std::vector<Base::Vector3f>& n) const
|
||||
{
|
||||
PlaneFit planeFit;
|
||||
planeFit.AddPoints(n);
|
||||
planeFit.Fit();
|
||||
return planeFit.GetNormal();
|
||||
}
|
||||
|
||||
void CylinderFit::SetInitialValues(const Base::Vector3f& b, const Base::Vector3f& n)
|
||||
{
|
||||
_vBase = b;
|
||||
_vAxis = n;
|
||||
_initialGuess = true;
|
||||
}
|
||||
|
||||
float CylinderFit::Fit()
|
||||
{
|
||||
if (CountPoints() < 7)
|
||||
@@ -1046,15 +1062,22 @@ float CylinderFit::Fit()
|
||||
[](const Base::Vector3f& v) { return Wm4::Vector3d(v.x, v.y, v.z); });
|
||||
|
||||
Wm4::Vector3d cnt, axis;
|
||||
if (_initialGuess) {
|
||||
cnt = Base::convertTo<Wm4::Vector3d>(_vBase);
|
||||
axis = Base::convertTo<Wm4::Vector3d>(_vAxis);
|
||||
}
|
||||
|
||||
double radius, height;
|
||||
Wm4::CylinderFit3<double> fit(input.size(), input.data(), cnt, axis, radius, height, false);
|
||||
Wm4::CylinderFit3<double> fit(input.size(), input.data(), cnt, axis, radius, height, _initialGuess);
|
||||
_initialGuess = false;
|
||||
|
||||
_vBase = Base::convertTo<Base::Vector3f>(cnt);
|
||||
_vAxis = Base::convertTo<Base::Vector3f>(axis);
|
||||
_fRadius = float(radius);
|
||||
|
||||
_fLastResult = double(fit);
|
||||
|
||||
#if defined(_DEBUG)
|
||||
#if defined(FC_DEBUG)
|
||||
Base::Console().Message(" WildMagic Cylinder Fit: Base: (%0.4f, %0.4f, %0.4f), Axis: (%0.6f, %0.6f, %0.6f), Radius: %0.4f, Std Dev: %0.4f\n",
|
||||
_vBase.x, _vBase.y, _vBase.z, _vAxis.x, _vAxis.y, _vAxis.z, _fRadius, GetStdDeviation());
|
||||
#endif
|
||||
@@ -1067,7 +1090,7 @@ float CylinderFit::Fit()
|
||||
if (result < FLOAT_MAX) {
|
||||
Base::Vector3d base = cylFit.GetBase();
|
||||
Base::Vector3d dir = cylFit.GetAxis();
|
||||
#if defined(_DEBUG)
|
||||
#if defined(FC_DEBUG)
|
||||
Base::Console().Message("MeshCoreFit::Cylinder Fit: Base: (%0.4f, %0.4f, %0.4f), Axis: (%0.6f, %0.6f, %0.6f), Radius: %0.4f, Std Dev: %0.4f, Iterations: %d\n",
|
||||
base.x, base.y, base.z, dir.x, dir.y, dir.z, cylFit.GetRadius(), cylFit.GetStdDeviation(), cylFit.GetNumIterations());
|
||||
#endif
|
||||
@@ -1183,6 +1206,29 @@ float CylinderFit::GetStdDeviation() const
|
||||
return sqrt((ulPtCt / (ulPtCt - 3.0f)) * ((1.0f / ulPtCt) * fSumXi2 - fMean * fMean));
|
||||
}
|
||||
|
||||
void CylinderFit::GetBounding(Base::Vector3f& bottom, Base::Vector3f& top) const
|
||||
{
|
||||
float distMin = FLT_MAX;
|
||||
float distMax = FLT_MIN;
|
||||
|
||||
std::list<Base::Vector3f>::const_iterator cIt;
|
||||
for (cIt = _vPoints.begin(); cIt != _vPoints.end(); ++cIt) {
|
||||
float dist = cIt->DistanceToPlane(_vBase, _vAxis);
|
||||
if (dist < distMin) {
|
||||
distMin = dist;
|
||||
bottom = *cIt;
|
||||
}
|
||||
if (dist > distMax) {
|
||||
distMax = dist;
|
||||
top = *cIt;
|
||||
}
|
||||
}
|
||||
|
||||
// Project the points onto the cylinder axis
|
||||
bottom = bottom.Perpendicular(_vBase, _vAxis);
|
||||
top = top.Perpendicular(_vBase, _vAxis);
|
||||
}
|
||||
|
||||
void CylinderFit::ProjectToCylinder()
|
||||
{
|
||||
Base::Vector3f cBase(GetBase());
|
||||
|
||||
@@ -393,11 +393,16 @@ public:
|
||||
virtual ~CylinderFit();
|
||||
float GetRadius() const;
|
||||
Base::Vector3f GetBase() const;
|
||||
void SetInitialValues(const Base::Vector3f&, const Base::Vector3f&);
|
||||
/**
|
||||
* Returns the axis of the fitted cylinder. If Fit() has not been called the null vector is
|
||||
* returned.
|
||||
*/
|
||||
Base::Vector3f GetAxis() const;
|
||||
/**
|
||||
* Returns an initial axis based on point normals.
|
||||
*/
|
||||
Base::Vector3f GetInitialAxisFromNormals(const std::vector<Base::Vector3f>& n) const;
|
||||
/**
|
||||
* Fit a cylinder into the given points. If the fit fails FLOAT_MAX is returned.
|
||||
*/
|
||||
@@ -416,11 +421,17 @@ public:
|
||||
* Projects the points onto the fitted cylinder.
|
||||
*/
|
||||
void ProjectToCylinder();
|
||||
/**
|
||||
* Get the bottom and top points of the cylinder. The distance of these
|
||||
* points gives the height of the cylinder.
|
||||
*/
|
||||
void GetBounding(Base::Vector3f& bottom, Base::Vector3f& top) const;
|
||||
|
||||
protected:
|
||||
Base::Vector3f _vBase; /**< Base vector of the cylinder. */
|
||||
Base::Vector3f _vAxis; /**< Axis of the cylinder. */
|
||||
float _fRadius; /**< Radius of the cylinder. */
|
||||
bool _initialGuess;
|
||||
};
|
||||
|
||||
// -------------------------------------------------------------------------------
|
||||
|
||||
@@ -1079,6 +1079,26 @@ std::vector<Base::Vector3f> MeshKernel::CalcVertexNormals() const
|
||||
return normals;
|
||||
}
|
||||
|
||||
std::vector<Base::Vector3f> MeshKernel::GetFacetNormals(const std::vector<unsigned long>& facets) const
|
||||
{
|
||||
std::vector<Base::Vector3f> normals;
|
||||
normals.reserve(facets.size());
|
||||
|
||||
for (std::vector<unsigned long>::const_iterator it = facets.begin(); it != facets.end(); ++it) {
|
||||
const MeshFacet& face = _aclFacetArray[*it];
|
||||
|
||||
const Base::Vector3f& p1 = _aclPointArray[face._aulPoints[0]];
|
||||
const Base::Vector3f& p2 = _aclPointArray[face._aulPoints[1]];
|
||||
const Base::Vector3f& p3 = _aclPointArray[face._aulPoints[2]];
|
||||
|
||||
Base::Vector3f n = (p2 - p1) % (p3 - p1);
|
||||
n.Normalize();
|
||||
normals.emplace_back(n);
|
||||
}
|
||||
|
||||
return normals;
|
||||
}
|
||||
|
||||
// Evaluation
|
||||
float MeshKernel::GetSurface() const
|
||||
{
|
||||
|
||||
@@ -112,6 +112,7 @@ public:
|
||||
* by summarizing the normals of the associated facets.
|
||||
*/
|
||||
std::vector<Base::Vector3f> CalcVertexNormals() const;
|
||||
std::vector<Base::Vector3f> GetFacetNormals(const std::vector<unsigned long>&) const;
|
||||
|
||||
/** Returns the facet at the given index. This method is rather slow and should be
|
||||
* called occasionally only. For fast access the MeshFacetIterator interface should
|
||||
|
||||
@@ -37,6 +37,7 @@
|
||||
#include <App/Application.h>
|
||||
#include <App/Document.h>
|
||||
#include <App/DocumentObjectGroup.h>
|
||||
#include <Gui/Command.h>
|
||||
#include <Gui/SelectionObject.h>
|
||||
#include <Mod/Mesh/App/Core/Approximation.h>
|
||||
#include <Mod/Mesh/App/Core/Segmentation.h>
|
||||
@@ -54,7 +55,7 @@ public:
|
||||
virtual std::vector<float> getParameter(FitParameter::Points pts) const {
|
||||
std::vector<float> values;
|
||||
MeshCore::PlaneFit fit;
|
||||
fit.AddPoints(pts);
|
||||
fit.AddPoints(pts.points);
|
||||
if (fit.Fit() < FLOAT_MAX) {
|
||||
Base::Vector3f base = fit.GetBase();
|
||||
Base::Vector3f axis = fit.GetNormal();
|
||||
@@ -77,9 +78,20 @@ public:
|
||||
virtual std::vector<float> getParameter(FitParameter::Points pts) const {
|
||||
std::vector<float> values;
|
||||
MeshCore::CylinderFit fit;
|
||||
fit.AddPoints(pts);
|
||||
fit.AddPoints(pts.points);
|
||||
if (!pts.normals.empty()) {
|
||||
Base::Vector3f base = fit.GetGravity();
|
||||
Base::Vector3f axis = fit.GetInitialAxisFromNormals(pts.normals);
|
||||
fit.SetInitialValues(base, axis);
|
||||
|
||||
#if defined(FC_DEBUG)
|
||||
Base::Console().Message("Initial axis: (%f, %f, %f)\n", axis.x, axis.y, axis.z);
|
||||
#endif
|
||||
}
|
||||
|
||||
if (fit.Fit() < FLOAT_MAX) {
|
||||
Base::Vector3f base = fit.GetBase();
|
||||
Base::Vector3f base, top;
|
||||
fit.GetBounding(base, top);
|
||||
Base::Vector3f axis = fit.GetAxis();
|
||||
float radius = fit.GetRadius();
|
||||
values.push_back(base.x);
|
||||
@@ -89,6 +101,25 @@ public:
|
||||
values.push_back(axis.y);
|
||||
values.push_back(axis.z);
|
||||
values.push_back(radius);
|
||||
|
||||
#if defined(FC_DEBUG)
|
||||
// Only for testing purposes
|
||||
try {
|
||||
float height = Base::Distance(base, top);
|
||||
Gui::Command::doCommand(Gui::Command::App,
|
||||
"cyl = App.ActiveDocument.addObject('Part::Cylinder', 'Cylinder')\n"
|
||||
"cyl.Radius = %f\n"
|
||||
"cyl.Height = %f\n"
|
||||
"cyl.Placement = App.Placement(App.Vector(%f,%f,%f), App.Rotation(App.Vector(0,0,1), App.Vector(%f,%f,%f)))\n",
|
||||
radius, height, base.x, base.y, base.z, axis.x, axis.y, axis.z);
|
||||
|
||||
Gui::Command::doCommand(Gui::Command::App,
|
||||
"axis = cyl.Placement.Rotation.multVec(App.Vector(0,0,1))\n"
|
||||
"print('Final axis: ({}, {}, {})'.format(axis.x, axis.y, axis.z))\n");
|
||||
}
|
||||
catch (...) {
|
||||
}
|
||||
#endif
|
||||
}
|
||||
return values;
|
||||
}
|
||||
@@ -102,7 +133,7 @@ public:
|
||||
virtual std::vector<float> getParameter(FitParameter::Points pts) const {
|
||||
std::vector<float> values;
|
||||
MeshCore::SphereFit fit;
|
||||
fit.AddPoints(pts);
|
||||
fit.AddPoints(pts.points);
|
||||
if (fit.Fit() < FLOAT_MAX) {
|
||||
Base::Vector3f base = fit.GetCenter();
|
||||
float radius = fit.GetRadius();
|
||||
@@ -230,14 +261,15 @@ void ParametersDialog::on_compute_clicked()
|
||||
{
|
||||
const Mesh::MeshObject& kernel = myMesh->Mesh.getValue();
|
||||
if (kernel.hasSelectedFacets()) {
|
||||
FitParameter::Points fitpts;
|
||||
std::vector<unsigned long> facets, points;
|
||||
kernel.getFacetsFromSelection(facets);
|
||||
points = kernel.getPointsFromFacets(facets);
|
||||
MeshCore::MeshPointArray coords = kernel.getKernel().GetPoints(points);
|
||||
fitpts.normals = kernel.getKernel().GetFacetNormals(facets);
|
||||
|
||||
// Copy points into right format
|
||||
FitParameter::Points fitpts;
|
||||
fitpts.insert(fitpts.end(), coords.begin(), coords.end());
|
||||
fitpts.points.insert(fitpts.points.end(), coords.begin(), coords.end());
|
||||
coords.clear();
|
||||
|
||||
values = fitParameter->getParameter(fitpts);
|
||||
|
||||
@@ -41,7 +41,10 @@ class Ui_SegmentationBestFit;
|
||||
class FitParameter
|
||||
{
|
||||
public:
|
||||
typedef std::vector<Base::Vector3f> Points;
|
||||
struct Points {
|
||||
std::vector<Base::Vector3f> points;
|
||||
std::vector<Base::Vector3f> normals;
|
||||
};
|
||||
virtual ~FitParameter() {}
|
||||
virtual std::vector<float> getParameter(Points) const = 0;
|
||||
};
|
||||
|
||||
Reference in New Issue
Block a user