WebSVN – gelsvn – Diff – /branches/cpp11-devel/src/HMesh/harmonics.cpp


/*
 *  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/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/mesh_optimization.h>
#include <HMesh/curvature.h>
#include <HMesh/triangulate.h>
#include <HMesh/load.h>
#include <HMesh/x3d_save.h>
 
using namespace CGLA;
using namespace std;
using namespace HMesh;
using namespace Geometry;
using namespace LinAlg;
 
 
void Harmonics::make_laplace_operator()
{
	Q.Resize(mani.no_vertices(), mani.no_vertices());
	
	for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)
		if(!boundary(mani, *v)){
			int i = vtouched[*v];
			double area_i = mixed_area(mani, *v);
			Vec3d vertex(mani.pos(*v));
			Vec3d curv_normal(0);
			double a_sum = 0;
            for(Walker wv = mani.walker(*v); !wv.full_circle(); wv = wv.circulate_vertex_cw())
			{
				int j = vtouched[wv.vertex()];
				double area_j = mixed_area(mani, wv.vertex());
				
				Vec3d nbr(mani.pos(wv.vertex()));
				Vec3d left(mani.pos(wv.next().vertex()));
				Vec3d right(mani.pos(wv.opp().prev().opp().vertex()));
				
				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);
}
 
vector<Vec3d> Harmonics::analyze_signal(const VertexAttributeVector<Vec3d>& sig) {
    vector<Vec3d> sig_proj(maximum_eigenvalue, Vec3d(0));
    
    for(int es=0; es<maximum_eigenvalue; ++es)
		for(VertexID v: mani.vertices())
			sig_proj[es] +=  sig[v] * Q[es][vtouched[v]];
    
    return sig_proj;
}
 
VertexAttributeVector<Vec3d> Harmonics::reconstruct_signal(const vector<Vec3d>& sig_proj, int max_es) {
    VertexAttributeVector<Vec3d> sig(mani.allocated_vertices(), Vec3d(0));
    for(int es=0; es<max_es; ++es)
        for(VertexID v: mani.vertices())
            sig[v] += sig_proj[es] * Q[es][vtouched[v]];
    return sig;
}
 
 
 
Harmonics::Harmonics(HMesh::Manifold& _mani):mani(_mani), vtouched(_mani.allocated_vertices(), 0)
{
  
    int i = 0;
    for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v, ++i)
        vtouched[*v] = i;
	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(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)
		{
			proj[es] +=  Vec3d(mani.pos(*v)) * Q[es][vtouched[*v]];
			max_eig_values[es] = max(max_eig_values[es], static_cast<float>(abs(Q[es][vtouched[*v]])));
		}
	}
}
 
void Harmonics::add_frequency(int es, float scale)
{
	if(es<maximum_eigenvalue)
		for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v){
			Vec3d p = Vec3d(proj[es]);
			double Qval = Q[es][vtouched[*v]];
			
			mani.pos(*v) += p * Qval * scale;
		}
}
 
void Harmonics::reset_shape()
{
	for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)
		mani.pos(*v) = Vec3d(0);
}
void Harmonics::partial_reconstruct(int E0, int E1, float scale)
{
	for(int es=E0;es<=E1;++es)
		add_frequency(es, scale);
}
 
 
double Harmonics::compute_adf(HMesh::VertexAttributeVector<double>& F, double t, double fiedler_boost)
{
	double F_max = 0;
    for(VertexID vid : mani.vertices())
        F[vid]=0;
    
	for(VertexID vid : mani.vertices()){
		for(int e = 1; e < V.Length(); ++e)
			F[vid] += sqr(Q[e][vtouched[vid]]) * exp(-t*V[e]/V[1]);
            if(fiedler_boost>0)
                for(int e = 1; e < V.Length(); ++e)
                    F[vid] += sqr(Q[e][vtouched[vid]]) * exp(-t) * fiedler_boost;
        
		F_max = max(F[vid], F_max);
	}
    return F_max;
}
 
double Harmonics::compute_esum(HMesh::VertexAttributeVector<double>& F, int e0, int e1)
{
	double F_max = 0;
    for(VertexID vid : mani.vertices())
        F[vid]=0;
    
	for(VertexID vid : mani.vertices()){
		for(int e = e0; e < e1; ++e)
			F[vid] += sqr(Q[e][vtouched[vid]]);
        
		F_max = max(F[vid], F_max);
	}
    return F_max;
}
 

Rev 663	Rev 664
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		11
12	`#include "harmonics.h"`	12	`#include "harmonics.h"`
13		13
14	`#if USE_SPARSE_MATRIX`	14	`#if USE_SPARSE_MATRIX`
15	`#include <arlgsym.h>`	15	`#include <arlgsym.h>`
16	`#endif`	16	`#endif`
17		17
18	`#include <CGLA/Vec3d.h>`	18	`#include <CGLA/Vec3d.h>`
19	`#include <CGLA/Mat2x2d.h>`	19	`#include <CGLA/Mat2x2d.h>`
20	`#include <LinAlg/Matrix.h>`	20	`#include <LinAlg/Matrix.h>`
21	`#include <LinAlg/Vector.h>`	21	`#include <LinAlg/Vector.h>`
22	`#include <LinAlg/LapackFunc.h>`	22	`#include <LinAlg/LapackFunc.h>`
23		23
24	`#include <GL/glew.h>`	24	`#include <GL/glew.h>`
25	`#include <GLGraphics/glsl_shader.h>`	25	`#include <GLGraphics/glsl_shader.h>`
26		26
27	`#include <HMesh/mesh_optimization.h>`	27	`#include <HMesh/mesh_optimization.h>`
28	`#include <HMesh/curvature.h>`	28	`#include <HMesh/curvature.h>`
29	`#include <HMesh/triangulate.h>`	29	`#include <HMesh/triangulate.h>`
30	`#include <HMesh/load.h>`	30	`#include <HMesh/load.h>`
31	`#include <HMesh/x3d_save.h>`	31	`#include <HMesh/x3d_save.h>`
32		32
33	`using namespace CGLA;`	33	`using namespace CGLA;`
34	`using namespace std;`	34	`using namespace std;`
35	`using namespace HMesh;`	35	`using namespace HMesh;`
36	`using namespace Geometry;`	36	`using namespace Geometry;`
37	`using namespace LinAlg;`	37	`using namespace LinAlg;`
38		38
39		39
40	`void Harmonics::make_laplace_operator()`	40	`void Harmonics::make_laplace_operator()`
41	`{`	41	`{`
42	`Q.Resize(mani.no_vertices(), mani.no_vertices());`	42	`Q.Resize(mani.no_vertices(), mani.no_vertices());`
43		43
44	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`	44	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`
45	`if(!boundary(mani, *v)){`	45	`if(!boundary(mani, *v)){`
46	`int i = vtouched[*v];`	46	`int i = vtouched[*v];`
47	`double area_i = mixed_area(mani, *v);`	47	`double area_i = mixed_area(mani, *v);`
48	`Vec3d vertex(mani.pos(*v));`	48	`Vec3d vertex(mani.pos(*v));`
49	`Vec3d curv_normal(0);`	49	`Vec3d curv_normal(0);`
50	`double a_sum = 0;`	50	`double a_sum = 0;`
51	`for(Walker wv = mani.walker(*v); !wv.full_circle(); wv = wv.circulate_vertex_cw())`	51	`for(Walker wv = mani.walker(*v); !wv.full_circle(); wv = wv.circulate_vertex_cw())`
52	`{`	52	`{`
53	`int j = vtouched[wv.vertex()];`	53	`int j = vtouched[wv.vertex()];`
54	`double area_j = mixed_area(mani, wv.vertex());`	54	`double area_j = mixed_area(mani, wv.vertex());`
55		55
56	`Vec3d nbr(mani.pos(wv.vertex()));`	56	`Vec3d nbr(mani.pos(wv.vertex()));`
57	`Vec3d left(mani.pos(wv.next().vertex()));`	57	`Vec3d left(mani.pos(wv.next().vertex()));`
58	`Vec3d right(mani.pos(wv.opp().prev().opp().vertex()));`	58	`Vec3d right(mani.pos(wv.opp().prev().opp().vertex()));`
59		59
60	`double d_left = dot(normalize(nbr-left),`	60	`double d_left = dot(normalize(nbr-left),`
61	`normalize(vertex-left));`	61	`normalize(vertex-left));`
62	`double d_right = dot(normalize(nbr-right),`	62	`double d_right = dot(normalize(nbr-right),`
63	`normalize(vertex-right));`	63	`normalize(vertex-right));`
64	`double a_left = acos(min(1.0, max(-1.0, d_left)));`	64	`double a_left = acos(min(1.0, max(-1.0, d_left)));`
65	`double a_right = acos(min(1.0, max(-1.0, d_right)));`	65	`double a_right = acos(min(1.0, max(-1.0, d_right)));`
66		66
67	`double w = 1.0/tan(a_left) + 1.0/tan(a_right);`	67	`double w = 1.0/tan(a_left) + 1.0/tan(a_right);`
68		68
69	`Q[i][j] = -w/sqrt(area_i*area_j);`	69	`Q[i][j] = -w/sqrt(area_i*area_j);`
70	`//Q[i][j] = -1;`	70	`//Q[i][j] = -1;`
71	`a_sum += Q[i][j];`	71	`a_sum += Q[i][j];`
72	`}`	72	`}`
73	`Q[i][i] = -a_sum;`	73	`Q[i][i] = -a_sum;`
74	`}`	74	`}`
75	`EigenSolutionsSym(Q,V);`	75	`EigenSolutionsSym(Q,V);`
76	`}`	76	`}`
77		77
-		78	`vector<Vec3d> Harmonics::analyze_signal(const VertexAttributeVector<Vec3d>& sig) {`
-		79	`vector<Vec3d> sig_proj(maximum_eigenvalue, Vec3d(0));`
-		80
-		81	`for(int es=0; es<maximum_eigenvalue; ++es)`
-		82	`for(VertexID v: mani.vertices())`
-		83	`sig_proj[es] += sig[v] * Q[es][vtouched[v]];`
-		84
-		85	`return sig_proj;`
-		86	`}`
-		87
-		88	`VertexAttributeVector<Vec3d> Harmonics::reconstruct_signal(const vector<Vec3d>& sig_proj, int max_es) {`
-		89	`VertexAttributeVector<Vec3d> sig(mani.allocated_vertices(), Vec3d(0));`
-		90	`for(int es=0; es<max_es; ++es)`
-		91	`for(VertexID v: mani.vertices())`
-		92	`sig[v] += sig_proj[es] * Q[es][vtouched[v]];`
-		93	`return sig;`
-		94	`}`
-		95
-		96
-		97
78	`Harmonics::Harmonics(HMesh::Manifold& _mani):mani(_mani), vtouched(_mani.allocated_vertices(), 0)`	98	`Harmonics::Harmonics(HMesh::Manifold& _mani):mani(_mani), vtouched(_mani.allocated_vertices(), 0)`
79	`{`	99	`{`
80		100
81	`int i = 0;`	101	`int i = 0;`
82	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v, ++i)`	102	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v, ++i)`
83	`vtouched[*v] = i;`	103	`vtouched[*v] = i;`
84	`maximum_eigenvalue = mani.no_vertices()-1;`	104	`maximum_eigenvalue = mani.no_vertices()-1;`
85	`make_laplace_operator();`	105	`make_laplace_operator();`
86		106
87	`proj.resize(maximum_eigenvalue);`	107	`proj.resize(maximum_eigenvalue);`
88	`max_eig_values.resize(maximum_eigenvalue, 1e-10f);`	108	`max_eig_values.resize(maximum_eigenvalue, 1e-10f);`
89		109
90	`cout << "Projection magnitude" << endl;`	110	`cout << "Projection magnitude" << endl;`
91	`for(int es=0; es<maximum_eigenvalue; ++es)`	111	`for(int es=0; es<maximum_eigenvalue; ++es)`
92	`{`	112	`{`
93	`proj[es] = Vec3d(0.0);`	113	`proj[es] = Vec3d(0.0);`
94	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`	114	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`
95	`{`	115	`{`
96	`proj[es] += Vec3d(mani.pos(v)) Q[es][vtouched[*v]];`	116	`proj[es] += Vec3d(mani.pos(v)) Q[es][vtouched[*v]];`
97	`max_eig_values[es] = max(max_eig_values[es], static_cast<float>(abs(Q[es][vtouched[*v]])));`	117	`max_eig_values[es] = max(max_eig_values[es], static_cast<float>(abs(Q[es][vtouched[*v]])));`
98	`}`	118	`}`
99	`}`	119	`}`
100	`}`	120	`}`
101		121
102	`void Harmonics::add_frequency(int es, float scale)`	122	`void Harmonics::add_frequency(int es, float scale)`
103	`{`	123	`{`
104	`if(es<maximum_eigenvalue)`	124	`if(es<maximum_eigenvalue)`
105	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v){`	125	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v){`
106	`Vec3d p = Vec3d(proj[es]);`	126	`Vec3d p = Vec3d(proj[es]);`
107	`double Qval = Q[es][vtouched[*v]];`	127	`double Qval = Q[es][vtouched[*v]];`
108		128
109	`mani.pos(v) += p Qval * scale;`	129	`mani.pos(v) += p Qval * scale;`
110	`}`	130	`}`
111	`}`	131	`}`
112		132
113	`void Harmonics::reset_shape()`	133	`void Harmonics::reset_shape()`
114	`{`	134	`{`
115	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`	135	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`
116	`mani.pos(*v) = Vec3d(0);`	136	`mani.pos(*v) = Vec3d(0);`
117	`}`	137	`}`
118	`void Harmonics::partial_reconstruct(int E0, int E1, float scale)`	138	`void Harmonics::partial_reconstruct(int E0, int E1, float scale)`
119	`{`	139	`{`
120	`for(int es=E0;es<=E1;++es)`	140	`for(int es=E0;es<=E1;++es)`
121	`add_frequency(es, scale);`	141	`add_frequency(es, scale);`
122	`}`	142	`}`
123		143
124		144
125	`double Harmonics::compute_adf(HMesh::VertexAttributeVector<double>& F, double t, double fiedler_boost)`	145	`double Harmonics::compute_adf(HMesh::VertexAttributeVector<double>& F, double t, double fiedler_boost)`
126	`{`	146	`{`
127	`double F_max = 0;`	147	`double F_max = 0;`
128	`for(VertexID vid : mani.vertices())`	148	`for(VertexID vid : mani.vertices())`
129	`F[vid]=0;`	149	`F[vid]=0;`
130		150
131	`for(VertexID vid : mani.vertices()){`	151	`for(VertexID vid : mani.vertices()){`
132	`for(int e = 1; e < V.Length(); ++e)`	152	`for(int e = 1; e < V.Length(); ++e)`
133	`F[vid] += sqr(Q[e][vtouched[vid]]) * exp(-t*V[e]/V[1]);`	153	`F[vid] += sqr(Q[e][vtouched[vid]]) * exp(-t*V[e]/V[1]);`
134	`if(fiedler_boost>0)`	154	`if(fiedler_boost>0)`
135	`for(int e = 1; e < V.Length(); ++e)`	155	`for(int e = 1; e < V.Length(); ++e)`
136	`F[vid] += sqr(Q[e][vtouched[vid]]) * exp(-t) * fiedler_boost;`	156	`F[vid] += sqr(Q[e][vtouched[vid]]) * exp(-t) * fiedler_boost;`
137		157
138	`F_max = max(F[vid], F_max);`	158	`F_max = max(F[vid], F_max);`
139	`}`	159	`}`
140	`return F_max;`	160	`return F_max;`
141	`}`	161	`}`
142		162
143	`double Harmonics::compute_esum(HMesh::VertexAttributeVector<double>& F, int e0, int e1)`	163	`double Harmonics::compute_esum(HMesh::VertexAttributeVector<double>& F, int e0, int e1)`
144	`{`	164	`{`
145	`double F_max = 0;`	165	`double F_max = 0;`
146	`for(VertexID vid : mani.vertices())`	166	`for(VertexID vid : mani.vertices())`
147	`F[vid]=0;`	167	`F[vid]=0;`
148		168
149	`for(VertexID vid : mani.vertices()){`	169	`for(VertexID vid : mani.vertices()){`
150	`for(int e = e0; e < e1; ++e)`	170	`for(int e = e0; e < e1; ++e)`
151	`F[vid] += sqr(Q[e][vtouched[vid]]);`	171	`F[vid] += sqr(Q[e][vtouched[vid]]);`
152		172
153	`F_max = max(F[vid], F_max);`	173	`F_max = max(F[vid], F_max);`
154	`}`	174	`}`
155	`return F_max;`	175	`return F_max;`
156	`}`	176	`}`
157		177

Subversion Repositories gelsvn

(root)/branches/cpp11-devel/src/HMesh/harmonics.cpp – Rev 663 → 664