SpringbackCorrection.cpp

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (c) 2007                                                    *
 *   Joachim Zettler <Joachim.Zettler@gmx.de>                  *
 *   Human Rezai <Human@web.de>                                            *
 *   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"

#include "SpringbackCorrection.h"
#include "best_fit.h"


#include <Mod/Mesh/App/Core/Builder.h>
#include <Mod/Mesh/App/WildMagic4/Wm4MeshCurvature.h>

#include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
#include <Base/Builder3D.h>
#include <GeomAPI_ProjectPointOnSurf.hxx>
#include <BRepTools.hxx>
#include <BRep_Builder.hxx>
#include <Poly_PolygonOnTriangulation.hxx>
#include <Poly_Triangulation.hxx>
#include <BRepExtrema_DistShapeShape.hxx>
#include <TopExp_Explorer.hxx>
#include <TopExp.hxx>
#include <TopoDS_Vertex.hxx>
#include <BRep_Tool.hxx>
#include <TColgp_Array1OfPnt2d.hxx>
#include <BRepAdaptor_Surface.hxx>
#include <BRepExtrema_DistShapeShape.hxx>
#include <TopoDS.hxx>
#include <TopTools_ListIteratorOfListOfShape.hxx>
#include <Geom_Surface.hxx>




#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
#include <boost/numeric/bindings/atlas/clapack.hpp>

using namespace boost::numeric::bindings;

using namespace boost::numeric;


SpringbackCorrection::SpringbackCorrection()
{
}

SpringbackCorrection::SpringbackCorrection(const TopoDS_Shape& aShape, const MeshCore::MeshKernel& aMesh)
        :m_Shape(aShape),m_Mesh(aMesh)
{
    m_EdgeStruct.clear();
    EdgeMap.clear();
    MeshMap.clear();
    //Remove any existing Triangulation on the Shape
    BRepTools::Clean(m_Shape);
    best_fit::Tesselate_Shape(m_Shape,m_CadMesh,(float) 0.1); // Basistriangulierung des ganzen Shapes

    MeshCore::MeshTopoAlgorithm algo(m_CadMesh);
    algo.HarmonizeNormals();

        // flip all normals so that they are oriented in positive z-direction 
    MeshCore::MeshGeomFacet facet = m_CadMesh.GetFacet(0);
    if (facet.GetNormal().z < 0.0)
        algo.FlipNormals();

    int n = m_CadMesh.CountFacets();
    Base::Vector3f normal;

    MeshCore::MeshPointArray points = m_CadMesh.GetPoints();
    MeshCore::MeshFacetArray facets = m_CadMesh.GetFacets();

    for (int i=0; i<n; ++i)
    {
        MeshCore::MeshGeomFacet face = m_CadMesh.GetFacet(i);
        normal = face.GetNormal();

        if (normal.z < 0.0)
        {
            facets[i].FlipNormal();
        }
    }

    m_CadMesh.Assign(points, facets);
}

SpringbackCorrection::~SpringbackCorrection()
{
}
bool SpringbackCorrection::Load(const MeshCore::MeshKernel& aMesh)
{
        m_MeshVec.push_back(aMesh);
        //m_CadMesh = MeshRef;

        return true;
}

bool SpringbackCorrection::Load(const TopoDS_Shape& aShape)
{
        m_Shape = aShape;
        return true;
}
bool SpringbackCorrection::Load(const TopoDS_Shape& aShape, const MeshCore::MeshKernel& aMesh)
{
    m_Shape = aShape;
    m_Mesh  = aMesh;

    m_EdgeStruct.clear();
    EdgeMap.clear();
    MeshMap.clear();

    //Remove any existing Triangulation on the Shape
    BRepTools::Clean(m_Shape);
    best_fit::Tesselate_Shape(m_Shape,m_CadMesh,(float) 0.1); // Basistriangulierung des ganzen Shapes

    MeshCore::MeshTopoAlgorithm algo(m_CadMesh);
    algo.HarmonizeNormals();

    MeshCore::MeshGeomFacet facet = m_CadMesh.GetFacet(0);
    if (facet.GetNormal().z < 0.0)
        algo.FlipNormals();

    int n = m_CadMesh.CountFacets();
    Base::Vector3f normal;

    //MeshCore::MeshPointArray points = m_CadMesh.GetPoints();
    //MeshCore::MeshFacetArray facets = m_CadMesh.GetFacets();

    //for (int i=0; i<n; ++i)
    //{
    //    MeshCore::MeshGeomFacet face = m_CadMesh.GetFacet(i);
    //    normal = face.GetNormal();

    //    if (normal.z < 0.0)
    //    {
    //        facets[i].FlipNormal();
    //    }
    //}

    //m_CadMesh.Assign(points, facets);

    return true;
}

bool SpringbackCorrection::CalcCurv()
{
    Base::Builder3D logo;
    BRep_Builder VertexBuild;
    BRepExtrema_DistShapeShape pnt2edge;
    pnt2edge.SetDeflection(0.1);

    TopExp_Explorer aExpFace;
    TopExp_Explorer aExpEdge;
    TopoDS_Vertex V;

    MeshPnt aMeshStruct;
    MeshCore::MeshPointIterator pntIt(m_CadMesh);
    MeshCore::MeshKernel FaceMesh;
    MeshCore::MeshPointArray MeshPnts, MeshPntsCop, MeshPntsCad;
    MeshCore::MeshPointArray MeshPntsCopy;
    MeshCore::MeshFacetArray MeshFacets, MeshFacets2;
    MeshCore::MeshFacetIterator facetIt(m_CadMesh);

    std::pair<Base::Vector3f, MeshPnt> inp;
    std::vector<int> MeanVec;
    std::vector<FacePnt> FacePntVector;
    TopLoc_Location aloc;
    Base::Vector3f mP;
    gp_Pnt e_pnt, m_pnt;
    double dist, distSum, revSum;;
    int n;

    std::map<TopoDS_Edge, std::vector<double>, Edge_Less>::iterator edge_it;
    std::map<Base::Vector3f,MeshPnt,MeshPntLess >::iterator meshIt;
    std::vector<MeshCore::MeshPoint>::iterator pIt;

    m_CurvPos.resize(m_CadMesh.CountPoints());
    m_CurvNeg.resize(m_CadMesh.CountPoints());
    m_MeshStruct.resize(m_CadMesh.CountPoints());
    MeanVec.resize(m_CadMesh.CountPoints(), 0);

    // Fülle Map mit Punkten (Key) und deren, zur Basistriangulierung korrespondierenden, Indizes
    MeshPnts = m_CadMesh.GetPoints();
        MeshPntsCop = MeshPnts;
        MeshFacets2 = m_CadMesh.GetFacets();

    for (unsigned int i=0; i<MeshPnts.size(); ++i)
    {
        aMeshStruct.pnt = MeshPnts[i];        // stores point
        aMeshStruct.index = i;                // stores index
        aMeshStruct.maxCurv = 1e+3;
        aMeshStruct.minCurv = -1e+3;
        inp.first  = aMeshStruct.pnt;
        inp.second = aMeshStruct;
        MeshMap.insert(inp);
    }

    // explores all faces  ------------  Hauptschleife
    for (aExpFace.Init(m_Shape,TopAbs_FACE);aExpFace.More();aExpFace.Next())
    {
        TopoDS_Face aFace = TopoDS::Face(aExpFace.Current());
        TransferFaceTriangulationtoFreeCAD(aFace, FaceMesh);
        MeshPnts.clear();
        MeshPnts = FaceMesh.GetPoints();
        MeshPntsCopy = MeshPnts;
        n = MeshPnts.size();
        TopLoc_Location aLocation;
        // berechne normalen -> m_MeshStruct
        Handle_Poly_Triangulation aTr = BRep_Tool::Triangulation(aFace,aLocation);
        const TColgp_Array1OfPnt& aNodes = aTr->Nodes();
        // create array of node points in absolute coordinate system
        TColgp_Array1OfPnt aPoints(1, aNodes.Length());
        for ( Standard_Integer i = 1; i <= aNodes.Length(); i++)
            aPoints(i) = aNodes(i).Transformed(aLocation);

        const TColgp_Array1OfPnt2d& aUVNodes = aTr->UVNodes();


        BRepAdaptor_Surface aSurface(aFace);
        Base::Vector3f Pnt1, Pnt2;
        gp_Pnt2d par;
        gp_Pnt P;
        gp_Vec D1U, D1V;

        for (int i=1; i<n+1; ++i)
        {
            par = aUVNodes.Value(i);
            aSurface.D1(par.X(),par.Y(),P,D1U,D1V);
            P = aPoints(i);
            Pnt1.x = (float) P.X();
            Pnt1.y = (float) P.Y();
            Pnt1.z = (float) P.Z();
            meshIt = MeshMap.find(Pnt1);
            if (meshIt == MeshMap.end())
            {
                cout << "error";
                return false;
            }

            D1U.Cross(D1V);
            D1U.Normalize();
            if (aFace.Orientation() == TopAbs_FORWARD) D1U.Scale(-1.0);

            Pnt2.Set(float(Pnt1.x+D1U.X()),float(Pnt1.y+D1U.Y()),float(Pnt1.z+D1U.Z()));
            logo.addSingleArrow(Pnt1, Pnt2);

            m_MeshStruct[((*meshIt).second).index].normal.x = (float) D1U.X();
            m_MeshStruct[((*meshIt).second).index].normal.y = (float) D1U.Y();
            m_MeshStruct[((*meshIt).second).index].normal.z = (float) D1U.Z();

        }

        logo.saveToFile("c:/norm.iv");



        // get Inner Mesh Points
        int innerpoints = GetBoundary(FaceMesh, MeshPnts);



        FacePntVector.resize(innerpoints);  // stores inner-points and sourrounding edges-distances

        // explores the edges of the face ----------------------
        for (aExpEdge.Init(aFace,TopAbs_EDGE);aExpEdge.More();aExpEdge.Next())
        {
            TopoDS_Edge aEdge = TopoDS::Edge(aExpEdge.Current());

            edge_it = EdgeMap.find(aEdge);
            if (edge_it == EdgeMap.end())
            {
                cout << "error1";
                return false;
            }

            pnt2edge.LoadS1(aEdge);

            int counter = 0;

            for (pIt = MeshPnts.begin(); pIt != MeshPnts.end(); ++pIt)
            {
                //check if there is no boundary
                if (pIt->_ulProp == 0)
                {
                    m_pnt.SetCoord(pIt->x, pIt->y, pIt->z);
                    VertexBuild.MakeVertex(V,m_pnt,0.001);
                    pnt2edge.LoadS2(V);
                    pnt2edge.Perform();

                    if (pnt2edge.IsDone() == false)
                    {
                        throw Base::Exception("couldn't perform distance calculation pnt2edge \n");
                    }

                    dist = pnt2edge.Value();

                    FacePntVector[counter].pnt = *pIt;
                    FacePntVector[counter].distances.push_back(dist);
                    FacePntVector[counter].MinEdgeOff.push_back((*edge_it).second[0]);
                    FacePntVector[counter].MaxEdgeOff.push_back((*edge_it).second[1]);
                    ++counter;
                }

            }

            // Knoten auf der EDGE die EdgeOffset-Werte zuweisen

            // Edge -> Polygon -> Nodes
            Handle(Poly_PolygonOnTriangulation) polyg;
            Handle_Poly_Triangulation aTrLoc;
            BRep_Tool::PolygonOnTriangulation(aEdge,polyg,aTrLoc,aloc);

            TColStd_Array1OfInteger IndArr(1,polyg->NbNodes());
            IndArr = polyg->Nodes();

            const TColgp_Array1OfPnt& Nodes = aTrLoc->Nodes();
            TColgp_Array1OfPnt TrLocPnts(1, Nodes.Length());
            for ( Standard_Integer i = 1; i <= Nodes.Length(); i++)
                TrLocPnts(i) = Nodes(i).Transformed(aloc);

            for (int k=0; k<polyg->Nodes().Upper(); ++k)
            {
                e_pnt = TrLocPnts(IndArr.Value(k+1));

                mP.x = (float) e_pnt.X();
                mP.y = (float) e_pnt.Y();
                mP.z = (float) e_pnt.Z();

                meshIt = MeshMap.find(mP);
                if (meshIt == MeshMap.end())
                {
                    cout << "error2" << endl;
                    return false;
                }


                ++MeanVec[meshIt->second.index];

                // nehme stets den minimalen wert
                if (((*edge_it).second[0])>((*meshIt).second).minCurv)
                    ((*meshIt).second).minCurv = (*edge_it).second[0];

                if ((*edge_it).second[1]<((*meshIt).second).maxCurv)
                    ((*meshIt).second).maxCurv = (*edge_it).second[1];
            }
        }

        // Knoten INNERHALB eines Faces die Offset-Werte zuweisen
        for (unsigned int k=0; k<FacePntVector.size(); ++k)
        {
            meshIt = MeshMap.find(FacePntVector[k].pnt);

            if (meshIt == MeshMap.end())
            {
                cout << "error3";
                return false;
            }

            distSum = 0.0;
            for (unsigned int l=0; l<FacePntVector[k].distances.size(); ++l)
                distSum +=  FacePntVector[k].distances[l];

            revSum = 0.0;
            for (unsigned int l=0; l<FacePntVector[k].distances.size(); ++l)
            {
                revSum += distSum/FacePntVector[k].distances[l];

                if (((*meshIt).second).minCurv == -1e+3)
                    ((*meshIt).second).minCurv = 0.0;

                if (((*meshIt).second).maxCurv == 1e+3)
                    ((*meshIt).second).maxCurv = 0.0;

                ((*meshIt).second).minCurv += distSum*(FacePntVector[k].MinEdgeOff[l])/FacePntVector[k].distances[l];
                ((*meshIt).second).maxCurv += distSum*(FacePntVector[k].MaxEdgeOff[l])/FacePntVector[k].distances[l];

            }

            ((*meshIt).second).minCurv /= revSum;
            ((*meshIt).second).maxCurv /= revSum;

        }
        FacePntVector.clear();
    }



    Base::Builder3D log,log3d;
    Base::Vector3f pa,pb;
    char text[10];
    int w;

    // übergebe krümmungswerte an m_MeshStruct
    for (meshIt = MeshMap.begin(); meshIt != MeshMap.end(); ++meshIt)
    {
        w = meshIt->second.index;
        m_MeshStruct[w].pnt = (*meshIt).first;

        m_MeshStruct[w].minCurv = (((*meshIt).second).minCurv);
        m_MeshStruct[w].maxCurv = (((*meshIt).second).maxCurv);

        // Ausgabe
        if (m_MeshStruct[w].maxCurv<10000000000)
        {
            snprintf(text,10,"%f",m_MeshStruct[w].minCurv);
            //snprintf(text,10,"%i",w);
            log.addText(m_MeshStruct[w].pnt.x, m_MeshStruct[w].pnt.y, m_MeshStruct[w].pnt.z,text);
        }

        if (MeanVec[w] == 1)
            log3d.addSinglePoint(m_MeshStruct[w].pnt,4,0,0,0); // innere punkte - > schwarz
        else if (MeanVec[w] == 2)
            log3d.addSinglePoint(m_MeshStruct[w].pnt,4,1,1,1); // edge-punkte   - > weiß
        else
            log3d.addSinglePoint(m_MeshStruct[w].pnt,4,1,0,0); // eck-punkte    - > rot

        //m_MeshStruct[w].minCurv = (((*meshIt).second).minCurv)/MeanVec[w];
        //m_MeshStruct[w].maxCurv = (((*meshIt).second).maxCurv)/MeanVec[w];

        //if(m_MeshStruct[w].maxCurv < 100 || fabs(m_MeshStruct[w].minCurv) < 100)
        //{
        // if(m_MeshStruct[w].maxCurv > fabs(m_MeshStruct[w].minCurv))
        // {
        //  snprintf(text,10,"%f",m_MeshStruct[w].minCurv);
        //  //m_MeshStruct[w].normal.Scale(m_MeshStruct[w].minCurv,m_MeshStruct[w].minCurv,m_MeshStruct[w].minCurv);
        //  logo.addSingleArrow(m_MeshStruct[w].pnt, m_MeshStruct[w].pnt+m_MeshStruct[w].normal);
        // }
        // else
        // {
        //  //snprintf(text,4,"%f",m_MeshStruct[w].maxCurv);
        //  m_MeshStruct[w].normal.Scale(m_MeshStruct[w].maxCurv,m_MeshStruct[w].maxCurv,m_MeshStruct[w].maxCurv);
        //  logo.addSingleArrow(m_MeshStruct[w].pnt, m_MeshStruct[w].pnt+m_MeshStruct[w].normal);
        // }
        //}

        //if(m_MeshStruct[w].maxCurv < 100)
        //{
        // m_MeshStruct[w].normal.Scale(0.5*m_MeshStruct[w].maxCurv,0.5*m_MeshStruct[w].maxCurv,0.5*m_MeshStruct[w].maxCurv);
        // logo.addSingleArrow(m_MeshStruct[w].pnt, m_MeshStruct[w].pnt+m_MeshStruct[w].normal);
        //}*/

        //snprintf(text,10,"%f",m_MeshStruct[w].maxCurv);
        //snprintf(text,10,"%i",w);


        /*if(m_MeshStruct[w].maxCurv > fabs(m_MeshStruct[w].minCurv))
        {
         snprintf(text,10,"%f",m_MeshStruct[w].minCurv);
         
        }
        else
        {
         snprintf(text,10,"%f",m_MeshStruct[w].maxCurv);
        }*/



    }

    log.saveToFile("c:/printCurv.iv");
    log3d.saveToFile("c:/triPnts.iv");
    logo.saveToFile("c:/NormalsCurv.iv");

        MeshPntsCad = m_CadMesh.GetPoints();
        
        
         for (unsigned int j=0; j<m_FixFaces.size(); ++j)
         {
        TopoDS_Face aFace = TopoDS::Face(m_FixFaces[j]);
        TransferFaceTriangulationtoFreeCAD(aFace, FaceMesh);
                MeshPnts = FaceMesh.GetPoints();
                MeshFacets2 = m_CadMesh.GetFacets();
        MeshPntsCopy = MeshPnts;
        n = MeshPnts.size();
        TopLoc_Location aLocation;
        // berechne normalen -> m_MeshStruct
        Handle_Poly_Triangulation aTr = BRep_Tool::Triangulation(aFace,aLocation);
        const TColgp_Array1OfPnt& aNodes = aTr->Nodes();
        // create array of node points in absolute coordinate system
        TColgp_Array1OfPnt aPoints(1, aNodes.Length());
        for ( Standard_Integer i = 1; i <= aNodes.Length(); i++)
            aPoints(i) = aNodes(i).Transformed(aLocation);

        const TColgp_Array1OfPnt2d& aUVNodes = aTr->UVNodes();

        BRepAdaptor_Surface aSurface(aFace);
        Base::Vector3f Pnt1, Pnt2;
        gp_Pnt2d par;
        gp_Pnt P;
        gp_Vec D1U, D1V;
                std::vector<bool> trans_check(m_MeshStruct.size(), false);

        for (int i=1; i<n+1; ++i)
        {
            par = aUVNodes.Value(i);
            aSurface.D1(par.X(),par.Y(),P,D1U,D1V);
            P = aPoints(i);
            Pnt1.x = (float) P.X();
            Pnt1.y = (float) P.Y();
            Pnt1.z = (float) P.Z();
            meshIt = MeshMap.find(Pnt1);
            if (meshIt == MeshMap.end())
            {
                cout << "error";
                return false;
            }

            D1U.Cross(D1V);
            D1U.Normalize();
            if (aFace.Orientation() == TopAbs_FORWARD) D1U.Scale(-1.0);

                        int g = ((*meshIt).second).index;

                        D1U.Scale(0.1);
                        if(trans_check[g] == false)
                        {
                                m_dist_vec[g].x = - D1U.X();
                                m_dist_vec[g].y = - D1U.Y();
                                m_dist_vec[g].z = - D1U.Z();

                                MeshPntsCad[g].Set(MeshPntsCad[g].x - D1U.X(),
                                                                   MeshPntsCad[g].y - D1U.Y(),
                                                                   MeshPntsCad[g].z - D1U.Z());
                                
                                (m_MeshStruct[g].pnt).Set(MeshPntsCad[g].x - D1U.X(),
                                                                                  MeshPntsCad[g].y - D1U.Y(),
                                                                          MeshPntsCad[g].z - D1U.Z());

                                trans_check[g] = true;
                        }
                }
         }
         m_CadMesh.Assign(MeshPntsCad, MeshFacets2);

     return true;
}


