#include "sqrt3subdiv.h"
#include <vcg/complex/trimesh/allocate.h>
#include <vcg/complex/trimesh/clean.h>
#include <QtScript>
#include <vector>
#include <stdio.h>
#include <math.h>
#include "vector.h"

using namespace vcg;
using namespace tri;
using namespace std;
// Constructor usually performs only two simple tasks of filling the two lists 
//  - typeList: with all the possible id of the filtering actions
//  - actionList with the corresponding actions. If you want to add icons to your filtering actions you can do here by construction the QActions accordingly

Sqrt3Subdiv::Sqrt3Subdiv()
{
    typeList << SQRT3_SUBDIV;
    foreach(FilterIDType tt , types()) {
        actionList << new QAction(filterName(tt), this);
    }
}


// ST() must return the very short string describing each filtering action 
// (this string is used also to define the menu entry)
QString Sqrt3Subdiv::filterName(FilterIDType filterId) const
{
    switch(filterId) {
        case SQRT3_SUBDIV :
                return QString("Sqrt(3) subdivision");
        default :
                assert(0);
    }
    return QString();
}


// Info() must return the longer string describing each filtering action 
// (this string is used in the About plugin dialog)
 QString Sqrt3Subdiv::filterInfo(FilterIDType filterId) const
{
    switch(filterId) {
        case SQRT3_SUBDIV :
        return QString("bbertka: Performs one iteration of the sqrt(3) subdivision algorithm.");
        default :
                assert(0);
    }
    return QString("bbertka: Unknown Filter");
}


// The FilterClass describes in which generic class of filters it fits. 
// This choice affect the submenu in which each filter will be placed 
// More than a single class can be choosen.
Sqrt3Subdiv::FilterClass Sqrt3Subdiv::getClass(QAction *a)
{
    switch(ID(a)) {
        case SQRT3_SUBDIV :
                return MeshFilterInterface::Smoothing;
        default :
                assert(0);
    }
    return MeshFilterInterface::Generic;
}


/*
  No parameters
  */
void Sqrt3Subdiv::initParameterSet(QAction *action,MeshModel &m, RichParameterSet & parlst)
{
    switch( ID(action) ) {
        case SQRT3_SUBDIV :
             break;
	}
}


