Subversion Repositories gelsvn

Rev

Rev 493 | Rev 499 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | RSS feed

/*
*  MeshEdit is a small application which allows you to load and edit a mesh.
*  The mesh will be stored in GEL's half edge based Manifold data structure.
*  A number of editing operations are supported. Most of these are accessible from the 
*  console that pops up when you hit 'esc'.
*
*  Created by J. Andreas Bærentzen on 15/08/08.
*  Copyright 2008 __MyCompanyName__. All rights reserved.
*
*/

#include <string>
#include <iostream>
#include <vector>

#include <GL/glew.h>

#include <GLConsole/GLConsole.h>


#include <CGLA/eigensolution.h>
#include <CGLA/Vec2d.h>
#include <CGLA/Vec3d.h>
#include <CGLA/Mat3x3d.h>
#include <CGLA/Mat2x2d.h>
#include <CGLA/Mat2x3d.h>

#include <LinAlg/Matrix.h>
#include <LinAlg/Vector.h>
#include <LinAlg/LapackFunc.h>

#include <GLGraphics/gel_glut.h>

#include <HMesh/Manifold.h>
#include <HMesh/HalfEdgeHandle.h>
#include <HMesh/VertexHandle.h>
#include <HMesh/FaceHandle.h>
#include <HMesh/VertexCirculator.h>
#include <HMesh/FaceCirculator.h>
#include <HMesh/mesh_optimization.h>
#include <HMesh/triangulate.h>
#include <HMesh/load.h>
#include <HMesh/quadric_simplify.h>
#include <HMesh/smooth.h>
#include <HMesh/x3d_save.h>
#include <HMesh/obj_save.h>
#include <HMesh/off_save.h>
#include <HMesh/mesh_optimization.h>
#include <HMesh/triangulate.h>
#include <HMesh/close_holes.h>
#include <HMesh/caps_and_needles.h>
#include <HMesh/refine_edges.h>
#include <HMesh/subdivision.h>

#include <Util/Timer.h>
#include <Util/ArgExtracter.h>

#include "harmonics.h"
#include "curvature.h"
#include "Renderer.h"

#include "VisObj.h"

using namespace std;
using namespace HMesh;
using namespace Geometry;
using namespace GLGraphics;
using namespace CGLA;
using namespace Util;
using namespace LinAlg;
using namespace CVarUtils;

inline VisObj& get_vis_obj(int i)
{
    static VisObj vo[9];
    return vo[i];
}

inline VisObj& avo()
{
    static int& active = CreateCVar("active_mesh",0);
    return get_vis_obj(active);
}

inline Manifold& active_mesh()
{
    return avo().mesh();
}

inline GLViewController& active_view_control()
{
    return avo().view_control();
}

// Single global instance so glut can get access
Trie CVarTrie;
GLConsole theConsole;


////////////////////////////////////////////////////////////////////////////////
bool MyConsoleHelp(std::vector<std::string> *args)
{
    theConsole.Printf("");
    theConsole.Printf("----------------- HELP -----------------");
    theConsole.Printf("Press ESC key to open and close console");
    theConsole.Printf("Press TAB to see the available commands and functions");
    theConsole.Printf("Functions are shown in green and variables in yellow");
    theConsole.Printf("Setting a value: [command] = value");
    theConsole.Printf("Getting a value: [command]");
    theConsole.Printf("Functions: [function] [arg1] [arg2] ...");
    theConsole.Printf("Entering arg1=? or arg1=help will give a description.");
    theConsole.Printf("History: Up and Down arrow keys move through history.");
    theConsole.Printf("Tab Completion: TAB does tab completion and makes suggestions.");
    theConsole.Printf("");
    theConsole.Printf("Keyboard commands (when console is not active):");
    theConsole.Printf("w   : switch to display.render_mode = wireframe");
    theConsole.Printf("i   : switch to display.render_mode = isophotes");
    theConsole.Printf("r   : switch to display.render_mode = reflection");
    theConsole.Printf("m   : switch to display.render_mode = metallic");
    theConsole.Printf("g   : switch to display.render_mode = glazed");
    theConsole.Printf("n   : switch to display.render_mode = normal");
    theConsole.Printf("h   : switch to display.render_mode = harmonics");
    theConsole.Printf("f   : toggle smooth/flat shading");
    theConsole.Printf("1-9 : switch between active meshes.");
    theConsole.Printf("d   : (display.render_mode = harmonics) diffuse light on and off");
    theConsole.Printf("h   : (display.render_mode = harmonics) highlight on and off ");
    theConsole.Printf("+/- : (display.render_mode = harmonics) which eigenvector to show");
    theConsole.Printf("q   : quit program");
    theConsole.Printf("ESC : open console");
    theConsole.Printf("");
    theConsole.Printf("Mouse: Left button rotates, middle zooms, right pans");
    theConsole.Printf("----------------- HELP -----------------");
    theConsole.Printf("");
    return "";
}

bool wantshelp(std::vector<std::string> *args)
{
    if(args->size() == 0) return false;
    string str = (*args)[0];
    if(str=="help" || str=="HELP" || str=="Help" || str=="?") return true;
    return false;
}

