| 646 |
janba |
1 |
/*
|
|
|
2 |
* harmonics.cpp
|
|
|
3 |
* GEL
|
|
|
4 |
*
|
|
|
5 |
* Created by J. Andreas Bærentzen on 01/09/08.
|
|
|
6 |
* Copyright 2008 __MyCompanyName__. All rights reserved.
|
|
|
7 |
*
|
|
|
8 |
*/
|
|
|
9 |
|
|
|
10 |
#include <iostream>
|
|
|
11 |
|
|
|
12 |
#include "harmonics.h"
|
|
|
13 |
|
|
|
14 |
#if USE_SPARSE_MATRIX
|
|
|
15 |
#include <arlgsym.h>
|
|
|
16 |
#endif
|
|
|
17 |
|
|
|
18 |
#include <CGLA/Vec3d.h>
|
|
|
19 |
#include <CGLA/Mat2x2d.h>
|
|
|
20 |
#include <LinAlg/Matrix.h>
|
|
|
21 |
#include <LinAlg/Vector.h>
|
|
|
22 |
#include <LinAlg/LapackFunc.h>
|
|
|
23 |
|
|
|
24 |
#include <GL/glew.h>
|
|
|
25 |
#include <GLGraphics/glsl_shader.h>
|
|
|
26 |
|
|
|
27 |
#include <HMesh/mesh_optimization.h>
|
|
|
28 |
#include <HMesh/curvature.h>
|
|
|
29 |
#include <HMesh/triangulate.h>
|
|
|
30 |
#include <HMesh/load.h>
|
|
|
31 |
#include <HMesh/x3d_save.h>
|
|
|
32 |
|
|
|
33 |
using namespace CGLA;
|
|
|
34 |
using namespace std;
|
|
|
35 |
using namespace HMesh;
|
|
|
36 |
using namespace Geometry;
|
|
|
37 |
using namespace LinAlg;
|
|
|
38 |
|
|
|
39 |
|
|
|
40 |
void Harmonics::make_laplace_operator()
|
|
|
41 |
{
|
|
|
42 |
Q.Resize(mani.no_vertices(), mani.no_vertices());
|
|
|
43 |
|
|
|
44 |
for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)
|
|
|
45 |
if(!boundary(mani, *v)){
|
|
|
46 |
int i = vtouched[*v];
|
|
|
47 |
double area_i = voronoi_area(mani, *v);
|
|
|
48 |
Vec3d vertex(mani.pos(*v));
|
|
|
49 |
Vec3d curv_normal(0);
|
|
|
50 |
double a_sum = 0;
|
|
|
51 |
for(Walker wv = mani.walker(*v); !wv.full_circle(); wv = wv.circulate_vertex_cw())
|
|
|
52 |
{
|
|
|
53 |
int j = vtouched[wv.vertex()];
|
|
|
54 |
double area_j = voronoi_area(mani, wv.vertex());
|
|
|
55 |
|
|
|
56 |
Vec3d nbr(mani.pos(wv.vertex()));
|
|
|
57 |
Vec3d left(mani.pos(wv.next().vertex()));
|
|
|
58 |
Vec3d right(mani.pos(wv.opp().prev().opp().vertex()));
|
|
|
59 |
|
|
|
60 |
double d_left = dot(normalize(nbr-left),
|
|
|
61 |
normalize(vertex-left));
|
|
|
62 |
double d_right = dot(normalize(nbr-right),
|
|
|
63 |
normalize(vertex-right));
|
|
|
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)));
|
|
|
66 |
|
|
|
67 |
double w = 1.0/tan(a_left) + 1.0/tan(a_right);
|
|
|
68 |
|
|
|
69 |
Q[i][j] = -w/sqrt(area_i*area_j);
|
|
|
70 |
//Q[i][j] = -1;
|
|
|
71 |
a_sum += Q[i][j];
|
|
|
72 |
}
|
|
|
73 |
Q[i][i] = -a_sum;
|
|
|
74 |
}
|
|
|
75 |
EigenSolutionsSym(Q,V);
|
|
|
76 |
}
|
|
|
77 |
|
|
|
78 |
Harmonics::Harmonics(HMesh::Manifold& _mani):mani(_mani), vtouched(_mani.allocated_vertices(), 0)
|
|
|
79 |
{
|
|
|
80 |
|
|
|
81 |
int i = 0;
|
|
|
82 |
for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v, ++i)
|
|
|
83 |
vtouched[*v] = i;
|
|
|
84 |
maximum_eigenvalue = mani.no_vertices()-1;
|
|
|
85 |
make_laplace_operator();
|
|
|
86 |
|
|
|
87 |
proj.resize(maximum_eigenvalue);
|
|
|
88 |
max_eig_values.resize(maximum_eigenvalue, 1e-10f);
|
|
|
89 |
|
|
|
90 |
cout << "Projection magnitude" << endl;
|
|
|
91 |
for(int es=0; es<maximum_eigenvalue; ++es)
|
|
|
92 |
{
|
|
|
93 |
proj[es] = Vec3d(0.0);
|
|
|
94 |
for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)
|
|
|
95 |
{
|
|
|
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]])));
|
|
|
98 |
}
|
|
|
99 |
}
|
|
|
100 |
}
|
|
|
101 |
|
|
|
102 |
void Harmonics::add_frequency(int es, float scale)
|
|
|
103 |
{
|
|
|
104 |
if(es<maximum_eigenvalue)
|
|
|
105 |
for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v){
|
|
|
106 |
Vec3d p = Vec3d(proj[es]);
|
|
|
107 |
double Qval = Q[es][vtouched[*v]];
|
|
|
108 |
|
|
|
109 |
mani.pos(*v) += p * Qval * scale;
|
|
|
110 |
}
|
|
|
111 |
}
|
|
|
112 |
|
|
|
113 |
void Harmonics::reset_shape()
|
|
|
114 |
{
|
|
|
115 |
for(VertexIDIterator v = mani.vertices_begin(); v != mani.vertices_end(); ++v)
|
|
|
116 |
mani.pos(*v) = Vec3d(0);
|
|
|
117 |
}
|
|
|
118 |
void Harmonics::partial_reconstruct(int E0, int E1, float scale)
|
|
|
119 |
{
|
|
|
120 |
for(int es=E0;es<=E1;++es)
|
|
|
121 |
add_frequency(es, scale);
|
|
|
122 |
}
|
|
|
123 |
|
|
|
124 |
|
|
|
125 |
double Harmonics::compute_adf(HMesh::VertexAttributeVector<double>& F, double t)
|
|
|
126 |
{
|
|
|
127 |
double F_max = 0;
|
|
|
128 |
for(VertexID vid : mani.vertices()){
|
|
|
129 |
for(int e = 1; e < V.Length(); ++e)
|
|
|
130 |
F[vid] += sqr(Q[e][vtouched[vid]]) * exp(-t*V[e]/V[1]);
|
|
|
131 |
F_max = max(F[vid], F_max);
|
|
|
132 |
}
|
|
|
133 |
return F_max;
|
|
|
134 |
}
|