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; 	
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])));
		}
	}
	
	reconstruct(mani, MAX_E);
}
 
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);
}
 
 
template<typename T>
T& get_CVar_ref(const std::string& s)
{
	return *reinterpret_cast<T*> (GetCVarData(s));
}
 
void draw_eigenvalues(Manifold& m)
{
	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)
	{
		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()
{
	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);
 
}
 

Rev 386	Rev 388
1	`/*`	1	`/*`
2	`* harmonics.cpp`	2	`* harmonics.cpp`
3	`* GEL`	3	`* GEL`
4	`*`	4	`*`
5	`* Created by J. Andreas Bærentzen on 01/09/08.`	5	`* Created by J. Andreas Bærentzen on 01/09/08.`
6	`* Copyright 2008 __MyCompanyName__. All rights reserved.`	6	`* Copyright 2008 __MyCompanyName__. All rights reserved.`
7	`*`	7	`*`
8	`*/`	8	`*/`
9		9
10	`#include <iostream>`	10	`#include <iostream>`
11	`#include <CGLA/Vec3f.h>`	11	`#include <CGLA/Vec3f.h>`
12	`#include <CGLA/Vec3d.h>`	12	`#include <CGLA/Vec3d.h>`
13	`#include <LinAlg/Matrix.h>`	13	`#include <LinAlg/Matrix.h>`
14	`#include <LinAlg/Vector.h>`	14	`#include <LinAlg/Vector.h>`
15	`#include <LinAlg/LapackFunc.h>`	15	`#include <LinAlg/LapackFunc.h>`
16		16
17	`#include <GL/glew.h>`	17	`#include <GL/glew.h>`
18	`#include <GLGraphics/glsl_shader.h>`	18	`#include <GLGraphics/glsl_shader.h>`
19		19
20	`#include <HMesh/Manifold.h>`	20	`#include <HMesh/Manifold.h>`
21	`#include <HMesh/VertexCirculator.h>`	21	`#include <HMesh/VertexCirculator.h>`
22	`#include <HMesh/FaceCirculator.h>`	22	`#include <HMesh/FaceCirculator.h>`
23	`#include <HMesh/build_manifold.h>`	23	`#include <HMesh/build_manifold.h>`
24	`#include <HMesh/mesh_optimization.h>`	24	`#include <HMesh/mesh_optimization.h>`
25	`#include <HMesh/triangulate.h>`	25	`#include <HMesh/triangulate.h>`
26	`#include <HMesh/load.h>`	26	`#include <HMesh/load.h>`
27	`#include <HMesh/x3d_save.h>`	27	`#include <HMesh/x3d_save.h>`
28		28
29	`#include <GLConsole/GLConsole.h>`	29	`#include <GLConsole/GLConsole.h>`
30		30
31	`#include "harmonics.h"`	31	`#include "harmonics.h"`
32	`#include "wireframe.h"`	32	`#include "wireframe.h"`
33		33
34	`using namespace CGLA;`	34	`using namespace CGLA;`
35	`using namespace std;`	35	`using namespace std;`
36	`using namespace HMesh;`	36	`using namespace HMesh;`
37	`using namespace Geometry;`	37	`using namespace Geometry;`
38	`using namespace GLGraphics;`	38	`using namespace GLGraphics;`
39	`using namespace LinAlg;`	39	`using namespace LinAlg;`
40		40
41	`int E = 1;`	41	`int E = 1;`
42	`int MAX_E;`	42	`int MAX_E;`
43		43
44		44
45	`namespace`	45	`namespace`
46	`{`	46	`{`
47	`GLuint prog_P0;`	47	`GLuint prog_P0;`
48		48
49	`vector<double> projx;`	49	`vector<double> projx;`
50	`vector<double> projy;`	50	`vector<double> projy;`
51	`vector<double> projz;`	51	`vector<double> projz;`
52		52
53	`CMatrix Q;`	53	`CMatrix Q;`
54	`CVector V;`	54	`CVector V;`
55	`vector<float> max_eig_values;`	55	`vector<float> max_eig_values;`
56	`vector<Vec3f> vertex_positions;`	56	`vector<Vec3f> vertex_positions;`
57		57
58	`double voronoi_area(VertexIter v)`	58	`double voronoi_area(VertexIter v)`
59	`{`	59	`{`
60	`double area_mixed = 0;`	60	`double area_mixed = 0;`
61	`//For each triangle T from the 1-ring neighborhood of x`	61	`//For each triangle T from the 1-ring neighborhood of x`
62	`for(VertexCirculator vc(v); !vc.end(); ++vc)`	62	`for(VertexCirculator vc(v); !vc.end(); ++vc)`
63	`{`	63	`{`
64	`FaceIter f = vc.get_face();`	64	`FaceIter f = vc.get_face();`
65	`double f_area = area(f);`	65	`double f_area = area(f);`
66		66
67	`HalfEdgeIter he = vc.get_halfedge();`	67	`HalfEdgeIter he = vc.get_halfedge();`
68	`Vec3d v1(he->vert->pos);`	68	`Vec3d v1(he->vert->pos);`
69	`Vec3d v2(he->next->vert->pos);`	69	`Vec3d v2(he->next->vert->pos);`
70	`Vec3d v0(he->next->next->vert->pos);`	70	`Vec3d v0(he->next->next->vert->pos);`
71		71
72	`double a0 = acos(dot(v1-v0, v2-v0)/(length(v1-v0)*length(v2-v0)));`	72	`double a0 = acos(dot(v1-v0, v2-v0)/(length(v1-v0)*length(v2-v0)));`
73	`double a1 = acos(dot(v2-v1, v0-v1)/(length(v2-v1)*length(v0-v1)));`	73	`double a1 = acos(dot(v2-v1, v0-v1)/(length(v2-v1)*length(v0-v1)));`
74	`double a2 = acos(dot(v0-v2, v1-v2)/(length(v0-v2)*length(v1-v2)));`	74	`double a2 = acos(dot(v0-v2, v1-v2)/(length(v0-v2)*length(v1-v2)));`
75		75
76	`if(a0>(M_PI/2.0) && a1>(M_PI/2.0) && a2>(M_PI/2.0)) // f is non-obtuse`	76	`if(a0>(M_PI/2.0) && a1>(M_PI/2.0) && a2>(M_PI/2.0)) // f is non-obtuse`
77	`{`	77	`{`
78	`// Add Voronoi formula (see Section 3.3)`	78	`// Add Voronoi formula (see Section 3.3)`
79	`area_mixed += (1.0/8) *`	79	`area_mixed += (1.0/8) *`
80	`((1.0/tan(a1)) * sqr_length(v2-v0) +`	80	`((1.0/tan(a1)) * sqr_length(v2-v0) +`
81	`(1.0/tan(a2)) * sqr_length(v1-v0));`	81	`(1.0/tan(a2)) * sqr_length(v1-v0));`
82	`}`	82	`}`
83	`else // Voronoi inappropriate`	83	`else // Voronoi inappropriate`
84	`{`	84	`{`
85	`// Add either area(f)/4 or area(f)/2`	85	`// Add either area(f)/4 or area(f)/2`
86	`if(a0>M_PI/2.0)// the angle of f at x is obtuse`	86	`if(a0>M_PI/2.0)// the angle of f at x is obtuse`
87	`area_mixed += f_area/2;`	87	`area_mixed += f_area/2;`
88	`else`	88	`else`
89	`area_mixed += f_area/4;`	89	`area_mixed += f_area/4;`
90	`}`	90	`}`
91	`}`	91	`}`
92	`return area_mixed;`	92	`return area_mixed;`
93	`}`	93	`}`
94		94
95	`double barycentric_area(VertexIter v)`	95	`double barycentric_area(VertexIter v)`
96	`{`	96	`{`
97	`double barea = 0;`	97	`double barea = 0;`
98	`//For each triangle T from the 1-ring neighborhood of x`	98	`//For each triangle T from the 1-ring neighborhood of x`
99	`for(VertexCirculator vc(v); !vc.end(); ++vc)`	99	`for(VertexCirculator vc(v); !vc.end(); ++vc)`
100	`{`	100	`{`
101	`FaceIter f = vc.get_face();`	101	`FaceIter f = vc.get_face();`
102	`barea += area(f)/3.0;`	102	`barea += area(f)/3.0;`
103	`}`	103	`}`
104	`return barea;`	104	`return barea;`
105	`}`	105	`}`
106		106
107	`}`	107	`}`
108		108
109	`void make_laplace_operator(Manifold& mani)`	109	`void make_laplace_operator(Manifold& mani)`
110	`{`	110	`{`
111	`Q.Resize(mani.no_vertices(), mani.no_vertices());`	111	`Q.Resize(mani.no_vertices(), mani.no_vertices());`
112		112
113	`for(VertexIter v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`	113	`for(VertexIter v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`
114	`if(!is_boundary(v))`	114	`if(!is_boundary(v))`
115	`{`	115	`{`
116	`int i = v->touched;`	116	`int i = v->touched;`
117	`double area_i = voronoi_area(v);`	117	`double area_i = voronoi_area(v);`
118	`Vec3d vertex(v->pos);`	118	`Vec3d vertex(v->pos);`
119	`Vec3d curv_normal(0);`	119	`Vec3d curv_normal(0);`
120	`double a_sum = 0;`	120	`double a_sum = 0;`
121	`for(VertexCirculator vc(v); !vc.end(); ++vc)`	121	`for(VertexCirculator vc(v); !vc.end(); ++vc)`
122	`{`	122	`{`
123	`int j = vc.get_vertex()->touched;`	123	`int j = vc.get_vertex()->touched;`
124	`double area_j = voronoi_area(vc.get_vertex());`	124	`double area_j = voronoi_area(vc.get_vertex());`
125	`HalfEdgeIter h = vc.get_halfedge();`	125	`HalfEdgeIter h = vc.get_halfedge();`
126	`Vec3d nbr(h->vert->pos);`	126	`Vec3d nbr(h->vert->pos);`
127	`Vec3d left(h->next->vert->pos);`	127	`Vec3d left(h->next->vert->pos);`
128	`Vec3d right(h->opp->prev->opp->vert->pos);`	128	`Vec3d right(h->opp->prev->opp->vert->pos);`
129		129
130	`double d_left = dot(normalize(nbr-left),`	130	`double d_left = dot(normalize(nbr-left),`
131	`normalize(vertex-left));`	131	`normalize(vertex-left));`
132	`double d_right = dot(normalize(nbr-right),`	132	`double d_right = dot(normalize(nbr-right),`
133	`normalize(vertex-right));`	133	`normalize(vertex-right));`
134	`double a_left = acos(min(1.0, max(-1.0, d_left)));`	134	`double a_left = acos(min(1.0, max(-1.0, d_left)));`
135	`double a_right = acos(min(1.0, max(-1.0, d_right)));`	135	`double a_right = acos(min(1.0, max(-1.0, d_right)));`
136		136
137	`double w = 1.0/tan(a_left) + 1.0/tan(a_right);`	137	`double w = 1.0/tan(a_left) + 1.0/tan(a_right);`
138		138
139	`Q[i][j] = -w/sqrt(area_i*area_j);`	139	`Q[i][j] = -w/sqrt(area_i*area_j);`
140	`//Q[i][j] = -1;`	140	`//Q[i][j] = -1;`
141	`a_sum += Q[i][j];`	141	`a_sum += Q[i][j];`
142	`}`	142	`}`
143	`Q[i][i] = -a_sum;`	143	`Q[i][i] = -a_sum;`
144	`}`	144	`}`
145	`EigenSolutionsSym(Q,V);`	145	`EigenSolutionsSym(Q,V);`
146		146
147	`}`	147	`}`
148		148
149		149
150	`void analyze_mesh(Manifold& mani)`	150	`void analyze_mesh(Manifold& mani)`
151	`{`	151	`{`
152	`static bool was_here = false;`	152	`static bool was_here = false;`
153	`if(!was_here)`	153	`if(!was_here)`
154	`{`	154	`{`
155	`triangulate(mani);`	155	`triangulate(mani);`
156	`mani.enumerate_vertices();`	156	`mani.enumerate_vertices();`
157	`MAX_E = mani.no_vertices()-1;`	157	`MAX_E = mani.no_vertices()-1;`
158	`E = min(E, MAX_E);`	158	`E = min(E, MAX_E);`
159	`was_here = true;`	159	`was_here = true;`
160	`make_laplace_operator(mani);`	160	`make_laplace_operator(mani);`
161	`vertex_positions.resize(mani.no_vertices());`	161	`vertex_positions.resize(mani.no_vertices());`
162	`}`	162	`}`
163	`projx.resize(MAX_E);`	163	`projx.resize(MAX_E);`
164	`projy.resize(MAX_E);`	164	`projy.resize(MAX_E);`
165	`projz.resize(MAX_E);`	165	`projz.resize(MAX_E);`
166	`max_eig_values.resize(MAX_E, 1e-10);`	166	`max_eig_values.resize(MAX_E, 1e-10);`
167	`for(int es=0; es<MAX_E; ++es)`	167	`for(int es=0; es<MAX_E; ++es)`
168	`{`	168	`{`
169	`for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)`	169	`for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)`
170	`{`	170	`{`
171	`projx[es] += v->pos[0] * Q[es][v->touched];`	171	`projx[es] += v->pos[0] * Q[es][v->touched];`
172	`projy[es] += v->pos[1] * Q[es][v->touched];`	172	`projy[es] += v->pos[1] * Q[es][v->touched];`
173	`projz[es] += v->pos[2] * Q[es][v->touched];`	173	`projz[es] += v->pos[2] * Q[es][v->touched];`
174	`max_eig_values[es] = max(max_eig_values[es], static_cast<float>(abs(Q[es][v->touched])));`	174	`max_eig_values[es] = max(max_eig_values[es], static_cast<float>(abs(Q[es][v->touched])));`
175	`}`	175	`}`
176	`}`	176	`}`
-		177
-		178	`reconstruct(mani, MAX_E);`
177	`}`	179	`}`
178		180
179	`void add_frequency(Manifold& mani, int es, float scale)`	181	`void add_frequency(Manifold& mani, int es, float scale)`
180	`{`	182	`{`
181	`if(es<MAX_E)`	183	`if(es<MAX_E)`
182	`for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)`	184	`for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)`
183	`{`	185	`{`
184	`Vec3f& p = vertex_positions[v->touched];`	186	`Vec3f& p = vertex_positions[v->touched];`
185	`p[0] += projx[es]* Q[es][v->touched] * scale;`	187	`p[0] += projx[es]* Q[es][v->touched] * scale;`
186	`p[1] += projy[es]* Q[es][v->touched] * scale;`	188	`p[1] += projy[es]* Q[es][v->touched] * scale;`
187	`p[2] += projz[es]* Q[es][v->touched] * scale;`	189	`p[2] += projz[es]* Q[es][v->touched] * scale;`
188	`}`	190	`}`
189		191
190	`}`	192	`}`
191		193
192	`void reset_shape(Manifold& mani)`	194	`void reset_shape(Manifold& mani)`
193	`{`	195	`{`
194	`for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)`	196	`for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)`
195	`vertex_positions[v->touched] = Vec3f(0);`	197	`vertex_positions[v->touched] = Vec3f(0);`
196	`}`	198	`}`
197		199
198	`void reconstruct(Manifold& mani, int E)`	200	`void reconstruct(Manifold& mani, int E)`
199	`{`	201	`{`
200	`reset_shape(mani);`	202	`reset_shape(mani);`
201	`for(int es=0;es<=E;++es)`	203	`for(int es=0;es<=E;++es)`
202	`add_frequency(mani,es);`	204	`add_frequency(mani,es);`
203	`}`	205	`}`
204		206
205	`void partial_reconstruct(Manifold& mani, int E0, int E1, float scale)`	207	`void partial_reconstruct(Manifold& mani, int E0, int E1, float scale)`
206	`{`	208	`{`
207	`for(int es=E0;es<=E1;++es)`	209	`for(int es=E0;es<=E1;++es)`
208	`add_frequency(mani,es, scale);`	210	`add_frequency(mani,es, scale);`
209	`}`	211	`}`
210		212
211		213
-		214	`template<typename T>`
-		215	`T& get_CVar_ref(const std::string& s)`
-		216	`{`
-		217	`return reinterpret_cast<T> (GetCVarData(s));`
-		218	`}`
212		219
213	`void draw_eigenvalues(Manifold& m)`	220	`void draw_eigenvalues(Manifold& m)`
214	`{`	221	`{`
215	`static CVar<int> display_eigen("display.eigenvalue",0);`	222	`int& display_eigen = get_CVar_ref<int>("display.harmonics.eigenvalue");`
216	`static CVar<int> display_eigen2("display.eigenvalue2",0);`	223	`int& display_eigen2 = get_CVar_ref<int>("display.harmonics.eigenvalue");`
217	`static CVar<int> do_highlight("display.highlight",1);`	224	`int& do_diffuse = get_CVar_ref<int>("display.harmonics.diffuse");`
218	`static CVar<int> do_diffuse("display.diffuse",1);`	225	`int& do_highlight = get_CVar_ref<int>("display.harmonics.highlight");`
-		226
219		227
220	`glUseProgram(prog_P0);`	228	`glUseProgram(prog_P0);`
221	`glUniform1f(glGetUniformLocation(prog_P0,"eig_max"),max_eig_values[display_eigen]);`	229	`glUniform1f(glGetUniformLocation(prog_P0,"eig_max"),max_eig_values[display_eigen]);`
222	`glUniform1f(glGetUniformLocation(prog_P0,"eig_max2"),max_eig_values[display_eigen2]);`	230	`glUniform1f(glGetUniformLocation(prog_P0,"eig_max2"),max_eig_values[display_eigen2]);`
223	`glUniform1i(glGetUniformLocation(prog_P0,"do_diffuse"),do_diffuse);`	231	`glUniform1i(glGetUniformLocation(prog_P0,"do_diffuse"),do_diffuse);`
224	`glUniform1i(glGetUniformLocation(prog_P0,"do_highlight"),do_highlight);`	232	`glUniform1i(glGetUniformLocation(prog_P0,"do_highlight"),do_highlight);`
225	`GLuint attrib = glGetAttribLocationARB(prog_P0, "eigenvalue");`	233	`GLuint attrib = glGetAttribLocationARB(prog_P0, "eigenvalue");`
226	`GLuint attrib2 = glGetAttribLocationARB(prog_P0, "eigenvalue2");`	234	`GLuint attrib2 = glGetAttribLocationARB(prog_P0, "eigenvalue2");`
227		235
228	`glFrontFace(GL_CW);`	236	`glFrontFace(GL_CW);`
229	`for(FaceIter f=m.faces_begin(); f != m.faces_end(); ++f)`	237	`for(FaceIter f=m.faces_begin(); f != m.faces_end(); ++f)`
230	`{`	238	`{`
231	`FaceCirculator fc(f);`	239	`FaceCirculator fc(f);`
232	`glBegin(GL_TRIANGLES);`	240	`glBegin(GL_TRIANGLES);`
233	`while(!fc.end())`	241	`while(!fc.end())`
234	`{`	242	`{`
235	`int i = fc.get_vertex()->touched;`	243	`int i = fc.get_vertex()->touched;`
236	`glVertexAttrib1f(attrib,Q[display_eigen][i]);`	244	`glVertexAttrib1f(attrib,Q[display_eigen][i]);`
237	`glVertexAttrib1f(attrib2,Q[display_eigen2][i]);`	245	`glVertexAttrib1f(attrib2,Q[display_eigen2][i]);`
238	`glNormal3fv(normal(fc.get_vertex()).get());`	246	`glNormal3fv(normal(fc.get_vertex()).get());`
239	`glVertex3fv(vertex_positions[i].get());`	247	`glVertex3fv(vertex_positions[i].get());`
240	`++fc;`	248	`++fc;`
241	`}`	249	`}`
242	`glEnd();`	250	`glEnd();`
243	`}`	251	`}`
244	`glFrontFace(GL_CCW);`	252	`glFrontFace(GL_CCW);`
245	`glUseProgram(0);`	253	`glUseProgram(0);`
246	`}`	254	`}`
247		255
248	`void init_harmonics(Manifold& mani)`	256	`void init_harmonics()`
249	`{`	257	`{`
250	`string shader_path = "/Users/jab/GEL/apps/MeshEdit/";`	258	`string shader_path = "/Users/jab/GEL/apps/MeshEdit/";`
251	`GLuint vs = create_glsl_shader(GL_VERTEX_SHADER, shader_path, "tri.vert");`	259	`GLuint vs = create_glsl_shader(GL_VERTEX_SHADER, shader_path, "tri.vert");`
252	`GLuint fs = create_glsl_shader(GL_FRAGMENT_SHADER, shader_path, "tri.frag");`	260	`GLuint fs = create_glsl_shader(GL_FRAGMENT_SHADER, shader_path, "tri.frag");`
253		261
254	`// Create the program`	262	`// Create the program`
255	`prog_P0 = glCreateProgram();`	263	`prog_P0 = glCreateProgram();`
256		264
257	`// Attach all shaders`	265	`// Attach all shaders`
258	`if(vs) glAttachShader(prog_P0, vs);`	266	`if(vs) glAttachShader(prog_P0, vs);`
259	`if(fs) glAttachShader(prog_P0, fs);`	267	`if(fs) glAttachShader(prog_P0, fs);`
260		268
261	`// Link the program object and print out the info log`	269	`// Link the program object and print out the info log`
262	`glLinkProgram(prog_P0);`	270	`glLinkProgram(prog_P0);`
263	`print_glsl_program_log(prog_P0);`	271	`print_glsl_program_log(prog_P0);`
264		272
265	`// Install program object as part of current state`	273	`// Install program object as part of current state`
266	`glUseProgram(prog_P0);`	274	`glUseProgram(prog_P0);`
267		275
268	`analyze_mesh(mani);`	276	`static CVar<int> display_eigen("display.harmonics.eigenvalue",0);`
-		277	`static CVar<int> display_eigen2("display.harmonics.eigenvalue2",0);`
-		278	`static CVar<int> do_highlight("display.harmonics.highlight",1);`
269	`reconstruct(mani, MAX_E);`	279	`static CVar<int> do_diffuse("display.harmonics.diffuse",1);`
270		280
271	`}`	281	`}`
272		282

Subversion Repositories gelsvn

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