bool SpringbackCorrection::Init()
{
        m_EdgeStruct.clear();
    EdgeMap.clear();
    MeshMap.clear();
    //Remove any existing Triangulation on the Shape
    BRepTools::Clean(m_Shape);
    best_fit::Tesselate_Shape(m_Shape,m_CadMesh,(float) 0.1); // Basistriangulierung des ganzen Shapes

    MeshCore::MeshTopoAlgorithm algo(m_CadMesh);
    algo.HarmonizeNormals();

    MeshCore::MeshGeomFacet facet = m_CadMesh.GetFacet(0);
    if (facet.GetNormal().z < 0.0)
        algo.FlipNormals();

    int n = m_CadMesh.CountFacets();
    Base::Vector3f normal;

    MeshCore::MeshPointArray points = m_CadMesh.GetPoints();
    MeshCore::MeshFacetArray facets = m_CadMesh.GetFacets();

    for (int i=0; i<n; ++i)
    {
        MeshCore::MeshGeomFacet face = m_CadMesh.GetFacet(i);
        normal = face.GetNormal();

        if (normal.z < 0.0)
        {
            facets[i].FlipNormal();
        }
    }

    m_CadMesh.Assign(points, facets);

    double max = cMin;
    TopTools_IndexedDataMapOfShapeListOfShape anIndex;
    anIndex.Clear();
    TopExp aMap;
    aMap.MapShapesAndAncestors(m_Shape,TopAbs_EDGE,TopAbs_FACE,anIndex);
    TopExp_Explorer anExplorer, anExplorer2;
    TopoDS_Face aFace;
    EdgeStruct tempEdgeStruct;
    std::vector<EdgeStruct> EdgeVec;
    MeshCore::MeshKernel FaceMesh;
    std::vector<double> curvature(2);
    std::pair<TopoDS_Edge, std::vector<double> > aPair;

    float maxOffset, minOffset;

    for (anExplorer.Init(m_Shape,TopAbs_EDGE);anExplorer.More();anExplorer.Next())
    {
        tempEdgeStruct.aFace.clear();
        tempEdgeStruct.MaxOffset = 1e+10;
        tempEdgeStruct.MinOffset = -1e+10;
        tempEdgeStruct.anEdge = TopoDS::Edge(anExplorer.Current()); // store current edge

        const TopTools_ListOfShape& aFaceList = anIndex.FindFromKey(tempEdgeStruct.anEdge);
        TopTools_ListIteratorOfListOfShape aListIterator;

        for (aListIterator.Initialize(aFaceList);aListIterator.More(); aListIterator.Next())
        {
            // store corresponding faces
            aFace = TopoDS::Face(aListIterator.Value());
            tempEdgeStruct.aFace.push_back(aFace);
            TransferFaceTriangulationtoFreeCAD(aFace, FaceMesh);
            curvature = MeshCurvature(aFace, FaceMesh);

            if (aFace.Orientation() == TopAbs_REVERSED)
            {
                maxOffset = (float) -curvature[1];
                minOffset = (float) -curvature[0];
            }
            else
            {
                maxOffset = (float) curvature[0];
                minOffset = (float) curvature[1];
            }

            if (maxOffset < tempEdgeStruct.MaxOffset)
                tempEdgeStruct.MaxOffset = maxOffset;

            if (minOffset > tempEdgeStruct.MinOffset)
                tempEdgeStruct.MinOffset = minOffset;
        }


        if (tempEdgeStruct.MaxOffset > max)
            curvature[1] = tempEdgeStruct.MaxOffset;
        else
            curvature[1] = cMin;

        if (-tempEdgeStruct.MinOffset > max)
            curvature[0] = tempEdgeStruct.MinOffset;
        else
            curvature[0] = -cMin;

        aPair.first = tempEdgeStruct.anEdge;
        aPair.second = curvature;
        EdgeMap.insert(aPair);

    }
    return true;
}

bool SpringbackCorrection::Init_Setting(struct CuttingToolsSettings& set)
{
    m_set = set;
    return true;
}

bool SpringbackCorrection::SetFixEdges()
{
    TopExp_Explorer anExplorer;
    TopoDS_Edge anEdge;
    std::map<TopoDS_Edge, std::vector<double>, Edge_Less>::iterator edge_it;
    std::vector<double> pair(2,0.0);

    for (unsigned int i=0; i<m_FixFaces.size(); ++i)
    {
        for (anExplorer.Init(m_FixFaces[i],TopAbs_EDGE); anExplorer.More(); anExplorer.Next())
        {
            anEdge  = TopoDS::Edge(anExplorer.Current());
            edge_it = EdgeMap.find(anEdge);
            edge_it->second = pair;
        }
    }

    return true;
}

bool SpringbackCorrection::TransferFaceTriangulationtoFreeCAD(const TopoDS_Face& aFace, MeshCore::MeshKernel& FaceMesh)
{
    FaceMesh.Clear();
    MeshCore::MeshBuilder builder(FaceMesh);
    builder.Initialize(1000);
    Base::Vector3f Points[3];

    TopLoc_Location aLocation;
    // takes the triangulation of the face aFace:
    Handle_Poly_Triangulation aTr = BRep_Tool::Triangulation(aFace,aLocation);
    if (!aTr.IsNull()) // if this triangulation is not NULL
    {
        // takes the array of nodes for this triangulation:
        const TColgp_Array1OfPnt& aNodes = aTr->Nodes();
        // takes the array of triangles for this triangulation:
        const Poly_Array1OfTriangle& triangles = aTr->Triangles();
        // create array of node points in absolute coordinate system
        TColgp_Array1OfPnt aPoints(1, aNodes.Length());
        for ( Standard_Integer i = 1; i <= aNodes.Length(); i++)
            aPoints(i) = aNodes(i).Transformed(aLocation);
        // Takes the node points of each triangle of this triangulation.
        // takes a number of triangles:
        Standard_Integer nnn = aTr->NbTriangles();
        Standard_Integer nt,n1,n2,n3;
        for ( nt = 1 ; nt < nnn+1 ; nt++)
        {
            // takes the node indices of each triangle in n1,n2,n3:
            triangles(nt).Get(n1,n2,n3);
            // takes the node points:
            gp_Pnt aPnt1 = aPoints(n1);
            Points[0].Set(float(aPnt1.X()),float(aPnt1.Y()),float(aPnt1.Z()));
            gp_Pnt aPnt2 = aPoints(n2);
            Points[1].Set((float) aPnt2.X(),(float) aPnt2.Y(),(float) aPnt2.Z());
            gp_Pnt aPnt3 = aPoints(n3);
            Points[2].Set((float) aPnt3.X(),(float) aPnt3.Y(),(float) aPnt3.Z());
            // give the occ faces to the internal mesh structure of freecad
            MeshCore::MeshGeomFacet Face(Points[0],Points[1],Points[2]);
            builder.AddFacet(Face);
        }

    }
    // if the triangulation of only one face is not possible to get
    else
    {
        throw Base::Exception("Empty face triangulation\n");
    }

    // finish FreeCAD Mesh Builder and exit with new mesh
    builder.Finish();
    return true;

}

