diff --git a/Ponca/Fitting b/Ponca/Fitting
index 4699cd34e..45b9cfcbd 100644
--- a/Ponca/Fitting
+++ b/Ponca/Fitting
@@ -31,6 +31,7 @@
 #include "src/Fitting/mongePatch.h"
 #include "src/Fitting/sphereFit.h"
 #include "src/Fitting/orientedSphereFit.h"
+#include "src/Fitting/rimlsPlaneFit.h"
 #ifndef __CUDACC__
 # include "src/Fitting/unorientedSphereFit.h"
 #endif
diff --git a/Ponca/src/Fitting/rimlsPlaneFit.h b/Ponca/src/Fitting/rimlsPlaneFit.h
new file mode 100644
index 000000000..31a36a738
--- /dev/null
+++ b/Ponca/src/Fitting/rimlsPlaneFit.h
@@ -0,0 +1,89 @@
+/*
+ This Source Code Form is subject to the terms of the Mozilla Public
+ License, v. 2.0. If a copy of the MPL was not distributed with this
+ file, You can obtain one at http://mozilla.org/MPL/2.0/.
+*/
+
+
+#pragma once
+
+#include "./defines.h"
+namespace Ponca
+{
+
+/*!
+    \brief Plane fitting procedure on oriented point sets using robust implicit mean local surfaces
+
+    Method published in \cite oztireli:2009:feature
+
+    A plane is fitted by iteratively minimizing the following equation
+
+    \f$ f^k(x) = argmin_{S_0} \sum ( S_0 + ( x_i - x )^T n_i )^2 Φ _i(x) w(r_i^{k-1} ) 
+    w_n( Δ n_i^k )\f$
+
+    with the residuals \f$ r_i^{k-1} = f^{k-1}(x) - ( x - x_i )^T n_i \f$
+
+    \f$w\f$ and \f$w_n\f$ are Gaussian functions
+
+    Φ is an arbitrary weighting function
+
+    \inherit Concept::FittingProcedureConcept
+
+    \see Plane
+
+ */
+template<class DataPoint, class _WFunctor, typename T>
+class RimlsPlaneFitImpl : public T
+{
+    PONCA_FITTING_DECLARE_DEFAULT_TYPES
+
+protected:
+    enum { Check = Base::PROVIDES_PLANE };
+
+protected:
+    /*!< \brief A parameter which controls the sharpness of the result, smaller value lead to sharper result */
+    Scalar m_sigmaN;
+
+    /*!< \brief convergence threshold ( default = \f$ 10^{-3} \f$ ) */
+    Scalar m_minConvergence;
+
+    /*!< \brief maximum number of iteration  ( default = \f$ 10 \f$ ) */
+    unsigned int m_maxIteration;
+
+    /*!< \brief neighbors are accumulated in a vector */
+    std::vector<DataPoint> m_neighbors;
+
+public:
+    PONCA_MULTIARCH inline RimlsPlaneFitImpl() : Base() {}
+
+    PONCA_FITTING_DECLARE_INIT_ADD_FINALIZE
+
+    PONCA_MULTIARCH inline const Scalar sigmaN() const { return m_sigmaN; }
+    PONCA_MULTIARCH inline Scalar& sigmaN() { return m_sigmaN; }
+
+    PONCA_MULTIARCH inline const Scalar minConvergence() const { return m_minConvergence; }
+    PONCA_MULTIARCH inline Scalar& minConvergence() { return m_minConvergence; }
+
+    PONCA_MULTIARCH inline const unsigned int maxIteration() const { return m_maxIteration; }
+    PONCA_MULTIARCH unsigned int& maxIteration() { return m_maxIteration; }
+
+private:
+    PONCA_MULTIARCH inline Scalar gaussian(Scalar x, Scalar sigma)
+    {
+        PONCA_MULTIARCH_STD_MATH(exp);
+        return exp(-(x) / (sigma * sigma));
+    }
+
+}; // class RimlsPlaneFitImpl
+
+/// @brief Helper alias for rimls plane fitting on points using RimlsPlaneFittingImpl
+//! [RimlsPlaneFit Definition]
+    template< class DataPoint, class _WFunctor, typename T >
+    using RimlsPlaneFit = 
+        RimlsPlaneFitImpl< DataPoint, _WFunctor, 
+        Plane< DataPoint, _WFunctor, T > >;
+//! [RimlsPlaneFit Definition]
+
+#include "rimlsPlaneFit.hpp"
+} //namespace Ponca
+
diff --git a/Ponca/src/Fitting/rimlsPlaneFit.hpp b/Ponca/src/Fitting/rimlsPlaneFit.hpp
new file mode 100644
index 000000000..5ea5fd596
--- /dev/null
+++ b/Ponca/src/Fitting/rimlsPlaneFit.hpp
@@ -0,0 +1,106 @@
+/*
+ This Source Code Form is subject to the terms of the Mozilla Public
+ License, v. 2.0. If a copy of the MPL was not distributed with this
+ file, You can obtain one at http://mozilla.org/MPL/2.0/.
+*/
+
+
+template<class DataPoint, class _WFunctor, typename T>
+void RimlsPlaneFitImpl<DataPoint, _WFunctor, T>::init(const VectorType& _evalPos)
+{
+    Base::init( _evalPos );
+
+    // Setup internal values
+    m_minConvergence = Scalar(1e-3);
+    m_maxIteration = 10;
+    m_sigmaN = Scalar( 0.5 );
+}
+
+template<class DataPoint, class _WFunctor, typename T>
+bool RimlsPlaneFitImpl<DataPoint, _WFunctor, T>::addLocalNeighbor(Scalar w,
+                                                                    const VectorType &localQ,
+                                                                    const DataPoint &attributes)
+{
+
+    //if valid, the neighbor is added to the neighbors vector
+    if( Base::addLocalNeighbor( w, localQ, attributes ) ){
+        m_neighbors.push_back( attributes ); 
+        return true;
+    }
+    return false;
+}
+
+template<class DataPoint, class _WFunctor, typename T>
+FIT_RESULT RimlsPlaneFitImpl<DataPoint, _WFunctor, T>::finalize()
+{
+
+    // handle UNDEFINED cases
+    if( Base::finalize() != STABLE ) {
+        return Base::m_eCurrentState = UNDEFINED;
+    }
+
+    // handle conflict error
+    if( Base::plane().isValid() ){
+        return Base::m_eCurrentState = CONFLICT_ERROR_FOUND;
+    }
+
+    unsigned int iteration{ 0 };
+    Scalar convergence = Scalar(1);
+
+    Scalar f = Scalar(0.);
+    VectorType gradF = VectorType::Zero();
+    VectorType prevGrad = VectorType::Zero();
+    Scalar alpha = Scalar(1.);
+
+    VectorType px;
+    Scalar fx;
+    Scalar w;
+    VectorType gradW;
+
+    while( iteration < m_maxIteration && convergence > m_minConvergence ) {
+
+        // initialise values
+        Scalar sumF = Scalar(0);
+        Scalar sumW = Scalar(0);
+        VectorType sumGF = VectorType::Zero();
+        VectorType sumGW = VectorType::Zero();
+        VectorType sumN = VectorType::Zero();
+
+        for(DataPoint& nei : m_neighbors) {
+
+            px = Base::m_w.evalPos() - nei.pos();
+            fx = px.dot( nei.normal() );
+
+            // use gaussian weights
+            if( iteration > 0 ) {
+                alpha = gaussian( ( nei.normal() - gradF ).squaredNorm(), m_sigmaN );
+                alpha *= gaussian( (fx - f) * (fx - f), 0.5 * Base::m_w.evalScale() );
+            }
+
+            w = alpha * Base::m_w.w( nei.pos(), nei ).first;
+            gradW = alpha * 2 * px * Base::m_w.scaledw(px, nei);
+
+            sumW += w;
+            sumGW += gradW;
+            sumF += w * fx;
+            sumGF += gradW * fx;
+            sumN += w * nei.normal();
+        }
+
+        if( sumW == 0 ){
+            return Base::m_eCurrentState = UNDEFINED;
+        }
+
+        prevGrad = gradF;
+
+        // update the actual potential and gradient
+        f = sumF / sumW;
+        gradF = ( sumGF - f * sumGW + sumN ) / sumW;
+        convergence = ( prevGrad - gradF ).squaredNorm();
+        ++iteration;
+		
+    }
+
+    Base::setPlane( -gradF, -f * gradF );
+    return Base::m_eCurrentState = STABLE;
+}
diff --git a/cmake/PoncaConfigureFitting.cmake b/cmake/PoncaConfigureFitting.cmake
index b947bbb78..a5e78d758 100644
--- a/cmake/PoncaConfigureFitting.cmake
+++ b/cmake/PoncaConfigureFitting.cmake
@@ -28,6 +28,8 @@ set(ponca_Fitting_INCLUDE
     "${PONCA_src_ROOT}/Ponca/src/Fitting/orientedSphereFit.hpp"
     "${PONCA_src_ROOT}/Ponca/src/Fitting/plane.h"
     "${PONCA_src_ROOT}/Ponca/src/Fitting/primitive.h"
+    "${PONCA_src_ROOT}/Ponca/src/Fitting/rimlsPlaneFit.h"
+    "${PONCA_src_ROOT}/Ponca/src/Fitting/rimlsPlaneFit.hpp"
     "${PONCA_src_ROOT}/Ponca/src/Fitting/sphereFit.h"
     "${PONCA_src_ROOT}/Ponca/src/Fitting/sphereFit.hpp"
     "${PONCA_src_ROOT}/Ponca/src/Fitting/unorientedSphereFit.h"
diff --git a/doc/src/ponca.bib b/doc/src/ponca.bib
index 1079b7bab..d79543c31 100644
--- a/doc/src/ponca.bib
+++ b/doc/src/ponca.bib
@@ -175,4 +175,14 @@ @article{Alexa:2009:Hermite
 articleno = {20},
 numpages = {10},
 keywords = {interpolation, Point-based modeling, Hermite data}
+}
+
+@inproceedings{oztireli:2009:feature,
+  title={Feature preserving point set surfaces based on non-linear kernel regression},
+  author={Oztireli, Cengiz and Guennebaud, Gael and Gross, Markus},
+  booktitle={Computer graphics forum},
+  volume={28},
+  number={2},
+  pages={493--501},
+  year={2009}
 }
\ No newline at end of file
diff --git a/tests/src/fit_plane.cpp b/tests/src/fit_plane.cpp
index 77dcc7ca2..423e2cab1 100644
--- a/tests/src/fit_plane.cpp
+++ b/tests/src/fit_plane.cpp
@@ -20,6 +20,7 @@
 #include <Ponca/src/Fitting/basket.h>
 #include <Ponca/src/Fitting/covariancePlaneFit.h>
 #include <Ponca/src/Fitting/meanPlaneFit.h>
+#include <Ponca/src/Fitting/rimlsPlaneFit.h>
 #include <Ponca/src/Fitting/weightFunc.h>
 #include <Ponca/src/Fitting/weightKernel.h>
 
@@ -143,6 +144,9 @@ void callSubTests()
     typedef Basket<Point, WeightSmoothFunc, MeanPlaneFit> MeanFitSmooth;
     typedef Basket<Point, WeightConstantFunc, MeanPlaneFit> MeanFitConstant;
 
+    typedef Basket<Point, WeightSmoothFunc, RimlsPlaneFit> RimlsFitSmooth;
+    typedef Basket<Point, WeightConstantFunc, RimlsPlaneFit> RimlsFitConstant;
+
     // test if conflicts are detected
     //! [Conflicting type]
     typedef Basket<Point, WeightConstantFunc, Plane,
@@ -161,6 +165,8 @@ void callSubTests()
         CALL_SUBTEST(( testFunction<Point, CovFitConstant, WeightConstantFunc, true>() ));
         CALL_SUBTEST(( testFunction<Point, MeanFitSmooth, WeightSmoothFunc, false>() ));
         CALL_SUBTEST(( testFunction<Point, MeanFitConstant, WeightConstantFunc, false>() ));
+        CALL_SUBTEST(( testFunction<Point, RimlsFitSmooth, WeightSmoothFunc, false>() ));
+        CALL_SUBTEST(( testFunction<Point, RimlsFitConstant, WeightConstantFunc, false>() ));
         // Check if fitting conflict is detected
         CALL_SUBTEST(( testFunction<Point, Hybrid1, WeightConstantFunc, false>(false, false, false, true) ));
         CALL_SUBTEST(( testFunction<Point, Hybrid2, WeightConstantFunc, false>(false, false, false, true) ));
@@ -172,6 +178,8 @@ void callSubTests()
     {
         CALL_SUBTEST(( testFunction<Point, CovFitSmooth, WeightSmoothFunc, true>(false, true, true) ));
         CALL_SUBTEST(( testFunction<Point, CovFitConstant, WeightConstantFunc, true>(false, true, true) ));
+        CALL_SUBTEST(( testFunction<Point, RimlsFitSmooth, WeightSmoothFunc, false>(false, true, true, false) ));
+        CALL_SUBTEST(( testFunction<Point, RimlsFitConstant, WeightConstantFunc, false>(false, true, true, false) ));
     }
     cout << "Ok!" << endl;
 }