/// Function that aligns two meshes.
bool console_align(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: align <dest> <src>");
        theConsole.Printf("This function aligns dest mesh with src");
        theConsole.Printf("In practice the GLViewController of src is copied to dst.");
        theConsole.Printf("both arguments are mandatory and must be numbers between 1 and 9.");
        theConsole.Printf("Note that results might be unexpexted if the meshes are not on the same scale");
        return "";
    }

    int dest = 0;
    if(args->size()>0)
    {
        istringstream a0((*args)[0]);
        a0 >> dest;
        --dest;
        if(dest <0 || dest>8) return "dest mesh out of range (1-9)";
    }
    else return "neither source nor destination mesh?!";
    int src = 0;
    if(args->size()>1)
    {
        istringstream a1((*args)[1]);
        a1 >> src;
        --src;
        if(src <0 || src>8) return "src mesh out of range (1-9)";
    }   
    else return "no src mesh?";

    get_vis_obj(dest).view_control() = get_vis_obj(src).view_control();

    return "";
}

bool console_flatten(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: flatten <floater|harmonic|barycentric>");
        theConsole.Printf("This function flattens a meshs with a simple boundary. It is mostly for showing mesh");
        theConsole.Printf("parametrization methods. The current mesh MUST have a SINGLE boundary loop");
        theConsole.Printf("This loop is mapped to the unit circle in a regular fashion (equal angle intervals).");
        theConsole.Printf("All non boundary vertices are placed at the origin. Then the system is relaxed iteratively");
        theConsole.Printf("using the weight scheme given as argument.");
        return "";
    }
    enum WeightScheme {FLOATER_W, HARMONIC_W, BARYCENTRIC_W};
    WeightScheme ws = BARYCENTRIC_W;
    if(args->size()>0)
    {
        if((*args)[0] == "floater")
            ws = FLOATER_W;
        else if((*args)[0] == "harmonic")
            ws = HARMONIC_W;
    }
    else
        return "";

    active_mesh().enumerate_vertices();
    active_mesh().enumerate_halfedges();

    vector<float> edge_weights(active_mesh().no_halfedges());
    for(HalfEdgeHandle h = active_mesh().halfedges_begin(); h != active_mesh().halfedges_end(); ++h)
        if(!is_boundary(h)){
            Vec3f p0 = h.vert().get().pos;
            Vec3f p1 = h.next().vert().get().pos;
            Vec3f p2 = h.opp().vert().get().pos;
            Vec3f p3 = h.opp().next().vert().get().pos;

            if(ws == FLOATER_W){
                float d = acos(min(1.0f, max(-1.0f, dot(normalize(p2-p0), normalize(p3-p0)))));
                float g = acos(min(1.0f, max(-1.0f, dot(normalize(p2-p0), normalize(p1-p0)))));
                edge_weights[h.opp().get().touched]  = (tan(d/2) + tan(g/2)) / (p0-p2).length();


                d = acos(min(1.0f, max(-1.0f, dot(normalize(p0-p2), normalize(p1-p2)))));
                g = acos(min(1.0f, max(-1.0f, dot(normalize(p0-p2), normalize(p3-p2)))));
                edge_weights[h.get().touched]  = (tan(d/2) + tan(g/2)) / (p0-p2).length();
            }
            else if(ws == HARMONIC_W)
            {

                float a = acos(min(1.0f, max(-1.0f, dot(normalize(p0-p3), normalize(p2-p3)))));
                float b = acos(min(1.0f, max(-1.0f, dot(normalize(p2-p1), normalize(p0-p1)))));
                float w=0;
                if(a+b < M_PI)
                    w = sin(a+b)/(sin(a)+sin(b));
                edge_weights[h.get().touched]  = w;
                edge_weights[h.opp().get().touched]  = w;
            }
            else
            {
                edge_weights[h.get().touched]  = valency(h.opp().vert());
                edge_weights[h.opp().get().touched]  = valency(h.vert());
            }   
        }

        ofstream ofs("parametrized.obj");

        ofs << "mtllib parametrized.mtl\nusemtl mat\n" << endl;

        VertexHandle v = active_mesh().vertices_begin();
        for(; v != active_mesh().vertices_end(); ++v)
            ofs << "v " << v.get().pos[0] << " " << v.get().pos[1] << " " << v.get().pos[2] << endl;
        ofs << endl;

        v = active_mesh().vertices_begin();
        for(; v != active_mesh().vertices_end(); ++v){
            if(is_boundary(v))
                break;
        }
        int n = 0;
        VertexHandle bv = v;
        do{
            ++n;
            bv = bv.out().vert();
        }
        while(bv != v);

        int i = 0;
        do{
            float a = 2.0*M_PI*float(i)/n;
            bv.get().pos = Vec3f(cos(a), sin(a), 0);
            ++i;
            bv = bv.out().vert();
        }
        while(bv != v);

        for(v = active_mesh().vertices_begin(); v != active_mesh().vertices_end(); ++v)
            if(!is_boundary(v))
                v.get().pos = Vec3f(0.0);

        for(int i = 0; i < 10000; ++i){
            for(v = active_mesh().vertices_begin(); v != active_mesh().vertices_end(); ++v){
                if(!is_boundary(v)){
                    Vec3f p_new(0);
                    float w_sum = 0;
                    for(VertexCirculator vc(v); !vc.end(); ++vc){
                        float w = edge_weights[vc.halfedge().get().touched];
                        p_new += vc.vertex().get().pos * w;
                        w_sum += w;
                    }
                    v.get().pos = p_new/w_sum;
                }
            }
        }

        for(v = active_mesh().vertices_begin(); v != active_mesh().vertices_end(); ++v)
            ofs << "vt " << (0.5*v.get().pos[0]+0.5) << " " << (0.5*v.get().pos[1]+0.5)  << endl;

        ofs << endl;

        for(FaceHandle f = active_mesh().faces_begin(); f != active_mesh().faces_end(); ++f)
        {
            ofs << "f ";
            for(FaceCirculator fc(f); !fc.end(); ++fc){
                int idx = fc.vertex().get().touched + 1;
                ofs << idx << "/" << idx <<" ";
            }
            ofs << endl;
        }

        return "";
}

