From 54a0f44b81b0d6af2ad0b4aa5a746537b9331ecc Mon Sep 17 00:00:00 2001 From: wmayer Date: Sat, 6 Aug 2022 16:08:20 +0200 Subject: [PATCH] Mesh: implement smoothing based on median filter --- src/Mod/Mesh/App/Core/Elements.h | 6 ++ src/Mod/Mesh/App/Core/Smoothing.cpp | 132 +++++++++++++++++++++++++++- src/Mod/Mesh/App/Core/Smoothing.h | 33 ++++++- src/Mod/Mesh/App/MeshPyImp.cpp | 16 +++- src/Mod/Mesh/Gui/DlgSmoothing.cpp | 10 +++ src/Mod/Mesh/Gui/DlgSmoothing.h | 3 +- 6 files changed, 191 insertions(+), 9 deletions(-) diff --git a/src/Mod/Mesh/App/Core/Elements.h b/src/Mod/Mesh/App/Core/Elements.h index b3c5c5aebf..ffddf36376 100644 --- a/src/Mod/Mesh/App/Core/Elements.h +++ b/src/Mod/Mesh/App/Core/Elements.h @@ -328,6 +328,12 @@ public: */ bool HasNeighbour (unsigned short usSide) const { return (_aulNeighbours[usSide] != FACET_INDEX_MAX); } + /** + * Checks if the given index is a neighbour facet. + */ + bool IsNeighbour(FacetIndex index) const { + return Side(index) < 3; + } /** Counts the number of edges without neighbour. */ inline unsigned short CountOpenEdges() const; /** Returns true if there is an edge without neighbour, otherwise false. */ diff --git a/src/Mod/Mesh/App/Core/Smoothing.cpp b/src/Mod/Mesh/App/Core/Smoothing.cpp index e0345cd9f1..05a0fa5dc2 100644 --- a/src/Mod/Mesh/App/Core/Smoothing.cpp +++ b/src/Mod/Mesh/App/Core/Smoothing.cpp @@ -31,6 +31,7 @@ #include "Elements.h" #include "Iterator.h" #include "Approximation.h" +#include using namespace MeshCore; @@ -38,7 +39,6 @@ using namespace MeshCore; AbstractSmoothing::AbstractSmoothing(MeshKernel& m) : kernel(m) - , tolerance(0) , component(Normal) , continuity(C0) { @@ -56,6 +56,7 @@ void AbstractSmoothing::initialize(Component comp, Continuity cont) PlaneFitSmoothing::PlaneFitSmoothing(MeshKernel& m) : AbstractSmoothing(m) + , maximum(FLT_MAX) { } @@ -102,7 +103,7 @@ void PlaneFitSmoothing::Smooth(unsigned int iterations) N.Scale(-1.0, -1.0, -1.0); // maximum value to move is distance to mean plane - float d = std::min(fabs(this->tolerance),fabs(N*L)); + float d = std::min(fabs(this->maximum),fabs(N*L)); N.Scale(d,d,d); PointArray[v_it.Position()].Set(v_it->x - N.x, v_it->y - N.y, v_it->z - N.z); @@ -156,7 +157,7 @@ void PlaneFitSmoothing::SmoothPoints(unsigned int iterations, const std::vector< N.Scale(-1.0, -1.0, -1.0); // maximum value to move is distance to mean plane - float d = std::min(fabs(this->tolerance),fabs(N*L)); + float d = std::min(fabs(this->maximum),fabs(N*L)); N.Scale(d,d,d); PointArray[v_it.Position()].Set(v_it->x - N.x, v_it->y - N.y, v_it->z - N.z); @@ -304,3 +305,128 @@ void TaubinSmoothing::SmoothPoints(unsigned int iterations, const std::vector; +inline Base::Vector3d find_median(std::vector &container) +{ + auto compare_angle_normal = [](const AngleNormal& an1, const AngleNormal& an2) { + return an1.first < an2.first; + }; + size_t n = container.size() / 2; + std::nth_element(container.begin(), container.begin() + n, container.end(), compare_angle_normal); + + if ((container.size() % 2) == 1) { + return container[n].second; + } + else { + // even sized vector -> average the two middle values + auto max_it = std::max_element(container.begin(), container.begin() + n, compare_angle_normal); + Base::Vector3d vec = (max_it->second + container[n].second) / 2.0; + vec.Normalize(); + return vec; + } +} +} + +MedianFilterSmoothing::MedianFilterSmoothing(MeshKernel& m) + : AbstractSmoothing(m), weights(1) +{ +} + +MedianFilterSmoothing::~MedianFilterSmoothing() +{ +} + +void MedianFilterSmoothing::Smooth(unsigned int iterations) +{ + std::vector point_indices(kernel.CountPoints()); + std::generate(point_indices.begin(), point_indices.end(), Base::iotaGen(0)); + MeshCore::MeshRefFacetToFacets ff_it(kernel); + MeshCore::MeshRefPointToFacets vf_it(kernel); + + for (unsigned int i=0; i& point_indices) +{ + MeshCore::MeshRefFacetToFacets ff_it(kernel); + MeshCore::MeshRefPointToFacets vf_it(kernel); + + for (unsigned int i=0; i& point_indices) +{ + const MeshCore::MeshPointArray& points = kernel.GetPoints(); + const MeshCore::MeshFacetArray& facets = kernel.GetFacets(); + + // Initialize the array with the real normals + std::vector faceNormals; + faceNormals.reserve(facets.size()); + MeshCore::MeshFacetIterator iter(kernel); + for (iter.Init(); iter.More(); iter.Next()) { + faceNormals.emplace_back(Base::toVector(iter->GetNormal())); + } + + // Step 1: determine face normals + for (FacetIndex pos = 0; pos < facets.size(); pos++) { + iter.Set(pos); + Base::Vector3d refNormal = Base::toVector(iter->GetNormal()); + const std::set& cv = ff_it[pos]; + const MeshCore::MeshFacet& facet = facets[pos]; + + std::vector anglesWithFaces; + for (auto fi : cv) { + iter.Set(fi); + Base::Vector3d faceNormal = Base::toVector(iter->GetNormal()); + double angle = refNormal.GetAngle(faceNormal); + + int absWeight = std::abs(weights); + if (absWeight > 1 && facet.IsNeighbour(fi)) { + if (weights < 0) { + angle = -angle; + } + for (int i = 0; i < absWeight; i++) { + anglesWithFaces.emplace_back(angle, faceNormal); + } + } + else { + anglesWithFaces.emplace_back(angle, faceNormal); + } + } + + faceNormals[pos] = find_median(anglesWithFaces); + } + + // Step 2: move vertices + for (auto pos : point_indices) { + Base::Vector3d P = Base::toVector(points[pos]); + const std::set& cv = vf_it[pos]; + + double totalArea = 0.0; + Base::Vector3d totalvT; + for (auto it : cv) { + iter.Set(it); + + double faceArea = iter->Area(); + totalArea += faceArea; + + Base::Vector3d C = Base::toVector(iter->GetGravityPoint()); + + Base::Vector3d PC = C - P; + Base::Vector3d mT = faceNormals[it]; + Base::Vector3d vT = (PC * mT) * mT; + totalvT += vT * faceArea; + } + + P = P + totalvT / totalArea; + kernel.SetPoint(pos, Base::toVector(P)); + } +} diff --git a/src/Mod/Mesh/App/Core/Smoothing.h b/src/Mod/Mesh/App/Core/Smoothing.h index 300d97b1d5..ed39aa0126 100644 --- a/src/Mod/Mesh/App/Core/Smoothing.h +++ b/src/Mod/Mesh/App/Core/Smoothing.h @@ -32,6 +32,7 @@ namespace MeshCore class MeshKernel; class MeshRefPointToPoints; class MeshRefPointToFacets; +class MeshRefFacetToFacets; /** Base class for smoothing algorithms. */ class MeshExport AbstractSmoothing @@ -60,7 +61,6 @@ public: protected: MeshKernel& kernel; - float tolerance; Component component; Continuity continuity; }; @@ -70,8 +70,14 @@ class MeshExport PlaneFitSmoothing : public AbstractSmoothing public: PlaneFitSmoothing(MeshKernel&); virtual ~PlaneFitSmoothing(); + void SetMaximum(float max) { + maximum = max; + } void Smooth(unsigned int); void SmoothPoints(unsigned int, const std::vector&); + +private: + float maximum; }; class MeshExport LaplaceSmoothing : public AbstractSmoothing @@ -107,6 +113,31 @@ protected: double micro; }; +/*! + * \brief The MedianFilterSmoothing class + * Smoothing based on median filter from the paper: + * Mesh Median Filter for Smoothing 3-D Polygonal Surfaces + */ +class MeshExport MedianFilterSmoothing : public AbstractSmoothing +{ +public: + MedianFilterSmoothing(MeshKernel&); + ~MedianFilterSmoothing() override; + void SetWeight(int w) { + weights = w; + } + void Smooth(unsigned int) override; + void SmoothPoints(unsigned int, const std::vector&) override; + +private: + void UpdatePoints(const MeshRefFacetToFacets&, + const MeshRefPointToFacets&, + const std::vector&); + +private: + int weights; +}; + } // namespace MeshCore diff --git a/src/Mod/Mesh/App/MeshPyImp.cpp b/src/Mod/Mesh/App/MeshPyImp.cpp index 85dff1f3a3..56f52997e7 100644 --- a/src/Mod/Mesh/App/MeshPyImp.cpp +++ b/src/Mod/Mesh/App/MeshPyImp.cpp @@ -1751,15 +1751,17 @@ PyObject* MeshPy::trimByPlane(PyObject *args) Py_Return; } -PyObject* MeshPy::smooth(PyObject *args, PyObject *kwds) +PyObject* MeshPy::smooth(PyObject *args, PyObject *kwds) { char* method = "Laplace"; int iter=1; double lambda = 0; double micro = 0; - static char* keywords_smooth[] = {"Method","Iteration","Lambda","Micro",nullptr}; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "|sidd",keywords_smooth, - &method, &iter, &lambda, µ)) + double maximum = 1000; + int weight = 1; + static char* keywords_smooth[] = {"Method", "Iteration", "Lambda", "Micro", "Maximum", "Weight", nullptr}; + if (!PyArg_ParseTupleAndKeywords(args, kwds, "|sidddi",keywords_smooth, + &method, &iter, &lambda, µ, &maximum, &weight)) return nullptr; PY_TRY { @@ -1781,6 +1783,12 @@ PyObject* MeshPy::smooth(PyObject *args, PyObject *kwds) } else if (strcmp(method, "PlaneFit") == 0) { MeshCore::PlaneFitSmoothing smooth(kernel); + smooth.SetMaximum(maximum); + smooth.Smooth(iter); + } + else if (strcmp(method, "MedianFilter") == 0) { + MeshCore::MedianFilterSmoothing smooth(kernel); + smooth.SetWeight(weight); smooth.Smooth(iter); } else { diff --git a/src/Mod/Mesh/Gui/DlgSmoothing.cpp b/src/Mod/Mesh/Gui/DlgSmoothing.cpp index ee89b4db3c..afcd142420 100644 --- a/src/Mod/Mesh/Gui/DlgSmoothing.cpp +++ b/src/Mod/Mesh/Gui/DlgSmoothing.cpp @@ -212,6 +212,16 @@ bool TaskSmoothing::accept() s.Smooth(widget->iterations()); } } break; + case MeshGui::DlgSmoothing::MedianFilter: + { + MeshCore::MedianFilterSmoothing s(mm->getKernel()); + if (widget->smoothSelection()) { + s.SmoothPoints(widget->iterations(), selection); + } + else { + s.Smooth(widget->iterations()); + } + } break; default: break; } diff --git a/src/Mod/Mesh/Gui/DlgSmoothing.h b/src/Mod/Mesh/Gui/DlgSmoothing.h index 2b64d6e88e..06eb1d37c0 100644 --- a/src/Mod/Mesh/Gui/DlgSmoothing.h +++ b/src/Mod/Mesh/Gui/DlgSmoothing.h @@ -45,7 +45,8 @@ public: enum Smooth { None, Taubin, - Laplace + Laplace, + MedianFilter }; DlgSmoothing(QWidget* parent = nullptr);