diff --git a/src/Mod/Mesh/App/CMakeLists.txt b/src/Mod/Mesh/App/CMakeLists.txt
index ee813be7b..97e285680 100644
--- a/src/Mod/Mesh/App/CMakeLists.txt
+++ b/src/Mod/Mesh/App/CMakeLists.txt
@@ -1,20 +1,26 @@
 if(WIN32)
-    add_definitions(-DFCAppMesh -DWM4_FOUNDATION_DLL_EXPORT)
+    add_definitions(-DFCAppMesh -DWM4_FOUNDATION_DLL_EXPORT -DEIGEN2_SUPPORT)
+else (Win32)
+    add_definitions(-DEIGEN2_SUPPORT)
 endif(WIN32)
 
 include_directories(
     ${CMAKE_CURRENT_BINARY_DIR}
     ${CMAKE_CURRENT_SOURCE_DIR}
+    ${CMAKE_SOURCE_DIR}/src/3rdParty
     ${Boost_INCLUDE_DIRS}
     ${PYTHON_INCLUDE_PATH}
     ${XERCESC_INCLUDE_DIR}
+    ${QT_QTCORE_INCLUDE_DIR}
     ${ZLIB_INCLUDE_DIR}
-    ${EIGEN2_INCLUDE_DIR}
+    ${EIGEN3_INCLUDE_DIR}
 )
 
 
 set(Mesh_LIBS
     ${Boost_LIBRARIES}
+    ${QT_QTCORE_LIBRARY}
+    ${QT_QTCORE_LIBRARY_DEBUG}
     FreeCADBase
     FreeCADApp
 )
