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 "harmonics.h"

#if USE_SPARSE_MATRIX

#include <arlgsym.h>

#endif

#include <CGLA/Vec3f.h>

#include <CGLA/Vec3d.h>

#include <CGLA/Mat2x2d.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 "curvature.h"

#include "CSCMatrixBuilder.h"

using namespace CGLA;

using namespace std;

using namespace HMesh;

using namespace Geometry;

using namespace GLGraphics;

using namespace LinAlg;

using namespace CVarUtils;

bool Harmonics::is_initialized=false;

GLuint Harmonics::prog_P0;

namespace

{

	string vss =

	"#version 120\n"

	"#extension GL_EXT_gpu_shader4 : enable\n"

	"	\n"

	"	\n"

	"	attribute float eigenvalue;\n"

	"	attribute float eigenvalue2;\n"

	"	varying vec3 normal;\n"

	"	varying float eig;\n"

	"	varying float eig2;\n"

	"	\n"

	"	void main(void)\n"

	"	{\n"

	"		gl_Position =  ftransform();\n"

	"		normal = normalize(gl_NormalMatrix * gl_Normal);\n"

	"		eig = eigenvalue;\n"

	"		eig2 = eigenvalue2;\n"

	"	}\n";

	string fss = 	

	"#version 120\n"

	"#extension GL_EXT_gpu_shader4 : enable\n"

	"	\n"

	"	varying vec3 normal;\n"

	"	varying float eig;\n"

	"	varying float eig2;\n"

	"	uniform float eig_max;\n"

	"	uniform float eig_max2;\n"

	"	uniform bool do_highlight;\n"

	"	uniform bool do_diffuse;\n"

	"	const vec3 light_dir = vec3(0,0,1);\n"

	"	\n"

	" float basef(float x) {return max(0.0,min(1.0,2.0-4.0*abs(x)));\n}" 

	"	void main()\n"

	"	{\n"

	"		float dot_ln = max(0.0,dot(light_dir, normal));\n"

	"		\n"

	"		float eig_norm = eig/eig_max;\n"

	"		float stripe_signal = 250 * eig_norm;\n"

	//"		vec4 stripe_col = abs(stripe_signal) < 3.14 ? vec4(1,1,0,0) : vec4(.4,.4,.4,0);\n"

	"		vec4 stripe_col = -vec4(.4,.4,.4,0);\n"

	"		\n"

	"       float alpha = (1.0-eig_norm) * 2.0 * 3.1415926;\n"

	"       float offs = 2.0*3.1415/3.0;\n"

	"		gl_FragColor = vec4(0,0,1,0)*basef(eig_norm)+vec4(0,1,0,0)*basef(eig_norm-0.5)+vec4(1,0,0,0)* basef(eig_norm-1.0);\n"

	"		if(do_diffuse)   gl_FragColor *= dot_ln;\n"

	"		if(do_highlight) gl_FragColor += dot_ln*dot_ln*dot_ln*dot_ln*dot_ln*dot_ln*dot_ln*vec4(.5,.5,.5,0);\n"

	"		gl_FragColor -= vec4(.4,.4,.4,.4)*smoothstep(0.2,0.6,cos(stripe_signal));\n"

	"	}\n";

}

#if USE_SPARSE_MATRIX

double calculate_laplace_entry(VertexCirculator vc,Vec3d vertex)

{

	//* this piece of code is taken from the orginal make_laplace_operator...

	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))); // w positive...?

	double a_right = acos(min(1.0, max(-1.0, d_right)));

	double w = 1.0/tan(a_left) + 1.0/tan(a_right);

	return -w;

}

void Harmonics::make_laplace_operator_sparse()

