Subversion Repositories gelsvn

//

//  glut_main.cpp

//  GEL

//

//  Created by J. Andreas Bærentzen on 04/10/13.

//

//

#include "glut_main.h"

/*

 *  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 <algorithm>

#include <queue>

#include <GL/glew.h>

#include <GLGraphics/Console.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 <CGLA/Mat4x4d.h>

#include <LinAlg/Matrix.h>

#include <LinAlg/Vector.h>

#include <LinAlg/LapackFunc.h>

#include <GLGraphics/gel_glut.h>

#include <HMesh/Manifold.h>

#include <HMesh/AttributeVector.h>

#include <HMesh/mesh_optimization.h>

#include <HMesh/curvature.h>

#include <HMesh/triangulate.h>

#include <HMesh/flatten.h>

#include <HMesh/dual.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/cleanup.h>

#include <HMesh/cleanup.h>

#include <HMesh/refine_edges.h>

#include <HMesh/subdivision.h>

#include <Util/Timer.h>

#include <Util/ArgExtracter.h>

#include "polarize.h"

#include "harmonics.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;

// Single global instance so glut can get access

Console theConsole;

bool console_visible = false;

inline VisObj& get_vis_obj(int i)

{

    static VisObj vo[9];

    return vo[i];

}

Console::variable<int> active(0);

inline VisObj& avo()

{

    return get_vis_obj(active);

}

inline Manifold& active_mesh()

{

    return avo().mesh();

}

inline GLViewController& active_view_control()

{

    return avo().view_control();

}

////////////////////////////////////////////////////////////////////////////////

bool MyConsoleHelp(const 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 true;

}

bool wantshelp(const 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.

void console_align(const 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");

    }

    int dest = 0;

    if(args.size()>0){

        istringstream a0(args[0]);

        a0 >> dest;

        --dest;

        if(dest <0 || dest>8)

        {

            theConsole.printf("dest mesh out of range (1-9)");

            return;

        }

    }

    else

    {

        theConsole.printf("neither source nor destination mesh?!");

        return;

    }

    int src = 0;

    if(args.size()>1){

        istringstream a1(args[1]);

        a1 >> src;

        --src;

        if(src <0 || src>8)

        {

            theConsole.printf("src mesh out of range (1-9)");

            return;

        }

    }

    else

    {

        theConsole.printf("no src mesh?");

        return;

    }

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

}

void console_polarize(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: polarize");

        return;

    }

    int divisions = 50;

    if(args.size() > 0){

        istringstream a0(args[0]);

        a0 >> divisions;

    }

    double t=1;

    if(args.size() > 1){

        istringstream a0(args[1]);

        a0 >> t;

    }

    avo().save_old();

	double vmin, vmax;

    VertexAttributeVector<double> fun;

    VertexAttributeVector<Vec2d> par;

    make_adf_fun(active_mesh(), t, fun, vmin, vmax);

    //    polarize_mesh_new(active_mesh(), fun, vmin, vmax, divisions, par);

    polarize_mesh(active_mesh(), fun, vmin, vmax, divisions, par);

}

void console_simplify_polar(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: simplify.polar <frac>");

        return;

    }

    double frac = 0.9;

    if(args.size() > 0){

        istringstream a0(args[0]);

        a0 >> frac;

    }

    int iter=1;

    if(args.size() > 1){

        istringstream a0(args[1]);

        a0 >> iter;

    }

    avo().save_old();

    simplify_polar_mesh(avo().mesh(), frac, iter);

}

void console_polar_segment(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: polar.segment");

        return;

    }

    int segments = 1;

    if(args.size() > 0){

        istringstream a0(args[0]);

        a0 >> segments;

    }

    polar_segment(avo().mesh(), segments);

}

void console_polar_skeleton(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: polar.skeleton <frac>");

        return;

    }

    double frac = 0.9;

    if(args.size() > 0){

        istringstream a0(args[0]);

        a0 >> frac;

    }

    avo().save_old();

    skeleton_retract(avo().mesh(), frac);

}

void console_polar_subdivide(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: polar.subdivide <iter>");

        return;

    }

    int iter=1;

    if(args.size() > 0){

        istringstream a0(args[0]);

        a0 >> iter;

    }

    avo().save_old();

    polar_subdivide (avo().mesh(), iter);

}

void console_ridge_lines(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: ridge_lines");

        return;

    }

    avo().save_old();

    Manifold& mani = avo().mesh();

    VertexAttributeVector<Mat3x3d> curvature_tensors(mani.allocated_vertices());

    VertexAttributeVector<Vec3d> min_curv_direction(mani.allocated_vertices());

    VertexAttributeVector<Vec3d> max_curv_direction(mani.allocated_vertices());

    VertexAttributeVector<Vec2d> curvature(mani.allocated_vertices());

    //    curvature_tensors_from_edges(mani, curvature_tensors);

    //    for(int i=0;i<3; ++i)

    //        smooth_curvature_tensors(mani,curvature_tensors);

    //    curvature_from_tensors(mani, curvature_tensors,

    //                        min_curv_direction,

    //                        max_curv_direction,

    //                        curvature);

    curvature_paraboloids(mani,

                          min_curv_direction,

                          max_curv_direction,

                          curvature);

    for(auto vid : mani.vertices())

    {

        Vec3d max_curv_dir = normalize(max_curv_direction[vid]);

        Vec3d min_curv_dir = normalize(min_curv_direction[vid]);

        double vid_min_pc = curvature[vid][0];

        double vid_max_pc = curvature[vid][1];

        bool ridge = true;

        bool ravine = true;

        Walker w = mani.walker(vid);

        double wsum =0;

        Vec3d r(0);

        for(; !w.full_circle();w = w.circulate_vertex_ccw())

        {

            Vec3d e = (mani.pos(w.vertex()) - mani.pos(vid));

            if(abs(dot(min_curv_dir,e)) > abs(dot(max_curv_dir,e)))

            {

                if(curvature[w.vertex()][0]<vid_min_pc+20)

                    ravine = false;

            }

            else

            {

                if(curvature[w.vertex()][1]>vid_max_pc-20)

                    ridge = false;

            }

        }

        DebugRenderer::vertex_colors[vid] = Vec3f(ridge,ravine,0.0);

    }

    for(auto fid : mani.faces())

        DebugRenderer::face_colors[fid] = Vec3f(.3,.3,.6);

    for(auto hid : mani.halfedges()) {

        Walker w = mani.walker(hid);

        Vec3f c0 = DebugRenderer::vertex_colors[w.opp().vertex()];

        Vec3f c1 = DebugRenderer::vertex_colors[w.vertex()];

        DebugRenderer::edge_colors[hid] = (c0==c1) ? c0 : Vec3f(0.1,0.1,0.3);

    }

}

void console_refit_polar(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: simplify.polar <mesh 1> <mesh 2> <iter>");

        return;

    }

    int m1=1;

    int m2=2;

    int iter=1;

    int dim = 64;

    if(args.size() > 0){

        istringstream a0(args[0]);

        a0 >> m1;

    }

    if(args.size() > 1){

        istringstream a0(args[1]);

        a0 >> m2;

    }

    if(args.size() > 2){

        istringstream a0(args[2]);

        a0 >> iter;

    }

    if(args.size() > 3){

        istringstream a0(args[3]);

        a0 >> dim;

    }

    avo().save_old();

    smooth_and_refit(get_vis_obj(m1-1).mesh() , get_vis_obj(m2-1).mesh(), iter, dim);

}

void transform_mesh(Manifold& mani, const Mat4x4d& m)

{

    for(VertexIDIterator vid = mani.vertices_begin(); vid != mani.vertices_end(); ++vid)

        mani.pos(*vid) = m.mul_3D_point(mani.pos(*vid));

}

void console_scale(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: scale sx sy sz");

        return;

    }

    Vec3d s;

    if(args.size() > 0){

        istringstream a0(args[0]);

        a0 >> s[0];

    }

    if(args.size() > 1){

        istringstream a0(args[0]);

        a0 >> s[1];

    }

    if(args.size() > 2){

        istringstream a0(args[0]);

        a0 >> s[2];

    }

    avo().save_old();

    transform_mesh(avo().mesh(),scaling_Mat4x4d(s));

    avo().refit();

}

void console_flatten(const 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;

    }

    avo().save_old();

    WeightScheme ws = BARYCENTRIC_W;

    if(args.size()>0){

        if(args[0] == "floater")

            ws = FLOATER_W;

        else if(args[0] == "harmonic")

            ws = HARMONIC_W;

        else if(args[0] == "lscm")

            ws = LSCM_W;

    }

    else

        return;

    flatten(active_mesh(), ws);

    return;

}

void console_save(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: save <name.x3d|name.obj> ");

        return;

    }

    const 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;

        }

        theConsole.printf("unknown format");

        return;

    }

    theConsole.printf("usage: save <name.x3d|name.obj> ");

}

void console_refine_edges(const 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;

    }

    avo().save_old();

    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;

}

void console_refine_faces(const 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;

    }

    avo().save_old();

    triangulate_by_vertex_face_split(active_mesh());

    return;

}

void console_cc_subdivide(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: subdivide.catmull_clark ");

        theConsole.printf("Does one step of Catmull-Clark subdivision");

        return;

    }

    avo().save_old();

    cc_split(active_mesh(),active_mesh());

    cc_smooth(active_mesh());

    return;

}

void console_loop_subdivide(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: subdivide.loop");

        theConsole.printf("Does one step of Loop subdivision");

        return;

    }

    avo().save_old();

    loop_split(active_mesh(),active_mesh());

    loop_smooth(active_mesh());

    return;

}

void console_root3_subdivide(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: subdivide.root3");

        theConsole.printf("Does one step of sqrt(3) subdivision");

        return;

    }

    avo().save_old();

    root3_subdivide(active_mesh(),active_mesh());

    return;

}

void console_doosabin_subdivide(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: subdivide.doo_sabin ");

        theConsole.printf("Does one step of Doo-Sabin Subdivision");

        return;

    }

    avo().save_old();

    cc_split(active_mesh(),active_mesh());

    dual(active_mesh());

    return;

}

void console_butterfly_subdivide(const std::vector<std::string> & args)

{

    if(wantshelp(args)) {

        theConsole.printf("usage: subdivide.butterfly ");

        theConsole.printf("Does one step of Modified Butterfly Subdivision");

        return;

    }

    avo().save_old();

    butterfly_subdivide(active_mesh(),active_mesh());

    return;

}

void console_dual(const 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;

    }

    avo().save_old();

    dual(active_mesh());

    return;

}

void console_minimize_curvature(const 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;

    }

    avo().save_old();

    bool anneal=false;

    if(args.size() > 0)

    {

        istringstream a0(args[0]);

        a0 >> anneal;

    }

    minimize_curvature(active_mesh(), anneal);

    avo().post_create_display_list();

    return;

}

void console_minimize_dihedral(const 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;

    }

    avo().save_old();

    int iter = 1000;

    if(args.size() > 0)

    {

        istringstream a0(args[0]);

        a0 >> iter;

    }

    bool anneal = false;

    if(args.size() > 1)

    {

        istringstream a0(args[1]);

        a0 >> anneal;

    }

    bool use_alpha = true;

    if(args.size() > 2)

    {

        istringstream a0(args[2]);

        a0 >> use_alpha;

    }

    float gamma = 4.0f;

    if(args.size() > 3)

    {

        istringstream a0(args[3]);

        a0 >> gamma;

    }

    minimize_dihedral_angle(active_mesh(), iter, anneal, use_alpha, gamma);

    return;

}

void console_maximize_min_angle(const 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;

    }

    avo().save_old();

    float thresh = 0.0f;

    if(args.size() > 0)

    {

        istringstream a0(args[0]);

        a0 >> thresh;

    }

    bool anneal = false;

    if(args.size() > 1)

    {

        istringstream a0(args[1]);

        a0 >> anneal;

    }

    maximize_min_angle(active_mesh(),thresh,anneal);

    return;

}

void console_optimize_valency(const 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;

    }

    avo().save_old();

    bool anneal = false;

    if(args.size() > 0)

    {

        istringstream a0(args[0]);

        a0 >> anneal;

    }

    optimize_valency(active_mesh(), anneal);

    return;

}

void console_analyze(const 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;

}

void console_partial_reconstruct(const 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.");

    }

    avo().save_old();

    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;

}

void console_reset_shape(const 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().save_old();

    avo().harmonics_reset_shape();

    return;

}

void console_close_holes(const 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;

    }

    avo().save_old();

    close_holes(active_mesh());

    return;

}

void console_reload(const 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;

    }

    avo().save_old();

    if(!avo().reload(args.size() > 0 ? args[0]:""))

        theConsole.printf("failed to load");

    return;

}

void console_add_mesh(const std::vector<std::string> & args)

{

    if(wantshelp(args))

    {