diff --git a/src/Mod/Part/App/BSplineSurfacePy.xml b/src/Mod/Part/App/BSplineSurfacePy.xml
index 413e94247..6bfebdb44 100644
--- a/src/Mod/Part/App/BSplineSurfacePy.xml
+++ b/src/Mod/Part/App/BSplineSurfacePy.xml
@@ -677,5 +677,10 @@
 				<UserDocu>Returns a reparametrized copy of this surface</UserDocu>
 			</Documentation>
 		</Methode>
+		<Methode Name="interpolate">
+			<Documentation>
+				<UserDocu>Replaces this B-Spline surface by interpolating a set of points.</UserDocu>
+			</Documentation>
+		</Methode>
 	</PythonExport>
 </GenerateModel>
diff --git a/src/Mod/Part/App/BSplineSurfacePyImp.cpp b/src/Mod/Part/App/BSplineSurfacePyImp.cpp
index c016f277a..f642662f4 100644
--- a/src/Mod/Part/App/BSplineSurfacePyImp.cpp
+++ b/src/Mod/Part/App/BSplineSurfacePyImp.cpp
@@ -31,6 +31,8 @@
 # include <TColStd_Array1OfInteger.hxx>
 # include <TColgp_Array1OfPnt.hxx>
 # include <TColgp_Array2OfPnt.hxx>
+# include <Precision.hxx>
+# include <GeomAPI_PointsToBSplineSurface.hxx>
 #endif
 
 #include <Base/GeometryPyCXX.h>
@@ -1242,6 +1244,55 @@ PyObject* BSplineSurfacePy::reparametrize(PyObject * args)
     }
 }
 
+PyObject* BSplineSurfacePy::interpolate(PyObject *args)
+{
+    PyObject* obj;
+    double tol3d = Precision::Approximation();
+    PyObject* closed = Py_False;
+    PyObject* t1=0; PyObject* t2=0;
+    if (!PyArg_ParseTuple(args, "O!",&(PyList_Type), &obj))
+        return 0;
+    try {
+        Py::List list(obj);
+        Standard_Integer lu = list.size();
+        Py::List col(list.getItem(0));
+        Standard_Integer lv = col.size();
+        TColgp_Array2OfPnt interpolationPoints(1, lu, 1, lv);
+
+        Standard_Integer index1 = 0;
+        Standard_Integer index2 = 0;
+        for (Py::List::iterator it1 = list.begin(); it1 != list.end(); ++it1) {
+            index1++;
+            index2=0;
+            Py::List row(*it1);
+            for (Py::List::iterator it2 = row.begin(); it2 != row.end(); ++it2) {
+                index2++;
+                Py::Vector v(*it2);
+                Base::Vector3d pnt = v.toVector();
+                gp_Pnt newPoint(pnt.x,pnt.y,pnt.z);
+                interpolationPoints.SetValue(index1, index2, newPoint);
+            }
+        }
+
+        if (interpolationPoints.RowLength() < 2 || interpolationPoints.ColLength() < 2) {
+            Standard_Failure::Raise("not enough points given");
+        }
+
+        GeomAPI_PointsToBSplineSurface surInterpolation;
+        surInterpolation.Interpolate (interpolationPoints);
+        Handle_Geom_BSplineSurface sur(surInterpolation.Surface());
+        this->getGeomBSplineSurfacePtr()->setHandle(sur);
+        Py_Return;
+    }
+    catch (Standard_Failure) {
+        Handle_Standard_Failure e = Standard_Failure::Caught();
+        std::string err = e->GetMessageString();
+        if (err.empty()) err = e->DynamicType()->Name();
+        PyErr_SetString(PyExc_Exception, err.c_str());
+        return 0;
+    }
+}
+
 Py::Int BSplineSurfacePy::getUDegree(void) const
 {
     Handle_Geom_BSplineSurface surf = Handle_Geom_BSplineSurface::DownCast
diff --git a/src/Mod/Part/App/Geometry.cpp b/src/Mod/Part/App/Geometry.cpp
index 7f3b9f71b..2f8450d3a 100644
--- a/src/Mod/Part/App/Geometry.cpp
+++ b/src/Mod/Part/App/Geometry.cpp
@@ -1362,6 +1362,11 @@ GeomBSplineSurface::~GeomBSplineSurface()
 {
 }
 
+void GeomBSplineSurface::setHandle(const Handle_Geom_BSplineSurface& s)
+{
+    mySurface = Handle_Geom_BSplineSurface::DownCast(s->Copy());
+}
+
 const Handle_Geom_Geometry& GeomBSplineSurface::handle() const
 {
     return mySurface;
diff --git a/src/Mod/Part/App/Geometry.h b/src/Mod/Part/App/Geometry.h
index f8a8efa5e..d04797083 100644
--- a/src/Mod/Part/App/Geometry.h
+++ b/src/Mod/Part/App/Geometry.h
@@ -448,6 +448,7 @@ public:
     // Base implementer ----------------------------
     virtual PyObject *getPyObject(void);
 
+    void setHandle(const Handle_Geom_BSplineSurface&);
     const Handle_Geom_Geometry& handle() const;
 
 private: