Rev 416 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | RSS feed
/*
* 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);
}