kindred/create
1
1
Fork 0
Code Issues 37 Pull Requests Actions Packages Projects 9 Releases 2 Wiki Activity

Mesh: implement smoothing based on median filter

Browse Source
This commit is contained in:
wmayer
2022-08-06 16:08:20 +02:00
parent ebda764bf1
commit 54a0f44b81
6 changed files with 191 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/Mod/Mesh/App/Core/Elements.h
View File

@@ -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. */

132
src/Mod/Mesh/App/Core/Smoothing.cpp
View File

@@ -31,6 +31,7 @@
#include "Elements.h"
#include "Iterator.h"
#include "Approximation.h"
#include <Base/Tools.h>
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<float>(fabs(this->tolerance),fabs(N*L));
float d = std::min<float>(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<float>(fabs(this->tolerance),fabs(N*L));
float d = std::min<float>(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<Po
Umbrella(vv_it, vf_it, -(lambda+micro), point_indices);
}
}
namespace {
using AngleNormal = std::pair<double, Base::Vector3d>;
inline Base::Vector3d find_median(std::vector<AngleNormal> &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<unsigned long> point_indices(kernel.CountPoints());
std::generate(point_indices.begin(), point_indices.end(), Base::iotaGen<unsigned long>(0));
MeshCore::MeshRefFacetToFacets ff_it(kernel);
MeshCore::MeshRefPointToFacets vf_it(kernel);
for (unsigned int i=0; i<iterations; i++) {
UpdatePoints(ff_it, vf_it, point_indices);
}
}
void MedianFilterSmoothing::SmoothPoints(unsigned int iterations, const std::vector<PointIndex>& point_indices)
{
MeshCore::MeshRefFacetToFacets ff_it(kernel);
MeshCore::MeshRefPointToFacets vf_it(kernel);
for (unsigned int i=0; i<iterations; i++) {
UpdatePoints(ff_it, vf_it, point_indices);
}
}
void MedianFilterSmoothing::UpdatePoints(const MeshRefFacetToFacets& ff_it,
const MeshRefPointToFacets& vf_it,
const std::vector<PointIndex>& point_indices)
{
const MeshCore::MeshPointArray& points = kernel.GetPoints();
const MeshCore::MeshFacetArray& facets = kernel.GetFacets();
// Initialize the array with the real normals
std::vector<Base::Vector3d> faceNormals;
faceNormals.reserve(facets.size());
MeshCore::MeshFacetIterator iter(kernel);
for (iter.Init(); iter.More(); iter.Next()) {
faceNormals.emplace_back(Base::toVector<double>(iter->GetNormal()));
}
// Step 1: determine face normals
for (FacetIndex pos = 0; pos < facets.size(); pos++) {
iter.Set(pos);
Base::Vector3d refNormal = Base::toVector<double>(iter->GetNormal());
const std::set<FacetIndex>& cv = ff_it[pos];
const MeshCore::MeshFacet& facet = facets[pos];
std::vector<AngleNormal> anglesWithFaces;
for (auto fi : cv) {
iter.Set(fi);
Base::Vector3d faceNormal = Base::toVector<double>(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<double>(points[pos]);
const std::set<FacetIndex>& 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<double>(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<float>(P));
}
}

33
src/Mod/Mesh/App/Core/Smoothing.h
View File

@@ -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<PointIndex>&);
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<PointIndex>&) override;
private:
void UpdatePoints(const MeshRefFacetToFacets&,
const MeshRefPointToFacets&,
const std::vector<PointIndex>&);
private:
int weights;
};
} // namespace MeshCore

16
src/Mod/Mesh/App/MeshPyImp.cpp
View File

@@ -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, &micro))
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, &micro, &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 {

10
src/Mod/Mesh/Gui/DlgSmoothing.cpp
View File

@@ -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;
}

3
src/Mod/Mesh/Gui/DlgSmoothing.h
View File

@@ -45,7 +45,8 @@ public:
enum Smooth {
None,
Taubin,
Laplace
Laplace,
MedianFilter
};
DlgSmoothing(QWidget* parent = nullptr);