Subversion Repositories gelsvn

/*

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

bool Harmonics::is_initialized=false;

GLuint Harmonics::prog_P0;

namespace

{

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

}

Harmonics::Harmonics(Manifold& _mani):mani(_mani)

{

	assert(is_initialized);

	triangulate(mani);

	mani.enumerate_vertices();

	maximum_eigenvalue = mani.no_vertices()-1;

	make_laplace_operator();

	proj.resize(maximum_eigenvalue);

	max_eig_values.resize(maximum_eigenvalue, 1e-10);

	for(int es=0; es<maximum_eigenvalue; ++es)

	{

		for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)

		{

			int i= 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][i])));

		}

	}

}

void Harmonics::add_frequency(int es, float scale)

{

	if(es<maximum_eigenvalue)

		for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)

		{

			Vec3f& p = v->pos; 

			v->pos += Vec3f(proj[es] * Q[es][v->touched] * scale);

		}

}

void Harmonics::reset_shape()

{

	for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)

		v->pos = Vec3f(0);	

}

void Harmonics::partial_reconstruct(int E0, int E1, float scale)

{

	for(int es=E0;es<=E1;++es)

		add_frequency(es, scale);

}

template<typename T>

T& get_CVar_ref(const std::string& s)

{

	return *reinterpret_cast<T*> (GetCVarData(s));

}

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=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(fc.get_vertex()->pos.get());

			++fc;

		}

		glEnd();

	}

	glFrontFace(GL_CCW);

	glUseProgram(0);

}

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

}