Subversion Repositories gelsvn

/*

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

#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 <Util/Timer.h>

#include <Util/ArgExtracter.h>

#include <GL/glew.h>

#include <GLGraphics/gel_glut.h>

#include <GLGraphics/draw.h>

#include <GLGraphics/glsl_shader.h>

#include <GLGraphics/GLViewController.h>

#include <HMesh/Manifold.h>

#include <HMesh/VertexCirculator.h>

#include <HMesh/FaceCirculator.h>

#include <HMesh/build_manifold.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/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 <GLConsole/GLConsole.h>

#include <Util/Timer.h>

#include "harmonics.h"

#include "wireframe.h"

using namespace std;

using namespace HMesh;

using namespace Geometry;

using namespace GLGraphics;

using namespace CGLA;

using namespace Util;

using namespace LinAlg;

int WINX=800, WINY=800;

class VisObj

{

	string file;

	GLViewController view_ctrl;

	GLuint display_list;

	bool create_display_list;

	Manifold mani;

	Harmonics* harmonics;

public:

	Manifold& mesh() {return mani;}

	GLViewController& view_control() {return view_ctrl;}

	bool reload(string _file)

	{

		if(_file != "") file = _file;

		mani.clear();

		if(!load(file, mani))

			return false;

		Vec3f c(0,0,0);

		float r = 5;

		mani.get_bsphere(c,r);

		view_ctrl.set_centre(c);

		view_ctrl.set_eye_dist(2*r);

		return true;

	}

	VisObj():

	file(""), view_ctrl(WINX,WINY, Vec3f(0), 1.0), display_list(glGenLists(1)), create_display_list(true), harmonics(0) 

	{

	}

	void display(bool wire, bool harm, bool flat)

	{

/*		static Timer tim;

		tim.start();*/

		if(create_display_list)

		{

			create_display_list = false;

			glNewList(display_list,GL_COMPILE);

			if(wire)

			{

				draw_as_wire(mani,false);

				glUseProgram(0);	

			}

			else if(harm)

				harmonics->draw();

			else 

				draw(mani,!flat);

			glEndList();

		}

		view_ctrl.set_gl_modelview();

		glCallList(display_list);

		// tim.get_secs();

	}

	void post_create_display_list()

	{

		create_display_list = true;

	}

	void harmonics_analyze_mesh()

	{

		delete harmonics;

		harmonics = new Harmonics(mani);

	}

	void harmonics_reset_shape()

	{

		if(harmonics)

			harmonics->reset_shape();

	}

	void harmonics_parse_key(unsigned char key)

	{

		harmonics->parse_key(key);

	}

	void harmonics_partial_reconstruct(int eig0, int eig1, float scale)

	{

		if(harmonics)

			harmonics->partial_reconstruct(eig0, eig1, scale);

	}

};

inline VisObj& get_vis_obj(int i)

{

	static VisObj vo[9];

	return vo[i];

}

inline VisObj& avo()

{

	static CVar<int> active("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

GLConsole theConsole;

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

char* ConsoleHelp(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   : toggle wireframe");

    theConsole.Printf("f   : toggle flatshading");

    theConsole.Printf("1-9 : switch between active meshes.");

    theConsole.Printf("d   : (display.show_harmonics = 1) diffuse light on and off");

    theConsole.Printf("h   : (display.show_harmonics = 1) highlight on and off ");

    theConsole.Printf("+/- : (display.show_harmonics = 1) 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.

char* 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 "";

}

char* 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())==".x3d")

		{

			x3d_save(file_name, active_mesh());

			return "";

		}

		return "unknown format";

	}

	return "usage: save <name.x3d|name.obj> ";

}

char* 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.5;

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

}

char* 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 "";

		}

	safe_triangulate(active_mesh());

	return "";

}

char* 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 "";

}

