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"
 
 
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
{
 
	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"
"	void main()\n"
"	{\n"
"		float dot_ln = dot(light_dir, normal);\n"
"		\n"
"		float eig_norm = eig/eig_max;\n"
"		float stripe_signal = 100 * eig_norm;\n"
"		vec4 stripe_col = abs(stripe_signal) < 3.14 ? vec4(1,1,0,0) : vec4(.1,.1,.1,0);\n"
"		\n"
"		gl_FragColor = eig_norm * vec4(-1,0,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 += stripe_col * smoothstep(0.8,1.0,cos(stripe_signal));\n"
"	}\n";
	
	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)
	{
		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])));
		}
	}
}
 
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]);
			float p0 = p[0];
			float p1 = p[1];
			float p2 = p[2];
			int idx = v->touched;
			double Qval = Q[es][v->touched];
			
			v->pos += p * Qval * scale; // 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, 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);
	
	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 397	Rev 399
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"`	-
33		32
34	`using namespace CGLA;`	33	`using namespace CGLA;`
35	`using namespace std;`	34	`using namespace std;`
36	`using namespace HMesh;`	35	`using namespace HMesh;`
37	`using namespace Geometry;`	36	`using namespace Geometry;`
38	`using namespace GLGraphics;`	37	`using namespace GLGraphics;`
39	`using namespace LinAlg;`	38	`using namespace LinAlg;`
40		39
41	`bool Harmonics::is_initialized=false;`	40	`bool Harmonics::is_initialized=false;`
42	`GLuint Harmonics::prog_P0;`	41	`GLuint Harmonics::prog_P0;`
43		42
44		43
45	`namespace`	44	`namespace`
46	`{`	45	`{`
47		46
48	`string vss =`	47	`string vss =`
49	`"#version 120\n"`	48	`"#version 120\n"`
50	`"#extension GL_EXT_gpu_shader4 : enable\n"`	49	`"#extension GL_EXT_gpu_shader4 : enable\n"`
51	`" \n"`	50	`" \n"`
52	`" \n"`	51	`" \n"`
53	`" attribute float eigenvalue;\n"`	52	`" attribute float eigenvalue;\n"`
54	`" attribute float eigenvalue2;\n"`	53	`" attribute float eigenvalue2;\n"`
55	`" varying vec3 normal;\n"`	54	`" varying vec3 normal;\n"`
56	`" varying float eig;\n"`	55	`" varying float eig;\n"`
57	`" varying float eig2;\n"`	56	`" varying float eig2;\n"`
58	`" \n"`	57	`" \n"`
59	`" void main(void)\n"`	58	`" void main(void)\n"`
60	`" {\n"`	59	`" {\n"`
61	`" gl_Position = ftransform();\n"`	60	`" gl_Position = ftransform();\n"`
62	`" normal = normalize(gl_NormalMatrix * gl_Normal);\n"`	61	`" normal = normalize(gl_NormalMatrix * gl_Normal);\n"`
63	`" eig = eigenvalue;\n"`	62	`" eig = eigenvalue;\n"`
64	`" eig2 = eigenvalue2;\n"`	63	`" eig2 = eigenvalue2;\n"`
65	`" }\n";`	64	`" }\n";`
66		65
67	`string fss =`	66	`string fss =`
68	`"#version 120\n"`	67	`"#version 120\n"`
69	`"#extension GL_EXT_gpu_shader4 : enable\n"`	68	`"#extension GL_EXT_gpu_shader4 : enable\n"`
70	`" \n"`	69	`" \n"`
71	`" varying vec3 normal;\n"`	70	`" varying vec3 normal;\n"`
72	`" varying float eig;\n"`	71	`" varying float eig;\n"`
73	`" varying float eig2;\n"`	72	`" varying float eig2;\n"`
74	`" uniform float eig_max;\n"`	73	`" uniform float eig_max;\n"`
75	`" uniform float eig_max2;\n"`	74	`" uniform float eig_max2;\n"`
76	`" uniform bool do_highlight;\n"`	75	`" uniform bool do_highlight;\n"`
77	`" uniform bool do_diffuse;\n"`	76	`" uniform bool do_diffuse;\n"`
78	`" const vec3 light_dir = vec3(0,0,1);\n"`	77	`" const vec3 light_dir = vec3(0,0,1);\n"`
79	`" \n"`	78	`" \n"`
80	`" void main()\n"`	79	`" void main()\n"`
81	`" {\n"`	80	`" {\n"`
82	`" float dot_ln = dot(light_dir, normal);\n"`	81	`" float dot_ln = dot(light_dir, normal);\n"`
83	`" \n"`	82	`" \n"`
84	`" float eig_norm = eig/eig_max;\n"`	83	`" float eig_norm = eig/eig_max;\n"`
85	`" float stripe_signal = 100 * eig_norm;\n"`	84	`" float stripe_signal = 100 * eig_norm;\n"`
86	`" vec4 stripe_col = abs(stripe_signal) < 3.14 ? vec4(1,1,0,0) : vec4(.1,.1,.1,0);\n"`	85	`" vec4 stripe_col = abs(stripe_signal) < 3.14 ? vec4(1,1,0,0) : vec4(.1,.1,.1,0);\n"`
87	`" \n"`	86	`" \n"`
88	`" gl_FragColor = eig_norm * vec4(-1,0,1,0);\n"`	87	`" gl_FragColor = eig_norm * vec4(-1,0,1,0);\n"`
89	`" if(do_diffuse) gl_FragColor *= dot_ln;\n"`	88	`" if(do_diffuse) gl_FragColor *= dot_ln;\n"`
90	`" if(do_highlight) gl_FragColor += dot_lndot_lndot_lndot_lndot_lndot_lndot_ln*vec4(.5,.5,.5,0);\n"`	89	`" if(do_highlight) gl_FragColor += dot_lndot_lndot_lndot_lndot_lndot_lndot_ln*vec4(.5,.5,.5,0);\n"`
91	`" gl_FragColor += stripe_col * smoothstep(0.8,1.0,cos(stripe_signal));\n"`	90	`" gl_FragColor += stripe_col * smoothstep(0.8,1.0,cos(stripe_signal));\n"`
92	`" }\n";`	91	`" }\n";`
93		92
94	`double voronoi_area(VertexIter v)`	93	`double voronoi_area(VertexIter v)`
95	`{`	94	`{`
96	`double area_mixed = 0;`	95	`double area_mixed = 0;`
97	`//For each triangle T from the 1-ring neighborhood of x`	96	`//For each triangle T from the 1-ring neighborhood of x`
98	`for(VertexCirculator vc(v); !vc.end(); ++vc)`	97	`for(VertexCirculator vc(v); !vc.end(); ++vc)`
99	`{`	98	`{`
100	`FaceIter f = vc.get_face();`	99	`FaceIter f = vc.get_face();`
101	`double f_area = area(f);`	100	`double f_area = area(f);`
102		101
103	`HalfEdgeIter he = vc.get_halfedge();`	102	`HalfEdgeIter he = vc.get_halfedge();`
104	`Vec3d v1(he->vert->pos);`	103	`Vec3d v1(he->vert->pos);`
105	`Vec3d v2(he->next->vert->pos);`	104	`Vec3d v2(he->next->vert->pos);`
106	`Vec3d v0(he->next->next->vert->pos);`	105	`Vec3d v0(he->next->next->vert->pos);`
107		106
108	`double a0 = acos(dot(v1-v0, v2-v0)/(length(v1-v0)*length(v2-v0)));`	107	`double a0 = acos(dot(v1-v0, v2-v0)/(length(v1-v0)*length(v2-v0)));`
109	`double a1 = acos(dot(v2-v1, v0-v1)/(length(v2-v1)*length(v0-v1)));`	108	`double a1 = acos(dot(v2-v1, v0-v1)/(length(v2-v1)*length(v0-v1)));`
110	`double a2 = acos(dot(v0-v2, v1-v2)/(length(v0-v2)*length(v1-v2)));`	109	`double a2 = acos(dot(v0-v2, v1-v2)/(length(v0-v2)*length(v1-v2)));`
111		110
112	`if(a0>(M_PI/2.0) && a1>(M_PI/2.0) && a2>(M_PI/2.0)) // f is non-obtuse`	111	`if(a0>(M_PI/2.0) && a1>(M_PI/2.0) && a2>(M_PI/2.0)) // f is non-obtuse`
113	`{`	112	`{`
114	`// Add Voronoi formula (see Section 3.3)`	113	`// Add Voronoi formula (see Section 3.3)`
115	`area_mixed += (1.0/8) *`	114	`area_mixed += (1.0/8) *`
116	`((1.0/tan(a1)) * sqr_length(v2-v0) +`	115	`((1.0/tan(a1)) * sqr_length(v2-v0) +`
117	`(1.0/tan(a2)) * sqr_length(v1-v0));`	116	`(1.0/tan(a2)) * sqr_length(v1-v0));`
118	`}`	117	`}`
119	`else // Voronoi inappropriate`	118	`else // Voronoi inappropriate`
120	`{`	119	`{`
121	`// Add either area(f)/4 or area(f)/2`	120	`// Add either area(f)/4 or area(f)/2`
122	`if(a0>M_PI/2.0)// the angle of f at x is obtuse`	121	`if(a0>M_PI/2.0)// the angle of f at x is obtuse`
123	`area_mixed += f_area/2;`	122	`area_mixed += f_area/2;`
124	`else`	123	`else`
125	`area_mixed += f_area/4;`	124	`area_mixed += f_area/4;`
126	`}`	125	`}`
127	`}`	126	`}`
128	`return area_mixed;`	127	`return area_mixed;`
129	`}`	128	`}`
130		129
131	`double barycentric_area(VertexIter v)`	130	`double barycentric_area(VertexIter v)`
132	`{`	131	`{`
133	`double barea = 0;`	132	`double barea = 0;`
134	`//For each triangle T from the 1-ring neighborhood of x`	133	`//For each triangle T from the 1-ring neighborhood of x`
135	`for(VertexCirculator vc(v); !vc.end(); ++vc)`	134	`for(VertexCirculator vc(v); !vc.end(); ++vc)`
136	`{`	135	`{`
137	`FaceIter f = vc.get_face();`	136	`FaceIter f = vc.get_face();`
138	`barea += area(f)/3.0;`	137	`barea += area(f)/3.0;`
139	`}`	138	`}`
140	`return barea;`	139	`return barea;`
141	`}`	140	`}`
142		141
143	`}`	142	`}`
144		143
145	`void Harmonics::make_laplace_operator()`	144	`void Harmonics::make_laplace_operator()`
146	`{`	145	`{`
147	`Q.Resize(mani.no_vertices(), mani.no_vertices());`	146	`Q.Resize(mani.no_vertices(), mani.no_vertices());`
148		147
149	`for(VertexIter v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`	148	`for(VertexIter v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`
150	`if(!is_boundary(v))`	149	`if(!is_boundary(v))`
151	`{`	150	`{`
152	`int i = v->touched;`	151	`int i = v->touched;`
153	`double area_i = voronoi_area(v);`	152	`double area_i = voronoi_area(v);`
154	`Vec3d vertex(v->pos);`	153	`Vec3d vertex(v->pos);`
155	`Vec3d curv_normal(0);`	154	`Vec3d curv_normal(0);`
156	`double a_sum = 0;`	155	`double a_sum = 0;`
157	`for(VertexCirculator vc(v); !vc.end(); ++vc)`	156	`for(VertexCirculator vc(v); !vc.end(); ++vc)`
158	`{`	157	`{`
159	`int j = vc.get_vertex()->touched;`	158	`int j = vc.get_vertex()->touched;`
160	`double area_j = voronoi_area(vc.get_vertex());`	159	`double area_j = voronoi_area(vc.get_vertex());`
161	`HalfEdgeIter h = vc.get_halfedge();`	160	`HalfEdgeIter h = vc.get_halfedge();`
162	`Vec3d nbr(h->vert->pos);`	161	`Vec3d nbr(h->vert->pos);`
163	`Vec3d left(h->next->vert->pos);`	162	`Vec3d left(h->next->vert->pos);`
164	`Vec3d right(h->opp->prev->opp->vert->pos);`	163	`Vec3d right(h->opp->prev->opp->vert->pos);`
165		164
166	`double d_left = dot(normalize(nbr-left),`	165	`double d_left = dot(normalize(nbr-left),`
167	`normalize(vertex-left));`	166	`normalize(vertex-left));`
168	`double d_right = dot(normalize(nbr-right),`	167	`double d_right = dot(normalize(nbr-right),`
169	`normalize(vertex-right));`	168	`normalize(vertex-right));`
170	`double a_left = acos(min(1.0, max(-1.0, d_left)));`	169	`double a_left = acos(min(1.0, max(-1.0, d_left)));`
171	`double a_right = acos(min(1.0, max(-1.0, d_right)));`	170	`double a_right = acos(min(1.0, max(-1.0, d_right)));`
172		171
173	`double w = 1.0/tan(a_left) + 1.0/tan(a_right);`	172	`double w = 1.0/tan(a_left) + 1.0/tan(a_right);`
174		173
175	`Q[i][j] = -w/sqrt(area_i*area_j);`	174	`Q[i][j] = -w/sqrt(area_i*area_j);`
176	`//Q[i][j] = -1;`	175	`//Q[i][j] = -1;`
177	`a_sum += Q[i][j];`	176	`a_sum += Q[i][j];`
178	`}`	177	`}`
179	`Q[i][i] = -a_sum;`	178	`Q[i][i] = -a_sum;`
180	`}`	179	`}`
181	`EigenSolutionsSym(Q,V);`	180	`EigenSolutionsSym(Q,V);`
182		181
183	`}`	182	`}`
184		183
185		184
186	`Harmonics::Harmonics(Manifold& _mani):mani(_mani)`	185	`Harmonics::Harmonics(Manifold& _mani):mani(_mani)`
187	`{`	186	`{`
188	`assert(is_initialized);`	187	`assert(is_initialized);`
189		188
190	`triangulate(mani);`	189	`triangulate(mani);`
191	`mani.enumerate_vertices();`	190	`mani.enumerate_vertices();`
192	`maximum_eigenvalue = mani.no_vertices()-1;`	191	`maximum_eigenvalue = mani.no_vertices()-1;`
193	`make_laplace_operator();`	192	`make_laplace_operator();`
194		193
195	`proj.resize(maximum_eigenvalue);`	194	`proj.resize(maximum_eigenvalue);`
196	`max_eig_values.resize(maximum_eigenvalue, 1e-10);`	195	`max_eig_values.resize(maximum_eigenvalue, 1e-10);`
-		196
197	`for(int es=0; es<maximum_eigenvalue; ++es)`	197	`for(int es=0; es<maximum_eigenvalue; ++es)`
198	`{`	198	`{`
-		199	`proj[es] = Vec3d(0.0);`
199	`for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)`	200	`for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)`
200	`{`	201	`{`
201	`int i= v->touched;`	-
-		202
202	`proj[es] += Vec3d(v->pos) * Q[es][i];`	203	`proj[es] += Vec3d(v->pos) * Q[es][v->touched];`
203	`max_eig_values[es] = max(max_eig_values[es], static_cast<float>(abs(Q[es][i])));`	204	`max_eig_values[es] = max(max_eig_values[es], static_cast<float>(abs(Q[es][v->touched])));`
204	`}`	205	`}`
205	`}`	206	`}`
206	`}`	207	`}`
207		208
208	`void Harmonics::add_frequency(int es, float scale)`	209	`void Harmonics::add_frequency(int es, float scale)`
209	`{`	210	`{`
210	`if(es<maximum_eigenvalue)`	211	`if(es<maximum_eigenvalue)`
211	`for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)`	212	`for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)`
212	`{`	213	`{`
213	`Vec3f& p = v->pos;`	214	`Vec3f p = Vec3f(proj[es]);`
-		215	`float p0 = p[0];`
-		216	`float p1 = p[1];`
-		217	`float p2 = p[2];`
-		218	`int idx = v->touched;`
-		219	`double Qval = Q[es][v->touched];`
-		220
214	`v->pos += Vec3f(proj[es] * Q[es][v->touched] * scale);`	221	`v->pos += p * Qval * scale; // Vec3f(proj[es] * Q[es][v->touched] * scale);`
215	`}`	222	`}`
216		-
217	`}`	223	`}`
218		224
219	`void Harmonics::reset_shape()`	225	`void Harmonics::reset_shape()`
220	`{`	226	`{`
221	`for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)`	227	`for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)`
222	`v->pos = Vec3f(0);`	228	`v->pos = Vec3f(0);`
223	`}`	229	`}`
224	`void Harmonics::partial_reconstruct(int E0, int E1, float scale)`	230	`void Harmonics::partial_reconstruct(int E0, int E1, float scale)`
225	`{`	231	`{`
226	`for(int es=E0;es<=E1;++es)`	232	`for(int es=E0;es<=E1;++es)`
227	`add_frequency(es, scale);`	233	`add_frequency(es, scale);`
228	`}`	234	`}`
229		235
230		236
231	`template<typename T>`	237	`template<typename T>`
232	`T& get_CVar_ref(const std::string& s)`	238	`T& get_CVar_ref(const std::string& s)`
233	`{`	239	`{`
234	`return reinterpret_cast<T> (GetCVarData(s));`	240	`return reinterpret_cast<T> (GetCVarData(s));`
235	`}`	241	`}`
236		242
237	`void Harmonics::parse_key(unsigned char key)`	243	`void Harmonics::parse_key(unsigned char key)`
238	`{`	244	`{`
239	`int& display_eigenvalue = get_CVar_ref<int>("display.harmonics.eigenvalue");`	245	`int& display_eigenvalue = get_CVar_ref<int>("display.harmonics.eigenvalue");`
240	`int& display_diffuse = get_CVar_ref<int>("display.harmonics.diffuse");`	246	`int& display_diffuse = get_CVar_ref<int>("display.harmonics.diffuse");`
241	`int& display_highlight = get_CVar_ref<int>("display.harmonics.highlight");`	247	`int& display_highlight = get_CVar_ref<int>("display.harmonics.highlight");`
242	`switch(key) {`	248	`switch(key) {`
243	`case '+':`	249	`case '+':`
244	`display_eigenvalue = min(display_eigenvalue+1, maximum_eigenvalue);`	250	`display_eigenvalue = min(display_eigenvalue+1, maximum_eigenvalue);`
245	`break;`	251	`break;`
246	`case '-':`	252	`case '-':`
247	`display_eigenvalue = max(display_eigenvalue-1, 0);`	253	`display_eigenvalue = max(display_eigenvalue-1, 0);`
248	`break;`	254	`break;`
249	`case 'd':`	255	`case 'd':`
250	`display_diffuse = !display_diffuse;`	256	`display_diffuse = !display_diffuse;`
251	`break;`	257	`break;`
252	`case 'h':`	258	`case 'h':`
253	`display_highlight = !display_highlight;`	259	`display_highlight = !display_highlight;`
254	`break;`	260	`break;`
255	`}`	261	`}`
256		262
257	`}`	263	`}`
258		264
259	`void Harmonics::draw()`	265	`void Harmonics::draw()`
260	`{`	266	`{`
261	`int& display_eigen = get_CVar_ref<int>("display.harmonics.eigenvalue");`	267	`int& display_eigen = get_CVar_ref<int>("display.harmonics.eigenvalue");`
262	`int& display_eigen2 = get_CVar_ref<int>("display.harmonics.eigenvalue");`	268	`int& display_eigen2 = get_CVar_ref<int>("display.harmonics.eigenvalue");`
263	`int& do_diffuse = get_CVar_ref<int>("display.harmonics.diffuse");`	269	`int& do_diffuse = get_CVar_ref<int>("display.harmonics.diffuse");`
264	`int& do_highlight = get_CVar_ref<int>("display.harmonics.highlight");`	270	`int& do_highlight = get_CVar_ref<int>("display.harmonics.highlight");`
265		271
266		272
267	`glUseProgram(prog_P0);`	273	`glUseProgram(prog_P0);`
268	`glUniform1f(glGetUniformLocation(prog_P0,"eig_max"),max_eig_values[display_eigen]);`	274	`glUniform1f(glGetUniformLocation(prog_P0,"eig_max"),max_eig_values[display_eigen]);`
269	`glUniform1f(glGetUniformLocation(prog_P0,"eig_max2"),max_eig_values[display_eigen2]);`	275	`glUniform1f(glGetUniformLocation(prog_P0,"eig_max2"),max_eig_values[display_eigen2]);`
270	`glUniform1i(glGetUniformLocation(prog_P0,"do_diffuse"),do_diffuse);`	276	`glUniform1i(glGetUniformLocation(prog_P0,"do_diffuse"),do_diffuse);`
271	`glUniform1i(glGetUniformLocation(prog_P0,"do_highlight"),do_highlight);`	277	`glUniform1i(glGetUniformLocation(prog_P0,"do_highlight"),do_highlight);`
272	`GLuint attrib = glGetAttribLocationARB(prog_P0, "eigenvalue");`	278	`GLuint attrib = glGetAttribLocationARB(prog_P0, "eigenvalue");`
273	`GLuint attrib2 = glGetAttribLocationARB(prog_P0, "eigenvalue2");`	279	`GLuint attrib2 = glGetAttribLocationARB(prog_P0, "eigenvalue2");`
274		280
275	`glFrontFace(GL_CW);`	281	`glFrontFace(GL_CW);`
276	`for(FaceIter f=mani.faces_begin(); f != mani.faces_end(); ++f)`	282	`for(FaceIter f=mani.faces_begin(); f != mani.faces_end(); ++f)`
277	`{`	283	`{`
278	`FaceCirculator fc(f);`	284	`FaceCirculator fc(f);`
279	`glBegin(GL_TRIANGLES);`	285	`glBegin(GL_TRIANGLES);`
280	`while(!fc.end())`	286	`while(!fc.end())`
281	`{`	287	`{`
282	`int i = fc.get_vertex()->touched;`	288	`int i = fc.get_vertex()->touched;`
283	`glVertexAttrib1f(attrib,Q[display_eigen][i]);`	289	`glVertexAttrib1f(attrib,Q[display_eigen][i]);`
284	`glVertexAttrib1f(attrib2,Q[display_eigen2][i]);`	290	`glVertexAttrib1f(attrib2,Q[display_eigen2][i]);`
285	`glNormal3fv(normal(fc.get_vertex()).get());`	291	`glNormal3fv(normal(fc.get_vertex()).get());`
286	`glVertex3fv(fc.get_vertex()->pos.get());`	292	`glVertex3fv(fc.get_vertex()->pos.get());`
287	`++fc;`	293	`++fc;`
288	`}`	294	`}`
289	`glEnd();`	295	`glEnd();`
290	`}`	296	`}`
291	`glFrontFace(GL_CCW);`	297	`glFrontFace(GL_CCW);`
292	`glUseProgram(0);`	298	`glUseProgram(0);`
293	`}`	299	`}`
294		300
295	`void Harmonics::init()`	301	`void Harmonics::init()`
296	`{`	302	`{`
297	`is_initialized = true;`	303	`is_initialized = true;`
298	`string shader_path = "/Users/jab/GEL/apps/MeshEdit/";`	304	`string shader_path = "/Users/jab/GEL/apps/MeshEdit/";`
299	`GLuint vs = create_glsl_shader(GL_VERTEX_SHADER, vss);`	305	`GLuint vs = create_glsl_shader(GL_VERTEX_SHADER, vss);`
300	`GLuint fs = create_glsl_shader(GL_FRAGMENT_SHADER, fss);`	306	`GLuint fs = create_glsl_shader(GL_FRAGMENT_SHADER, fss);`
301		307
302	`// Create the program`	308	`// Create the program`
303	`prog_P0 = glCreateProgram();`	309	`prog_P0 = glCreateProgram();`
304		310
305	`// Attach all shaders`	311	`// Attach all shaders`
306	`if(vs) glAttachShader(prog_P0, vs);`	312	`if(vs) glAttachShader(prog_P0, vs);`
307	`if(fs) glAttachShader(prog_P0, fs);`	313	`if(fs) glAttachShader(prog_P0, fs);`
308		314
309	`// Link the program object and print out the info log`	315	`// Link the program object and print out the info log`
310	`glLinkProgram(prog_P0);`	316	`glLinkProgram(prog_P0);`
311	`print_glsl_program_log(prog_P0);`	317	`print_glsl_program_log(prog_P0);`
312		318
313	`// Install program object as part of current state`	319	`// Install program object as part of current state`
314	`glUseProgram(0);`	320	`glUseProgram(0);`
315		321
316	`static CVar<int> display_eigen("display.harmonics.eigenvalue",0);`	322	`static CVar<int> display_eigen("display.harmonics.eigenvalue",0);`
317	`static CVar<int> display_eigen2("display.harmonics.eigenvalue2",0);`	323	`static CVar<int> display_eigen2("display.harmonics.eigenvalue2",0);`
318	`static CVar<int> do_highlight("display.harmonics.highlight",1);`	324	`static CVar<int> do_highlight("display.harmonics.highlight",1);`
319	`static CVar<int> do_diffuse("display.harmonics.diffuse",1);`	325	`static CVar<int> do_diffuse("display.harmonics.diffuse",1);`
320		326
321	`}`	327	`}`
322		328

Subversion Repositories gelsvn

(root)/trunk/apps/MeshEdit/harmonics.cpp – Rev 397 → 399