WebSVN – gelsvn – Diff – /trunk/apps/MeshEdit/harmonics.cpp

 /*
  *  harmonics.cpp
  *  GEL
+ *
  *  Created by J. Andreas Bærentzen on 01/09/08.
  *  Copyright 2008 __MyCompanyName__. All rights reserved.
+ *
  */
 #include <iostream>
 #include <CGLA/Vec3f.h>
 #include <CGLA/Vec3d.h>
 #include <LinAlg/Matrix.h>
 #include <LinAlg/Vector.h>
 #include <LinAlg/LapackFunc.h>
 #include <GL/glew.h>
 #include <GLGraphics/glsl_shader.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/x3d_save.h>
 #include <GLConsole/GLConsole.h>
 #include "harmonics.h"
 #include "wireframe.h"
 using namespace CGLA;
 using namespace std;
 using namespace HMesh;
 using namespace Geometry;
 using namespace GLGraphics;
 using namespace LinAlg;
-int E = 1;
+bool Harmonics::is_initialized=false;
-int MAX_E;
+GLuint Harmonics::prog_P0;
 namespace
+{
-	GLuint prog_P0;
-	vector<double> projx;
-	vector<double> projy;
-	vector<double> projz;
-	CMatrix Q;
-	CVector V;
-	vector<float> max_eig_values;
-	vector<Vec3f> vertex_positions;
 	double voronoi_area(VertexIter v)
+	{
 	double area_mixed = 0;
 	//For each triangle T from the 1-ring neighborhood of x
 	for(VertexCirculator vc(v); !vc.end(); ++vc)
+	{
 		FaceIter f = vc.get_face();
 		double f_area = area(f);
 		HalfEdgeIter he = vc.get_halfedge();
 		Vec3d v1(he->vert->pos);
 		Vec3d v2(he->next->vert->pos);
 		Vec3d v0(he->next->next->vert->pos);
 		double a0 = acos(dot(v1-v0, v2-v0)/(length(v1-v0)*length(v2-v0)));
 		double a1 = acos(dot(v2-v1, v0-v1)/(length(v2-v1)*length(v0-v1)));
 		double a2 = acos(dot(v0-v2, v1-v2)/(length(v0-v2)*length(v1-v2)));
 		if(a0>(M_PI/2.0) && a1>(M_PI/2.0) && a2>(M_PI/2.0)) // f is non-obtuse
+		{
 			// Add Voronoi formula (see Section 3.3)
 			area_mixed += (1.0/8) *
 			((1.0/tan(a1)) * sqr_length(v2-v0) +
 			 (1.0/tan(a2)) * sqr_length(v1-v0));
+		}
 		else // Voronoi inappropriate
+		{
 			// Add either area(f)/4 or area(f)/2
 			if(a0>M_PI/2.0)// the angle of f at x is obtuse
 				area_mixed += f_area/2;
 			else
 				area_mixed += f_area/4;
+		}
+	}
 	return area_mixed;
+}
 	double barycentric_area(VertexIter v)
+	{
 	double barea = 0;
 	//For each triangle T from the 1-ring neighborhood of x
 	for(VertexCirculator vc(v); !vc.end(); ++vc)
+	{
 		FaceIter f = vc.get_face();
 		barea += area(f)/3.0;
+	}
 	return barea;
+}
+}
-void make_laplace_operator(Manifold& mani)
+void Harmonics::make_laplace_operator()
+{
 	Q.Resize(mani.no_vertices(), mani.no_vertices());
 	for(VertexIter v = mani.vertices_begin(); v != mani.vertices_end(); ++v)
 		if(!is_boundary(v))
+		{
 			int i = v->touched;
 			double area_i = voronoi_area(v);
 			Vec3d vertex(v->pos);
 			Vec3d curv_normal(0);
 			double a_sum = 0;
 			for(VertexCirculator vc(v); !vc.end(); ++vc)
+			{
 				int j = vc.get_vertex()->touched;
 				double area_j = voronoi_area(vc.get_vertex());
 				HalfEdgeIter h = vc.get_halfedge();
 				Vec3d nbr(h->vert->pos);
 				Vec3d left(h->next->vert->pos);
 				Vec3d right(h->opp->prev->opp->vert->pos);
 				double d_left = dot(normalize(nbr-left),
 									normalize(vertex-left));
 				double d_right = dot(normalize(nbr-right),
 									 normalize(vertex-right));
 				double a_left  = acos(min(1.0, max(-1.0, d_left)));
 				double a_right = acos(min(1.0, max(-1.0, d_right)));
 				double w = 1.0/tan(a_left) + 1.0/tan(a_right);
 				Q[i][j] = -w/sqrt(area_i*area_j);
 				//Q[i][j] = -1;
 				a_sum += Q[i][j];
+			}
 			Q[i][i] = -a_sum;
+		}
 	EigenSolutionsSym(Q,V);
+}
-void analyze_mesh(Manifold& mani)
+Harmonics::Harmonics(Manifold& _mani):mani(_mani)
+{
-	static bool was_here = false;
+	assert(is_initialized);
-	if(!was_here)
-	{
-		triangulate(mani);
+	triangulate(mani);
-		mani.enumerate_vertices();
+	mani.enumerate_vertices();
-		MAX_E = mani.no_vertices()-1;
+	maximum_eigenvalue = mani.no_vertices()-1;
-		E = min(E, MAX_E);
-		was_here = true;
-		make_laplace_operator(mani);
+	make_laplace_operator();
-		vertex_positions.resize(mani.no_vertices());
-	}
-	projx.resize(MAX_E);
-	projy.resize(MAX_E);
-	projz.resize(MAX_E);
+	proj.resize(maximum_eigenvalue);
-	max_eig_values.resize(MAX_E, 1e-10);
+	max_eig_values.resize(maximum_eigenvalue, 1e-10);
-	for(int es=0; es<MAX_E; ++es)
+	for(int es=0; es<maximum_eigenvalue; ++es)
+	{
 		for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)
+		{
-			projx[es] += v->pos[0] * Q[es][v->touched];
-			projy[es] += v->pos[1] * Q[es][v->touched];
+			int i= v->touched;
-			projz[es] += v->pos[2] * Q[es][v->touched];
+			proj[es] += Vec3d(v->pos) * Q[es][i];
-			max_eig_values[es] = max(max_eig_values[es], static_cast<float>(abs(Q[es][v->touched])));
+			max_eig_values[es] = max(max_eig_values[es], static_cast<float>(abs(Q[es][i])));
+		}
+	}
-	reconstruct(mani, MAX_E);
+}
-void add_frequency(Manifold& mani, int es, float scale)
+void Harmonics::add_frequency(int es, float scale)
+{
-	if(es<MAX_E)
+	if(es<maximum_eigenvalue)
 		for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)