/*
MeshCore::MeshKernel SpringbackCorrection::BuildMesh(Handle_Poly_Triangulation aTri, std::vector<Base::Vector3f> TrPoints)
{
    MeshCore::MeshKernel FaceMesh;
    MeshCore::MeshBuilder builder(FaceMesh);
    builder.Initialize(1000);
    Base::Vector3f Points[3];

    const Poly_Array1OfTriangle& triangles = aTri->Triangles();
    Standard_Integer nnn = aTri->NbTriangles();
    Standard_Integer nt,n1,n2,n3;

    for ( nt = 1 ; nt < nnn+1 ; nt++)
    {
        // takes the node indices of each triangle in n1,n2,n3:
        triangles(nt).Get(n1,n2,n3);
        // takes the node points:
        Points[0] = TrPoints[n1];
        Points[1] = TrPoints[n2];
        Points[2] = TrPoints[n3];
        // give the occ faces to the internal mesh structure of freecad
        MeshCore::MeshGeomFacet Face(Points[0],Points[1],Points[2]);
        builder.AddFacet(Face);
    }

    // finish FreeCAD Mesh Builder and exit with new mesh
    builder.Finish();

    return FaceMesh;
}
*/

int SpringbackCorrection::GetBoundary(const MeshCore::MeshKernel &mesh, MeshCore::MeshPointArray &meshPnts)
{
    // get mesh-boundary
    // zweck: nur für die inneren punkte wird der curvature wert über die edge abstände berechnet
    std::list< std::vector <unsigned long> >  BoundariesIndex;
    std::list< std::vector <unsigned long> >::iterator bInd;


    MeshCore::MeshAlgorithm algo(mesh);
    algo.GetMeshBorders(BoundariesIndex);
    //Set all MeshPoints to Property 0
    for (unsigned int i=0;i<meshPnts.size();++i)
        meshPnts[i].SetProperty(0);

    //Set Boundary Points to 1
    int inner_points = 0;
    for (bInd = BoundariesIndex.begin(); bInd != BoundariesIndex.end(); ++bInd)
    {
        for (unsigned int i=0;i< bInd->size();++i)
        {
            meshPnts[(*bInd).at(i)].SetProperty(1);
        }

    }
    for (unsigned int i=0;i<meshPnts.size();++i)
    {
        if (meshPnts[i]._ulProp == 0)
            inner_points++;
    }

    return (inner_points);
}
/*
bool SpringbackCorrection::SmoothCurvature()
{
    MeshCore::MeshPointIterator v_it(m_Mesh);
    MeshCore::MeshRefPointToPoints vv_it(m_Mesh);
    std::set<MeshCore::MeshPointArray::_TConstIterator> PntNei;
    std::set<MeshCore::MeshPointArray::_TConstIterator>::iterator pnt_it;

    int n = m_MeshStruct.size();

    std::vector<double> SCMax(n);
    std::vector<double> SCMin(n);
    double maxCurv = 0.0, minCurv = 0.0;

    for (v_it.Begin(); v_it.More(); v_it.Next())
    {
        PntNei = vv_it[v_it.Position()];
        maxCurv += m_MeshStruct[v_it.Position()].maxCurv;
        minCurv += m_MeshStruct[v_it.Position()].minCurv;

        for (pnt_it = PntNei.begin(); pnt_it !=PntNei.end(); ++pnt_it)
        {
            maxCurv += m_MeshStruct[(*pnt_it)[0]._ulProp].maxCurv;
            minCurv += m_MeshStruct[(*pnt_it)[0]._ulProp].minCurv;
        }

        SCMax[v_it.Position()] = maxCurv/(PntNei.size()+1);
        SCMin[v_it.Position()] = minCurv/(PntNei.size()+1);

        maxCurv = 0.0;
        minCurv = 0.0;
    }

    for (int i=0; i<n; ++i)
    {
        m_MeshStruct[i].maxCurv = SCMax[i];
        m_MeshStruct[i].minCurv = SCMin[i];
    }

    return true;
}
*/

bool SpringbackCorrection::SmoothMesh(MeshCore::MeshKernel &Mesh, double d_max)
{
    Base::Builder3D log;
    MeshCore::MeshKernel localMesh, Meshtmp;
    MeshCore::MeshPoint mpnt, spnt;
    MeshCore::MeshPointArray locPointArray, PointArray = Mesh.GetPoints();
    MeshCore::MeshFacetArray locFacetArray, FacetArray = Mesh.GetFacets();

    MeshCore::MeshPointIterator v_it(Mesh);
    MeshCore::MeshRefPointToPoints vv_it(Mesh);
    std::set<unsigned long>::const_iterator pnt_it;
    MeshCore::MeshPointArray::_TConstIterator v_beg = Mesh.GetPoints().begin();

    Base::Vector3f N, L, coor;
    int n = m_MeshStruct.size();
    int g;

    for (v_it.Begin(); v_it.More(); v_it.Next())
    {

        g = v_it.Position();
        spnt.Set(0.0, 0.0, 0.0);
        locPointArray.push_back(*v_it);
        spnt += *v_it;
        const std::set<unsigned long>& PntNei = vv_it[(*v_it)._ulProp];

        if (PntNei.size() < 3)
            continue;

        for (pnt_it = PntNei.begin(); pnt_it !=PntNei.end(); ++pnt_it)
        {
            locPointArray.push_back(v_beg[*pnt_it]);
            spnt += v_beg[*pnt_it];
        }

        spnt.Scale((float) (1.0/(double(PntNei.size()) + 1.0)),(float) (1.0/(double(PntNei.size()) + 1.0)),(float) (1.0/(double(PntNei.size()) + 1.0)));

        localMesh.Assign(locPointArray, locFacetArray);


        // need at least one facet to perform the PCA
        MeshCore::MeshGeomFacet face;
        face._aclPoints[0] = locPointArray[0];
        face._aclPoints[1] = locPointArray[1];
        face._aclPoints[2] = locPointArray[2];
        localMesh.AddFacet(face);

        MeshCore::MeshEigensystem pca(localMesh);
        pca.Evaluate();
        Base::Matrix4D T1 = pca.Transform();

        N.x = (float) T1[0][0];
        N.y = (float) T1[0][1];
        N.z = (float) T1[0][2];
        L.Set(v_it->x - spnt.x, v_it->y - spnt.y, v_it->z - spnt.z);
        N.Normalize();


        if (N*L < 0.0)
            N.Scale(-1.0, -1.0, -1.0);

        if ((*v_it)._ulProp == 2286)
        {
            log.addSinglePoint(spnt,4,1,0,0);
            log.addSinglePoint(*v_it,4,0,0,0);
            for (unsigned int i=0; i<locPointArray.size(); ++i)
                log.addSinglePoint(locPointArray[i],4);

            for (int i=0; i<3; ++i)
            {
                coor.x = (float) T1[i][0];
                coor.y = (float) T1[i][1];
                coor.z = (float) T1[i][2];

                coor.Normalize();

                log.addSingleArrow(*v_it, *v_it + coor);
            }

            /*log.addSingleArrow(*v_it, *v_it - N);
            snprintf(text,10,"%i",(*v_it)._ulProp);
            log.addText(v_it->x, v_it->y, v_it->z,text);*/
        }

        double d = d_max;
        if (fabs(N*L)< d_max) d = fabs(N*L);
        N.Scale((float) d,(float) d,(float) d);

        PointArray[v_it.Position()].Set(v_it->x - N.x, v_it->y - N.y, v_it->z - N.z);
        locPointArray.clear();
        localMesh = Meshtmp;
    }

    Mesh.Assign(PointArray, FacetArray);
    log.saveToFile("c:/checkNormals.iv");

    return true;
}

bool SpringbackCorrection::SmoothMesh(MeshCore::MeshKernel &Mesh, std::vector<int> ind, double d_max)
{
    Base::Builder3D log;
    MeshCore::MeshKernel localMesh, Meshtmp;
    MeshCore::MeshPoint mpnt, spnt;
    MeshCore::MeshPointArray locPointArray, PointArray = Mesh.GetPoints();
    MeshCore::MeshFacetArray locFacetArray, FacetArray = Mesh.GetFacets();

    MeshCore::MeshPointIterator v_it(Mesh);
    MeshCore::MeshRefPointToPoints vv_it(Mesh);
    std::set<unsigned long>::const_iterator pnt_it;
    MeshCore::MeshPointArray::_TConstIterator v_beg = Mesh.GetPoints().begin();

    Base::Vector3f N, L, coor;

    int n = ind.size();

    for (int i=0; i<n; ++i)
    {
        v_it.Set(ind[i]);
        spnt.Set(0.0, 0.0, 0.0);
        locPointArray.push_back(*v_it);
        spnt += *v_it;
        const std::set<unsigned long>& PntNei = vv_it[(*v_it)._ulProp];

        if (PntNei.size() < 3)
            continue;

        for (pnt_it = PntNei.begin(); pnt_it !=PntNei.end(); ++pnt_it)
        {
            locPointArray.push_back(v_beg[*pnt_it]);
            spnt += v_beg[*pnt_it];
        }

        spnt.Scale((float) (1.0/(double(PntNei.size()) + 1.0)),(float) (1.0/(double(PntNei.size()) + 1.0)),(float) (1.0/(double(PntNei.size()) + 1.0)));

        localMesh.Assign(locPointArray, locFacetArray);


        // need at least one facet to perform the PCA
        MeshCore::MeshGeomFacet face;
        face._aclPoints[0] = locPointArray[0];
        face._aclPoints[1] = locPointArray[1];
        face._aclPoints[2] = locPointArray[2];
        localMesh.AddFacet(face);

        MeshCore::MeshEigensystem pca(localMesh);
        pca.Evaluate();
        Base::Matrix4D T1 = pca.Transform();

        N.x = (float) T1[0][0];
        N.y = (float) T1[0][1];
        N.z = (float) T1[0][2];
        L.Set(v_it->x - spnt.x, v_it->y - spnt.y, v_it->z - spnt.z);
        N.Normalize();


        if (N*L < 0.0)
            N.Scale(-1.0, -1.0, -1.0);

        if ((*v_it)._ulProp == 2286)
        {
            log.addSinglePoint(spnt,4,1,0,0);
            log.addSinglePoint(*v_it,4,0,0,0);
            for (unsigned int i=0; i<locPointArray.size(); ++i)
                log.addSinglePoint(locPointArray[i],4);

            for (int i=0; i<3; ++i)
            {
                coor.x = (float) T1[i][0];
                coor.y = (float) T1[i][1];
                coor.z = (float) T1[i][2];

                coor.Normalize();

                log.addSingleArrow(*v_it, *v_it + coor);
            }

            /*log.addSingleArrow(*v_it, *v_it - N);
            snprintf(text,10,"%i",(*v_it)._ulProp);
            log.addText(v_it->x, v_it->y, v_it->z,text);*/
        }

        double d = d_max;
        if (fabs(N*L)< d_max) d = fabs(N*L);
        N.Scale((float) d,(float) d,(float) d);

        PointArray[v_it.Position()].Set(v_it->x - N.x, v_it->y - N.y, v_it->z - N.z);
        locPointArray.clear();
        localMesh = Meshtmp;
    }

    Mesh.Assign(PointArray, FacetArray);
    log.saveToFile("c:/checkNormals.iv");

    return true;
}

/*
bool SpringbackCorrection::GetFaceAng(MeshCore::MeshKernel &mesh, int deg_tol)
{
    Base::Vector3f normal, base;
    double deg;
    int n = mesh.CountFacets();
    bool b = true;

    for (int i=0; i<n; ++i)
    {
        MeshCore::MeshGeomFacet face = mesh.GetFacet(i);
        face.CalcNormal();
        normal = face.GetNormal();
        normal.Normalize();
        base = normal;
        base.z = 0.0;
        base.Normalize();

        if (normal.z < 0.0)
            deg = 90.0 + acos(normal*base)*180.0/PI;
        else
            deg = 90.0 - acos(normal*base)*180.0/PI;

        if (deg > deg_tol)
            b = false;

        m_FaceAng.push_back(deg);
    }

    return b;
}
*/

