#include <features/CVertexInfoGarland.h>

//=====================================================================================================================
CVertexInfoGarland::CVertexInfoGarland( vctl::MCVertex * vertex ) :
        CVertexInfoMV( vertex, 3, 3, 3 )
{
}

//=====================================================================================================================
CVertexInfoGarland::~CVertexInfoGarland()
{
}

//=====================================================================================================================
CVertexInfoGeneratorGarland::CVertexInfoGeneratorGarland( unsigned type, unsigned source, int size, bool cont )
{
    bCont = cont;
    iType = type;
    iSource = source;
    iSize = size;
}

//=====================================================================================================================
CVertexInfoGeneratorGarland::~CVertexInfoGeneratorGarland()
{
}

//=====================================================================================================================
bool CVertexInfoGeneratorGarland::generate( CSmoothingMesh * mesh )
{
    vctl::MCVerticeS * vertices    = mesh->GetVerticeS();
    vctl::MCVertex *     vertex    = vertices->GetFirst();

    CVertexInfoGarland   * info;
    CVertexInfoGarland::clearMax();

    double  absmax = 0;

    while ( vertex )
    {
        mds::math::CDMatrix X;

        if      ( this->iType == TYPE_COVARIANCE && this->iSource == FROM_VERTICES  ) this->xMatrixVertCov( vertex, mesh, X );
        else if ( this->iType == TYPE_COVARIANCE && this->iSource == FROM_TRIANGLES ) this->xMatrixTriCov( vertex, mesh, X );
        else if ( this->iType == TYPE_GARLAND && this->iSource == FROM_VERTICES )     this->xMatrixVertGar( vertex, mesh, X );
        else this->xMatrixTriGar( vertex, mesh, X );

        // prepare transposition of the data matrix
        mds::math::CDMatrix XT( X.getNumOfCols(), 3 );
        XT.transpose( X );        

        // compute covariance matrix
        mds::math::CDMatrix     COV( 3, 3 );
        COV.mult( X, XT );

        // prepare eigenvector matrix and eigenvalue vector
        mds::math::CDMatrix     eigenvectors( 3, 3 );
        mds::math::CDVector     eigenvalues( 3 );

        // mds::math::inverse( COV );
        eigenvectors = COV;

        mds::math::eig( eigenvectors, eigenvalues );
         // trim extreme values
        if ( eigenvalues.get(0) < 0.05 ) eigenvalues.get(0) = 0.05;
        eigenvalues.get(0) = 1.0 / eigenvalues.get(0);

        if ( eigenvalues.get(1) < 0.05 ) eigenvalues.get(1) = 0.05;
        eigenvalues.get(1) = 1.0 / eigenvalues.get(1);

        if ( eigenvalues.get(2) < 0.05 ) eigenvalues.get(2) = 0.05;
        eigenvalues.get(2) = 1.0 / eigenvalues.get(2);

        if ( eigenvalues.get(0) > absmax ) absmax = eigenvalues.get(0);
        if ( eigenvalues.get(1) > absmax ) absmax = eigenvalues.get(1);
        if ( eigenvalues.get(2) > absmax ) absmax = eigenvalues.get(2);

        if ( bCont )
        {
            double nn = eigenvalues.get(0) * eigenvalues.get(0);
            nn = nn   + eigenvalues.get(1) * eigenvalues.get(1);
            nn = nn   + eigenvalues.get(2) * eigenvalues.get(2);
            nn = sqrt( nn );

            eigenvalues.get(0)  = 2 * eigenvalues.get(0)  / nn;
            eigenvalues.get(1)  = 2 * eigenvalues.get(1)  / nn;
            eigenvalues.get(2)  = 2 * eigenvalues.get(2)  / nn;
        }

        CVertexInfoGarland::tryNewMax( absmax );

        info = new CVertexInfoGarland( vertex );
        info->setMatrix( eigenvectors ),
        info->setVector( eigenvalues );

        vertex->SetValuePtr( static_cast< void* >( info ) );

        // move to the next vertex
        vertex = vertex->GetNext();
    }
	return true;
}

//=====================================================================================================================
void CVertexInfoGeneratorGarland::xMatrixTriCov( vctl::MCVertex * vertex, CSmoothingMesh * mesh, mds::math::CDMatrix & X )
{

    tTriPtrList neigh;
    mesh->getNeighbourTri( vertex, neigh, this->iSize );

    int size = neigh.size();

    vctl::MCVector3D avg_normal( 0.0, 0.0, 0.0 );

    for ( int i = 0; i < size; i++ )
    {
        vctl::MCTri * tri  = neigh[ i ];
        vctl::MCVector3D n = tri->GetNormal();

        avg_normal.SetX( avg_normal.GetX() + n.GetX() );
        avg_normal.SetY( avg_normal.GetY() + n.GetY() );
        avg_normal.SetZ( avg_normal.GetZ() + n.GetZ() );
    }

    avg_normal.SetX( avg_normal.GetX() / size );
    avg_normal.SetY( avg_normal.GetY() / size );
    avg_normal.SetZ( avg_normal.GetZ() / size );

    X.create( 3, size );
    X.fill( 0 );

    // get normal of current vertex
    vctl::MCVector3D vertex_normal = mesh->getNormal( vertex );

    // through the whole neighbourhood
    for ( int i = 0; i < size; i++ )
    {
        vctl::MCTri * tri = neigh[ i ];

        // get normal of neighbouring vertex
        vctl::MCVector3D normal =    tri->GetNormal();
        double a = tri->GetArea();

        // and add it to data vector
        X.get( 0, i )   =  normal.GetX() - avg_normal.GetX();
        X.get( 1, i )   =  normal.GetY() - avg_normal.GetY();
        X.get( 2, i )   =  normal.GetZ() - avg_normal.GetZ();
    }
}

