diff --git a/src/hole_filling.cpp b/src/hole_filling.cpp
index 3a980e3..df641dd 100644
--- a/src/hole_filling.cpp
+++ b/src/hole_filling.cpp
@@ -54,5 +54,248 @@ void fillHolesDumb(MyMesh &mesh) {
 }
 
 
+using namespace std ;
+using namespace MeshReconstruction;
+
+/* ************** Implicit RBF ************** */
+
+Implicit_RBF::Implicit_RBF(vector<float> alpha, vector<float> beta, vector<MyMesh::Point> center) : _alpha(alpha), _beta(beta), _center(center)
+{
+    _n = _alpha.size() ;
+    _d = _beta.size() ;
+    if (_center.size() != _n)
+        throw runtime_error("Inconsistent size of alpha and centers in Implicit_RBF constructor.");
+}
+
+// Computation of the value of the implicit surface at point X
+float Implicit_RBF::val(MyMesh::Point X) const
+{
+    float res = 0 ;
+    // Computation of the sum of RBF at centers
+    for(int i=0; i<_n; i++)
+    {
+        res += _alpha.at(i) * myphi((X-_center.at(i)).norm()) ;
+    }
+    // Computation of the polynomial part
+    for(int j=0; j<_d; j++)
+    {
+        res += _beta.at(j) * myp(j, X) ;
+    }
+    return res ;
+}
+
+/* ************** Hole Filling ************** */
+
+Hole_Filling::Hole_Filling(MyMesh &mesh) : _mesh(mesh)
+{
+    _mesh.add_property(_vprop,  "vprop_flag");
+    cout << "Starting hole filling ..." << endl ;
+}
+
+// ***** Computation of boundary and its neighborhood
+
+MyMesh::HalfedgeHandle Hole_Filling::find_boundary_edge()
+{
+    MyMesh::HalfedgeIter he_it = _mesh.halfedges_begin();
+    while ( (he_it != _mesh.halfedges_end()) && (!_mesh.is_boundary(*he_it)))
+    {
+        ++he_it ;
+    }
+    if (he_it != _mesh.halfedges_end())
+        return *he_it ;
+    else
+        throw std::runtime_error("Boundary HE does not exist") ;
+}
+
+vector<MyMesh::VertexHandle> Hole_Filling::find_boundary(MyMesh::HalfedgeHandle heh)
+{
+    MyMesh::HalfedgeHandle heh_ini = heh ;
+    vector<MyMesh::VertexHandle> boundary ;
+
+    // Follow (and memorize) boundary edges
+    do
+    {
+        boundary.push_back(_mesh.to_vertex_handle(heh));
+        heh = _mesh.next_halfedge_handle(heh);
+    } while (heh != heh_ini) ;
+    return boundary ;
+}
+
+void Hole_Filling::init_mark_boundary(const vector<MyMesh::VertexHandle> & bnd)
+{
+    for (MyMesh::VertexIter v_it = _mesh.vertices_begin() ; v_it != _mesh.vertices_end(); ++v_it)
+        _mesh.property(_vprop, *v_it) = false ;
+    for (int i=0; i<bnd.size(); ++i)
+        _mesh.property(_vprop, bnd.at(i)) = true ;
+}
+
+vector<MyMesh::VertexHandle> Hole_Filling::next_neighbors(const vector<MyMesh::VertexHandle> & bnd)
+{
+    // Visit bnd vertices to find and mark next circle
+    vector<MyMesh::VertexHandle> next_bnd ;
+    for (int i=0; i<bnd.size(); i++)
+    {
+        for (MyMesh::VertexVertexIter vv_it = _mesh.vv_iter(bnd.at(i));vv_it.is_valid();++vv_it)
+        {
+            if (_mesh.property(_vprop, *vv_it) == false) // new vertex
+            {
+                _mesh.property(_vprop, *vv_it) = true ;
+                next_bnd.push_back(*vv_it);
+            }
+        }
+    }
+    return next_bnd ;
+}
+
+// ***** Computation of RBF
+
+pair<pair<Eigen::MatrixXd &, Eigen::VectorXd &>, vector<MyMesh::Point> &> Hole_Filling::compute_approx_mat(vector<MyMesh::VertexHandle> vlist)
+{
+    const int n(vlist.size()), d(10) ;
+    Eigen::MatrixXd & A = *(new Eigen::MatrixXd(3*n+d,3*n+d)) ;
+    Eigen::MatrixXd Phi(3*n,3*n);
+    Eigen::MatrixXd P(3*n,d);
+    Eigen::VectorXd & B = *(new Eigen::VectorXd(3*n+d));
+
+    vector<MyMesh::Point> & pts_list = *(new vector<MyMesh::Point>);
+    //Append vertices to pts_list
+    for (int i=0; i<n; i++)
+    {
+        pts_list.push_back(_mesh.point(vlist.at(i))) ;
+    }
+    //Append vertices+normals to pts_list
+    for (int i=0; i<n; i++)
+    {
+        pts_list.push_back(_mesh.point(vlist.at(i)) + _mesh.normal(vlist.at(i))) ;
+    }
+    //Append vertices-normals to pts_list
+    for (int i=0; i<n; i++)
+    {
+        pts_list.push_back(_mesh.point(vlist.at(i)) - _mesh.normal(vlist.at(i))) ;
+    }
+
+    int nn = pts_list.size() ;
+    // Compute corresponding B vector (0 / 1 / -1 )
+    B << Eigen::VectorXd::Zero(n), Eigen::VectorXd::Ones(n), -Eigen::VectorXd::Ones(n), Eigen::VectorXd::Zero(d) ;
+
+    // Fill Phi matrix
+    // TODO
+
+    // Fill P matrix
+    //TODO
+
+    // Set final A matrix
+    /* A = Phi | P
+     *      P' | 0    */
+
+    // TODO
+
+    cout << "size of pts_list : " << nn << endl ;
+    cout << "End computation of matrices" << endl ;
+
+    return {{A,B},pts_list} ;
+}
+
+pair<vector<float>&, vector<float>&> Hole_Filling::solve_approx(const pair<Eigen::MatrixXd &, Eigen::VectorXd &> &p, int n, int d)
+{
+    Eigen::MatrixXd & A = p.first ;
+    Eigen::VectorXd & B = p.second ;
+
+    Eigen::VectorXd res = A.householderQr().solve(B) ;
+
+    cout << "End of solver" << endl ;
+    cout << "res : " << res.head(10) << endl ;
+    vector<float> & alpha = *(new vector<float>);
+    vector<float> & beta = *(new vector<float>);
+
+    if (res.size() != (n+d))
+    {
+        cout << "taille du res : " << res.size() << endl ;
+        throw std::runtime_error("Error in solve_approx") ;
+    }
+    for (int i=0; i<n; i++)
+    {
+        alpha.push_back(res(i)) ;
+    }
+
+    for (int j=0; j<d; j++)
+    {
+        beta.push_back(res(n+j)) ;
+    }
+
+    return {alpha, beta} ;
+}
+
+// ***** Hole filling
+
+void Hole_Filling::fill_hole(string out)
+{
+	// TODO !!!
+	cout << "Au travail !" << endl ;
+}
+
+// ***** IO
+
+void Hole_Filling::colorize_prop()
+{
+    for (MyMesh::VertexIter v_it = _mesh.vertices_begin(); v_it != _mesh.vertices_end() ; ++v_it)
+    {
+        if(_mesh.property(_vprop, *v_it) == true)
+            _mesh.set_color(*v_it, MyMesh::Color(255, 0, 0)) ;
+        else
+           _mesh.set_color(*v_it, MyMesh::Color(200, 200, 200)) ;
+    }
+}
+
+void Hole_Filling::colorize_verts(const vector<MyMesh::VertexHandle> &vlist)
+{
+    for (MyMesh::VertexIter v_it = _mesh.vertices_begin(); v_it != _mesh.vertices_end() ; ++v_it)
+    {
+           _mesh.set_color(*v_it, MyMesh::Color(200, 200, 200)) ;
+    }
+
+    for(int i=0; i<vlist.size(); i++)
+    {
+        _mesh.set_color(vlist.at(i), MyMesh::Color(255, 0, 0)) ;
+    }
+}
+
+Rect3 Hole_Filling::estimate_BB(const vector<MyMesh::VertexHandle> &vlist)
+{
+    MyMesh::Point minp = _mesh.point(vlist.at(0)), maxp=minp ;
+    for (int i=0; i<vlist.size(); i++)
+    {
+        minp = min(minp, _mesh.point(vlist.at(i))) ;
+        maxp = max(maxp, _mesh.point(vlist.at(i))) ;
+    }
+    MyMesh::Point sizep(maxp-minp), centerp = (minp+maxp)/2 ;
+    minp = centerp - _scale*sizep/2 ;
+    sizep *= _scale ;
+
+    Rect3 domain ;
+    domain.min = {minp[0], minp[1], minp[2]} ;
+    domain.size = {sizep[0], sizep[1], sizep[2]} ;
+    return domain ;
+}
+
+void Hole_Filling::poly_n_out(const Implicit_RBF &implicit, Rect3 domain, string filename)
+{
+    auto implicitEq = [&implicit](Vec3 const& pos)
+    {
+        MyMesh::Point X(pos.x, pos.y, pos.z) ;
+        return implicit.val(X) ;
+    };
+
+
+    Vec3 cubeSize(domain.size*(_discr)) ;
+    cout << "mind " << domain.min << endl ;
+    cout << "sized" << domain.size << endl ;
+    cout << "cubesize "<< cubeSize << endl ;
+
+    auto mesh = MarchCube(implicitEq, domain, cubeSize);
+    WriteObjFile(mesh, filename);
+}
+
+
 void fillHolesImplicit(MyMesh &mesh) {
 }
