diff --git a/src/hole_filling.cpp b/src/hole_filling.cpp
index d522d16..1aa4402 100644
--- a/src/hole_filling.cpp
+++ b/src/hole_filling.cpp
@@ -2,6 +2,7 @@
 #include "IO.h"
 #include "util.h"
 #include <MeshReconstruction.h>
+#include <OpenMesh/Core/Utils/PropertyManager.hh>
 
 
 static std::vector<std::vector<HalfedgeHandle>> findHoles(MyMesh &mesh) {
@@ -153,7 +154,7 @@ vector<MyMesh::VertexHandle> Hole_Filling::next_neighbors(const vector<MyMesh::V
 
 // ***** Computation of RBF
 
-pair<pair<Eigen::MatrixXd &, Eigen::VectorXd &>, vector<MyMesh::Point> &> Hole_Filling::compute_approx_mat(vector<MyMesh::VertexHandle> vlist)
+pair<pair<Eigen::MatrixXd &, Eigen::VectorXd &>, vector<MyMesh::Point> &> Hole_Filling::compute_approx_mat(vector<MyMesh::VertexHandle> vlist, double normal_scale)
 {
     const int n(vlist.size()), d(10) ;
     Eigen::MatrixXd & A = *(new Eigen::MatrixXd(3*n+d,3*n+d)) ;
@@ -170,12 +171,12 @@ pair<pair<Eigen::MatrixXd &, Eigen::VectorXd &>, vector<MyMesh::Point> &> Hole_F
     //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))) ;
+        pts_list.push_back(_mesh.point(vlist.at(i)) + _mesh.normal(vlist.at(i)) * normal_scale) ;
     }
     //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))) ;
+        pts_list.push_back(_mesh.point(vlist.at(i)) - _mesh.normal(vlist.at(i)) * normal_scale) ;
     }
 
     int nn = pts_list.size() ;
@@ -302,15 +303,37 @@ Mesh Hole_Filling::poly_n_out(const Implicit_RBF &implicit, Rect3 domain)
 }
 
 
+/* Computes a mesh's bounding box and stores it in a mesh property
+ * named "bounding_box". */
+static void computeMeshBoundingBox(MyMesh &mesh, Hole_Filling &hf) {
+	try {
+		auto mesh_bb = OpenMesh::getProperty<void, Rect3>
+			(mesh, "bounding_box");
+	} catch (const std::runtime_error &e) {
+		auto mesh_bb = OpenMesh::getOrMakeProperty<void, Rect3>
+			(mesh, "bounding_box");
+		std::vector<VertexHandle> verts;
+		for (VertexHandle vh : mesh.vertices()) {
+			verts.push_back(vh);
+		}
+		*mesh_bb = hf.estimate_BB(verts);
+	}
+}
+
+
 MyMesh fillHoleImplicit(MyMesh &mesh, Hole_Filling &hf,
 						std::vector<HalfedgeHandle> &hole) {
+	computeMeshBoundingBox(mesh, hf);
+	Rect3 mesh_bb = *OpenMesh::getProperty<void, Rect3>(mesh, "bounding_box");
+	double diag = mesh_bb.size.Norm();
+
 	std::vector<VertexHandle> verts;
 	for (HalfedgeHandle hh : hole) {
 		verts.push_back(mesh.to_vertex_handle(hh));
 	}
 	auto bb = hf.estimate_BB(verts) ;
 	verts = hf.next_neighbors(verts);
-	auto [system, pts_list] = hf.compute_approx_mat(verts);
+	auto [system, pts_list] = hf.compute_approx_mat(verts, diag * .1);
 	auto [alpha, beta] = hf.solve_approx(system, pts_list.size(), 10);
 	Implicit_RBF rbf(alpha, beta, pts_list);
 
diff --git a/src/hole_filling.h b/src/hole_filling.h
index d636c69..939054e 100644
--- a/src/hole_filling.h
+++ b/src/hole_filling.h
@@ -103,7 +103,7 @@ public:
     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<pair<Eigen::MatrixXd &,Eigen::VectorXd &>,vector<MyMesh::Point> &> compute_approx_mat(vector<MyMesh::VertexHandle> vlist, double normal_scale=1) ;
     pair<vector<float>&, vector<float>&> solve_approx(const pair<Eigen::MatrixXd &, Eigen::VectorXd &> &p, int n, int d) ;
 
     // IO