diff --git a/src/Mod/ReverseEngineering/App/AppReverseEngineering.cpp b/src/Mod/ReverseEngineering/App/AppReverseEngineering.cpp
index a12a6da7f..da79627bc 100644
--- a/src/Mod/ReverseEngineering/App/AppReverseEngineering.cpp
+++ b/src/Mod/ReverseEngineering/App/AppReverseEngineering.cpp
@@ -160,6 +160,10 @@ private:
 
             throw Py::RuntimeError("Computation of B-Spline surface failed");
         }
+        catch (const Py::Exception&) {
+            // re-throw
+            throw;
+        }
         catch (Standard_Failure &e) {
             std::string str;
             Standard_CString msg = e.GetMessageString();
diff --git a/src/Mod/ReverseEngineering/App/ApproxSurface.cpp b/src/Mod/ReverseEngineering/App/ApproxSurface.cpp
index 5c9f4d060..17a74a170 100644
--- a/src/Mod/ReverseEngineering/App/ApproxSurface.cpp
+++ b/src/Mod/ReverseEngineering/App/ApproxSurface.cpp
@@ -26,9 +26,15 @@
 #include <math_Householder.hxx>
 #include <Geom_BSplineSurface.hxx>
 
+#include <QFuture>
+#include <QFutureWatcher>
+#include <QtConcurrentMap>
+#include <boost/bind.hpp>
+
 #include <Mod/Mesh/App/Core/Approximation.h>
 #include <Base/Sequencer.h>
 #include <Base/Tools2D.h>
+#include <Base/Tools.h>
 
 #include "ApproxSurface.h"
 
@@ -174,6 +180,23 @@ void BSplineBasis::AllBasisFunctions(double fParam, TColStd_Array1OfReal& vFuncV
     }
 }
 
+BSplineBasis::ValueT BSplineBasis::LocalSupport(int iIndex, double fParam)
+{
+    int m = _vKnotVector.Length()-1;
+    int p = _iOrder-1;
+
+    if ((iIndex == 0 && fParam == _vKnotVector(0)) || 
+        (iIndex == m-p-1 && fParam == _vKnotVector(m))) {
+        return BSplineBasis::Full;
+    }
+
+    if (fParam < _vKnotVector(iIndex) || fParam >= _vKnotVector(iIndex+p+1)) {
+        return BSplineBasis::Zero;
+    }
+
+    return BSplineBasis::Other;
+}
+
 double BSplineBasis::BasisFunction(int iIndex, double fParam)
 {
     int m = _vKnotVector.Length()-1;
@@ -889,10 +912,29 @@ bool BSplineParameterCorrection::SolveWithoutSmoothing()
         double fV = (*_pvcUVParam)(i).Y();
         unsigned long ulIdx=0;
 
+        // Vorberechnung der Werte der Basis-Funktionen
+        std::vector<double> basisU(_usUCtrlpoints);
         for (unsigned short j=0; j<_usUCtrlpoints; j++) {
-            for (unsigned short k=0; k<_usVCtrlpoints; k++) {
-                M(i,ulIdx) = _clUSpline.BasisFunction(j,fU)*_clVSpline.BasisFunction(k,fV);
-                ulIdx++;
+            basisU[j] = _clUSpline.BasisFunction(j,fU);
+        }
+        std::vector<double> basisV(_usVCtrlpoints);
+        for (unsigned short k=0; k<_usVCtrlpoints; k++) {
+            basisV[k] = _clVSpline.BasisFunction(k,fV);
+        }
+
+        for (unsigned short j=0; j<_usUCtrlpoints; j++) {
+            double valueU = basisU[j];
+            if (valueU == 0.0) {
+                for (unsigned short k=0; k<_usVCtrlpoints; k++) {
+                    M(i,ulIdx) = 0.0;
+                    ulIdx++;
+                }
+            }
+            else {
+                for (unsigned short k=0; k<_usVCtrlpoints; k++) {
+                    M(i,ulIdx) = valueU * basisV[k];
+                    ulIdx++;
+                }
             }
         }
     }
@@ -925,52 +967,116 @@ bool BSplineParameterCorrection::SolveWithoutSmoothing()
     return true;
 }
 