// The Real Core Function doing the actual mesh processing.
// Subdivision sqrt(3) method
bool Sqrt3Subdiv::applyFilter( QAction */*filter*/, MeshDocument &md,
                               RichParameterSet & par, vcg::CallBackPos *cb )
{
    MeshModel &m=*md.mm();

//=========================================================================
// Step 0: Save the neighborhoods of the original vertices
//=========================================================================
//this is for saving neighbors on the original mesh
    vector<Point3f> vertexNeighbors;
    vector<vector<Point3f> > neighbors;  // stores each neighboring vertex

    m.cm.face.EnableFFAdjacency();
    UpdateTopology<CMeshO>::FaceFace(m.cm);
    m.cm.vert.EnableVFAdjacency();
    m.cm.face.EnableVFAdjacency();
    UpdateTopology<CMeshO>::VertexFace(m.cm);
    for(unsigned int vi = 0; vi < m.cm.vert.size(); vi++){
        CVertexO *vI = &m.cm.vert[vi];  // vertex pointer to start vertex
        face::Pos<CFaceO> pos(vI->VFp(), vI);   //pos for current vertex
        CVertexO *firstV = pos.VFlip();
        vertexNeighbors.clear();
        vertexNeighbors.push_back(firstV->P());
        do{
            pos.NextE();
            vI = pos.VFlip();
            vertexNeighbors.push_back( vI->P() );  // save the neighboring point
        }while(vI != firstV);
            neighbors.push_back( vertexNeighbors );
     }
    // unlock
    m.cm.face.DisableFFAdjacency();
    m.cm.vert.DisableVFAdjacency();
    m.cm.face.DisableVFAdjacency();


//=========================================================================
// Step 1: Divide each face into three, about a point in the middle
//=========================================================================
    vector<Point3f> facePoint;
    vector<vector<Point3f> > startFacePoints;
    vector<CVertexO*> faceVertex_ref;
    vector<vector<CVertexO*> > subdivFaceVertex_ref;
    // P() references the point, and we can do math with it
    // V() references a pointer to the point, we can reuse it, for other face

    CMeshO::FaceIterator faceitor;
    unsigned len = m.cm.face.size();        //number of original faces
    for(unsigned int i = 0; i < len; i++){
        cb( 100 * i / len, "Computing initial split...");
        if(m.cm.face[i].IsD())  //skip any face which may be hanging around
            continue;

        //compute a center vertex on the face, and add to master vertex list
        Allocator<CMeshO>::AddVertices(m.cm, 1);
        m.cm.vert[m.cm.vert.size() - 1].P() = (m.cm.face[i].P(0) +
                                               m.cm.face[i].P(1) +
                                               m.cm.face[i].P(2))/ 3.0;
        //make new faces forming the initial split
        faceitor =  Allocator<CMeshO>::AddFaces(m.cm, 3);
        faceitor->V(0) = m.cm.face[i].V(0);
        faceitor->V(1) = m.cm.face[i].V(1);
        faceitor->V(2) = &m.cm.vert[m.cm.vert.size() - 1];
        faceitor++;
        faceitor->V(0) = m.cm.face[i].V(2);
        faceitor->V(1) = m.cm.face[i].V(0);
        faceitor->V(2) = &m.cm.vert[m.cm.vert.size() - 1];
        faceitor++;;
        faceitor->V(0) = m.cm.face[i].V(1);
        faceitor->V(1) = m.cm.face[i].V(2);
        faceitor->V(2) = &m.cm.vert[m.cm.vert.size() - 1];

        //remove original face, and save the original face points
        facePoint.clear();
        facePoint.push_back(m.cm.face[i].P(0));
        facePoint.push_back(m.cm.face[i].P(1));
        facePoint.push_back(m.cm.face[i].P(2));
        startFacePoints.push_back(facePoint);
        Allocator<CMeshO>::DeleteFace(m.cm,m.cm.face[i]);
    }

//=========================================================================
// Step 2: Find the spoke edge for each vertex, and check if in original
//=========================================================================

    // Set up adjacency lists.  NOTE: there is a locking mechanism which makes
    // it impossible to add new verticecs or edges directly on mesh (crashes),
    // so I save the face vertex pointers for processing later
    //lock the mesh for adjacency and one ring traversal operations
    m.cm.face.EnableFFAdjacency();
    UpdateTopology<CMeshO>::FaceFace(m.cm);
    m.cm.vert.EnableVFAdjacency();
    m.cm.face.EnableVFAdjacency();
    UpdateTopology<CMeshO>::VertexFace(m.cm);

    subdivFaceVertex_ref.clear();       // reset subdiv face vertex pointers

    // tests each spoke eminating from vertex to see if also in original faces
    // when an edge is found, a new face is made
    for(unsigned int vi = 0; vi < m.cm.vert.size(); vi++){
        Point3f v0, v1,v2;              // actual points on current face
        CVertexO *vI = &m.cm.vert[vi];  // vertex pointer to start vertex
        CVertexO *vJ, *vA, *vB;         //  vertex pointer to opposite vertex
                                        // vI and vJ spokes, opposite v ptr
                                        // vA is first face diagonal v pointer
                                        // vB is second face diagonal v pointer
        face::Pos<CFaceO> pos(vI->VFp(), vI);   //pos for current vertex
        face::Pos<CFaceO> A_pos, B_pos;         //pos for adjacent faces

        CVertexO *firstV = pos.VFlip(); // check each spoke around the vertex until the
                                        // opposite spoke's vertex is the original vertex
        do{
            v0 = vI->P();               // spoke starting vertex (vI starting vertex pointer)
            pos.NextE();
            vJ = pos.VFlip();           // vertex pointer to opposite spoke verte

            A_pos = pos;                // pos for first adjacency
            A_pos.FlipE();
            vA = A_pos.VFlip();         // first diagonal vertex pointer
            v1 = vA->P();               // first diagonal vertex

            B_pos = pos;                // pos for second adjacency
            B_pos.FlipF();
            B_pos.FlipE();
            vB = B_pos.VFlip();         // second diagonal vertex pointer
            v2 = vB->P();               // second diagonal vertex

            // check to see if this spoke forms an an original edge by comparing the
            // spoke pointers to locations of original face edges
            if(hasEdgePointers(vI, vJ, startFacePoints) >= 0){                
                // Add new sqrt(3) subdivision face vertices
                faceVertex_ref.clear();
                faceVertex_ref.push_back(vA);   // right diagonal face vertex pointer
                faceVertex_ref.push_back(vB);   // left diagonal face vertex pointer
                faceVertex_ref.push_back(vI);   // starting diagonal face vertex pointer
                                                // Add new sqrt(3) subdivision face
                subdivFaceVertex_ref.push_back(faceVertex_ref);
            }
            cb( 100 * vi / m.cm.vert.size(), "Computing edge flip on vertices...");
        }while(vJ != firstV);   // check each spoke around the vertex until the
                                // opposite spoke's vertex is the original vertex
    }

    // unlock
    m.cm.face.DisableFFAdjacency();
    m.cm.vert.DisableVFAdjacency();
    m.cm.face.DisableVFAdjacency();

    //Delete the current faces which form the initital three-split
    int newLen = m.cm.face.size();
    for(int i = len; i < newLen; i++){
        cb( 100 * i / (newLen - len), "Removing old faces...");
        Allocator<CMeshO>::DeleteFace(m.cm,m.cm.face[i]);
    }
    // Draw the sqrt(3) subdivision faces from the pos operation re-using vertices from,
    // re-using vertex pointers from the first split
    for(unsigned int i = 0; i < subdivFaceVertex_ref.size(); i++){
        cb( 100 * i / subdivFaceVertex_ref.size(), "Initializing sqrt(3) subdivision faces...");
        faceitor =  Allocator<CMeshO>::AddFaces(m.cm, 1);
        faceitor->V(0) = subdivFaceVertex_ref[i][0];
        faceitor->V(1) = subdivFaceVertex_ref[i][1];
        faceitor->V(2) = subdivFaceVertex_ref[i][2];
    }
//=========================================================================
// Step 3: Relaxation, using Step 0 neighbors
//=========================================================================
    //only relax the old vertices, found at beginning of vertex list (3 * num faces)
    for( unsigned int i = 0; i < neighbors.size(); i++) {  //ith neighborhood for ith vertex
        // sum around the vertex neighborhood
        Vector3 sum;
        Vector3 pv;
        sum.Set(0,0,0);
        float an = 0.0;
        Point3f p = m.cm.vert[i].P();
        pv.Set(p.X(), p.Y(), p.Z());
        int valence = int(neighbors[i].size());
        vertexNeighbors = neighbors[i];
        for(int k = 0; k < valence; k++){
            Vector3 tmp;
            tmp.Set(vertexNeighbors[k].X(),
                    vertexNeighbors[k].Y(),
                    vertexNeighbors[k].Z() );
            sum.Add(sum, tmp);
        }
        an = (4.0 - 2.0*cos( (2.0*M_PI)/float(valence) ) ) / 9.0;
        float k = float(1)/float(valence);
        sum.Scale(an*k, sum);
        pv.Scale( (1.0-an), pv);
        pv.Add(pv, sum);
        p.X() = pv.x;
        p.Y() = pv.y;
        p.Z() = pv.z;
        m.cm.vert[i].P() = p;
    }


    // some clean up just in case
    Clean<CMeshO>::RemoveDuplicateVertex(m.cm);
    Clean<CMeshO>::RemoveDuplicateFace(m.cm);
    UpdateFlags<CMeshO>::VertexClear(m.cm);
    UpdateFlags<CMeshO>::FaceClear(m.cm);
    Clean<CMeshO>::RemoveUnreferencedVertex( m.cm );

    // set the normals (this is the better than per-face of the two)
    UpdateNormals<CMeshO>::PerVertexNormalizedPerFace(m.cm);

    UpdateBounding<CMeshO>::Box(m.cm);
    return true;
}

/*
  Returns the index in faces which the edge vertices are found in
  Returns -1 on fail
  */
int Sqrt3Subdiv::hasEdgePointers(CVertexO *vI, CVertexO *vJ,
                                 vector<vector<Point3f> > startfaceList){
    bool foundI = false;
    bool foundJ = false;
    for(unsigned int i = 0; i < startfaceList.size(); i++){
        foundI = false;
        foundJ = false;
        if(startfaceList[i][0] == vI->P() ||
           startfaceList[i][1] == vI->P() ||
           startfaceList[i][2] == vI->P())
            foundI = true;

        if(startfaceList[i][0] == vJ->P() ||
           startfaceList[i][1] == vJ->P() ||
           startfaceList[i][2] == vJ->P())
            foundJ = true;

        if(foundI && foundJ)
            return i;
    }
    return -1;
}


QString Sqrt3Subdiv::filterScriptFunctionName( FilterIDType filterID )
{
    switch(filterID) {
        case SQRT3_SUBDIV :
                return QString("squareRootThreeSubdivision");
        default :
                assert(0);
	}
	return QString();
}



Q_EXPORT_PLUGIN(Sqrt3Subdiv)