std::vector<int> SpringbackCorrection::InitFaceCheck(MeshCore::MeshKernel &mesh, int degree)
{
    std::vector<int> faceInd;
    double tol = degree;
    double alpha;
    MeshCore::MeshFacetIterator fIt(mesh);
    MeshCore::MeshPointIterator mIt(mesh);
    MeshCore::MeshRefPointToFacets p2fIt(mesh);
    Base::Vector3f normal, base, gpnt;
    Base::Builder3D log;

    MeshCore::MeshPointArray mPnts   = mesh.GetPoints();
    MeshCore::MeshFacetArray mFacets = mesh.GetFacets();
    MeshCore::MeshFacetArray::_TConstIterator f_beg = mesh.GetFacets().begin();

    int n = mesh.CountFacets();


    for (int i=0; i<n; ++i)
    {
        MeshCore::MeshGeomFacet face = mesh.GetFacet(i);
        //face.CalcNormal();
        normal = face.GetNormal();
        normal.Normalize();
        base = normal;
        base.z = 0.0;
        base.Normalize();

        gpnt = face.GetGravityPoint();


        if (normal.z < 0.0)
        {
            alpha = 90.0 + acos(normal*base)*180.0/PI;
            //log.addSingleArrow(gpnt,gpnt+normal,6,1,0,0);
        }
        else
        {
            alpha = 90.0 - acos(normal*base)*180.0/PI;
            //log.addSingleArrow(gpnt,gpnt+normal);
        }

        if (alpha > tol)
        {
            // markiere kritisches face
            mFacets[i].SetProperty(0);

            for (int j=0; j<3; ++j)
            {
                const std::set<unsigned long>& faceSet = p2fIt[mFacets[i]._aulPoints[j]];

                for (std::set<unsigned long>::const_iterator it = faceSet.begin(); it != faceSet.end(); ++it)
                {
                    f_beg[*it].SetProperty(5);
                }
            }

            log.addSingleArrow(gpnt,gpnt+normal,4,1,0,0);
        }
        else
        {
            mFacets[i].SetProperty(1);
            faceInd.push_back(i);
        }
    }

    log.saveToFile("c:/normalschecker.iv");

    MeshCore::MeshFacetArray mFacets2 = mesh.GetFacets();

    for (unsigned int i=0; i<mFacets2.size(); ++i)
    {
        if (mFacets2[i]._ulProp == 5)
        {
            mFacets[i].SetProperty(0);
        }
    }

    mesh.Assign(mPnts, mFacets);

    log.saveToFile("c:/angles.iv");
    return faceInd;
}

std::vector<int> SpringbackCorrection::FaceCheck(MeshCore::MeshKernel &mesh, int degree)
{
    std::vector<int> faceInd;
    double tol = degree;
    double alpha;
    MeshCore::MeshFacetIterator fIt(mesh);
    MeshCore::MeshPointIterator mIt(mesh);
    Base::Vector3f normal, base, gpnt;
    Base::Builder3D log;


    MeshCore::MeshPointArray mPnts   = mesh.GetPoints();
    MeshCore::MeshFacetArray mFacets = mesh.GetFacets();

    int n = mesh.CountFacets();

    for (int i=0; i<n; ++i)
    {
        MeshCore::MeshGeomFacet face = mesh.GetFacet(i);
        face.CalcNormal();
        normal = face.GetNormal();
        normal.Normalize();
        base = normal;
        base.z = 0.0;
        base.Normalize();

        gpnt = face.GetGravityPoint();

        if (normal.z < 0.0)
        {
            alpha = 90.0 + acos(normal*base)*180.0/PI;
            log.addSingleArrow(gpnt,gpnt+normal,2,1,0,0);
        }
        else
            alpha = 90.0 - acos(normal*base)*180.0/PI;

        if (alpha > tol)
        {
            faceInd.push_back(i);
                        cout << "i: " << i << ", angle: " << alpha << endl;
        }
    }

    log.saveToFile("c:/angles.iv");
    return faceInd;
}

/*
bool SpringbackCorrection::CorrectScale(double zLim)
{
    // VERSCHIEBUNGSVEKTOREN ...
    // berechnung der skalierungsfaktoren bzgl. der Cad-Normalen
    int n = m_CadMesh.CountPoints();
    double zLev;
    for (int i=0; i<n; ++i)
    {
        m_MeshStruct[i].normal.Normalize();
        zLev = m_MeshStruct[i].pnt.z + m_Offset[i]*m_MeshStruct[i].normal.z;

        if (zLev>zLim)
            m_Offset[i] = (zLim - m_MeshStruct[i].pnt.z)/m_MeshStruct[i].normal.z;
    }

    return true;
}
*/

bool SpringbackCorrection::Perform(int deg_Tol, bool out)
{
    unsigned int NumRef;
        std::vector<double> dist_vec;

    //GetFaceAng(m_CadMesh, deg_Tol+1);
    //cout << "Init" << endl;
    //Init();      // tesseliere shape -> Basistriangulierung CAD-Mesh
    // EdgeMap wird hier gefüllt

        Base::Vector3f nullvec(0.0,0.0,0.0);
        m_dist_vec.resize(m_CadMesh.CountPoints(), nullvec); // fülle mit Nullvektoren

    cout << "SetFixEdges" << endl;
    
        SetFixEdges();

        const MeshCore::MeshKernel RefMesh23 = m_CadMesh;

    cout << "CalcCurv" << endl;

    if (!CalcCurv())    
           return false;
        
        // berechne Krümmungswerte über Edgekrümmungen
    // MeshMap wird hier gefüllt

    const MeshCore::MeshKernel RefMesh = m_CadMesh;  // übegebe CAD-Mesh vor der Verformung

    /*Base::Builder3D log;
    Base::Vector3f gpnt, normal;
    int numb = m_CadMesh.CountFacets();

    for (int i=0; i<numb; ++i)
       {
           MeshCore::MeshGeomFacet face = m_CadMesh.GetFacet(i);
           normal = face.GetNormal();


     if(normal.z < 0.0)
      log.saveToFile("c:/didid.iv");

     gpnt = face.GetGravityPoint();
     normal.Scale(10,10,10);
     log.addSingleArrow(gpnt,gpnt+normal);
    }

    log.saveToFile("c:/Facenormals.iv");*/

    // speichere normalen seperat
    m_normals.clear();
    m_normals.resize(m_MeshStruct.size());

    //Base::Builder3D logo;

    for (unsigned int i=0; i<m_normals.size(); ++i)
    {
        m_normals[i] = m_MeshStruct[i].normal;  // Normale am Punkt i der Basistriangulierung
        //logo.addSingleArrow(m_MeshStruct[i].pnt, m_MeshStruct[i].pnt + m_normals[i]);
    }

    //logo.saveToFile("c:/normals.iv");

    // übergebe Normalen und CAD-Mesh für Fehlerberechnng
    best_fit befi;
    befi.m_normals = m_normals;
    befi.m_CadMesh = m_CadMesh;

        m_error.resize(m_CadMesh.CountPoints(), 0.0);

        for(unsigned int i=0; i<m_MeshVec.size(); ++i)
        {
                befi.CompTotalError(m_MeshVec[i]);  // Fehlerberechnung
        }

    // VERSCHIEBUNGSVEKTOREN ...
    // berechnung der skalierungsfaktoren bzgl. der Cad-Normalen
    int n = m_CadMesh.CountPoints();
    for (int i=0; i<n; ++i)
    {
        m_MeshStruct[i].error = befi.m_error[i];

        if (m_MeshStruct[i].error < 0)
        {
            if ((m_MeshStruct[i].maxCurv - m_set.master_radius) < 0.0)
                m_Offset.push_back(0.0);
            else if (-m_MeshStruct[i].error > (m_MeshStruct[i].maxCurv - m_set.master_radius))
                m_Offset.push_back(m_MeshStruct[i].maxCurv - m_set.master_radius);
            else
                m_Offset.push_back(-m_MeshStruct[i].error);
        }
        else
        {
            if ((m_MeshStruct[i].minCurv + m_set.slave_radius) > 0.0)
                m_Offset.push_back(0.0);
            else if (m_MeshStruct[i].error > -(m_MeshStruct[i].minCurv + m_set.slave_radius))
                m_Offset.push_back(m_MeshStruct[i].minCurv + m_set.slave_radius);
            else
                m_Offset.push_back(-m_MeshStruct[i].error);
        }
    }

    //Base::BoundBox3f bbox = m_CadMesh.GetBoundBox();  // hole bounding-box
    //CorrectScale(bbox.MaxZ);                          // korrigiere verschiebungsvektoren...

        //m_set.correction_factor = -0.1;
    // skalierung der Cad-Normalen
    for (int i=0; i<n; ++i)
    {
        m_normals[i].Scale((float) (m_set.correction_factor*m_Offset[i]),
                           (float) (m_set.correction_factor*m_Offset[i]),
                           (float) (m_set.correction_factor*m_Offset[i]));
    }

        for(int i=0; i<m_normals.size(); ++i)
        {
                 m_dist_vec[i] = m_normals[i];
        }

        if(out==true) //hier werden die Outputvektoren für Catia geschrieben
        {
                MeshCore::MeshPointArray mpts = RefMesh23.GetPoints();

                Base::Builder3D loo;
                ofstream anOutputFile;
                anOutputFile.open("c:/catia_offset.txt");
                anOutputFile.precision(7);

                for(int i=0; i< mpts.size(); ++i)
                {
                        loo.addSingleArrow(mpts[i], mpts[i] + m_dist_vec[i]);
                        anOutputFile <<  mpts[i].x << " " << mpts[i].y << " " << mpts[i].z << " " << m_dist_vec[i].x << " " << m_dist_vec[i].y << " " << m_dist_vec[i].z << endl;
                }

                loo.saveToFile("c:/prpopo.iv");
            anOutputFile.close();
                return true; //mehr braucne wir auch nicht....
        }



    std::vector<int> tmpVec;
    tmpVec = FaceCheck(m_CadMesh, deg_Tol+1);
        
    NumRef = (int) tmpVec.size();

    MeshCore::MeshPointArray pntAr = m_CadMesh.GetPoints();
    MeshCore::MeshFacetArray facAr = m_CadMesh.GetFacets();

    // Spiegelung
    Base::Builder3D log3d;
    for (int i=0; i<n; ++i)
    {
        log3d.addSingleArrow(pntAr[i], pntAr[i] + m_normals[i]);
        pntAr[i].Set(pntAr[i].x + m_normals[i].x, pntAr[i].y + m_normals[i].y, pntAr[i].z + m_normals[i].z);
    }

    log3d.saveToFile("c:/project2result.iv");

    m_CadMesh.Assign(pntAr, facAr);

    //std::vector<int> tmpVec;
    std::vector<Base::Vector3f> normalsRef = m_normals;
    int cc = 0;



    // ********************************** GLOBALE KORREKTUR ****************************************
    int cm = 49;
    InitFaceCheck(m_CadMesh, deg_Tol+1);  // kritische dreiecke -> _ulProp = 0
    MeshCore::MeshFacetArray facAr2 = m_CadMesh.GetFacets();

    while (true)
    {
        tmpVec = FaceCheck(m_CadMesh, deg_Tol+1);
        cout << "remaining: " << tmpVec.size() << endl;
        if (tmpVec.size() == NumRef || cm == 0)
        {
            cout << cm << "left" << endl;
            break;
        }

        --cm;
        MeshCore::MeshPointArray pntAr2 = RefMesh.GetPoints();

        for (unsigned int i=0; i<pntAr2.size(); ++i)
        {
            m_normals[i].Normalize();
            m_normals[i].Scale((cm*normalsRef[i].Length())/50,(cm*normalsRef[i].Length())/50, (cm*normalsRef[i].Length())/50);
            pntAr2[i].Set(pntAr2[i].x + m_normals[i].x, pntAr2[i].y + m_normals[i].y, pntAr2[i].z + m_normals[i].z);
        }

        m_CadMesh.Assign(pntAr2, facAr2);
    }

        for(int i=0;  i< m_normals.size(); ++i)
        {
                 m_dist_vec[i] += m_normals[i];
        }


        
    // *********************** REGION-GROWING **********************

    MeshCore::MeshFacetArray mFacets = m_CadMesh.GetFacets();
    MeshCore::MeshPointArray mPoints = m_CadMesh.GetPoints();
    int num = mFacets.size();

    for (int i=0; i<num; ++i) mFacets[i].ResetFlag(MeshCore::MeshFacet::VISIT);   // lösche alle VISIT-Flags

    m_RingCurrent = 0;
    std::vector<unsigned long> aRegion;
    std::vector< std::pair<unsigned long, double> > skals;
    std::vector< std::vector< std::pair<unsigned long, double> > > RegionSkals;

        MeshCore::MeshBuilder builder(m_Mesh_vis);
        builder.Initialize(10000);
        MeshCore::MeshBuilder builder2(m_Mesh_vis2);
        builder2.Initialize(10000);

    for (int i=0; i<num; ++i)
    {

        cout << i << " von " << num << endl;
        for (int m=0; m<num; ++m)
        {
            if (mFacets[m]._ulProp == 0)
                        {
                mFacets[m].SetFlag(MeshCore::MeshFacet::VISIT);  // kritische faces kriegen ein visit-flag gesetzt
                        }
        }

                if (mFacets[i]._ulProp == 0)
                {
                                MeshCore::MeshGeomFacet Face(mPoints[mFacets[i]._aulPoints[0]],
                                                                             mPoints[mFacets[i]._aulPoints[1]],
                                                                                 mPoints[mFacets[i]._aulPoints[2]]);
                        
                            builder.AddFacet(Face);
                }
                else
                {
                            MeshCore::MeshGeomFacet Face2(mPoints[mFacets[i]._aulPoints[0]],
                                                                              mPoints[mFacets[i]._aulPoints[1]],
                                                                                  mPoints[mFacets[i]._aulPoints[2]]);
                        
                            builder2.AddFacet(Face2);
                }


        m_CadMesh.Assign(mPoints, mFacets);

        cout << "a" << endl;
        if (mFacets[i]._ulProp == 1 && !mFacets[i].IsFlag(MeshCore::MeshFacet::VISIT)) // falls unkritisches facet
            // welches noch nicht einer region zugewiesen wurde
        {

            cout << "b" << endl;
            m_RingVector.clear();
            m_RingVector.push_back(i);
            m_RegionVector.push_back(m_RingVector);
            m_RingVector.clear();

            m_RingCurrent = 0;

            m_CadMesh.VisitNeighbourFacets(*this,i);

            for (unsigned int j=0; j<m_RegionVector.size(); ++j)
            {
                cout << "c" << endl;
                for (unsigned int k=0; k<m_RegionVector[j].size(); ++k)
                {
                    aRegion.push_back(m_RegionVector[j][k]);

                }
            }

                        /*Base::Vector3f pnt1(0,0,0);
                        Base::Vector3f pnt2(1,0,0);
                        Base::Vector3f pnt3(0,1,0);
                        
                        MeshCore::MeshGeomFacet Face(pnt1,pnt2,pnt3);
                        builder.AddFacet(Face);*/

            m_Regions.push_back(aRegion);



            for (int l=0; l<num; ++l)            mFacets[l]._ucFlag = 0;

            cout << "d" << endl;

            for (unsigned int l=0; l<m_Regions.size(); ++l)
                for (unsigned int m=0; m<m_Regions[l].size(); ++m)
                    mFacets[m_Regions[l][m]].SetFlag(MeshCore::MeshFacet::VISIT);

    

            skals = RegionEvaluate(m_CadMesh, aRegion, normalsRef);
            RegionSkals.push_back(skals);
            skals.clear();

            m_RegionVector.clear();
            aRegion.clear();

            cout << "e" << endl;

        }
    }

        builder.Finish();
        builder2.Finish();

    cout << "fast" << endl;

    Base::Builder3D logg;
    double d;

    // stelle die regionen dar

    for (unsigned int i=0; i<RegionSkals.size(); ++i)
    {
        cout << i << endl;
        d = double(i)/double(RegionSkals.size()-1.0);
        cout << d << endl;
        for (unsigned int j=0; j<RegionSkals[i].size(); ++j)
        {

            logg.addSinglePoint((float) mPoints[RegionSkals[i][j].first].x,
                                (float) mPoints[RegionSkals[i][j].first].y,
                                (float) mPoints[RegionSkals[i][j].first].z,3,(float) d,(float) d,(float) d);


        }
    }

    logg.saveToFile("c:/regions.iv");





    /*
       // hole rand der regionen
       MeshCore::MeshRefPointToFacets p2fIt(tmpMesh);

       for(int i=0; i<m_Regions[0].size(); ++i)
       {
        for(int j=0; j<3; ++j)
        {
         std::set<MeshCore::MeshFacetArray::_TConstIterator>& faceSet = p2fIt[m_Regions[0][i]._aulPoints[j]];

                  for (std::set<MeshCore::MeshFacetArray::_TConstIterator>::const_iterator it = faceSet.begin(); it != faceSet.end(); ++it)
         {
          if((*it)[0]._ulProp == 0)
           m_RegionBounds.push_back((*it)[0]);
         }
        }
       }*/



    tmpVec = InitFaceCheck(m_CadMesh, deg_Tol+1);

    cout << " glob: " << tmpVec.size() << endl;

    // ********************************** LOKALE KORREKTUR ****************************************
    unsigned long ind;

    MeshCore::MeshPointArray tmpPnts = m_CadMesh.GetPoints();
    MeshCore::MeshFacetArray tmpFact = m_CadMesh.GetFacets();

    Base::Builder3D logic;
    int cm_ref = 50-cm;

    for (unsigned int i=0; i<RegionSkals.size(); ++i)
    {
        cm = cm_ref;

        while (true)
        {
            tmpVec = FaceCheck(m_CadMesh, deg_Tol+1);
            cout << "remaining: " << tmpVec.size() << endl;

            if (tmpVec.size() > NumRef || cm == 0)
            {
                tmpPnts = m_CadMesh.GetPoints();

                for (unsigned int k=0; k<tmpVec.size(); ++k)
                {
                    for (int l=0; l<3; ++l)
                    {
                        logic.addSinglePoint(tmpPnts[tmpFact[tmpVec[k]]._aulPoints[l]].x,
                                             tmpPnts[tmpFact[tmpVec[k]]._aulPoints[l]].y,
                                             tmpPnts[tmpFact[tmpVec[k]]._aulPoints[l]].z,5,0,0,1);
                    }
                }


                logic.saveToFile("c:/tired.iv");

                for (unsigned int j=0; j<RegionSkals[i].size(); ++j)
                {
                    ind = RegionSkals[i][j].first;
                    m_normals[ind].Normalize();
                    m_normals[ind].Scale(-normalsRef[ind].Length()/50,
                                         -normalsRef[ind].Length()/50,
                                         -normalsRef[ind].Length()/50);

                    tmpPnts[ind].Set(tmpPnts[ind].x + m_normals[ind].x,
                                     tmpPnts[ind].y + m_normals[ind].y,
                                     tmpPnts[ind].z + m_normals[ind].z);


                                        m_dist_vec[ind] += m_normals[ind];


                }

                m_CadMesh.Assign(tmpPnts, tmpFact);
                break;
            }

            tmpPnts = m_CadMesh.GetPoints();

            for (unsigned int j=0; j<RegionSkals[i].size(); ++j)
            {
                ind = RegionSkals[i][j].first;

                m_normals[ind].Normalize();
                m_normals[ind].Scale((normalsRef[ind].Length())/50,
                                     (normalsRef[ind].Length())/50,
                                     (normalsRef[ind].Length())/50);

                tmpPnts[ind].Set(tmpPnts[ind].x + m_normals[ind].x,
                                 tmpPnts[ind].y + m_normals[ind].y,
                                 tmpPnts[ind].z + m_normals[ind].z);

                                m_dist_vec[ind] += m_normals[ind];
            }

            m_CadMesh.Assign(tmpPnts, tmpFact);
            --cm;
        }
    }

    //SmoothMesh(m_CadMesh, 50);
        //SmoothMesh(m_CadMesh, 50);
        /*SmoothMesh(m_CadMesh, 50);
        SmoothMesh(m_CadMesh, 50);
        SmoothMesh(m_CadMesh, 50);
        SmoothMesh(m_CadMesh, 50);
        SmoothMesh(m_CadMesh, 50);
        SmoothMesh(m_CadMesh, 50);*/

//      MeshCore::MeshPointArray mpts = RefMesh23.GetPoints();

        //Base::Builder3D loooo;
        //ofstream anOutputFile;
        //anOutputFile.open("c:/catia_offset.txt");
 //   anOutputFile.precision(7);

        //for(int i=0; i< mpts.size(); ++i)
        //{
        //      loooo.addSingleArrow(mpts[i], mpts[i] + m_dist_vec[i]);
        //      anOutputFile <<  mpts[i].x << " " << mpts[i].y << " " << mpts[i].z << " " << m_dist_vec[i].x << " " << m_dist_vec[i].y << " " << m_dist_vec[i].z << endl;
        //}

        //loooo.saveToFile("c:/prpopo.iv");
        //  anOutputFile.close();




        //anOutputFile << 
        



    return true;
}