bool console_save(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: save <name.x3d|name.obj> ");
        return "";
    }
    string& file_name = (*args)[0];
    if(args->size() == 1)
    {
        if(file_name.substr(file_name.length()-4,file_name.length())==".obj")
        {
            obj_save(file_name, active_mesh());
            return "";
        }
        else if(file_name.substr(file_name.length()-4,file_name.length())==".off")
        {
        off_save(file_name, active_mesh());
        return "";
        }
        else if(file_name.substr(file_name.length()-4,file_name.length())==".x3d")
        {
            x3d_save(file_name, active_mesh());
            return "";
        }
        return "unknown format";
    }
    return "usage: save <name.x3d|name.obj> ";
}

bool MyConsoleLoad( std::vector<std::string> *vArgs )
{
    std::string sFile = "cvars.xml";
    std::vector< std::string > vAcceptedSubstrings;

    if( vArgs->size() > 0 ) {
        sFile = (*vArgs)[0]; 
        for( size_t i = 1; i < vArgs->size(); i++ ) {
            vAcceptedSubstrings.push_back( (*vArgs)[i] );
        }
    }
    theConsole.Printf("Loading file from \"%s\".", sFile.c_str() );
    if( !CVarUtils::Load( sFile, vAcceptedSubstrings) ) {
        theConsole.Printf( "Error loading file.\n" );
    }
    return "";
}


bool console_refine_edges(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: refine.split_edges <length>");
        theConsole.Printf("splits edges longer than <length>; default is 0.5 times average length");
        return "";
    }

    float thresh = 0.5f;
    if(args->size() > 0)
    {
        istringstream a0((*args)[0]);
        a0 >> thresh;
    }
    float avg_length = average_edge_length(active_mesh());
    refine_edges(active_mesh(), thresh * avg_length);
    return "";

}

bool console_refine_faces(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: refine.split_faces ");
        theConsole.Printf("usage:  Takes no arguments. Inserts a vertex at the centre of each face.");
        return "";
    }

    triangulate_by_vertex_face_split(active_mesh());
    return "";

}

bool console_cc_subdivide(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: refine.catmull_clark ");
        theConsole.Printf("Splits each polygon into four (Catmull Clark style)");
        return "";
    }
    cc_split(active_mesh(),active_mesh());
    return "";
}

bool console_dual(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: dual ");
        theConsole.Printf("Produces the dual by converting each face to a vertex placed at the barycenter.");
        return "";
    }

    Manifold& m = active_mesh();

    // make sure every face knows its number
    m.enumerate_faces();

    vector<Vec3f> vertices(m.no_faces());
    vector<unsigned int> faces;
    vector<unsigned int> indices;

    // Create new vertices. Each face becomes a vertex whose position
    // is the centre of the face
    unsigned int i=0;
    for(FaceHandle f = m.faces_begin(); f != m.faces_end(); ++f,++i)
        vertices[i] = centre(f);

    // Create new faces. Each vertex is a new face with N=valency of vertex
    // edges.
    for(VertexHandle v = m.vertices_begin(); v!= m.vertices_end(); ++v){
        if(!is_boundary(v)){
            VertexCirculator vc(v);
            vector<size_t> index_tmp;
            for(; !vc.end(); ++vc){
                index_tmp.push_back(vc.face().get().touched);
            }

            // Push vertex indices for this face onto indices vector.
            // The circulator moves around the face in a clockwise fashion
            // so we just reverse the ordering.
            indices.insert(indices.end(), index_tmp.rbegin(), index_tmp.rend());

            // Insert face valency in the face vector.
            faces.push_back(vc.no_steps());
        }
    }

    // Clear the manifold before new geometry is inserted.
    m.clear();

    // And build
    m.build(static_cast<unsigned int>(vertices.size()), 
                        reinterpret_cast<float*>(&vertices[0]), 
                        static_cast<unsigned int>(faces.size()), 
                        static_cast<const unsigned int*>(&faces[0]),
                        static_cast<const unsigned int*>(&indices[0]));

    return "";
}


bool console_minimize_curvature(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: optimize.minimize_curvature <anneal>");
        theConsole.Printf("Flip edges to minimize mean curvature.");
        theConsole.Printf("If anneal is true, simulated annealing (slow) is used rather than a greedy scheme");
        return "";
    }
    bool anneal=false;
    if(args->size() > 0)
    {
        istringstream a0((*args)[0]);
        a0 >> anneal;
    }

    minimize_curvature(active_mesh(), anneal);
    avo().post_create_display_list();
    return "";
}