{

	bool is_generalized=true;

	CSCMatrixBuilder<double> mb_M;

	//+GENERALIZED

	CSCMatrixBuilder<double> mb_S;

	// S has the same type as V

	S.Resize(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);

			double a_sum = 0;

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

			{

				int j = vc.get_vertex()->touched;

				double entry = calculate_laplace_entry(vc,vertex);

				if(!is_generalized)

				{

					double area_j = voronoi_area(vc.get_vertex());

					entry /= sqrt(area_i*area_j);

				}

				if(j > i)

					mb_M.insert_entry(j,entry);

				a_sum += entry;

			}

			//cout << a_sum << " ";

			mb_M.insert_entry(i,-a_sum);

			mb_M.next_column();

			if(!is_generalized)

				area_i = 1; // if standard S is an identity matrix;

			mb_S.insert_entry(i,area_i);

			mb_S.next_column_nonsort();

			S[i] = area_i;

		}     

	cout << "solving generalized problem... i = " << mani.no_vertices()<< endl;

	//+STANDARD

	//ArpackPP::ARluSymStdEig<double> dprob(50L, mb_M.get_Matrix(), "SA");

	//+GENERALIZED (shifted inv mode)

	ARluSymGenEig<double> dprob('S',maximum_eigenvalue, mb_M.get_Matrix(), mb_S.get_Matrix(), 0.0);

	//ArpackPP::ARluSymGenEig<double> dprob(number_of_eigenvalues, mb_M.get_Matrix(), mb_S.get_Matrix());

	dprob.FindEigenvectors();

	int conv = dprob.ConvergedEigenvalues();

	cout << conv << " eigenvectors found" << endl;

	Q.Resize(conv, dprob.GetN());

	V.Resize(conv);

	qnorm.Resize(conv);

	for(int i = 0; i < conv; i++)

	{

		V[i] = dprob.Eigenvalue(i);

		qnorm[i] = 0;

		for(int j = 0; j < dprob.GetN(); j++)

		{

			Q[i][j] = dprob.Eigenvector(i,j);

			qnorm[i] += Q[i][j]*Q[i][j];

		}

		cout << "(" << i << ":" << sqrt(V[i]/V[1]

										) << ")" << V[i]/V[1] << "V= "<< V[i] << " exp " << exp(-0.01*V[i]/V[1]) << " Qnorm " << qnorm[i]<< endl;

	}

	cout  << endl;

}

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

{

	assert(is_initialized);

	maximum_eigenvalue = min(size_t(500),mani.no_vertices());

	triangulate(mani);

	mani.enumerate_vertices();

	make_laplace_operator_sparse();

	if(maximum_eigenvalue == -1)

	{

		cout << "not found" << endl;

		return;

	}

	proj.resize(maximum_eigenvalue);

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

	cout << endl << "Proj" << endl;

	for(int es=0; es<maximum_eigenvalue; ++es)  //o(n^2)

	{

		proj[es] = Vec3d(0.0);

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

		{

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

			max_eig_values[es] = max(max_eig_values[es], static_cast<float>(abs(Q[es][v->touched])));     

		}     

	}

	cout << endl;

}

#else

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

	cout << "Projection magnitude" << endl;

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

	{

		proj[es] = Vec3d(0.0);

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

		{

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

			max_eig_values[es] = max(max_eig_values[es], static_cast<float>(abs(Q[es][v->touched])));

		}

	}

}

#endif

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 = Vec3f(proj[es]);

			double Qval = Q[es][v->touched];

			v->pos += p * Qval * 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);

}

void Harmonics::parse_key(unsigned char key)

{

	float& display_t = GetCVarRef<float>("display.harmonics.time");

	int& display_diffuse = GetCVarRef<int>("display.harmonics.diffuse");

	int& display_highlight = GetCVarRef<int>("display.harmonics.highlight");

	switch(key) {

		case '+': 

			display_t = display_t+0.001; 

			break;

		case '-': 

			display_t = display_t-0.001; 

			break;

		case 'd':	

			display_diffuse = !display_diffuse; 

			break;

		case 'h':

			display_highlight = !display_highlight;

			break;			

	}

}

void Harmonics::draw_adf()

{

	float& t = GetCVarRef<float>("display.harmonics.time");

	int& do_diffuse = GetCVarRef<int>("display.harmonics.diffuse");

	int& do_highlight = GetCVarRef<int>("display.harmonics.highlight");

	vector<double> F(mani.no_vertices(),0);

	double F_max = 0;

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

	{

		int v_idx = v->touched;

		for(int e=1;e<V.Length();++e)

			F[v_idx] += sqr(Q[e][v_idx])*exp(-(t)*V[e]/V[1]);

		F_max = max(F[v_idx], F_max);

	}

	cout << "F max" <<  F_max << endl;

	glUseProgram(prog_P0);

	glUniform1f(glGetUniformLocation(prog_P0,"eig_max"),F_max);//2*M_PI);

	glUniform1i(glGetUniformLocation(prog_P0,"do_diffuse"),do_diffuse);

   	glUniform1i(glGetUniformLocation(prog_P0,"do_highlight"),do_highlight);

	GLuint attrib = glGetAttribLocationARB(prog_P0, "eigenvalue");

	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,F[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, vss);

	GLuint fs = create_glsl_shader(GL_FRAGMENT_SHADER, fss);

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

	CreateCVar("display.harmonics.time",0.0f);

	CreateCVar("display.harmonics.eigenvalue2",0);

	CreateCVar("display.harmonics.highlight",1);

	CreateCVar("display.harmonics.diffuse",1);

	CreateCVar("display.harmonics.diffusion_distance_init_vertex",0);

}