std::vector< std::pair<unsigned long, double> > SpringbackCorrection::RegionEvaluate(const MeshCore::MeshKernel &mesh, std::vector<unsigned long> &RegionFacets, std::vector<Base::Vector3f> &normals)
{
    std::list< std::vector  <Base::Vector3f> > Borders;
    std::list< std::vector  <Base::Vector3f> >::iterator bIt;
    std::map <unsigned long, Base::Vector3f> RegionMap;
    std::map <unsigned long, Base::Vector3f>::iterator rIt;
    std::map <double, unsigned long> DistMap;
    std::map <double, unsigned long>::iterator dIt;
    std::pair<unsigned long, Base::Vector3f> Ind2Pnt;
    std::pair<double, unsigned long> Dist2Ind;

    std::vector< std::pair<double, unsigned long> >  dists;
    std::vector< std::pair< unsigned long, double> > skals;

    MeshCore::MeshAlgorithm algo(mesh);  // übergebe mesh
    algo.GetFacetBorders(RegionFacets, Borders); // speichert ränder der region in Borders

    MeshCore::MeshFacetArray facets = mesh.GetFacets();
    MeshCore::MeshPointArray points = mesh.GetPoints();

    // um doppelte indizes zu vermeiden push die punkte der region in eine map
    for (unsigned int i=0; i<RegionFacets.size(); ++i)
    {
        for (int j=0; j<3; ++j)
        {
            Ind2Pnt.first  = facets[RegionFacets[i]]._aulPoints[j];
            Ind2Pnt.second = points[Ind2Pnt.first];

            RegionMap.insert(Ind2Pnt);
        }
                
    }

    // berechne distanzen der regionenpunkte zu den rändern

    Base::Vector3f distVec;   // speichert abstandsvektor
    double distVal;     // speichert maximalen abstandswert zum rand

    for (rIt = RegionMap.begin(); rIt != RegionMap.end(); ++rIt)
    {
        distVal = 1e+3;
        for (bIt = Borders.begin(); bIt != Borders.end(); ++bIt)
        {
            for (unsigned int i=0; i< bIt->size(); ++i)
            {
                distVec =  (*rIt).second - (*bIt).at(i);

                if (distVec.Length() < distVal)
                    distVal = distVec.Length();
            }
        }

        Dist2Ind.first   = distVal;
        Dist2Ind.second  = (*rIt).first;

        dists.push_back(Dist2Ind);
        DistMap.insert(Dist2Ind);
    }

    dIt = DistMap.end();
    --dIt;
    double d_max = (*dIt).first;
    std::pair<unsigned long, double> pair;

    for (unsigned int i=0; i<dists.size(); ++i)
    {
        pair.first  = dists[i].second;

        if (dists[i].first != 0.0)
            pair.second = ( (cos( ((dists[i].first * PI) / d_max) - PI) + 1) / 2);  // 0.5 * [ cos( (dist * PI / dist_max) - PI ) + 1 ]
        else
            pair.second = 0.0;

        skals.push_back(pair);
        normals[pair.first].Scale((float) pair.second, (float) pair.second, (float) pair.second);
    }

    return skals;
}



bool SpringbackCorrection::FacetRegionGrowing(MeshCore::MeshKernel &mesh,
        std::vector<MeshCore::MeshFacet> &FacetRegion,
        MeshCore::MeshFacetArray &mFacets)
{
    MeshCore::MeshRefFacetToFacets ff_It(mesh);

    MeshCore::MeshFacet facet = FacetRegion.back();
    const std::set<unsigned long>& FacetNei = ff_It[facet._ulProp];
    MeshCore::MeshFacetArray::_TConstIterator f_beg = mesh.GetFacets().begin();

    std::set<unsigned long>::const_iterator f_it;
    for (f_it = FacetNei.begin(); f_it != FacetNei.end(); ++f_it)
    {
        if (f_beg[*f_it]._ucFlag == MeshCore::MeshFacet::VISIT)
        {
            FacetRegion.push_back(f_beg[*f_it]);
            mFacets[f_beg[*f_it]._ulProp].SetFlag(MeshCore::MeshFacet::INVALID); // markiere als schon zugewiesen
            FacetRegionGrowing(mesh, FacetRegion, mFacets);
        }
    }

    return true;
}

//
//bool SpringbackCorrection::Init()
//{
// Base::Builder3D log;
// Base::Vector3f  p1,p2;
// double curv;
// // Triangulate Shape (vorläufig)
// best_fit::Tesselate_Shape(m_Shape, m_Mesh, 1);
//
// MeshCurvature(m_Mesh);
//
// std::vector<Base::Vector3f> normals = best_fit::Comp_Normals(m_Mesh);
// MeshCore::MeshPointIterator p_it(m_Mesh);
// int i = 0;
//
// for(p_it.Begin(); p_it.More(); p_it.Next()){
//  normals[i].Scale(m_CurvMax[i], m_CurvMax[i], m_CurvMax[i]);
//  //p1.Set(PointArray[i].x,       PointArray[i].y,       PointArray[i].z);
//  p2.Set(p_it->x + normals[i].x,   p_it->y + normals[i].y,   p_it->z + normals[i].z);
//
//
//  if(m_CurvMax[i] < 100)
//  log.addSingleLine(*p_it, p2, 1 ,0, 0, 0 );
//
//  //cout << p_it->x << ", " << p_it->y << ", " << p_it->z << endl;
//  //cout <<normals[i].x << ", " << normals[i].y << ", " << normals[i].z << endl;
//  //cout << m_CurvMax[i] << endl;
//  ++i;
//
// }
//
// log.saveToFile("c:/servolation.iv");
//}