bool console_minimize_dihedral(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: optimize.minimize_dihedral <iter> <anneal> <use_alpha> <gamma> ");
        theConsole.Printf("Flip edges to minimize dihedral angles.");
        theConsole.Printf("Iter is the max number of iterations. anneal tells us whether to use ");
        theConsole.Printf("simulated annealing and not greedy optimization. use_alpha (default=true) ");
        theConsole.Printf("means to use angle and not cosine of anglegamma (default=4) is the power ");
        theConsole.Printf("to which we raise the dihedral angle");
        return "";
    }
    int iter = 1000;
    if(args->size() > 0)
    {
        istringstream a0((*args)[0]);
        a0 >> iter;
    }

    bool anneal = false;
    if(args->size() > 1)
    {
        istringstream a0((*args)[0]);
        a0 >> anneal;
    }

    bool use_alpha = true;
    if(args->size() > 2)
    {
        istringstream a0((*args)[0]);
        a0 >> use_alpha;
    }

    float gamma = 4.0f;
    if(args->size() > 3)
    {
        istringstream a0((*args)[0]);
        a0 >> gamma;
    }


    minimize_dihedral_angle(active_mesh(), iter, anneal, use_alpha, gamma);
    return "";
}

bool console_maximize_min_angle(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: optimize.maximize_min_angle <thresh> <anneal>");
        theConsole.Printf("Flip edges to maximize min angle - to make mesh more Delaunay.");
        theConsole.Printf("If the dot product of the normals between adjacent faces < thresh");
        theConsole.Printf("no flip will be made. anneal selects simulated annealing rather ");
        theConsole.Printf("nthan greedy optimization.");
        return "";
    }
    float thresh = 0.0f;
    if(args->size() > 0)
    {
        istringstream a0((*args)[0]);
        a0 >> thresh;
    }
    bool anneal = false;
    if(args->size() > 1)
    {
        istringstream a0((*args)[0]);
        a0 >> anneal;
    }
    maximize_min_angle(active_mesh(),thresh,anneal);
    return "";
}


bool console_optimize_valency(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: optimize.valency <anneal> ");
        theConsole.Printf("Optimizes valency for triangle meshes. Anneal selects simulated annealing rather than greedy optim.");
        return "";
    }
    bool anneal = false;
    if(args->size() > 0)
    {
        istringstream a0((*args)[0]);
        a0 >> anneal;
    }
    optimize_valency(active_mesh(), anneal);
    return "";
}

bool console_analyze(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage:  harmonics.analyze");
        theConsole.Printf("Creates the Laplace Beltrami operator for the mesh and finds all eigensolutions.");
        theConsole.Printf("It also projects the vertices onto the eigenvectors - thus transforming the mesh");
        theConsole.Printf("to this basis.");
        theConsole.Printf("Note that this will stall the computer for a large mesh - as long as we use Lapack.");
        return "";
    }
    avo().harmonics_analyze_mesh();
    return "";
}


bool console_partial_reconstruct(std::vector<std::string> *args)
{
    if(args->size() != 3)
        theConsole.Printf("usage: haramonics.partial_reconstruct <e0> <e1> <s>");

    if(wantshelp(args)) {
        theConsole.Printf("Reconstruct from projections onto eigenvectors. The two first arguments indicate");
        theConsole.Printf("the eigenvector interval that we reconstruct from. The last argument is the ");
        theConsole.Printf("scaling factor. Thus, for a vertex, v, the formula for computing the position, p, is:");
        theConsole.Printf("for (i=e0; i<=e1;++i) p += proj[i] * Q[i][v] * s;");
        theConsole.Printf("where proj[i] is the 3D vector containing the x, y, and z projections of the mesh onto");
        theConsole.Printf("eigenvector i. Q[i][v] is the v'th coordinate of the i'th eigenvector.");
        theConsole.Printf("Note that if vertex coordinates are not first reset, the result is probably unexpected.");
    }

    if(args->size() != 3)
        return "";

    int E0,E1;
    float scale;
    istringstream a0((*args)[0]);
    a0 >> E0;
    istringstream a1((*args)[1]);
    a1 >> E1;
    istringstream a2((*args)[2]);
    a2 >> scale;
    avo().harmonics_partial_reconstruct(E0,E1,scale);
    return "";
}

bool console_reset_shape(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: harmonics.reset_shape ");
        theConsole.Printf("Simply sets all vertices to 0,0,0. Call this before doing partial_reconstruct");
        theConsole.Printf("unless you know what you are doing.");
        return "";
    }
    avo().harmonics_reset_shape();
    return "";
}


bool console_close_holes(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: cleanup.close_holes");
        theConsole.Printf("This function closes holes. It simply follows the loop of halfvectors which");
        theConsole.Printf("enclose the hole and add a face to which they all point.");
        return "";
    }
    close_holes(active_mesh());
    return "";
}

bool console_reload(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage:  load <file>");
        theConsole.Printf("(Re)loads the current file if no argument is given, but");
        theConsole.Printf("if an argument is given, then that becomes the current file");
        return "";
    }
    if(!avo().reload(args->size() > 0 ? (*args)[0]:""))
        return "failed to load";

    return "";
}