char* 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<int> faces(m.no_vertices());

	vector<int> indices;

	// Create new vertices. Each face becomes a vertex whose position

	// is the centre of the face

	int i=0;

	for(FaceIter 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.

	i=0;

	for(VertexIter v=m.vertices_begin(); v!= m.vertices_end(); ++v,++i)

	{

		VertexCirculator vc(v);

		vector<int> index_tmp;

		for(; !vc.end(); ++vc)

			index_tmp.push_back(vc.get_face()->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[i] = vc.no_steps();

	}

	// Clear the manifold before new geometry is inserted.

	m.clear();

	// And build

	build_manifold(m, vertices.size(), &vertices[0], faces.size(),

				   &faces[0],&indices[0]);

	return "";

}

char* 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 "";

}

char* 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.0;

	if(args.size()>3)

	{

		istringstream a0(args[0]);

		a0 >> gamma;

	}

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

	return "";

}

char* 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.0;

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

}

char* 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 "";

}

char* 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 "";

}

char* console_partial_reconstruct(std::vector<std::string> &args)

{

	if(wantshelp(args)) 

		{

			theConsole.Printf("usage: haramonics.partial_reconstruct <e0> <e1> <s>");

			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.");

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

}

char* 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 "";

}

char* 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 "";

}

char* console_reload(std::vector<std::string> &args)

{

	if(wantshelp(args)) 

		{

			theConsole.Printf("usage:  reload <file>");

			theConsole.Printf("Reloads 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 "";

}

char* 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().get_bbox(p0, p7);

	Vec3f d = p7-p0;

	float s = 1.0/d.max_coord();

	Vec3f pcentre = (p7+p0)/2.0;

	for(VertexIter vi = active_mesh().vertices_begin(); vi != active_mesh().vertices_end(); ++vi)

		vi->pos = (vi->pos - pcentre) * s;

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

	for(VertexIter vi = active_mesh().vertices_begin(); vi != active_mesh().vertices_end(); ++vi)

		vi->pos = vi->pos*d.max_coord() + pcentre;

	return "";

}

char* 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. To ensure uniformness, the vector must lie in the");

			theConsole.Printf("unit sphere. The length of the vector is multiplied by the average edge length and then amplitude");

			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;

	}

	srand(0);

	for(VertexIter vi = active_mesh().vertices_begin(); vi != active_mesh().vertices_end(); ++vi)

	{

		Vec3f v;

		do {

			v = Vec3f(rand(),rand(),rand());

			v /= RAND_MAX;

		} while(sqr_length(v) > 1.0);

		v -= Vec3f(0.5);

		v *= 2.0;

		v *= noise_amplitude;

		v *= avg_length;

		vi->pos += v;

	}		

	return "";

}

char* 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 "";

}

char* console_laplacian_smooth(std::vector<std::string> &args)

{

	if(wantshelp(args)) 

		{

			theConsole.Printf("usage:  smooth.laplacian <weight>");

			theConsole.Printf("Perform Laplacian smoothing. weight is the scaling factor for the Laplacian.");

			return "";

		}

	float t=1.0;

	if(args.size()>0)

	{

		istringstream a0(args[0]);

		a0 >> t;

	}

	/// Simple laplacian smoothing with an optional weight.

	laplacian_smooth(active_mesh(), t);

	return "";

}

char* 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 "";

}

char* 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 "";

}

char* 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());

	return "";

}

char* 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.85;

	if(args.size()>0)

	{

		istringstream a0(args[0]);

		a0 >> t;

	}

	remove_caps_from_trimesh(active_mesh(), static_cast<float>(M_PI) *t);

	return "";

}

char* 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.1;

	if(args.size()>0)

	{

		istringstream a0(args[0]);

		a0 >> thresh;

	}

	float avg_length = average_edge_length(active_mesh());

	remove_needles_from_trimesh(active_mesh(), thresh * avg_length);

	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 CVar<int> display_wireframe("display.wireframe",0);

	static CVar<int> display_eigenmodes("display.show_harmonics",0);

	static CVar<int> display_flat("display.flatshading",0);

	glPushMatrix();

	avo().display(display_wireframe, display_eigenmodes, display_flat);

	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();