+		{
-			Vec3f& p = vertex_positions[v->touched];
+			Vec3f& p = v->pos;
-			p[0] += projx[es]* Q[es][v->touched] * scale;
+			v->pos += Vec3f(proj[es] * Q[es][v->touched] * scale);
-			p[1] += projy[es]* Q[es][v->touched] * scale;
-			p[2] += projz[es]* Q[es][v->touched] * scale;
+		}
+}
-void reset_shape(Manifold& mani)
+void Harmonics::reset_shape()
+{
 	for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)
-		vertex_positions[v->touched] = Vec3f(0);
+		v->pos = Vec3f(0);
-}
-void reconstruct(Manifold& mani, int E)
-{
-	reset_shape(mani);
-	for(int es=0;es<=E;++es)
-		add_frequency(mani,es);
+}
-void partial_reconstruct(Manifold& mani, int E0, int E1, float scale)
+void Harmonics::partial_reconstruct(int E0, int E1, float scale)
+{
 	for(int es=E0;es<=E1;++es)
-		add_frequency(mani,es, scale);
+		add_frequency(es, scale);
+}
 template<typename T>
 T& get_CVar_ref(const std::string& s)
+{
 	return *reinterpret_cast<T*> (GetCVarData(s));
+}
-void draw_eigenvalues(Manifold& m)
+void Harmonics::parse_key(unsigned char key)
+{
+		int& display_eigenvalue = get_CVar_ref<int>("display.harmonics.eigenvalue");
+		int& display_diffuse = get_CVar_ref<int>("display.harmonics.diffuse");
+		int& display_highlight = get_CVar_ref<int>("display.harmonics.highlight");
+		switch(key) {
+			case '+':
+				display_eigenvalue = min(display_eigenvalue+1, maximum_eigenvalue);
+				break;
+			case '-':
+				display_eigenvalue = max(display_eigenvalue-1, 0);
+				break;
+			case 'd':
+				display_diffuse = !display_diffuse;
+				break;
+			case 'h':
+				display_highlight = !display_highlight;
+				break;
+		}
+}
+void Harmonics::draw()
+{
 	int& display_eigen = get_CVar_ref<int>("display.harmonics.eigenvalue");
 	int& display_eigen2 = get_CVar_ref<int>("display.harmonics.eigenvalue");
 	int& do_diffuse = get_CVar_ref<int>("display.harmonics.diffuse");
 	int& do_highlight = get_CVar_ref<int>("display.harmonics.highlight");
 	glUseProgram(prog_P0);
 	glUniform1f(glGetUniformLocation(prog_P0,"eig_max"),max_eig_values[display_eigen]);
 	glUniform1f(glGetUniformLocation(prog_P0,"eig_max2"),max_eig_values[display_eigen2]);
 	glUniform1i(glGetUniformLocation(prog_P0,"do_diffuse"),do_diffuse);
    	glUniform1i(glGetUniformLocation(prog_P0,"do_highlight"),do_highlight);
 	GLuint attrib = glGetAttribLocationARB(prog_P0, "eigenvalue");
 	GLuint attrib2 = glGetAttribLocationARB(prog_P0, "eigenvalue2");
 	glFrontFace(GL_CW);
-	for(FaceIter f=m.faces_begin(); f != m.faces_end(); ++f)
+	for(FaceIter f=mani.faces_begin(); f != mani.faces_end(); ++f)
+	{
 		FaceCirculator fc(f);
 		glBegin(GL_TRIANGLES);
 		while(!fc.end())
+		{
 			int i = fc.get_vertex()->touched;
 			glVertexAttrib1f(attrib,Q[display_eigen][i]);
 			glVertexAttrib1f(attrib2,Q[display_eigen2][i]);
 			glNormal3fv(normal(fc.get_vertex()).get());
-			glVertex3fv(vertex_positions[i].get());
+			glVertex3fv(fc.get_vertex()->pos.get());
 			++fc;
+		}
 		glEnd();
+	}
 	glFrontFace(GL_CCW);
 	glUseProgram(0);
+}
-void init_harmonics()
+void Harmonics::init()
+{
+	is_initialized = true;
 	string shader_path = "/Users/jab/GEL/apps/MeshEdit/";
 	GLuint vs = create_glsl_shader(GL_VERTEX_SHADER, shader_path, "tri.vert");
 	GLuint fs = create_glsl_shader(GL_FRAGMENT_SHADER, shader_path, "tri.frag");
 	// Create the program
 	prog_P0 = glCreateProgram();
 	// Attach all shaders
 	if(vs) glAttachShader(prog_P0, vs);
 	if(fs) glAttachShader(prog_P0, fs);
 	// Link the program object and print out the info log
 	glLinkProgram(prog_P0);
 	print_glsl_program_log(prog_P0);
 	// Install program object as part of current state
 	glUseProgram(prog_P0);
 	static CVar<int> display_eigen("display.harmonics.eigenvalue",0);
 	static CVar<int> display_eigen2("display.harmonics.eigenvalue2",0);
 	static CVar<int> do_highlight("display.harmonics.highlight",1);
 	static CVar<int> do_diffuse("display.harmonics.diffuse",1);
+}

Subversion Repositories gelsvn

(root)/trunk/apps/MeshEdit/harmonics.cpp – Rev 388 → 391