@@ -41,6 +47,8 @@ SET(Core_SRCS
     Core/Approximation.h
     Core/Builder.cpp
     Core/Builder.h
+    Core/Curvature.cpp
+    Core/Curvature.h
     Core/Definitions.cpp
     Core/Definitions.h
     Core/Degeneration.cpp
@@ -61,6 +69,8 @@ SET(Core_SRCS
     Core/MeshKernel.h
     Core/Projection.cpp
     Core/Projection.h
+    Core/Segmentation.cpp
+    Core/Segmentation.h
     Core/SetOperations.cpp
     Core/SetOperations.h
     Core/Smoothing.cpp
diff --git a/src/Mod/Mesh/App/Core/Algorithm.cpp b/src/Mod/Mesh/App/Core/Algorithm.cpp
index 40d4907d4..281fba7f7 100644
--- a/src/Mod/Mesh/App/Core/Algorithm.cpp
+++ b/src/Mod/Mesh/App/Core/Algorithm.cpp
@@ -1156,7 +1156,6 @@ void MeshAlgorithm::CheckFacets(const Base::ViewProjMethod* pclProj, const Base:
 {
     const MeshPointArray& p = _rclMesh.GetPoints();
     const MeshFacetArray& f = _rclMesh.GetFacets();
-    Base::SequencerLauncher seq("Check facets", f.size());
     Base::Vector3f pt2d;
     unsigned long index=0;
     for (MeshFacetArray::_TConstIterator it = f.begin(); it != f.end(); ++it,++index) {
@@ -1167,7 +1166,6 @@ void MeshAlgorithm::CheckFacets(const Base::ViewProjMethod* pclProj, const Base:
                 break;
             }
         }
-        seq.next();
     }
 }
 
@@ -1743,47 +1741,38 @@ std::set<unsigned long> MeshRefPointToFacets::NeighbourPoints(const std::vector<
     return nb;
 }
 
-// ermittelt alle Nachbarn zum Facet deren Schwerpunkt unterhalb der mpx. Distanz befindet. 
-// Facet deren VISIT-Flag gesetzt ist werden nicht beruecksichtig. 
-/// @todo
-void MeshRefPointToFacets::Neighbours (unsigned long ulFacetInd, float fMpxDist, std::vector<MeshFacetArray::_TConstIterator> &rclNb)
+void MeshRefPointToFacets::Neighbours (unsigned long ulFacetInd, float fMaxDist, MeshCollector& collect) const
 {
-    rclNb.clear();
+    std::set<unsigned long> visited;
     Base::Vector3f  clCenter = _rclMesh.GetFacet(ulFacetInd).GetGravityPoint();
 
     const MeshFacetArray& rFacets = _rclMesh.GetFacets();
-    SearchNeighbours(rFacets.begin() + ulFacetInd, clCenter, fMpxDist * fMpxDist, rclNb);
-
-    for (std::vector<MeshFacetArray::_TConstIterator>::iterator i = rclNb.begin(); i != rclNb.end(); i++)
-        (*i)->ResetFlag(MeshFacet::VISIT);  
+    SearchNeighbours(rFacets, ulFacetInd, clCenter, fMaxDist * fMaxDist, visited, collect);
 }
 
-/// @todo
-void MeshRefPointToFacets::SearchNeighbours(MeshFacetArray::_TConstIterator f_it, const Base::Vector3f &rclCenter, float fMpxDist, std::vector<MeshFacetArray::_TConstIterator> &rclNb)
+void MeshRefPointToFacets::SearchNeighbours(const MeshFacetArray& rFacets, unsigned long index, const Base::Vector3f &rclCenter,
+                                            float fMaxDist2, std::set<unsigned long>& visited, MeshCollector& collect) const
 {
-    if (f_it->IsFlag(MeshFacet::VISIT) == true)
+    if (visited.find(index) != visited.end())
         return;
 
-    if (Base::DistanceP2(rclCenter, _rclMesh.GetFacet(*f_it).GetGravityPoint()) > fMpxDist)
+    const MeshFacet& face = rFacets[index];
+    if (Base::DistanceP2(rclCenter, _rclMesh.GetFacet(face).GetGravityPoint()) > fMaxDist2)
         return;
 
-    rclNb.push_back(f_it);
-    f_it->SetFlag(MeshFacet::VISIT);
-
-    MeshPointArray::_TConstIterator p_beg = _rclMesh.GetPoints().begin();
-    MeshFacetArray::_TConstIterator f_beg = _rclMesh.GetFacets().begin();
-
+    visited.insert(index);
+    collect.Append(_rclMesh, index);
     for (int i = 0; i < 3; i++) {
-        const std::set<unsigned long> &f = (*this)[f_it->_aulPoints[i]];
+        const std::set<unsigned long> &f = (*this)[face._aulPoints[i]];
 
         for (std::set<unsigned long>::const_iterator j = f.begin(); j != f.end(); ++j) {
-            SearchNeighbours(f_beg+*j, rclCenter, fMpxDist, rclNb);
+            SearchNeighbours(rFacets, *j, rclCenter, fMaxDist2, visited, collect);
         }
     }
 }
 
 MeshFacetArray::_TConstIterator
-MeshRefPointToFacets::getFacet (unsigned long index) const
+MeshRefPointToFacets::GetFacet (unsigned long index) const
 {
     return _rclMesh.GetFacets().begin() + index;
 }
diff --git a/src/Mod/Mesh/App/Core/Algorithm.h b/src/Mod/Mesh/App/Core/Algorithm.h
index 9baa667e2..78a4caf5b 100644
--- a/src/Mod/Mesh/App/Core/Algorithm.h
+++ b/src/Mod/Mesh/App/Core/Algorithm.h
@@ -310,6 +310,43 @@ protected:
   const MeshKernel      &_rclMesh; /**< The mesh kernel. */
 };
 
+class MeshExport MeshCollector
+{
+public:
+    MeshCollector(){}
+    virtual void Append(const MeshCore::MeshKernel&, unsigned long index) = 0;
+};
+
+class MeshExport PointCollector : public MeshCollector
+{
+public:
+    PointCollector(std::vector<unsigned long>& ind) : indices(ind){}
+    virtual void Append(const MeshCore::MeshKernel& kernel, unsigned long index)
+    {
+        unsigned long ulP1, ulP2, ulP3;
+        kernel.GetFacetPoints(index, ulP1, ulP2, ulP3);
+        indices.push_back(ulP1);
+        indices.push_back(ulP2);
+        indices.push_back(ulP3);
+    }
+
+private:
+    std::vector<unsigned long>& indices;
+};
+
+class MeshExport FacetCollector : public MeshCollector
+{
+public:
+    FacetCollector(std::vector<unsigned long>& ind) : indices(ind){}
+    void Append(const MeshCore::MeshKernel&, unsigned long index)
+    {
+        indices.push_back(index);
+    }
+
+private:
+    std::vector<unsigned long>& indices;
+};
+
 /**
  * The MeshRefPointToFacets builds up a structure to have access to all facets indexing
  * a point.
@@ -329,14 +366,14 @@ public:
     /// Rebuilds up data structure
     void Rebuild (void);
     const std::set<unsigned long>& operator[] (unsigned long) const;
-    MeshFacetArray::_TConstIterator getFacet (unsigned long) const;
+    MeshFacetArray::_TConstIterator GetFacet (unsigned long) const;
     std::set<unsigned long> NeighbourPoints(const std::vector<unsigned long>& , int level) const;
-    void Neighbours (unsigned long ulFacetInd, float fMaxDist, std::vector<MeshFacetArray::_TConstIterator> &rclNb);
+    void Neighbours (unsigned long ulFacetInd, float fMaxDist, MeshCollector& collect) const;
     Base::Vector3f GetNormal(unsigned long) const;
 
 protected:
-    void SearchNeighbours(MeshFacetArray::_TConstIterator pFIter, const Base::Vector3f &rclCenter, 
-        float fMaxDist, std::vector<MeshFacetArray::_TConstIterator> &rclNb);
+    void SearchNeighbours(const MeshFacetArray& rFacets, unsigned long index, const Base::Vector3f &rclCenter, 
+        float fMaxDist, std::set<unsigned long> &visit, MeshCollector& collect) const;
 
 protected:
     const MeshKernel  &_rclMesh; /**< The mesh kernel. */
diff --git a/src/Mod/Mesh/App/Core/Approximation.cpp b/src/Mod/Mesh/App/Core/Approximation.cpp
index ace1a796e..15b0c522f 100644
--- a/src/Mod/Mesh/App/Core/Approximation.cpp
+++ b/src/Mod/Mesh/App/Core/Approximation.cpp
@@ -37,6 +37,7 @@
 #include <Mod/Mesh/App/WildMagic4/Wm4ApprPolyFit3.h>
 
 //#define FC_USE_EIGEN
+//#define FC_USE_BOOST
 #if defined(FC_USE_BOOST)
 #include <boost/numeric/ublas/io.hpp>
 #include <boost/numeric/bindings/traits/ublas_matrix.hpp>
@@ -54,25 +55,26 @@ extern "C" void LAPACK_DGESV (int const* n, int const* nrhs,
 #elif defined(FC_USE_EIGEN)
 # include <Eigen/LeastSquares>
 #endif
+# include <Eigen/LeastSquares>
 
 using namespace MeshCore;
 
 
-void Approximation::Convert( const Wm4::Vector3<float>& Wm4, Base::Vector3f& pt)
+void Approximation::Convert( const Wm4::Vector3<double>& Wm4, Base::Vector3f& pt)
 {
-    pt.Set( Wm4.X(), Wm4.Y(), Wm4.Z() );
+    pt.Set( (float)Wm4.X(), (float)Wm4.Y(), (float)Wm4.Z() );
 }
 
-void Approximation::Convert( const Base::Vector3f& pt, Wm4::Vector3<float>& Wm4)
+void Approximation::Convert( const Base::Vector3f& pt, Wm4::Vector3<double>& Wm4)
 {
     Wm4.X() = pt.x; Wm4.Y() = pt.y; Wm4.Z() = pt.z;
 }
 
-void Approximation::GetMgcVectorArray(std::vector< Wm4::Vector3<float> >& rcPts) const
+void Approximation::GetMgcVectorArray(std::vector< Wm4::Vector3<double> >& rcPts) const
 {
     std::list< Base::Vector3f >::const_iterator It;
     for (It = _vPoints.begin(); It != _vPoints.end(); ++It) {
-        Wm4::Vector3<float> pt( (*It).x, (*It).y, (*It).z );
+        Wm4::Vector3<double> pt( (*It).x, (*It).y, (*It).z );
         rcPts.push_back( pt );
     }
 }
@@ -176,7 +178,7 @@ float PlaneFit::Fit()
     int r = lapack::syev('V','U',A,eigenval,lapack::optimal_workspace());
     if (r) {
     }
-    float sigma;
+    float sigma = 0;
 #elif defined(FC_USE_EIGEN)
     Eigen::Matrix3d covMat = Eigen::Matrix3d::Zero();
     covMat(0,0) = sxx;
@@ -360,17 +362,15 @@ void PlaneFit::ProjectToPlane ()
 
 // -------------------------------------------------------------------------------
 
-float FunctionContainer::dKoeff[]; // Koeffizienten der Quadrik
-
-bool QuadraticFit::GetCurvatureInfo(float x, float y, float z,
-                                    float &rfCurv0, float &rfCurv1,
-                                    Base::Vector3f &rkDir0, Base::Vector3f &rkDir1, float &dDistance)
+bool QuadraticFit::GetCurvatureInfo(double x, double y, double z,
+                                    double &rfCurv0, double &rfCurv1,
+                                    Base::Vector3f &rkDir0, Base::Vector3f &rkDir1, double &dDistance)
 {
     assert( _bIsFitted );
     bool bResult = false;
 
     if (_bIsFitted) {
-        Wm4::Vector3<float> Dir0, Dir1;
+        Wm4::Vector3<double> Dir0, Dir1;
         FunctionContainer  clFuncCont( _fCoeff );
         bResult = clFuncCont.CurvatureInfo( x, y, z, rfCurv0, rfCurv1, Dir0, Dir1, dDistance );
 
@@ -382,7 +382,7 @@ bool QuadraticFit::GetCurvatureInfo(float x, float y, float z,
     return bResult;
 }
 
-bool QuadraticFit::GetCurvatureInfo(float x, float y, float z, float &rfCurv0, float &rfCurv1)
+bool QuadraticFit::GetCurvatureInfo(double x, double y, double z, double &rfCurv0, double &rfCurv1)
 {
     bool bResult = false;
 
@@ -394,12 +394,12 @@ bool QuadraticFit::GetCurvatureInfo(float x, float y, float z, float &rfCurv0, f
     return bResult;
 }
 
-const float& QuadraticFit::GetCoeffArray() const
+const double& QuadraticFit::GetCoeffArray() const
 {
     return _fCoeff[0];
 }
 
-float QuadraticFit::GetCoeff(unsigned long ulIndex) const
+double QuadraticFit::GetCoeff(unsigned long ulIndex) const
 {
     assert( ulIndex >= 0 && ulIndex < 10 );
 
@@ -414,9 +414,9 @@ float QuadraticFit::Fit()
     float fResult = FLOAT_MAX;
 
     if (CountPoints() > 0) {
-        std::vector< Wm4::Vector3<float> > cPts;
+        std::vector< Wm4::Vector3<double> > cPts;
         GetMgcVectorArray( cPts );
-        fResult = Wm4::QuadraticFit3<float>( CountPoints(), &(cPts[0]), _fCoeff );
+        fResult = Wm4::QuadraticFit3<double>( CountPoints(), &(cPts[0]), _fCoeff );
         _fLastResult = fResult;
 
         _bIsFitted = true;
@@ -425,7 +425,7 @@ float QuadraticFit::Fit()
     return fResult;
 }
 
-void QuadraticFit::CalcEigenValues(float &dLambda1, float &dLambda2, float &dLambda3,
+void QuadraticFit::CalcEigenValues(double &dLambda1, double &dLambda2, double &dLambda3,
                                    Base::Vector3f &clEV1, Base::Vector3f &clEV2, Base::Vector3f &clEV3) const
 {
     assert( _bIsFitted );
@@ -451,16 +451,16 @@ void QuadraticFit::CalcEigenValues(float &dLambda1, float &dLambda2, float &dLam
      *
      */
 
-    Wm4::Matrix3<float>  akMat(_fCoeff[4],       _fCoeff[7]/2.0f, _fCoeff[8]/2.0f,
+    Wm4::Matrix3<double>  akMat(_fCoeff[4],       _fCoeff[7]/2.0f, _fCoeff[8]/2.0f,
                                _fCoeff[7]/2.0f,  _fCoeff[5],      _fCoeff[9]/2.0f,
                                _fCoeff[8]/2.0f,  _fCoeff[9]/2.0f, _fCoeff[6]       );
 
-    Wm4::Matrix3<float> rkRot, rkDiag;
+    Wm4::Matrix3<double> rkRot, rkDiag;
     akMat.EigenDecomposition( rkRot, rkDiag );
 
-    Wm4::Vector3<float> vEigenU = rkRot.GetColumn(0);
-    Wm4::Vector3<float> vEigenV = rkRot.GetColumn(1);
-    Wm4::Vector3<float> vEigenW = rkRot.GetColumn(2);
+    Wm4::Vector3<double> vEigenU = rkRot.GetColumn(0);
+    Wm4::Vector3<double> vEigenV = rkRot.GetColumn(1);
+    Wm4::Vector3<double> vEigenW = rkRot.GetColumn(2);
 
     Convert( vEigenU, clEV1 );
     Convert( vEigenV, clEV2 );
@@ -471,11 +471,11 @@ void QuadraticFit::CalcEigenValues(float &dLambda1, float &dLambda2, float &dLam
     dLambda3 = rkDiag[2][2];
 }
 
-void QuadraticFit::CalcZValues( float x, float y, float &dZ1, float &dZ2 ) const
+void QuadraticFit::CalcZValues( double x, double y, double &dZ1, double &dZ2 ) const
 {
     assert( _bIsFitted );
 
-    float dDisk = _fCoeff[3]*_fCoeff[3]+2*_fCoeff[3]*_fCoeff[8]*x+2*_fCoeff[3]*_fCoeff[9]*y+
+    double dDisk = _fCoeff[3]*_fCoeff[3]+2*_fCoeff[3]*_fCoeff[8]*x+2*_fCoeff[3]*_fCoeff[9]*y+
                   _fCoeff[8]*_fCoeff[8]*x*x+2*_fCoeff[8]*x*_fCoeff[9]*y+_fCoeff[9]*_fCoeff[9]*y*y-
                   4*_fCoeff[6]*_fCoeff[0]-4*_fCoeff[6]*_fCoeff[1]*x-4*_fCoeff[6]*_fCoeff[2]*y-
                   4*_fCoeff[6]*_fCoeff[7]*x*y-4*_fCoeff[6]*_fCoeff[4]*x*x-4*_fCoeff[6]*_fCoeff[5]*y*y;
@@ -494,8 +494,8 @@ void QuadraticFit::CalcZValues( float x, float y, float &dZ1, float &dZ2 ) const
     else
         dDisk = sqrt( dDisk );
 
-    dZ1 = 0.5f * ( ( -_fCoeff[3] - _fCoeff[8]*x - _fCoeff[9]*y + dDisk ) / _fCoeff[6] );
-    dZ2 = 0.5f * ( ( -_fCoeff[3] - _fCoeff[8]*x - _fCoeff[9]*y - dDisk ) / _fCoeff[6] );
+    dZ1 = 0.5 * ( ( -_fCoeff[3] - _fCoeff[8]*x - _fCoeff[9]*y + dDisk ) / _fCoeff[6] );
+    dZ2 = 0.5 * ( ( -_fCoeff[3] - _fCoeff[8]*x - _fCoeff[9]*y - dDisk ) / _fCoeff[6] );
 }
 
 // -------------------------------------------------------------------------------
@@ -529,13 +529,13 @@ float SurfaceFit::Fit()
     return fResult;
 }
 
-bool SurfaceFit::GetCurvatureInfo(float x, float y, float z, float &rfCurv0, float &rfCurv1,
-                                  Base::Vector3f &rkDir0, Base::Vector3f &rkDir1, float &dDistance )
+bool SurfaceFit::GetCurvatureInfo(double x, double y, double z, double &rfCurv0, double &rfCurv1,
+                                  Base::Vector3f &rkDir0, Base::Vector3f &rkDir1, double &dDistance )
 {
     bool bResult = false;
 
     if (_bIsFitted) {
-        Wm4::Vector3<float> Dir0, Dir1;
+        Wm4::Vector3<double> Dir0, Dir1;
         FunctionContainer  clFuncCont( _fCoeff );
         bResult = clFuncCont.CurvatureInfo( x, y, z, rfCurv0, rfCurv1, Dir0, Dir1, dDistance );
 
@@ -547,7 +547,7 @@ bool SurfaceFit::GetCurvatureInfo(float x, float y, float z, float &rfCurv0, flo
     return bResult;
 }
 
-bool SurfaceFit::GetCurvatureInfo(float x, float y, float z, float &rfCurv0, float &rfCurv1)
+bool SurfaceFit::GetCurvatureInfo(double x, double y, double z, double &rfCurv0, double &rfCurv1)
 {
     assert( _bIsFitted );
     bool bResult = false;
@@ -560,18 +560,17 @@ bool SurfaceFit::GetCurvatureInfo(float x, float y, float z, float &rfCurv0, flo
     return bResult;
 }
 
-float SurfaceFit::PolynomFit()
+double SurfaceFit::PolynomFit()
 {
     if (PlaneFit::Fit() == FLOAT_MAX)
         return FLOAT_MAX;
 
-#if 0
-#if defined(FC_USE_BOOST)
     Base::Vector3d bs(this->_vBase.x,this->_vBase.y,this->_vBase.z);
     Base::Vector3d ex(this->_vDirU.x,this->_vDirU.y,this->_vDirU.z);
     Base::Vector3d ey(this->_vDirV.x,this->_vDirV.y,this->_vDirV.z);
     Base::Vector3d ez(this->_vDirW.x,this->_vDirW.y,this->_vDirW.z);
 
+#if defined(FC_USE_BOOST)
     ublas::matrix<double> A(6,6);
     ublas::vector<double> b(6);
     for (int i=0; i<6; i++) {
@@ -580,6 +579,11 @@ float SurfaceFit::PolynomFit()
         }
         b(i) = 0.0;
     }
+#else
+    Eigen::Matrix<double,6,6> A = Eigen::Matrix<double,6,6>::Zero();
+    Eigen::Matrix<double,6,1> b = Eigen::Matrix<double,6,1>::Zero();
+    Eigen::Matrix<double,6,1> x = Eigen::Matrix<double,6,1>::Zero();
+#endif
 
     for (std::list<Base::Vector3f>::const_iterator it = _vPoints.begin(); it != _vPoints.end(); it++) {
         Base::Vector3d clPoint(it->x,it->y,it->z);
@@ -648,29 +652,34 @@ float SurfaceFit::PolynomFit()
 
     A(5,4) = A(4,5);
 
-
+#if defined(FC_USE_BOOST)
     namespace lapack= boost::numeric::bindings::lapack;
     //std::clog << A << std::endl;
     //std::clog << b << std::endl;
-    lapack::gesv(A,b);
+    //lapack::gesv(A,b);
+    ublas::vector<double> x(6);
+    x = b;
     //std::clog << b << std::endl;
+#else
+    // A.llt().solve(b,&x); // not sure if always positive definite
+    A.qr().solve(b,&x);
+#endif
 
-    _fCoeff[0] = (float)(-b(5));
-    _fCoeff[1] = (float)(-b(3));
-    _fCoeff[2] = (float)(-b(4));
+    _fCoeff[0] = (float)(-x(5));
+    _fCoeff[1] = (float)(-x(3));
+    _fCoeff[2] = (float)(-x(4));
     _fCoeff[3] = 1.0f;
-    _fCoeff[4] = (float)(-b(0));
-    _fCoeff[5] = (float)(-b(1));
+    _fCoeff[4] = (float)(-x(0));
+    _fCoeff[5] = (float)(-x(1));
     _fCoeff[6] = 0.0f;
-    _fCoeff[7] = (float)(-b(2));
+    _fCoeff[7] = (float)(-x(2));
     _fCoeff[8] = 0.0f;
     _fCoeff[9] = 0.0f;
-#endif
-#endif
+
     return 0.0f;
 }
 
-float SurfaceFit::Value(float x, float y) const
+double SurfaceFit::Value(double x, double y) const
 {
     float z = 0.0f;
     if (_bIsFitted) {
diff --git a/src/Mod/Mesh/App/Core/Approximation.h b/src/Mod/Mesh/App/Core/Approximation.h
index 74124f35e..623a17dac 100644
--- a/src/Mod/Mesh/App/Core/Approximation.h
+++ b/src/Mod/Mesh/App/Core/Approximation.h
@@ -27,6 +27,7 @@
 #include <Mod/Mesh/App/WildMagic4/Wm4Vector3.h>
 #include <Mod/Mesh/App/WildMagic4/Wm4QuadricSurface.h>
 #include <Mod/Mesh/App/WildMagic4/Wm4Eigen.h>
+#include <Mod/Mesh/App/WildMagic4/Wm4ImplicitSurface.h>
 #include <list>
 #include <set>
 #include <vector>
@@ -34,8 +35,64 @@
 #include <Base/Vector3D.h>
 #include <Base/Matrix.h>
 
-namespace MeshCore {
+namespace Wm4
+{
 
+/**
+ * An implicit surface is defined by F(x,y,z) = 0.
+ * This polynomial surface is actually defined as z = f(x,y) = ax^2 + by^2 + cx + dy + exy + g.
+ * To use Wm3 routines for implicit surfaces we can write the surface also as F(x,y,z) = f(x,y) - z = 0.
+ * @author Werner Mayer
+ */
+template <class Real>
+class PolynomialSurface : public ImplicitSurface<Real>
+{
+public:
+  PolynomialSurface (const Real afCoeff[6])
+  { for (int i=0; i<6; i++) m_afCoeff[i] = afCoeff[i]; }
+
+  virtual ~PolynomialSurface () {}
+
+  // the function
+  virtual Real F (const Vector3<Real>& rkP) const
+  { 
+    return ( m_afCoeff[0]*rkP.X()*rkP.X() + 
+             m_afCoeff[1]*rkP.Y()*rkP.Y() + 
+             m_afCoeff[2]*rkP.X()         + 
+             m_afCoeff[3]*rkP.Y()         + 
+             m_afCoeff[4]*rkP.X()*rkP.Y() + 
+             m_afCoeff[5]-rkP.Z())        ;
+  }
+
+  // first-order partial derivatives
+  virtual Real FX (const Vector3<Real>& rkP) const
+  { return (Real)(2.0*m_afCoeff[0]*rkP.X() + m_afCoeff[2] + m_afCoeff[4]*rkP.Y()); }
+  virtual Real FY (const Vector3<Real>& rkP) const
+  { return (Real)(2.0*m_afCoeff[1]*rkP.Y() + m_afCoeff[3] + m_afCoeff[4]*rkP.X()); }
+  virtual Real FZ (const Vector3<Real>& rkP) const
+  { return (Real)-1.0; }
+
+  // second-order partial derivatives
+  virtual Real FXX (const Vector3<Real>& rkP) const
+  { return (Real)(2.0*m_afCoeff[0]); }
+  virtual Real FXY (const Vector3<Real>& rkP) const
+  { return (Real)(m_afCoeff[4]); }
+  virtual Real FXZ (const Vector3<Real>& rkP) const
+  { return (Real)0.0; }
+  virtual Real FYY (const Vector3<Real>& rkP) const
+  { return (Real)(2.0*m_afCoeff[1]); }
+  virtual Real FYZ (const Vector3<Real>& rkP) const
+  { return (Real)0.0; }
+  virtual Real FZZ (const Vector3<Real>& rkP) const
+  { return (Real)0.0; }
+
+protected:
+  Real m_afCoeff[6];
+};
+
+}
+
+namespace MeshCore {
 
 /**
  * Abstract base class for approximation of a geometry to a given set of points.
@@ -104,15 +161,15 @@ protected:
     /**
      * Converts point from Wm4::Vector3 to Base::Vector3f.
      */
-    static void Convert( const Wm4::Vector3<float>&, Base::Vector3f&);
+    static void Convert( const Wm4::Vector3<double>&, Base::Vector3f&);
     /**
      * Converts point from Base::Vector3f to Wm4::Vector3.
      */
-    static void Convert( const Base::Vector3f&, Wm4::Vector3<float>&);
+    static void Convert( const Base::Vector3f&, Wm4::Vector3<double>&);
     /**
      * Creates a vector of Wm4::Vector3 elements.
      */
-    void GetMgcVectorArray( std::vector< Wm4::Vector3<float> >& rcPts ) const;
+    void GetMgcVectorArray( std::vector< Wm4::Vector3<double> >& rcPts ) const;
 
 protected:
     std::list< Base::Vector3f > _vPoints; /**< Holds the points for the fit algorithm.  */
@@ -201,22 +258,22 @@ public:
     /**
      * Übertragen der Quadric-Koeffizienten
      * @param ulIndex Nummer des Koeffizienten (0..9)
-     * @return float Wert des Koeffizienten
+     * @return double Wert des Koeffizienten
      */
-    float GetCoeff(unsigned long ulIndex) const;
+    double GetCoeff(unsigned long ulIndex) const;
     /**
      * Übertragen der Koeffizientan als Referenz
      * auf das interne Array
-     * @return const float& Referenz auf das float-Array
+     * @return const double& Referenz auf das double-Array
      */
-    const float& GetCoeffArray() const;
+    const double& GetCoeffArray() const;
     /**
      * Aufruf des Fit-Algorithmus
      * @return float Qualität des Fits.
      */
     float Fit();
 
-    void CalcZValues(float x, float y, float &dZ1, float &dZ2) const;
+    void CalcZValues(double x, double y, double &dZ1, double &dZ2) const;
     /**
      * Berechnen der Krümmungswerte der Quadric in einem bestimmten Punkt.
      * @param x X-Koordinate
@@ -229,12 +286,12 @@ public:
      * @param dDistance
      * @return bool Fehlerfreie Ausfürhung = true, ansonsten false
      */
-    bool GetCurvatureInfo(float x, float y, float z,
-                          float &rfCurv0, float &rfCurv1,
-                          Base::Vector3f &rkDir0, Base::Vector3f &rkDir1, float &dDistance);
+    bool GetCurvatureInfo(double x, double y, double z,
+                          double &rfCurv0, double &rfCurv1,
+                          Base::Vector3f &rkDir0, Base::Vector3f &rkDir1, double &dDistance);
 
-    bool GetCurvatureInfo(float x, float y, float z,
-                          float &rfCurv0, float &rfcurv1);
+    bool GetCurvatureInfo(double x, double y, double z,
+                          double &rfCurv0, double &rfcurv1);
     /**
      * Aufstellen der Formanmatrix A und Berechnen der Eigenwerte.
      * @param dLambda1 Eigenwert 1
@@ -244,11 +301,11 @@ public:
      * @param clEV2    Eigenvektor 2
      * @param clEV3    Eigenvektor 3
      */
-    void CalcEigenValues(float &dLambda1, float &dLambda2, float &dLambda3,
+    void CalcEigenValues(double &dLambda1, double &dLambda2, double &dLambda3,
                          Base::Vector3f &clEV1, Base::Vector3f &clEV2, Base::Vector3f &clEV3) const;
 
 protected:
-    float _fCoeff[ 10 ];  /**< Ziel der Koeffizienten aus dem Fit */
+    double _fCoeff[ 10 ];  /**< Ziel der Koeffizienten aus dem Fit */
 };
 
 // -------------------------------------------------------------------------------
@@ -275,15 +332,15 @@ public:
      */
     virtual ~SurfaceFit(){};
 
-    bool GetCurvatureInfo(float x, float y, float z, float &rfCurv0, float &rfCurv1,
-                          Base::Vector3f &rkDir0, Base::Vector3f &rkDir1, float &dDistance);
-    bool GetCurvatureInfo(float x, float y, float z, float &rfCurv0, float &rfcurv1);
+    bool GetCurvatureInfo(double x, double y, double z, double &rfCurv0, double &rfCurv1,
+                          Base::Vector3f &rkDir0, Base::Vector3f &rkDir1, double &dDistance);
+    bool GetCurvatureInfo(double x, double y, double z, double &rfCurv0, double &rfcurv1);
     float Fit();
-    float Value(float x, float y) const;
+    double Value(double x, double y) const;
 
 protected:
-    float PolynomFit();
-    float _fCoeff[ 10 ];  /**< Ziel der Koeffizienten aus dem Fit */
+    double PolynomFit();
+    double _fCoeff[ 10 ];  /**< Ziel der Koeffizienten aus dem Fit */
 };
 
 // -------------------------------------------------------------------------------
@@ -300,42 +357,24 @@ public:
      * Die MGC-Algorithmen arbeiten mit Funktionen dieses
      * Types
      */
-    typedef float (*Function)(float,float,float);
+    typedef double (*Function)(double,double,double);
     /**
      * Der parametrisierte Konstruktor. Erwartet ein Array
      * mit den Quadric-Koeffizienten.
      * @param pKoef Zeiger auf die Quadric-Parameter
-     *        (float [10])
+     *        (double [10])
      */
-    FunctionContainer(const float *pKoef)
+    FunctionContainer(const double *pKoef)
     {
         Assign( pKoef );
-/*
-        Function oF;     
-        Function aoDF[3];
-        Function aoD2F[6];
-
-        oF = &F;
-        aoDF[0] = &Fx;
-        aoDF[1] = &Fy;
-        aoDF[2] = &Fz;
-        aoD2F[0] = &Fxx;
-        aoD2F[1] = &Fxy;
-        aoD2F[2] = &Fxz;
-        aoD2F[3] = &Fyy;
-        aoD2F[4] = &Fyz;
-        aoD2F[5] = &Fzz;
-
-        pImplSurf = new Wm4::QuadricSurface<float>( oF, aoDF, aoD2F );*/
-
-        pImplSurf = new Wm4::QuadricSurface<float>( dKoeff );
+        pImplSurf = new Wm4::QuadricSurface<double>( dKoeff );
     }
     /**
      * Übernehmen der Quadric-Parameter
      * @param pKoef Zeiger auf die Quadric-Parameter
-     *        (doube [10])
+     *        (double [10])
      */
-    void Assign( const float *pKoef )
+    void Assign( const double *pKoef )
     {
         for (long ct=0; ct < 10; ct++)
             dKoeff[ ct ] = pKoef[ ct ];
@@ -344,13 +383,13 @@ public:
      * Destruktor. Löscht die ImpicitSurface Klasse
      * der MGC-Bibliothek wieder
      */
-    virtual ~FunctionContainer(){ delete pImplSurf; }
+    ~FunctionContainer(){ delete pImplSurf; }
     /** 
      * Zugriff auf die Koeffizienten der Quadric
      * @param idx Index des Parameters
-     * @return float& Der Koeffizient
+     * @return double& Der Koeffizient
      */
-    float& operator[](int idx){ return dKoeff[ idx ]; }
+    double& operator[](int idx){ return dKoeff[ idx ]; }
     /**
      * Redirector auf eine Methode der MGC Bibliothek. Ermittelt
      * die Hauptkrümmungen und ihre Richtungen im angegebenen Punkt.
@@ -364,22 +403,22 @@ public:
      * @param dDistance Ergebnis das die Entfernung des Punktes von der Quadrik angibt.
      * @return bool Fehlerfreie Ausfürhung = true, ansonsten false
      */
-    bool CurvatureInfo(float x, float y, float z, 
-                       float &rfCurv0, float &rfCurv1,
-                       Wm4::Vector3<float> &rkDir0,  Wm4::Vector3<float> &rkDir1, float &dDistance)
+    bool CurvatureInfo(double x, double y, double z, 
+                       double &rfCurv0, double &rfCurv1,
+                       Wm4::Vector3<double> &rkDir0,  Wm4::Vector3<double> &rkDir1, double &dDistance)
     {
-        return pImplSurf->ComputePrincipalCurvatureInfo( Wm4::Vector3<float>(x, y, z),rfCurv0, rfCurv1, rkDir0, rkDir1 );
+        return pImplSurf->ComputePrincipalCurvatureInfo( Wm4::Vector3<double>(x, y, z),rfCurv0, rfCurv1, rkDir0, rkDir1 );
     }
 
-    Base::Vector3f GetGradient( float x, float y, float z ) const
+    Base::Vector3f GetGradient( double x, double y, double z ) const
     {
-        Wm4::Vector3<float> grad = pImplSurf->GetGradient( Wm4::Vector3<float>(x, y, z) );
+        Wm4::Vector3<double> grad = pImplSurf->GetGradient( Wm4::Vector3<double>(x, y, z) );
         return Base::Vector3f( grad.X(), grad.Y(), grad.Z() );
     }
 
-    Base::Matrix4D GetHessian( float x, float y, float z ) const
+    Base::Matrix4D GetHessian( double x, double y, double z ) const
     {
-        Wm4::Matrix3<float> hess = pImplSurf->GetHessian( Wm4::Vector3<float>(x, y, z) );
+        Wm4::Matrix3<double> hess = pImplSurf->GetHessian( Wm4::Vector3<double>(x, y, z) );
         Base::Matrix4D cMat; cMat.setToUnity();
         cMat[0][0] = hess[0][0]; cMat[0][1] = hess[0][1]; cMat[0][2] = hess[0][2];
         cMat[1][0] = hess[1][0]; cMat[1][1] = hess[1][1]; cMat[1][2] = hess[1][2];
@@ -387,23 +426,23 @@ public:
         return cMat;
     }
 
-    bool CurvatureInfo(float x, float y, float z,
-                       float &rfCurv0, float &rfCurv1)
+    bool CurvatureInfo(double x, double y, double z,
+                       double &rfCurv0, double &rfCurv1)
     {
-        float dQuot = Fz(x,y,z);
-        float zx = - ( Fx(x,y,z) / dQuot );
-        float zy = - ( Fy(x,y,z) / dQuot );
+        double dQuot = Fz(x,y,z);
+        double zx = - ( Fx(x,y,z) / dQuot );
+        double zy = - ( Fy(x,y,z) / dQuot );
         
-        float zxx = - ( 2.0f * ( dKoeff[5] + dKoeff[6] * zx * zx + dKoeff[8] * zx ) ) / dQuot;
-        float zyy = - ( 2.0f * ( dKoeff[5] + dKoeff[6] * zy * zy + dKoeff[9] * zy ) ) / dQuot;
-        float zxy = - ( dKoeff[6] * zx * zy + dKoeff[7] + dKoeff[8] * zy + dKoeff[9] * zx ) / dQuot;
+        double zxx = - ( 2.0f * ( dKoeff[5] + dKoeff[6] * zx * zx + dKoeff[8] * zx ) ) / dQuot;
+        double zyy = - ( 2.0f * ( dKoeff[5] + dKoeff[6] * zy * zy + dKoeff[9] * zy ) ) / dQuot;
+        double zxy = - ( dKoeff[6] * zx * zy + dKoeff[7] + dKoeff[8] * zy + dKoeff[9] * zx ) / dQuot;
 
-        float dNen = 1 + zx*zx + zy*zy;
-        float dNenSqrt = (float)sqrt( dNen );
-        float K = ( zxx * zyy - zxy * zxy ) / ( dNen * dNen );
-        float H = 0.5f * ( ( 1.0f+zx*zx - 2*zx*zy*zxy + (1.0f+zy*zy)*zxx ) / ( dNenSqrt * dNenSqrt * dNenSqrt ) ) ;
+        double dNen = 1 + zx*zx + zy*zy;
+        double dNenSqrt = (double)sqrt( dNen );
+        double K = ( zxx * zyy - zxy * zxy ) / ( dNen * dNen );
+        double H = 0.5f * ( ( 1.0f+zx*zx - 2*zx*zy*zxy + (1.0f+zy*zy)*zxx ) / ( dNenSqrt * dNenSqrt * dNenSqrt ) ) ;
 
-        float dDiscr = (float)sqrt(fabs(H*H-K));
+        double dDiscr = (double)sqrt(fabs(H*H-K));
         rfCurv0 = H - dDiscr;
         rfCurv1 = H + dDiscr;
 
@@ -411,7 +450,7 @@ public:
     }
 
     //+++++++++ Quadric +++++++++++++++++++++++++++++++++++++++
-    static float F  ( float x, float y, float z ) 
+    double F  ( double x, double y, double z ) 
     {
         return (dKoeff[0] + dKoeff[1]*x + dKoeff[2]*y + dKoeff[3]*z +
                 dKoeff[4]*x*x + dKoeff[5]*y*y + dKoeff[6]*z*z +
@@ -419,48 +458,48 @@ public:
     }
   
     //+++++++++ 1. derivations ++++++++++++++++++++++++++++++++
-    static float Fx ( float x, float y, float z )
+    double Fx ( double x, double y, double z )
     {
         return( dKoeff[1] + 2.0f*dKoeff[4]*x + dKoeff[7]*y + dKoeff[8]*z );
     }
-    static float Fy ( float x, float y, float z ) 
+    double Fy ( double x, double y, double z ) 
     {
         return( dKoeff[2] + 2.0f*dKoeff[5]*y + dKoeff[7]*x + dKoeff[9]*z );
     }
-    static float Fz ( float x, float y, float z ) 
+    double Fz ( double x, double y, double z ) 
     {
         return( dKoeff[3] + 2.0f*dKoeff[6]*z + dKoeff[8]*x + dKoeff[9]*y );
     }
 
     //+++++++++ 2. derivations ++++++++++++++++++++++++++++++++
-    static float Fxx( float x, float y, float z ) 
+    double Fxx( double x, double y, double z ) 
     {
         return( 2.0f*dKoeff[4] );
     }
-    static float Fxy( float x, float y, float z ) 
+    double Fxy( double x, double y, double z ) 
     {
         return( dKoeff[7] );
     }
-    static float Fxz( float x, float y, float z ) 
+    double Fxz( double x, double y, double z ) 
     {
         return( dKoeff[8] );
     }
-    static float Fyy( float x, float y, float z ) 
+    double Fyy( double x, double y, double z ) 
     {
         return( 2.0f*dKoeff[5] );
     }
-    static float Fyz( float x, float y, float z ) 
+    double Fyz( double x, double y, double z ) 
     {
         return( dKoeff[9] );
     }
-    static float Fzz( float x, float y, float z ) 
+    double Fzz( double x, double y, double z ) 
     {
         return( 2.0f*dKoeff[6] );
     }
    
 protected:
-    static float dKoeff[ 10 ];     /**< Koeffizienten der Quadric */
-    Wm4::ImplicitSurface<float> *pImplSurf;  /**< Zugriff auf die MGC-Bibliothek */
+    double dKoeff[ 10 ];     /**< Koeffizienten der Quadric */
+    Wm4::ImplicitSurface<double> *pImplSurf;  /**< Zugriff auf die MGC-Bibliothek */
 
 private:
     /**
diff --git a/src/Mod/Mesh/App/Core/Curvature.cpp b/src/Mod/Mesh/App/Core/Curvature.cpp
new file mode 100644
index 000000000..d0371da4b
--- /dev/null
+++ b/src/Mod/Mesh/App/Core/Curvature.cpp
@@ -0,0 +1,236 @@
+/***************************************************************************
+ *   Copyright (c) 2012 Imetric 3D GmbH                                    *
+ *                                                                         *
+ *   This file is part of the FreeCAD CAx development system.              *
+ *                                                                         *
+ *   This library is free software; you can redistribute it and/or         *
+ *   modify it under the terms of the GNU Library General Public           *
+ *   License as published by the Free Software Foundation; either          *
+ *   version 2 of the License, or (at your option) any later version.      *
+ *                                                                         *
+ *   This library  is distributed in the hope that it will be useful,      *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU Library General Public License for more details.                  *
+ *                                                                         *
+ *   You should have received a copy of the GNU Library General Public     *
+ *   License along with this library; see the file COPYING.LIB. If not,    *
+ *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
+ *   Suite 330, Boston, MA  02111-1307, USA                                *
+ *                                                                         *
+ ***************************************************************************/
+
+#include "PreCompiled.h"
+#ifndef _PreComp_
+# include <algorithm>
+#endif
+
+#include <QFuture>
+#include <QFutureWatcher>
+#include <QtConcurrentMap>
+#include <boost/bind.hpp>
+
+#include <Mod/Mesh/App/WildMagic4/Wm4Vector3.h>
+#include <Mod/Mesh/App/WildMagic4/Wm4MeshCurvature.h>
+
+#include "Curvature.h"
+#include "Algorithm.h"
+#include "Approximation.h"
+#include "MeshKernel.h"
+#include "Iterator.h"
+#include "Tools.h"
+#include <Base/Sequencer.h>
+#include <Base/Tools.h>
+
+using namespace MeshCore;
+
+MeshCurvature::MeshCurvature(const MeshKernel& kernel)
+  : myKernel(kernel), myMinPoints(20), myRadius(0.5f)
+{
+    mySegment.resize(kernel.CountFacets());
+    std::generate(mySegment.begin(), mySegment.end(), Base::iotaGen<unsigned long>(0));
+}
+
+MeshCurvature::MeshCurvature(const MeshKernel& kernel, const std::vector<unsigned long>& segm)
+  : myKernel(kernel), myMinPoints(20), myRadius(0.5f), mySegment(segm)
+{
+}
+
+void MeshCurvature::ComputePerFace(bool parallel)
+{
+    Base::Vector3f rkDir0, rkDir1, rkPnt;
+    Base::Vector3f rkNormal;
+    myCurvature.clear();
+    MeshRefPointToFacets search(myKernel);
+    FacetCurvature face(myKernel, search, myRadius, myMinPoints);
+
+    if (!parallel) {
+        Base::SequencerLauncher seq("Curvature estimation", mySegment.size());
+        for (std::vector<unsigned long>::iterator it = mySegment.begin(); it != mySegment.end(); ++it) {
+            CurvatureInfo info = face.Compute(*it);
+            myCurvature.push_back(info);
+            seq.next();
+        }
+    }
+    else {
+        QFuture<CurvatureInfo> future = QtConcurrent::mapped
+            (mySegment, boost::bind(&FacetCurvature::Compute, &face, _1));
+        QFutureWatcher<CurvatureInfo> watcher;
+        watcher.setFuture(future);
+        watcher.waitForFinished();
+        for (QFuture<CurvatureInfo>::const_iterator it = future.begin(); it != future.end(); ++it) {
+            myCurvature.push_back(*it);
+        }
+    }
+}
+
+void MeshCurvature::ComputePerVertex()
+{
+    myCurvature.clear();
+
+    // get all points
+    std::vector< Wm4::Vector3<double> > aPnts;
+    aPnts.reserve(myKernel.CountPoints());
+    MeshPointIterator cPIt(myKernel);
+    for (cPIt.Init(); cPIt.More(); cPIt.Next()) {
+        Wm4::Vector3<double> cP(cPIt->x, cPIt->y, cPIt->z);
+        aPnts.push_back(cP);
+    }
+
+    // get all point connections
+    std::vector<int> aIdx;
+    aIdx.reserve(3*myKernel.CountFacets());
+    const MeshFacetArray& raFts = myKernel.GetFacets();
+    for (MeshFacetArray::const_iterator jt = raFts.begin(); jt != raFts.end(); ++jt) {
+        for (int i=0; i<3; i++) {
+            aIdx.push_back((int)jt->_aulPoints[i]);
+        }
+    }
+
+    // compute vertex based curvatures
+    Wm4::MeshCurvature<double> meshCurv(myKernel.CountPoints(), &(aPnts[0]), myKernel.CountFacets(), &(aIdx[0]));
+
+    // get curvature information now
+    const Wm4::Vector3<double>* aMaxCurvDir = meshCurv.GetMaxDirections();
+    const Wm4::Vector3<double>* aMinCurvDir = meshCurv.GetMinDirections();
+    const double* aMaxCurv = meshCurv.GetMaxCurvatures();
+    const double* aMinCurv = meshCurv.GetMinCurvatures();
+
+    myCurvature.reserve(myKernel.CountPoints());
+    for (unsigned long i=0; i<myKernel.CountPoints(); i++) {
+        CurvatureInfo ci;
+        ci.cMaxCurvDir = Base::Vector3f((float)aMaxCurvDir[i].X(), (float)aMaxCurvDir[i].Y(), (float)aMaxCurvDir[i].Z());
+        ci.cMinCurvDir = Base::Vector3f((float)aMinCurvDir[i].X(), (float)aMinCurvDir[i].Y(), (float)aMinCurvDir[i].Z());
+        ci.fMaxCurvature = (float)aMaxCurv[i];
+        ci.fMinCurvature = (float)aMinCurv[i];
+        myCurvature.push_back(ci);
+    }
+}
+
+// --------------------------------------------------------
+
+namespace MeshCore {
+class FitPointCollector : public MeshCollector
+{
+public:
+    FitPointCollector(std::set<unsigned long>& ind) : indices(ind){}
+    virtual void Append(const MeshCore::MeshKernel& kernel, unsigned long index)
+    {
+        unsigned long ulP1, ulP2, ulP3;
+        kernel.GetFacetPoints(index, ulP1, ulP2, ulP3);
+        indices.insert(ulP1);
+        indices.insert(ulP2);
+        indices.insert(ulP3);
+    }
+
+private:
+    std::set<unsigned long>& indices;
+};
+}
+
+// --------------------------------------------------------
+
+FacetCurvature::FacetCurvature(const MeshKernel& kernel, const MeshRefPointToFacets& search, float r, unsigned long pt)
+  : myKernel(kernel), mySearch(search), myRadius(r), myMinPoints(pt)
+{
+}
+
+CurvatureInfo FacetCurvature::Compute(unsigned long index) const
+{
+    Base::Vector3f rkDir0, rkDir1, rkPnt;
+    Base::Vector3f rkNormal;
+
+    MeshGeomFacet face = myKernel.GetFacet(index);
+    Base::Vector3f face_gravity = face.GetGravityPoint();
+    Base::Vector3f face_normal = face.GetNormal();
+    std::set<unsigned long> point_indices;
+    FitPointCollector collect(point_indices);
+
+    float searchDist = myRadius;
+    int attempts=0;
+    do {
+        mySearch.Neighbours(index, searchDist, collect);
+        if (point_indices.empty())
+            break;
+        float min_points = myMinPoints;
+        float use_points = point_indices.size();
+        searchDist = searchDist * sqrt(min_points/use_points);
+    }
+    while((point_indices.size() < myMinPoints) && (attempts++ < 3));
+
+    std::vector<Base::Vector3f> fitPoints;
+    const MeshPointArray& verts = myKernel.GetPoints();
+    fitPoints.reserve(point_indices.size());
+    for (std::set<unsigned long>::iterator it = point_indices.begin(); it != point_indices.end(); ++it) {
+        fitPoints.push_back(verts[*it] - face_gravity);
+    }
+
+    float fMin, fMax;
+    if (fitPoints.size() >= myMinPoints) {
+        SurfaceFit surf_fit;
+        surf_fit.AddPoints(fitPoints);
+        surf_fit.Fit();
+        rkNormal = surf_fit.GetNormal();
+        double dMin, dMax, dDistance;
+        if (surf_fit.GetCurvatureInfo(0.0, 0.0, 0.0, dMin, dMax, rkDir1, rkDir0, dDistance)) {
+            fMin = (float)dMin;
+            fMax = (float)dMax;
+        }
+        else {
+            fMin = FLT_MAX;
+            fMax = FLT_MAX;
+        }
+    }
+    else {
+        // too few points => cannot calc any properties
+        fMin = FLT_MAX;
+        fMax = FLT_MAX;
+    }
+
+    CurvatureInfo info;
+    if (fMin < fMax) {
+        info.fMaxCurvature = fMax;
+        info.fMinCurvature = fMin;
+        info.cMaxCurvDir = rkDir1;
+        info.cMinCurvDir = rkDir0;
+    }
+    else {
+        info.fMaxCurvature = fMin;
+        info.fMinCurvature = fMax;
+        info.cMaxCurvDir = rkDir0;
+        info.cMinCurvDir = rkDir1;
+    }
+
+    // Reverse the direction of the normal vector if required
+    // (Z component of "local" normal vectors should be opposite in sign to the "local" view vector)
+    if (rkNormal * face_normal < 0.0) {
+        // Note: Changing the normal directions is similar to flipping over the object.
+        // In this case we must adjust the curvature information as well.
+        std::swap(info.cMaxCurvDir,info.cMinCurvDir);
+        std::swap(info.fMaxCurvature,info.fMinCurvature);
+        info.fMaxCurvature *= (-1.0);
+        info.fMinCurvature *= (-1.0);
+    }
+
+    return info;
+}
diff --git a/src/Mod/Mesh/App/Core/Curvature.h b/src/Mod/Mesh/App/Core/Curvature.h
new file mode 100644
index 000000000..7c200946a
--- /dev/null
+++ b/src/Mod/Mesh/App/Core/Curvature.h
@@ -0,0 +1,75 @@
+/***************************************************************************
+ *   Copyright (c) 2012 Imetric 3D GmbH                                    *
+ *                                                                         *
+ *   This file is part of the FreeCAD CAx development system.              *
+ *                                                                         *
+ *   This library is free software; you can redistribute it and/or         *
+ *   modify it under the terms of the GNU Library General Public           *
+ *   License as published by the Free Software Foundation; either          *
+ *   version 2 of the License, or (at your option) any later version.      *
+ *                                                                         *
+ *   This library  is distributed in the hope that it will be useful,      *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU Library General Public License for more details.                  *
+ *                                                                         *
+ *   You should have received a copy of the GNU Library General Public     *
+ *   License along with this library; see the file COPYING.LIB. If not,    *
+ *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
+ *   Suite 330, Boston, MA  02111-1307, USA                                *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef MESHCORE_CURVATURE_H
+#define MESHCORE_CURVATURE_H
+
+#include <vector>
+#include <Base/Vector3D.h>
+
+namespace MeshCore {
+
+class MeshKernel;
+class MeshRefPointToFacets;
+
+/** Curvature information. */
+struct MeshExport CurvatureInfo
+{
+    float fMaxCurvature, fMinCurvature;
+    Base::Vector3f cMaxCurvDir, cMinCurvDir;
+};
+
+class MeshExport FacetCurvature
+{
+public:
+    FacetCurvature(const MeshKernel& kernel, const MeshRefPointToFacets& search, float, unsigned long);
+    CurvatureInfo Compute(unsigned long index) const;
+
+private:
+    const MeshKernel& myKernel;
+    const MeshRefPointToFacets& mySearch;
+    unsigned long myMinPoints;
+    float myRadius;
+};
+
+class MeshExport MeshCurvature
+{
+public:
+    MeshCurvature(const MeshKernel& kernel);
+    MeshCurvature(const MeshKernel& kernel, const std::vector<unsigned long>& segm);
+    float GetRadius() const { return myRadius; }
+    void SetRadius(float r) { myRadius = r; }
+    void ComputePerFace(bool parallel);
+    void ComputePerVertex();
+    const std::vector<CurvatureInfo>& GetCurvature() const { return myCurvature; }
+
+private:
+    const MeshKernel& myKernel;
+    unsigned long myMinPoints;
+    float myRadius;
+    std::vector<unsigned long> mySegment;
+    std::vector<CurvatureInfo> myCurvature;
+};
+
+} // MeshCore
+
+#endif // MESHCORE_CURVATURE_H
diff --git a/src/Mod/Mesh/App/Core/Segmentation.cpp b/src/Mod/Mesh/App/Core/Segmentation.cpp
new file mode 100644
index 000000000..4326457e5
--- /dev/null
+++ b/src/Mod/Mesh/App/Core/Segmentation.cpp
@@ -0,0 +1,228 @@
+/***************************************************************************
+ *   Copyright (c) 2012 Imetric 3D GmbH                                    *
+ *                                                                         *
+ *   This file is part of the FreeCAD CAx development system.              *
+ *                                                                         *
+ *   This library is free software; you can redistribute it and/or         *
+ *   modify it under the terms of the GNU Library General Public           *
+ *   License as published by the Free Software Foundation; either          *
+ *   version 2 of the License, or (at your option) any later version.      *
+ *                                                                         *
+ *   This library  is distributed in the hope that it will be useful,      *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU Library General Public License for more details.                  *
+ *                                                                         *
+ *   You should have received a copy of the GNU Library General Public     *
+ *   License along with this library; see the file COPYING.LIB. If not,    *
+ *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
+ *   Suite 330, Boston, MA  02111-1307, USA                                *
+ *                                                                         *
+ ***************************************************************************/
+
+#include "PreCompiled.h"
+#ifndef _PreComp_
+#include <algorithm>
+#endif
+
+#include "Segmentation.h"
+#include "Algorithm.h"
+#include "Approximation.h"
+
+using namespace MeshCore;
+
+void MeshSurfaceSegment::PrepareFacet(unsigned long)
+{
+}
+
+void MeshSurfaceSegment::AddSegment(const std::vector<unsigned long>& segm)
+{
+    if (segm.size() >= minFacets) {
+        segments.push_back(segm);
+    }
+}
+
+// --------------------------------------------------------
+
+MeshDistancePlanarSegment::MeshDistancePlanarSegment(const MeshKernel& mesh, unsigned long minFacets, float tol)
+  : MeshDistanceSurfaceSegment(mesh, minFacets, tol), fitter(new PlaneFit)
+{
+}
+
+MeshDistancePlanarSegment::~MeshDistancePlanarSegment()
+{
+    delete fitter;
+}
+
+void MeshDistancePlanarSegment::PrepareFacet(unsigned long index)
+{
+    fitter->Clear();
+
+    MeshGeomFacet triangle = kernel.GetFacet(index);
+    basepoint = triangle.GetGravityPoint();
+    normal = triangle.GetNormal();
+    fitter->AddPoint(triangle._aclPoints[0]);
+    fitter->AddPoint(triangle._aclPoints[1]);
+    fitter->AddPoint(triangle._aclPoints[2]);
+}
+
+bool MeshDistancePlanarSegment::TestFacet (const MeshFacet& face) const
+{
+    if (!fitter->Done())
+        fitter->Fit();
+    MeshGeomFacet triangle = kernel.GetFacet(face);
+    for (int i=0; i<3; i++) {
+        if (fabs(fitter->GetDistanceToPlane(triangle._aclPoints[i])) > tolerance)
+            return false;
+    }
+
+    fitter->AddPoint(triangle.GetGravityPoint());
+    return true;
+}
+
+// --------------------------------------------------------
+
+bool MeshCurvaturePlanarSegment::TestFacet (const MeshFacet &rclFacet) const
+{
+    for (int i=0; i<3; i++) {
+        const CurvatureInfo& ci = info[rclFacet._aulPoints[i]];
+        if (fabs(ci.fMinCurvature) > tolerance)
+            return false;
+        if (fabs(ci.fMaxCurvature) > tolerance)
+            return false;
+    }
+
+    return true;
+}
+
+bool MeshCurvatureCylindricalSegment::TestFacet (const MeshFacet &rclFacet) const
+{
+    for (int i=0; i<3; i++) {
+        const CurvatureInfo& ci = info[rclFacet._aulPoints[i]];
+        if (ci.fMaxCurvature > ci.fMinCurvature) {
+            // convexe
+            if (fabs(ci.fMinCurvature) > tolerance)
+                return false;
+            float diff = ci.fMaxCurvature - curvature;
+            if (fabs(diff) > tolerance)
+                return false;
+        }
+        else {
+            // concave
+            if (fabs(ci.fMaxCurvature) > tolerance)
+                return false;
+            float diff = ci.fMinCurvature + curvature;
+            if (fabs(diff) > tolerance)
+                return false;
+        }
+    }
+
+    return true;
+}
+
+bool MeshCurvatureSphericalSegment::TestFacet (const MeshFacet &rclFacet) const
+{
+    for (int i=0; i<3; i++) {
+        const CurvatureInfo& ci = info[rclFacet._aulPoints[i]];
+        if (ci.fMaxCurvature * ci.fMinCurvature < 0)
+            return false;
+        float diff;
+        diff = fabs(ci.fMinCurvature) - curvature;
+        if (fabs(diff) > tolerance)
+            return false;
+        diff = fabs(ci.fMaxCurvature) - curvature;
+        if (fabs(diff) > tolerance)
+            return false;
+    }
+
+    return true;
+}
+
+bool MeshCurvatureFreeformSegment::TestFacet (const MeshFacet &rclFacet) const
+{
+    for (int i=0; i<3; i++) {
+        const CurvatureInfo& ci = info[rclFacet._aulPoints[i]];
+        if (fabs(ci.fMinCurvature-c2) > tolerance)
+            return false;
+        if (fabs(ci.fMaxCurvature-c1) > tolerance)
+            return false;
+    }
+
+    return true;
+}
+
+// --------------------------------------------------------
+
+MeshSurfaceVisitor::MeshSurfaceVisitor (const MeshSurfaceSegment& segm, std::vector<unsigned long> &indices)
+  : indices(indices), segm(segm)
+{
+}
+
+MeshSurfaceVisitor::~MeshSurfaceVisitor ()
+{
+}
+
+bool MeshSurfaceVisitor::AllowVisit (const MeshFacet& face, const MeshFacet&, 
+                                     unsigned long, unsigned long, unsigned short)
+{
+    return segm.TestFacet(face);
+}
+
+bool MeshSurfaceVisitor::Visit (const MeshFacet & face, const MeshFacet &,
+                                unsigned long ulFInd, unsigned long)
+{
+    indices.push_back(ulFInd);
+    return true;
+}
+
+// --------------------------------------------------------
+
+void MeshSegmentAlgorithm::FindSegments(std::vector<MeshSurfaceSegment*>& segm)
+{
+    // reset VISIT flags
+    unsigned long startFacet;
+    MeshCore::MeshAlgorithm cAlgo(myKernel);
+    cAlgo.ResetFacetFlag(MeshCore::MeshFacet::VISIT);
+
+    const MeshCore::MeshFacetArray& rFAry = myKernel.GetFacets();
+    MeshCore::MeshFacetArray::_TConstIterator iCur = rFAry.begin();
+    MeshCore::MeshFacetArray::_TConstIterator iBeg = rFAry.begin();
+    MeshCore::MeshFacetArray::_TConstIterator iEnd = rFAry.end();
+
+    // start from the first not visited facet
+    cAlgo.CountFacetFlag(MeshCore::MeshFacet::VISIT);
+    std::vector<unsigned long> resetVisited;
+
+    for (std::vector<MeshSurfaceSegment*>::iterator it = segm.begin(); it != segm.end(); ++it) {
+        cAlgo.ResetFacetsFlag(resetVisited, MeshCore::MeshFacet::VISIT);
+        resetVisited.clear();
+
+        iCur = std::find_if(iBeg, iEnd, std::bind2nd(MeshCore::MeshIsNotFlag<MeshCore::MeshFacet>(),
+            MeshCore::MeshFacet::VISIT));
+        startFacet = iCur - iBeg;
+        while (startFacet != ULONG_MAX) {
+            // collect all facets of the same geometry
+            std::vector<unsigned long> indices;
+            indices.push_back(startFacet);
+            (*it)->PrepareFacet(startFacet);
+            MeshSurfaceVisitor pv(**it, indices);
+            myKernel.VisitNeighbourFacets(pv, startFacet);
+
+            // add or discard the segment
+            if (indices.size() == 1) {
+                resetVisited.push_back(startFacet);
+            }
+            else {
+                (*it)->AddSegment(indices);
+            }
+
+            // search for the next start facet
+            iCur = std::find_if(iCur, iEnd, std::bind2nd(MeshCore::MeshIsNotFlag<MeshCore::MeshFacet>(),
+                MeshCore::MeshFacet::VISIT));
+            if (iCur < iEnd)
+                startFacet = iCur - iBeg;
+            else
+                startFacet = ULONG_MAX;
+        }
+    }
+}
diff --git a/src/Mod/Mesh/App/Core/Segmentation.h b/src/Mod/Mesh/App/Core/Segmentation.h
new file mode 100644
index 000000000..e45ecff4c
--- /dev/null
+++ b/src/Mod/Mesh/App/Core/Segmentation.h
@@ -0,0 +1,161 @@
+/***************************************************************************
+ *   Copyright (c) 2012 Imetric 3D GmbH                                    *
+ *                                                                         *
+ *   This file is part of the FreeCAD CAx development system.              *
+ *                                                                         *
+ *   This library is free software; you can redistribute it and/or         *
+ *   modify it under the terms of the GNU Library General Public           *
+ *   License as published by the Free Software Foundation; either          *
+ *   version 2 of the License, or (at your option) any later version.      *
+ *                                                                         *
+ *   This library  is distributed in the hope that it will be useful,      *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU Library General Public License for more details.                  *
+ *                                                                         *
+ *   You should have received a copy of the GNU Library General Public     *
+ *   License along with this library; see the file COPYING.LIB. If not,    *
+ *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
+ *   Suite 330, Boston, MA  02111-1307, USA                                *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef MESHCORE_SEGMENTATION_H
+#define MESHCORE_SEGMENTATION_H
+
+#include "MeshKernel.h"
+#include "Curvature.h"
+#include "Visitor.h"
+#include <vector>
+
+namespace MeshCore {
+
+class PlaneFit;
+class MeshFacet;
+typedef std::vector<unsigned long> MeshSegment;
+
+class MeshExport MeshSurfaceSegment
+{
+public:
+    MeshSurfaceSegment(unsigned long minFacets)
+        : minFacets(minFacets) {}
+    virtual ~MeshSurfaceSegment() {}
+    virtual bool TestFacet (const MeshFacet &rclFacet) const = 0;
+    virtual void PrepareFacet(unsigned long);
+    void AddSegment(const std::vector<unsigned long>&);
+    const std::vector<MeshSegment> GetSegments() const { return segments; }
+
+protected:
+    std::vector<MeshSegment> segments;
+    unsigned long minFacets;
+};
+
+// --------------------------------------------------------
+
+class MeshExport MeshDistanceSurfaceSegment : public MeshSurfaceSegment
+{
+public:
+    MeshDistanceSurfaceSegment(const MeshKernel& mesh, unsigned long minFacets, float tol)
+        : MeshSurfaceSegment(minFacets), kernel(mesh), tolerance(tol) {}
+
+protected:
+    const MeshKernel& kernel;
+    float tolerance;
+};
+
+class MeshExport MeshDistancePlanarSegment : public MeshDistanceSurfaceSegment
+{
+public:
+    MeshDistancePlanarSegment(const MeshKernel& mesh, unsigned long minFacets, float tol);
+    virtual ~MeshDistancePlanarSegment();
+    bool TestFacet (const MeshFacet &rclFacet) const;
+    void PrepareFacet(unsigned long);
+
+protected:
+    Base::Vector3f basepoint;
+    Base::Vector3f normal;
+    PlaneFit* fitter;
+};
+
+// --------------------------------------------------------
+
+class MeshExport MeshCurvatureSurfaceSegment : public MeshSurfaceSegment
+{
+public:
+    MeshCurvatureSurfaceSegment(const std::vector<CurvatureInfo>& ci, unsigned long minFacets, float tol)
+        : MeshSurfaceSegment(minFacets), info(ci), tolerance(tol) {}
+
+protected:
+    const std::vector<CurvatureInfo>& info;
+    float tolerance;
+};
+
+class MeshExport MeshCurvaturePlanarSegment : public MeshCurvatureSurfaceSegment
+{
+public:
+    MeshCurvaturePlanarSegment(const std::vector<CurvatureInfo>& ci, unsigned long minFacets, float tol)
+        : MeshCurvatureSurfaceSegment(ci, minFacets, tol) {}
+    virtual bool TestFacet (const MeshFacet &rclFacet) const;
+};
+
+class MeshExport MeshCurvatureCylindricalSegment : public MeshCurvatureSurfaceSegment
+{
+public:
+    MeshCurvatureCylindricalSegment(const std::vector<CurvatureInfo>& ci, unsigned long minFacets, float tol, float radius)
+        : MeshCurvatureSurfaceSegment(ci, minFacets, tol) { curvature = 1/radius;}
+    virtual bool TestFacet (const MeshFacet &rclFacet) const;
+
+private:
+    float curvature;
+};
+
+class MeshExport MeshCurvatureSphericalSegment : public MeshCurvatureSurfaceSegment
+{
+public:
+    MeshCurvatureSphericalSegment(const std::vector<CurvatureInfo>& ci, unsigned long minFacets, float tol, float radius)
+        : MeshCurvatureSurfaceSegment(ci, minFacets, tol) { curvature = 1/radius;}
+    virtual bool TestFacet (const MeshFacet &rclFacet) const;
+
+private:
+    float curvature;
+};
+
+class MeshExport MeshCurvatureFreeformSegment : public MeshCurvatureSurfaceSegment
+{
+public:
+    MeshCurvatureFreeformSegment(const std::vector<CurvatureInfo>& ci, unsigned long minFacets, float tol, float c1, float c2)
+        : MeshCurvatureSurfaceSegment(ci, minFacets, tol), c1(c1), c2(c2) {}
+    virtual bool TestFacet (const MeshFacet &rclFacet) const;
+
+private:
+    float c1, c2;
+};
+
+class MeshExport MeshSurfaceVisitor : public MeshFacetVisitor
+{
+public:
+    MeshSurfaceVisitor (const MeshSurfaceSegment& segm, std::vector<unsigned long> &indices);
+    virtual ~MeshSurfaceVisitor ();
+    bool AllowVisit (const MeshFacet& face, const MeshFacet&, 
+                     unsigned long, unsigned long, unsigned short neighbourIndex);
+    bool Visit (const MeshFacet & face, const MeshFacet &,
+                unsigned long ulFInd, unsigned long);
+
+protected:
+    std::vector<unsigned long>  &indices;
+    const MeshSurfaceSegment& segm;
+};
+
+class MeshExport MeshSegmentAlgorithm
+{
+public:
+    MeshSegmentAlgorithm(const MeshKernel& kernel) : myKernel(kernel) {}
+    void FindSegments(std::vector<MeshSurfaceSegment*>&);
+
+private:
+    const MeshKernel& myKernel;
+};
+
+} // MeshCore
+
+#endif // MESHCORE_SEGMENTATION_H
diff --git a/src/Mod/Mesh/App/FeatureMeshCurvature.cpp b/src/Mod/Mesh/App/FeatureMeshCurvature.cpp
index f618f2dea..f7bc01205 100644
--- a/src/Mod/Mesh/App/FeatureMeshCurvature.cpp
+++ b/src/Mod/Mesh/App/FeatureMeshCurvature.cpp
@@ -29,18 +29,16 @@
 #include <Base/Exception.h>
 #include <Base/Matrix.h>
 #include <Base/Sequencer.h>
-#include <Mod/Mesh/App/WildMagic4/Wm4Vector3.h>
-#include <Mod/Mesh/App/WildMagic4/Wm4MeshCurvature.h>
 #include "FeatureMeshCurvature.h"
 #include "MeshFeature.h"
 
+#include "Core/Curvature.h"
 #include "Core/Elements.h"
 #include "Core/Iterator.h"
 
 
 
 using namespace Mesh;
-using namespace MeshCore;
 
 PROPERTY_SOURCE(Mesh::Curvature, App::DocumentObject)
 
@@ -68,42 +66,20 @@ App::DocumentObjectExecReturn *Curvature::execute(void)
     }
  
     // get all points
-    const MeshKernel& rMesh = pcFeat->Mesh.getValue().getKernel();
-    std::vector< Wm4::Vector3<float> > aPnts;
-    aPnts.reserve(rMesh.CountPoints());
-    MeshPointIterator cPIt( rMesh );
-    for (cPIt.Init(); cPIt.More(); cPIt.Next()) {
-        Wm4::Vector3<float> cP(cPIt->x, cPIt->y, cPIt->z);
-        aPnts.push_back(cP);
-    }
+    const MeshCore::MeshKernel& rMesh = pcFeat->Mesh.getValue().getKernel();
+    MeshCore::MeshCurvature meshCurv(rMesh);
+    meshCurv.ComputePerVertex();
+    const std::vector<MeshCore::CurvatureInfo>& curv = meshCurv.GetCurvature();
 
-    // get all point connections
-    std::vector<int> aIdx;
-    aIdx.reserve(3*rMesh.CountFacets());
-    const std::vector<MeshFacet>& raFts = rMesh.GetFacets();
-    for (std::vector<MeshFacet>::const_iterator jt = raFts.begin(); jt != raFts.end(); ++jt) {
-        for (int i=0; i<3; i++) {
-            aIdx.push_back((int)jt->_aulPoints[i]);
-        }
-    }
-
-    // compute vertex based curvatures
-    Wm4::MeshCurvature<float> meshCurv(rMesh.CountPoints(), &(aPnts[0]), rMesh.CountFacets(), &(aIdx[0]));
-
-    // get curvature information now
-    const Wm4::Vector3<float>* aMaxCurvDir = meshCurv.GetMaxDirections();
-    const Wm4::Vector3<float>* aMinCurvDir = meshCurv.GetMinDirections();
-    const float* aMaxCurv = meshCurv.GetMaxCurvatures();
-    const float* aMinCurv = meshCurv.GetMinCurvatures();
-
-    std::vector<CurvatureInfo> values(rMesh.CountPoints());
-    for (unsigned long i=0; i<rMesh.CountPoints(); i++) {
+    std::vector<CurvatureInfo> values;
+    values.reserve(curv.size());
+    for (std::vector<MeshCore::CurvatureInfo>::const_iterator it = curv.begin(); it != curv.end(); ++it) {
         CurvatureInfo ci;
-        ci.cMaxCurvDir = Base::Vector3f(aMaxCurvDir[i].X(), aMaxCurvDir[i].Y(), aMaxCurvDir[i].Z());
-        ci.cMinCurvDir = Base::Vector3f(aMinCurvDir[i].X(), aMinCurvDir[i].Y(), aMinCurvDir[i].Z());
-        ci.fMaxCurvature = aMaxCurv[i];
-        ci.fMinCurvature = aMinCurv[i];
-        values[i] = ci;
+        ci.cMaxCurvDir = it->cMaxCurvDir;
+        ci.cMinCurvDir = it->cMinCurvDir;
+        ci.fMaxCurvature = it->fMaxCurvature;
+        ci.fMinCurvature = it->fMinCurvature;
+        values.push_back(ci);
     }
 
     CurvInfo.setValues(values);
diff --git a/src/Mod/Mesh/App/Makefile.am b/src/Mod/Mesh/App/Makefile.am
index 6309f0d91..27537e060 100644
--- a/src/Mod/Mesh/App/Makefile.am
+++ b/src/Mod/Mesh/App/Makefile.am
@@ -22,6 +22,8 @@ libMesh_la_SOURCES=\
 		Core/Approximation.h \
 		Core/Builder.cpp \
 		Core/Builder.h \
+		Core/Curvature.cpp \
+		Core/Curvature.h \
 		Core/Definitions.cpp \
 		Core/Definitions.h \
 		Core/Degeneration.cpp \
@@ -42,6 +44,8 @@ libMesh_la_SOURCES=\
 		Core/MeshIO.h \
 		Core/Projection.cpp \
 		Core/Projection.h \
+		Core/Segmentation.cpp \
+		Core/Segmentation.h \
 		Core/SetOperations.cpp \
 		Core/SetOperations.h \
 		Core/Smoothing.cpp \
@@ -324,9 +328,9 @@ nobase_include_HEADERS = \
 
 
 # the library search path.
-libMesh_la_LDFLAGS = -L../../../Base -L../../../App $(all_libraries) $(GTS_LIBS) \
+libMesh_la_LDFLAGS = -L../../../Base -L../../../App $(QT4_CORE_LIBS) $(all_libraries) $(GTS_LIBS) \
 		-version-info @LIB_CURRENT@:@LIB_REVISION@:@LIB_AGE@
-libMesh_la_CPPFLAGS = -DMeshExport= 
+libMesh_la_CPPFLAGS = -DMeshExport= -DEIGEN2_SUPPORT
 
 libMesh_la_LIBADD   = \
 		@BOOST_FILESYSTEM_LIB@ @BOOST_REGEX_LIB@ @BOOST_SYSTEM_LIB@ \
@@ -354,7 +358,8 @@ Mesh_la_DEPENDENCIES = libMesh.la
 #--------------------------------------------------------------------------------------
 
 # set the include path found by configure
-AM_CXXFLAGS = -I$(top_srcdir)/src/3rdParty -I$(top_srcdir)/src -I$(top_builddir)/src $(GTS_CFLAGS) $(all_includes) 
+AM_CXXFLAGS = -I$(top_srcdir)/src/3rdParty -I$(top_srcdir)/src -I$(top_builddir)/src $(GTS_CFLAGS) \
+		$(all_includes) $(QT4_CORE_CXXFLAGS)
 
 includedir = @includedir@/Mod/Mesh/App
 libdir = $(prefix)/Mod/Mesh
diff --git a/src/Mod/Mesh/App/Mesh.cpp b/src/Mod/Mesh/App/Mesh.cpp
index bd7be06d0..89f85856a 100644
--- a/src/Mod/Mesh/App/Mesh.cpp
+++ b/src/Mod/Mesh/App/Mesh.cpp
@@ -1422,10 +1422,11 @@ MeshObject* MeshObject::meshFromSegment(const std::vector<unsigned long>& indice
     return new MeshObject(kernel, _Mtrx);
 }
 
-std::vector<Segment> MeshObject::getSegmentsFromType(MeshObject::Type type, const Segment& aSegment, float dev) const
+std::vector<Segment> MeshObject::getSegmentsFromType(MeshObject::Type type, const Segment& aSegment,
+                                                     float dev, unsigned long minFacets) const
 {
     std::vector<Segment> segm;
-    unsigned long startFacet, visited;
+    unsigned long startFacet;
     if (this->_kernel.CountFacets() == 0)
         return segm;
 
@@ -1445,7 +1446,7 @@ std::vector<Segment> MeshObject::getSegmentsFromType(MeshObject::Type type, cons
     MeshCore::MeshFacetArray::_TConstIterator iEnd = rFAry.end();
 
     // start from the first not visited facet
-    visited = cAlgo.CountFacetFlag(MeshCore::MeshFacet::VISIT);
+    cAlgo.CountFacetFlag(MeshCore::MeshFacet::VISIT);
     iTri = std::find_if(iTri, iEnd, std::bind2nd(MeshCore::MeshIsNotFlag<MeshCore::MeshFacet>(),
         MeshCore::MeshFacet::VISIT));
     startFacet = iTri - iBeg;
@@ -1455,7 +1456,7 @@ std::vector<Segment> MeshObject::getSegmentsFromType(MeshObject::Type type, cons
         std::vector<unsigned long> indices;
         indices.push_back(startFacet);
         MeshCore::MeshPlaneVisitor pv(this->_kernel, startFacet, dev, indices);
-        visited += this->_kernel.VisitNeighbourFacets(pv, startFacet);
+        this->_kernel.VisitNeighbourFacets(pv, startFacet);
 
         iTri = std::find_if(iTri, iEnd, std::bind2nd(MeshCore::MeshIsNotFlag<MeshCore::MeshFacet>(),
             MeshCore::MeshFacet::VISIT));
@@ -1463,7 +1464,8 @@ std::vector<Segment> MeshObject::getSegmentsFromType(MeshObject::Type type, cons
             startFacet = iTri - iBeg;
         else
             startFacet = ULONG_MAX;
-        segm.push_back(Segment(const_cast<MeshObject*>(this), indices, false));
+        if (indices.size() > minFacets)
+            segm.push_back(Segment(const_cast<MeshObject*>(this), indices, false));
     }
 
     return segm;
diff --git a/src/Mod/Mesh/App/Mesh.h b/src/Mod/Mesh/App/Mesh.h
index dcd1ca19d..3d884517f 100644
--- a/src/Mod/Mesh/App/Mesh.h
+++ b/src/Mod/Mesh/App/Mesh.h
@@ -266,7 +266,7 @@ public:
     const Segment& getSegment(unsigned long) const;
     Segment& getSegment(unsigned long);
     MeshObject* meshFromSegment(const std::vector<unsigned long>&) const;
-    std::vector<Segment> getSegmentsFromType(Type, const Segment& aSegment, float dev) const;
+    std::vector<Segment> getSegmentsFromType(Type, const Segment& aSegment, float dev, unsigned long minFacets) const;
     //@}
 
     /** @name Primitives */
diff --git a/src/Mod/Mesh/App/MeshPy.xml b/src/Mod/Mesh/App/MeshPy.xml
index 6070d9fc1..f29a59730 100644
--- a/src/Mod/Mesh/App/MeshPy.xml
+++ b/src/Mod/Mesh/App/MeshPy.xml
@@ -370,13 +370,26 @@ an empty dictionary if there is no intersection.
 </UserDocu>
 			</Documentation>
 		</Methode>
-		<Methode Name="getPlanes" Const="true">
+		<Methode Name="getPlanarSegments" Const="true">
 			<Documentation>
-				<UserDocu>Get all planes of the mesh as segment.
+				<UserDocu>getPlanarSegments(dev,[min faces=0]) -> list
+Get all planes of the mesh as segment.
 In the worst case each triangle can be regarded as single
 plane if none of its neighours is coplanar.</UserDocu>
 			</Documentation>
 		</Methode>
+		<Methode Name="getSegmentsByCurvature" Const="true">
+			<Documentation>
+				<UserDocu>getSegmentsByCurvature(list) -> list
+The argument list gives a list if tuples where it defines the preferred maximum curvature,
+the preferred minumum curvature, the tolerance and the number of minimum faces for the segment.
+Example:
+c=(1.0, 0.0, 0.1, 500) # search for a cylinder with radius 1.0
+p=(0.0, 0.0, 0.1, 500) # search for a plane
+mesh.getSegmentsByCurvature([c,p])
+				</UserDocu>
+			</Documentation>
+		</Methode>
 		<Attribute Name="Points" ReadOnly="true">
 			<Documentation>
 				<UserDocu>A collection of the mesh points
diff --git a/src/Mod/Mesh/App/MeshPyImp.cpp b/src/Mod/Mesh/App/MeshPyImp.cpp
index 143ff61f9..1e32bf2f3 100644
--- a/src/Mod/Mesh/App/MeshPyImp.cpp
+++ b/src/Mod/Mesh/App/MeshPyImp.cpp
@@ -41,6 +41,8 @@
 #include "Core/Grid.h"
 #include "Core/MeshKernel.h"
 #include "Core/Triangulation.h"
+#include "Core/Segmentation.h"
+#include "Core/Curvature.h"
 
 using namespace Mesh;
 
@@ -1356,15 +1358,16 @@ PyObject* MeshPy::nearestFacetOnRay(PyObject *args)
     }
 }
 
-PyObject*  MeshPy::getPlanes(PyObject *args)
+PyObject*  MeshPy::getPlanarSegments(PyObject *args)
 {
     float dev;
-    if (!PyArg_ParseTuple(args, "f",&dev))
+    unsigned long minFacets=0;
+    if (!PyArg_ParseTuple(args, "f|k",&dev,&minFacets))
         return NULL;
 
     Mesh::MeshObject* mesh = getMeshObjectPtr();
     std::vector<Mesh::Segment> segments = mesh->getSegmentsFromType
-        (Mesh::MeshObject::PLANE, Mesh::Segment(mesh,false), dev);
+        (Mesh::MeshObject::PLANE, Mesh::Segment(mesh,false), dev, minFacets);
 
     Py::List s;
     for (std::vector<Mesh::Segment>::iterator it = segments.begin(); it != segments.end(); ++it) {
@@ -1379,6 +1382,48 @@ PyObject*  MeshPy::getPlanes(PyObject *args)
     return Py::new_reference_to(s);
 }
 
+PyObject*  MeshPy::getSegmentsByCurvature(PyObject *args)
+{
+    PyObject* l;
+    if (!PyArg_ParseTuple(args, "O!",&PyList_Type,&l))
+        return NULL;
+
+    const MeshCore::MeshKernel& kernel = getMeshObjectPtr()->getKernel();
+    MeshCore::MeshSegmentAlgorithm finder(kernel);
+    MeshCore::MeshCurvature meshCurv(kernel);
+    meshCurv.ComputePerVertex();
+
+    Py::List func(l);
+    std::vector<MeshCore::MeshSurfaceSegment*> segm;
+    //segm.push_back(new MeshCore::MeshCurvatureCylindricalSegment(meshCurv.GetCurvature(), minFacets, dev, 4.75f));
+    //segm.push_back(new MeshCore::MeshCurvaturePlanarSegment(meshCurv.GetCurvature(), minFacets, dev));
+    for (Py::List::iterator it = func.begin(); it != func.end(); ++it) {
+        Py::Tuple t(*it);
+        float c1 = (float)Py::Float(t[0]);
+        float c2 = (float)Py::Float(t[1]);
+        float tol = (float)Py::Float(t[2]);
+        int num = (int)Py::Int(t[3]);
+        segm.push_back(new MeshCore::MeshCurvatureFreeformSegment(meshCurv.GetCurvature(), num, tol, c1, c2));
+    }
+
+    finder.FindSegments(segm);
+
+    Py::List list;
+    for (std::vector<MeshCore::MeshSurfaceSegment*>::iterator segmIt = segm.begin(); segmIt != segm.end(); ++segmIt) {
+        std::vector<MeshCore::MeshSegment> data = (*segmIt)->GetSegments();
+        delete (*segmIt);
+        for (std::vector<MeshCore::MeshSegment>::iterator it = data.begin(); it != data.end(); ++it) {
+            Py::List ary;
+            for (MeshCore::MeshSegment::const_iterator jt = it->begin(); jt != it->end(); ++jt) {
+                ary.append(Py::Int((int)*jt));
+            }
+            list.append(ary);
+        }
+    }
+
+    return Py::new_reference_to(list);
+}
+
 Py::Int MeshPy::getCountPoints(void) const
 {
     return Py::Int((long)getMeshObjectPtr()->countPoints());