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);
}
 
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 662	Rev 663
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	`Harmonics::Harmonics(HMesh::Manifold& _mani):mani(_mani), vtouched(_mani.allocated_vertices(), 0)`	78	`Harmonics::Harmonics(HMesh::Manifold& _mani):mani(_mani), vtouched(_mani.allocated_vertices(), 0)`
79	`{`	79	`{`
80		80
81	`int i = 0;`	81	`int i = 0;`
82	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v, ++i)`	82	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v, ++i)`
83	`vtouched[*v] = i;`	83	`vtouched[*v] = i;`
84	`maximum_eigenvalue = mani.no_vertices()-1;`	84	`maximum_eigenvalue = mani.no_vertices()-1;`
85	`make_laplace_operator();`	85	`make_laplace_operator();`
86		86
87	`proj.resize(maximum_eigenvalue);`	87	`proj.resize(maximum_eigenvalue);`
88	`max_eig_values.resize(maximum_eigenvalue, 1e-10f);`	88	`max_eig_values.resize(maximum_eigenvalue, 1e-10f);`
89		89
90	`cout << "Projection magnitude" << endl;`	90	`cout << "Projection magnitude" << endl;`
91	`for(int es=0; es<maximum_eigenvalue; ++es)`	91	`for(int es=0; es<maximum_eigenvalue; ++es)`
92	`{`	92	`{`
93	`proj[es] = Vec3d(0.0);`	93	`proj[es] = Vec3d(0.0);`
94	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`	94	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`
95	`{`	95	`{`
96	`proj[es] += Vec3d(mani.pos(v)) Q[es][vtouched[*v]];`	96	`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]])));`	97	`max_eig_values[es] = max(max_eig_values[es], static_cast<float>(abs(Q[es][vtouched[*v]])));`
98	`}`	98	`}`
99	`}`	99	`}`
100	`}`	100	`}`
101		101
102	`void Harmonics::add_frequency(int es, float scale)`	102	`void Harmonics::add_frequency(int es, float scale)`
103	`{`	103	`{`
104	`if(es<maximum_eigenvalue)`	104	`if(es<maximum_eigenvalue)`
105	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v){`	105	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v){`
106	`Vec3d p = Vec3d(proj[es]);`	106	`Vec3d p = Vec3d(proj[es]);`
107	`double Qval = Q[es][vtouched[*v]];`	107	`double Qval = Q[es][vtouched[*v]];`
108		108
109	`mani.pos(v) += p Qval * scale;`	109	`mani.pos(v) += p Qval * scale;`
110	`}`	110	`}`
111	`}`	111	`}`
112		112
113	`void Harmonics::reset_shape()`	113	`void Harmonics::reset_shape()`
114	`{`	114	`{`
115	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`	115	`for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)`
116	`mani.pos(*v) = Vec3d(0);`	116	`mani.pos(*v) = Vec3d(0);`
117	`}`	117	`}`
118	`void Harmonics::partial_reconstruct(int E0, int E1, float scale)`	118	`void Harmonics::partial_reconstruct(int E0, int E1, float scale)`
119	`{`	119	`{`
120	`for(int es=E0;es<=E1;++es)`	120	`for(int es=E0;es<=E1;++es)`
121	`add_frequency(es, scale);`	121	`add_frequency(es, scale);`
122	`}`	122	`}`
123		123
124		124
125	`double Harmonics::compute_adf(HMesh::VertexAttributeVector<double>& F, double t, double fiedler_boost)`	125	`double Harmonics::compute_adf(HMesh::VertexAttributeVector<double>& F, double t, double fiedler_boost)`
126	`{`	126	`{`
127	`double F_max = 0;`	127	`double F_max = 0;`
-		128	`for(VertexID vid : mani.vertices())`
-		129	`F[vid]=0;`
-		130
128	`for(VertexID vid : mani.vertices()){`	131	`for(VertexID vid : mani.vertices()){`
129	`for(int e = 1; e < V.Length(); ++e)`	132	`for(int e = 1; e < V.Length(); ++e)`
130	`F[vid] += sqr(Q[e][vtouched[vid]]) * exp(-t*V[e]/V[1]);`	133	`F[vid] += sqr(Q[e][vtouched[vid]]) * exp(-t*V[e]/V[1]);`
131	`if(fiedler_boost>0)`	134	`if(fiedler_boost>0)`
132	`for(int e = 1; e < V.Length(); ++e)`	135	`for(int e = 1; e < V.Length(); ++e)`
133	`F[vid] += sqr(Q[e][vtouched[vid]]) * exp(-t) * fiedler_boost;`	136	`F[vid] += sqr(Q[e][vtouched[vid]]) * exp(-t) * fiedler_boost;`
134		137
-		138	`F_max = max(F[vid], F_max);`
-		139	`}`
-		140	`return F_max;`
-		141	`}`
-		142
-		143	`double Harmonics::compute_esum(HMesh::VertexAttributeVector<double>& F, int e0, int e1)`
-		144	`{`
-		145	`double F_max = 0;`
-		146	`for(VertexID vid : mani.vertices())`
-		147	`F[vid]=0;`
-		148
-		149	`for(VertexID vid : mani.vertices()){`
-		150	`for(int e = e0; e < e1; ++e)`
-		151	`F[vid] += sqr(Q[e][vtouched[vid]]);`
-		152
135	`F_max = max(F[vid], F_max);`	153	`F_max = max(F[vid], F_max);`
136	`}`	154	`}`
137	`return F_max;`	155	`return F_max;`
138	`}`	156	`}`
139		157

Subversion Repositories gelsvn

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