diff --git a/src/hole_filling.h b/src/hole_filling.h
index aa213be..8e47367 100644
--- a/src/hole_filling.h
+++ b/src/hole_filling.h
@@ -3,11 +3,125 @@
 
 #include "my_mesh.h"
 #include <vector>
+#include <iostream>
+#include <cmath>
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+#include <Eigen/Dense>
+#include <MeshReconstruction.h>
+#include <IO.h>
 
 
 void fillHoleDumb(MyMesh &mesh, std::vector<HalfedgeHandle> &hole);
 void fillHolesDumb(MyMesh &mesh);
 void fillHolesImplicit(MyMesh &mesh);
 
+using namespace std ;
+using namespace MeshReconstruction;
+
+// ******************************
+// Function computing polynomial P_i (degree 2 polynomials in 3 variables)
+inline float myp(int i, MyMesh::Point X) {
+    float x = X[0] ;
+    float y = X[1] ;
+    float z = X[2] ;
+    switch (i) {
+    case 0:
+        return x*x ;
+        break;
+    case 1:
+        return y*y ;
+        break;
+    case 2:
+        return z*z ;
+        break;
+    case 3:
+        return x*y ;
+        break;
+    case 4:
+        return y*z ;
+        break;
+    case 5:
+        return x*z ;
+        break;
+    case 6:
+        return x ;
+        break;
+    case 7:
+        return y ;
+        break;
+    case 8:
+        return z;
+        break;
+    case 9:
+        return 1;
+        break;
+    default:
+        throw std::runtime_error("Error on indice i : unknown polynomial P.");
+    }
+}
+
+// ******************************
+// Function computing thin spline RBF
+inline float myphi(float r) { if (r == 0) return 0 ; else return r*r*log(r) ; }
+
+// ******************************
+// Class encoding an implicit function built from the basis of RBF
+// f(X) = Sum_i=0^n-1 alpha_i * phi(|| X-c_i ||) + Sum_j=0^9 beta_j * P_j(X)
+class Implicit_RBF {
+private:
+    int _n ;
+    int _d ;
+    vector<float> _alpha, _beta ;
+    vector<MyMesh::Point> _center ;
+
+public:
+    Implicit_RBF (vector<float> alpha, vector<float> beta, vector<MyMesh::Point> center) ;
+    float val(MyMesh::Point) const ;
+};
+
+// ******************************
+// Class encoding all the functions required for implicit hole filling
+class Hole_Filling
+{
+private:
+    MyMesh &_mesh ;
+    OpenMesh::VPropHandleT<bool> _vprop ;
+
+    // Constantes servant au rendu
+    // TODO : à mettre dans des boutons ...
+    const float _scale = 4 ; // Boîte englobante de rendu = _scale * BB des centres
+    const float _discr = 1./40. ; // BB discrétisée en 1/_discr voxels
+public:
+    Hole_Filling(MyMesh &mesh);
+    inline ~Hole_Filling() { std::cout << "Ending hole filling" << std::endl ; }
+
+    // Computation of boundary and its neighborhood
+    MyMesh::HalfedgeHandle find_boundary_edge() ;
+    vector<MyMesh::VertexHandle> find_boundary(MyMesh::HalfedgeHandle heh) ;
+    void init_mark_boundary(const vector<MyMesh::VertexHandle> & bnd) ;
+    vector<MyMesh::VertexHandle> next_neighbors(const vector<MyMesh::VertexHandle> & bnd) ;
+
+    // Computation of RBF
+    pair<pair<Eigen::MatrixXd &,Eigen::VectorXd &>,vector<MyMesh::Point> &> compute_approx_mat(vector<MyMesh::VertexHandle> vlist) ;
+    pair<vector<float>&, vector<float>&> solve_approx(const pair<Eigen::MatrixXd &, Eigen::VectorXd &> &p, int n, int d) ;
+
+    // Hole filling
+    void fill_hole(string out) ;
+
+    // IO
+    void colorize_prop() ;
+    void colorize_verts(const vector<MyMesh::VertexHandle> &vlist) ;
+    Rect3 estimate_BB(const vector<MyMesh::VertexHandle> &vlist) ;
+    void poly_n_out (const Implicit_RBF & implicit, Rect3 domain, std::string filename) ;
+};
+
+// Other ....
+
+inline MyMesh::Point min (MyMesh::Point X, MyMesh::Point Y) { MyMesh::Point P ; for (int i=0; i<3; i++) P[i] = min(X[i], Y[i]) ; return P ;}
+inline MyMesh::Point max (MyMesh::Point X, MyMesh::Point Y) { MyMesh::Point P ; for (int i=0; i<3; i++) P[i] = max(X[i], Y[i]) ; return P ;}
+
+inline ostream & operator<< (ostream & out, Vec3 v) { out << v.x << ", " << v.y << ", " << v.z ; return out ;}
+
 
 #endif
\ No newline at end of file