diff --git a/src/Mod/Mesh/App/Core/Approximation.cpp b/src/Mod/Mesh/App/Core/Approximation.cpp index 068cdf81ca..41b607810b 100644 --- a/src/Mod/Mesh/App/Core/Approximation.cpp +++ b/src/Mod/Mesh/App/Core/Approximation.cpp @@ -25,6 +25,7 @@ #ifndef _PreComp_ # include +# include #endif #include "Approximation.h" @@ -795,6 +796,120 @@ void SurfaceFit::GetCoefficients(double& a,double& b,double& c,double& d,double& f = _fCoeff[0]; } +// ------------------------------------------------------------------------------- + +CylinderFit::CylinderFit() + : _vBase(0,0,0) + , _vAxis(0,0,1) + , _fRadius(0) +{ +} + +CylinderFit::~CylinderFit() +{ +} + +float CylinderFit::Fit() +{ + if (CountPoints() < 7) + return FLOAT_MAX; + _bIsFitted = true; + + // TODO + + _fLastResult = 0; + return _fLastResult; +} + +float CylinderFit::GetRadius() const +{ + return _fRadius; +} + +Base::Vector3f CylinderFit::GetBase() const +{ + if (_bIsFitted) + return _vBase; + else + return Base::Vector3f(); +} + +Base::Vector3f CylinderFit::GetAxis() const +{ + if (_bIsFitted) + return _vAxis; + else + return Base::Vector3f(); +} + +float CylinderFit::GetDistanceToCylinder(const Base::Vector3f &rcPoint) const +{ + float fResult = FLOAT_MAX; + if (_bIsFitted) + fResult = rcPoint.DistanceToLine(_vBase, _vAxis) - _fRadius; + return fResult; +} + +float CylinderFit::GetStdDeviation() const +{ + // Mean: M=(1/N)*SUM Xi + // Variance: VAR=(N/N-3)*[(1/N)*SUM(Xi^2)-M^2] + // Standard deviation: SD=SQRT(VAR) + // Standard error of the mean: SE=SD/SQRT(N) + if (!_bIsFitted) + return FLOAT_MAX; + + float fSumXi = 0.0f, fSumXi2 = 0.0f, + fMean = 0.0f, fDist = 0.0f; + + float ulPtCt = (float)CountPoints(); + std::list< Base::Vector3f >::const_iterator cIt; + + for (cIt = _vPoints.begin(); cIt != _vPoints.end(); ++cIt) { + fDist = GetDistanceToCylinder( *cIt ); + fSumXi += fDist; + fSumXi2 += ( fDist * fDist ); + } + + fMean = (1.0f / ulPtCt) * fSumXi; + return (float)sqrt((ulPtCt / (ulPtCt - 3.0)) * ((1.0 / ulPtCt) * fSumXi2 - fMean * fMean)); +} + +void CylinderFit::ProjectToCylinder() +{ + Base::Vector3f cBase(GetBase()); + Base::Vector3f cAxis(GetAxis()); + + for (std::list< Base::Vector3f >::iterator it = _vPoints.begin(); it != _vPoints.end(); ++it) { + Base::Vector3f& cPnt = *it; + if (cPnt.DistanceToLine(cBase, cAxis) > 0) { + Base::Vector3f proj; + cBase.ProjectToPlane(cPnt, cAxis, proj); + Base::Vector3f diff = cPnt - proj; + diff.Normalize(); + cPnt = proj + diff * _fRadius; + } + else { + // Point is on the cylinder axis, so it can be moved in + // any direction perpendicular to the cylinder axis + Base::Vector3f cMov(cPnt); + do { + float x = ((float)rand() / (float)RAND_MAX); + float y = ((float)rand() / (float)RAND_MAX); + float z = ((float)rand() / (float)RAND_MAX); + cMov.Move(x,y,z); + } + while (cMov.DistanceToLine(cBase, cAxis) == 0); + + Base::Vector3f proj; + cMov.ProjectToPlane(cPnt, cAxis, proj); + Base::Vector3f diff = cPnt - proj; + diff.Normalize(); + cPnt = proj + diff * _fRadius; + } + } +} + // ----------------------------------------------------------------------------- PolynomialFit::PolynomialFit() diff --git a/src/Mod/Mesh/App/Core/Approximation.h b/src/Mod/Mesh/App/Core/Approximation.h index 272c85c900..310c08ed0e 100644 --- a/src/Mod/Mesh/App/Core/Approximation.h +++ b/src/Mod/Mesh/App/Core/Approximation.h @@ -342,7 +342,7 @@ public: /** * Destruction */ - virtual ~SurfaceFit(){}; + virtual ~SurfaceFit(){} bool GetCurvatureInfo(double x, double y, double z, double &rfCurv0, double &rfCurv1, Base::Vector3f &rkDir0, Base::Vector3f &rkDir1, double &dDistance); @@ -358,6 +358,54 @@ protected: // ------------------------------------------------------------------------------- +/** + * Approximation of a cylinder into a given set of points. + */ +class MeshExport CylinderFit : public Approximation +{ +public: + /** + * Construction + */ + CylinderFit(); + /** + * Destruction + */ + virtual ~CylinderFit(); + float GetRadius() const; + Base::Vector3f GetBase() const; + /** + * Returns the axis of the fitted cylinder. If Fit() has not been called the null vector is + * returned. + */ + Base::Vector3f GetAxis() const; + /** + * Fit a cylinder into the given points. If the fit fails FLOAT_MAX is returned. + */ + float Fit(); + /** + * Returns the distance from the point \a rcPoint to the fitted cylinder. If Fit() has not been + * called FLOAT_MAX is returned. + */ + float GetDistanceToCylinder(const Base::Vector3f &rcPoint) const; + /** + * Returns the standard deviation from the points to the fitted cylinder. If Fit() has not been + * called FLOAT_MAX is returned. + */ + float GetStdDeviation() const; + /** + * Projects the points onto the fitted cylinder. + */ + void ProjectToCylinder(); + +protected: + Base::Vector3f _vBase; /**< Base vector of the cylinder. */ + Base::Vector3f _vAxis; /**< Axis of the cylinder. */ + float _fRadius; /**< Radius of the cylinder. */ +}; + +// ------------------------------------------------------------------------------- + /** * Helper class for the quadric fit. Includes the * partial derivates of the quadric and serves for diff --git a/src/Mod/Mesh/App/Core/Segmentation.cpp b/src/Mod/Mesh/App/Core/Segmentation.cpp index e862fc32cb..785f13f624 100644 --- a/src/Mod/Mesh/App/Core/Segmentation.cpp +++ b/src/Mod/Mesh/App/Core/Segmentation.cpp @@ -35,6 +35,11 @@ void MeshSurfaceSegment::Initialize(unsigned long) { } +bool MeshSurfaceSegment::TestInitialFacet(unsigned long) const +{ + return true; +} + void MeshSurfaceSegment::AddFacet(const MeshFacet&) { } @@ -122,6 +127,11 @@ void PlaneSurfaceFit::Initialize(const MeshCore::MeshGeomFacet& tria) fitter->AddPoint(tria._aclPoints[2]); } +bool PlaneSurfaceFit::TestTriangle(const MeshGeomFacet&) const +{ + return true; +} + void PlaneSurfaceFit::AddTriangle(const MeshCore::MeshGeomFacet& tria) { fitter->AddPoint(tria.GetGravityPoint()); @@ -137,13 +147,150 @@ float PlaneSurfaceFit::Fit() return fitter->Fit(); } -float PlaneSurfaceFit::GetDistanceToSurface(const Base::Vector3f& pnt) +float PlaneSurfaceFit::GetDistanceToSurface(const Base::Vector3f& pnt) const { return fitter->GetDistanceToPlane(pnt); } // -------------------------------------------------------- +CylinderSurfaceFit::CylinderSurfaceFit() + : fitter(new CylinderFit) +{ + axis.Set(0,0,0); + radius = FLOAT_MAX; +} + +/*! + * \brief CylinderSurfaceFit::CylinderSurfaceFit + * Set a pre-defined cylinder. Internal cylinder fits are not done, then. + */ +CylinderSurfaceFit::CylinderSurfaceFit(const Base::Vector3f& b, const Base::Vector3f& a, float r) + : basepoint(b) + , axis(a) + , radius(r) + , fitter(nullptr) +{ +} + +CylinderSurfaceFit::~CylinderSurfaceFit() +{ + delete fitter; +} + +void CylinderSurfaceFit::Initialize(const MeshCore::MeshGeomFacet& tria) +{ + if (fitter) { + fitter->Clear(); + fitter->AddPoint(tria._aclPoints[0]); + fitter->AddPoint(tria._aclPoints[1]); + fitter->AddPoint(tria._aclPoints[2]); + } +} + +void CylinderSurfaceFit::AddTriangle(const MeshCore::MeshGeomFacet& tria) +{ + if (fitter) { + fitter->AddPoint(tria._aclPoints[0]); + fitter->AddPoint(tria._aclPoints[1]); + fitter->AddPoint(tria._aclPoints[2]); + } +} + +bool CylinderSurfaceFit::TestTriangle(const MeshGeomFacet& tria) const +{ + float dot = axis.Dot(tria.GetNormal()); + return fabs(dot) < 0.01f; +} + +bool CylinderSurfaceFit::Done() const +{ + if (fitter) { + return fitter->Done(); + } + + return true; +} + +float CylinderSurfaceFit::Fit() +{ + if (!fitter) + return 0; + + float fit = fitter->Fit(); + if (fit < FLOAT_MAX) { + basepoint = fitter->GetBase(); + axis = fitter->GetAxis(); + radius = fitter->GetRadius(); + } + return fit; +} + +float CylinderSurfaceFit::GetDistanceToSurface(const Base::Vector3f& pnt) const +{ + if (fitter && !fitter->Done()) { + // collect some points + return 0; + } + float dist = pnt.DistanceToLine(basepoint, axis); + return (dist - radius); +} + +// -------------------------------------------------------- + +SphereSurfaceFit::SphereSurfaceFit() +{ + center.Set(0,0,0); + radius = FLOAT_MAX; +} + +SphereSurfaceFit::SphereSurfaceFit(const Base::Vector3f& c, float r) + : center(c) + , radius(r) +{ + +} + +SphereSurfaceFit::~SphereSurfaceFit() +{ +} + +void SphereSurfaceFit::Initialize(const MeshCore::MeshGeomFacet& tria) +{ + //FIXME +} + +void SphereSurfaceFit::AddTriangle(const MeshCore::MeshGeomFacet& tria) +{ + //FIXME +} + +bool SphereSurfaceFit::TestTriangle(const MeshGeomFacet& tria) const +{ + // Already handled by GetDistanceToSurface + return true; +} + +bool SphereSurfaceFit::Done() const +{ + //FIXME + return true; +} + +float SphereSurfaceFit::Fit() +{ + //FIXME + return 0; +} + +float SphereSurfaceFit::GetDistanceToSurface(const Base::Vector3f& pnt) const +{ + float dist = Base::Distance(pnt, center); + return (dist - radius); +} + +// -------------------------------------------------------- + MeshDistanceGenericSurfaceSegment::MeshDistanceGenericSurfaceSegment(AbstractSurfaceFit* fit, const MeshKernel& mesh, unsigned long minFacets, @@ -164,6 +311,16 @@ void MeshDistanceGenericSurfaceSegment::Initialize(unsigned long index) fitter->Initialize(triangle); } +bool MeshDistanceGenericSurfaceSegment::TestInitialFacet(unsigned long index) const +{ + MeshGeomFacet triangle = kernel.GetFacet(index); + for (int i=0; i<3; i++) { + if (fabs(fitter->GetDistanceToSurface(triangle._aclPoints[i])) > tolerance) + return false; + } + return fitter->TestTriangle(triangle); +} + bool MeshDistanceGenericSurfaceSegment::TestFacet (const MeshFacet& face) const { if (!fitter->Done()) @@ -174,7 +331,7 @@ bool MeshDistanceGenericSurfaceSegment::TestFacet (const MeshFacet& face) const return false; } - return true; + return fitter->TestTriangle(triangle); } void MeshDistanceGenericSurfaceSegment::AddFacet(const MeshFacet& face) @@ -297,13 +454,14 @@ void MeshSegmentAlgorithm::FindSegments(std::vector& segm) while (startFacet != ULONG_MAX) { // collect all facets of the same geometry std::vector indices; - indices.push_back(startFacet); (*it)->Initialize(startFacet); + if ((*it)->TestInitialFacet(startFacet)) + indices.push_back(startFacet); MeshSurfaceVisitor pv(**it, indices); myKernel.VisitNeighbourFacets(pv, startFacet); // add or discard the segment - if (indices.size() == 1) { + if (indices.size() <= 1) { resetVisited.push_back(startFacet); } else { diff --git a/src/Mod/Mesh/App/Core/Segmentation.h b/src/Mod/Mesh/App/Core/Segmentation.h index 33a749ec0a..1f124ea35e 100644 --- a/src/Mod/Mesh/App/Core/Segmentation.h +++ b/src/Mod/Mesh/App/Core/Segmentation.h @@ -31,6 +31,7 @@ namespace MeshCore { class PlaneFit; +class CylinderFit; class MeshFacet; typedef std::vector MeshSegment; @@ -43,6 +44,7 @@ public: virtual bool TestFacet (const MeshFacet &rclFacet) const = 0; virtual const char* GetType() const = 0; virtual void Initialize(unsigned long); + virtual bool TestInitialFacet(unsigned long) const; virtual void AddFacet(const MeshFacet& rclFacet); void AddSegment(const std::vector&); const std::vector& GetSegments() const { return segments; } @@ -88,10 +90,11 @@ public: AbstractSurfaceFit(){} virtual ~AbstractSurfaceFit(){} virtual void Initialize(const MeshGeomFacet&) = 0; + virtual bool TestTriangle(const MeshGeomFacet&) const = 0; virtual void AddTriangle(const MeshGeomFacet&) = 0; virtual bool Done() const = 0; virtual float Fit() = 0; - virtual float GetDistanceToSurface(const Base::Vector3f&) = 0; + virtual float GetDistanceToSurface(const Base::Vector3f&) const = 0; }; class MeshExport PlaneSurfaceFit : public AbstractSurfaceFit @@ -100,10 +103,11 @@ public: PlaneSurfaceFit(); ~PlaneSurfaceFit(); void Initialize(const MeshGeomFacet&); + bool TestTriangle(const MeshGeomFacet&) const; void AddTriangle(const MeshGeomFacet&); bool Done() const; float Fit(); - float GetDistanceToSurface(const Base::Vector3f&); + float GetDistanceToSurface(const Base::Vector3f&) const; private: Base::Vector3f basepoint; @@ -111,6 +115,44 @@ private: PlaneFit* fitter; }; +class MeshExport CylinderSurfaceFit : public AbstractSurfaceFit +{ +public: + CylinderSurfaceFit(); + CylinderSurfaceFit(const Base::Vector3f& b, const Base::Vector3f& a, float r); + ~CylinderSurfaceFit(); + void Initialize(const MeshGeomFacet&); + bool TestTriangle(const MeshGeomFacet&) const; + void AddTriangle(const MeshGeomFacet&); + bool Done() const; + float Fit(); + float GetDistanceToSurface(const Base::Vector3f&) const; + +private: + Base::Vector3f basepoint; + Base::Vector3f axis; + float radius; + CylinderFit* fitter; +}; + +class MeshExport SphereSurfaceFit : public AbstractSurfaceFit +{ +public: + SphereSurfaceFit(); + SphereSurfaceFit(const Base::Vector3f& c, float r); + ~SphereSurfaceFit(); + void Initialize(const MeshGeomFacet&); + bool TestTriangle(const MeshGeomFacet&) const; + void AddTriangle(const MeshGeomFacet&); + bool Done() const; + float Fit(); + float GetDistanceToSurface(const Base::Vector3f&) const; + +private: + Base::Vector3f center; + float radius; +}; + class MeshExport MeshDistanceGenericSurfaceSegment : public MeshDistanceSurfaceSegment { public: @@ -120,6 +162,7 @@ public: bool TestFacet (const MeshFacet& rclFacet) const; const char* GetType() const { return "Generic"; } void Initialize(unsigned long); + bool TestInitialFacet(unsigned long) const; void AddFacet(const MeshFacet& rclFacet); protected: diff --git a/src/Mod/Mesh/App/Mesh.cpp b/src/Mod/Mesh/App/Mesh.cpp index 56abad9663..71c5062933 100644 --- a/src/Mod/Mesh/App/Mesh.cpp +++ b/src/Mod/Mesh/App/Mesh.cpp @@ -1699,8 +1699,8 @@ MeshObject* MeshObject::meshFromSegment(const std::vector& indice return new MeshObject(kernel, _Mtrx); } -std::vector MeshObject::getSegmentsFromType(MeshObject::GeometryType type, - float dev, unsigned long minFacets) const +std::vector MeshObject::getSegmentsOfType(MeshObject::GeometryType type, + float dev, unsigned long minFacets) const { std::vector segm; if (this->_kernel.CountFacets() == 0) @@ -1714,10 +1714,13 @@ std::vector MeshObject::getSegmentsFromType(MeshObject::GeometryType ty surf.reset(new MeshCore::MeshDistanceGenericSurfaceSegment(new MeshCore::PlaneSurfaceFit, this->_kernel, minFacets, dev)); break; - // todo! case CYLINDER: + surf.reset(new MeshCore::MeshDistanceGenericSurfaceSegment(new MeshCore::CylinderSurfaceFit, + this->_kernel, minFacets, dev)); break; case SPHERE: + surf.reset(new MeshCore::MeshDistanceGenericSurfaceSegment(new MeshCore::SphereSurfaceFit, + this->_kernel, minFacets, dev)); break; default: break; diff --git a/src/Mod/Mesh/App/Mesh.h b/src/Mod/Mesh/App/Mesh.h index 229af20efb..3418899914 100644 --- a/src/Mod/Mesh/App/Mesh.h +++ b/src/Mod/Mesh/App/Mesh.h @@ -297,7 +297,7 @@ public: const Segment& getSegment(unsigned long) const; Segment& getSegment(unsigned long); MeshObject* meshFromSegment(const std::vector&) const; - std::vector getSegmentsFromType(GeometryType, float dev, unsigned long minFacets) const; + std::vector getSegmentsOfType(GeometryType, float dev, unsigned long minFacets) const; //@} /** @name Primitives */ diff --git a/src/Mod/Mesh/App/MeshPy.xml b/src/Mod/Mesh/App/MeshPy.xml index ed549875dd..34d1f256ac 100644 --- a/src/Mod/Mesh/App/MeshPy.xml +++ b/src/Mod/Mesh/App/MeshPy.xml @@ -452,7 +452,14 @@ In the worst case each triangle can be regarded as single plane if none of its neighours is coplanar. - + + + getSegmentsOfType(type, dev,[min faces=0]) -> list +Get all segments of type. +Type can be Plane, Cylinder or Sphere + + + getSegmentsByCurvature(list) -> list The argument list gives a list if tuples where it defines the preferred maximum curvature, diff --git a/src/Mod/Mesh/App/MeshPyImp.cpp b/src/Mod/Mesh/App/MeshPyImp.cpp index 0c7c10c4a4..1ec738f450 100644 --- a/src/Mod/Mesh/App/MeshPyImp.cpp +++ b/src/Mod/Mesh/App/MeshPyImp.cpp @@ -1759,7 +1759,7 @@ PyObject* MeshPy::getPlanarSegments(PyObject *args) return NULL; Mesh::MeshObject* mesh = getMeshObjectPtr(); - std::vector segments = mesh->getSegmentsFromType + std::vector segments = mesh->getSegmentsOfType (Mesh::MeshObject::PLANE, dev, minFacets); Py::List s; @@ -1779,6 +1779,50 @@ PyObject* MeshPy::getPlanarSegments(PyObject *args) return Py::new_reference_to(s); } +PyObject* MeshPy::getSegmentsOfType(PyObject *args) +{ + char* type; + float dev; + unsigned long minFacets=0; + if (!PyArg_ParseTuple(args, "sf|k",&type,&dev,&minFacets)) + return NULL; + + Mesh::MeshObject::GeometryType geoType; + if (strcmp(type, "Plane") == 0) { + geoType = Mesh::MeshObject::PLANE; + } + else if (strcmp(type, "Cylinder") == 0) { + geoType = Mesh::MeshObject::CYLINDER; + } + else if (strcmp(type, "Sphere") == 0) { + geoType = Mesh::MeshObject::SPHERE; + } + else { + PyErr_SetString(PyExc_ValueError, "Unsupported surface type"); + return nullptr; + } + + Mesh::MeshObject* mesh = getMeshObjectPtr(); + std::vector segments = mesh->getSegmentsOfType + (geoType, dev, minFacets); + + Py::List s; + for (std::vector::iterator it = segments.begin(); it != segments.end(); ++it) { + const std::vector& segm = it->getIndices(); + Py::List ary; + for (std::vector::const_iterator jt = segm.begin(); jt != segm.end(); ++jt) { +#if PY_MAJOR_VERSION >= 3 + ary.append(Py::Long((int)*jt)); +#else + ary.append(Py::Int((int)*jt)); +#endif + } + s.append(ary); + } + + return Py::new_reference_to(s); +} + PyObject* MeshPy::getSegmentsByCurvature(PyObject *args) { PyObject* l;