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 <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));
+    }
+}
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<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
 
 
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, &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 {
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);