bool console_simplify(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: simplify <fraction> ");
        theConsole.Printf("Performs Garland Heckbert (quadric based) mesh simplification.");
        theConsole.Printf("The only argument is the fraction of vertices to keep.");
        return "";
    }
    float keep_fraction;
    if(args->size() == 0) return "you must specify fraction of vertices to keep";
    istringstream a0((*args)[0]);
    a0 >> keep_fraction;

    Vec3f p0, p7;
    active_mesh().bbox(p0, p7);
    Vec3f d = p7-p0;
    float s = 1.0/d.max_coord();
    Vec3f pcentre = (p7+p0)/2.0;
    for(VertexHandle vi = active_mesh().vertices_begin(); vi != active_mesh().vertices_end(); ++vi){
        vi.get().pos = (vi.get().pos - pcentre) * s;
    }
    cout << "Timing the Garland Heckbert (quadric based) mesh simplication..." << endl;
    Timer timer;
    timer.start();

    //simplify
    quadric_simplify(active_mesh(),keep_fraction,0.0001f,true);

    cout << "Simplification complete, process time: " << timer.get_secs() << " seconds" << endl;

    //clean up the mesh, a lot of edges were jus collapsed 
        active_mesh().cleanup();

    for(VertexHandle vi = active_mesh().vertices_begin(); vi != active_mesh().vertices_end(); ++vi)
        vi.get().pos = vi.get().pos*d.max_coord() + pcentre;
    return "";
}

bool console_vertex_noise(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: noise.perturb_vertices <amplitude>");
        theConsole.Printf("adds a random vector to each vertex. A random vector in the unit cube is generated and");
        theConsole.Printf("to ensure an isotropic distribution, vectors outside the unit ball are discarded.");
        theConsole.Printf("The vector is multiplied by the average edge length and then by the amplitude specified.");
        theConsole.Printf("If no amplitude is specified, the default (0.5) is used.");
        return "";
    }
    float avg_length = average_edge_length(active_mesh());

    float noise_amplitude = 0.5f;
    if(args->size() > 0) 
    {
        istringstream a0((*args)[0]);
        a0 >> noise_amplitude;
    }

    gel_srand(0);
    for(VertexHandle vi = active_mesh().vertices_begin(); vi != active_mesh().vertices_end(); ++vi)
    {
        Vec3f v;
        do {
            v = Vec3f(gel_rand(),gel_rand(),gel_rand());
            v /= (float)(GEL_RAND_MAX);
        } while(sqr_length(v) > 1.0);
        v -= Vec3f(0.5);
        v *= 2.0;
        v *= noise_amplitude;
        v *= avg_length;
        vi.get().pos += v;
    }           
    return "";
}

bool console_perpendicular_vertex_noise(std::vector<std::string> *args)
{
    if(wantshelp(args)) 
    {
        theConsole.Printf("usage: noise.perturb_vertices_perpendicular <amplitude>");
        theConsole.Printf("adds the normal times a random scalar times amplitude times");
        theConsole.Printf("times average edge length to the vertex. (default amplitude=0.5)");
        return "";
    }
    float avg_length = average_edge_length(active_mesh());

    float noise_amplitude = 0.5;
    if(args->size() > 0) 
    {
        istringstream a0((*args)[0]);
        a0 >> noise_amplitude;
    }

    vector<Vec3f> normals(active_mesh().no_vertices());
    int i=0;
    for(VertexHandle vi = active_mesh().vertices_begin(); vi != active_mesh().vertices_end(); ++vi,++i)
        normals[i] = normal(vi);
    i=0;
    gel_srand(0);
    for(VertexHandle vi = active_mesh().vertices_begin(); vi != active_mesh().vertices_end(); ++vi,++i)
    {
        float rval = 0.5-gel_rand() / float(GEL_RAND_MAX);
        vi.get().pos += normals[i]*rval*noise_amplitude*avg_length*2.0;
    }
    return "";
}

bool console_noisy_flips(std::vector<std::string> *args)
{
    if(wantshelp(args)){
        theConsole.Printf("usage:  noise.perturb_topology <iter>");
        theConsole.Printf("Perform random flips. iter (default=1) is the number of iterations.");
        theConsole.Printf("mostly for making nasty synthetic test cases.");
        return "";
    }
    int iter = 1;
    if(args->size() > 0){
        istringstream a0((*args)[0]);
        a0 >> iter;
    }

    randomize_mesh(active_mesh(),  iter);
    return "";
}


bool console_laplacian_smooth(std::vector<std::string> *args)
{
    if(wantshelp(args)) {
        theConsole.Printf("usage:  smooth.laplacian <weight> <iter>");
        theConsole.Printf("Perform Laplacian smoothing. weight is the scaling factor for the Laplacian.");
        theConsole.Printf("default weight = 1.0. Default number of iterations = 1");
        return "";
    }
    float t=1.0;
    if(args->size() > 0){
        istringstream a0((*args)[0]);
        a0 >> t;
    }
    int iter = 1;
    if(args->size()>1){
        istringstream a0((*args)[1]);
        a0 >> iter;
    }
    /// Simple laplacian smoothing with an optional weight.
    for(int i=0;i<iter;++i) laplacian_smooth(active_mesh(), t);
    return "";
}