// return true;
bool SpringbackCorrection::GetCurvature(TopoDS_Face aFace)
{
    Base::Builder3D log;
    Handle_Geom_Surface geom_surface;
    geom_surface = BRep_Tool::Surface(aFace);
    gp_Pnt proPnt;
    double u_par,v_par,k1,k2;
    gp_Vec D1U,D1V,D2U,D2V,D2UV;
    Base::Vector3f p1,p2;

    MeshCore::MeshPointArray PointArray = m_Mesh.GetPoints();
    MeshCore::MeshFacetArray FacetArray = m_Mesh.GetFacets();

    MeshCore::MeshPointIterator p_it(m_Mesh);

    m_CurvMax.resize(PointArray.size());

    for (unsigned int i=0; i<PointArray.size(); ++i)
    {

        proPnt.SetX(PointArray[i].x);
        proPnt.SetY(PointArray[i].y);
        proPnt.SetZ(PointArray[i].z);

        GeomAPI_ProjectPointOnSurf aProjection(proPnt,geom_surface);
        aProjection.LowerDistanceParameters(u_par,v_par);

        // Berechne Krümmung
        geom_surface->D2(u_par,v_par,proPnt,D1U,D1V,D2U,D2V,D2UV);

        // erste & zweite Hauptkrümmung
        k1 = D1U.CrossMagnitude(D2U) / D1U.Magnitude();
        k2 = D1V.CrossMagnitude(D2V) / D1V.Magnitude();

        /*if(k1>k2) PointArray[i].SetProperty(1/k1);
        else      PointArray[i].SetProperty(1/k2);*/

        if (abs(k1)>abs(k2))
        {
            m_CurvMax[i] = 1/k1;
        }
        else
            if (k2!=0)
                m_CurvMax[i] = 1/k2;
            else
                m_CurvMax[i] = 100;


        D1U.Cross(D1V);
        D1U.Normalize();
        D1U.Scale(m_CurvMax[i]);

        p1.x = (float) (proPnt.X() + D1U.X());
        p1.y = (float) (proPnt.Y() + D1U.Y());
        p1.z = (float) (proPnt.Z() + D1U.Z());

        p2.x = (float) proPnt.X();
        p2.y = (float) proPnt.Y();
        p2.z = (float) proPnt.Z();

        log.addSingleLine(p2, p1, 1 ,0, 0, 0 );
    }


    MeshCore::MeshPointIterator v_it(m_Mesh);
    MeshCore::MeshRefPointToPoints vv_it(m_Mesh);
    MeshCore::MeshPointArray::_TConstIterator v_beg = m_Mesh.GetPoints().begin();
    std::set<unsigned long> PntNei;
    std::set<unsigned long> PntNei2;
    std::set<unsigned long> PntNei3;
    std::set<unsigned long> PntNei4;
    std::set<unsigned long>::iterator pnt_it1;
    std::set<unsigned long>::iterator pnt_it2;
    std::set<unsigned long>::iterator pnt_it3;
    std::set<unsigned long>::iterator pnt_it4;
    std::vector<unsigned long> nei;
    double curv;

    std::vector<double> CurvCop = m_CurvMax;

    for (v_it.Begin(); v_it.More(); v_it.Next())
    {
        PntNei = vv_it[v_it.Position()];
        curv = m_CurvMax[v_it.Position()];

        for (pnt_it1 = PntNei.begin(); pnt_it1 !=PntNei.end(); ++pnt_it1)
        {
            if (m_CurvMax[v_beg[*pnt_it1]._ulProp] < curv)
                curv = m_CurvMax[v_beg[*pnt_it1]._ulProp];

            PntNei2 = vv_it[v_beg[*pnt_it1]._ulProp];
            for (pnt_it2 = PntNei2.begin(); pnt_it2 !=PntNei2.end(); ++pnt_it2)
            {
                if (m_CurvMax[v_beg[*pnt_it2]._ulProp] < curv)
                    curv = m_CurvMax[v_beg[*pnt_it2]._ulProp];


                PntNei3 = vv_it[v_beg[*pnt_it2]._ulProp];
                for (pnt_it3 = PntNei3.begin(); pnt_it3 !=PntNei3.end(); ++pnt_it3)
                {
                    if (m_CurvMax[v_beg[*pnt_it3]._ulProp] < curv)
                        curv = m_CurvMax[v_beg[*pnt_it3]._ulProp];

                    PntNei4 = vv_it[v_beg[*pnt_it3]._ulProp];
                    for (pnt_it4 = PntNei4.begin(); pnt_it4 !=PntNei4.end(); ++pnt_it4)
                    {
                        if (m_CurvMax[v_beg[*pnt_it4]._ulProp] < curv)
                            curv = m_CurvMax[v_beg[*pnt_it4]._ulProp];
                    }
                }
            }
        }

        CurvCop[v_it.Position()] = curv;
    }

    m_CurvMax = CurvCop;


    return true;
}