//=====================================================================================================================
void CVertexInfoGeneratorGarland::xMatrixTriGar( vctl::MCVertex * vertex, CSmoothingMesh * mesh, mds::math::CDMatrix & X )
{


    tTriPtrList neigh;
    mesh->getNeighbourTri( vertex, neigh, this->iSize );

    int size = neigh.size();

    X.create( 3, size );
    X.fill( 0 );

    // get normal of current vertex
    vctl::MCVector3D vertex_normal = mesh->getNormal( vertex );

    // through the whole neighbourhood
    for ( int i = 0; i < size; i++ )
    {
        vctl::MCTri * tri = neigh[ i ];

        // get normal of neighbouring vertex
        vctl::MCVector3D normal =    tri->GetNormal();
        double a = tri->GetArea();

        // and add it to data vector
        X.get( 0, i )   =  normal.GetX();
        X.get( 1, i )   =  normal.GetY();
        X.get( 2, i )   =  normal.GetZ();
    }
}

//=====================================================================================================================
void CVertexInfoGeneratorGarland::xMatrixVertCov( vctl::MCVertex * vertex, CSmoothingMesh * mesh, mds::math::CDMatrix & X )
{

    tVertPtrList    neigh;
    mesh->getNeighbourVert( vertex, neigh, this->iSize );

    int size = neigh.size();

    vctl::MCVector3D avg_normal( 0.0, 0.0, 0.0 );

    for ( int i = 0; i < size; i++ )
    {
        vctl::MCVertex * vert  = neigh[ i ];
        vctl::MCVector3D n = mesh->getNormal( vert );

        avg_normal.SetX( avg_normal.GetX() + n.GetX() );
        avg_normal.SetY( avg_normal.GetY() + n.GetY() );
        avg_normal.SetZ( avg_normal.GetZ() + n.GetZ() );
    }

    avg_normal.SetX( avg_normal.GetX() / size );
    avg_normal.SetY( avg_normal.GetY() / size );
    avg_normal.SetZ( avg_normal.GetZ() / size );

    X.create( 3, size );
    X.fill( 0 );

    // through the whole neighbourhood
    for ( int i = 0; i < size; i++ )
    {
        vctl::MCVertex * vert = neigh[ i ];

        // get normal of neighbouring vertex
        vctl::MCVector3D normal =    mesh->getNormal( vert );

        // and add it to data vector
        X.get( 0, i )   =  normal.GetX() - avg_normal.GetX();
        X.get( 1, i )   =  normal.GetY() - avg_normal.GetY();
        X.get( 2, i )   =  normal.GetZ() - avg_normal.GetZ();
    }
}

//=====================================================================================================================
void CVertexInfoGeneratorGarland::xMatrixVertGar( vctl::MCVertex * vertex, CSmoothingMesh * mesh, mds::math::CDMatrix & X )
{
    tVertPtrList    neigh;
    mesh->getNeighbourVert( vertex, neigh, this->iSize );

    int size = neigh.size();

    X.create( 3, size );
    X.fill( 0 );

    // get normal of current vertex
    vctl::MCVector3D vertex_normal = mesh->getNormal( vertex );

    // through the whole neighbourhood
    for ( int i = 0; i < size; i++ )
    {
        vctl::MCVertex * vert = neigh[ i ];

        // get normal of neighbouring vertex
        vctl::MCVector3D normal =    mesh->getNormal( vert );

        // and add it to data vector
        X.get( 0, i )   =  normal.GetX();
        X.get( 1, i )   =  normal.GetY();
        X.get( 2, i )   =  normal.GetZ();
    }
}

//=====================================================================================================================
CVertexInfoGenerator * CVertexInfoGeneratorGarland::create()
{
    return new CVertexInfoGeneratorGarland();
}

//=====================================================================================================================
bool CVertexInfoGeneratorGarland::init()
{
	return true;
}

//=====================================================================================================================
void CVertexInfoGeneratorGarland::normalize()
{
}

//=====================================================================================================================
bool CVertexInfoGeneratorGarland::operator!= ( const CVertexInfoGeneratorGarland & g )
{
    if ( this->iSize    != g.iSize )    return true;
    if ( this->iType    != g.iType )    return true;
    if ( this->iSource  != g.iSource )  return true;
    return false;
}