bool console_mean_curvature_smooth(std::vector<std::string> *args){
    if(wantshelp(args)) {
        theConsole.Printf("usage:  smooth.mean_curvature <weight> <iter>");
        theConsole.Printf("Perform mean curvature smoothing. weight is the scaling factor for the");
        theConsole.Printf("mean curvature vector which has been normalized by dividing by edge lengths");
        theConsole.Printf("this allows for larger steps as suggested by Desbrun et al.");
        theConsole.Printf("default weight = 1.0. Default number of iterations = 1");
        return "";
    }
    float t=1.0;
    if(args->size() > 0){
        istringstream a0((*args)[0]);
        a0 >> t;
    }
    int iter=1;
    if(args->size() > 1){
        istringstream a0((*args)[1]);
        a0 >> iter;
    }   
    vector<Vec3d> new_pos(active_mesh().no_vertices());
    for(int j = 0; j < iter; ++j){
        int i = 0;
        for(VertexHandle v = active_mesh().vertices_begin(); v != active_mesh().vertices_end(); ++v,++i) {
            Vec3d m;
            double w_sum;
            unnormalized_mean_curvature_normal(v, m, w_sum);
            new_pos[i] = Vec3d(v.get().pos)  - (t * m/w_sum);
        }
        i = 0;
        for(VertexHandle v = active_mesh().vertices_begin(); v != active_mesh().vertices_end(); ++v,++i)
            v.get().pos = Vec3f(new_pos[i]);
    }
    return "";
}

double angle_defect(const VertexHandle& v)
{
    if(!is_boundary(v)){
        Vec3f vertex(v.get().pos);
        vector<Vec3d> edges;
        for(VertexCirculator vc(v); !vc.end(); ++vc){
            Vec3d e(normalize(vc.vertex().get().pos - vertex));
            edges.push_back(e);
        }
        int N=edges.size();
        double angle_sum = 0;
        for(int i=0;i<N;++i){
            double dot_prod = 
                s_max(-1.0, s_min(1.0, dot(edges[i],edges[(i+1)%N])));
            angle_sum += acos(dot_prod);
        }
        return fabs(2*M_PI - angle_sum);
    }
    return 1000;
}

bool console_experimental_smooth(std::vector<std::string> *args)
{
    if(wantshelp(args)){
        theConsole.Printf("usage:  smooth.experimental <weight> <iter>");
        theConsole.Printf("Perform experimental smoothing. weight is the scaling factor for the");
        theConsole.Printf("experimetnal vector which has been normalized by dividing by edge lengths");
        theConsole.Printf("this allows for larger steps as suggested by Desbrun et al.");
        theConsole.Printf("default weight = 1.0. Default number of iterations = 1");
        return "";
    }
    float t = 1.0f;
    if(args->size() > 0){
        istringstream a0((*args)[0]);
        a0 >> t;
    }
    int iter = 1;
    if(args->size() > 1){
        istringstream a0((*args)[1]);
        a0 >> iter;
    }   
    Vec3f p0_before, p7_before;
    active_mesh().bbox(p0_before, p7_before);

    double avg_area = 0.0;
    for(FaceHandle f = active_mesh().faces_begin(); f != active_mesh().faces_end(); ++f){
        avg_area += area(f);
    }
    avg_area /= active_mesh().no_faces();

    for(int j=0;j<iter;++j){
        vector<Vec3d> new_pos(active_mesh().no_vertices());

        int i=0;
        for(VertexHandle v = active_mesh().vertices_begin(); v != active_mesh().vertices_end(); ++v,++i) {
            double w_sum=0;
            Vec3d m(0);
            VertexCirculator vc(v);
            for(; !vc.end(); ++vc){
                HalfEdgeHandle h = vc.halfedge();
                double area_left  = area(h.face());
                double area_right = area(h.opp().face());

                double w = 0.5 * (area_left + area_right);
                m += w * Vec3d(vc.vertex().get().pos - v.get().pos);
                w_sum += w;
            }
            double t2 = pow((w_sum/vc.no_steps())/avg_area,.25);
            if(t2 > 1e-10)
                new_pos[i] = Vec3d(v.get().pos)  + (t * t2 * m/w_sum);
        }
        i=0;
        for(VertexHandle v = active_mesh().vertices_begin(); v != active_mesh().vertices_end(); ++v,++i)
            v.get().pos = Vec3f(new_pos[i]);

        for(HalfEdgeHandle h=active_mesh().halfedges_begin(); h!=active_mesh().halfedges_end(); ++h){

            Vec3d mcna,mcnb;
            double ws;

            unnormalized_mean_curvature_normal(h.vert(), mcna, ws);
            unnormalized_mean_curvature_normal(h.opp().vert(), mcnb, ws);
            if(sqr_length(mcna) < 1 && sqr_length(mcnb) < 1){
                if(active_mesh().collapse_halfedge_test(h)){
                    active_mesh().collapse_halfedge(h, false);
                }
            }
        }
        //new manifold cleanup call
        active_mesh().cleanup();

        maximize_min_angle(active_mesh(), 0.95f, false);
    }

    Vec3f p0, p7;
    active_mesh().bbox(p0, p7);
    for(VertexHandle vi = active_mesh().vertices_begin(); vi != active_mesh().vertices_end(); ++vi)
        vi.get().pos = (p7_before-p0_before)*(vi.get().pos - p0)/(p7-p0)+p0_before;

    return "";
}