+namespace Reen {
+class ScalarProduct
+{
+public:
+    ScalarProduct(const math_Matrix& mat) : mat(mat)
+    {
+    }
+    std::vector<double> multiply(int col) const
+    {
+        math_Vector vec = mat.Col(col);
+        std::vector<double> out(mat.ColNumber());
+        for (int n=mat.LowerCol(); n<=mat.UpperCol(); n++) {
+            out[n] = vec * mat.Col(n);
+        }
+        return out;
+    }
+
+private:
+    const math_Matrix& mat;
+};
+}
+
 bool BSplineParameterCorrection::SolveWithSmoothing(double fWeight)
 {
     unsigned long ulSize = _pvcPoints->Length();
     unsigned long ulDim  = _usUCtrlpoints*_usVCtrlpoints;
-    math_Matrix M  (0, ulSize-1, 0,_usUCtrlpoints*_usVCtrlpoints-1);
-    math_Matrix MTM(0, _usUCtrlpoints*_usVCtrlpoints-1, 
-                  0, _usUCtrlpoints*_usVCtrlpoints-1);
-    math_Vector Xx (0, _usUCtrlpoints*_usVCtrlpoints-1);
-    math_Vector Xy (0, _usUCtrlpoints*_usVCtrlpoints-1);
-    math_Vector Xz (0, _usUCtrlpoints*_usVCtrlpoints-1);
+    math_Matrix M  (0, ulSize-1, 0, ulDim-1);
+    math_Vector Xx (0, ulDim-1);
+    math_Vector Xy (0, ulDim-1);
+    math_Vector Xz (0, ulDim-1);
     math_Vector bx (0, ulSize-1);
     math_Vector by (0, ulSize-1);
     math_Vector bz (0, ulSize-1);
-    math_Vector Mbx(0, _usUCtrlpoints*_usVCtrlpoints-1);
-    math_Vector Mby(0, _usUCtrlpoints*_usVCtrlpoints-1);
-    math_Vector Mbz(0, _usUCtrlpoints*_usVCtrlpoints-1);
+    math_Vector Mbx(0, ulDim-1);
+    math_Vector Mby(0, ulDim-1);
+    math_Vector Mbz(0, ulDim-1);
 
     //Bestimmung der Koeffizientenmatrix des ueberbestimmten LGS
     for (unsigned long i=0; i<ulSize; i++) {
         double fU = (*_pvcUVParam)(i).X();
         double fV = (*_pvcUVParam)(i).Y();
-        unsigned long ulIdx=0;
+        int ulIdx=0;
+
+        // Vorberechnung der Werte der Basis-Funktionen
+        std::vector<double> basisU(_usUCtrlpoints);
+        for (unsigned short j=0; j<_usUCtrlpoints; j++) {
+            basisU[j] = _clUSpline.BasisFunction(j,fU);
+        }
+        std::vector<double> basisV(_usVCtrlpoints);
+        for (unsigned short k=0; k<_usVCtrlpoints; k++) {
+            basisV[k] = _clVSpline.BasisFunction(k,fV);
+        }
 
         for (unsigned short j=0; j<_usUCtrlpoints; j++) {
-            for (unsigned short k=0; k<_usVCtrlpoints; k++) {
-                M(i,ulIdx) = _clUSpline.BasisFunction(j,fU)*_clVSpline.BasisFunction(k,fV);
-                ulIdx++;
+            double valueU = basisU[j];
+            if (valueU == 0.0) {
+                for (unsigned short k=0; k<_usVCtrlpoints; k++) {
+                    M(i,ulIdx) = 0.0;
+                    ulIdx++;
+                }
+            }
+            else {
+                for (unsigned short k=0; k<_usVCtrlpoints; k++) {
+                    M(i,ulIdx) = valueU * basisV[k];
+                    ulIdx++;
+                }
             }
         }
     }
 
-    //Das Produkt aus ihrer Transformierten und ihr selbst ergibt die quadratische Systemmatrix
+    //Das Produkt aus ihrer Transformierten und ihr selbst ergibt die quadratische Systemmatrix (langsam)
+#if 0
+    math_Matrix MTM = M.TMultiply(M);
+#elif 0
+    math_Matrix MTM(0, ulDim-1, 0, ulDim-1);
     for (unsigned long m=0; m<ulDim; m++) {
+        math_Vector Mm = M.Col(m);
         for (unsigned long n=m; n<ulDim; n++) {
-            MTM(m,n)=MTM(n,m)=M.Col(m)*M.Col(n);
+            MTM(m,n) = MTM(n,m) = Mm * M.Col(n);
         }
     }
+#else // multi-threaded
+    std::vector<int> columns(ulDim);
+    std::generate(columns.begin(), columns.end(), Base::iotaGen<int>(0));
+    ScalarProduct scalar(M);
+    QFuture< std::vector<double> > future = QtConcurrent::mapped
+        (columns, boost::bind(&ScalarProduct::multiply, &scalar, _1));
+    QFutureWatcher< std::vector<double> > watcher;
+    watcher.setFuture(future);
+    watcher.waitForFinished();
+
+    math_Matrix MTM(0, ulDim-1, 0, ulDim-1);
+    int rowIndex=0;
+    for (QFuture< std::vector<double> >::const_iterator it = future.begin(); it != future.end(); ++it) {
+        int colIndex=0;
+        for (std::vector<double>::const_iterator jt = it->begin(); jt != it->end(); ++jt, colIndex++)
+            MTM(rowIndex, colIndex) = *jt;
+        rowIndex++;
+    }
+#endif
 
     //Bestimmen der rechten Seite
     for (int ii=_pvcPoints->Lower(); ii<=_pvcPoints->Upper(); ii++) {
         bx(ii) = (*_pvcPoints)(ii).X(); by(ii) = (*_pvcPoints)(ii).Y(); bz(ii) = (*_pvcPoints)(ii).Z();
     }
     for (unsigned long i=0; i<ulDim; i++) {
-        Mbx(i) = M.Col(i) * bx;
-        Mby(i) = M.Col(i) * by;
-        Mbz(i) = M.Col(i) * bz;
+        math_Vector Mi = M.Col(i);
+        Mbx(i) = Mi * bx;
+        Mby(i) = Mi * by;
+        Mbz(i) = Mi * bz;
     }
 
     // Loese das LGS mit der LU-Zerlegung
@@ -1048,7 +1154,7 @@ void BSplineParameterCorrection::CalcSecondSmoothMatrix(Base::SequencerLauncher&
                 for (unsigned short j=0; j<_usVCtrlpoints; j++) {
                     _clSecondMatrix(m,n) =  _clUSpline.GetIntegralOfProductOfBSplines(i,k,2,2) *
                                             _clVSpline.GetIntegralOfProductOfBSplines(j,l,0,0) +
-                                            2*_clUSpline.GetIntegralOfProductOfBSplines(i,k,1,1) *
+                                          2*_clUSpline.GetIntegralOfProductOfBSplines(i,k,1,1) *
                                             _clVSpline.GetIntegralOfProductOfBSplines(j,l,1,1) +
                                             _clUSpline.GetIntegralOfProductOfBSplines(i,k,0,0) *
                                             _clVSpline.GetIntegralOfProductOfBSplines(j,l,2,2);
diff --git a/src/Mod/ReverseEngineering/App/ApproxSurface.h b/src/Mod/ReverseEngineering/App/ApproxSurface.h
index 0ca340d4b..0fc6f92d2 100644
--- a/src/Mod/ReverseEngineering/App/ApproxSurface.h
+++ b/src/Mod/ReverseEngineering/App/ApproxSurface.h
@@ -45,6 +45,11 @@ namespace Reen {
 class ReenExport SplineBasisfunction
 {
 public:
+    enum ValueT {
+        Zero = 0,
+        Full,
+        Other
+    };
     /**
      * Konstruktor
      * @param iSize Length of Knots vector
@@ -71,6 +76,16 @@ public:
 
     virtual ~SplineBasisfunction();
 
+    /**
+     * Gibt an, ob der Funktionswert Nik(t) an der Stelle fParam
+     * 0, 1 oder ein Wert dazwischen ergibt.
+     * Dies dient dazu, um die Berechnung zu u.U. zu beschleunigen.
+     *
+     * @param iIndex Index
+     * @param fParam Parameterwert
+     * @return ValueT
+     */
+    virtual ValueT LocalSupport(int iIndex, double fParam)=0;
     /**
      * Berechnet den Funktionswert Nik(t) an der Stelle fParam
      * (aus: Piegl/Tiller 96 The NURBS-Book)
@@ -167,6 +182,17 @@ public:
      */
     virtual void AllBasisFunctions(double fParam, TColStd_Array1OfReal& vFuncVals);
 
+    /**
+     * Gibt an, ob der Funktionswert Nik(t) an der Stelle fParam
+     * 0, 1 oder ein Wert dazwischen ergibt.
+     * Dies dient dazu, um die Berechnung zu u.U. zu beschleunigen.
+     *
+     * @param iIndex Index
+     * @param fParam Parameterwert
+     * @return ValueT
+     */
+    virtual ValueT LocalSupport(int iIndex, double fParam);
+
     /**
      * Berechnet den Funktionswert Nik(t) an der Stelle fParam
      * (aus: Piegl/Tiller 96 The NURBS-Book)
diff --git a/src/Mod/ReverseEngineering/App/CMakeLists.txt b/src/Mod/ReverseEngineering/App/CMakeLists.txt
index 60d54a75a..f0d97254b 100644
--- a/src/Mod/ReverseEngineering/App/CMakeLists.txt
+++ b/src/Mod/ReverseEngineering/App/CMakeLists.txt
@@ -18,6 +18,7 @@ include_directories(
     ${EIGEN3_INCLUDE_DIR}
     ${PCL_INCLUDE_DIRS}
     ${FLANN_INCLUDE_DIRS}
+    ${QT_QTCORE_INCLUDE_DIR}
 )
 
 link_directories(${OCC_LIBRARY_DIR})
@@ -32,6 +33,7 @@ set(Reen_LIBS
     ${PCL_FEATURES_LIBRARIES}
     ${PCL_SEARCH_LIBRARIES}
     ${PCL_SURFACE_LIBRARIES}
+    ${QT_QTCORE_LIBRARY}
 )
 
 SET(Reen_SRCS