Rev 388 | 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 <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;
int MAX_E;
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)
{
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)
{
static bool was_here = false;
if(!was_here)
{
triangulate(mani);
mani.enumerate_vertices();
MAX_E = mani.no_vertices()-1;
E = min(E, MAX_E);
was_here = true;
make_laplace_operator(mani);
vertex_positions.resize(mani.no_vertices());
}
projx.resize(MAX_E);
projy.resize(MAX_E);
projz.resize(MAX_E);
max_eig_values.resize(MAX_E, 1e-10);
for(int es=0; es<MAX_E; ++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];
projz[es] += v->pos[2] * Q[es][v->touched];
max_eig_values[es] = max(max_eig_values[es], static_cast<float>(abs(Q[es][v->touched])));
}
}
}
void add_frequency(Manifold& mani, int es, float scale)
{
if(es<MAX_E)
for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)
{
Vec3f& p = vertex_positions[v->touched];
p[0] += projx[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)
{
for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)
vertex_positions[v->touched] = 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)
{
for(int es=E0;es<=E1;++es)
add_frequency(mani,es, scale);
}
void draw_eigenvalues(Manifold& m)
{
static CVar<int> display_eigen("display.eigenvalue",0);
static CVar<int> display_eigen2("display.eigenvalue2",0);
static CVar<int> do_highlight("display.highlight",1);
static CVar<int> do_diffuse("display.diffuse",1);
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)
{
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());
++fc;
}
glEnd();
}
glFrontFace(GL_CCW);
glUseProgram(0);
}
void init_harmonics(Manifold& mani)
{
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);
analyze_mesh(mani);
reconstruct(mani, MAX_E);
}