bool console_taubin_smooth(std::vector<std::string> *args)
{
    if(wantshelp(args)){
        theConsole.Printf("usage:  smooth.taubin <iter>");
        theConsole.Printf("Perform Taubin smoothing. iter (default=1) is the number of iterations.");
        return "";
    }
    int iter = 1;
    if(args->size() > 0){
        istringstream a0((*args)[0]);
        a0 >> iter;
    }
    /// Taubin smoothing is similar to laplacian smoothing but reduces shrinkage
    taubin_smooth(active_mesh(),  iter);

    return "";
}

bool console_fvm_smooth(std::vector<std::string> *args)
{       
    if(wantshelp(args)){
        theConsole.Printf("usage: smooth.fuzzy_vector_median <iter>");
        theConsole.Printf("Smooth normals using fuzzy vector median smoothing. iter (default=1) is the number of iterations");
        theConsole.Printf("This function does a very good job of preserving sharp edges.");
        return "";
    }
    int iter=1;
    if(args->size() > 0){
        istringstream a0((*args)[0]);
        a0 >> iter;
    }
    // Fuzzy vector median smoothing is effective when it comes to preserving sharp edges. 
    fvm_smooth(active_mesh(),  iter);

    return "";
}

bool console_triangulate(std::vector<std::string> *args)
{       
    if(wantshelp(args)) {
        theConsole.Printf("usage:  triangulate");
        theConsole.Printf("This function triangulates all non triangular faces of the mesh.");
        theConsole.Printf("you may want to call it after hole closing. For a polygon it simply connects");
        theConsole.Printf("the two closest vertices in a recursive manner until only triangles remain");
        return "";
    }
    shortest_edge_triangulate(active_mesh());
    active_mesh().cleanup();

    return "";
}


bool console_remove_caps(std::vector<std::string> *args)
{       
    if(wantshelp(args)) {
        theConsole.Printf("usage:  cleanup.remove_caps thresh");
        theConsole.Printf("Remove caps (triangles with one very big angle). The thresh argument is the fraction of PI to");
        theConsole.Printf("use as threshold for big angle. Default is 0.85. Caps are removed by flipping.");
        return "";
    }
    float t = 0.85f;
    if(args->size() > 0){
        istringstream a0((*args)[0]);
        a0 >> t;
    }
    remove_caps(active_mesh(), static_cast<float>(M_PI) *t);
    active_mesh().cleanup();

    return "";
}

bool console_remove_needles(std::vector<std::string> *args)
{       
    if(wantshelp(args)){
        theConsole.Printf("usage: cleanup.remove_needles <thresh>");
        theConsole.Printf("Removes very short edges by collapse. thresh is multiplied by the average edge length");
        theConsole.Printf("to get the length shorter than which we collapse. Default = 0.1");
        return "";
    }
    float thresh = 0.1f;
    if(args->size() > 0){
        istringstream a0((*args)[0]);
        a0 >> thresh;
    }
    float avg_length = average_edge_length(active_mesh());
    remove_needles(active_mesh(), thresh * avg_length);
    active_mesh().cleanup();

    return "";
}

void reshape(int W, int H)
{
    active_view_control().reshape(W,H);
}

void display() 
{
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    static string& display_render_mode = CreateCVar<string>("display.render_mode","");
    static int& display_smooth= CreateCVar<int>("display.smooth_shading",1);

    glPushMatrix();

    //cout << display_render_mode << endl;
    avo().display(display_render_mode, display_smooth);

    glPopMatrix();

    glUseProgram(0);
    theConsole.RenderConsole();

    glutSwapBuffers();
}

void animate() 
{       
    //usleep( (int)1e4 );
    active_view_control().try_spin();
    glutPostRedisplay();
}


void mouse(int button, int state, int x, int y) 
{
    Vec2i pos(x,y);
    if (state==GLUT_DOWN) 
    {
        if (button==GLUT_LEFT_BUTTON) 
            active_view_control().grab_ball(ROTATE_ACTION,pos);
        else if (button==GLUT_MIDDLE_BUTTON) 
            active_view_control().grab_ball(ZOOM_ACTION,pos);
        else if (button==GLUT_RIGHT_BUTTON) 
            active_view_control().grab_ball(PAN_ACTION,pos);
    }
    else if (state==GLUT_UP)
        active_view_control().release_ball();
}

void motion(int x, int y) {
    Vec2i pos(x,y);
    active_view_control().roll_ball(Vec2i(x,y));
}

void keyboard_spec(int key, int x, int y)
{
    int mod = glutGetModifiers();
    if( theConsole.IsOpen() ) {
        // If shift held, scroll the console
        if( mod == GLUT_ACTIVE_SHIFT ) {
            switch (key){
                case GLUT_KEY_UP:
                    theConsole.ScrollDownLine();
                    break;
                case GLUT_KEY_DOWN: 
                    theConsole.ScrollUpLine();
                    break;
            }
        } else {
            theConsole.SpecialFunc( key );
        }
    }
}