// 1. Computes the curvature of the mesh at the verticies
// 2. Stores the minimum curvatures along the neighbours for all verticies in a vector
std::vector<double> SpringbackCorrection::MeshCurvature(const TopoDS_Face& aFace, const MeshCore::MeshKernel& mesh)
{

        TopLoc_Location aLocation;
    Handle_Poly_Triangulation aTr = BRep_Tool::Triangulation(aFace,aLocation);
        const TColgp_Array1OfPnt2d& aUVNodes = aTr->UVNodes();
        int n = aUVNodes.Length();
        BRepAdaptor_Surface aSurface(aFace);
    
    gp_Pnt2d par;
    gp_Pnt P;
    gp_Vec D1U, D1V, D2U, D2V, D2UV, nor, xvv, xuv, xuu;
        double H,K;   // Gaußsche- und mittlere Krümmung

        std::vector<double> aMaxCurve, aMinCurve;
        std::vector<double> curv(2);
    
        for (int i=0; i<n; ++i)
    {
        par = aUVNodes.Value(i+1);
                aSurface.D2(par.X(),par.Y(),P,D1U,D1V,D2U,D2V,D2UV);

                //berechne Hilfsnormale
                nor = D1U;
                nor.Cross(D1V);
                nor.Normalize();

                xvv = D2V;
                xuv = D2UV;
                xuu = D2U;

                xvv.Multiply(D1U*D1U);
                xuv.Multiply(2*(D1U*D1V));
                xuu.Multiply(D1V*D1V);

                H = ((xvv - xuv + xuu)*nor);
                H /= 2*((D1U*D1U)*(D1V*D1V) - (D1U*D1V)*(D1U*D1V));
                K = ((nor*D2U)*(nor*D2V) - (nor*D2UV)*(nor*D2UV))/((D1U*D1U)*(D1V*D1V) - (D1U*D1V)*(D1U*D1V));


                aMaxCurve.push_back(-(H - sqrt(H*H - K)));
                aMinCurve.push_back(-(H + sqrt(H*H - K)));

        }
                
                
                
        double maxCurv = 0.0;
        double avgCurv = 0.0;
        double tmp1 =  1e+10;
        double tmp2 = -1e+10;

        double mean_curvMax = 0.0, mean_curvMin = 0.0;
    int c1 = 0, c2 = 0;

        m_CurvPos.clear();
        m_CurvNeg.clear();
        m_CurvMax.clear();
        m_CurvPos.resize(mesh.CountPoints());
        m_CurvNeg.resize(mesh.CountPoints());
        m_CurvMax.resize(mesh.CountPoints());

        for ( unsigned long i=0; i<n; i++ )
        {
                if (aMaxCurve[i] > 0) m_CurvPos[i] = 1 / aMaxCurve[i];
                else                  m_CurvPos[i] = 1e+10;

                if (aMinCurve[i] < 0) m_CurvNeg[i] = 1 / aMinCurve[i];
                else                  m_CurvNeg[i] = -1e+10;

                if (m_CurvPos[i] < tmp1)
                        tmp1 = m_CurvPos[i];

                if (m_CurvNeg[i] > tmp2)
                        tmp2 = m_CurvNeg[i];

                if (aMaxCurve[i] > 0)
                {
                        ++c1;
                        mean_curvMax += 1/aMaxCurve[i];
                }

                if (aMinCurve[i] < 0)
                {
                        ++c2;
                        mean_curvMin += 1/aMinCurve[i];
                }
        }

        if (c1==0)
                mean_curvMax = 1e+3;
        else
                mean_curvMax /= mesh.CountPoints();

        if (c2==0)
                mean_curvMin = -1e+3;
        else
                mean_curvMin /= mesh.CountPoints();

        curv[0] = tmp1;
        curv[1] = tmp2;
        

        
        
        //Krümmung auf Netzbasis
        
        for(int i=0; i<n; ++i)
        {

                // get all points
            
                //MeshCore::MeshKernel& rMesh = mesh;
                std::vector< Wm4::Vector3<float> > aPnts;
                std::vector<Base::Vector3f> aPnts_tmp;
                MeshCore::MeshPointIterator cPIt( mesh );
                for ( cPIt.Init(); cPIt.More(); cPIt.Next() )
                {
                        Wm4::Vector3<float> cP( cPIt->x, cPIt->y, cPIt->z );
                        Base::Vector3f cP_tmp( cPIt->x, cPIt->y, cPIt->z );
                        aPnts.push_back( cP );
                        aPnts_tmp.push_back(cP_tmp);
                }

                // get all point connections
                std::vector<int> anIdx;
                const std::vector<MeshCore::MeshFacet>& MeshFacetArray = mesh.GetFacets();
                for ( std::vector<MeshCore::MeshFacet>::const_iterator jt = MeshFacetArray.begin(); jt != MeshFacetArray.end(); ++jt )
                {
                        for (int i=0; i<3; i++)
                        {
                                anIdx.push_back( (int)jt->_aulPoints[i] );
                        }
                }

                // compute vertex based curvatures
                Wm4::MeshCurvature<float> meshCurv(mesh.CountPoints(), &(aPnts[0]), mesh.CountFacets(), &(anIdx[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();

                double maxCurv = 0.0;
                double avgCurv = 0.0;
                double tmp1 =  1e+10;
                double tmp2 = -1e+10;

                double mean_curvMax = 0.0, mean_curvMin = 0.0;
                int c1 = 0, c2 = 0;

                m_CurvPos.clear();
                m_CurvNeg.clear();
                m_CurvMax.clear();
                m_CurvPos.resize(mesh.CountPoints());
                m_CurvNeg.resize(mesh.CountPoints());
                m_CurvMax.resize(mesh.CountPoints());

                for ( unsigned long i=0; i<mesh.CountPoints(); i++ )
                {

                        if (aMaxCurv[i] > 0) m_CurvPos[i] = 1 / aMaxCurv[i];
                        else                 m_CurvPos[i] = 1e+10;

                        if (aMinCurv[i] < 0) m_CurvNeg[i] = 1 / aMinCurv[i];
                        else                 m_CurvNeg[i] = -1e+10;

                        if (m_CurvPos[i] < tmp1)
                                tmp1 = m_CurvPos[i];

                        if (m_CurvNeg[i] > tmp2)
                                tmp2 = m_CurvNeg[i];

                        if (aMaxCurv[i] > 0)
                        {
                                ++c1;
                                mean_curvMax += 1/aMaxCurv[i];
                        }

                        if (aMinCurv[i] < 0)
                        {
                                ++c2;
                                mean_curvMin += 1/aMinCurv[i];
                        }
                }

                if (c1==0)
                        mean_curvMax = 1e+3;
                else
                        mean_curvMax /= mesh.CountPoints();

                if (c2==0)
                        mean_curvMin = -1e+3;
                else
                        mean_curvMin /= mesh.CountPoints();

                curv[0] = tmp1;
                curv[1] = tmp2;
        }
        

    return curv;
}

/*
bool SpringbackCorrection::MirrorMesh(std::vector<double> error)
{
    // Flags: 0 - not yet checked
    //       1 - no correction applied, but neccessary
    //        2 - first correction applied and still neccessary
    //        3 - done successfully
    //        4 - done without success
    //        5 - boundary point

    int n = error.size();
    m_error = error;

    std::vector<int> InterPnts;
    std::vector<unsigned long> InterFace;
    std::vector<unsigned long> IFCopy;

    bool boo = true;

    std::vector<unsigned long> ind;
    std::vector<unsigned long> copy;

    MeshRef = m_Mesh;
    std::vector<Base::Vector3f> NormalsRef = m_normals;
    std::vector<Base::Vector3f> NormalOrig = m_normals;

    PointArray = MeshRef.GetPoints();
    MeshCore::MeshPointArray PntArray   = m_Mesh.GetPoints();
    MeshCore::MeshFacetArray FctArray   = m_Mesh.GetFacets();
    MeshCore::MeshFacetArray FacetArray = MeshRef.GetFacets();

    MeshCore::MeshFacetIterator f_it(MeshRef);
    MeshCore::MeshPointIterator v_it(MeshRef);
    MeshCore::MeshRefPointToPoints vv_it(MeshRef);
    MeshCore::MeshRefPointToFacets vf_it(MeshRef);

    std::set<MeshCore::MeshPointArray::_TConstIterator> PntNei;
    std::set<MeshCore::MeshFacetArray::_TConstIterator> FacetNei;

    // Iterates in facets in form of MeshGeomFacets
    std::set<MeshCore::MeshFacetArray::_TConstIterator>::iterator face_it;
    std::set<MeshCore::MeshPointArray::_TConstIterator>::iterator pnt_it;

    std::vector<Base::Vector3f> NeiLoc;
    std::vector<Base::Vector3f> NorLoc;
    std::vector<Base::Vector3f> TmpPntVec;
    std::vector<unsigned long>  TmpIndVec;
    std::vector<unsigned long>  nei;

    bool k=true;
    int count = 0;

    while (k==true)
    {
        // ------ Self-Intersection Check + Correction ------------------------------------------------------

        if (count ==2) break;
        count++;

        Base::Builder3D log3d;
        Base::Builder3D loga;
        Base::Vector3f Pnt, NeiPnt, NeiPntNext, NeiPntPrev, TmpPnt, TmpNor;

        std::vector<Base::Vector3f> Neis;
        double minScale, tmp;
        bool sign;

        ublas::matrix<double> A(3, 3);
        ublas::matrix<double> b(1, 3);

        ind.clear();
        for (unsigned int i=0; i<m_Mesh.CountPoints(); ++i)
            ind.push_back(i);


        // do correction for selected points ---------------------------------
        for (unsigned int j=0; j<ind.size(); ++j)
        {
            sign = (error[ind[j]] > 0);  // (+) -> true , (-) -> false

            MeshRef.Assign(PointArray, FacetArray);

            v_it.Set(ind[j]);

            Pnt      = *v_it;
            PntNei   = vv_it[ind[j]];
            FacetNei = vf_it[ind[j]];

            nei.clear();

            // falls boundary-point tue nichts bzw. gehe zum nächsten punkt
            if (FacetNei.size() != PntNei.size())
            {
                PointArray[ind[j]]._ucFlag = 5;
                m_CurvMax[ind[j]] = abs(error[ind[j]]);   // vorübergehend
                continue;
            }

            // ordne nachbarn
            ReorderNeighbourList(PntNei,FacetNei,nei,ind[j]);

            NeiLoc.clear();
            NorLoc.clear();
            NeiLoc.push_back(Pnt);
            NorLoc.push_back(NormalsRef[ind[j]]);

            for (unsigned int r=0; r<nei.size(); ++r)
            {
                v_it.Set(nei[r]);
                NeiLoc.push_back(*v_it);
                NorLoc.push_back(NormalsRef[nei[r]]);
            }


            if (count == 2)
            {
                if (sign ==true) minScale = -(1e+10);
                else      minScale =   1e+10;
                minScale = LocalCorrection(NeiLoc,NorLoc,sign,minScale);

                if ((abs(0.8*minScale) - abs(error[ind[j]])) > 0)
                {
                    PointArray[ind[j]]._ucFlag = 3;
                    m_CurvMax[ind[j]] = abs(minScale);
                    //cout << "t, ";
                    continue;
                }
                else
                {
                    PointArray[ind[j]]._ucFlag = 4;
                    m_CurvMax[ind[j]] = abs(minScale);
                    cout << "f, ";
                    continue;
                }
            }


            if (sign ==true) minScale = -(1e+10);
            else      minScale =   1e+10;

            minScale = LocalCorrection(NeiLoc,NorLoc,sign,minScale);

            if ((abs(0.8*minScale) - abs(error[ind[j]])) > 0)
            {
                PointArray[ind[j]]._ucFlag = 3;
                m_CurvMax[ind[j]] = abs(minScale);
                continue;
            }

            //---------------------------- 1. Korrektur -------------------------
            tmp = 0;
            std::vector<Base::Vector3f> ConvComb = FillConvex(NeiLoc,NorLoc,20);

            for (unsigned int l=0; l<ConvComb.size(); ++l)
            {
                if (sign ==true) minScale = -(1e+10);
                else      minScale =   1e+10;

                NeiLoc[0] = ConvComb[l];
                minScale = LocalCorrection(NeiLoc,NorLoc,sign,minScale);

                if (abs(minScale) > abs(tmp))
                {
                    tmp = minScale;
                    TmpPnt = NeiLoc[0];
                }
            }

            NeiLoc[0] = PointArray[ind[j]];
            cout << "new minScale: " << tmp << endl;

            if ((abs(0.8*tmp) - abs(error[ind[j]])) > 0)
            {
                PointArray[ind[j]].Set(TmpPnt.x, TmpPnt.y, TmpPnt.z);
                PointArray[ind[j]]._ucFlag = 3;
                m_CurvMax[ind[j]] = abs(tmp);
                continue;
            }
            else
            {
                PointArray[ind[j]]._ucFlag = 2;
                m_CurvMax[ind[j]] = abs(tmp);
            }

            //------------------------------- END -------------------------------

            // update
            for (unsigned int r=0; r<nei.size(); ++r)
            {
                v_it.Set(nei[r]);
                NeiLoc[r+1] = *v_it;
            }

            MeshRef.Assign(PointArray, FacetArray);

            //---------------------------- 2. Korrektur -------------------------
            tmp = GlobalCorrection(NeiLoc,NorLoc,sign,ind[j]);

            //MeshRef.Assign(PointArray, FacetArray);
            //m_Mesh = MeshRef;
            //return true;

            cout << "2: " << tmp << endl;

            if ((abs(0.8*tmp) - abs(error[ind[j]])) > 0)
            {
                PointArray[ind[j]].Set(NeiLoc[0].x, NeiLoc[0].y, NeiLoc[0].z);
                NormalsRef[ind[j]].Set(NorLoc[0].x, NorLoc[0].y, NorLoc[0].z);
                PointArray[ind[j]]._ucFlag = 3;
                m_CurvMax[ind[j]] = abs(tmp);
                continue;
            }
            else
            {
                if (m_CurvMax[ind[j]] < abs(tmp))
                {
                    m_CurvMax[ind[j]] = abs(tmp);
                    PointArray[ind[j]].Set(NeiLoc[0].x, NeiLoc[0].y, NeiLoc[0].z);
                    NormalsRef[ind[j]].Set(NorLoc[0].x, NorLoc[0].y, NorLoc[0].z);
                }
                // flag bleibt 2
            }

            //------------------------------- END -------------------------------

        } // end for

        cout << "Replace Mesh..." << endl;
        PntArray   = PointArray;
        m_normals  = NormalsRef;
        m_Mesh = MeshRef;

    } // end while


    double ver, mx=100;
    for (unsigned int i=0; i<error.size(); ++i)
    {
        ver = m_CurvMax[i] / abs(error[i]);
        if (ver < mx)
            mx = ver;
    }

    //cout << "skalierung: " << mx << endl;
    //for(int i=0; i<error.size(); ++i)
    //error[i] *= mx;



    // deform CAD-Mesh with regard to the scaled normals
    MeshCore::MeshPointIterator p_it(m_Mesh);
    Base::Builder3D log,log2;
    std::vector<int> IndexVec,NumVec, NumVecCopy;
    int c, lastNum, lastInd;

    // kritische punkte plus anzahl seiner kritischen nachbarn speichern
    //ind.clear();
    for (unsigned int j=0; j<ind.size(); ++j)
    {
        if (PointArray[ind[j]]._ucFlag == 3 || PointArray[ind[j]]._ucFlag == 5)
        {
            //ind.push_back(ind[j]);
            continue;
        }
        else
        {
            c=0;
            PntNei = vv_it[ind[j]];
            for (pnt_it = PntNei.begin(); pnt_it != PntNei.end(); ++pnt_it)
            {
                if (PointArray[(*pnt_it)->_ulProp]._ucFlag == 4)
                    c++;
            }

            IndexVec.push_back(ind[j]);
            NumVec.push_back(c);
        }
    }

    // normalenskalierung
    for (unsigned int j=0; j<ind.size(); ++j)
    {
        if ( ( abs(0.8*m_CurvMax[ind[j]]) - abs(error[ind[j]]) ) < 0)
        {
            if (error[ind[j]] < 0)
                m_normals[ind[j]].Scale((float) (0.8*m_CurvMax[ind[j]]),(float) (0.8*m_CurvMax[ind[j]]),(float) (0.8*m_CurvMax[ind[j]]));
            else
                m_normals[ind[j]].Scale((float) (-0.8*m_CurvMax[ind[j]]),(float) (-0.8*m_CurvMax[ind[j]]), (float) (-0.8*m_CurvMax[ind[j]]));
        }
        else
            m_normals[ind[j]].Scale((float) -error[ind[j]],(float)  -error[ind[j]],(float)  -error[ind[j]]);
    }


    for (unsigned int j=0; j<ind.size(); ++j)
    {
        log.addSingleLine(PointArray[ind[j]],PointArray[ind[j]] + m_normals[ind[j]],2,1,1,1);  // white arrows
        PntArray[ind[j]].Set(PntArray[ind[j]].x + m_normals[ind[j]].x,
                             PntArray[ind[j]].y + m_normals[ind[j]].y,
                             PntArray[ind[j]].z + m_normals[ind[j]].z);
    }

    m_Mesh.Assign(PntArray, FctArray);
    MeshRef = m_Mesh;
    log.saveToFile("c:/deformation.iv");


    // kritische punkte nach der anzahl seiner kritischen nachbarn sortieren
    NumVecCopy = NumVec;
    std::vector<int> SortNumVec(NumVec.size(), -1);
    std::vector<int> SortIndVec(NumVec.size(), -1);
    cout << NumVec.size() << endl;

    //cout << "vor dem sortieren: ";
    //for(int j=0; j<NumVec.size(); ++j)
    //{
    //cout << NumVec[j] << ", ";
    //}
    //cout << endl;

    for (unsigned int j=0; j<NumVecCopy.size(); ++j)
    {

        lastNum = NumVec.back();
        lastInd = IndexVec.back();

        NumVec.pop_back();
        IndexVec.pop_back();

        c=0;
        for (unsigned int i=0; i<NumVecCopy.size(); ++i)
        {
            if (NumVecCopy[i] < lastNum)
                ++c;
        }

        while (SortNumVec[c] == lastNum)
            ++c;

        SortNumVec[c] = lastNum;
        SortIndVec[c] = lastInd;
    }

    //cout << "danach: ";
    //for(int j=0; j<SortNumVec.size(); ++j)
    //{
    //cout << SortNumVec[j] << ", ";
    //}
    //cout << endl;

    NeiLoc.clear();
    NorLoc.clear();

    Base::Builder3D log1;

    MeshCore::MeshKernel mesh;
    MeshCore::MeshPoint mpnt;
    MeshCore::MeshPointArray mPointArray;
    MeshCore::MeshFacetArray mFacetArray;

    float w;

    for (unsigned int j=0; j<SortIndVec.size(); ++j)
    {
        Base::Vector3f bvec(0.0, 0.0, 0.0);
        NeiLoc.clear();
        NeiLoc.push_back(bvec);

        PntNei = vv_it[SortIndVec[j]];
        FacetNei = vf_it[SortIndVec[j]];
        nei.clear();

        // ordne nachbarn
        ReorderNeighbourList(PntNei,FacetNei,nei,SortIndVec[j]);

        mPointArray.clear();
        for (unsigned int i=0; i<nei.size(); ++i)
        {
            mpnt = PntArray[nei[i]];
            mPointArray.push_back(mpnt);
            if (PntArray[nei[i]]._ucFlag != 4)
            {
                NeiLoc.push_back(PntArray[nei[i]]);
                log1.addSinglePoint(PntArray[nei[i]].x, PntArray[nei[i]].y, PntArray[nei[i]].z, 2,1,1,1);
                mPointArray.push_back(mpnt);
            }
            else
            {
                bvec += PntArray[nei[i]];
                log1.addSinglePoint(PntArray[nei[i]].x, PntArray[nei[i]].y, PntArray[nei[i]].z, 2,1,0,0);
            }
        }

        log1.saveToFile("c:/nachbarn.iv");

        mesh.Assign(mPointArray, mFacetArray);

        // need at least one facet to perform the PCA
        MeshCore::MeshGeomFacet face;
        mesh.AddFacet(face);

        // berechne hauptachsen
        MeshCore::MeshEigensystem pca(mesh);
        pca.Evaluate();
        Base::Matrix4D T1 = pca.Transform();
        std::vector<Base::Vector3f> v(3);

        v[0].x = (float) T1[0][0];
        v[0].y = (float) T1[0][1];
        v[0].z = (float) T1[0][2];
        v[1].x = (float) T1[1][0];
        v[1].y = (float) T1[1][1];
        v[1].z = (float) T1[1][2];
        v[2].x = (float) T1[2][0];
        v[2].y = (float) T1[2][1];
        v[2].z = (float) T1[2][2];

        v[0].Normalize();
        v[1].Normalize();
        v[2].Normalize();

        double deg,temp=0;
        int s;

        for (int i=0; i<3; ++i)
        {
            deg = v[i]*NormalOrig[SortIndVec[j]];
            if (abs(deg)>temp)
            {
                temp = abs(deg);
                s=i;
            }
        }

        if (v[s].z < 0.0)
            v[s] *= -1;

        for (unsigned int i=0; i<NeiLoc.size(); ++i)
            NorLoc.push_back(v[s]);

        // gewichtung (im verhältnis 1:2)
        w = float(1.0) /(float(NeiLoc.size()) + float(PntNei.size()) - float(1.0));
        bvec.Scale(w,w,w);

        for (unsigned int i=1; i<NeiLoc.size(); ++i)
        {
            bvec.x = bvec.x + 2*w*NeiLoc[i].x;
            bvec.y = bvec.y + 2*w*NeiLoc[i].y;
            bvec.z = bvec.z + 2*w*NeiLoc[i].z;
        }

        NeiLoc[0] = bvec;
        log1.addSinglePoint(NeiLoc[0],2,0,1,0);
        std::vector<Base::Vector3f> Convs = FillConvex(NeiLoc, NorLoc, 20);

        for (unsigned int i=0; i<Convs.size(); ++i)
            log1.addSinglePoint(Convs[i],2,0,0,0);

        log1.saveToFile("c:/nachbarn.iv");

        PntArray[SortIndVec[j]].Set(Convs[0].x, Convs[0].y, Convs[0].z);
        PntArray[SortIndVec[j]]._ucFlag = 3;

    }

    m_Mesh.Assign(PntArray, FctArray);

    return true;


    for (p_it.Begin(); p_it.More(); p_it.Next())
    {
        if (PointArray[p_it.Position()]._ucFlag == 3)
        {
            log.addSingleLine(*p_it,*p_it + m_normals[p_it.Position()],2,1,1,1);  // red arrows
            PntArray[p_it.Position()].Set(PntArray[p_it.Position()].x + m_normals[p_it.Position()].x,
                                          PntArray[p_it.Position()].y + m_normals[p_it.Position()].y,
                                          PntArray[p_it.Position()].z + m_normals[p_it.Position()].z);
        }

        if (PointArray[p_it.Position()]._ucFlag == 4)
        {
            log.addSingleLine(*p_it,*p_it + m_normals[p_it.Position()],2,1,0,0);  // red arrows
            PntArray[p_it.Position()].Set(PntArray[p_it.Position()].x + m_normals[p_it.Position()].x,
                                          PntArray[p_it.Position()].y + m_normals[p_it.Position()].y,
                                          PntArray[p_it.Position()].z + m_normals[p_it.Position()].z);
        }
    }

    m_Mesh.Assign(PntArray, FctArray);
    log.saveToFile("c:/deformation.iv");

    return true;
}
*/

/*
double SpringbackCorrection::LocalCorrection(std::vector<Base::Vector3f> Neib ,
        std::vector<Base::Vector3f> Normal, bool sign, double minScale)
{
    Base::Vector3f Pnt, NeiPnt, NeiPntNext, NeiPntPrev;

    ublas::matrix<double> A(3, 3);
    ublas::matrix<double> b(1, 3);

    // 2. Spalte bleibt unverändert
    A(0,2) = -Normal[0].x;
    A(1,2) = -Normal[0].y;
    A(2,2) = -Normal[0].z;

    // forward
    for (unsigned int k=1; k<Neib.size(); ++k)
    {
        NeiPnt = Neib[k];

        if ((k+1) == Neib.size())
        {
            NeiPntNext = Neib[1];
        }
        else
        {
            NeiPntNext = Neib[k+1];
        }

        A(0,2) = -Normal[0].x;
        A(0,0) = Normal[k].x;
        A(0,1) = (NeiPntNext.x - NeiPnt.x);
        A(1,2) = -Normal[0].y;
        A(1,0) = Normal[k].y;
        A(1,1) = (NeiPntNext.y - NeiPnt.y);
        A(2,2) = -Normal[0].z;
        A(2,0) = Normal[k].z;
        A(2,1) = (NeiPntNext.z - NeiPnt.z);

        b(0,0) = Neib[0].x - NeiPnt.x;
        b(0,1) = Neib[0].y - NeiPnt.y;
        b(0,2) = Neib[0].z - NeiPnt.z;

        //cout << "A-Matrix: " << endl;
        //cout << A(0,0) << ", " << A(0,1) << ", " << A(0,2) << endl;
        //cout << A(1,0) << ", " << A(1,1) << ", " << A(1,2) << endl;
        //cout << A(2,0) << ", " << A(2,1) << ", " << A(2,2) << endl;
        //cout << endl;

        //cout << "b-Vector: " << endl;
        //cout << b(0,0) << ", " << b(0,1) << ", " << b(0,2) << endl;
        //cout << endl;

        atlas::gesv(A,b);

        //cout << "solution: " << endl;
        //cout << b(0,0) << ", " << b(0,1) << ", " << b(0,2) << endl;
        //cout << endl;

        if ((b(0,2)<=0) && (sign == true) && (b(0,1)>=0) && (b(0,1)<=1))
        {
            if (minScale < b(0,2))
                minScale = b(0,2);
        }
        else if (b(0,2)>=0 && sign == false && (b(0,1)>=0) && (b(0,1)<=1))
        {
            if (minScale > b(0,2))
                minScale = b(0,2);
        }
    }

    // backward
    for (unsigned int k=1; k<Neib.size(); ++k)
    {
        NeiPnt = Neib[k];

        if (k == 1)
        {
            NeiPntPrev = Neib[Neib.size()-1];
        }
        else
        {
            NeiPntPrev = Neib[k-1];
        }

        A(0,2) = -Normal[0].x;
        A(0,0) = Normal[k].x;
        A(0,1) = (NeiPntPrev.x - NeiPnt.x);
        A(1,2) = -Normal[0].y;
        A(1,0) = Normal[k].y;
        A(1,1) = (NeiPntPrev.y - NeiPnt.y);
        A(2,2) = -Normal[0].z;
        A(2,0) = Normal[k].z;
        A(2,1) = (NeiPntPrev.z - NeiPnt.z);

        b(0,0) = Neib[0].x - NeiPnt.x;
        b(0,1) = Neib[0].y - NeiPnt.y;
        b(0,2) = Neib[0].z - NeiPnt.z;

        //cout << "A-Matrix: " << endl;
        //cout << A(0,0) << ", " << A(0,1) << ", " << A(0,2) << endl;
        //cout << A(1,0) << ", " << A(1,1) << ", " << A(1,2) << endl;
        //cout << A(2,0) << ", " << A(2,1) << ", " << A(2,2) << endl;
        //cout << endl;

        //cout << "b-Vector: " << endl;
        //cout << b(0,0) << ", " << b(0,1) << ", " << b(0,2) << endl;
        //cout << endl;

        atlas::gesv(A,b);

        //cout << "solution: " << endl;
        //cout << b(0,0) << ", " << b(0,1) << ", " << b(0,2) << endl;
        //cout << endl;

        if (b(0,2)<0 && sign == true && (b(0,1)>=0) && (b(0,1)<=1))
        {
            if (minScale < b(0,2))
                minScale = b(0,2);
        }
        else if (b(0,2)>0 && sign == false && (b(0,1)>=0) && (b(0,1)<=1))
        {
            if (minScale > b(0,2))
                minScale = b(0,2);
        }
    }

    return minScale;
}
*/

/*
std::vector<Base::Vector3f> SpringbackCorrection::FillConvex(std::vector<Base::Vector3f> Neib, std::vector<Base::Vector3f> Normals, int n)
{
    std::vector<Base::Vector3f> ConvComb;
    Base::Vector3f tmp;
    for (unsigned int i=1; i<Neib.size(); ++i)
    {
        for (int j=0; j<n; ++j)
        {
            //cout << "t, " ;
            tmp = Neib[i] - Neib[0];
            tmp.Scale(float(j)/float(n), float(j)/float(n), float(j)/float(n));
            if (InsideCheck(Neib[0] + tmp, Normals[0], Neib) == true)
                ConvComb.push_back(Neib[0] + tmp);
            else;
            //cout << "f, " ;
        }
    }
    cout << endl;
    return ConvComb;
}
*/

/*
bool SpringbackCorrection::InsideCheck(Base::Vector3f pnt, Base::Vector3f normal, std::vector<Base::Vector3f> Neib)
{
    Base::Vector3f Pnt, NeiPnt, NeiPntNext, Cross;

    ublas::matrix<double> A(3, 3);
    ublas::matrix<double> b(1, 3);

    bool b1 = true, b2;

    for (unsigned int k=1; k<Neib.size(); ++k)
    {
        NeiPnt = Neib[k];

        if ((k+1) == Neib.size())
        {
            NeiPntNext = Neib[1];
        }
        else
        {
            NeiPntNext = Neib[k+1];
        }

        Cross = normal%(NeiPntNext - NeiPnt);

        A(0,0) = normal.x;
        A(0,1) = (NeiPntNext.x - NeiPnt.x);
        A(0,2) = -Cross.x;
        A(1,0) = normal.y;
        A(1,1) = (NeiPntNext.y - NeiPnt.y);
        A(1,2) = -Cross.y;
        A(2,0) = normal.z;
        A(2,1) = (NeiPntNext.z - NeiPnt.z);
        A(2,2) = -Cross.z;

        b(0,0) = pnt.x - NeiPnt.x;
        b(0,1) = pnt.y - NeiPnt.y;
        b(0,2) = pnt.z - NeiPnt.z;

        atlas::gesv(A,b);

        if (b1==true)
        {
            if ( b(0,2) < 0)
            {
                b1 = false;
                b2 = true;
            }
            else if (b(0,2)>0)
            {
                b1 = false;
                b2 = false;
            }
        }

        if (b1==false)
        {
            if ((b2 == true)  && (b(0,2) > 0)) return false;
            if ((b2 == false) && (b(0,2) < 0)) return false;
        }
    }
    return true;
}
*/

/*
double SpringbackCorrection::GlobalCorrection(std::vector<Base::Vector3f> NeiLoc, std::vector<Base::Vector3f> NorLoc,
        bool sign, int ind)
{
    double minScale,mScRef,tmp = 0;
    Base::Vector3f TmpPnt, TmpNor;

    if (sign ==true) mScRef = -(1e+10);
    else      mScRef =   1e+10;

    minScale = mScRef;
    minScale = LocalCorrection(NeiLoc,NorLoc,sign,minScale);

    if ((abs(0.8*minScale) - abs(m_error[ind])) > 0)
    {
        PointArray[ind]._ucFlag = 3;
        m_CurvMax[ind] = abs(minScale);
        return minScale;
    }

    std::vector<Base::Vector3f> ConvComb = FillConvex(NeiLoc,NorLoc,20);

    // füge schwerpunkte der umgebenden faces hinzu
    std::set<MeshCore::MeshFacetArray::_TConstIterator> FacetNei;
    std::set<MeshCore::MeshFacetArray::_TConstIterator>::iterator face_it;
    MeshCore::MeshFacetIterator    f_it (MeshRef);
    MeshCore::MeshRefPointToFacets vf_it(MeshRef);

    FacetNei = vf_it[ind];

    for (face_it = FacetNei.begin(); face_it != FacetNei.end(); ++face_it)
    {
        TmpPnt.Set(0.0, 0.0, 0.0);
        TmpPnt += PointArray[(*face_it)->_aulPoints[0]];
        TmpPnt += PointArray[(*face_it)->_aulPoints[1]];
        TmpPnt += PointArray[(*face_it)->_aulPoints[2]];
        TmpPnt.Scale((float) (1.0/3.0),(float) (1.0/3.0),(float) (1.0/3.0));
        ConvComb.push_back(TmpPnt);
    }

    //Base::Builder3D logg;
    //for(int i=0; i<NeiLoc.size(); ++i)
    // logg.addSinglePoint(NeiLoc[i], 2,1,0,0);

    //for(int i=0; i<ConvComb.size(); ++i)
    // logg.addSinglePoint(ConvComb[i]);

    //logg.saveToFile("c:/convex.iv");

    for (unsigned int i=0; i<ConvComb.size(); ++i)
    {
        for (unsigned int j=1; j<NorLoc.size()-1; ++j)
        {

            minScale  = mScRef;
            NeiLoc[0] = ConvComb[i];
            NorLoc[0] = NorLoc[j];
            minScale  = LocalCorrection(NeiLoc,NorLoc,sign,minScale);


            if (minScale > tmp)
            {
                tmp = minScale;
                TmpPnt = NeiLoc[0];
                TmpNor = NorLoc[0];
            }
        }
    }

    NeiLoc[0] = TmpPnt;
    NorLoc[0] = TmpNor;

    if ((abs(0.8*tmp) - abs(m_error[ind])) > 0)
    {
        PointArray[ind]._ucFlag = 3;
        m_CurvMax[ind] = abs(tmp);
    }
    else
    {
        PointArray[ind]._ucFlag = 4;
        m_CurvMax[ind] = abs(tmp);
    }

    return tmp;
}
/*

/*
void SpringbackCorrection::ReorderNeighbourList(std::set<MeshCore::MeshPointArray::_TConstIterator> &pnt,
        std::set<MeshCore::MeshFacetArray::_TConstIterator> &face, std::vector<unsigned long> &nei, unsigned long CurInd)
{
    std::set<MeshCore::MeshPointArray::_TConstIterator>::iterator pnt_it;
    std::set<MeshCore::MeshFacetArray::_TConstIterator>::iterator face_it;
    std::vector<unsigned long>::iterator vec_it;
    std::vector<unsigned long>::iterator ulong_it;
    unsigned long PrevIndex;

    if (face.size() != pnt.size())  //this means boundary point
    {
        for (pnt_it = pnt.begin(); pnt_it != pnt.end(); ++pnt_it)
        {
            int c = 0;
            for (face_it = face.begin(); face_it != face.end(); ++face_it)
            {
                if ((*pnt_it)[0]._ulProp == (*face_it)[0]._aulPoints[0])
                    ++c;
                if ((*pnt_it)[0]._ulProp == (*face_it)[0]._aulPoints[1])
                    ++c;
                if ((*pnt_it)[0]._ulProp == (*face_it)[0]._aulPoints[2])
                    ++c;
            }

            if (c==1)
                break;
        }
    }
    else
        pnt_it = pnt.begin();



    face_it = face.begin();

    nei.push_back((*pnt_it)->_ulProp);  //push back first neighbour
    vec_it = nei.begin();
    PrevIndex = nei[0];  //Initialize PrevIndex
    for (unsigned int i = 1; i < pnt.size(); i++) //Start
    {
        while (true)
        {
            std::vector<unsigned long> facetpnt;
            facetpnt.push_back((*face_it)->_aulPoints[0]); //push back into a vector for easier iteration
            facetpnt.push_back((*face_it)->_aulPoints[1]);
            facetpnt.push_back((*face_it)->_aulPoints[2]);
            if ((ulong_it = find(facetpnt.begin(), facetpnt.end(), PrevIndex)) == facetpnt.end())  //if PrevIndex not found
            {
                //in current facet
                if (face_it != face.end())
                    ++face_it; //next face please
                else
                    break;

                continue;
            }
            else
            {
                for (unsigned int k = 0 ; k < 3; k++)
                {
                    //If current facetpnt[k] is not yet in nei_list AND it is not the CurIndex
                    if (((vec_it = find(nei.begin(), nei.end(), facetpnt[k])) == nei.end()) && facetpnt[k] != CurInd)
                    {
                        face.erase(face_it);   //erase this face
                        nei.push_back(facetpnt[k]);  //push back the index
                        PrevIndex = facetpnt[k];   //this index is now PrevIndex
                        face_it = face.begin(); //reassign the iterator
                        break;   //end this for-loop
                    }
                }
                break;  //end this while loop
            }
        }
    }
}
*/


bool SpringbackCorrection::Visit(const MeshCore::MeshFacet &rclFacet, const MeshCore::MeshFacet &rclFrom, unsigned long ulFInd, unsigned long ulLevel)
{
    if (ulLevel != m_RingCurrent)
    {
        ++m_RingCurrent;

        if (m_RingVector.size() == 0)
        {
            return false;
        }
        //else if(m_RingVector.size() < 2 && m_RingCurrent > 2)
        //{
        // m_RingVector.clear();
        // return false;
        //}
        else
        {
            m_RegionVector.push_back(m_RingVector);
            m_RingVector.clear();
        }
    }

    if (rclFacet._ulProp == 1 && rclFrom._ulProp == 1)
        m_RingVector.push_back(ulFInd);

    return true;
}