#ifndef MIQ_STIFFENING
#define MIQ_STIFFENING


template <class ScalarType>
ScalarType Gauss(ScalarType &value)
{
    const ScalarType E_NEPER=2.71828;
    const ScalarType den=sqrt(2.0*M_PI);
    ScalarType exponent=-0.5*pow(value,2);
    ScalarType res=((1.0/den)*pow(E_NEPER,exponent));
    return res;
}

template <class MeshType>
class StiffeningInitializer
{
    typedef typename MeshType::ScalarType ScalarType;
    typedef typename MeshType::FaceType FaceType;
    typedef typename MeshType::VertexType VertexType;
    typedef typename MeshType::CoordType CoordType;

public:

    static void colorByStiffening(MeshType & mesh,
                                  typename MeshType::ScalarType MaxVal=16)
    {
        bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
        assert(hasStiffness);
        typename MeshType::template PerFaceAttributeHandle<float> Handle_Stiffness=vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<float>(mesh,std::string("Stiffness"));

        for (unsigned int i=0;i<mesh.face.size();i++)
        {
            //CMesh::ScalarType val=MaxVal-mesh.face[i].stiffening+1;
            ScalarType val=MaxVal-Handle_Stiffness[i]+1;
            if (val<1)val=1;
            mesh.face[i].C()=vcg::Color4b::ColorRamp(1.0,MaxVal,val);
        }
    }

    static void AddGaussStiffening(MeshType & mesh,ScalarType C)
    {
        int radius=floor(C);
        if (C<4)radius=4;
        typename MeshType::template PerFaceAttributeHandle<float> Handle_Stiffness;

        bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
        if(!hasStiffness)
            Handle_Stiffness=vcg::tri::Allocator<MeshType>::template AddPerFaceAttribute<float>(mesh,std::string("Stiffness"));
        else
            Handle_Stiffness=vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<float>(mesh,std::string("Stiffness"));

        bool hasSingular = vcg::tri::HasPerVertexAttribute(mesh,std::string("Singular"));
        assert(hasSingular);

        typename MeshType::template PerVertexAttributeHandle<bool> Handle_Singular;
        Handle_Singular=vcg::tri::Allocator<CMesh>::GetPerVertexAttribute<bool>(mesh,std::string("Singular"));

        std::vector<VertexType*> to_stiff;
        for(unsigned int i=0;i<mesh.vert.size();i++)
        {
            VertexType *v=&mesh.vert[i];
            if (v->IsD())continue;
            //if (!v->IsSingular())continue;
            if (!Handle_Singular[v])continue;
            to_stiff.push_back(v);
        }
        for(unsigned int i=0;i<mesh.face.size();i++)
        {

            FaceType *f=&(mesh.face[i]);
            if (f->IsD())continue;
            if (!f->IsV())continue;
            to_stiff.push_back(f->V(0));
            to_stiff.push_back(f->V(1));
            to_stiff.push_back(f->V(2));
        }
        std::sort(to_stiff.begin(),to_stiff.end());
        typename std::vector<VertexType*>::iterator new_end;
        new_end=std::unique(to_stiff.begin(),to_stiff.end());
        int dist=distance(to_stiff.begin(),new_end);
        to_stiff.resize(dist);
        for (unsigned int i=0;i<to_stiff.size();i++)
        {
            VertexType *v=to_stiff[i];
            for (int r=0;r<radius;r++)
            {
                CMesh::ScalarType stiffVal=((CMesh::ScalarType)r)/(CMesh::ScalarType)radius;//((ScalarType)(radius-r))/(ScalarType)radius;
                stiffVal*=3.0;
                stiffVal=Gauss(stiffVal)/0.4;
                stiffVal=1+(stiffVal*C);
                std::vector<FaceType*> ring;
                //mesh.GetConnectedFaces(v,r,ring);
                VFExtendedStarVF(v,r,ring);
                ///then set stiffening
                for (unsigned int k=0;k<ring.size();k++)
                {
                    FaceType* f=ring[k];
                    //if (f->stiffening<stiffVal)
                    //    f->stiffening=stiffVal;
                    if (Handle_Stiffness[f]<stiffVal)
                        Handle_Stiffness[f]=stiffVal;
                }
            }
        }
    }

    static bool updateStiffeningJacobianDistorsion(MeshType & mesh,ScalarType grad_size)
    {

        typename MeshType::template PerFaceAttributeHandle<float> Handle_Stiffness;

        bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
        if(!hasStiffness)
            Handle_Stiffness=vcg::tri::Allocator<MeshType>::template AddPerFaceAttribute<float>(mesh,std::string("Stiffness"));
        else
            Handle_Stiffness=vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<float>(mesh,std::string("Stiffness"));

        bool flipped = NumFlips(mesh)>0;
        //if (h == 0.0)
        //    return flipped;
        //
        //assert(h != 0.0);

        if (!flipped)
            return false;
       ScalarType maxL=0;
       ScalarType maxD=0;
        if (flipped)
        {
            const double c = 1.0;
            const double d = 5.0;

            for (unsigned int i = 0; i < mesh.face.size(); ++i)
            {
                ScalarType dist=Distortion(mesh.face[i],grad_size);
                if (dist>maxD)maxD=dist;
                ScalarType absLap=fabs(LaplaceDistortion(mesh.face[i], grad_size));
                if (absLap>maxL)maxL=absLap;
                ScalarType stiffDelta = std::min(c * absLap, d);
                //mesh.face[i].stiffening+=stiffDelta;
                Handle_Stiffness[i]+=stiffDelta;
            }
        }
        printf("Maximum Distorsion %4.4f \n",maxD);
        printf("Maximum Laplacian %4.4f \n",maxL);
        return flipped;
    }

    static void InitDefaultStiffening(MeshType & mesh)
    {
        typename MeshType::template PerFaceAttributeHandle<float> Handle_Stiffness;

        bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
        if(!hasStiffness)
            Handle_Stiffness=vcg::tri::Allocator<MeshType>::template AddPerFaceAttribute<float>(mesh,std::string("Stiffness"));

        for(unsigned int i=0;i<mesh.face.size();i++)
        {
            FaceType *f=&(mesh.face[i]);
            //f->stiffening=1;
            Handle_Stiffness[f]=1;
        }
    }

};

#endif