void keyboard(unsigned char key, int x, int y) 
{       
    if(theConsole.IsOpen())
    {
        switch(key) {
            case '\033': 
                theConsole.ToggleConsole();
            default:
                theConsole.KeyboardFunc(key);
                break;
        }
        if(key == 13)   avo().post_create_display_list();

    }   
    else {
        string& display_render_mode = GetCVarRef<string>("display.render_mode");
        int& display_smooth = GetCVarRef<int>("display.smooth_shading");
        int& active  = GetCVarRef<int>("active_mesh");


        switch(key) {
            case 'q': exit(0);
            case '\033':
                theConsole.ToggleConsole();
                break;
            case '1':
            case '2':
            case '3':
            case '4':
            case '5':
            case '6':
            case '7':
            case '8':
            case '9':
                active = key - '1'; break;
            case 'f': display_smooth = !display_smooth; break;
            case 'w':
                display_render_mode = "wire"; break;
            case 'n':
                display_render_mode = "normal"; break;
            case 'i':
                display_render_mode = "isophotes"; break;
            case 'r':
                display_render_mode = "reflection"; break;
            case 'h':
                display_render_mode = "harmonics"; break;
            case 't':
                display_render_mode = "toon"; break;
            case 'g':
                display_render_mode = "glazed"; break;
            case 'a':
                display_render_mode = "ambient_occlusion"; break;
            case 'c':
                display_render_mode = "copper"; break;
            case 'C':
                display_render_mode = "curvature_lines"; break;
            case 'M':
                display_render_mode = "mean_curvature"; break;
            case 'G':
                display_render_mode = "gaussian_curvature"; break;
        }

        if(display_render_mode.substr(0,3) == "har")
            avo().harmonics_parse_key(key);

        if(key != '\033') avo().post_create_display_list();             
    }
}

void init_glut(int argc, char** argv)
{  
    glutInitDisplayMode(GLUT_RGBA|GLUT_DOUBLE|GLUT_DEPTH|GLUT_ALPHA);
    glutInitWindowSize(WINX, WINY);
    glutInit(&argc, argv);
    glutCreateWindow("MeshEdit");
    glutDisplayFunc(display);
    glutKeyboardFunc(keyboard);
    glutSpecialFunc(keyboard_spec);
    glutReshapeFunc(reshape);
    glutMouseFunc(mouse);
    glutMotionFunc(motion);
    glutIdleFunc(animate);
}
void init_gl()
{
    glewInit();
    glEnable(GL_LIGHTING);
    glEnable(GL_LIGHT0);
    glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, 1);

    // Set the value of a uniform
    //glUniform2f(glGetUniformLocation(prog_P0,"WIN_SCALE"), win_size_x/2.0, win_size_y/2.0);

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glClearColor(0.50f, 0.50f, 0.50f, 0.f);
    glColor4f(1.0f, 1.0f, 1.0f, 0.f);
    float material[4] = {1,1,1,1};
    glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, material);
    glEnable(GL_DEPTH_TEST);

//    CreateCVar("help", MyConsoleHelp );
    CreateCVar("harmonics.reset_shape", console_reset_shape);
    CreateCVar("harmonics.analyze", console_analyze);
    CreateCVar("harmonics.partial_reconstruct", console_partial_reconstruct);
    CreateCVar("simplify", console_simplify);
    CreateCVar("smooth.mean_curvature", console_mean_curvature_smooth);
    CreateCVar("smooth.experimental", console_experimental_smooth);
    CreateCVar("smooth.laplacian", console_laplacian_smooth);
    CreateCVar("smooth.taubin", console_taubin_smooth);
    CreateCVar("smooth.fuzzy_vector_median", console_fvm_smooth);

    CreateCVar("optimize.valency", console_optimize_valency);
    CreateCVar("optimize.minimize_dihedral_angles", console_minimize_dihedral);
    CreateCVar("optimize.minimize_curvature", console_minimize_curvature);
    CreateCVar("optimize.maximize_min_angle", console_maximize_min_angle);
    CreateCVar("cleanup.close_holes", console_close_holes);
    CreateCVar("load_mesh", console_reload);

    CreateCVar("cleanup.remove_caps", console_remove_caps);
    CreateCVar("cleanup.remove_needles", console_remove_needles);
    CreateCVar("triangulate", console_triangulate);
    CreateCVar("refine.split_edges", console_refine_edges);
    CreateCVar("refine.split_faces", console_refine_faces);
    CreateCVar("refine.catmull_clark", console_cc_subdivide);
    CreateCVar("save_mesh", console_save);
    CreateCVar("noise.perturb_vertices", console_vertex_noise);
    CreateCVar("noise.perturb_vertices_perpendicular", console_perpendicular_vertex_noise);
    CreateCVar("noise.perturb_topology", console_noisy_flips);

    CreateCVar("dual", console_dual);
    CreateCVar("flatten", console_flatten);

    CreateCVar("align", console_align);

}

int main(int argc, char** argv)
{

    ArgExtracter ae(argc, argv);

    init_glut(argc, argv);
    init_gl();

    Harmonics::init();

    if(argc>1){         
        string file = ae.get_last_arg();
        avo().reload(file);
    }
    glutMainLoop();
    return 0;
}